Bone quantitative susceptibility mapping using a chemical species-specific R2* signal model with ultrashort and conventional echo data.

Purpose: To develop quantitative susceptibility mapping (QSM) of bone using an ultrashort echo time (UTE) gradient echo (GRE) sequence for signal acquisition and a bone-specific effective transverse relaxation rate (R (cid:2) 2 ) to model water–fat MR signals for field mapping. Methods: Three-dimensional radial UTE data (echo times (cid:3) 40 m s) was acquired on a 3 Tesla scanner and fitted with a bone-specific signal model to map the chemical species and susceptibility field. Experiments were performed ex vivo on a porcine hoof and in vivo on healthy human subjects (n ¼ 7). For water–fat separation, a bone-specific model assigning R (cid:2) 2 decay mostly to water was compared with the standard models that assigned the same decay for both fat and water. In the ex vivo experiment, bone QSM was correlated with CT. Results: Compared with standard models, the bone-specific R (cid:2) 2 method significantly reduced errors in the fat fraction within the cortical bone in all tested data sets, leading to reduced artifacts in QSM. Good correlation was found between bone CT and QSM values in the porcine hoof (R 2 ¼ 0.77). Bone QSM was successfully generated in all subjects. Conclusions: The QSM of bone is feasible using UTE with a conventional echo time GRE acquisition and a bone-specific R (cid:2) 2 signal model. Magn Reson Med 79:121–128, 2018. 2017 Society for

