Lorentzian-Corrected Apparent Exchange-Dependent Relaxation (LAREX) Ω-Plot Analysis—An Adaptation for qCEST in a Multi-Pool System: Comprehensive In Silico, In Situ, and In Vivo Studies

Based on in silico, in situ, and in vivo studies, this study aims to develop a new method for the quantitative chemical exchange saturation transfer (qCEST) technique considering multi-pool systems. To this end, we extended the state-of-the-art apparent exchange-dependent relaxation (AREX) method with a Lorentzian correction (LAREX). We then validated this new method with in situ and in vivo experiments on human intervertebral discs (IVDs) using the Kendall-Tau correlation coefficient. In the in silico experiments, we observed significant deviations of the AREX method as a function of the underlying exchange rate (kba) and fractional concentration (fb) compared to the ground truth due to the influence of other exchange pools. In comparison to AREX, the LAREX-based Ω-plot approach yielded a substantial improvement. In the subsequent in situ and in vivo experiments on human IVDs, no correlation to the histological reference standard or Pfirrmann classification could be found for the fb (in situ: τ = −0.17 p = 0.51; in vivo: τ = 0.13 p = 0.30) and kba (in situ: τ = 0.042 p = 0.87; in vivo: τ = −0.26 p = 0.04) of Glycosaminoglycan (GAG) with AREX. In contrast, the influence of interfering pools could be corrected by LAREX, and a moderate to strong correlation was observed for the fractional concentration of GAG for both in situ (τ = −0.71 p = 0.005) and in vivo (τ = −0.49 p < 0.001) experiments. The study presented here is the first to introduce a new qCEST method that enables qCEST imaging in systems with multiple proton pools.


