The selection and parameterization of gradients that are concurrently applied to the RF pulse, forming the excitation k-space trajectory [26], play a vital role in achieving an appropriate pulse design [21, 23, 27,28,29,30,31]. To address this, we conduct an offline trajectory optimization with six volunteer data sets prior to the scans and evaluate the resulting trajectory with six additional volunteer data sets (in total: 6 female, 2 diverse, median age 47 years, range 19–80 years). The data sets include information of the sagittal cervical spine. All in vivo MRI measurements are performed in accordance with institutional guidelines and all healthy volunteers provided informed consent prior to the examination. This optimization involves considering corresponding field maps, masks, scanner properties, and other relevant factors. The resulting trajectory parameterization is then utilized for subject-specific online calculations of the SPSP fat saturation RF pulse. The optimization process presented can be used as a method for improving existing parameter values of an already initialized trajectory, but also as a structured initial parameterization when creating a new trajectory. If a parameterized trajectory that is well suited to the demands (e.g. target application, region) already exists, this process can be skipped.

For the purposes of SPSP fat saturation focused on sagittal imaging of the cervical spine cord, a 2D spiral trajectory with six parameters, enabling in-plane frequency modulation through temporally varying radius and velocity in k-space, is employed. In contrast to a standard spiral trajectory [32], this implementation offers more flexibility to cover k-space and exploit scanner limits. In combination with a matching RF pulse, this trajectory aims to compensate for the existing field inhomogeneities and aims to provide uniform fat saturation without exciting water. A more detailed description of the optimization procedure and an illustration of the workflow can be found in the Supporting Information. The final trajectory is visualized in Fig. 1a, where \(x\) and \(y\) of the trajectory definition correspond to readout and phase encoding directions, respectively. The minimum spatial resolution of the trajectory is 3.13 cm for both directions. The inner spiral (which covers the last 4.2 ms) enables spatial variations with a resolution of 16.0 cm or more.

a Excitation k-space trajectory and gradients for the final parameter set. These parameters are optimized with six volunteer data sets (sagittal cervical spine) and ten initial points following MATLABs pattern search algorithm (part of the Global Optimization Toolbox in R2019b) and individual pulse calculations with subsequent Bloch simulations. The shown trajectory is used for all further experiments. b Complex RF pulse shapes for the proposed fat saturation method (SPSP) in standard and accelerated configurations for an exemplary volunteer measurement. The first column shows the RF magnitude and the second column represents the corresponding phase. As sampling compression factors (sc) increase, the RF pulse shape becomes more discretized

The RF pulse design process involves two phases: the configuration phase and the optimization phase. The configuration phase is performed offline and is the prerequisite of the subsequent online optimization phase.

Within the configuration phase, all settings and parameters required for the optimization procedure are determined. This includes defining the target flip angles (FAs) of 110° for fat (− 3.4 ppm) and 0° for water (0.0 ppm) [33], as well as the desired spatial pattern for both frequencies to achieve homogeneous (non-)excitation. The gradient trajectory is optimized beforehand as described in the previous section and the Supporting Information. The total pulse duration is determined by the number of trajectory samples, considering a sample duration of 10 μs. To ensure compliance with hardware long-term requirements, the RF voltage per sample is limited. In the standard SPSP configuration, an upper boundary of 100 iterations is set for the optimization algorithm. However, for two of the four presented SPSP fat saturation configurations, the number of iterations is increased to 500 to compensate for the potential accuracy loss introduced by the applied accelerated RF calculation approach discussed below. To reduce the size of the optimization problem, the resolutions of the field maps are scaled down by a factor of four.

The optimization phase represents the core element of the SPSP RF pulse design, which is based on an interior point solver [34, 35] implemented by Majewski [25]. In principle, this solver enables a simultaneous calculation of the RF pulse shape, sample duration, and gradient shape. For this work, the sampling duration and gradient shape are not determined during online optimization, but only the RF pulse shape. The optimization is based on explicit first and second order derivatives and aims to minimize the deviation between the resulting magnetization vector \(^\) and a given target magnetization vector \(^\) for all voxels \(n\in (1,\dots ,_)\) and frequency bands \(f\in (1,\dots ,F)\) (here: \(F= 2\)) within a single slice using full Bloch simulations. The objective function for the optimization, as expressed in Eq. 1, minimizes the L2-norm of the difference between the magnetization vectors:

