Functional Magnetic Resonance Imaging (fMRI) is typically acquired using gradient-echo sequences with a long echo time at high temporal resolution. Gradient-echo sequences inherently encode information about the magnetic field in the often discarded image phase. We demonstrate a method for processing the phase of reconstructed fMRI data to isolate temporal fluctuations in the harmonic fields associated with respiration by solving a blind source separation problem. The fMRI-derived field fluctuations are shown to be in strong agreement with breathing belt data acquired during the same scan. This work presents a concurrent, hardware-free measurement of respiration-induced field fluctuations, providing a respiratory regressor for fMRI analysis which is independent of local contrast changes, and with potential applications in image reconstruction and fMRI analysis.

Functional magnetic resonance imaging (fMRI) images time-varying fluctuations in MRI signal associated with changes in blood oxygenation and blood oxygenation level dependent (BOLD) contrast (Ogawa et al., 1992). It is used to identify and map regions of increased and decreased blood oxygenation across the brain, which have been correlated with brain function (Ogawa et al., 2000). Its popularity has led it to become commonplace in medical imaging research investigating brain function at rest and in response to stimulus, under a variety of conditions and states. For BOLD fMRI, data are usually acquired with a T2*-weighted single-shot echo-planar scan. Neural activity causes a local increase in cerebral blood flow (CBF). Due to the resulting increase in venous blood oxygenation, veins become less paramagnetic and the field inhomogeneities around the vessels decrease (Rauscher et al., 2005). Consequently, the loss of signal coherence in the vicinity of blood vessels is slowed down, resulting in an increase in the T2* weighted signal (Ogawa et al., 1992). Changes in CBF are not only due to neural activity, for example blood CO2 acts as a strong vasodilator. Adding 5% CO2 to the inhaled air or oxygen causes an increase in CBF by around 50% (Rauscher et al., 2005). With the brain’s oxygen demand remaining unchanged, the venous blood oxygenation increases, even in the absence of neural activity.

Fluctuations in blood CO2 also occur due to variations in the depth and frequency of breathing, which have a direct influence on cerebral blood flow and consequently influence the BOLD signal in the fMRI scan (Chang & Glover, 2009). A breath hold of 30 seconds, for instance, was shown to increase the BOLD signal by 3-5% (Kastrup et al., 1999). Regular respiratory frequency is around 0.2 to 0.4 Hz, with deep breaths occurring at a rate of 0.05-0.15 Hz (Power et al., 2020). Low-frequency fluctuations are particularly problematic in resting-state fMRI (rs-fMRI), as this technique infers functional connections within the brain from correlations in low-frequency fluctuations in the BOLD signal, which are in the range of 0.1 Hz (Biswal et al., 1995). End tidal fluctuations in CO2 occur at around 0.05 Hz and are correlated with BOLD fMRI signal fluctuations (Wise et al., 2004). Birn et al. found that respiration depth varied by 19.1% ± 8.6% during a lexical task and by 17.9% ± 4.2% when subjects were cued to breath at a constant rate and these changes were significantly correlated to BOLD signal changes (Birn et al., 2006).

Respiratory information can be used to model respiratory contributions to the fMRI signal. Regressing out breathing fluctuations during fMRI analysis requires a respiratory measure. External monitoring of respiration using a breathing belt remains the gold standard, and robust methods have been developed to use the measurements from such sensors to remove noise from fMRI acquisitions (Glover et al., 2000; Kasper et al., 2017). Power et al. analyzed 440 rs-fMRI scans from the Human Connectome Project data set and found a wide variety of intra- and inter-subject breathing periodicity (Power et al., 2020; Van Essen et al., 2013). They observed that over the 14.4 minutes scan, respiratory rate decreased along with depth of breathing, whereas the variance increased (Power et al., 2020).

Single-shot gradient echo (GRE) echo-planar imaging (EPI) sequences with low spatial resolution and high temporal resolution are commonly used for fMRI (Glover, 2011; Ogawa et al., 1992). The signal evolution in a GRE sequence can be described by a magnitude which changes according to the aforementioned tissue-specific decay constant T2*, as well as a phase which accumulates with TE and scales with the local resonant frequency offset from the main magnet field B0 (Rauscher et al., 2006; Reichenbach & Haacke, 2011). The local field offset is determined by the distribution of magnetic susceptibility within and outside of the reconstructed image, and as such, GRE sequences are widely used as input to processing techniques including quantitative susceptibility mapping (QSM) and susceptibility weighted imaging (SWI) (Rauscher et al., 2006, 2011; Reichenbach & Haacke, 2011; Wang & Liu, 2015). These image processing techniques use the information about the local magnetic field encoded in the image phase to quantify the local susceptibility distribution and enhance the depiction of vasculature (Ayaz et al., 2011; Reichenbach et al., 1997).

Considering the phase, each reconstructed fMRI volume is a complex valued depiction of the magnetization sampled around the echo time (TE). It has been shown that respiration induces changes in the magnetic field inside the brain of 2-5 Hz (which is ~0.01 ppm), contributing to physiological noise in the fMRI measurement (Van de Moortele et al., 2002; Vannesjo et al., 2018). This field fluctuation amplitude is within the sensitivity range of phase processing routines employed in QSM (Deistung et al., 2017). Local changes in susceptibility have been observed by processing the phase of fMRI acquisitions in an approach known as functional QSM, an emerging technique which provides motivation for the observability of small field fluctuations in the fMRI scan (Sun et al., 2017). While QSM is sensitive to changes in local susceptibility, the measured phase contains contributions from both the local tissue fields and the harmonic background field, which can confound the measurement of susceptibility values in the volume of interest. It is necessary in QSM to subtract the background field to isolate the local susceptibility distribution (Sun & Wilman, 2014). One of the useful properties of the background field is that it is harmonic and therefore does not contain contributions from signals within the imaging volume. For the fMRI scan, this means that any temporal changes in the harmonic field must comprise field changes which are independent of the BOLD signal.

Although temporal fluctuations in the harmonic field are independent of BOLD contrast, they cannot be attributed to a specific source of physiological noise without additional constraints. In this work, the identification of harmonic field fluctuations which result from respiratory motion is approached using a blind source separation (BSS) strategy (Golub & Reinsch, 1970). BSS attempts to isolate signals from multiple sources which have been measured with a single channel. In the case of respiratory field estimation, the fluctuations in the background field are a superposition of a set of fields arising from different susceptibility distributions in the imaging environment and their behavior.

We can model susceptibility sources as dipoles, which have a magnetic field which decays with the cube of the distance from the dipole B(r)1r3, and so any harmonic field existing within the imaging volume is attributable to nearby susceptibility sources. In the case of fMRI and neuroimaging, this is the body tissue, which is relatively static during the imaging protocol and finite in extent, with motion that can be well-described by a small number of terms (Welch et al., 2002). Using a truncated singular value decomposition (SVD), a small set of spatial basis functions and time-varying coefficient vectors can be found which approximate the spatiotemporal evolution of the harmonic field.

