Trace-based correction of breathing-induced field fluctuations in T2*-weighted imaging of the spinal cord

Purpose Spinal cord MRI at ultra-high field is hampered by time-varying magnetic fields associated with the breathing cycle, giving rise to ghosting artifacts in multi-shot acquisitions. Here, we suggest a correction approach based on linking the signal from a respiratory bellows to field changes inside the spinal cord. The information is used to correct the data at the image reconstruction level. Methods The correction was demonstrated in the context of multi-shot T2*-weighted imaging of the cervical spinal cord at 7T. A respiratory trace was acquired during a high-resolution multi-echo gradient-echo sequence, used for structural imaging and quantitative T2* mapping, and a multi-shot EPI time series, as would be suitable for fMRI. The coupling between the trace and the breathing-induced fields was determined by a short calibration scan in each individual. Images were reconstructed with and without trace-based correction. Results In the multi-echo GRE, breathing-induced fields caused severe ghosting in the long echo time images, which led to a systematic underestimation of T2* in the spinal cord. The trace-based correction reduced the ghosting and increased the estimated T2* values. Breathing-related ghosting was also observed in the multi-shot EPI images. The correction largely removed the ghosting, thereby improving the temporal signal-to-noise ratio of the time series. Conclusions Trace-based retrospective correction of breathing-induced field variations can reduce ghosting and improve quantitative metrics in multi-shot structural and functional T2*-weighted imaging of the spinal cord. The method is straightforward to implement and does not rely on sequence modifications or additional hardware beyond a respiratory bellows.


Introduction
Spatial encoding in MRI relies on the assumption that the background magnetic field is homogenous and stable over time. However, the presence of a subject in the scanner gives rise to both static field inhomogeneity and dynamic field fluctuations. Magnetic susceptibility differences between tissue and air, cause local field distortions around interfaces (1). Because breathing changes the air-tissue distribution of the thorax and abdomen, the surrounding field distribution varies periodically with the breathing cycle (2). The time-varying fields cause mis-localisation of signal, which can manifest as image distortion, apparent motion, blurring or ghosting, depending on the sequence. In neuroimaging, breathing-related field fluctuations can be measured as far away as in the brain (3), but are particularly prominent in the spine because of the proximity to the lungs (4,5). Indeed, breathing-induced fields have been identified as one of the major challenges to overcome to achieve robust functional MRI (fMRI) of the spinal cord (6).
One way to address the problem is to track the field variations over time, and use the information for prospective or retrospective correction approaches. One potential tracker is the signal from a respiratory bellows, which indicates the state of the breathing cycle. It has recently been shown that the respiratory trace can predict over 90% of the breathing-induced time-variance of the field in the spinal cord during normal shallow breathing (5). A respiratory trace has previously been used as basis for realtime shim updating in ultra-high field brain imaging (7), and has more recently been explored for realtime shimming of the spinal cord at 3T using a custom-built 24-channel shim coil (8,9). However, realtime shim updating demands specialized hardware, which is not available at most sites.
In this work, we investigate using the signal from a respiratory bellows to retrospectively correct the acquired MR data for breathing-induced field fluctuations in spinal cord imaging. A linear model is used to transform the acquired respiratory trace to field variations inside the spine along the superiorinferior (z) axis. The model parameters are determined in each individual subject using a short calibration set of fast phase-sensitive FLASH acquisitions. The breathing-induced field variations are then estimated from the trace alone during subsequent scans, and the acquired MR signal is demodulated by the corresponding phase variations before image reconstruction. The method can be implemented on standard MR systems, and for any type of sequence. Here we explore the correction in the context of T2*-weighted imaging at 7T. Ultra-high field accentuates the T2*-weighted contrast and allows for higher resolution.
However, T2*-weighted acquisitions are particularly vulnerable to the effects of field fluctuations.
Furthermore, the amplitude of the field variations scale with the strength of the background field. Hence, accurate correction is crucial for optimal image quality in ultra-high field T2*-weighted imaging. While T2* is therefore challenging at ultra-high field, it is also of particular interest in the spine. In structural imaging, T2*-weighting provides excellent gray/white matter contrast in the spine (10)(11)(12) and it is also the contrast underlying blood-oxygen-level-dependent functional imaging. We therefore implement the correction for high-resolution multi-echo gradient-echo acquisitions, used for structural imaging and quantitative T2* mapping, and time series of multi-shot EPI images, intended for functional imaging. We focus on multi-shot acquisitions as they are more robust against static B 0 field inhomogeneity compared to single-shot acquisitions, while being especially affected by dynamic field variations.

