Fitting of Lorentzians
Any zspectrum can be modeled by the sum of several Lorentzian functions described by the following function:
$$L\left({a}_{k},{\omega }_{k}^{c},{\sigma }_{k}\right)= 1\frac{I}{{I}_{0}}=\sum _{k=1}^{K}\frac{{a}_{k}}{1+4{\left(\frac{\omega {\omega }_{k}^{c}}{{\sigma }_{k}}\right)}^{2}} \left(1\right)$$
where \(I\) is the image intensity, \({I}_{0}\) is the image intensity without presaturation, \(K\) is the number of Lorentzian components, \(\omega\) is the frequency and \({a}_{k}\), \({\omega }_{k}^{c}\) and \({\sigma }_{k}\) are the amplitude, frequency offset and width of the \(k\)th CEST proton pool.
According to probability theory, to obtain maximum likelihood estimates, the following functional \(\)must be minimized:
where \({y}={\left[{y}_{0},{y}_{1},{y}_{2},\dots {y}_{N1}\right]}^{{T}}\) is the signal vector, \(N\) is the number of frequency offsets measured, \({a}={\left[{a}_{0},{a}_{1},{a}_{2},\dots {a}_{K}\right]}^{{T}}\) is the amplitude vector, \({{\omega }}^{{c}}={\left[{\omega }_{1}^{c},{\omega }_{2}^{c},{\omega }_{2}^{c},\dots {\omega }_{K}^{c}\right]}^{{T}}\) is the frequency vector and \({\sigma }={\left[{\sigma }_{1},{\sigma }_{2},{\sigma }_{3},\dots {\sigma }_{K}\right]}^{{T}}\) is the width vector of the Lorentzian components. Ψ is the following \(N x K\) matrix of full rank:
The superscript T denotes the transpose and \(‖\dots ‖\) the Euclidian norm.
The functional Φ consists of the sum of squared residuals and leads to a typical nonlinear least squared problem. To solve such problems, Dennis et al. developed an adaptive NLLS algorithm that can handle large residuals or very nonlinear problems better than a LevenbergMarquardt algorithm [25, 26]. In our implementation, we used the subroutine dn2gb from the port library of netlib [27], an established collection of mathematical software commonly used in science and engineering. By applying upper and lower bounds in dn2gb, prior knowledge such as large linewidths for magnetization transfer components or known frequencies of CEST metabolites can be included in the fitting process. Adding such physical requirements as additional bounds leads to maximum accuracy and robustness [28]. As a secant method, the NLLS algorithm applied here uses the evaluation of the residuals and the Jacobi matrix\(J\), which consists of the first derivatives of the residuals according to the unknown parameters:
$$\frac{{\partial J}_{i,k}}{\partial {a}_{k}}= \frac{1}{1+4{\left(\frac{{\omega }_{i}{\omega }_{k}^{c}}{{\sigma }_{k}}\right)}^{2}}, \frac{{\partial J}_{i,k}}{\partial {\omega }_{k}}= \frac{8{a}_{k}({\omega }_{i}{\omega }_{k}^{c})}{{{\sigma }_{k}^{2}\left[1+4{\left(\frac{{\omega }_{i}{\omega }_{k}^{c}}{{\sigma }_{k}}\right)}^{2}\right]}^{2}}, \frac{{\partial J}_{i,k}}{\partial {\sigma }_{k}}=\frac{8{a}_{k}{({\omega }_{i}{\omega }_{k}^{c})}^{2}}{{{\sigma }_{k}^{3}\left[1+4{\left(\frac{{\omega }_{i}{\omega }_{k}^{c}}{{\sigma }_{k}}\right)}^{2}\right]}^{2}} \left(4\right)$$
Software development
To create an easy to use tool for CEST analysis, the proposed software was developed under Windows 10 operating system in Microsoft Visual C + + 2022 using the MFC library as graphical user interface. The usual design elements of a typical Windows software have been used so that users can easily and quickly find their way around. Since the netlib library was developed in the Fortran programming language, dn2gb was compiled with the Intel® Fortran compiler ifort v.2021.6.0. Then the corresponding libfile was created using Microsoft Library Manager, which is part of Visual Studio. Fortran passes all arguments of subroutines by reference. Therefore, pointers must be used when calling Fortran subroutines from C++. Data structures like vectors, matrices and arrays are represented in different ways in Fortran and C++. Consequently, when calling the external Fortran subroutines from C++, special care must be taken with the construction and order of the variables.
Filtering and pyramidal approach
The SNR of CEST data is a relevant aspect to be considered. For this reason, two approaches were proposed to improve the stability of the results. The first one uses image filtering prior to zspectra generation and analysis, whereas the second one is based on a pyramidal approach similar to the Image Dowsampling Expedited Adaptive Leastsquares (IDEAL) method presented by Zhou et al. [12]. For noise filtering in image space, a non local means (NLM) filter [29] was implemented using the OpenCV 4.5.1 library [30]. Compared to simple Gaussian filters, which cause significant image smoothing, NLM filters preserve image details. In the pyramidal method, image downsampling is used to increase SNR. First, the image size of the original image data is strongly reduced by a factor of 2N, where N is a small integer number. Next, the averaged signal of the downsampled images is analyzed by Lorentzian fitting and the obtained results are then extrapolated by a factor of two and serve as initial values for analyzing the upsampled images. The upsampling steps are repeated until the original image resolution is reached. Zhou et al. demonstrated that this approach leads to more robust results both in phantom and rat brain data [12].
CEST analysis
Due to the inhomogeneity of the static magnetic field B0, the center frequency of the water peak, which is normally at 0 ppm, is often shifted. Depending on the inhomogeneity of B0, the shift can exceed 1 ppm [31]. When analyzing the zspectra by Lorentzian fitting, B0 shift is less relevant as long as the fitted Lorentzian components can be assigned to the individual proton pools. This means that corrections of the B0 inhomogeneity by pixelwise shifting of the zspectrum, e.g. with the watershift referencing (WASSR) method [14], are not necessary for the Lorentzian fitting because no asymmetry analysis is performed. Nevertheless, our CEST analysis performed frequency correction by prefitting to obtain reliable frequency offsets for the individual components. In this initial step, a sum of six Lorentzian’s is fitted to each zspectrum to determine the center of the water resonance. The fitting is loosely constrained and starts with a large peak for water at 0 ppm, two peaks above and below the water peak at ± 1.0 ppm and ± 2.0 ppm and a broad peak at 2.0 ppm representing possible semisolid magnetization transfer (MT) contrast in the spectrum. The resulting frequency of the fitted water peak is then used to shift the original zspectra, and the actual fitting analysis starts. This allows easy assignment of the individual components to the known frequencies.
In the actual analysis of the zspectra, a set of starting values is introduced to the fitting process. This set contains the number of Lorentzian components, the frequency range swept by RF saturation pulses, and three parameters for each component, namely the amplitude, frequency offset, and line width. For the latter three parameters, constraints are given by lower and upper bounds, restricting the possible solutions for each individual Lorentzian line. In this way, prior physical knowledge is incorporated into the fitting process. When analyzing CEST images with the pyramid method, it is possible to select the number of pyramid levels and a tolerance value in percent that specifies a range within which the parameters are allowed to change from one pyramid level to the next.
The fitting process of the zspectra yields the optimal amplitude, frequency and width for each individual Lorentzian line, the integral of each component, and the sum of squares between the fitting result and the original zspectra as a quality measure.
Creation of digital phantoms
To evaluate the newly developed algorithms and software, two digital phantoms with different compositions of proton pools were created. The advantage of a digital phantom is the known Lorentzian composition as a ground truth. The phantoms were overlaid pixel by pixel with Ricedistributed noise of varying amplitude. The noise amplitude ranged from 0–20% in the image without presaturation, resulting in ideal zspectra without noise and "more realistic” zspectra of different quality
A cylindrical structure with two inner circles was chosen as the first phantom, as shown in Fig. 1a. The outer region of the phantom consisted entirely of a water pool, and the left and right inner regions contained a composition of water (0 ppm), creatine (1.9 ppm), amide (3.5 ppm), OH (1.0 ppm) and MT (2.0 ppm) pools. The amide and creatine amplitudes varied between the left and right inner regions, while all other components remained constant. The complete composition is shown in Table 1.
Table 1
Composition of digital phantom #1. Amide and creatine amplitude varied between the left and right regions, while all other pools remained constant. The outer region consisted of water only
Proton pool