We propose to process the phase of fMRI sequences to capture respiration-induced changes in the harmonic field within the imaging volume. The time-course of the harmonic field is approximated using SVD to separate breathing-induced field modulations from other sources of noise (Agrawal et al., 2020). The resulting field modulation can be further projected onto the set of solid harmonic basis functions, which are used to describe the field fluctuations during acquisition that stem from respiration.

We acquired a set of 10 fMRI scans in 7 healthy volunteers (weight: 79.7 ± 12.3 kg, height: 1.80 ± 0.08 m) at 3 T (Philips Elition) using a gradient-echo EPI sequence (duration ~5 minutes) and TR/TE of 1150/30 ms. An isotropic voxel size of 2.5 mm, a parallel imaging factor of 2, and a multi-band factor of 3 (overall R = 6) were employed to achieve a short TR (Pruessmann et al., 1999). Complex data were exported at the console. The size of the reconstructed data was 96 × 96 × 57 × 260 corresponding to 3 spatial dimensions and one temporal dimension respectively. No zero-filling or partial Fourier acquisition was performed. Nine (9) of the scans were acquired during the imaging portion of a study on concussion in ice hockey players, and a control subject was acquired separately. Three of the subjects were scanned twice. The study was approved by the University of British Columbia Clinical Research Ethics Board (H21-00400). All volunteers gave written, informed consent.

Respiration of the subjects was monitored using the scanner's respiratory physiology sensor (500 Hz sampling rate, positioned around lower chest). The control subject was instructed to perform a series of controlled breathing patterns during the scan: A) 10 short breaths B) 10 deep breaths, C) free-breathing, D) breath-hold of 10 seconds to 20 seconds. All other subjects were given no specific instructions regarding respiration.

After acquisition, the resting-state fMRI data were sliced along the temporal dimension to produce 260 volumes which were individually processed (Fig. 1) using the following steps: Laplacian-based phase unwrapping was performed to remove phase discontinuities in the image (Schofield & Zhu, 2003). The local field was estimated from the unwrapped phase with regularized SHARP (Sun & Wilman, 2014). Then, the background harmonic field was estimated by subtracting the local field from the total field.

Fig. 1.

Illustration of the pipeline from raw image phase to harmonic background field. All images are of a single slice taken from the first timepoint of the fMRI time-series. (A) Raw image phase. (B) Image phase after unwrapping with Laplacian unwrapping. (C) The projection of the total field onto the local susceptibility sources within the volume of interest. Note the absence of low-frequency phase modulation across the slice. (D) The isolated harmonic field, obtained by subtraction of the local field in (C) from the total field depicted in (B).

Fig. 1.

Illustration of the pipeline from raw image phase to harmonic background field. All images are of a single slice taken from the first timepoint of the fMRI time-series. (A) Raw image phase. (B) Image phase after unwrapping with Laplacian unwrapping. (C) The projection of the total field onto the local susceptibility sources within the volume of interest. Note the absence of low-frequency phase modulation across the slice. (D) The isolated harmonic field, obtained by subtraction of the local field in (C) from the total field depicted in (B).

Close modal

The estimated harmonic field data from all volumes was reshaped into a M×N matrix P where M is the number of spatial voxels and N is the number of time samples during the acquisition. The field was then decomposed using Singular-Value Decomposition (SVD) and a temporally low-rank representation of the background field was reconstructed by selecting elements on the diagonal of the matrix S to obtain SL (Algorithm 1) (Wold et al., 1987).

Discrimination of the harmonic field component corresponding to respiration v˜r(t) was performed by finding the scaled right singular vector v˜i(t):i[1:5] of the low-rank representation PL with maximal instantaneous variation. This was accomplished by defining v˜r(t) to be the v˜i(t) for which |dv˜i(t)dt| was maximized. A visualization of the separated harmonic field components is shown for an arbitrary region of interest in Figure 2.

Fig. 2.

(A) An example region of interest, depicted in red, overlaid onto the mean of the 1st component of the harmonic field decomposition in an example slice from the control subject. (B) Time-varying magnetic field within the region of interest shown as a dashed line. The colored lines show the first 5 right singular vectors in the time domain, which sum to the raw field evolution (dashed). The black (third) component is clearly periodic and contains most of the temporal variation in the measured magnetic field.

Fig. 2.

(A) An example region of interest, depicted in red, overlaid onto the mean of the 1st component of the harmonic field decomposition in an example slice from the control subject. (B) Time-varying magnetic field within the region of interest shown as a dashed line. The colored lines show the first 5 right singular vectors in the time domain, which sum to the raw field evolution (dashed). The black (third) component is clearly periodic and contains most of the temporal variation in the measured magnetic field.

Close modal

After discrimination of the correct v˜r(t), the resulting field estimated from v˜r(t) at each timepoint was projected onto the set of solid harmonics (up to third order) to produce a matrix of solid harmonic coefficients which describes the time evolution of the respiration-induced field modulation. The zeroth coefficient was chosen as the respiratory regressor for comparison with the breathing belt measurement. Pseudocode for the proposed method is given in Algorithm 1. Processing of the data was performed in MATLAB and took 150 seconds on a 32-core workstation.

graphic
graphic

Algorithm 1: Pseudocode outlining the proposed method. Variables U and V are the lefthand side and righthand side of the SVD of the measurement matrix P respectively, and S is a diagonal matrix containing the singular values of P on its diagonal. V˜ is a matrix containing the right singular vectors of P scaled by the corresponding singular values.

Processing of the breathing-belt respiratory trace was also performed in MATLAB using standard tools available in the PhysIO toolbox (Kasper et al., 2017). The trace was lowpass filtered and de-spiked to remove spurious noise. The breathing-belt respiratory trace was sampled at ~500 Hz, while the field evolution derived using the proposed method was sampled every volume TR (1150 ms, or 0.87 Hz). An offset correction step was performed to ensure that recorded sampling times for both measurements were synchronized. The temporal evolution of the field modulations measured using the proposed method was compared to the breathing-belt respiratory trace by measuring the error in respiratory period and peak/trough locations between the signals. Peak locations were identified and compared across both traces, and overlapping peaks were identified as those present in both traces within a time interval of half of the volume TR (575 ms). As well, visual comparison of the respiratory trace and the field modulations was performed for all subjects, and the timepoints associated with breathing pattern changes were observed and compared with the instructions given to the control subject.

Synthesis of respiratory phase was performed for the control subject using a histogram-based method proposed for RETROICOR, as well as a Hilbert transform method, and applied to both the zeroth order field modulation and the respiratory sensor signal (Glover et al., 2000; Harrison et al., 2021; Noto et al., 2018). The reference respiratory phase and the respiratory phase obtained from the proposed method were compared by measuring the error in respiratory period and peak/trough locations.