Please note that the sampling duration is set to 10 μs and the excitation k-space trajectory is optimized offline before the actual RF pulse calculation. Therefore, the minimization process in Eq. 1 pertains to the optimization of the complex RF shape only. This results in \(}_}\) degrees of freedom (DOFs) [25], where \(}_}\) represents the number of RF pulse time samples to be optimized. For example, a pulse to be optimized with a total duration of 10 ms and sample duration of 10 μs consists of 1000 complex pulse samples and has a maximum of 2000 DOFs, independent from the number of voxels \(_}}\). To extend the slice-specific optimization to a multi-slice application, neighboring slices are automatically compared according to the similarity of their B0 maps to identify slice blocks. Slices \(_\) and \(_\) (where \(j,i\in \left\\), \(j\ne i\) with S the total number of slices) which are considered within one block need to fulfill both of the following conditions based on the voxel wise B0 difference \(\Delta ^__} = \sqrt^_\right)}-^_\right)})^}\) between voxels at the same in-plane position.

$$(1) \, }\left(\left\^__}: x < 50\, }\right\}\right) < 20\, }$$

$$(2) \, |\left\^__}: x > 25\, }\right\}| < 0.3|\Delta ^__}|$$

(2)

with \(|x|\) denoting the cardinality of a given set. Threshold values are found empirically. These conditions ensure that the mean B0 difference for reasonable values (\(< 50 }\)) is below 20 Hz and that less than 30% of all values differ by more than 25 Hz. A common SPSP pulse is optimized for each slice block based on the center slice. The overall scheme of the TSE sequence remains unchanged, only the fat pre-saturation pulses are replaced by SPSP pulses. Neighboring slices within one block share the same SPSP pulse, but are acquired independently. For example, a measurement with 15 slices requires a maximum of 15 individual SPSP pulses. If five neighboring slices each have a similar B0 distribution, only three SPSP pulses need to be calculated. In the applied TSE sequence, these three pulses replace the standard fat saturation pulses of the corresponding slices, i.e. SPSP pulse #1 is applied to slices 1–5, SPSP pulse #2 is applied to slices 6–10, SPSP pulse #3 is applied to slices 11–15. Due to the selection of slice blocks for the pulse calculation, the overall preparation time is reduced. The algorithm is capable of minimizing the magnetization deviation by taking into account, among other aspects, B0 and B1 field inhomogeneities as well as system properties, such as the adjust transmitter voltage or the maximum applicable RF-voltage. The solver presented by Majewski, therefore, offers a well-fitting answer for the demands of a SPSP selective pulse design for online subject-specific fat saturation.

The actual implementation of the SPSP pulse design in an online workflow includes the following steps. After starting the TSE adapted with the proposed method, B1 and B0 mapping sequences are acquired. The results are processed (e.g. actual mapping, interpolation) and stored in binary MATLAB files together with the scanner properties. The actual pulse optimization uses this information and another file with the settings from the configuration phase to calculate the SPSP RF pulse in a MATLAB framework on the scanner computer as described in the optimization phase. The output of this optimization replaces the standard spectral fat pre-saturation pulse slice by slice.

To decrease the calculation time and enhance the clinical acceptance of the proposed fat saturation method, different acceleration techniques for RF calculation are developed. One technique involves treating \(_}}\) (sampling compression, sc) neighboring pulse samples as a single combined sample, resulting in a reduction of degrees of freedom (DOFs) by a factor of \(_}}\). This can be seen as down sampling the RF time axis. The resulting smaller optimization problem reduces the calculation time, with higher values of \(_}}\) providing stronger acceleration. Given that the number of RF pulse samples in our scenarios is 990, a configuration with a small \(_}}\) still offers sufficient DOFs for effective optimization. Another technique for accelerating the RF calculation is integrated into the interior point solver itself. In this approach, the second order derivatives are not explicitly determined but approximated based on the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm [36]. The application of BFGS is expected to require more iterations, but the individual iteration steps are performed faster compared to using the exact Hessians. These two acceleration techniques (RF time axis down sampling and BFGS) are independent of each other and can be applied simultaneously. Due to the trajectory setting, the optimization procedure uses \(}_}\) = 990, 445, 165 for SPSP configurations with different sampling compression factors.

To assess the impact of the proposed accelerated RF calculation approaches on calculation time and fat saturation quality, four different SPSP configurations are considered. The "SPSP standard" configuration does not employ any acceleration techniques and is limited to a maximum of 100 iterations. "SPSP bfgs500 sc1" and "SPSP bfgs500 sc2" introduce moderate acceleration by activating the BFGS approximation, combining one or two neighboring samples (\(_}}\)=1, \(_}}\)=2), and allowing a maximum of 500 iteration steps. The fourth configuration, "SPSP bfgs100 sc6," employs the BFGS algorithm, a sampling compression factor of \(_}}\)=6 and a maximum of 100 iterations. This configuration is expected to reveal the limitations of the proposed acceleration methods.