Introduction
Low back pain is one of the most common health concerns worldwide, with a lifetime prevalence of up to 80% and a huge impact on patients' quality of life and socioeconomic status [1]. Degenerative disc disease (DDD) is believed to be the main cause of low back pain [2]. However, the underlying degenerative processes of DDD are not yet fully understood. While clinical-standard magnetic resonance imaging (MRI) is the most sensitive imaging method to assess DDD by depicting structural damage, early alterations cannot be quantified [3,4]. Consequently, more sensitive biochemical techniques are needed to assess the underlying biochemical changes that develop during the early stages of cartilage degeneration.
Chemical exchange saturation transfer (CEST) imaging has emerged as a promising bio-sensitive MRI technique alongside Na-, T1ρ-, and diffusion tensor imaging [3,[5][6][7]. In musculoskeletal imaging, CEST has been used to assess the biochemical composition of the intervertebral discs (IVDs) [3,4,8]. The early detection of degenerative IVD alterations using CEST imaging may allow for the timely diagnosis of degenerative diseases, the initiation of targeted therapy, and a better understanding of degenerative processes. Glycosaminoglycans are an important component of the intervertebral discs and form the side chains of the complex proteoglycan molecule [9]. They are linear polysaccharides, which are composed of disaccharide units. Each of these disaccharide units has three hydroxyl and one NH group [10]. The intervertebral discs themselves consist of the collagen-fibrous annulus fibrosus and the nucleus pulposus, which is predominantly composed of proteoglycans [11]. They serve to soften shocks as well as to distribute pressure to adjacent vertebral bodies. Glycosaminoglycans are essential for this function [11].
By assessing changes in the water volume signal after selective saturation, CEST imaging is sensitive to the detection of small amounts of labile protons [12,13], and thus provides essential information that can complement conventional morphologic MRI methods. Particularly, the CEST contrast is approximately proportional to the concentration and exchange rate of the observed labile protons [14]. In fact, the measured CEST effects vary with the fractional concentration of labile protons and the exchange rate and depend on other parameters such as radiofrequency (RF) power, B 0 field strength, regional pH value, temperature, and T 1 and T 2 relaxation times [15][16][17].
CEST imaging techniques have further evolved over the last few years due to contributions from various technical fields, including the development of mathematical models and hardware improvement [18,19]. Furthermore, methods for the simultaneous determination of fractional concentrations and exchange rate have been introduced [17,20,21]. Based on these new methods, CEST imaging has developed into quantitative CEST analysis (qCEST) [18]. Among others, Dixon et al. indicated that the CEST effect can be plotted as a linear function of 1/B 1 (Ω-plot) and that, consequently, the proton exchange rate (k ba ) from Pool B to Pool A and the labile proton ratio (f b ) of Pool B can be determined by linear regression [20]. However, this method is limited to paramagnetic CEST agents, which can be considered independent of direct water saturation (so-called spillover effect) and the magnetization transfer (MT) effect due to the large chemical shift to water [22]. For the qCEST imaging of endogenous CEST agents with smaller chemical shifts (<5 ppm), Meissner et al. proposed a novel quantitative CEST-MRI method termed AREX (apparent exchange dependent relaxation)-based Ω-plots [21]. This metric uses the MTR REX calculation introduced by Zaiss et al., eliminates the spillover and MT effects, and extends it with a T 1 relaxation component [23]. However, while this approach is based on the assumption of a two-pool system, most human biochemical systems interact with more than one labile proton pool. Consequently, many studies have extended the conventional asymmetry analysis (MTR asym ) to a multi-Lorentzian analysis [19,[24][25][26], allowing the separate determination of different pools' effects. Assuming that each saturation transfer and semi-solid MT signal can be approximated as a Lorentzian lineshape, these various pool effects are considered and corrected [27]. Recently, the first studies have demonstrated an increased accuracy in detecting both the MTR asym and AREX signals based on multi-pool Lorentzian fits in simulations and animal tumor models [25]. Yet, to our knowledge, no prior study has used the separation of individual exchange processes via Lorentzian analysis of different B 1 values to extend AREX-based Ω-plot analysis. Therefore, the present study aimed to (a) develop a CEST quantification algorithm applicable to multi-pool systems called Lorentzian-corrected apparent exchange-dependent relaxation (LAREX)-based Ω-plot analysis, and (b) validate the presented algorithm by in situ and in vivo experiments for several medically relevant multi-pool systems, such as amide proton transfer (APT, frequency offset ∆s = 3.5 ppm) determination of the white matter in the visual cortex and histological determination of glycosaminoglycan (GAG, ∆s = 1 ppm) concentration in human IVDs. Additionally, in vivo experiments are used to illustrate the transferability of the in situ results to real MR measurements. To this end, we hypothesized that (a) the AREX-based Ω-Plot method cannot be readily transferred from a two-pool system to a multi-pool system, (b) the introduced LAREX method can be used to separate the different pools and achieve comparable accuracy concerning the parameters f b and k ba as AREX in an ideal two-pool system, and (c) that conventional AREX-based evaluation cannot detect degeneration both in situ and in vivo due to the precision of different pools in human IVDs, but that these effects can be separated and subsequently evaluated using our LAREX approach.

In Silico Study
The Bloch-McConnell simulation (please refer to Table 1 for an overview of all the parameters used) over an extended physiological range of fractional concentrations (f b ) and exchange rates (k ba ) demonstrated substantial differences in accuracy and feasibility between AREX-and LAREX-based Ω-plots ( Figure 1). To better illustrate the accuracy of the methods and to contextualize them in the AREX-based Ω-plot in the two-pool system, ∆f b ( Figure 1A) and ∆k ba maps ( Figure 1C) were calculated with a color-coded error range of ∆f b = ±1‰ ( Figure 1B) and ∆k ba = ±300 Hz ( Figure 1D). The AREX approach did not accurately display the full analyzed range of f b and k ba . Our proposed LAREX approach allows for widely sufficient accuracy for both k ba value and fractional concentration f b determination. However, analogous to the AREX approach for two-proton pools, the exchange rate k ba and the fractional concentration f b tend to be underestimated for fast exchange rates. Similar findings were observed for our APT-qCEST models of the white and gray matter of the visual cortex (Appendix A Figures A1 and A2), our amine-qCEST model of ex-vivo blood ( Figure A3), and for our creatine qCEST brain model ( Figure A4). Table 1. qCEST parameters with references for the six-pool IVD numerical simulation with solute concentration (f), the solute-water exchange rate (k ba ), longitudinal relaxation time (T 1 ), transversal relaxation time (T 2 ), and solute resonance frequency offset (∆). Figure A7 also shows chondroitin sulfate as structural formula, showing the 3:1 ratio of hydroxyl to amide protons.

Water [8]
Hydroxyl [  (GAG) concentrations f b and exchange rates k ba , as a function of f b and k ba , determined by apparent exchange-dependent relaxation (AREX) Ω-plots on two-and six-pool systems, LAREX Ω-plots on a six pool system, as well as the ground truth maps with a red box with marked the clinically relevant physiological range of IVD changes. The physiological parameters used are based on the exchange rates and relaxation times of human IVDs at 3 Tesla observed in previous studies [8,25,28], and are listed separately in Table 1. In this case, we only used the LAREX-based evaluation for the six-pool system because the separation of different pools characterizes this approach. The calculations are then based on the two-pool system separated in this way. Corresponding figures for further multi-pool systems investigated, such as APT qCEST imaging in the human brain, are provided in Appendix A.

Discussion
This study demonstrates the feasibility of qCEST using a new LAREX-based Ω-plot approach to study fractional concentration in a multi-pool system. We validated our proposed LAREX-based Ω-plot approach for GAG concentration, in silico, in situ (human IVDs from body donors), and in vivo (using IVDs from volunteers). In addition, we have demonstrated the applicability of our method using in silico experiments for APT-qCEST imaging in the brain, amine-qCEST imaging of human blood, and creatine-qCEST analysis in the human brain, considering additional exchanged proton pools in Appendix A.
The results indicate that our proposed method substantially improves the accuracy of determining the fractional concentration and the exchange rate.
In contrast to morphological MRI, compositional MRI quantifies tissue composition beyond mere morphology and structure. However, human tissues are complex biochemical systems consisting of numerous components. Thus, while tissue composition can be adequately quantified using compositional MRI, alterations of a specific metabolite within the tissue are far more difficult to assess. Consequently, simple analysis methods like classical MTR asym analysis or quantitative AREX-based Ω-plot analysis reach their limits here. In contrast, our introduced LAREX-based Ω-plot approach overcomes these challenges and allows for the separate analysis of the tissue's individual functional groups.
We were able to demonstrate that the presence of additional proton pools hampered the classical qCEST approach based on an AREX-based Ω-plot analysis and showed substantial deviations. Our simulations for IVDs showed that AREX substantially and systematically underestimated the fractional concentrations. This occurs in the physiological range from f B = 0.5-5‰ and k sw = 100-500 Hz for the GAG-OH pool, which is clinically relevant. Comparable results were observed by Zhou et al. [31] in a comparison of ideal two-pool phantoms and in situ IVDs. In their study, systematic deviations between in vitro and in situ experiments as a function of pH values were also observed, which were attributed to the non-analyzed pools. In contrast, the LAREX-based evaluation allows for accuracy comparable to that achieved with an ideal system consisting solely of two pools. In addition, we validated the LAREX method by in situ experiments for the quantitative assessment of APT in white and gray matter, the quantitative assessment of amines in human blood, and the assessment of creatine in the brain (Appendix A). Here, we observed comparable results for IVDs. Thus, in future studies, multi-pool methods can be analyzed equivalently using the LAREX approach to two-pool systems. The simulations used in the Appendix A were performed at a field strength of 7T. With the help of the simulations we performed, the LAREX approach we proposed could be successfully validated at 3T (in silico IVD study) and 7T (studies in the Appendix A).
In the in situ experiments, we observed that the AREX-based Ω-plot approach specifically for the determination of fractional concentration did not yield physiological results, and the results analogous to the T 1 and T 2 measurements showed no correlation to the histological scoring. Using the LAREX-based Ω-plot approach, we observed a strong and significant correlation (τ = −0.71 p = 0.005) of fractional GAG concentration to the histological reference and lowered fractional concentrations compared to the subject collective. Studies of our group demonstrated an age dependence for the MTR asym effect [32]. Since the classical AREX method is based on the MTR asym under different B 1 field strengths, a decrease in the GAG-induced MTR asym effect such as advanced age leads to an increase in the influence of noise and the other pools. In particular, the NOE#1 effect at −1.8 ppm significantly affects the evaluation of GAG OH protons at one ppm. This may lead to non-physiological results for IVD of advanced age (age of body donors for the in situ measurements: 88.6 ± 8.7 years). The LAREX-based approach allows, due to the Lorentzian correction, not only for the correction of other pools but also a smoothing of the Z-spectrum optimized for CEST experiments [33], which makes this approach suitable for low fractional concentrations. However, it should be considered that Lorentzian analyses have a number of 3 × n (n = number of assumed pools) free parameters. Moreover, CEST exchange problems are nonlinear, which results in the optimization function having not only one minimum but several local minima. If the initial values are chosen unfavorably compared to the input parameters, jumps between minima may occur. However, this is primarily relevant for artificially generated data with only one evaluated pixel. For in situ or in vivo measurements, this effect is compensated by the slight differences between the local pixels and the number of acquired pixels.
In the in vivo experiments performed using the LAREX approach, we observed a moderate and significant decrease in both the fractional concentration (τ = −0.49, p ≤ 0.0001) and the exchange rate (τ = −0.48, p ≤ 0.0001) of GAGs depending on the IVDs' degenerative stage. The AREX-based evaluation, on the other hand, showed no significant correlations and determined a non-physiologically increase in GAG concentrations with progressive degeneration. With LAREX, a fractional concentration of 4.96 ± 2.54‰ was observed in IVDs without any degeneration (Pfirrmann grade 1) and 0.78 ± 0.26‰ in IVDs with advanced degeneration (Pfirrmann grade 4). Considering three exchanging OH protons at one GAG [10] and a water concentration of about 80% in IVDs [34], the corresponding concentrations are 150 mM (Pfirrmann grade 1) and 25 mM of GAG concentration (Pfirrmann grade 4). Iatridis et al. determined GAG concentrations of 250 µg ± 134 GAG/mg dry IVD tissue [34], which relates to a 50-150 mM GAG concentration, which is consistent with our results. Moreover, in our in vivo experiments, we observed about two times higher fractional concentrations than in the in situ experiments, which agrees with previous studies [10,35], which showed a decrease in GAG as a function of age and associated progression of IVD degeneration. Furthermore, in previous studies, 66-100 mM GAG concentrations were found in articular cartilage [10,35]. Therefore, the cutoff and baseline values of the Lorentzian analysis we presented could be adopted to examine articular cartilage for alterations in GAG concentration in future studies. Moreover, regional differences in healthy IVDs were observed in our in vivo experiments (Pfirrmann ≤ 3, Figure 4), which might be due to the differences between nucleus pulposus and annulus fibrosus. Previous studies have also observed these regional differences using MTR asym [3,4,32].
Numerous previous studies have demonstrated the potential for deep learning in medical imaging and post-processing [36][37][38]. Among others, Zaiss et al. showed the potential for deep learning in qCEST imaging [39]. Recently, Huang et al. published two fast and accurate ways to generate CEST/AREX contrast maps using DeepCEST and DeepAREX in animal experiments [40]. Based on those studies and our results, neural networks could be developed to eliminate the influence of other pools. Thus, our results suggest a promising new approach to qCEST imaging. Further longitudinal studies could periodically examine subjects or animal models to determine if the early deterioration can be detected.
However, our study has some limitations. First, the long measurement time of 84:35 min (in situ) and 36:15 min (in vivo) for a one-slice sequence must be mentioned. This is foremost because regular CEST experiments require long TR times, thus guaranteeing T 1 relaxation time decay. In addition, qCEST analysis requires an RF saturation time (10 s in our study) to reach a steady-state and at least three different CEST experiments with different B 1 amplitudes to perform the Ω-plots. By using dual-gradient echo strategies [41], simulated multi-slice gradient echo (GRE) sequences [42], phase-shifted multiplanar CEST-fast imaging with steady-state free precession (FISP) sequences [43], ultrafast one-shot acquisition [44], compressed sensing [45], and other imaging techniques, this limitation could be addressed in future studies [46]. Second, the two IVD MR studies we performed (in situ and in vivo) have only limited comparability. As these experiments were conducted at different temperatures (20 • C in situ and 37 • C in vivo), there may be systematic differences between the experiments due to temperature differences, which might affect the accuracy of LAREX and AREX. However, it should be noted that temperature changes do not change the linear dependence of the CEST effect on concentration [13]. Therefore, the increased accuracy of the LAREX approach compared to the AREX approach would not change by adjusting the temperature. Third, the present method relies on the fact that Lorentz functions can describe the various CEST effects. Previous studies have shown that Lorentz functions cannot represent rapidly exchanging metabolites [47]. Therefore, the accuracy of the LAREX method depends on the exchange rate of the metabolites of the multi-pool system, and as a result, LAREX is not necessarily better than AREX for all CEST metabolites [47]. Fourth, the application of the Lorentz fitting not only enables a reduction of the spurious effects of other pools but also has a CEST adaptive denois-ing effect. Contrary to this, the AREX-based application only used the NLM denoising filter [48], which is not optimized for CEST. Fifth, our in situ assessments were only semiquantitative concerning GAG concentrations. Recently, Kubaski et al. demonstrated a rapid, sensitive, and accurate measurement of GAG and all isomers by mass spectrometric detection analysis [49]. Thus, in the subsequent studies, the sensitivity of the LAREX approach could be investigated in more detail than the laborious spectroscopic methods. Sixth, the in situ and in vivo studies we performed were performed at only a single time point. In further studies, detailed reproducibility measurements and longitudinal studies over a long period of time are needed to detect degenerative processes at multiple time points and to investigate the degenerative changes. Seventh, in the in situ and in vivo studies we performed using LAREX, we observed a moderate change in k ba as a function of the Pfirrmann grade. Previous studies have shown that pH changes in the IVD due to degeneration [31,50]. Subsequent in situ or in vitro studies will need to infer changes of k ba as function of pH-values. Eighth, the application of the LAREX approach at 7T and its use outside IVDs were validated in our study only by in silico experiments. Further in situ and in vivo studies, guided by the IVD studies we performed, are needed to validate the promising LAREX approach further and prepare it for clinical applicability.

Study Design
This prospective feasibility study included sequential in silico, in situ, and in vivo CEST MRI assessments and was, thus, conducted in three consecutive steps: (1) implementation and validation of the novel LAREX approach for quantitative multipool CEST-MRI evaluation based on a Lorentzian adaptation and AREX-based Ω-plot analyses using in silico studies (2) evaluation of the newly implemented LAREX approach using in situ experiments and histological referencing of GAG content in human IVDs, and (3) demonstration of the in vivo applicability in eight subjects with different stages of DDD.
Written informed consent was obtained from all body donors as well as all volunteers. The study was approved by the local ethics committee (Ethics Committee of the Medical Faculty of Heinrich Heine University, Düsseldorf, Germany, study numbers: in situ 2021-1528 and in vivo 2019-551).

MR Imaging
All MRI measurements were performed on a clinical 3T MRI scanner (MAGNETOM Prisma, Siemens Healthineers, Erlangen, Germany) with a dedicated 15-channel knee coil (Tx/Rx Knee 15 Flare Coil, Siemens Healthineers) for the in situ experiments or with the integrated 24-channel spine coil (Spine Matrix coil, Siemens Healthineers) for the in vivo experiments.
For morphologic reference imaging, a sagittal T1-weighted (T1w) turbo spin-echo (TSE) sequence for the alignment of the compositional sequences as well as a sagittal T2-weighted (T2w) TSE sequence were used. For bio-sensitive MRI, T 1 and T 2 mapping sequences, as well as CEST sequences, were obtained; for T 1 mapping, an inversion recovery TSE sequence with seven inversion times (TIs: 25-3000 ms) (Table 3), while for T 2 mapping, a spin-echo (SE) sequence with 20 different echo times was obtained (TEs: 9.7-194 ms) (Table 3). Furthermore, seven (in situ) or three (in vivo) CEST sequences with different high-frequency pulse amplitudes (B 1 ) were acquired ( Table 3). The saturation module of our CEST sequence used Gaussian RF pulses. To this end, a total of 64 saturated images were acquired at different saturation frequencies around the water resonance between −5 and 5 ppm, with a reference image at 300 ppm, and we systematically increased the B 1 pulse amplitude from 0.6 to 1.2 µT in 0.1 µT (in situ) or 0.3 µT (in vivo) steps. In addition, the pulse duration (t p = 100 ms), pulse delay (t d = 100 ms), and number of presaturation pulses (n p = 40) were chosen to achieve a steady-state condition. For the B 0 inhomogeneity correction [51], a Water Saturation Shift Referencing (WASSR) sequence was acquired with t p = 25 ms, t d = 25 ms, n p = 1, B 1 = 0.2 µT, 22 dynamics, and frequency offsets between −1 and 1 ppm [52]; otherwise, the same sequence parameters were used for the CEST sequence.

In Silico Study
In accordance with previous studies, we simulated the Z-spectra based on the Bloch-  [13,19,[53][54][55]. We evaluated the accuracy of the AREX-and LAREX-based Ω-plot analyses as a function of fractional GAG concentration f b and hydroxyl proton exchange rate k ba . Previously, Stabinska et al. showed that AREX-based Ω-plot analyses are varyingly accurate as a function of f b and k ba [17]. Therefore, we performed additional ideal simulations for a two-pool system consisting solely of water and the substrate to compare our multi-pool results. The same MR (TR = 2500 ms, TE = 3.5 ms) and presaturation parameters (n p = 40, t p = 100 ms, t d = 100 ms, presaturation pulse shape = Gaussian, ∆ω = 5 ppm) as well as a field strength of 3T were used for the simulations and the subsequent in situ and in vivo reference studies. GAG concentrations were systematically varied from 10 mM to 250 mM (f b = 0.23-8.52‰) in an extended physiological range [19,31,56]. Accordingly, OH proton exchange rates k ba were varied between 50 and 1000 Hz as described in previous studies [31,57]. A water concentration of 80% in IVDs was assumed, as shown by Baldoni et al. [29]; the other parameters used are listed with references in Table 1.
We also performed further in silico analyses, which can be found in Appendix A. In these analyses, we investigated our proposed LAREX-based Ω-plot approach for APT qCEST imaging in white matter (Appendix A Figure A1) and gray matter (Appendix A Figure A2), APT qCEST imaging of blood (Appendix A Figure A3), and creatine qCEST imaging in the human brain (Appendix A Figure A4) to illustrate its clinical applicability to a wide range of clinically relevant multi-pool systems.

In Situ Study
Human lumbar intervertebral disc cadavers: The local Institute of Anatomy I (Heinrich Heine University, Düsseldorf, Germany) provided eleven freshly frozen human IVDs from three body donors for the in situ measurements. At the time of death, the mean age of the body donors was 88.6 ± 8.7 years (range 82-101 years), one was female, and two were male. Before each MRI examination, specimens were thawed and warmed to room temperature for at least 24 hours.
In situ MR imaging: All specimens were positioned centrally in the dedicated knee coil. In addition, mechanical positioning devices such as sandbags were used to fix the specimens in the coil. A baseline morphological T1w sequence was acquired to align the bio-sensitive imaging sequences. Subsequently, the specimens were examined with seven different B 1 field strengths as described above. Consistent with previous studies, a WASSR sequence was used for B 0 inhomogeneity correction [3,8,19].
Histological preparation: Following the MRI examination, the specimens were subjected to a standard histological procedure [58]. For this purpose, the IVDs and adjacent vertebrae were decalcified and fixed in Ossa fixona (Diagonal, Münster, Germany), dehydrated, and embedded in paraffin. IVDs were then cut along the mid-sagittal plane. The initial condition of the discs was semi-quantitatively scored according to Thompson et al. [29]. Thompson scoring. as the gold standard, allows for a full assessment of the degenerative stage of discs (score 1-no degeneration, score 5-fully degenerative IVD). T.J.F., who has 33 years of experience in musculoskeletal histopathology, assessed each tissue sample individually. In addition, IVDs were cut into 5-µm-thick slices and stained with Safranin O to validate the Thompson assessments [59]. A conventional light microscope (Motic Easy Scan Infinity 100, MoticEurope, Barcelona, Spain) and dedicated software (Motic ® Images Devices MoticEurope, Barcelona, Spain) were used for validation; deviating staining compared to the Thompson score would lead to IVD exclusion.

In Vivo Study
The in vivo study was carried out on eight volunteers (mean age: 29 ± 4 years; age range: 24 to 35 years; six male, two female) with various stages of disc degenerations. All subjects were positioned headfirst and supine in the center of the MRI scanner. Morphological sequences (i.e., T1w-and T2w-sequences) were obtained for clinical assessment, such as Pfirrmann scoring of IVDs [30]. To determine the Pfirrmann degree of IVD degeneration, signal intensity on T2-weighted MR images was used to estimate water content with morphologic parameters on a scale of 1 (no degeneration) to 5 (completely degenerated). To reduce scanning time, only three of the CEST sequences used for the in vivo experiments (B 1 = 0.6, 0.9, and 1.2 µT, t p = 100 ms, n p = 40, pulse shape = Gauss) were acquired. For the B 0 inhomogeneity correction, a WASSR sequence was performed according to the in situ measurements as previously described.

LAREX: Lorentzian-Corrected Apparent Exchange-Dependent Relaxation
The Lorentzian-corrected apparent exchange-dependent relaxation Ω-analysis is an extension of the classical AREX approach for systems with multiple pools, using Lorentzian analyses to eliminate the influence of other proton pools ( Figure 5). Analogous to the AREX approach, Z-spectra were first normalized. Subsequently, the different overlapping proton pools affecting the CEST effect are fitted and corrected by Lorentzian analyses ( Figure 5). Analogous to the AREX approach, Z-spectra were first normalized (step 1). Next, the different overlapping proton pools affecting the CEST effect are fitted by Lorentzian analyses (step 2, Figure 5A). Subsequently, 2-pool Z-spectra were calculated, which are composed of a superposition of the calculated water pool and the investigated metabolite pool, and MTR asym curves were calculated (step 3, Figure 5B). Finally, the Z-spectra were evaluated analogously to the classical AREX approach, where the measured Z-spectra were replaced by the corrected Z-spectra (step 4, Figure 5C). Therefore, for the IVDs experiments, the Z spectra were fitted considering the nuclear Overhauser enhancement-NOE #1 (−1.6 ppm), NOE #2 (−3.5 ppm), MT, spillover, OH, and NH effects. Due to the large number of free hyperparameters that arise in a multi-pool system, we used a user-defined error function (Equation (1)). For curve fitting, we used the MATLAB optimization solver fmincon (find minimum of constrained nonlinear multivariable function) with a state-ofthe-art sequential quadratic (SQP) algorithm. [60]. Following this correction mechanism, the Z-spectrum corresponds only to a two-pool system, which can be evaluated analogously to previous studies using AREX [17,22].  Table 3. (A) Plot of an exemplary Z-spectra (blue circles) at B1 = 0.7 µT and the fitted Z-spectra (red line), and the Lorentz fits of each pool (colored dashed lines). (B) The plot of the MTR asym curve for the simulated raw Z-spectra (blue circles, Raw MTR asym -Curve), the MTR asym -Curves based on the Lorentz corrected Z-spectra to a two-pool system of water and Gag-OH (red line, Corrected MTR asym -Curve), and the ground truth MTR asym -Curve with based on a Bloch-McConnell simulation with only 2-pools (yellow circles, MTRasym -Curve 2-Pool-System). (C) Visualization of the AREX (blue line) or LAREX-based (red line) Ω-plots. The gray dashed lines illustrate the 95% confidence interval of the linear omega plots.
The error function we proposed is based on the arctangent, which had proven itself as an error function under the minimization of numerous hyperparameters in machine learning [61].
where y corresponds to the simulated or measured values of the Z-spectra, f Lorentzian corresponds to the Lorentzian function, and → X corresponds to the hyperparameters to be optimized. Compared to the commonly used Root-Mean-Squared-Error, the error function differs by an adaptive gradient around zero (Appendix B Figure A5). Since the Z-spectra are normalized to y-values between 0 and 1, just the range smaller than 1 is decisive for the behavior of the error functions. Table 4 shows the fitting parameters we used. For the Lorentzian analysis of the Z spectrum, the different pools are not fitted to a fixed frequency offset. Instead, they are determined in a frequency range (Table 4); these deviations from the presumed frequency offset were corrected where necessary by replacing the frequency in the Lorentzian function with the assumed frequency. Consequently, no deviations from the assumed and evaluated frequency offsets are present in the corrected Z spectra of the LAREX evaluation. The subsequent calculations of the MTR Rex maps, AREX maps, and Ω-plots were consistently performed. Table 4. Starting points and boundaries of the amplitude (A), width (W), and offset (∆) of the six-pool Lorentzian fit, which we used for our in silico, in situ, and in vivo studies of IVDs. Furthermore, we used a custom loss function (Equation (1)

MR Image Analysis
For further data analysis, the acquired in situ and in vivo MR data were segmented independently by two experienced radiologists, L.M.W. (5 years of musculoskeletal imaging experience) and D.B.A. (6 years of musculoskeletal imaging experience), using ITK Snap software (v3.8.0, Cognitica, Philadelphia, PA, USA) [62]. L.M.W. segmented the ROIs twice (i.e., six weeks apart) to assess the intra-reader reliability of the techniques studied, whereas D.B.A. segmented the data only once to determine inter-reader reliability. In addition, D.B.A. assessed the Pfirrmann grade for the in vivo measurements [30]. Additionally, tissue properties were further analyzed using in-house developed MATLAB scripts. T 1 maps were created by fitting the IR measurement data with the non-linear least-squares method as a function of inversion delay (TI): S(TI)~S 0 (1 − 2 × exp(−TI/T1)), where S 0 is the equilibrium signal. T 2 maps were consistently calculated pixel-wise as a function of echo time (TE): S(TE)~S 0 × exp(−TE/T2). Based on the WASSR measurements, the pixel-wise B 1 inhomogeneity offset maps were calculated. The CEST data were normalized to the signal of the first acquired frame with a frequency offset of 300 ppm. Using a custom MATLAB script validated in previous in silico and in vitro studies, the subsequent Ω-plots were performed [17]. To this end, the MR images were first denoised using the non-local means (NLM) filter [48], and Z-spectra were frequency corrected using the offset maps; the MTR Rex maps with inverse asymmetry; and the AREX maps determined with the frequency shift of 1 ppm specific to hydroxyl protons. Subsequently, the relaxation-compensated Ω-plot analysis was calculated. In addition, an R 2b of 100 Hz was assumed for the hydroxyl protons bounded on GAG in accordance with Singh et al. [28].
Due to the experimental design of this study, the significance level was set from p ≤ 0.05 to an adjusted p ≤ 0.00625 according to the conservative alpha adjustment method Bonferroni [66]. This "low" significance level prevented alpha error inflation while simultaneously maintaining statistical power.

Conclusions
In our study of qCEST imaging in multi-pool models, we validated our proposed LAREX method using in silico experiments for a wide range of medically relevant systems such as human IVDs, the white and gray matter of the visual cortex, and human exvivo blood. Moreover, we have therefore shown that Lorentz analyses can be used to extract 2-pool spectra from multi-pool spectra. These calculated 2-pool spectra yield results comparable to real 2-pool spectra in an omega plot analysis. Furthermore, we could transfer the results to in situ and in vivo studies. Therefore, for the first time, it is possible to quantitatively investigate multi-pool systems by means of fractional concentration and exchange rate and by using qCEST with the LAREX approach.  Parameter maps (A,C) and color-coded error maps (B,D) for fractional APT concentrations f b and exchange rates k ba as a function of f b and k ba . The study was performed using a 4-pool system and a field strength of 7 Tesla. The parameters of the exchanging proton pools were: Water (T 1 = 1200 ms, T 2 = 40 ms, and ∆ = 0 ppm); APT (T 1 = 1000 ms, T 2 = 10 ms, f = variable between 0.001-0.2, k ba = variable between 1-300 Hz, and ∆ = 3.5); NOE (T 1 = 1000 ms, T 2 = 0.3 ms, f s = 0.06, k ba = 10 Hz, and ∆ s = −3.5 ppm); MT (T 1 = 1000 ms, T 2 = 10 µs, f s = 0.11, k ba = 50 Hz, and ∆ s = 0 ppm), and were based on the results of Mougin et al. [67]. Furthermore, analogous to the in silico, in situ, and in vivo experiments, the same MR parameters (TR = 2500 ms, TE = 3.5 ms) and presaturation parameters (n p = 40, t p = 100 ms, t d = 100 ms, presaturation pulse shape = Gaussian, ∆ω = 5 ppm) were used. Figure A2. In silico Appendix Study 2-APT qCEST analyses of gray matter of the visual cortex. Parameter maps (A,C) and color-coded error maps (B,D) for fractional APT concentrations f b and exchange rates k ba as a function of f b and k ba . The study was performed using a 4-pool system and a field strength of 7 Tesla. The parameters of the exchanging proton pools were: Water (T 1 = 2000 ms, T 2 = 55 ms, and ∆ = 0 ppm); APT (T 1 = 1000 ms, T 2 = 10 ms, f = variable between 0.001-0.01, k ba = variable between 0-300 Hz, and ∆ = 3.5); NOE (T 1 = 1000 ms, T 2 = 0.3 ms, f = 0.06, k ba = 10 Hz, and ∆ s = −3.5 ppm); MT (T 1 = 1000 ms, T 2 = 10 µs, f = 0.11, k ba = 50 Hz, and ∆ = 0 ppm), and were based on the results of Mougin et al. [67]. Furthermore, analogous to the in silico, in situ, and in vivo experiments, the same MR parameters (TR = 2500 ms, TE = 3.5 ms) and presaturation parameters (n p = 40, t p = 100 ms, t d = 100 ms, presaturation pulse shape = Gaussian, ∆ω = 5 ppm) were used. Figure A3. In silico Appendix Study 3-Amine qCEST analyses of the human ex-vivo blood. Parameter maps (A,C) and color-coded error maps (B,D) for fractional Amine concentrations f b and APT exchange rates k ba as a function of f b and k ba . The study was performed using a 6-pool system and a field strength of 7 Tesla. The parameters of the exchanging proton pools were: Water (T 1 = 2249 ms, T 2 = 48.5 ms, and ∆ s = 0 ppm); Amine (T 1 = 1000 ms, T 2 = 10 ms, f = variable 0.0001-0.01, k ba = variable 50-750 Hz and ∆ = 2.2); APT (T 1 = 1000 ms, T 2 = 0.5 ms, f = 0.0078, k ba = 338.7 and ∆ s = 3.5); NOE #1 (T 1 = 1000 ms, T 2 = 1.2 ms, f = 0.018, k ba = 10 Hz and ∆ = −1.7 ppm); NOE #2 (T 1 = 1000 ms, T 2 = 0.5 ms, f = 0.0144, k ba = 5 Hz and ∆ s = −3.5); MT (T 1 = 1000 ms, T 2 = 0.01 ms, f = 0.0105, k ba = 15 and ∆ = −2.34), and were based on the results of Shah et al. [68]. Furthermore, analogous to the in silico, in situ, and in vivo experiments, the same MR parameters (TR = 2500 ms, TE = 3.5 ms) and presaturation parameters (n p = 40, t p = 100 ms, t d = 100 ms, presaturation pulse shape = Gaussian, ∆ω = 5 ppm) were used. Figure A4. In silico Appendix Study 4-Creatine qCEST analyses of human brain. Parameter maps (A,C) and color-coded error maps (B,D) for fractional Creatine concentrations f b and Creatine exchange rates k ba as a function of f b and k ba . The study was performed using a 6-pool system and a field strength of 7 Tesla. The parameters of the exchanging proton pools were: Water (T 1 = 2000 ms, T 2 = 115 ms, and ∆ = 0 ppm); Creatine (T 1 = 1000 ms, T 2 = 10 ms, f = variable 0.00001-0.001, k sw = variable 50-1000 Hz and ∆ = 1.8); Amide (T 1 = 1000 ms, T 2 = 33 ms, f = 0.0008, k ba = 30 and ∆ = 3.5); Glutamate (T 1 = 1000 ms, T 2 = 10 ms, f s = 0.000125, k sw = 2000 and ∆ s = 3.0); Myo-inositol (T 1 = 1000 ms, T 2 = 10 ms, f = 0.0007, k ba = 600 and ∆ = 0.6); NOE (T 1 = 1000 ms, T 2 = 0.01 ms, f = 0.0795, k ba = 10 and ∆ = −2.4), and were based on the results of Sing et al. [69]. Furthermore, analogous to the in silico, in situ, and in vivo experiments, the same MR parameters (TR = 2500 ms, TE = 3.5 ms) and presaturation parameters (n p = 40, t p = 100 ms, t d = 100 ms, presaturation pulse shape = Gaussian, ∆ω = 5 ppm) were used.