Amplitude (I/I0) left

Amplitude (I/I0) right

Amplitude (I/I0) outside

Frequency [ppm]

Width [ppm]

H2O

0.95

0.95

0.95

0.0

0.5

Amide

0.30

0.10

0.00

3.5

0.5

Creatine

0.10

0.30

0.00

1.9

0.5

OH

0.03

0.03

0.00

1.0

0.5

MT

0.03

0.03

0.00

2.0

15.0

In addition, a second digital phantom was created with a structure similar to that of a sagittal spine image in which only the intervertebral discs (IVDs) were preserved (see Fig. 1b). Since this model is typical for CEST studies of glycosaminoglycans (GAG), two regions were defined in each disc, representing the inner nucleus pulposus (NP) and the outer annulus fibrosus (AF). The composition of the two regions was water (0 ppm), GAG (OH) (1.0 ppm), NH (3.2 ppm), MT (2.43 ppm), NOE (2.6 ppm) and NOE (1.0 ppm). The amplitude of the GAG component varied between NP and AF and was greater in NP. The complete composition is shown in table 2.
Table 2
Composition of digital phantom #2. The GAG (OH) amplitude varied between AF and NP, while all other pools remained constant
Proton pool

Amplitude (I/I0) AF

Amplitude (I/I0) NP