All measurements are conducted on a 1.5 T whole-body system (MAGNETOM Sola, Siemens Healthcare, Erlangen, Germany) using a Head/Neck Coil with 20 receive channels. 2-point Dixon measurements (later referred to as Dixon) are performed to acquire fat/water images for comparison purposes. For the SPSP fat saturation experiments, B0 and B1 maps are obtained prior to the actual measurements using default settings provided by the vendor (additional time ≈ 30 s). The acquired maps are subsequently automatically masked to eliminate any background influences. A clinical T2-weighted TSE sequence is utilized with the standard fat pre-saturation pulse, referred to as "Gaussian fat sat," and the subject-specific SPSP pulse, which replaces the Gaussian pulse. More sophisticated pre-saturation approaches like SHARP pulses [37, 38] or B1-insensitive WET pulses [39, 40] have been presented in the past but are not available for this work. Therefore, the Gaussian pulse is taken as comparison to the proposed SPSP method. The SPSP pulse optimization serves as “push-button” application that does not involve any transfer to a networked computer, but is triggered automatically (including B0 and B1 acquisition) by starting the sequence on the scanner computer. The primary focus of this study is to assess the feasibility of the proposed method in sagittal cervical spine imaging, where significant changes in B0 along the slice direction are not expected. Consequently, non-selective SPSP pulses calculated based on single slices, as described earlier, are applied. Furthermore, experiments involving sagittal foot/ankle imaging and simulations for coronal and transversal cervical spine imaging are conducted to demonstrate the potential extension of the method to other body parts and image orientations. Additionally, the necessity of B1 information in sagittal cervical spine imaging is investigated. The protocol parameters for all utilized sequences are provided in Table 1. Since the proposed SPSP pulses directly replace the Gaussian pre-saturation pulse, TR remains unchanged.

Seven configurations are employed for the measurements, where SPSP pulses use the previously tailored universal trajectory based on in vivo data: (A) with no fat saturation, (B) with 2-point Dixon, (C) with Gaussian fat saturation, (D) with the proposed SPSP fat saturation standard setting ("SPSP standard"), (E/F) with two moderately accelerated settings ("SPSP bfgs500 sc1" and "SPSP bfgs500 sc2"), and (G) with a strongly accelerated setting ("SPSP bfgs100 sc6").

The phantom and volunteer data are evaluated based on the remaining signal and signal inhomogeneity within the selected regions of interest (ROIs). The metric for the "remaining signal" is determined by calculating the overall mean and standard deviation of the signal intensities from individual measurements. This value is ideally 0 for fat only regions and 1 for water only regions. The "inhomogeneity" is quantified by calculating the mean and standard deviation of the signal variations within the single measurements and ROIs. The lower this value is, the more homogeneous the considered ROI appears.

The phantom used in this study has a foot-like shape and consists of two liquid compartments representing water and fat tissue. The water compartment is filled with a nickel sulfate solution (1.25 g NiSO4 × 6H2O per 1000 g H2O), while the fat compartment contains MARCOL 82 oil mixed with MACROLEX blue dye (0.011 g MACROLEX blue per 1000 ml MARCOL 82 oil).

To simulate different acquisitions, all measurements are repeated three times, with adjustment invalidation and re-shimming performed between each run. For evaluation, fat and water ROIs are determined manually.

We performed cervical spine (3 female, median age: 62 years, range 39–66 years) and foot/ankle (2 female, median age: 54 years, range 39–66 years) acquisitions in five healthy volunteers. All measurements are conducted in accordance with institutional guidelines and the volunteers provided informed consent prior examination. This sample size aligns with the median sample size typically employed in MRI technical development studies [41] and should be an adequate number to explore technical feasibility.

Three ROIs (subcutaneous fat dorsal, spine and subcutaneous fat ventral for spine; subcutaneous fat dorsal, heads of metatarsals and tibia distal for foot/ankle) are defined manually. For quantitative analysis the ratio to an image with no fat saturation reveals the relative residual signal. The imaging protocols and SPSP RF pulse optimization settings used for the phantom acquisitions (configurations (A)–(G)) are replicated to ensure comparability in the measurements.