Methods
All acquisitions were performed on a whole-body 7T system (Magnetom, Siemens Healthineers, Erlangen, Germany). Imaging was performed with a volume-transmit, 16 (20.4-29.0) kg/m 2 ), in compliance with local ethics guidelines. The volunteers were instructed to breathe regularly at a comfortable pace and to avoid swallowing during the scans. During all acquisitions, the signal from a respiratory bellows placed just below the thorax was acquired. The respiratory trace, ( ), was synchronized with the imaging data via a trigger in the sequence. No low-pass filter was applied to the trace, but the mean offset was subtracted for each scan to remove slow baseline drifts between scans.

Trace-based correction
The trace-based correction pipeline is summarized in Figure 1. Calibration of the individual breathinginduced spatial field profiles was performed as described in Vannesjo et al (5). In brief, a time series of FLASH images (13) (resolution 3.4x2.3x3.0 mm 3 , FOV 144x144 mm 2 , TR 8 ms, TE 4.08 ms, bandwidth 240 Hz/pixel, FA 6˚) of a single sagittal slice through the center of the spinal cord was acquired during normal breathing (Fig. 1a). The volume TR of the acquisition was 344 ms, yielding sufficient temporal resolution to capture the breathing cycle. The phase in each voxel, ( , ), was unwrapped over time and the time-averaged phase, ( ), was subtracted. The resulting phase evolution was divided by the echo time to yield a measure of the field changes over time: A mask covering the spinal cord was manually defined on the FLASH magnitude image. The field within the mask was averaged in the transverse plane, yielding a time-series of measured field offsets, Δ ! ( , ), along the superior-inferior (z) axis (Fig. 1b). A linear model linking the respiratory trace R(t) to the estimated field offset, Δ ! , , at any given z location was assumed: A least-squares fit based on 30 seconds of the calibration data was used to determine Δ !"# from Eq. [2] in each subject. Time points affected by swallowing were excluded from the fit. The measured Δ !"# was then used to estimate the field fluctuations from the respiratory trace acquired during subsequent scans (Fig. 1c).
For the correction, a spatially homogenous field offset within each transverse plane was assumed.
The field offset was assumed to be static during the length of the readout train. In this case the MR signal from the object is modulated by a phase offset, Δ ! , given by ( Fig. 1d): where !" denotes time since the last RF excitation, and represents time over the full length of the MR sequence. The acquired imaging data, , , !" , was demodulated with the estimated phase offset, Δ ! , , !" , for each time point in the acquisition: After demodulation, image reconstruction was performed with an iterative conjugate gradient optimization algorithm with SENSE unfolding (14). For comparison, image reconstruction was also performed without prior demodulation of the phase offset, to yield uncorrected images.

Acquisitions
Data for correction were acquired using two sequences: a multi-echo GRE sequence using a single-line readout at each echo time, suitable for T2*-weighted structural imaging and quantitative T2* mapping, and a single-TE GRE sequence using a segmented EPI readout, suitable for high-resolution functional imaging.
The multi-echo GRE sequence had the following imaging parameters: 24 axial slices, resolution . In order to minimize the achievable TE in light of very short T2* in the spine, the EPI sequences did not include phase correction lines for static ghost removal; this correction was performed using corresponding phase correction lines acquired separately in a phantom. A mask covering the spinal cord was extracted for the multi-shot EPI data, as described above. The temporal signal-tonoise ratio (tSNR) was calculated inside the spinal cord mask with and without correction. The ghosting systematically affects the T2* quantification as shown in Figure 3. The upper half of the plot (Fig 3a-b) shows the quantification results in two different slices from one subject. The mean signal within the spinal cord decays faster in the uncorrected case (Fig. 3b), as the ghosting increasingly scatters the signal from the spinal cord at longer echo times. This results in a systematic under-estimation of the local T2*, as displayed in the T2*-maps (Fig. 3a) and the histogram of measured T2* values within the spinal cord mask (Fig. 3b). All subjects showed a systematic shift towards higher T2* values inside the spinal cord with the correction (Fig. 3c). The median T2* value within the spinal cord mask was 15-23 ms in the uncorrected case and 23-35 ms with correction (Fig. 3d). In locations where the T2* was intrinsically low, the correction did not affect the quantification. This is evident from the upper slice in   Figure 5 shows a sagittal view of the tSNR without and with correction in three subjects (for completeness, we show the subjects not included in Figure 4). The improvement was larger towards lower levels of the spinal cord, where the field fluctuations are stronger. The mean tSNR within the whole spinal cord mask increased by 32% on average over the subjects (range 6-59%), and the tSNR at around level C6 increased by on average 69% (range 10-135%) (Fig. 5b).