Measured field modulations correlated strongly with the respiratory trace obtained from the breathing belt for all subjects. Correspondence with the reference respiratory phase was maintained across shallow, deep, free-breathing, and breath-hold periods, and for the range of respiratory frequencies observed across the subject cohort (Fig. 3). Peaks of the respiration-induced field modulation were well-resolved, with no obvious aliasing in the reconstructed coefficient vectors. For the control subject (Fig. 3, #1), distinct periods of controlled breathing were seen that corresponded directly with subject instruction. The root mean-square error in the measured respiratory period was 0.3 seconds between the ground-truth and proposed method for the control subject.

Fig. 3.

Respiration induced zeroth order field modulation derived from fMRI time-series data using the proposed method (blue) and externally monitored respiratory trace (red, dashed) for all fMRI scans (numbered #1-10). Strong agreement is shown between both the proposed method and the externally monitored respiratory trace. Each trace depicted (for both methods) was normalized to the interval [-1, 1] to facilitate visual comparison. The temporal range and sampling interval (x-axis) were identical for all fMRI scans. Controlled breathing intervals are indicated for subject #1 using the letters A, B, C, and D. Label A denotes shallow breathing, label B denotes deep breathing, label C denotes free-breathing, and label D denotes a breath hold. Breathing was unregulated for the other 9 scans.

Fig. 3.

Respiration induced zeroth order field modulation derived from fMRI time-series data using the proposed method (blue) and externally monitored respiratory trace (red, dashed) for all fMRI scans (numbered #1-10). Strong agreement is shown between both the proposed method and the externally monitored respiratory trace. Each trace depicted (for both methods) was normalized to the interval [-1, 1] to facilitate visual comparison. The temporal range and sampling interval (x-axis) were identical for all fMRI scans. Controlled breathing intervals are indicated for subject #1 using the letters A, B, C, and D. Label A denotes shallow breathing, label B denotes deep breathing, label C denotes free-breathing, and label D denotes a breath hold. Breathing was unregulated for the other 9 scans.

Close modal

The overlap in detected peaks across trace obtained with the proposed method and the breathing belt trace was 94% over all scans. Across all subjects, the root mean-square error in measured respiratory period between the ground-truth and the proposed method was 0.68 seconds or 17% of the average respiratory period. The mean error in peak locations between the ground-truth and proposed method was 0.57 seconds. The mean of the respiratory period measured across all subjects using the proposed method was 4.79 seconds ± 1.6 seconds, and the mean measured with the breathing belt was 4.79 seconds ± 1.6 seconds.

The spectrum of singular values (Fig. 4) obtained from the SVD of P was consistent with the assumption that fluctuations in the harmonic field were low rank, indicating that PL was a good approximation of P. Isolation of the v˜i(t) corresponding to respiration was reliably performed using our method. The L1-norm distribution of the derivatives of the columns of V˜ (Algorithm 1) was found to be consistently dominated by one component for all subjects and was not biased towards the component associated with the largest singular value of P.

Fig. 4.

Singular values of the matrix P demonstrating significant decay beyond the first 5 singular values. The first 5 singular vectors describe 97% of the variance in the data, validating the assumption in the proposed work that the matrix P can be well-characterized by a low-rank approximation.

Fig. 4.

Singular values of the matrix P demonstrating significant decay beyond the first 5 singular values. The first 5 singular vectors describe 97% of the variance in the data, validating the assumption in the proposed work that the matrix P can be well-characterized by a low-rank approximation.

Close modal

Determination of the solid harmonic coefficients describing the temporal and spatial behavior of the respiration-induced field during respiration was possible in all subjects, and the solid harmonic fitting procedure did not result in significant spurious coefficient values across the whole scan duration. This is seen in Figure 5, where the patterns of controlled breathing are evident in all solid harmonic basis function coefficients.

Fig. 5.

Projection of the temporal evolution of the field onto solid harmonics (see legend for definition) for the control subject, obtained using the proposed method. The primary fluctuation observed is in the zeroth order coefficient, however strong correlations are seen across all coefficients, indicating that respiratory induced field fluctuations can likely be described by a single mode. Instructed breathing periods are seen in all coefficients. The zeroth order fluctuations demonstrate a maximal field excursion of 2.5 ppb, however local field excursions due to higher order coefficients can be up to 15 ppb.

Fig. 5.

Projection of the temporal evolution of the field onto solid harmonics (see legend for definition) for the control subject, obtained using the proposed method. The primary fluctuation observed is in the zeroth order coefficient, however strong correlations are seen across all coefficients, indicating that respiratory induced field fluctuations can likely be described by a single mode. Instructed breathing periods are seen in all coefficients. The zeroth order fluctuations demonstrate a maximal field excursion of 2.5 ppb, however local field excursions due to higher order coefficients can be up to 15 ppb.

Close modal

Maximal field excursions due to respiration were observed to be between -5 ppb and +15 ppb within the imaged volume for the control subject, which corresponds to an approximate range of 2.2 Hz at 3 T, consistent with previous measurement of field perturbations due to breathing (Van de Moortele et al., 2002; Vannesjo et al., 2018). The average maximal field excursion over all subjects within the imaging volume was 3.05 Hz ± 1.2 Hz. These excursions were most significant in the periphery of the imaging volume, where the higher order terms in the solid harmonic expansion become dominant, as seen in Figure 6.

Fig. 6.

Illustration of the respiration-induced field modulation during a single respiratory period for the control subject. Each image is a masked depiction of the field in the central axial slice of the imaging volume, sampled during the scan at time points indicated in the lower sub-figure. The maximal field excursion over the respiratory period in the imaging volume is [-5 ppb, +15 ppb].

Fig. 6.

Illustration of the respiration-induced field modulation during a single respiratory period for the control subject. Each image is a masked depiction of the field in the central axial slice of the imaging volume, sampled during the scan at time points indicated in the lower sub-figure. The maximal field excursion over the respiratory period in the imaging volume is [-5 ppb, +15 ppb].

Close modal

Synthesized respiratory phase from the zeroth order solid harmonic fluctuations agreed strongly with the respiratory phase obtained from the breathing belt (Fig. 7). Distinct periods of controlled breathing were seen in both synthesized and reference traces that corresponded directly with subject instruction. Peaks of the respiratory phase were well-resolved, with no obvious aliasing in the reconstructed signal. For the control subject, root mean square error between respiratory period measured with the breathing belt and derived from the fMRI acquisition was 0.31 seconds. Correspondence with the reference respiratory phase was maintained across shallow, deep, and free-breathing periods. Respiratory phase was not accurately determined during the breath-hold period using either the proposed method or the physiology sensor.

Fig. 7.

Respiratory phase derived from fMRI time-series data (blue) and externally monitored respiratory phase (dashed red). The top image was generated using a histogram-based respiratory phase calculation method, and the bottom image was generated from instantaneous phase calculated using the Hilbert transform. Structured breathing patterns are denoted with the letters A through D. Shallow breathing is denoted with the letter A. Deep, slow breathing is denoted with the letter B. Free breathing is denoted with the letter C. A breath hold was conducted during the period denoted with the letter D. Strong agreement is shown between both the fMRI-derived respiratory phase and the externally monitored respiratory phase, with peak and trough locations showing good correspondence.

Fig. 7.

Respiratory phase derived from fMRI time-series data (blue) and externally monitored respiratory phase (dashed red). The top image was generated using a histogram-based respiratory phase calculation method, and the bottom image was generated from instantaneous phase calculated using the Hilbert transform. Structured breathing patterns are denoted with the letters A through D. Shallow breathing is denoted with the letter A. Deep, slow breathing is denoted with the letter B. Free breathing is denoted with the letter C. A breath hold was conducted during the period denoted with the letter D. Strong agreement is shown between both the fMRI-derived respiratory phase and the externally monitored respiratory phase, with peak and trough locations showing good correspondence.

Close modal

We demonstrate the feasibility of hardware-free estimation of respiration-induced field modulations from fMRI image data on a cohort of 7 healthy volunteers. The proposed method is concurrent with the fMRI experiment, requiring no additional or interleaved reference or navigator acquisitions. The method can be applied retrospectively if the image phase is available and facilitates physiological noise correction by providing an independent respiratory regressor in fMRI studies for which external physiology monitoring is not available or corrupted (Birn et al., 2006; Golestani & Chen, 2020; Power et al., 2020).

The proposed method was shown to be applicable to subjects with a wide range of respiratory rates, as seen in Figure 3. However, an important limitation is that the TR must be shorter than half of the shortest respiration period, a necessary condition for alias-free respiratory trace reconstruction. While this limitation was not encountered in the scope of the present work, the proposed method would not be applicable to populations for whom the resting respiratory rate is greater than twice the volume acquisition rate (1/TR). This limitation also indicates that the proposed method cannot be used to measure cardiac pulsation, which is normally an order of magnitude higher in frequency than respiration. Nevertheless, there is significant effort in fMRI research to increase temporal resolution with novel image acquisition and reconstruction strategies, such as 3D spiral readouts for fMRI (Kasper et al., 2022; Moeller et al., 2010; Engel et al., 2021; Narsude et al., 2014, 2016). These strategies have been demonstrated in vivo with a TR as short as 200 ms for the whole brain volume, which is sufficiently short to adequately sample respiratory fluctuations even in populations with very high respiratory frequency, such as newborns or very young children.

Obtaining breathing information using the fMRI data has several advantages over recording respiration with a belt. Using the belt is prone to errors, for example in obese or small participants. In situations of low chest movement, for example in the case of abdominal breathing, respiration may not be adequately detected. Furthermore, in some populations, the belt may be distracting and wearing it may interfere with the fMRI experiment. The breathing information obtained with the proposed method is automatically synchronized with the fMRI experiment and the processing of the data requires no user interaction. The measured error in respiratory period between the ground-truth belt recordings and the proposed method is less than the volume TR of the fMRI scan, and respiratory period group mean and error are identical across both methods. If the fMRI data are usable for fMRI analysis, they are automatically usable for the tracking of breathing, provide the phase information is saved.

The impact of respiratory period measurement error on physiological noise removal performance has not been specifically studied at the timescale corresponding to the 17% (680 ms) error reported for the proposed method, however there is significant literature on the impact of regressor shifts on fMRI analysis. Birn et al. describe a method by which optimally delayed respiratory regressors can be found for each voxel (Birn et al., 2006). The regressor is shifted with respect to the fMRI signal by -10 seconds to +15 seconds, and while an optimal shift does exist for each voxel, strong correlations were found to be relatively insensitive to shifts of a few seconds (Birn et al., 2006). As well, the respiration response function (RRF), which models the response of the fMRI signal in the brain to respiratory changes, has a duration of 25-40 seconds and varies across different regions of the brain (Birn et al., 2008). Given that changes in fMRI signal with respiration are the result of a convolution with the broad RRF, the accuracy of the proposed method is unlikely to be a limiting factor in resting-state analysis. This likely holds for task-based fMRI as well, as the haemodynamic response function (HRF) is 5-7 seconds in duration (Buxton et al., 2004).

The proposed method is physics driven, using the phase of the acquired fMRI images as the input data. A drawback is that it cannot be applied retrospectively to existing fMRI scans if phase information is unavailable. In such situations, estimation of breathing using approaches that work with the magnitude data have been proposed as a solution, including machine learning (ML) approaches such as a neural network, or deep learning (DL) methods (Addeh et al., 2023; Salas et al., 2020). The training of such approaches requires considerable quality assurance and data cleaning by hand, such as the removal of spike artifacts in the breathing belt data. Only about 14% of the data from the Developing Human Connectome Project could be used for training (Addeh et al., 2023; Somerville et al., 2018). In addition, the same authors identified 52.1% of the Young Adult HCP and 18.9% of the Aging HCP breathing data as good quality respiratory signals (Harms et al., 2018; Smith et al., 2013). Power et al. report that any of three measures typically derived from breathing belt data (respiration variation, the envelope of the respiratory trace, and respiratory volume per time) may miss deep breaths (Power et al., 2020). These findings further underline the notion that field fluctuations may be a more reliable measure of breathing than using a belt, especially in pediatric and aging populations. Retrospective application of ML or DL approaches may require the acquisition of additional training data to tune the network to the specific parameters of the fMRI scan. Furthermore, ML or DL approaches which are applied to magnitude data must extract a respiratory regressor from the BOLD time-series, which is the signal of interest in fMRI. This presents the potential for the respiratory regressor to attribute true BOLD contrast fluctuations to respiration, which would be incorrectly removed upon physiological noise correction. While training on a large, augmented dataset can improve generalizability of ML or DL methods, the acquisition of such a dataset can be difficult for some populations, and the possibility of confounding respiratory and true fMRI signals cannot be eliminated.

Higher temporal resolution field estimation has been reported, with a significant compromise with regards to imaging resolution (Zahneisen et al., 2014). It is, in principle, feasible to extend the proposed method to take advantage of the slice groups acquired in simultaneous multi-slice encoding schemes. This allows for the fitting of solid harmonics at the repetition time of each slice group (slice group TR), presenting an increase in temporal sampling rate by at least an order of magnitude (a factor of 19 in the protocol used in the present work). Fitting to the slice group TR was investigated, however it was not found to yield significant improvements over fitting to the volume TR. Further refinement to the multi-slice parallel imaging strategy would likely be needed to achieve different results, for example increasing the number of slice groups and reducing the gap between simultaneously acquired slices.

The choice of basis expansion to describe the field changes during fMRI acquisition was chosen as the set of solid harmonics up to 3rd order, consistent with existing field estimation work, however use of the proposed method for respiratory phase synthesis would only need an expansion of 1st order. The proposed method is fundamentally similar to a physiological noise measurement approach which uses the phase evolution of external field probes to measure harmonic field perturbations (Gross et al., 2017). Such a technique allows for estimation of the same solid harmonic field expansion as the proposed method, and is capable of measurements concurrent with the imaging experiment, albeit at a much higher temporal resolution than the volume or slice TR. The field probe approach has been shown to be accurate and amenable to fMRI sequences, and thus could serve as a reference measurement for further validation of the proposed method.

Global changes in blood oxygenation, for example due to the inhalation of CO2, can result in field changes of zeroth order. The susceptibility difference between fully oxygenated and fully deoxygenated blood is 264 ppb (Spees et al., 2001). Breathing of 5% CO2 changes venous blood oxygenation from ~0.5 to ~0.7 (Sedlacik et al., 2008). The venous volume fraction in the brain is ~1.7 to 2.5%, because around two thirds of total blood volume resides in veins (Leenders et al., 1990). Then, the global change in susceptibility due to the inhalation of 5% CO2 is between 0.75 and 1.5 ppb, which is about an order of magnitude smaller than the zeroth solid harmonic term.

Recently, a method for physiological regressor estimation on a slice-TR timescale has been proposed (PREPAIR), which identifies time-series signals from magnitude and phase of the fMRI images and performs power spectrum analysis and filtering to derive both respiratory and cardiac regressors (Bancelin et al., 2023). While this method utilizes the phase of the reconstructed image to aid regressor estimation, it does not comprise estimation of the magnetic field changes due to respiration but uses the phase as a source of physiological noise information. While it has been demonstrated as a robust method for regressor estimation and is notably capable of extracting a cardiac regressor sampled at slice-TR, it does not produce an estimate of the magnetic field during the fMRI acquisition. Furthermore, the regressors obtained from PREPAIR are determined to be independent of the BOLD signal by a filtering process in the frequency domain, whereas in the proposed method the estimated respiratory regressor does not require frequency domain filtering.

The computational overhead of the proposed method is minimal, only requiring the phase to be saved at the console and exported to a standard imaging format such as DICOM. The phase must be reconstructed regardless of whether it is exported, so no additional reconstruction time is required during the imaging session. Processing of the exported phase images can be done offline in ~2 minutes for the whole fMRI time series of 260 volumes (Kames et al., 2018). Standard processing pipelines for fMRI data, such as fMRIprep, already have significant computational overhead, and adding the proposed method as an additional step would not be a disruption to fMRI workflows (Esteban et al., 2019). The proposed approach is not limited to fMRI, but can be employed in all scans that acquire gradient echo data at sufficiently high temporal resolution, such as dynamic susceptibility contrast MRI or measurements of cerebrovascular reactivity (Ostergaard, Sorensen, et al., 1996; Ostergaard, Weisskoff, et al., 1996; Fisher et al., 2018). In conclusion, the phase of the fMRI scan is sensitive to subject breathing and can be used for the tracking of respiration, without the need for additional hardware, nor the recording of additional data.

The current implementation of the proposed method is available in MATLAB. Implementations of the phase processing algorithms used in this work are part of the open-source QSM toolboxes and are available in MATLAB and Julia programming languages (https://github.com/kamesy/QSM.jl, https://github.com/kamesy/QSM.m). The code to reproduce the results from this paper is available at https://github.com/alexjaffray/fMRI-phase-toolbox and the data is available upon request.

Alexander Jaffray: Conceptualization, study design, software and algorithm development, data analysis, figure generation, manuscript writing, and interpretation of results. Christian Kames: Software development, manuscript editing, and interpretation of results. Michelle Medina: Figure generation, manuscript editing. Christina Graf: Manuscript editing, data analysis guidance, and interpretation of results. Adam Clansey: Data acquisition, manuscript editing. Alexander Rauscher: Idea and initial concept, study supervision, manuscript writing, and interpretation of results.

The authors have no conflicts of interest to declare.

All work was conducted at the University of British Columbia, located on the traditional, ancestral and unceded territory of the Musqueam people. AR acknowledges support from Canada Research Chairs (950-230363), the Natural Sciences and Engineering Council of Canada (016-05371), and the Canadian Institutes of Health Research (RN382474-418628 and 479213-MPI-CAAA-18330).

Addeh
,
A.
,
Vega
,
F.
,
Medi
,
P. R.
,
Williams
,
R. J.
,
Pike
,
G. B.
, &
MacDonald
,
M. E.
(
2023
).
Direct machine learning reconstruction of respiratory variation waveforms from resting state fMRI data in a pediatric population
.
NeuroImage
,
269
,
119904
. https://doi.org/10.1016/j.neuroimage.2023.119904
Agrawal
,
U.
,
Brown
,
E. N.
, &
Lewis
,
L. D.
(
2020
).
Model-based physiological noise removal in fast fMRI
.
NeuroImage
,
205
,
116231
. https://doi.org/10.1016/j.neuroimage.2019.116231
Ayaz
,
M.
,
Boikov
,
A.
,
McAuley
,
G.
,
Schrag
,
M.
,
Kido
,
D. K.
,
Haacke
,
E. M.
, &
Kirsch
,
W.
(
2011
).
Imaging cerebral microbleeds with SWI
. In
Susceptibility weighted imaging in MRI
(pp.
191
214
).
John Wiley & Sons, Ltd
. https://doi.org/10.1002/9780470905203.ch12
Bancelin
,
D.
,
Bachrata
,
B.
,
Bollmann
,
S.
,
de Lima Cardoso
,
P.
,
Szomolanyi
,
P.
,
Trattnig
,
S.
, &
Robinson
,
S. D.
(
2023
).
Unsupervised physiological noise correction of functional magnetic resonance imaging data using phase and magnitude information (PREPAIR)
.
Human Brain Mapping
,
44
(
3
),
1209
1226
. https://doi.org/10.1002/hbm.26152
Birn
,
R. M.
,
Diamond
,
J. B.
,
Smith
,
M. A.
, &
Bandettini
,
P. A.
(
2006
).
Separating respiratory-variation-related fluctuations from neuronal-activity-related fluctuations in fMRI
.
NeuroImage
,
31
(
4
),
1536
1548
. https://doi.org/10.1016/j.neuroimage.2006.02.048
Birn
,
R. M.
,
Smith
,
M. A.
,
Jones
,
T. B.
, &
Bandettini
,
P. A.
(
2008
).
The respiration response function: The temporal dynamics of fMRI signal fluctuations related to changes in respiration
.
NeuroImage
,
40
(
2
),
644
654
. https://doi.org/10.1016/j.neuroimage.2007.11.059
Biswal
,
B.
,
Yetkin
,
F. Z.
,
Haughton
,
V. M.
, &
Hyde
,
J. S.
(
1995
).
Functional connectivity in the motor cortex of resting human brain using echo-planar MRI
.
Magnetic Resonance in Medicine
,
34
(
4
),
537
541
. https://doi.org/10.1002/mrm.1910340409
Buxton
,
R. B.
,
Uludağ
,
K.
,
Dubowitz
,
D. J.
, &
Liu
,
T. T.
(
2004
).
Modeling the hemodynamic response to brain activation
.
NeuroImage
,
23
,
S220
S233
. https://doi.org/10.1016/j.neuroimage.2004.07.013
Chang
,
C.
, &
Glover
,
G. H.
(
2009
).
Relationship between respiration, end-tidal CO2, and BOLD signals in resting-state fMRI
.
NeuroImage
,
47
(
4
),
1381
1393
. https://doi.org/10.1016/j.neuroimage.2009.04.048
Deistung
,
A.
,
Schweser
,
F.
, &
Reichenbach
,
J. R.
(
2017
).
Overview of quantitative susceptibility mapping: Overview of quantitative susceptibility mapping
.
NMR in Biomedicine
,
30
(
4
),
e3569
. https://doi.org/10.1002/nbm.3569
Engel
,
M.
,
Kasper
,
L.
,
Patzig
,
F.
,
Bianchi
,
S.
,
Prüssmann
,
K. P.
(
2021
).
Whole-brain fMRI at 5 frames per second using T-Hex spiral acquisition
. In
Proceedings of the International Society for Magnetic Resonance in Medicine (ISMRM) 2021
(Vol.
0886
).
Esteban
,
O.
,
Markiewicz
,
C. J.
,
Blair
,
R. W.
,
Moodie
,
C. A.
,
Isik
,
A. I.
,
Erramuzpe
,
A.
,
Kent
,
J. D.
,
Goncalves
,
M.
,
DuPre
,
E.
,
Snyder
,
M.
,
Oya
,
H.
,
Ghosh
,
S. S.
,
Wright
,
J.
,
Durnez
,
J.
,
Poldrack
,
R. A.
, &
Gorgolewski
,
K. J.
(
2019
).
fMRIPrep: A robust preprocessing pipeline for functional MRI
.
Nature Methods
,
16
(
1
),
111
116
. https://doi.org/10.1038/s41592-018-0235-4
Fisher
,
J. A.
,
Venkatraghavan
,
L.
, &
Mikulis
,
D. J.
(
2018
).
Magnetic resonance imaging-based cerebrovascular reactivity and hemodynamic reserve
.
Stroke
,
49
(
8
),
2011
2018
. https://doi.org/10.1161/STROKEAHA.118.021012
Glover
,
G. H.
(
2011
).
Overview of functional magnetic resonance imaging
.
Neurosurgery Clinics of North America
,
22
(
2
),
133
139
. https://doi.org/10.1016/j.nec.2010.11.001
Glover
,
G. H.
,
Li
,
T. Q.
, &
Ress
,
D.
(
2000
).
Image-based method for retrospective correction of physiological motion effects in fMRI: RETROICOR
.
Magnetic Resonance in Medicine
,
44
(
1
),
162
167
. https://doi.org/10.1002/1522-2594(200007)44:1<162::aid-mrm23>3.0.co;2-e
Golestani
,
A. M.
, &
Chen
,
J. J.
(
2020
).
Controlling for the effect of arterial-CO2 fluctuations in resting-state fMRI: Comparing end-tidal CO2 clamping and retroactive CO2 correction
.
NeuroImage
,
216
,
116874
. https://doi.org/10.1016/j.neuroimage.2020.116874
Golub
,
G. H.
, &
Reinsch
,
C.
(
1970
).
Singular value decomposition and least squares solutions
.
Numerische Mathematik
,
14
(
5
),
403
420
. https://doi.org/10.1007/BF02163027
Gross
,
S.
,
Vionnet
,
L.
,
Kasper
,
L.
,
Dietrich
,
B. E.
, &
Pruessmann
,
K. P.
(
2017
).
Physiology recording with magnetic field probes for fMRI denoising
.
NeuroImage
,
154
,
106
114
. https://doi.org/10.1016/j.neuroimage.2017.01.022
Harms
,
M. P.
,
Somerville
,
L. H.
,
Ances
,
B. M.
,
Andersson
,
J.
,
Barch
,
D. M.
,
Bastiani
,
M.
,
Bookheimer
,
S. Y.
,
Brown
,
T. B.
,
Buckner
,
R. L.
,
Burgess
,
G. C.
,
Coalson
,
T. S.
,
Chappell
,
M. A.
,
Dapretto
,
M.
,
Douaud
,
G.
,
Fischl
,
B.
,
Glasser
,
M. F.
,
Greve
,
D. N.
,
Hodge
,
C.
,
Jamison
,
K. W.
, …
Yacoub
,
E.
(
2018
).
Extending the Human Connectome Project across ages: Imaging protocols for the Lifespan Development and Aging projects
.
NeuroImage
,
183
,
972
984
. https://doi.org/10.1016/j.neuroimage.2018.09.060
Harrison
,
S. J.
,
Bianchi
,
S.
,
Heinzle
,
J.
,
Stephan
,
K. E.
,
Iglesias
,
S.
, &
Kasper
,
L.
(
2021
).
A Hilbert-based method for processing respiratory timeseries
.
NeuroImage
,
230
,
117787
. https://doi.org/10.1016/j.neuroimage.2021.117787
Kames
,
C.
,
Wiggermann
,
V.
, &
Rauscher
,
A.
(
2018
).
Rapid two-step dipole inversion for susceptibility mapping with sparsity priors
.
NeuroImage
,
167
,
276
283
. https://doi.org/10.1016/j.neuroimage.2017.11.018
Kasper
,
L.
,
Bollmann
,
S.
,
Diaconescu
,
A. O.
,
Hutton
,
C.
,
Heinzle
,
J.
,
Iglesias
,
S.
,
Hauser
,
T. U.
,
Sebold
,
M.
,
Manjaly
,
Z. M.
,
Pruessmann
,
K. P.
, &
Stephan
,
K. E.
(
2017
).
The PhysIO toolbox for modeling physiological noise in fMRI data
.
Journal of Neuroscience Methods
,
276
,
56
72
. https://doi.org/10.1016/j.jneumeth.2016.10.019
Kasper
,
L.
,
Engel
,
M.
,
Heinzle
,
J.
,
Mueller-Schrader
,
M.
,
Graedel
,
N. N.
,
Reber
,
J.
,
Schmid
,
T.
,
Barmet
,
C.
,
Wilm
,
B. J.
,
Stephan
,
K. E.
, &
Pruessmann
,
K. P.
(
2022
).
Advances in spiral fMRI: A high-resolution study with single-shot acquisition
.
NeuroImage
,
246
,
118738
. https://doi.org/10.1016/j.neuroimage.2021.118738
Kastrup
,
A.
,
Krueger
,
G.
,
Glover
,
G. H.
, &
Moseley
,
M. E.
(
1999
).
Assessment of cerebral oxidative metabolism with breath holding and fMRI
.
Magnetic Resonance in Medicine
,
42
(
3
),
608
611
. https://doi.org/10.1002/(SICI)1522-2594(199909)42:3<608::AID-MRM26>3.0.CO;2-I
Leenders
,
K. L.
,
Perani
,
D.
,
Lammertsma
,
A. A.
,
Heather
,
J. D.
,
Buckingham
,
P.
,
Jones
,
T.
,
Healy
,
M. J. R.
,
Gibbs
,
J. M.
,
Wise
,
R. J. S.
,
Hatazawa
,
J.
,
Herold
,
S.
,
Beaney
,
R. P.
,
Brooks
,
D. J.
,
Spinks
,
T.
,
Rhodes
,
C.
, &
Frackowiak
,
R. S. J.
(
1990
).
Cerebral blood flow, blood volume and oxygen utilization. Normal values and effect of age
.
Brain
,
113
(
1
),
27
47
. https://doi.org/10.1093/brain/113.1.27
Moeller
,
S.
,
Yacoub
,
E.
,
Olman
,
C. A.
,
Auerbach
,
E.
,
Strupp
,
J.
,
Harel
,
N.
, &
Ugurbil
,
K.
(
2010
).
Multiband multislice GE-EPI at 7 tesla, with 16-fold acceleration using partial parallel imaging with application to high spatial and temporal whole-brain fMRI
.
Magnetic Resonance in Medicine
,
63
(
5
),
1144
1153
. https://doi.org/10.1002mrm.22361
Narsude
,
M.
,
Gallichan
,
D.
,
van der Zwaag
,
W.
,
Gruetter
,
R.
, &
Marques
,
J. P.
(
2016
).
Three-dimensional echo planar imaging with controlled aliasing: A sequence for high temporal resolution functional MRI
.
Magnetic Resonance in Medicine
,
75
(
6
),
2350
2361
. https://doi.org/10.1002/mrm.25835
Narsude
,
M.
,
van der Zwaag
,
W.
,
Kober
,
T.
,
Gruetter
,
R.
, &
Marques
,
J. P.
(
2014
).
Improved temporal resolution for functional studies with reduced number of segments with three-dimensional echo planar imaging
.
Magnetic Resonance in Medicine
,
72
(
3
),
786
792
. https://doi.org/10.1002/mrm.24975
Noto
,
T.
,
Zhou
,
G.
,
Schuele
,
S.
,
Templer
,
J.
, &
Zelano
,
C.
(
2018
).
Automated analysis of breathing waveforms using BreathMetrics: A respiratory signal processing toolbox
.
Chemical Senses
,
43
(
8
),
583
597
. https://doi.org/10.1093/chemse/bjy045
Ogawa
,
S.
,
Lee
,
T. M.
,
Stepnoski
,
R.
,
Chen
,
W.
,
Zhu
,
X. H.
, &
Ugurbil
,
K.
(
2000
).
An approach to probe some neural systems interaction by functional MRI at neural time scale down to milliseconds
.
Proceedings of the National Academy of Sciences
,
97
(
20
),
11026
11031
. https://doi.org/10.1073/pnas.97.20.11026
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
.
Proceedings of the National Academy of Sciences
,
89
(
13
),
5951
5955
. https://doi.org/10.1073/pnas.89.13.5951
Ostergaard
,
L.
,
Sorensen
,
A. G.
,
Kwong
,
K. K.
,
Weisskoff
,
R. M.
,
Gyldensted
,
C.
, &
Rosen
,
B. R.
(
1996
).
High resolution measurement of cerebral blood flow using intravascular tracer bolus passages. Part II: Experimental comparison and preliminary results
.
Magnetic Resonance in Medicine
,
36
(
5
),
726
736
. https://doi.org/10.1002/mrm.1910360511
Ostergaard
,
L.
,
Weisskoff
,
R. M.
,
Chesler
,
D. A.
,
Gyldensted
,
C.
, &
Rosen
,
B. R.
(
1996
).
High resolution measurement of cerebral blood flow using intravascular tracer bolus passages. Part I: Mathematical approach and statistical analysis
.
Magnetic Resonance in Medicine
,
36
(
5
),
715
725
. https://doi.org/10.1002/mrm.1910360510
Power
,
J. D.
,
Lynch
,
C. J.
,
Dubin
,
M. J.
,
Silver
,
B. M.
,
Martin
,
A.
, &
Jones
,
R. M.
(
2020
).
Characteristics of respiratory measures in young adults scanned at rest, including systematic changes and “missed” deep breaths
.
NeuroImage
,
204
,
116234
. https://doi.org/10.1016/j.neuroimage.2019.116234
Pruessmann
,
K. P.
,
Weiger
,
M.
,
Scheidegger
,
M. B.
, &
Boesiger
,
P.
(
1999
).
SENSE: Sensitivity encoding for fast MRI
.
Magnetic Resonance in Medicine
,
42
(
5
),
952
962
. https://doi.org/10.1002/(SICI)1522-2594(199911)42:5<952::AID-MRM16>3.0.CO;2-S
Rauscher
,
A.
,
Haacke
,
E. M.
,
Neelavalli
,
J.
, &
Reichenbach
,
J. R.
(
2011
).
Phase and its relationship to imaging parameters and susceptibility
. In
E. M.
Haacke
&
J. R.
Reichenbach
(Eds.),
Susceptibility weighted imaging in MRI
(pp.
47
71
).
John Wiley & Sons, Ltd
. https://doi.org/10.1002/9780470905203.ch4
Rauscher
,
A.
,
Sedlacik
,
J.
,
Barth
,
M.
,
Haacke
,
E. M.
, &
Reichenbach
,
J. R.
(
2005
).
Noninvasive assessment of vascular architecture and function during modulated blood oxygenation using susceptibility weighted magnetic resonance imaging
.
Magnetic Resonance in Medicine
,
54
(
1
),
87
95
. https://doi.org/10.1002/mrm.20520
Rauscher
,
A.
,
Sedlacik
,
J.
,
Deistung
,
A.
,
Mentzel
,
H. J.
, &
Reichenbach
,
J. R.
(
2006
).
Susceptibility weighted imaging: Data acquisition, image reconstruction and clinical applications
.
Zeitschrift Für Medizinische Physik
,
16
(
4
),
240
250
. https://doi.org/10.1078/0939-3889-00322
Reichenbach
,
J. R.
, &
Haacke
,
E. M.
(
2011
).
Gradient echo imaging
. In
E.
Mark Haacke
&
J.
Reichenbach
(Eds.),
Susceptibility weighted imaging in MRI
(pp.
33
46
).
John Wiley & Sons, Ltd
. https://doi.org/10.1002/9780470905203.ch3
Reichenbach
,
J. R.
,
Venkatesan
,
R.
,
Schillinger
,
D. J.
,
Kido
,
D. K.
, &
Haacke
,
E. M.
(
1997
).
Small vessels in the human brain: MR venography with deoxyhemoglobin as an intrinsic contrast agent
.
Radiology
,
204
(
1
),
272
277
. https://doi.org/10.1148/radiology.204.1.9205259
Salas
,
J. A.
,
Bayrak
,
R. G.
,
Huo
,
Y.
, &
Chang
,
C.
(
2020
).
Reconstruction of respiratory variation signals from fMRI data
.
NeuroImage
,
225
,
117459
. https://doi.org/10.1016/j.neuroimage.2020.117459
Schofield
,
M. A.
, &
Zhu
,
Y.
(
2003
).
Fast phase unwrapping algorithm for interferometric applications
.
Optics Letters
,
28
(
14
),
1194
1196
. https://doi.org/10.1364/OL.28.001194
Sedlacik
,
J.
,
Kutschbach
,
C.
,
Rauscher
,
A.
,
Deistung
,
A.
, &
Reichenbach
,
J. R.
(
2008
).
Investigation of the influence of carbon dioxide concentrations on cerebral physiology by susceptibility-weighted magnetic resonance imaging (SWI)
.
NeuroImage
,
43
(
1
),
36
43
. https://doi.org/10.1016/j.neuroimage.2008.07.008
Smith
,
S. M.
,
Beckmann
,
C. F.
,
Andersson
,
J.
,
Auerbach
,
E. J.
,
Bijsterbosch
,
J.
,
Douaud
,
G.
,
Duff
,
E.
,
Feinberg
,
D. A.
,
Griffanti
,
L.
,
Harms
,
M. P.
,
Kelly
,
M.
,
Laumann
,
T.
,
Miller
,
K. L.
,
Moeller
,
S.
,
Petersen
,
S.
,
Power
,
J.
,
Salimi-Khorshidi
,
G.
,
Snyder
,
A. A.
,
Vu
,
A. T.
, …
Glasser
,
M. F.
(
2013
).
Resting-state fMRI in the Human Connectome Project
.
NeuroImage
,
80
,
144
168
. https://doi.org/10.1016/j.neuroimage.2013.05.039
Somerville
,
L. H.
,
Bookheimer
,
S. Y.
,
Buckner
,
R. L.
,
Burgess
,
G. C.
,
Curtiss
,
S. W.
,
Dapretto
,
M.
,
Elam
,
J. S.
,
Gaffrey
,
M. S.
,
Harms
,
M. P.
,
Hodge
,
C.
,
Kandala
,
S.
,
Kastman
,
E. K.
,
Nichols
,
T. E.
,
Schlaggar
,
B. L.
,
Smith
,
S. M.
,
Thomas
,
K. M.
,
Yacoub
,
E.
,
Van Essen
,
D. C.
, &
Barch
,
D. M.
(
2018
).
The Lifespan Human Connectome Project in Development: A large-scale study of brain connectivity development in 5–21 year olds
.
NeuroImage
,
183
,
456
468
. https://doi.org/10.1016/j.neuroimage.2018.08.050
Spees
,
W. M.
,
Yablonskiy
,
D. A.
,
Oswood
,
M. C.
, &
Ackerman
,
J. J. H.
(
2001
).
Water proton MR properties of human blood at 1.5 Tesla: Magnetic susceptibility, T1, T2, T*2, and non-Lorentzian signal behavior
.
Magnetic Resonance in Medicine
,
45
(
4
),
533
542
. https://doi.org/10.1002/mrm.1072
Sun
,
H.
,
Seres
,
P.
, &
Wilman
,
A. H.
(
2017
).
Structural and functional quantitative susceptibility mapping from standard fMRI studies: QSM from fMRI studies
.
NMR in Biomedicine
,
30
(
4
),
e3619
. https://doi.org/10.1002/nbm.3619
Sun
,
H.
, &
Wilman
,
A. H.
(
2014
).
Background field removal using spherical mean value filtering and Tikhonov regularization
.
Magnetic Resonance in Medicine
,
71
(
3
),
1151
1157
. https://doi.org/10.1002/mrm.24765
Van de Moortele
,
P. F.
,
Pfeuffer
,
J.
,
Glover
,
G. H.
,
Ugurbil
,
K.
, &
Hu
,
X.
(
2002
).
Respiration-induced B0 fluctuations and their spatial distribution in the human brain at 7 Tesla
.
Magnetic Resonance in Medicine
,
47
(
5
),
888
895
. https://doi.org/10.1002/mrm.10145
Van Essen
,
D. C.
,
Smith
,
S. M.
,
Barch
,
D. M.
,
Behrens
,
T. E. J.
,
Yacoub
,
E.
, &
Ugurbil
,
K.
(
2013
).
The WU-Minn Human Connectome Project: An overview
.
NeuroImage
,
80
,
62
79
. https://doi.org/10.1016/j.neuroimage.2013.05.041
Vannesjo
,
S. J.
,
Miller
,
K. L.
,
Clare
,
S.
, &
Tracey
,
I.
(
2018
).
Spatiotemporal characterization of breathing-induced B0 field fluctuations in the cervical spinal cord at 7T
.
NeuroImage
,
167
,
191
202
. https://doi.org/10.1016/j.neuroimage.2017.11.031
Wang
,
Y.
, &
Liu
,
T.
(
2015
).
Quantitative susceptibility mapping (QSM): Decoding MRI data for a tissue magnetic biomarker
.
Magnetic Resonance in Medicine
,
73
(
1
),
82
101
. https://doi.org/10.1002/mrm.25358
Welch
,
E. B.
,
Manduca
,
A.
,
Grimm
,
R. C.
,
Ward
,
H. A.
, &
Clifford
,
R. J.
, Jr
. (
2002
).
Spherical navigator echoes for full 3D rigid body motion measurement in MRI
.
Magnetic Resonance in Medicine
,
47
(
1
),
32
41
. https://doi.org/10.1002/mrm.10012
Wise
,
R. G.
,
Ide
,
K.
,
Poulin
,
M. J.
, &
Tracey
,
I.
(
2004
).
Resting fluctuations in arterial carbon dioxide induce significant low frequency variations in BOLD signal
.
NeuroImage
,
21
(
4
),
1652
1664
. https://doi.org/10.1016/j.neuroimage.2003.11.025
Wold
,
S.
,
Esbensen
,
K.
, &
Geladi
,
P.
(
1987
).
Principal component analysis
.
Chemometrics and Intelligent Laboratory Systems
,
2
(
1
),
37
52
. https://doi.org/10.1016/0169-7439(87)80084-9
Zahneisen
,
B.
,
Assländer
,
J.
,
LeVan
,
P.
,
Hugger
,
T.
,
Reisert
,
M.
,
Ernst
,
T.
, &
Hennig
,
J.
(
2014
).
Quantification and correction of respiration induced dynamic field map changes in fMRI using 3D single shot techniques: Respiration Induced Field Map Dynamics
.
Magnetic Resonance in Medicine
,
71
(
3
),
1093
1102
. https://doi.org/10.1002/mrm.24771
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International (CC BY 4.0) license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. For a full description of the license, please visit https://creativecommons.org/licenses/by/4.0/legalcode.