* Сепнет: Deep Neural Network for Magnetic Susputable Security Serviction.

Introduction

In the brain, iron and myelin are the primary magnetic susceptibility sources, predominantly determining tissue susceptibility. These two substances have been shown to be involved in normal brain functions, and changes in their concentrations are linked to the pathophysiology of various neurodegenerative disorders. For instance, iron deposition in substantia nigra is implicated in patients with Parkinson's disease (PD) (Dexter et al.  1987 ), and iron accumulation is also observed in the brains with Alzheimer's disease (AD) (Connor et al.  1992 ). Multiple sclerosis (MS) and Neuromyelitis Optica are known to induce demyelination (Martino and Hartung  1999 ; Wingerchuk et al.  2007 ).

Quantitative susceptibility mapping (QSM) is an MRI approach for quantifying bulk tissue magnetic susceptibility from magnetic field shift (or resonance frequency shift) (de Rochefort et al.  2008 ; Shmueli et al.  2009 ; Wang and Liu  2015 ). Due to its sensitivity to paramagnetic iron and diamagnetic myelin, QSM has been utilized to monitor the brain pathology of neurodegenerative disorders associated with these susceptibility sources. In many brain substructures, however, iron and myelin co‐exist within a voxel (Fukunaga et al.  2010 ), and influence QSM signal collectively. This makes it challenging to estimate their individual contributions to QSM, leading to a limited specificity for each substance.

To address this limitation, magnetic susceptibility source separation has been proposed, enabling separate estimation of the paramagnetic and diamagnetic susceptibility concentrations in brain tissue (Lee et al.  2017 ; Shin et al.  2021 ). The work by Shin et al. introduced χ‐ separation (or chi‐separation) technique to disentangle the contributions of para‐ and diamagnetic susceptibility sources ( χ para and χ dia ) using both frequency shift and reversible transverse relaxation R 2 ′
( = R 2 * − R 2 ) . This approach is based on the model where R 2 ′ is considered proportional to the absolute sum of paramagnetic and diamagnetic susceptibility concentrations (Lee et al.  2017 ; Shin et al.  2021 ). By successfully generating susceptibility source‐separated maps, χ‐ separation offers unique insights beyond traditional QSM, enabling the identification of MS lesions with iron and myelin alterations (Ji et al.  2024 ; Kim et al.  2023 ; Müller et al.  2024 ; Zhu et al.  2024 ) as well as visualization of cortical iron and myelin profiles (Lee et al.  2023 ). Similar susceptibility source separation approaches have been developed (Chen et al.  2021 ; Dimov et al.  2022 ; Emmerich et al.  2021 ; Kan et al.  2024 ; Li et al.  2023 ) and utilized in various applications (Jang et al.  2024 ; Nakashima et al.  2024 ). Chen et al. developed a multi‐exponential model that only required multi‐echo GRE data for susceptibility source separation (Chen et al.  2021 ). Dimov et al. proposed to approximate R 2 ′ from a linearly scaled R 2 * to separate susceptibility sources without R 2 measurement (Dimov et al.  2022 ). Li et al. presented a method that suggested a voxel‐specific spatially adaptive parameter relating R 2 ′ to susceptibility (Li et al.  2023 ). Recently, Kan et al. proposed a spatial optimization for the adaptive relaxometric constant relating R 2 * to susceptibility (Kan et al.  2024 ).

Most of these studies either require multi‐echo spin‐echo (SE) for R 2 measurements, which is time‐consuming, or rely on approximations that can introduce estimation errors (Ji et al.  2024 ). Additionally, similar to QSM, the susceptibility source‐separated maps reconstructed from single‐head orientation data may suffer from streaking artifact due to the ill‐conditioned dipole inversion process (de Rochefort et al.  2008 ; Shmueli et al.  2009 ), even though regularization algorithms can reduce the artifacts to some extent (Li et al.  2015 ; Liu et al.  2011 ). Recently, studies have demonstrated that the application of deep learning algorithms can effectively suppresse such artifacts and create high‐quality maps (Bollmann et al.  2019 ; Yoon et al.  2018 ). Furthermore, deep learning has been successfully applied to MRI contrast synthesis when input contains information, such as structural details and image contrast, required to synthesize output (e.g., synthesizing STIR images from T 1 ‐ and T 2 ‐weighted images) (Hagiwara et al.  2019 ; Moya‐Sáez et al.  2021 ; Tanenbaum et al.  2023 ).

In this work, we present two χ ‐separation pipelines, χ‐ sepnet‐ R 2 ′ and χ‐ sepnet‐ R 2 * , that generate χ ‐separation maps using a deep neural network ( χ‐ sepnet). The difference between the two pipelines is input: χ‐ sepnet‐ R 2 ′ demands both multi‐echo GRE data and multi‐echo SE data (or TSE), while χ‐ sepnet‐ R 2 * requires only multi‐echo GRE data. The missing R 2 ′ information in the χ‐ sepnet‐ R 2 * pipeline is generated by another deep neural network, R2PRIMEnet, reconstructing an R 2 ′ map from an R 2 * map. The proposed pipelines are validated using in vivo data from both healthy subjects and MS patients to assess their feasibility for clinical utility.

Materials and Methods

Datasets for Neural Network Training and Evaluation

Data Acquisition and Processing for Healthy Subjects

To construct χ‐ sepnet, a total of 72 scans was acquired from 12 healthy subjects (age = 25.5 ± 3.6 years; 8 males and 4 females) using a 3 T MRI system (Siemens Tim Trio, Erlangen, Germany) equipped with a 32‐channel phased‐array head coil. All participants signed written consent approved by the institutional review board. Each subject underwent three‐plane localization, manual B 0 shimming, magnetization‐prepared rapid gradient echo (MPRAGE), 2D multi‐echo SE, and then six 3D multi‐echo GRE scans for different head orientations, starting with natural head orientation. For each GRE scan, the three‐plane localization and manual B 0 shimming scan were re‐conducted to ensure field homogeneity after head rotation. For all scans, the imaging slab was axial so that the z‐axis was oriented along the B 0 field. Raw k‐space data were collected from the scanner and reconstructed using in‐house software. The total scan time was approximately 1 h 10 min.