Quantitative susceptibility mapping of the bone has been challenging because it requires complete measurements of phase everywhere within the region of interest (ROI), and cortical bone typically has very low signal at conventional echo times in gradient echo (GRE) imaging. Although water is abundant in cortical bone ($15% by volume (26)), it mostly exists in the "bound" form, that is, connected to the crystalline mineral structures or the collagen matrix. As a result, bound bone water has an ultrashort apparent transverse relaxation time (T Ã 2 $300 ms (27)), resulting in no meaningful phase for QSM reconstruction on conventional MRI. Because of these limitations, previous work in musculoskeletal applications of QSM was either focused on cartilage, or used piece-wise estimations of bone susceptibility (2,(28)(29)(30)(31)(32)(33)(34). An additional problem arises from intermingling of fat and water protons in the bone marrow, necessitating the application of water-fat separation techniques for field mapping.
The purpose of this preliminary study was to investigate the feasibility of using QSM for measuring bone MRI signal and to highlight the inherent technical issues involved in this application. To that end, we measured bone MRI signal using an ultrashort echo time (UTE) pulse sequence with an echo time (TE) of ( 1 ms (35)(36)(37), and investigated chemical shift and effective transverse relaxation rate (R Ã 2 ) components to properly model bone MRI signal in both UTE and conventional GRE.

THEORY
Magnetic field estimation is an essential initial step for QSM. For QSM of anatomical regions where magnetic susceptibility is the predominant contributor to the proton phase accrual such as in the brain, the field can be estimated from the MRI signal phase. For QSM on regions with high fat content, such as in bone with marrow, a complex MRI signal model is needed for estimating the inhomogeneous field generated by susceptibility through robust separation of water and fat signals (38)(39)(40)(41).
In a multiecho GRE sequence, the temporal behavior of the signal originating from the voxelr that contains multiple species can be expressed in general form as 1 sðr ; tÞ ¼ X k r k ðr Þ Á X n a k n e Ài2pf k n t e ÀR Ã 2 k n t Á e Ài2pfsðr Þt [1] Here, r k is the complex signal originating from the kth chemical species within the voxel at time t ¼ 0; a k n is the relative amplitude of the nth peak in the chemical spectrum of the kth species; f k n and R Ã 2 k n are the corresponding chemical shift and R Ã 2 decay rate for the nth spectral peak of the kth species, respectively; and f s ðr Þ is a spatially varying field induced by susceptibility sources, or the susceptibility field. Extracting multiple parameters from the MRI signal is highly sensitive to noise propagation, and reducing the number of parameters in Equation [1] is critical for robust estimation. We propose the following reduction in parameter numbers. In many clinical applications, water and fat are the dominant species, and Equation [1] becomes Here, a f n is the relative amplitude of the nth spectral peak of fat; r f and r w are fat and water complex amplitude at time t ¼ 0; and R Ã 2w and R Ãf 2n are water and the nth spectral component of fat transverse decay rates.
Although there is some variation in R Ã 2 decay rates among spectral peaks, an approximation is made here that it is constant. This simplifies Equation [2] to [3] In many applications, the difference in water and fat R Ã 2 within a voxelr is commonly neglected (39,40): sðr ; tÞ ¼ r w ðr Þ þ r f ðr Þ X n a f n e Ài2pf f n t Á e ÀR Ã 2 t Á e Ài2pfsðr Þt [4] Unfortunately, in bone tissue, water is distributed through a porous mineral matrix, and fat in the marrow tends to aggregate, and we have R Ã 2bone ) R Ã 2fat . Consequently, the use of Equation [4] to fit the bone MRI signal may lead to very large errors in the estimation of fat fraction, susceptibility field, and QSM. Therefore, we propose to avoid the R Ã 2 simplification in Equation [4]. The water signal is the main component acquired in cortical bone voxels. In our experiments (see below), the largest echo time is much shorter than the T Ã 2 of surrounding soft tissues, including marrow. Therefore, the bone water is the only species experiencing significant decay during acquisition. Accordingly, we further propose a bone-specific R Ã 2 model with the following reduction of parameters:

Pulse sequence
A radial 3D GRE-UTE sequence with TE TE ! 40 ms was implemented on a clinical 3 Tesla (T) scanner (GE Excite HD, GE Healthcare, Milwaukee, WI, USA). The sequence used a nonselective hard pulse (rectangular pulse of 100 ms duration) to achieve volumetric excitation and two readouts per repetition time (TR) to accelerate acquisition ( Fig. 1). On successive TRs, the TEs were shifted to achieve four unique echo times, two of which were considered ultrashort.

Porcine Specimen Experiment
A phantom was constructed from a porcine hoof specimen (length 16 cm, thickness 6.5 cm) embedded in 1% agarose gel (container height 20 cm, average diameter 11 cm) and, for comparison purposes, was imaged using an MR scanner (3T, GE Excite HD) and a clinical CT system ( The CT data were acquired using the following parameters: 120 kVp, 200 mA, 0.625-mm slice thickness, 512 x 512 matrix, achieving nearly isotropic resolution.

Image Processing and Analysis
Each radial MR data set was reconstructed using regridding (42), which interpolates the measured signal from radial spokes onto a Cartesian grid. In this work, NUFFT (43) with Kaiser-Bessel kernel and min-max interpolation was implemented. Because of the nonuniform sampling density of our radial trajectory, density compensation (44) was applied to the measured signal before regridding. Each data set was reconstructed on a coil-by-coil and echo-by echo basis. For each coil, the phase of the first echo was subtracted from the phase of each of the following echoes in order. This allowed a complex coil combination by simply summing the complex data across all coils for each echo. Finally, the phase of this summed signal was used as the phase for the corresponding echo.
Water, fat, and inhomogeneity components were then obtained using iterative decomposition of water and fat with echo asymmetry and least-squares estimation-based techniques (39)(40)(41). For this iterative method, it is important to use an initial guess for the inhomogeneity field that is reasonably close to the true solution, to avoid convergence to a local minimum. A preliminary field map estimation was carried out using simultaneous phase unwrapping and removal of chemical shift using graph cuts (45), assuming a single chemical shift peak with value f ¼ À3:5 ppm Á g B 0 Hz. An initial estimation of the distribution of R Ã 2 values was produced by mono-exponential fitting (46). The complex data were then fitted to Equations [4] and [5], assuming different models for the fat chemical spectrum as follows: single peak model (Df ¼ À3:46 ppm) and multipeak model (Df ¼ fÀ3:82; À3:46; À2:74; À1:86; À0:5; 0:5g ppm) (38). From the resulting inhomogeneity field f s , a susceptibility map was obtained using the morphology-enabled dipole inversion pipeline (16). Projection onto dipole fields (47) was used for background field removal. All susceptibility values reported in this work were referenced to adjacent homogenious muscle tissue.
The CT images of the porcine hoof were resampled and registered to the reconstructed susceptibility map using the FLIRT algorithm in the FMRIB Software Library toolbox (48). The ROI analysis was performed on coregistered volumes to correlate the CT signal (in Hounsfield units, HU) with the MR signal (calculated susceptibility values). For this, 58 in-slice ROIs were manually drawn in flexor tendon, and trabecular and cortical areas of metacarpal and phalanx bones. The QSM and CT volumetric averages for the ROIs were recorded and used for correlation analysis.

Porcine Phantom Results
A susceptibility map of the specimen was successfully reconstructed. Comparison of CT and QSM images is shown in Figure 3, along with the results of the ROI analysis ( Fig. 2) (49). It should be noted that the good correspondence between diamagnetic regions in QSM and regions of high HU values in CT images is supported by a strong linear correlation (correlation coefficient R 2 ¼ 0.77).

Volunteer Results
QSM reconstruction using the proposed method with water-specific R Ã 2 and multipeak fat spectrum modeling was successful in all volunteers (Fig. 4). Figure 5 shows a comparison of the field maps estimated with different techniques and the corresponding calculated susceptibility distributions. Systematic overestimations of the fat fraction within the bone and tendon (up to 80% above the negligible lipid content expected in healthy bone) areas were observed in maps calculated using conventional estimators. This led to errors in the field and, subsequently, the susceptibility map with significant errors in bone and tendon susceptibility values. This was most notable when assuming either a single fat peak or a common R Ã 2 between fat and water. In this case, the cortical bone region of the femur erroneously appeared to be paramagnetic (Fig. 5). It should also be noted that the proposed technique yielded the best visualization of trabecular bone (Figs. 4 and 5). Figure 6 shows a minimum intensity projection of the reconstructed QSM of the knee joint in one healthy subject. Homogeneous diamagnetic appearance of thick cortical areas of the femoral (top yellow arrow) and tibial (bottom yellow arrow) shafts are visible. Throughout the entire FOV, trabeculation is well-depicted by QSM, with a clear appearance of the epiphyseal plate (red arrow) and the area of transition from diaphyseal to metaphyseal bone (blue arrow). The observed nearly homogeneous appearance of the susceptibility throughout the bone is expected, given the likely absence of red marrow in the distal femur in this 28-year-old male volunteer.
The results of the volunteer scans (Table 1) reveal good intra-and intersubject agreement of susceptibility values for cortical bone (intersubject average x bone ¼ À1:4 60:2 ppm).

DISCUSSION
Our data demonstrate the feasibility of mapping bone magnetic susceptibility including cortical bone, trabecular bone, and marrow. This method uses a radial UTE-GRE sequence to acquire the rapidly decaying signal of bone water, and a novel bone-specific R Ã 2 in signal modeling to achieve accurate water-fat separation. Results in an ex vivo porcine hoof show excellent correlation with CT. By comparing QSM reconstructions using the current standard signal model with a common R Ã 2 for both fat and water components, the bone-specific R Ã 2 approach reduces fat fraction estimation errors within cortical bone and generates susceptibility maps in which cortical bone is correctly identified as diamagnetic (7,8). As the highly concentrated calcification is the dominant susceptibility source, QSM allows noninvasive quantitative mapping of bone mineralization without the need for radiation.
Field estimation from UTE and conventional echo GRE data using the classic T Ã 2 iterative decomposition of water and fat with echo asymmetry and least-squares estimation is nontrivial, and is prone to noise propagation and R Ã 2 bias. The robust extraction of biophysical parameters from GRE data require proper signal modeling with a minimal number of parameters. Standard GRE signal models that assign identical R Ã 2 to both water and fat species (39,40) fail to correctly estimate the field, consequently generating bone QSM with grossly erroneous values (such as paramagnetic bone), because rapid signal decay of bone is misinterpreted as dephasing as a result of the high fat fraction. From the echo times acquired in this study, there is negligible R Ã 2 decay of the fat signal, but rapid R Ã 2 decay of signal in cortical bone, because of the collagen-bound water (50). This led us to the bonespecific R Ã 2 model that accurately reflects the underlying tissue physics, attributing signal decay exclusively to the FIG. 5. Comparison of fat and water fraction maps, field maps, QSM, and estimated tissue signal decay rates for different signal models, including single-and multipeak fat spectrum and common and water-only R Ã water component. Both the porcine hoof and human volunteer data show that this bone-specific R Ã 2 modelwhen combined with iterative least-squares fitting-provides reliable estimation of chemical species distribution as well as a biologically plausible field map. Additionally, the application of the multipeak fat spectrum model further improves the quality of the water-fat separation and, consequently, the susceptibility map.
Previous studies show that bone magnetic susceptibility can be estimated by a highly regularized piece-wise inversion of the local magnetic field (2), which may be problematic when there is local variation in tissue susceptibility. A masked iterative calculation of susceptibility can also be used to estimate distribution of strong magnetic sources in tissue (4), but the iteration tends to enforce both uniformity and underestimation of highsusceptibility structures. Above all, these attempts do not properly account for the underlying fat-water biophysics in the MRI signal generation, nor do they acquire signal within bone using UTE, as illustrated in this work. For trabecular bones with densities much lower than those of cortical bones, QSM may be straightforwardly applied using conventional TEs (51).
Bone mineral density (BMD) assessment is central to the diagnosis of osteoporosis (52,53). Currently, BMD is measured using dual-energy x-ray absorptiometry (DXA), which may be confounded by factors including degenerative changes, small bone size, and overlaying anatomic structures (25). Quantitative CT has been proposed to overcome these limitations of DXA, but it requires a much higher xray dose (25). A number of noninvasive, nonionizing radiation-based methods for characterizing bone tissue are also available: Ultrasound can be used to characterize bone and assess osteoporosis, but may only provide correlativenot absolute-bone-density quantification (54). MRI-based techniques, including ZTE and UTE, can be used to detect water bound to the organic matrix of the bone, but bonebound water quantification requires separation of bound and free water signal components of different relaxation rates and signal scaling (50,(55)(56)(57)(58)(59). Conventional MRI may be used to assess bone quality but not mineral density (25).
The MRI-based bone QSM presented here suggests the possibility of BMD assessment without x-ray radiation. Therefore, our results warrant future research in comparing bone QSM with DXA on patients-including those suffering from osteoporosis-to establish an accurate non-x-ray BMD assessment for predicting fracture risk and guiding therapy. A comprehensive comparison of QSM as a biomarker with other nonionizing radiationbased measures falls outside of the scope of this preliminary study, which has focused on establishing the feasibility of using QSM. Additional investigation is required to assess the clinical utility of bone QSM.
The main limitation of including UTE in QSM is the relatively long acquisition time (approximately 13 min in this study). Scan time may be reduced substantially using data acquisition-acceleration strategies in Bayesian MRI (1,(60)(61)(62)(63)(64). Another potential limitation stems from the assumption that water is the only species experiencing significant decay by the time the last echo is acquired. Although this is correct in many practical situations when bone-matrix-bound water is present, this assumption may break down under conditions of strong field inhomogeneity, such as imaging near the abdominal cavity. These field variations can lead to strong spin dephasing and, as a result, rapid R Ã 2 decay in soft tissue regions. In this situation, high-order shimming might be required. An additional limiting factor not taken in the account in the present work is the image blurring inherent in radial imaging induced by field inhomogeneity and chemical shift (65). If Cartesian imaging had been used, the bandwidth chosen in this work would correspond to a shift for fat slightly below one pixel. The image blurring here is expected to affect final QSM results, and the application of off-resonance correction methods appropriate for non-Cartesian acquisition (65,66) may be required. Finally, as has previously been reported, QSM is prone to underestimation of strong FIG. 6. Local field (left) and thin slab maximum intensity projection of knee joint QSM (right), reconstructed using the proposed technique. A minimum intensity projection of the reconstructed QSM of the knee joint is shown in one healthy subject. Note the delineation of cortical areas of the femur and tibia (yellow arrows), the depiction of trabeculation, the epiphyseal line (red arrow), and the area of transition from diaphyseal to metaphyseal bone (blue arrow) in the femur, and reduced contrast between the bone and surrounding joint tissue. susceptibility sources (67), and further adjustments of the imaging protocol, such as matrix size and resolution, may be necessary to minimize error. A comparison of the results obtained with the proposed technique with those previously reported in the literature (6)(7)(8)68) demonstrates possible underestimation of bone magnetic susceptibility. For example, the authors in (7,8) report bone susceptibility (x bone % À2:4 ppm), obtained in an in vitro experiment. An in vivo estimation performed in (2,68) yielded a fairly wide range of values (x bone % À1:8 $ À2:3 ppm), depending on the reconstruction method. This underestimation is most notably observed around the joints (Fig. 6). The exact nature of this underestimation may be related to partialvolume effects (eg, volume fraction of water in the cortical bone) and the aforementioned limitations of the technique, although this issue requires further investigation.
A future direction of this work can be an investigation of the anisotropic properties of bone, tendon, and cartilage susceptibilities. Although strong magnetic anisotropy in bone and cartilage has been previously reported (6,33,34), it has not been accounted for in the signal model in this preliminary study. The use of a scalar susceptibility model even in the magnetically anisotropic cartilage may contribute to errors in the estimated susceptibility of nearby tissue regions.

CONCLUSIONS
Quantitative magnetic susceptibility maps across the entire bone cross section are feasible using a combination of UTE, a conventional TE gradient echo acquisition, and a bone-specific R Ã 2 signal model.