Frequency [ppm]

Width [ppm]

H2O

0.75

0.75

0.0

2.0

GAG (OH)

0.03

0.06

1.0

1.0

NOE2.6

0.05

0.05

2.6

1.0

NOE1.0

0.001

0.001

1.0

2.0

NH

0.03

0.03

3.2

0.5

MT

0.15

0.15

2.4

10.0

Validation of the software
To validate our software, the digital phantoms were overlaid with simulated noise and systematically analyzed. Noise with amplitudes from 0 to 20% for Phantom 1 with large geometric structures and from 0 to 10% for Phantom 2 with small geometric structures was superimposed in steps of 1%. Since the image noise severely reduces the visibility of small structures, a smaller amount of noise was added in Phantom 2. The resulting CEST images were analyzed in three ways: without any preprocessing (normal), after applying the NLM filter, and using the pyramidal approach (pyramid). In Phantom 2, a combined method of NLM and pyramidal approach was also used to further improve the CEST analysis. The pixelbypixel Lorentzian fitting resulted in parametric images of amplitudes, frequency offsets and widths. Further, integral images for each individual proton pool as well as the sum of squares maps were created. All the calculated parameter maps were analyzed within the predefined regions with different proton pool compositions. The mean and standard deviation of the regions were analyzed as a function of the added image noise. In addition, the standard deviations of the predefined regions were compared between the three approaches: normal, NLM and pyramid. Besides, the regional sum of squares for the three different methods were compared as a quality measure for the fit analysis.
Phantom
To validate the algorithms and software for analyzing measured MRI data, a phantom containing different amounts of creatine (Cr) and nicotinamide (NAD) was created. This waterbased phantom consisted of four inner tubes filled with NAD and Cr dissolved in phosphatebuffered saline (ROTI®Cell PBS, Carl ROTH) at concentrations ranging from 50 mM to 100 mM (Table 3), and physiological pH of 7.3.
Table 3
Composition of the phantom for real MRI measurements. The phantom consisted of 4 inner tubes with different concentrations of amide (NAD) and creatine (Cr)
Tube

#1

#2

#3

#4

NAD

100 mM

100 mM

50 mM

50 mM

Cr

50 mM

100 mM

100 mM

50 mM

The acquired CEST data from the phantom were preprocessed by NLM filter to increase the SNR before analysis by Lorentzian fitting.
Invivo study
To test the utility of the developed software, an invivo study was performed on a healthy volunteer. CEST data were acquired from the lumbar spine. Based on the results of the digital phantom #2, the CEST data of the spine were processed using the NLM and pyramidal approach. The study was approved by the local ethics committee (Ethics Committee of the Medical Faculty of Heinrich Heine University, Düsseldorf, Germany). Written informed consent was obtained from healthy volunteers who participated this study.
MRI measurements
All acquisitions were performed on a Siemens Magnetom Prisma MRI (Siemens Healthineers, Erlangen, Germany) at a magnetic field strength of 3T. The phantom experiment was performed with an eight channel knee coil. A 24channel spine coil was used for the lumbar spine acquisition. CEST was acquired with an inhouse developed spoiled gradient echo sequence using the following parameters: repetition time TR = 7.2 ms, echo time TE = 3.5 ms, excitation flip angle = 15°, image matrix = 128 x 128, slice thickness = 10 mm (phantom), 6 mm (invivo), field of view = 80 x 80 mm2 (phantom), 200 x 200 mm2 (invivo), 1 presaturation pulse (phantom), 40 presaturation pulses (invivo), B1 = 1.5 µT (phantom), 0.9 mT (invivo), pulse duration of PD = 900 ms (phantom), 100 ms (invivo), frequency offsets = 63, frequency range = 7.0–7.0 ppm (phantom), 5.0–5.0 ppm (spine).