The MPRAGE sequence had the following parameters: repetition time (TR) = 2400 ms, echo time (TE) = 2.10 ms, field of view (FOV) = 256 × 256 × 224 mm 3 , voxel size = 1 × 1 × 1 mm 3 , flip angle = 8°, bandwidth = 210 Hz/pixel, GRAPPA factor = 2, and time of acquisition = 5.45 min. The 2D multi‐echo SE data acquisition was performed using the following scan parameters: TR = 7800 ms, TEs = 15, 30, 45, 60, 75, and 90 ms, FOV = 256 × 204 mm 2 , number of slices = 76, slice gap = 0%, voxel size = 1 × 1 mm 2 , slice thickness = 2 mm, flip angle = 90° with refocusing angle = 180°, bandwidth = 110 Hz/pixel, GRAPPA factor = 2, and time of acquisition = 13.25 min. The 3D multi‐echo GRE data were acquired with the following parameters: TR = 38 ms, TEs = 7.70 ms, 12.73 ms, 17.76 ms, 22.79 ms, 27.82 ms, and 32.85 ms, FOV = 256 × 224 × 176 mm 3 , voxel size = 1 × 1 × 1 mm 3 , flip angle = 15°, bandwidth = 290 Hz/pixel, GRAPPA factor = 2, and time of acquisition = 6.42 min.

From each multi‐echo GRE data, a local field map was calculated adhering to the QSM consensus guidelines (QSM Consensus Organization Committee et al.  2024 ). First, the multi‐channel data were combined using sum‐of‐squares for the magnitude, and MCPC‐3D for the phase (Robinson et al.  2011 ). A brain mask was then extracted from the combined magnitude image using FSL BET (Smith et al.  2004 ). Each phase image was unwrapped using ROMEO (Dymerska et al.  2021 ), followed by background field removal with V‐SHARP, utilizing a spherical mean value size of 20 (Schweser et al.  2011 ; Wu et al.  2012b ). Then, the background‐removed multi‐echo phase images were combined using the weighted echo averaging method to generate a local field map (Wu et al.  2012a ). To create a QSM map from the multi‐orientation GRE data, the local field maps from each orientation were registered to the first orientation data using FSL FLIRT (Smith et al.  2004 ), and the COSMOS QSM reconstruction was performed (Liu et al.  2009 ).

