Fitting of Lorentzians
Any z-spectrum 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 pre-saturation, \(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}_{N-1}\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 Levenberg-Marquardt 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 lib-file 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 z-spectra generation and analysis, whereas the second one is based on a pyramidal approach similar to the Image Dowsampling Expedited Adaptive Least-squares (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 down-sampled 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 up-sampled images. The up-sampling 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 z-spectra 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 pixel-wise shifting of the z-spectrum, e.g. with the water-shift 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 pre-fitting to obtain reliable frequency offsets for the individual components. In this initial step, a sum of six Lorentzian’s is fitted to each z-spectrum 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 z-spectra, and the actual fitting analysis starts. This allows easy assignment of the individual components to the known frequencies.
In the actual analysis of the z-spectra, 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 z-spectra 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 z-spectra 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 Rice-distributed noise of varying amplitude. The noise amplitude ranged from 0–20% in the image without pre-saturation, resulting in ideal z-spectra without noise and "more realistic” z-spectra 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
|
NOE-2.6
|
0.05
|
0.05
|
-2.6
|
1.0
|
NOE-1.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 pre-processing (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 pixel-by-pixel 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 water-based phantom consisted of four inner tubes filled with NAD and Cr dissolved in phosphate-buffered 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.
In-vivo study
To test the utility of the developed software, an in-vivo 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 24-channel spine coil was used for the lumbar spine acquisition. CEST was acquired with an in-house 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 (in-vivo), field of view = 80 x 80 mm2 (phantom), 200 x 200 mm2 (in-vivo), 1 pre-saturation pulse (phantom), 40 pre-saturation pulses (in-vivo), B1 = 1.5 µT (phantom), 0.9 mT (in-vivo), pulse duration of PD = 900 ms (phantom), 100 ms (in-vivo), frequency offsets = 63, frequency range = -7.0–7.0 ppm (phantom), -5.0–5.0 ppm (spine).