Discussion
In this work, we have presented a correction for breathing-induced field fluctuations in spinal cord T2*weighted imaging. The correction utilizes a respiratory trace to remove breathing-related phase instabilities in the acquired MR signal, and can reduce ghosting in multi-shot anatomical and functional acquisitions. The method relies on a short calibration scan, but does not require any additional hardware beyond a respiratory bellows for tracking the breathing state. Dynamic B 0 fields are a considerable challenge in ultra-high field spinal cord imaging (17), and the method could therefore help to make greater use of the benefits of higher field strengths.
Quantitative T2* mapping was one of the demonstrated use cases for the proposed correction. In the uncorrected case, the T2* values in the spinal cord were systematically underestimated. This happens because the phase offset corresponding to a given field offset scales with the echo time, leading to increased ghosting in later echoes. The ghosting scatters the signal from the spinal cord over the image, thereby mimicking T2* signal loss inside the cord. The measured T2* values with correction were in good agreement with values previously reported at 7T (18). The previous study used a navigator echo for phase stabilisation (18,19). Navigators, however, take up time in the sequence and may prolong the minimum achievable TE and/or TR. In the context of T2* mapping, a navigator readout either replaces the first imaging echo or is appended after the last. In the former case, the early part of the decay, which carries information about the proton density of the tissue, is lost. In the latter case, at long TE the navigator may have too low SNR for robust phase extraction, which may produce an unstable correction that can even increase artefacts. Both cases limit the feasible echo times of the imaging readouts, and our navigator-less protocol is thus able to extend the observable part of the T2*-decay. Measuring the full decay provides more information for the exponential fit and could potentially allow for more complicated signal models as compared to a mono-exponential decay.
The second application investigated in this work was multi-shot EPI time series for functional imaging. In brain imaging, functional acquisitions are routinely performed with single-shot EPI. Singleshot EPI is relatively insensitive to time-varying field offsets, which translate into apparent motion in the time series that can be addressed retrospectively with motion correction algorithms. However, fMRI of the spinal cord at 7T has to date been conducted with 3D multi-shot EPI acquisitions (20,21), to reduce distortion and signal dropout due to local static field inhomogeneity. Multi-shot sequences are more susceptible to time-varying field effects, as this causes phase inconsistencies between shots. We here demonstrated that the tSNR of multi-shot EPI can be improved with correction for the breathing-induced fields, especially towards lower levels of the cervical spinal cord. This may improve the sensitivity to detect neural activation in fMRI of the spinal cord, especially in combination with post-processing methods to further reduce the impact of signal variations of physiological origin (22).
The proposed correction method was able to greatly reduce ghosting artefacts, but did not completely eliminate them. Residual artefacts after correction are also frequently observed with phase navigators. A number of potential reasons for incomplete correction can be identified. Firstly, the proposed correction relies on a reproducible and linear relationship between respiratory trace signal and the field state. This is a good approximation during regular shallow breathing, but is less reliable during The trace-based correction method could potentially be expanded to account for the above additional perturbations. For example, slice-dependent anterior-posterior field gradients could be estimated from the sagittal FLASH images, and the phase encoding gradient could be adjusted accordingly in the reconstruction. A non-linear model linking the trace to the field variations may be able to capture a larger range of fluctuations and breathing modes. Furthermore, additional external devices, such as NMR field probes (23) or motion-tracking optical devices (24) could yield more information about both the field state and the motion, which could then be incorporated into the reconstruction model (25). A further improvement to the current implementation of the correction method would be to eliminate the manual steps in the processing of the field calibration data, i.e. the spinal cord mask creation and the exclusion of time points affected by swallowing. This would be a crucial step to allow for integration into standard acquisition protocols. Potentially this could be achieved with the Spinal Cord Toolbox (16) for masking, combined with automatic outlier detection in the field time-courses.