For R 2 * mapping, auto‐regression on linear operations (ARLO) (Pei et al.  2015 ) was utilized for multi‐echo GRE data in each head orientation, enabling fast mono‐exponential fitting. For R 2 mapping, a dictionary of spin‐echo decay was simulated and matched for the multi‐echo SE magnitude image with the stimulated echo correction (McPhee and Wilman  2015 ). For the simulation, the StimFit toolbox ( https://github.com/rmlebel/StimFit ) was utilized. Subsequently, the R 2 * map of each orientation and the R 2 map were registered to the first head orientation R 2 * map using FSL FLIRT (Smith et al.  2004 ). The R 2 ′ maps were generated by subtracting the registered R 2 map from each orientational R 2 * map. Any negative values were set to zero.

The registered multi‐orientation R 2 ′ and local field maps were inputted to the multi‐orientation χ‐ separation algorithm ( χ‐ sep‐COSMOS) (Shin et al.  2022 ), generating a paramagnetic susceptibility map ( χ para ) and a diamagnetic susceptibility map with absolute values ( χ dia ) while the bulk magnetic susceptibility map equals to χ para − χ dia . In this algorithm, any negative values in χ para and χ dia maps were set to zero because they were defined to have non‐negative values The relaxometric constant ( D r ) was estimated by performing linear regression between the R 2 ′ and COSMOS QSM values in the five deep gray matter regions: caudate nucleus, putamen, globus pallidus, red nucleus, and substantia nigra. The resulting relaxometric contrast was 114 Hz/ppm with R 2 of 0.93, slightly deviating from the previously measured value of 137 Hz/ppm potentially due to different QSM algorithms (Shin et al.  2021 ).

Data Acquisition and Processing for Patients

To evaluate the two deep learning pipelines in pathological conditions, MRI data from 12 MS patients (age: 32.3 ± 8.3 years; disease duration: 56.7 ± 35.6 months; five males and seven females; 11 patients treated with immunomodulatory drugs, one untreated; see Table  S1 for details), utilized in a previous study (Ji et al.  2024 ) were reprocessed.

The imaging protocol for the patient datasets included single‐head orientation 3D multi‐echo GRE, 2D dual‐echo turbo spin‐echo (DE‐TSE), and 3D T 2 ‐weighted FLAIR, acquired using a 3 T MRI system (Magnetom Vida, Siemens Healthineers, Erlangen, Germany) equipped with a 64‐channel phased‐array head and neck coil. The 3D multi‐echo GRE data were acquired with the following parameters: TR = 30 ms, TEs = 6.2, 11.8, 17.3, and 22.9 ms, FOV = 187 × 230 × 144 mm 3 , voxel size = 0.7 × 0.7 × 1 mm 3 , flip angle = 18°, bandwidth = 270 Hz/pixel, GRAPPA factor = 2, and time of acquisition = 9.08 min. The 2D DE‐TSE data acquisition was performed using the following scan parameters: TR = 11,000 ms, TE = 10, and 100 ms, FOV = 210 × 210 mm 2 , number of slices = 120, slice gap = 0%, voxel size = 0.7 × 0.7 mm 2 , slice thickness = 2 mm, flip angles = 90°, 165°, and 150°, bandwidth = 222 Hz/pixel, GRAPPA factor = 2, and time of acquisition = 6.05 min. The 3D T 2 ‐weighted FLAIR sequence had the following parameters: TR = 7600 ms, TE = 431 ms, FOV = 230 × 230 × 180 mm 3 , voxel size = 0.7 × 0.7 × 0.7 mm 3 , flip angle = 90° with variable refocusing angles, bandwidth = 651 Hz/pixel, GRAPPA factor = 2 × 2, and time of acquisition = 2.56 min.

The multi‐echo GRE data were processed to produce a local field map and R 2 * map in accordance with the methods outlined in Section  2.1.1 . The DE‐TSE data were processed for an R 2 map using dictionary matching with the StimFit toolbox ( https://github.com/rmlebel/StimFit ), assuming a normalized B1+ value of one throughout the brain (Lebel and Wilman  2010 ; McPhee and Wilman  2015 ). The R 2 maps of the patients were registered to the multi‐echo GRE using FSL FLIRT (Smith et al.  2004 ). The R 2 ′ map was calculated by subtracting the R 2 map from R 2 * map and subsequently setting negative values to zero.

For non‐axial input data, the orientation was rotated such that the B 0 orientation is along the z‐direction.

χ ‐sepnet‐ R 2 ′ and χ ‐sepnet‐ R 2 * Pipelines

The two processing pipelines, χ‐ sepnet‐ R 2 ′ and χ‐ sepnet‐ R 2 * , are shown in Figure  1 . The two pipelines used the same data processing and deep learning‐based χ ‐separation neural network ( χ‐ sepnet), except for the procedures to estimate R 2 ′ . The χ‐ sepnet‐ R 2 ′ pipeline consisted of the pre‐processing steps of the multi‐echo GRE and multi‐echo SE data, producing local field, R 2 * and R 2 maps. Then, the local field map was converted to a QSM map via a neural network, QSMnet ME (see Section  2.2.2 ). Finally, the QSM, local field, and R 2 ′ ( = R 2 * − R 2 ) maps were processed by χ‐ sepnet (see Section  2.2.1 ) to generate paramagnetic ( χ para ) and diamagnetic susceptibility maps ( χ dia ; positive‐valued as mentioned in Section  2.1.1 ) as the output. For the χ‐ sepnet‐ R 2 * pipeline, only the GRE data were used as input, eliminating the need for the multi‐echo SE data (Figure  1b ). After the pre‐processing for local field and R 2 * maps, QSM and R 2 ′ maps were generated by QSMnet ME and R2PRIMEnet, respectively. The latter was a newly designed network to extract R 2 ′ from R 2 * (see Section  2.2.3 ). Finally, the QSM, local field, and R 2 ′ maps were processed via χ‐ sepnet, creating χ para and χ dia .

Overview of the proposed data processing pipelines for χ ‐separation with two different input options: χ‐ sepnet‐ R 2 ′ (a) and χ‐ sepnet‐ R 2 * (b). (a) The χ‐ sepnet‐ R 2 ′ pipeline requires multi‐echo GRE data for local field and R 2 * maps, and multi‐echo SE (or TSE) data for an R 2 map. (b) The χ‐ sepnet‐ R 2 * pipeline is designed to operate solely on multi‐echo GRE data. An R 2 ′ map is synthesized from R 2 * using an additional neural network, R2PRIMEnet.

χ ‐sepnet

χ‐ sepnet was designed to reconstruct χ para and χ dia from the inputs of QSM, local field, and R 2 ′ maps. The three inputs instead of two inputs (i.e., local field and R 2 ′ maps) helped to improve the generalization of the network for untrained resolution (see ablation study in Section  2.3 ). The three input maps were concatenated along the channel dimension. The output maps were designed to produce χ para and χ dia maps. For the network, a 3D U‐net architecture, as utilized in QSMnet (Yoon et al.  2018 ), was employed with modifications to set the input and output channel numbers to 3 and 2, respectively (see Section  S1 for more details). Additionally, the output values were forced to be non‐negative, setting any negative to zero.

For the training of χ‐ sepnet, the χ‐ sep‐COSMOS maps from Section  2.1.1 . were utilized as the label. Data acquired from six healthy subjects (five for training, and one for validation) were utilized. From each of the six head orientations, an input and output pair was generated by registering the χ ‐sep‐COSMOS maps to the orientation of the input maps. Then, the pair was cropped into 3D patches with a size of 64 × 64 × 64 for patch‐based training. Each patch had 40 voxels of overlap with its adjacent neighbours. Data were augmented by rotating the input and output pair at six different angles per subject, creating 20,160 total patches for training after data augmentation. Unlike the augmentation in QSMnet, the rotation was only applied within the plane perpendicular to B 0 (degree: −90° to +90°) to avoid complication from B 0 orientation dependent R 2 * in white matter (Bender and Klose  2010 ; Lee et al.  2011 ; Oh et al.  2013 ).

Before the training, the input and label data (e.g., local field in QSMnet ME and χ para in χ‐ sepnet) were normalized by calculating the mean and standard deviation values of all the training subjects within the brain mask. For each data, the normalization was performed by subtracting the mean value and then dividing by the standard deviation. The process was applied for all training, validation, and test datasets utilized for the neural networks.

The loss functions of the network were defined as follows: reconstruction loss ( L recon χ sep ), gradient loss ( L gradient χ sep ), and model loss ( L model χ sep ). To compute L recon χ sep , L1 loss was calculated between outputs and labels:
(1)
L recon χ sep = g χ φ , R 2 ′ para − χ para 1 + g χ φ , R 2 ′ dia − χ dia 1 ,
where g · para and g · dia denote the paramagnetic and diamagnetic susceptibility maps generated by χ‐ sepnet, respectively. χ , φ , and R 2 ′ refer to the QSM, local field, and R 2 ′ maps, respectively. χ para and χ dia are the labels that correspond to para‐ and diamagnetic susceptibility maps reconstructed by χ‐ sep‐COSMOS, respectively. The gradient loss was designed to enhance the edge information, which was defined as:
(2)
L gradient χ sep = ∇ g χ φ , R 2 ′ para x − ∇ χ para x 1 + ∇ g χ φ , R 2 ′ dia x − ∇ χ dia x 1 + ∇ g χ φ , R 2 ′ para y − ∇ χ para y 1 + ∇ g χ φ , R 2 ′ dia y − ∇ χ dia y 1 + ∇ g χ φ , R 2 ′ para z − ∇ χ para z 1 + ∇ g χ φ , R 2 ′ dia z − ∇ χ dia z 1 .

Lastly, the model loss enforced χ‐ sepnet to incorporate the physics information, represented by three losses ( L QSM χ sep , L field χ sep , and L R 2 ′ χ sep ). These losses have similar ranges because the maps were normalized. The three losses were defined as:
(3)
L QSM χ sep = g χ φ , R 2 ′ para − g χ φ , R 2 ′ dia − χ 1 ,


(4)
L field χ sep = d * g χ φ , R 2 ′ para − g χ φ , R 2 ′ dia − φ 1 ,


(5)
L R 2 ′ χ sep = D r · g χ φ , R 2 ′ para + g χ φ , R 2 ′ dia − R 2 ′ 1 ,


(6)
L model χ sep = L QSM χ sep + L field χ sep + L R 2 ′ χ sep ,
where d represents a dipole kernel, * convolution operation, and D r a relaxometric constant. The final loss function was the weighted sum of the reconstruction, gradient, and model losses with the weights set to 1, 0.1, and 1, respectively.

For training, the RMSprop optimizer was utilized with a learning rate of 0.0003. The stepLR learning rate scheduler was employed with a step size of 1000 and a gamma value of 0.98. The network was trained with a batch size of 12. The training was terminated at 60 epochs because the performance reached a plateau. The final model was selected as the one with the lowest normalized root‐mean‐squared error (NRMSE) on the validation data. The training process was completed in approximately 47 h.

QSMnet ME

QSMnet ME is a neural network to generate a QSM map from a local field map. The network was designed using the same structure as QSMnet (Yoon et al.  2018 ) (see Section  S1 ). The primary difference between the QSMnet ME and QSMnet was the training data. QSMnet ME was trained with COSMOS QSM maps (label) and local field maps (input) from the multi‐echo GRE data described in Section  2.1.1 , which allows for phase offset compensation and provides a higher signal‐to‐noise ratio (Biondetti et al.  2020 ), whereas QSMnet was trained using single‐echo GRE data.

The same subjects and approaches used for the χ‐ sepnet training were employed for the training of QSMnet ME . Only the data augmentation was modified to match that of the QSMnet approach, by rotating the input and label maps at angles relative to B 0 .

Three loss functions were utilized: L1 loss, gradient loss, and model loss, as in QSMnet+ (Jung et al.  2020 ). The loss weights for L1 loss, gradient loss, and model loss were set to 1, 0.1, and 1, respectively.

For training, the batch size was set to 12 and the learning rate was set to 0.001. The RMSprop optimizer and an exponential decay learning rate scheduler with a gamma value of 0.9999 were utilized. The training was conducted for 50 epochs, sufficient for the network's performance to converge, with the final model selected based on the minimum NRMSE on the validation data. The training time was approximately 31 h.

When inferencing an input other than 1 mm isotropic resolution, which was the training data resolution for QSMnet ME , a recently developed resolution generalization method was employed (Ji et al.  2023 ). This approach overcame the issue of resolution generalization in deep learning QSM (Jung, Bollmann, and Lee  2022 ), and enabled us to generate a QSM map up to the resolution of approximately 0.6 mm.

R2PRIMEnet

R2PRIMEnet was designed to estimate an R 2 ′ map from an R 2 * map and was motivated by the previous contrast synthesis studies (Hagiwara et al.  2019 ; Moya‐Sáez et al.  2021 ; Tanenbaum et al.  2023 ). The same network architecture as QSMnet ME was deployed while setting the input and output to R 2 * and R 2 ′ , respectively (see Section  S1 ). The negative values of the output were set to zero.

The same subjects and approaches used for the χ‐ sepnet training were employed for the training of R2PRIMEnet.

The loss of R2PRIMEnet ( L r 2 p ) consists of L1 loss ( L 1 r 2 p ) and gradient loss ( L grad r 2 p ):
(7)
L r 2 p = L 1 r 2 p + L grad r 2 p ,


(8)
L 1 r 2 p = h R 2 * − R 2 ′ 1 ,


(9)
L grad r 2 p = ∇ h R 2 * x − ∇ R 2 ′ x + ∇ h R 2 * y − ∇ R 2 ′ y + ∇ h R 2 * z − ∇ R 2 ′ z .
where h · is R2PRIMEnet, R 2 * is an R 2 * map, and R 2 ′ is an R 2 ′ map. The total loss function was the weighted sum of L 1 r 2 p and L grad r 2 p with the weights set to 1 and 0.1 for L 1 r 2 p and L grad r 2 p , respectively.

For training, the RMSprop optimizer was utilized with a learning rate set at 0.001. A stepLR learning rate scheduler was employed, featuring a step size of 1000 and a gamma value of 0.98. The network was trained with a batch size of 12, and the training was conducted for 50 epochs, which was sufficient for convergence. The total training time was approximately 30 h. The final model was selected for the minimum NRMSE loss on the validation data.

Implementation Environments

The data pre‐processing methods for multi‐echo GRE and multi‐echo SE were conducted using MATLAB (MATLAB 2020a, MathWorks Inc., Natick, MA, USA). The networks were trained using PyTorch (Paszke et al.  2019 ) with two NVIDIA Quadro RTX‐8000 GPUs.

Experiments

To evaluate the proposed pipelines in healthy subjects, the remaining six of the twelve healthy subjects, which were not seen by the neural networks during training, were utilized. Two additional R 2 ′ ‐based susceptibility source‐separation methods ( χ‐ sep‐MEDI and χ‐ sep‐iLSQR), which utilized MEDI regularization with a regularization parameter of 1 (Liu et al.  2011 ) and the iLSQR algorithm (Li et al.  2015 ), were processed (Shin et al.  2021 ) for comparison with the proposed pipelines. The evaluation of different algorithms was conducted by comparing the results to the χ‐ sep‐COSMOS maps as the reference. For quantitative evaluation, peak signal‐to‐noise (pSNR), normalized root‐mean‐squared error (NRMSE), high‐frequency error norm (HFEN), and structural similarity (SSIM) were measured. All the metrics were calculated within the brain mask excluding cerebrospinal fluid (CSF) (Liu et al.  2018 ) and vessels (Kim et al.  2024 ).

A region of interest (ROI) analysis was conducted in a total of 24 brain regions, including gray matter and white matter regions. These ROIs were segmented from the χ‐ separation atlas (see Section  S2 ) (Min et al.  2023 ).

An ablation study was designed to compare the performance of χ ‐sepnet with that of a model excluding the QSM input (i.e., χ ‐sepnet trained with field and R 2 ′ only) or the field input (i.e., χ ‐sepnet trained with QSM and R 2 ′ only). All the networks had the same architecture as χ‐ sepnet, except for the number of input channels (i.e., two input channels instead of three). The model loss term for the local field ( L field χ sep ) or for the QSM ( L QSM χ sep ) was not used for the network trained without local field or QSM, respectively.

An additional study was designed to assess the resolution generalization capability of χ ‐sepnet, which was known to be an important issue in deep learning QSM (Jung, Bollmann, and Lee  2022 ). Since χ ‐sepnet was trained using the 1 mm isotropic resolution data, inputs with different resolutions (2 × 2 × 2 mm 3 , 1 × 1 × 3 mm 3 , and 0.7 × 0.7 × 0.7 mm 3 ) were generated and tested. To generate the data with different resolutions, the k‐space of the 1 × 1 × 1 mm 3 input and label data used for χ ‐sepnet test were cropped (for lower resolution) or zero‐padded (for higher resolution). A 3D Tukey window with a filter size of 0.2 was then used to minimize artifacts in the resolution‐modified maps. A resolution generalization method (Ji et al.  2023 ) was employed for QSMnet ME . The performance of χ ‐sepnet was compared with the two models utilized in the ablation study.

For the 12 MS patients, a qualitative assessment of lesions was conducted by a neuroradiologist (J. J.). Unlike the healthy subject data, the patient data were limited to a single‐head orientation acquisition and, therefore, gold standard reference acquired from multi‐orientation data did not exist. Hence, we only compared the consistency between the χ‐ sepnet‐ R 2 ′ and χ‐ sepnet‐ R 2 * pipelines. MS lesions were segmented using a region growth algorithm from FLAIR images available at SPM12 software ( https://www.fil.ion.ucl.ac.uk/spm/software/spm12/ ) (Egger et al.  2017 ) and fine‐tuned by the neuroradiologist. A total of 250 lesions, each exceeding a volume of 20 mm 3 , were evaluated for visual characteristics. The para‐ and diamagnetic susceptibility maps of the lesions were categorized into five classifications: hypointense, subtly hypointense, isointense, subtly hyperintense, or hyperintense relative to surrounding normal‐appearing white matter (NAWM) (Kim et al.  2023 ). Additionally, paramagnetic rim lesion (PRL) was identified on the paramagnetic susceptibility maps as a subtype of the hyperintense paramagnetic susceptibility lesion.

Results

QSMnet ME
and R2PRIMEnet in Healthy Subjects

The QSM maps generated by QSMnet ME , along with the differences relative to the COSMOS QSM maps, are displayed in Figure  2a . The difference maps excluded CSF and vessels. As demonstrated in the figure, QSMnet ME successfully reconstructed COSMOS‐quality QSM maps. The quantitative metrics for QSMnet ME (pSNR: 33.3 ± 0.9 dB, NRMSE: 48.0% ± 4.3%, HFEN: 36.8% ± 2.8%, SSIM: 0.950 ± 0.011; Figure  2b ) are slightly better than the literature values of QSMnet and QSMnet+ except for pSNR (see Discussion).

(a) QSM maps reconstructed by the COSMOS algorithm and QSMnet ME are shown alongside the absolute difference images. QSMnet ME successfully reconstructed COSMOS‐quality QSM maps. (b) Quantitative metrics of QSMnet ME (pSNR, NRMSE, HFEN, and SSIM), calculated with respect to the COSMOS QSM maps. (c) R 2 ′ maps from MRI measurements (i.e., R 2 * − R 2 ) and R2PRIMEnet, along with the difference images. R2PRIMEnet effectively reconstructed the R 2 ′ maps, showing little differences from the reference. (d) Four metrics of R2PRIMEnet (pSNR, NRMSE, HFEN, and SSIM). pSNR and NRMSE of R2PRIMEnet are better than those from QSMnet ME , while being comparable in SSIM. For pSNR and SSIM, higher values indicate superior performance, whereas for NRMSE and HFEN, lower values are indicative of better performance.

The R 2 ′ maps reconstructed by R2PRIMEnet and the difference images excluding CSF and vessels are displayed in Figure  2c . R2PRIMEnet effectively reconstructed the R 2 ′ maps from R 2 * maps, showing little difference from the reference. The quantitative metrics for R2PRIMEnet (pSNR: 42.7 ± 4.2 dB, NRMSE: 25.3% ± 1.9%, HFEN: 44.9% ± 5.2%, SSIM: 0.949 ± 0.006; Figure  2d ) are better than those of QSMnet ME in terms of pSNR and NRMSE, while exhibiting comparable in SSIM.

χ ‐sepnet‐ R 2 ′ and χ ‐sepnet‐ R 2 * in Healthy Subjects

The χ para and χ dia maps reconstructed by five different χ‐ separation reconstruction algorithms ( χ ‐sep‐COSMOS, χ‐ sepnet‐ R 2 ′ , χ‐ sepnet‐ R 2 * , χ ‐sep‐MEDI, and χ ‐sep‐iLSQR), along with the differences from the χ‐ sep‐COSMOS maps (as the reference), are displayed in Figure  3 . The difference images are scaled by a factor of five and excluded CSF and vessels. Both χ para and χ dia maps from χ‐ sepnet‐ R 2 ′ and χ‐ sepnet‐ R 2 * reveal less noise and higher correspondence to χ ‐sep‐COSMOS than those from χ‐ sep‐MEDI and χ‐ sep‐iLSQR, illustrating reduced discrepancies in the difference images. The quantitative metrics in Table  1 further consolidate these observations, with χ‐ sepnet‐ R 2 ′ achieving the best outcomes, followed by χ‐ sepnet‐ R 2 * .

χ para and χ dia maps from the five χ‐ separation methods (first row for χ para and third row for χ dia ) and their difference images to reference (second and fourth rows; scaled by a factor of five with CSF and vessels removed). The χ para and χ dia maps from the two proposed pipelines ( χ‐ sepnet‐ R 2 ′ and χ‐ sepnet‐ R 2 * ) display smaller differences compared to those from the other conventional algorithms ( χ‐ sep‐MEDI and χ‐ sep‐iLSQR).

TABLE 1 pSNR, NRMSE, HFEN, and SSIM of the four χ‐ separation algorithms, calculated with respect to χ ‐sep‐COSMOS.
χ‐ sepnet‐ R 2 ′ (proposed)
χ‐ sepnet‐ R 2 * (proposed)
χ‐ sep‐MEDI
χ‐ sep‐iLSQR
χ para
pSNR (dB) 39.1 ± 1.0
a

38.4 ± 0.8
b

35.2 ± 1.2 35.8 ± 1.0 NRMSE (%) 38.1 ± 2.8
a

41.6 ± 2.4
b

60.0 ± 6.8 55.8 ± 5.5 HFEN (%) 47.2 ± 4.2
a

55.3 ± 5.1
b

74.2 ± 8.0 63.9 ± 6.9 SSIM 0.940 ± 0.007
a

0.927 ± 0.007
b

0.868 ± 0.015 0.877 ± 0.014
χ dia
pSNR (dB) 39.5 ± 0.9
a

38.6 ± 0.8
b

35.2 ± 1.2 35.9 ± 1.0 NRMSE (%) 35.7 ± 2.2
a

39.4 ± 1.8
b

60.0 ± 6.8 54.2 ± 5.0 HFEN (%) 48.4 ± 2.7
a

57.7 ± 2.1
b

78.0 ± 4.7 67.1 ± 4.4 SSIM 0.937 ± 0.008
a

0.922 ± 0.008
b

0.854 ± 0.017 0.864 ± 0.016
a
Best metric among the four χ‐ separation algorithms.
b
Second best metric among the four χ‐ separation algorithms.

pSNR, NRMSE, HFEN, and SSIM of the four χ‐ separation algorithms, calculated with respect to χ ‐sep‐COSMOS.

Best metric among the four χ‐ separation algorithms.

Second best metric among the four χ‐ separation algorithms.

When the χ para and χ dia maps are zoomed‐in for details in basal ganglia region, including caudate nucleus (CN), putamen (Put), globus pallidus (GP), pulvinar (Pul), dorsomedial nucleus (DN), posterior limb of internal capsule (PLIC) and anterior limb of internal capsule (ALIC), they reveal hyperintense χ para in iron‐rich regions (CN, Put, GP, and Pul), hyperintense χ dia in white matter regions (PLIC and ALIC), and hypointense χ dia in regions known to have little myelin (CN, Put, GP, Pul, and DN), agreeing with previous histological observations (Naidich et al.  2012 ; Schaltenbrand  1977 ) (Figure  4a ). Figure  4b shows a zoomed‐in axial slice of the χ ‐separation maps in the central sulcus region. The hand knob region clearly demonstrates a higher χ para intensity in the motor cortex (MC) than that in the sensory cortex (SC), which is consistent with the well‐known iron distribution observed in histology (Stuber et al.  2014 ). The cortical boundary appears blurred in the results from χ ‐sep‐MEDI and χ ‐sep‐iLSQR, potentially due to the streaking artifacts. Vessels show up in both χ para and χ dia , which is an artifact in χ ‐separation (Lee, Ji, and Oh  2024 ).

Zoom‐in axial slice of the χ‐ separation methods, displaying basal ganglia (a) and central sulcus (b). (a) Deep gray matter regions (Put: Putamen, CN: Caudate nucleus, GP: Globus pallidus, Pul: Pulvinar, and DN: Dorsomedial nucleus) and white matter regions (ALIC: Anterior limb of internal capsule, and PLIC: Posterior limb of internal capsule) are clearly delineated in the χ para and χ dia maps. (b) In the hand‐knob region, maps from χ‐ sepnet‐ R 2 ′ and χ‐ sepnet‐ R 2 * manifest a higher concentration of positive susceptibility sources (most likely iron) in the motor cortex (MC) than in the sensory cortex (SC). Vessels show up in both χ para and χ dia , which is an artifact in χ ‐separation.

The results of the ROI analysis in the healthy subjects are presented in Figure  5 . The χ‐ sepnet‐ R 2 ′ and χ‐ sepnet‐ R 2 * pipelines exhibit strong quantitative correspondence to the reference ( χ‐ sep‐COSMOS), reporting slopes of 1.043 and 1.051 for χ para , respectively, and 0.939 and 0.923 for χ dia , respectively. The results from χ‐ sep‐MEDI also exhibit good consistency with the reference, showing slopes of 0.938 for χ para and 0.934 for χ dia , whereas χ‐ sep‐iLSQR reveals the lowest correspondence (slopes of 0.835 and 0.813 for χ para and χ dia , respectively). Similarly, the R 2 values demonstrate the best correlations of χ‐ sepnet‐ R 2 ′ with respect to the reference, with R 2 of 0.994 and 0.974 for χ para and χ dia , respectively. χ‐ sepnet‐ R 2 * also shows good correlations with the reference with R 2 of 0.993 for χ para and 0.966 for χ dia . In contrast, χ‐ sep‐MEDI and χ‐ sep‐iLSQR reveal inferior performances, reporting R 2 of 0.989, and 0.984 for χ para , respectively, and 0.955, and 0.914 for χ dia , respectively. All measurements of the 24 ROIs are listed in Section  S2 .

ROI analysis results. Each point in the plot represents the mean value of one of the 26 ROIs in one participant. The linear regression slope and R 2 value are displayed in each panel, and reveal that χ‐ sepnet‐ R 2 ′ has the best performance, closely followed by χ ‐sepnet‐ R 2 * . χ ‐sep‐iLSQR shows the poorest performance.

All the results so far suggest that χ‐ sepnet‐ R 2 ′ provides the best performance among the four methods and is closely followed by χ‐ sepnet‐ R 2 * . On the other hand, the other two methods reveal qualitative and quantitative performance degradation.

The results of the ablation study and the resolution generalization study show that the proposed model with all three inputs performs the best while the model, trained with QSM and R 2 ′ maps, matches the performance (see Sections  S3 and S4 ). On the other hand, the model, trained with local field and R 2 ′ maps, demonstrates poor performance in the resolution generalization study for anisotropic resolution (1 × 1 × 3 mm 3 ) (see Section  S4 ).

χ ‐sepnet‐ R 2 ′ and χ ‐sepnet‐ R 2 * in MS Patients

When χ para and χ dia are reconstructed by χ‐ sepnet‐ R 2 ′ and χ‐ sepnet‐ R 2 * , and their lesion appearance are compared in the MS patients, majority of lesions exhibit similar contrasts. In terms of lesion categorization based on χ para and χ dia contrasts, χ para maps show consistent results in 249 out of the 250, while χ dia maps do in 246 lesions (Tables  2 and 3 ). Most lesions appear isointense or hyperintense compared to surrounding NAWM in χ para (isointense: 137, hyperintense: 96; Table  2 ). Out of the 96 hyperintense lesions in χ para , 26 showed the PRL shape (Table  2 ). χ dia typically shows hypointense or subtle hypointense lesions (hypointense: 159, subtle hypointense: 44; Table  3 ). No lesion exhibits hyperintense or even subtle hyperintense in χ dia .

TABLE 2 Total of 250 MS lesions are assessed for visual characteristics in the χ para maps from the χ‐ sepnet‐ R 2 ′ and χ‐ sepnet‐ R 2 * pipelines. Within the 96 hyperintense χ para lesions, 26 was classified as PRL.
χ para

χ‐ sepnet‐ R 2 *
Hyper (PRL) Subtle hyper Iso Subtle hypo Hypo Total
χ‐ sepnet‐ R 2 ′
Hyper (PRL) 96 (26) 0 0 0 0 96 Subtle hyper 0 1 0 0 0 1 Iso 0 0 137 0 0 137 Subtle hypo 0 0 0 5 1 6 Hypo 0 0 0 0 10 10 Total 96 1 137 5 11 250

Total of 250 MS lesions are assessed for visual characteristics in the χ para maps from the χ‐ sepnet‐ R 2 ′ and χ‐ sepnet‐ R 2 * pipelines. Within the 96 hyperintense χ para lesions, 26 was classified as PRL.

TABLE 3 Total of 250 MS lesions are evaluated for their visual characteristics using the χ dia maps generated by the χ‐ sepnet‐ R 2 ′ and χ‐ sepnet‐ R 2 * pipelines.
χ dia

χ‐ sepnet‐ R 2 *
Hyper Subtle hyper Iso Subtle hypo Hypo Total
χ‐ sepnet‐ R 2 ′
Hyper 0 0 0 0 0 0 Subtle hyper 0 0 0 0 0 0 Iso 0 0 43 0 0 43 Subtle hypo 0 0 1 44 3 48 Hypo 0 0 0 0 159 159 Total 0 0 44 44 162 250

Total of 250 MS lesions are evaluated for their visual characteristics using the χ dia maps generated by the χ‐ sepnet‐ R 2 ′ and χ‐ sepnet‐ R 2 * pipelines.

In Figure  6 , representative MS lesions from each category are visualized in T 2 ‐weighted FLAIR image, and χ para and χ dia maps, with yellow arrows indicating lesion locations. The blue boxes reveal the results with good consistency between χ‐ sepnet‐ R 2 ′ and χ‐ sepnet‐ R 2 * . The red boxes specify the results with the discrepancy, where the lesions appear subtly hypointense in the χ para map of χ‐ sepnet‐ R 2 ′ , but hypointense in χ‐ sepnet‐ R 2 * (Figure  6g ), and subtly hypointense in the χ dia maps of χ‐ sepnet‐ R 2 ′ , but hypointense (Figure  6l ) and isointense (Figure  6m ) in χ‐ sepnet‐ R 2 * .

Sample MS lesion appearances in the χ para and χ dia maps generated by χ‐ sepnet‐ R 2 ′ and χ‐ sepnet‐ R 2 * . The lesion locations are marked by yellow arrows. The maps exhibiting similar visual characteristics are indicated by blue boxes, while those with different characteristics are indicated by red boxes. In χ para , most of the MS lesions exhibit similar contrasts in both maps, reporting (a) hyperintense, (b) PRL, (c) subtle hyperintense, (d) hypointense, (e) subtle hypointense, and (f) isointense. Only one lesion out of 250 reveals slightly different visual characteristics (g). In the χ dia maps, most MS lesions appear as (h) hypointense, (j) subtle hypointense, or (k) isointense. Similar visual patterns were observed in the perilesional rim (i). Only four lesions out of 250 were categorized as different lesion types, two of which are displayed in (l) and (m).

Discussion

In this study, we developed two χ ‐separation reconstruction pipelines, χ‐ sepnet‐ R 2 ′ and χ‐ sepnet‐ R 2 * , that successfully produced COSMOS‐quality susceptibility source‐separated maps in healthy subjects and MS patients. Both pipelines revealed high‐quality images with good quantitative metrics. In particular, χ‐ sepnet‐ R 2 * achieved good image quality and quantitative accuracy even without R 2 measurement, substantially improving the applicability of χ ‐separation in clinical and neuroscientific research.

In this work, multi‐echo data trained QSMnet (QSMnet ME ) is developed to replace the original QSMnet, which utilized single‐echo data. The metrics of QSMnet ME within the brain mask excluding CSF and vessels (pSNR: 33.3  ±  0.9 dB, NRMSE: 48.0%  ±  4.3%, HFEN: 36.8%  ±  2.8%, SSIM: 0.95  ±  0.01; Figure  2b ) and within the entire brain mask (pSNR: 32.5  ±  0.9 dB, NRMSE: 47.1%  ±  3.9%, HFEN: 36.1%  ±  2.6%, SSIM: 0.95  ±  0.01) were slightly better than those of QSMnet (pSNR: 37.4  ±  1.1 dB, NRMSE: 50.0%  ±  4.5%, HFEN: 48.5%  ±  5.9%, SSIM: 0.91  ±  0.02) (Yoon et al.  2018 ) and QSMnet+ (pSNR: 37.3  ±  1.0 dB, NRMSE: 50.7%  ±  4.8%, HFEN: 49.2%  ±  6.1%, SSIM: 0.91  ±  0.02) (Jung et al.  2020 ) except for pSNR. The slight improvement may be from the larger number of head orientations used in the training dataset of QSMnet ME .

R2PRIMEnet was proposed to infer R 2 ′ from R 2 * , which was motivated by the fact that R 2 * contains information related to R 2 ′ . Similar contrast conversion tasks using deep learning have been proposed previously (Hagiwara et al.  2019 ; Moya‐Sáez et al.  2021 ; Tanenbaum et al.  2023 ). Furthermore, such a network was recently approved by a regulatory agency (Tanenbaum et al.  2023 ), supporting the reliability of such networks. Our quantitative metrics (pSNR: 42.7  ±  4.2 dB, NRMSE: 25.3%  ±  1.9%, HFEN: 44.9%  ±  5.2%, SSIM: 0.949  ±  0.006; Figure  2d ), which are better than those of QSMnet ME in pSNR and NRMSE, suggest that estimating R 2 ′ from R 2 * is feasible, although a direct comparison between R2PRIMEnet and QSMnet ME is cautious due to their different contrasts.

When compared to the χ ‐sep‐COSMOS results, maps from χ ‐sepnet based pipelines seem to provide less noisy results (Figures  3 and 4 ). This may be because a network might focus on structural information rather than noise. Additionally, χ‐ sep‐COSMOS may have encountered registration issues when processing data with multiple head orientations and multi‐echo SE data that had thicker slices than the multi‐echo GRE data.

When comparing the metrics of χ‐ sepnet‐ R 2 ′ ( χ para pSNR: 39.1  ±  1.0 dB, NRMSE: 38.1%  ±  2.8%, HFEN: 47.2%  ±  4.2%, SSIM: 0.940  ±  0.007; χ dia pSNR: 39.5  ±  0.9 dB, NRMSE: 35.7%  ±  2.2%, HFEN: 48.4%  ±  2.7%, SSIM: 0.937  ±  0.008) and χ‐ sepnet‐ R 2 * ( χ para pSNR: 38.4  ±  0.8 dB, NRMSE: 41.6%  ±  2.4%, HFEN: 55.3%  ±  5.1%, SSIM: 0.927  ±  0.007; χ dia pSNR: 38.6  ±  0.8 dB, NRMSE: 39.4%  ±  1.8%, HFEN: 57.7%  ±  2.1%, SSIM: 0.922  ±  0.008) to those of QSMnet ME (pSNR: 33.3  ±  0.9 dB, NRMSE: 48.0%  ±  4.3%, HFEN: 36.8%  ±  2.8%, SSIM: 0.950  ±  0.011), they demonstrated improved performances in pSNR and NRMSE, although HFEN and SSIM of χ‐ sepnet‐ R 2 ′ and χ‐ sepnet‐ R 2 * were somewhat inferior to those of QSMnet ME . Overall, these suggested high‐quality of the results of χ‐ sepnet‐ R 2 ′ and χ‐ sepnet‐ R 2 * . Still, χ‐ sepnet and QSMnet ME have different contrasts so direct comparison is not straightforward.

In the MS lesion study, χ‐ sepnet‐ R 2 * produced comparable lesion characteristic results to χ‐ sepnet‐ R 2 ′ , confirming its utility (Tables  2 , 3 , and Figure  6 ). Only five out of 500 lesions (both χ para and χ dia lesions combined) were categorized differently between the two methods. The differences were subtle as demonstrated in Figure  6g,l,m . These outcomes suggest that χ‐ sepnet‐ R 2 * can successfully generate χ ‐separation results without an R 2 map when evaluating χ‐ separation characteristics in MS lesions. Furthermore, there was no lesion exhibiting hyperintense or even subtle hyperintense in χ dia , which is consistent with the fact that demyelination is one of the primary pathophysiology of MS (Kim et al.  2023 ).

One limitation of lesion studies is that no gold‐standard map exists for the patient data. Hence, one has to be cautious in interpreting the results. Additionally, the R 2 maps were obtained using DE‐TSE with no B1 information instead of multi‐echo SE, introducing errors. Lastly, the proposed method only showed the results for MS lesions and further investigation is necessary to demonstate validity of the method in various diseases that present different susceptibility characteristics.

Although the proposed pipelines demonstrate superior performance compared to the other conventional χ ‐separation methods, they are subject to the intrinsic limitations of χ ‐separation: large vessels introduce artifacts due to flow and non‐local effects (Lee, Ji, and Oh  2024 ), and the models ignore susceptibility anisotropy (Lee et al. 2010 ), impeding the accurate reconstruction. For example, central sulcus in χ dia (Figure  4b ) showed hyperintense, likely caused by artifacts from vessels and pia mater. Further efforts are necessary to improve the χ ‐separation model.

In our study, the networks were trained with six orientation data from five subjects, which may seem to be small. However, the networks were trained by 64 × 64 × 64 size patches and the total number of patches for training were 20,160 for χ‐ sepnet. Hence, the number of training data were sufficient to cover diverse structural and susceptibility characteristics. To further validate this point, an experiment was conducted by increasing the training data up to 10 subjects to investigate the impacts of data size on performance (see Section  S5 ). The results showed that performance saturated when the training data reached five subjects (Figure  S7 ), although limitations can still exist for unseen input characteristics.

The cascading of multiple networks, as utilized in our pipelines, can potentially propagate errors. For example, χ ‐sepnet in the χ ‐sepnet‐ R 2 * pipeline has inputs that are outputs of QSMnet ME and R2PRIMEnet. Therefore, errors from QSMnet ME and R2PRIMEnet can be propagated into the reconstruction results of χ ‐sepnet. While the design and evaluation of a single neural network that integrates these multiple networks could address this issue, it would require careful consideration of factors such as network capacity, interpretability, and the availability of sufficient training data. Additionally, QSMnet ME is limited by its tendency to learn a dipole kernel specific to the resolution of the training dataset, which restricts its performance to that particular resolution (Jung, Bollmann, and Lee  2022 ). Consequently, errors arising from resolutions other than the training data may be transmitted to χ ‐sepnet. However, by using a resolution generalization method, this limitation of QSMnet ME can be partially overcome, and χ ‐sepnet can effectively extend its resolution generalization ability as demonstrated in our results (Section  S4 ).

The proposed pipelines will encounter difficulties when processing out‐of‐distribution data from the trained datasets. For instance, in lesions with high susceptibility values, the pipelines may produce maps with underestimated values as demonstrated in QSMnet+ (Jung et al.  2020 ). Nevertheless, translating the scaling augmentation approach of QSMnet+ to R2PRIMEnet is not straightforward because there is no existing model that directly produces R 2 ′ from R 2 * . Hence, further research is necessary to develop a generalized χ ‐sepnet.

This work has been partially presented at the 2021 Annual Meeting of the International Society of Magnetic Resonance in Medicine (Kim et al.  2022 ). It is important to note that χ‐ sepnet mentioned in the abstract differs from the proposed pipelines in terms of the number of datasets used for the training and the inclusion of R2PRIMEnet.

Conclusion

In conclusion, we proposed two deep learning‐powered susceptibility source separation pipelines, χ‐ sepnet‐ R 2 ′ and χ‐ sepnet‐ R 2 * . Leveraging the capabilities of the deep neural networks, χ‐ sepnet‐ R 2 ′ successfully reconstructed COSMOS‐quality paramagnetic and diamagnetic susceptibility source maps using single‐head orientation multi‐echo GRE data and multi‐echo SE data. Furthermore, χ‐ sepnet‐ R 2 * produced χ ‐separation maps with comparable quality using only multi‐echo GRE data. This advancement significantly lowers the burden of data acquisition, thereby enhancing the clinical availability and impact of χ ‐separation.

Conflicts of Interest

The authors declare no conflicts of interest.

Supporting information

Data S1. Supporting Information.