Norbert W Lutz
email: norbert.lutz@univ-amu.fr

Monique Bernard

Multiparametric quantification of heterogeneity of metal ion concentrations, as demonstrated for [Mg 2+ ] by way of 31 P MRS

Keywords: magnesium quantitation, tissue heterogeneity, 31 P MRS, lineshape analysis, heterogeneity descriptors, statistical distributions, histograms, free intracellular Mg 2+

Purpose: To design an approach for quantifying tissue heterogeneity in regard to concentrations of free intracellular magnesium(II), by way of multiple 31 P MRS-derived descriptors characterizing the statistical distribution of [Mg 2+ ] values. Theory: The basic paradigm underlying our novel approach exploits the fact that, for a tissue sample that is heterogeneous with respect to intracellular [Mg 2+ ], the lineshape of the ATP-β 31 P MRS resonance reflects the statistical distribution of [Mg 2+ ] values within the observed sample volume. Appropriate lineshape analysis reveals multiple quantitative statistical parameters (descriptors) characterizing the [Mg 2+ ] distribution. Methods: First, the ATP-β 31 P MRS resonance is transformed into a [Mg 2+ ] curve. Then, the digital points of this [Mg 2+ ] profile are used to construct a histogram using our specially developed algorithms. From this histogram, at least eight quantitative parameters describing the underlying statistical [Mg 2+ ] distribution are computed: weighted mean and median, standard deviation, range, mode(s), kurtosis, skewness, and entropy. Results: Comprehensive evaluation based on in-silico and experimental models demonstrates the validity of this new method, which is particularly sensitive and accurate for [Mg 2+ ] gradients > 10-20 µM. Conclusion: This basic feasibility study should open new avenues for future in-vivo studies in physiology and medicine, notably on conditions involving hypo or hypermagnesia.

INTRODUCTION

Magnesium(II) is the second most abundant intracellular cation in mammals after potassium (K + ) [START_REF] Swaminathan | Magnesium metabolism and its disorders[END_REF]. A number of magnesium deficiency disorders, among these congenital disorders, are known, which are frequently due to decreased gastrointestinal or renal magnesium(II) absorption. In addition, diabetes and endocrine causes may result in hypomagnesemia, while magnesium(II) decline has been shown to play a critical role in the development of irreversible tissue damage following direct or indirect neurotrauma [START_REF] Vink | Regulation of intracellular free magnesium in central nervous system injury[END_REF]. Moreover, drugs including antibiotics, anticancer chemotherapeutics and immunosuppressive agents have been linked to magnesium(II) depletion [START_REF] Swaminathan | Magnesium metabolism and its disorders[END_REF]. In such contexts, determining Mg 2+ levels present in body tissue (rather than blood) allows a direct assessment of variability of intracellular tissue Mg 2+ concentrations.

More than 30 years ago, a non-destructive and non-invasive method to measure intracellular free Mg 2+ levels by 31 P MRS has been introduced [START_REF] Gupta | Rose On the noninvasive measurement of intracellular free magnesium by 31P NMR spectroscopy[END_REF]. (Note that in the following, we exclusively deal with the determination of the concentration of free magnesium(II) cations; these are in fast to intermediate exchange with ATP.) This technique uses naturally occurring, endogenous ATP as an intrinsic Mg 2+ reporter molecule, thus rendering any administration of an exogenous Mg 2+ probe superfluous. Predominantly based on Gupta et al.'s seminal work on Mg 2+ quantitation in human blood and cancer cells [START_REF] Gupta | The determination of the free magnesium level in the human red blood cell by 31P NMR[END_REF][START_REF] Gupta | Noninvasive 31P NMR probes of free Mg2+, MgATP, and MgADP in intact Ehrlich ascites tumor cells[END_REF], early applications also included Mg 2+ measurements in intact frog skeletal muscle [START_REF] Gupta | 31P NMR studies of intracellular free Mg2+ in intact frog skeletal muscle[END_REF] and perfused heart muscle [START_REF] Gupta | Rose On the noninvasive measurement of intracellular free magnesium by 31P NMR spectroscopy[END_REF][START_REF] Borchgrevink | Magnesium and reperfusion of ischemic rat heart as assessed by 31P-NMR[END_REF], pointing to the in-vivo potential of this 31 P MRS technique. Since the 1990ies, this method has been applied to determine Mg 2+ levels in vivo in experimental animals and humans, notably in the brain [START_REF] Lodi | Deficient energy metabolism is associated with low free magnesium in the brains of patients with migraine and cluster headache[END_REF][START_REF] Taylor | Free magnesium levels in normal human brain and brain tumors: 31P chemical-shift imaging measurements at 1.5 T[END_REF][START_REF] Halvorson | Assessment of magnesium concentrations by 31P NMR in vivo[END_REF][START_REF] Resnick | Cellular-free magnesium depletion in brain and muscle of normal and preeclamptic pregnancy: a nuclear magnetic resonance spectroscopic study[END_REF][START_REF] Wary | Effect of 13. Iotti S, Malucelli E. In vivo assessment of Mg2+ in human brain and skeletal muscle by 31P-MRS[END_REF](13) and in skeletal muscle [START_REF] Resnick | Cellular-free magnesium depletion in brain and muscle of normal and preeclamptic pregnancy: a nuclear magnetic resonance spectroscopic study[END_REF][START_REF] Wary | Effect of 13. Iotti S, Malucelli E. In vivo assessment of Mg2+ in human brain and skeletal muscle by 31P-MRS[END_REF](13), but also in tumors [START_REF] Taylor | Free magnesium levels in normal human brain and brain tumors: 31P chemical-shift imaging measurements at 1.5 T[END_REF][START_REF] Stubbs | Metabolic consequences of a reversed pH gradient in rat tumors[END_REF] (an overview of Mg 2+ measurements in these and other tissues has been presented elsewhere [START_REF] Günther | Concentration, compartmentation and metabolic function of intracellular free Mg2+[END_REF]). In all these studies, one [Mg 2+ ] value has been determined for the tissue volume (or each voxel) examined, based on the chemical shift of the maximum of the ATP-β (or NTP-β) resonance of the measured 31 P MR spectrum. However, this approach does not reflect intra-volume (or intra-voxel) tissue heterogeneity with respect to intracellular Mg 2+ levels, as it only provides a single value that may or may not be close to the mean. This is also true for 31 P MRSI, where spatial resolution is rather coarse [START_REF] Barker | Magnesium and pH imaging of the human brain at 3.0 Tesla[END_REF]. In fact, experimental evidence of heterogeneous distribution of intracellular Mg 2+ (subcellular microheterogeneity) has been provided by fluorescence microscopy for a number of different cells [START_REF] Günther | Concentration, compartmentation and metabolic function of intracellular free Mg2+[END_REF]. In addition, intracellular [Mg 2+ ] is known to depend on malignancy and growth [START_REF] Günther | Concentration, compartmentation and metabolic function of intracellular free Mg2+[END_REF], and the latter two characteristics often vary considerably across a given tumor volume [START_REF] Sun | Intra-tumor heterogeneity of cancer cells and its implications for cancer treatment[END_REF]. ] ranges. We assess the validity of our method based on (i) simulated ATP-β 31 P MR spectral lines, and (ii) measured ATP-β 31 P MR spectra of well-defined aqueous test solutions. The rationale for our simulations has been described previously in great detail for an analogous application to thermal heterogeneity [START_REF] Lutz | Thermal heterogeneity within aqueous materials quantified by 1H NMR spectroscopy: Multiparametric validation in silico and in vitro[END_REF][START_REF] Lutz | Multiparametric quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: paradigms and algorithms[END_REF]. Briefly, the main purpose of our in-silico modeling is to examine questions such as: distribution descriptors obtained from simulated spectra compare with those from measured MR spectra of model solutions? Note that empirical rather than ab-initio simulations are employed; both chemical shift and linewidth values are taken from actual 31 P MRS experiments. This bestows high practical relevance upon our calculated models.

In addition to simulations, experiments performed on a series of specially designed phantoms providing well-defined [Mg 2+ ] gradients are presented, along with a preliminary example of an actual in-vivo application. In-vitro and in-silico validation results are presented for field strengths commonly used in present-day animal experiments (9.4 and 11.7 T); similar field strengths (e.g., 7 T) are increasingly used in human studies, which renders our validations relevant for future applications in both animal models and patients. All calculations are performed employing a preprogrammed EXCEL spreadsheet provided as Supporting Information. Our findings are expected to have significant implications for studying the relationship between [Mg 2+ ] variations and other biochemical, biophysical, and physiological properties in a variety of biological tissues in vivo. With future projects aiming at a comprehensive in-vivo validation being planned, a broad scope of research fields in biochemistry, pharmaceutical and medical chemistry may benefit from our novel approach, with potential applications to cardiology, myology, neurology and oncology.

THEORY

The use of 31

Mg !! = ! ! !"#$% ! ! -1 [1] with ! = [!"#] ! [!"#] !"! = ! !"#!! ! ! !"#!! !"#$% ! !"#!! !"# ! ! !"#!! !"#$% [2]
where ! ! !"#$% is the dissociation constant of the MgATP complex; [ATP] f , concentration of free (uncomplexed) ATP; [ATP] tot , total concentration of ATP (free and complexed with Mg 2+ );

! !"#!! , measured ATP-β chemical shift; ! !"#!! !"# , ATP-β chemical shift for uncomplexed ATP; ! !"#!! !"#$% , ATP-β chemical shift for ATP complexed with Mg 2+ ). ! ! !"#$% , ! !"#!! !"# and ! !"#!! !"#$% are empirical
values that are determined by way of calibration measurements; consequently, ! ! !"#$% is also denoted as K app (apparent K d ). These values may vary slightly as a function of intracellular parameters such as pH and ionic strength. Besides the basic relationship involving three parameters obtainable through appropriate calibration (eqs. 1 and 2), more complex relationships have been developed for specific tissue types, e.g., for brain and exercising muscle. Calibrations for these tissues resulted in five and ten fitted parameters, respectively (13). The dedicated EXCEL template provided in Supporting Information currently provides the option to choose between these three and five-parameter algorithms for ppm-to-µM Mg 2+

conversions.

The method for quantitative characterization of distribution curves with multiple modes (= maxima), the positions ([Mg 2+ ]) and relative intensities of these modes can be obtained. Additional details of our algorithms are provided as Supporting

Information, together with all equations used.

METHODS

In this section, the in-silico and in-vitro procedures employed in the validation of our new approach are described; additional technical details are presented in Supporting Information.

In-silico calculations

The purpose of our in-silico calculations was twofold: (i) to test the behavior of statistical Further specifics on the rationale for the choice of our in-silico parameters are given in Supporting Information, in conjunction with details concerning the statistical background. therefore, no new animal experiments were performed (see [START_REF] Lutz | Quantitative in-vivo characterization of intracellular and extracellular pH profiles in heterogeneous tumors: a novel method enabling multiparametric pH analysis[END_REF] for the approval of our animal care and experimental protocol).

RESULTS

The most relevant data of our study are presented in this section. Additional results are included in Supporting Information that also comprises Figures S1 through S9 and Tables S1 through S6. quantitation, the accuracy of the ATP/Mg 2+ method depends on the accuracy of the dissociation constant of the Mg 2+ complex [START_REF] Gupta | Rose On the noninvasive measurement of intracellular free magnesium by 31P NMR spectroscopy[END_REF][START_REF] Günther | Concentration, compartmentation and metabolic function of intracellular free Mg2+[END_REF][START_REF] Mottet | Experimental 31P NMR study of the influence of ionic strength on the apparent dissociation constant of MgATP[END_REF][START_REF] Grubbs | Intracellular magnesium and magnesium buffering[END_REF][START_REF] Trapani | Intracellular magnesium detection: imaging a brighter future[END_REF]. This is equally true for the classical ATP/Mg 2+ method and for our new approach introduced in this work. Although other phosphates with different dissociation constants may be developed that cover higher [Mg 2+ ] values with more precision, ATP has the huge advantage of being a native intracellular compound that renders injection of any extrinsic reporter molecules superfluous.

Studies of simulated multimodal [Mg

Subsequently, we evaluated 31 P MR spectra obtained from two combined coaxial NMR tubes, herein called phantoms. Each phantom contained two separate compartments with solutions of distinct, physiologically relevant Mg 2+ concentrations. Phantoms were designed to generate two modes of unequal heights in 31 P MR spectra (Fig. 5 bottom traces) and, consequently, in the resulting [Mg 2+ ] distribution curves. The relative areas under the modes correctly reflected the expected contributions of the respective samples volumes to the overall spectrum: since in our phantoms the outer coaxial volume measured was three times as large as the inner coaxial volume, the outer solution contributed three times as much signal to the spectrum than the inner solution (with [ATP] being identical in both solutions; see footnote to Table 1). This is shown in Figure 5 for the cases where (i) the 0.05-mM Mg 2+ solution was contained in the outer, and the 0.12-mM Mg 2+ solution in the inner compartment, and (ii) the 0.12-mM Mg 2+ solution was contained in the outer, and the 0.05-mM Mg 2+ solution in the inner compartment. While the spectra obtained from the phantoms show moderate but unambiguous separation of the two underlying modes, the separation is better visible if the two modes are of the same size, as shown in the combined spectrum resulting from a superimposition of two separately acquired spectra (Fig. 5 top traces). A detailed discussion of deconvolution effects on statistical descriptors of experimentally obtained [ATP] heterogeneity is presented in Supporting Information.

DISCUSSION

Magnesium(II) plays a very versatile role in body physiology because of its ability to form reversible bonds with a large number of chelating anions, proteins and membranes.

Consequently, magnesium(II) is an important cofactor involved in hundreds of enzymatic reactions, and is essential for the synthesis of nucleic acids and proteins [START_REF] Emila | Role of magnesium in health and disease[END_REF]. Its best known and perhaps most universal function relates to energy metabolism, where magnesium(II) is known to bind to ligands such as ATP and ATP-requiring enzymes, and to cause conformational changes in Na + ,K + -ATPase [START_REF] Swaminathan | Magnesium metabolism and its disorders[END_REF][START_REF] Lodi | Deficient energy metabolism is associated with low free magnesium in the brains of patients with migraine and cluster headache[END_REF]. The mechanisms behind the numerous known magnesium deficiency disorders as are complex, but are critically dependent on the availability of free bioactive magnesium. Nonetheless, many studies have up until now focused on measurements of total magnesium rather than its free form, making correlations to disease states difficult [START_REF] Glasdam | The Importance of Magnesium in the Human Body: A Systematic Literature Review[END_REF].

Ion-selective electrodes were used to measure intracellular free Mg 2+ levels in isolated ferret ventricular muscle [START_REF] Buri | Intracellular free magnesium and its regulation, studied in isolated ferret ventricular muscle with ion-selective microelectrodes[END_REF] and in other in-vitro settings. However, this method is invasive because it is based on tissue biopsies, and therefore of very limited use in patients. In other ex-vivo samples (usually isolated blood cells) and in experimental applications, fluorescence probes have been used to determine intracellular free Mg 2+ concentrations [START_REF] Arnaud | Update on the assessment of magnesium status[END_REF]. Based on this method, free Mg 2+ concentrations in lymphocytes were found to be decreased in migraine patients vs.

normal controls, while there was no difference in total magnesium(II) levels [START_REF] Thomas | Free and total magnesium in lymphocytes of migraine patients -effect of magnesium-rich mineral water intake[END_REF]. However, it is impossible to perform in-vivo measurements based on fluorescence probes because of limited tissue penetration due to light absorption and scattering, except for some of the most superficial tissues, and in the special case of near infrared fluorescence that requires dedicated dyes (up to 10 cm light penetration) [START_REF] Bai | Recent advances in receptor-targeted fluorescent probes for in vivo cancer imaging[END_REF].

As mentioned in the Introduction section, currently existing noninvasive in-vivo based on ATP solutions containing varying amounts of added MgCl 2 . Experimental conditions as described in footnote to Table 1.

  Necessarily, both macroscopic and microscopic intracellular-[Mg 2+ ] heterogeneity within a given tissue volume measured by31 P MRS is reflected in the overall lineshape of the resulting ATP-β resonance. Therefore, information on the underlying Mg 2+ heterogeneity can be extracted by an appropriate analysis of the ATP-β lineshape. By contrast, the classical method of exclusively evaluating the position of the lineshape maximum disregards most of the information contained in the overall resonance.We present here the first technique designed to quantitatively characterize the veryheterogeneity of tissue [Mg 2+ ].This method is based on multiparametric statistical analysis of the resonance of the ATP-β phosphate moiety. The new paradigm underlying our technique is grounded in the fact that (i) the chemical shift of the ATP-β 31 P MR resonance, δ ATP-β , is a function of the Mg 2+ concentration (based on the dissociation constant, K d , of the ATP-Mg 2+ complex (18,19)), and (ii) the overall ATP-β MRS lineshape reflects the statistical distribution of [Mg 2+ ] values in the measured sample volume or voxel, within certain limits. To analyze [Mg 2+ ] heterogeneity with this method, ATP-β 31 P MRS resonances are first converted to [Mg 2+ ] distribution curves. Then, these [Mg 2+ ] profiles are transformed into histograms by an appropriate algorithm. Based on these histograms, further mathematical algorithms developed by us are used to calculate at least eight quantitative parameters (descriptors) defining the underlying [Mg 2+ ] distribution, i.e., the [Mg 2+ ] heterogeneity within the measured sample volume. We exploit, for the first time, the entire [Mg 2+ ] curve (corresponding to an MRS lineshape) to quantitatively describe a distribution of [Mg 2+ ] values. We present the paradigm and algorithms underlying our method, together with the proof of principle for the following statistical descriptors of [Mg 2+ ] distribution: weighted mean and weighted median concentrations (20), standard deviation, [Mg 2+ ] range, one or multiple [Mg 2+ ] modes (= [Mg 2+ ] curve maxima), kurtosis (peakedness (20)), skewness (asymmetry (21)) and entropy (smoothness (22,23)), ratios of [Mg 2+ ] mode heights and areas under individual [Mg 2+

  What are the ranges of the statistical descriptor values for [Mg 2+ ] distribution curves to be expected from actual experiments? How sensitive are the resulting [Mg 2+ ] distribution curves, and their associated descriptor values, to actual [Mg 2+ ] gradients vs. spurious effects, i.e., curve shape contributions unrelated to [Mg 2+ ] distribution? What are the effects of ATP-β 31 P MRS lineshape correction by our empirical deconvolution method on the resulting [Mg 2+ ] distribution curve shapes and, consequently, on the derived descriptors of [Mg 2+ ] profiles? How do [Mg 2+ ]

  P MR spectroscopy for intracellular [Mg 2+ ] measurement is based on the [Mg 2+ ] dependence of the chemical shift of the ATP-β resonance, δ ATP-β . The relationship between δ ATP-β and [Mg 2+ ] is given by the following equation (26) :

  descriptors characterizing [Mg 2+ ] distributions resulting from addition, in varying proportions, of up to three simulated [Mg 2+ ] curves, and (ii) to enable quantitative comparisons between heterogeneity parameter values derived from models of theoretical (i.e., simulated) and experimentally measured [Mg 2+ ] distributions. A range of well-defined [Mg 2+ ] distributions were computer simulated to calculate at least eight descriptors for each of these distributions. The [Mg 2+ ] distributions investigated were chosen as to represent different degrees of heterogeneity; they comprise distributions about up to three modes (= dominating [Mg 2+ ] values) of varying intensities. The Mg 2+ concentrations studied in our in-silico examples lied between 28 to 48 µM.

Figure 4 .

 4 Figure 4. Calibration curve ( 31 P NMR chemical shifts of ATP-β vs. free-Mg 2+ concentration)

Figure 5 .

 5 Figure 5. Processed ATP-β resonances from experimental 31 P NMR spectra of ATP/MgCl 2 solutions (see footnote to Table 1 for chemical composition). Top traces: ATP-β resonances from spectra of solutions containing 0.05 mM (red trace) or 0.12 mM (green trace) free Mg 2+ , measured individually in 10-mm tubes, and the sum of these resonances (blue trace) multiplied by a factor of two for better visibility. Bottom traces: ATP-β resonances from spectra of coaxial phantoms, each containing two different solutions with two different concentrations of free Mg 2+ :0.05 mM in the inner (5-mm) NMR tube, combined with 0.12 mM in the outer (10-mm) tube (green trace); and, alternatively, 0.05 mM in the outer tube, combined with 0.12 mM in the inner tube (red trace). Arrows indicate contributions to resonances from ATP-β in the presence of 0.05 mM (red arrows) and 0.12 mM (green arrows) free Mg 2+ .

Figure 2 .

 2 Figure 2.

24

 24

Figure 3 .

 3 Figure 3.

Figure 4 .

 4 Figure 4.

Figure 5 .

 5 Figure 5.

  Mg 2+ ]-dependent δ !"#!! may influence experimentally determined [Mg 2+ ] curves, in particular for very small [Mg 2+ ] gradients covering up to 10-20 µM. The consequences of contributions from factors other than [Mg 2+ ] effects on δ !"#!! are, of course, reduced precision of our [Mg 2+ ] curves and, consequently, reduced precision of the statistical parameters extracted from these curves. The limitations imposed by these imprecisions become significant when the [Mg 2+ ] range covered by the [Mg 2+ ] distribution curve is of the same order as, or smaller than, the uncertainty introduced by the spurious effects. However, a deconvolution procedure developed by us is able to largely compensate for these spurious effects on [Mg 2+ ] curves (see Supporting Information); further improvements taking into account [Mg 2+ ]-dependent linewidths due to intermediate Mg 2+ -ATP exchange rates are currently We have now adapted the underlying concepts and equations for application to [Mg 2+ ] heterogeneity. Briefly, we suggest a new strategy providing a number of statistical parameters describing [Mg 2+ ] heterogeneity, each derived by analyzing an entire [Mg 2+ ] curve. This strategy is based on the concept that the [Mg 2+ ] curve calculated from a ATP-β 31 P MR resonance can be approximated as a histogram. Such a histogram is formed by using the intensities of the digital points of a [Mg 2+ ] curve as heights of the corresponding histogram bars. Besides accurate, weighted-average [Mg 2+ ] values (= weighted mean concentrations, Mg !! ), also weighted median concentrations ( Mg !! ) are obtained in this way, as well as skewness (G1), kurtosis (G2), normalized entropy (H n ) and standard entropy (H S ) of [Mg 2+ ] distributions. In addition, distinct areas under [Mg 2+ ] curves can be calculated, each area representing the relative prevalence of a particular range of [Mg 2+ ] values of interest. Furthermore, for [Mg 2+ ]

	in progress.

[Mg 2+

] distributions presented here is based on the following general paradigm, originally developed by us for applications to two other measurement parameters, viz., pH and temperature

[START_REF] Lutz | Thermal heterogeneity within aqueous materials quantified by 1H NMR spectroscopy: Multiparametric validation in silico and in vitro[END_REF][START_REF] Lutz | Multiparametric quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: paradigms and algorithms[END_REF][START_REF] Lutz | Quantitative in-vivo characterization of intracellular and extracellular pH profiles in heterogeneous tumors: a novel method enabling multiparametric pH analysis[END_REF]

: Suppose that a sample is heterogeneous with respect to a measurement parameter, p. Further suppose that the chemical shift of an MR resonance is a function of p. Then, the statistical distribution of p values within the sample can be obtained by appropriate lineshape analysis of said MR resonance.

To derive a [Mg 2+ ] distribution from an ATP-β 31 P MR resonance, the resonance needs to be converted to a [Mg 2+ ] curve using eqs. 1 and 2 above, and corrected for the non-linearity between the chemical-shift and concentration scales. This [Mg 2+ ] profile then represents, in principle, the [Mg 2+ ] distribution within the measured volume. However, the resulting curve may also be influenced by other contributions to the lineshape: (i) by magnetic-field inhomogeneity and transverse relaxation (T 2 ) processes, in particular in the case of intermediate Mg 2+ exchange; (ii) by line broadening due to 31 P-31 P and 1 H-31 P coupling; (iii) potentially by uneven ATP distribution across microscopic regions of varying [Mg 2+ ]; and (iv) by the spectral processing parameters chosen, in particular filter parameters. 1 H-31 P coupling effects can be easily removed by broadband or CPD decoupling, and 31 P-31 P coupling by selective irradiation of the ATP-α and ATP-γ resonances. For fast exchange, the T 2 effect (natural line width) is often much smaller than the T 2 * effect (line width dominated by magnetic field inhomogeneity). Therefore it is suitable to focus on optimizing the magnetic field homogeneity to minimize T 2 * effects on line width, and to use adequate, consistent filter parameters to compare statistical descriptors of [Mg 2+ ] heterogeneity between measured samples. Even with these optimizations, lineshape contributions attributable to factors other than [In summary, a [Mg 2+ ] curve obtained with eqs. 1 and 2 accurately depicts the underlying [Mg 2+ ] heterogeneity, within the limits described above. Such a curve is, in fact, an envelope representing the sum of all microenvironments existing within the measured volume. Thus, all microenvironments combined cover the range of [Mg 2+ ] values given by the resulting [Mg 2+ ] curve. However, to this date no effort has been made to characterize and interpret such [Mg 2+ ] envelope curves, judging by published literature. In the following we propose to analyze [Mg 2+ ] curves by multiple statistical methods in order to extract a number of parameters (statistical descriptors) providing quantitative information on the nature of intracellular [Mg 2+ ] heterogeneity.

The most important innovation in our approach is its capacity to calculate, for the first time, quantitative parameters characterizing the statistical distribution of intracellular [Mg 2+ ] values in a given tissue volume. The general philosophy and algorithmic tools of similar methods have been presented by us in the context of the analysis of tissue pH and temperature heterogeneity

[START_REF] Lutz | Thermal heterogeneity within aqueous materials quantified by 1H NMR spectroscopy: Multiparametric validation in silico and in vitro[END_REF][START_REF] Lutz | Multiparametric quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: paradigms and algorithms[END_REF][START_REF] Lutz | Quantitative in-vivo characterization of intracellular and extracellular pH profiles in heterogeneous tumors: a novel method enabling multiparametric pH analysis[END_REF]

.

Experimental 31 P MR spectroscopy Aqueous Mg 2+ /ATP phantoms.

  Several test samples, designed to provide [Mg 2+ ]dependent 31 P MR spectra from ATP solutions, were generated as follows. Various amounts of phantom designed to mimic a bimodal [Mg 2+ ] distribution. Sample pH was adjusted to the physiological (intracellular) value of 7.2, and 31 P MR spectroscopic measurements were performed at 310 K (physiological body temperature), save for a few low-

	ATP (3.3 to 10 mM) were added to aqueous solutions containing 150 mM KCl buffered with 20 mM Tris, and MgCl 2 (up to 36 mM) was dissolved in the resulting solutions. The final solutions were used (i) for calibration purposes, and (ii) to generate phantoms representing two distinct exchange with ATP. 31 P MR acquisition and processing. Test samples were measured in 5 or 10-mm NMR tubes (Wilmad, Vineland, NJ, USA) at 162 MHz resonance frequency by way of a quadronuclear probe in a 9.4 T Avance NMR spectrometer (Bruker, Rheinstetten, Germany). Hard pulses (9.5 µs) were employed, and the resulting FID (32 k) was recorded for AQ=2.019 s with a sweep free-Mg temperature measurements at 278 K designed to study line width effects due to slow Mg 2+ width of 50 ppm (8.117 kHz), with the receiver gain set to RG=16 k. FIDs were zero-filled once

2+

concentrations at physiologically

[START_REF] Grubbs | Intracellular magnesium and magnesium buffering[END_REF]

relevant levels (combinations of 0.05 and 0.12 mM). The two solutions were transferred separately into the inner and outer coaxial NMR tubes of a combined and Fourier transformed after Lorentzian-Gaussian lineshape transformation (typical filter parameters for calibration spectra: Lorentzian broadening, LB=-4 Hz; Gaussian broadening, GB=0.01; for phantom spectra: LB=-30 Hz, GB=0.002). After phase and baseline corrections of the resulting spectra, post-processing (deconvolution), chemical-shift-to-[Mg 2+ ] conversion and calculation of the statistical descriptors of [Mg 2+ ] heterogeneity were performed by using our dedicated pre-programmed EXCEL spreadsheet template in conjunction with the StatPlus:mac module, version v5 (AnalystSoft, Walnut, CA, USA). A detailed user manual is embedded in said template (Supporting Information). The mouse tumor 31 P MR spectrum used as a first in-vivo example of our approach is a re-evaluated spectrum from a previously published paper

[START_REF] Lutz | Quantitative in-vivo characterization of intracellular and extracellular pH profiles in heterogeneous tumors: a novel method enabling multiparametric pH analysis[END_REF]

;

  field profiles, but to -1.1 units for the high-field profiles. The same qualitative trends hold for kurtosis increase and standard entropy decrease.Real-world [Mg2+ ] distribution curves may have rather complex shapes. Furthermore, even perfectly Gaussian or Lorentzian ATP-β 31 P MRS resonances would not be translated to equally shaped [Mg 2+ ] profiles because of the non-linearity between the chemical-shift (ppm) and Mg 2+ (molar) scales. Validation examples for simulated bimodal [Mg 2+ ] profiles composed of asymmetric functions are presented in Supporting Information, based on Figures S1 to S3.

2+ ] profiles -overview. Most tissues containing cells with unequal Mg 2+ concentrations are likely characterized by complex [Mg 2+ ] distributions, due to the low-Measured profiles of [Mg 2+ ] heterogeneity phantoms. First, we generated a calibration curve for the determination of free-Mg 2+ concentration by ATP-β 31 P MR spectra. This was achieved by acquiring spectra for a series of ATP solutions containing between 0 and 32.4 mM free Mg 2+ (see Methods for details). The calibration curve (Fig. 4) indicates that the optimum range of free [Mg 2+ ] for measurement by the ATP-β 31 P resonance extends up to about 350 µM; above 400 µM, chemical-shift changes as a function of [Mg 2+ ] variations become rather small and, therefore, yield more imprecise results. complex used for [Mg 2+ ]

  31 P MRS methods based on [Mg 2+ ] effects on the ATP-β chemical shift, do not take into account the [Mg 2+ ] heterogeneity existing within a tissue volume or voxel of interest. In the present study, we have comprehensively validated the first concept for measuring statistical distributions of free Mg 2+ concentration values in a heterogeneous environment. The following quantitative descriptors of [Mg 2+ ] heterogeneity were successfully evaluated by numerical simulation (in-silico experiments) and 31 P MR spectroscopy (experiments with well-defined model solutions): [Mg 2+ ] mean, median, mode(s), range and standard deviation; skewness, kurtosis and entropy; as well as ratios of peak heights and areas under [Mg 2+ ] distribution curves. The algorithms underlying our new method are now confirmed for uni-, bi-and trimodal [Mg 2+ ] distributions. The suggested approach provided consistent and robust results for [Mg 2+ ] ranges down to about 10-20 µM. The deconvolution-based curve shape correction concept for compensating spurious influences, Once confirmed by in-vivo experiments, potential medical applications of our quantitative [Mg 2+ ] heterogeneity measurements may include a range of conditions affecting the concentration of free magnesium(II) in tissue, notably in muscle, brain and cancer tissue. The new method can easily be integrated with 31 P MRSI, a method already being used in humans. This approach should help clarify the extent and nature of pathological processes involving major perturbations of [Mg 2+ ] regulation. [Mg 2+ ] distribution curves (= [Mg 2+ ] profiles). Second to fourth rows: changes in statistical descriptors of [Mg 2+ ] distributions with increasing mode 1 contribution (see alsoTable 2, top, for quantitative results). : changes in statistical descriptors of [Mg 2+ ] distributions with increasing mode 1 contribution (see also Table 3, top, for quantitative results).

	ACKNOWLEDGMENT
	This study was supported by Centre National de la Recherche Scientifique (UMR7339).
	originally introduced by us for evaluating temperature profiles, was validated for measured
	bimodal [Mg 2+ ] distributions.
	Figure 3. [Mg 2+ ] distributions as a function of relative contributions of mode 1 (mg 1 = 28.18 µM)
	and mode 2 (mg 2 = 48.18 µM), based on numerically simulated Lorentzian functions ATP-β 31 P
	NMR resonances (line width of corresponding ATP-β 31 P NMR resonances: 0.2 ppm for high-

CONCLUSION

We have validated, both in silico and in vitro, [Mg 2+ ] heterogeneity quantitation by

31

P MRS, and presented a preliminary example of an in-vivo application. Future work will comprise the adjustment of our deconvolution method to [Mg 2+ ]-dependent response functions, as well as advanced in-vivo validations.

FIGURE CAPTIONS

Figure 1. [Mg 2+ ] distributions as a function of relative contributions of mode 1 (mg 1 = 28.18 µM) and mode 2 (mg 2 = 48.18 µM), based on numerically simulated Gaussian functions (line width of corresponding ATP-β 31 P NMR resonances: 0.2 ppm for high-field spectra at 14.1 T). Top row: [Mg 2+ ] distribution curves (= [Mg 2+ ] profiles). Second to fourth rows: changes in statistical descriptors of [Mg 2+ ] distributions with increasing mode 1 (mg 1 ) contribution (see also Table 1, top, for quantitative results).

Figure 2. [Mg 2+ ] distributions as a function of relative contributions of mode 1 (mg 1 = 28.18 µM) and mode 2 (mg 2 = 48.18 µM), based on numerically simulated Gaussian functions (line width of corresponding ATP-β 31 P NMR resonances: 0.35 ppm for low-field spectra at 9.4 T). Top row: field spectra at 14.1 T). Top row: [Mg 2+ ] distribution curves (= [Mg 2+ ] profiles). Second to fourth rows

TABLE CAPTION (= title)Table 1 .

 CAPTION1 Experimental conditions for 31 P NMR spectra of aqueous ATP/MgCl 2 test solutions [Mg 2+ ] distributions as a function of relative contributions of mode 1 (mg1 = 28.18 μM) and mode 2 (mg2 = 48.18 μM), calculated based on numerically simulated Gaussian ATP-β 31 P MR resonances (line width, 0.2 ppm) for high-field spectra (at 14.1 T). Top

	SUPPORTING FIGURE CAPTIONS
	Supporting Figure S1.

row: [Mg 2+ ] distribution curves (= [Mg 2+ ] profiles). Second to fourth rows: changes in statistical

Table 1 .

 1 Experimental conditions for 31 P NMR spectra of aqueous ATP/MgCl 2 test solutions [MgCl 2 ] [Mg 2+ ]f δ ATPβ -δ ATPα Concentration units: mM. Samples were based on aqueous solutions of 3.3 mM ATP, 150 mM KCl and 20 Tris buffer, and pH was adjusted to pH 7.2. Concentrations of free Mg 2+ ([Mg 2+ ]f) were calculated based on added [MgCl 2 ], and eq. 1 with K d =0.0385 mM. Spectra were acquired at 9.4 T magnetic field strength and 310 K sample temperature.

	0	0	-10.9115
	1.019	0.05	-10.1887
	1.502	0.09	-9.5997
	1.768	0.12	-9.3960
	2.553	0.27	-9.0906
	2.869	0.30	-8.7578
	3.168	0.40	-8.5958
	3.402	0.50	-8.5558
	3.772	0.70	-8.4278
	3.600		-8.5211
	3.929	0.80	-8.4258
	7.200	4.00	-8.2900

genetic, physiological, metabolic and other growth conditions varying throughout the studied tissue volume. However, for the purpose of validating our method we initially focused on numerically simulated unimodal, bimodal and trimodal [Mg 2+ ] distributions, because the concept of quantitative parameters of [Mg 2+ ] heterogeneity is best tested on the basis of well-defined distribution functions. In fact, most quantitative parameters suggested in this report are S1 to S3). For symmetric distributions corresponding to high-field spectra, the locations of the modes, the calculated ratios of peak heights and areas under peaks as well as skewness are in agreement with theoretical values (top parts of Tables S1 andS3), and with previous validations obtained for analogous algorithms [START_REF] Lutz | Thermal heterogeneity within aqueous materials quantified by 1H NMR spectroscopy: Multiparametric validation in silico and in vitro[END_REF][START_REF] Lutz | Quantitative in-vivo characterization of intracellular and extracellular pH profiles in heterogeneous tumors: a novel method enabling multiparametric pH analysis[END_REF]. Note that for the Gaussian-based bimodal distributions discussed, standardized entropy, H S , is systematically increased by ca. 0.5 units for low-field vs. high-field curves. This is due to the fact that identical linewidths (on the frequency scale) were chosen for all simulations of individual peaks; this translates to relatively broad resonances on the ppm scale for 9.4 vs. 14.1 T. Consequently, low-field [Mg 2+ ] curves are made up of more digital points than the corresponding high-field curves (for identical sweep widths and spectrum sizes), and the higher number of points in turn results in systematically increased H S values, as explained previously [START_REF] Lutz | Multiparametric quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: paradigms and algorithms[END_REF]. By contrast, when normalized for the actual number of points, entropy (H n = normalized entropy) is very similar for 9.4 and 14.1-T [Mg 2+ ] curves, and the somewhat lower values for 9.4 vs. 14.1 T reflect the higher overlap of the modes for 9.4 T. For the simulated bimodal [Mg 2+ ] distributions chosen, it is obvious by visual inspection that the [Mg 2+ ] effects on the lineshape would be more pronounced for high-field (Fig. 1, top row) than by low-field spectroscopy (Fig. 2, top row). In the former, the two modes can be separately recognized even when one of the modes (mg 1 ) is just 1/8th the height of the other mode (mg 1 ).

Validation of potential and limitations based on simulated bimodal [Mg

By contrast, in simulated [Mg 2+ ] profiles for low-field MRS the mg 1 mode could be clearly discerned only when it had the same height as the mg 2 mode (Fig. 2, top right curve). Thus, the sensitivity of our method to [Mg 2+ ] gradients is evidently higher for high-field than for low-field profiles, as the latter profiles are largely dominated by spurious effects (see also top parts of Tables S1 andS2). For instance, the nominal value for normalized entropy, H n , decreases only by about 0.02 units for the low-field profiles, but by about 0.04 units for the high-field profiles when going from the symmetric distribution (height ratio 1:1) to increasingly asymmetric distributions (down to 1: 8). In addition, skewness decreases from 0 to only about -0.4 units for descriptors of [Mg 2+ ] distributions with increasing mode 1 (mg1) contribution (see also Table 1 for quantitative results).

Supporting

S12)

. In each of the profiles C to G, two partially overlapping areas could be discerned sufficiently well for separate integration (see also Table S5). S6.

Supporting


Norbert W Lutz
email: norbert.lutz@univ-amu.fr

Monique Bernard
email: monique.bernard@univ-amu.fr

Thermal heterogeneity within aqueous materials quantified by 1 H NMR spectroscopy: multiparametric validation in silico and in vitro

We recently suggested a new paradigm for statistical analysis of thermal heterogeneity in (semi-)aqueous materials by 1 H NMR spectroscopy, using water as a temperature probe. Here, we present a comprehensive insilico and in-vitro validation that demonstrates the ability of this new technique to provide accurate quantitative parameters characterizing the statistical distribution of temperature values in a volume of (semi-)aqueous matter.

First, line shape parameters of numerically simulated water 1 H NMR spectra are systematically varied to study a range of mathematically well-defined temperature distributions. Then, corresponding models based on measured 1 H NMR spectra of agarose gel are analyzed. In addition, dedicated samples based on hydrogels or biological tissue are designed to produce temperature gradients changing over time, and dynamic NMR spectroscopy is employed to analyze the resulting temperature profiles at sub-second temporal resolution. Accuracy and consistency of the previously introduced statistical descriptors of temperature heterogeneity are determined: weighted median and mean temperature, standard deviation, temperature range, temperature mode(s), kurtosis, skewness, entropy, and relative areas under temperature curves. Potential and limitations of this method for quantitative analysis of thermal heterogeneity in (semi-)aqueous materials are discussed in view of prospective applications in materials science as well as biology and medicine.

Introduction

We have previously presented a theoretical concept for the quantification of thermal heterogeneity in aqueous and semi-aqueous materials [START_REF] Lutz | Multiparametric Quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: Paradigms and Algorithms[END_REF]. This new method is based on statistical analyses of the 1 H NMR line shape of water, a "chemical temperature probe". The new paradigm underlying our technique is grounded in the fact that (i) the chemical shift of a water 1 H NMR resonance,  H2O , is a linear function of the water temperature (with a slope of ca. 0.01 ppm/°C [START_REF] Muller | Temperature dependence of chemical shifts of protons in hydrogen bonds[END_REF]), and (ii) in the presence of temperature gradients, the overall water NMR line shape reflects the statistical distribution of temperature values within the measured sample volume, within certain limits discussed previously in detail [START_REF] Lutz | Multiparametric Quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: Paradigms and Algorithms[END_REF]. To analyze thermal heterogeneity with this method, water 1 H NMR spectra are first converted to temperature distribution curves. Then, these temperature profiles are transformed into histograms by an appropriate algorithm. Based on these histograms, further mathematical algorithms developed by us derive at least eight quantitative parameters (descriptors) defining the underlying temperature distribution, i.e., the thermal heterogeneity within the measured sample volume.

The specificity of our new paradigm lies in the fact that we exploit, for the first time, the entire temperature curve (corresponding to an NMR line shape), to quantitatively describe a distribution of temperature values. Here, we present the validation of our method for the following descriptors of temperature distribution: weighted mean and weighted median temperatures [START_REF] Lutz | Multiparametric Quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: Paradigms and Algorithms[END_REF], standard deviation, temperature range, one or multiple temperature modes (= temperature curve maxima), kurtosis (peakedness [START_REF] Lutz | Multiparametric Quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: Paradigms and Algorithms[END_REF][START_REF] Sokal | Biometry[END_REF]), skewness (asymmetry [START_REF] Doane | Measuring Skewness: A Forgotten Statistic?[END_REF]) and entropy (smoothness [START_REF] Cover | Elements of information theory[END_REF][START_REF] Lesne | Shannon entropy: a rigorous notion at the crossroads between probability, information theory, dynamical systems and statistical physics[END_REF]). Note that temperature skewness, kurtosis and entropy directly provide actual temperature information, akin to temperature mean and standard deviation. From a statistical perspective, the latter two parameters are 1st and 2nd moment descriptors, respectively, of a given (temperature) distribution, whereas skewness and kurtosis are 3rd and 4th moment descriptors, respectively. In other words, skewness and kurtosis provide temperature information in the same way as mean and standard deviation, albeit on a "higher level", i.e., in the form of "higher moments".

In addition, ratios of mode heights, and ratios of areas under individual temperature modes and/or ranges are validated as quantitative measures of the relative sizes of volumes characterized by different temperature values. We assess the validity of our method based on (i) numerically simulated water 1 H NMR spectral lines; (ii) experimen-tally measured 1 H NMR spectra of well-defined aqueous test solutions and hydrogels (agarose gels [START_REF] Lutz | Phantom material for quantitative evaluation of MR images[END_REF]); (iii) dynamic 1 H NMR spectroscopy of hydrogels in the presence of temperature gradients changing over time, using a dedicated experimental design; and (iv) in-vitro 1 H NMR experiments on biological tissue.

Our simulations are entirely unrelated to conventional NMR lineshape simulations based on spin dynamics or Bloch equations. The purpose of our numerical simulations is to examine the following questions: 1) How are the quantitative descriptors of statistical temperature distribution obtained through our algorithm interpreted in terms of distribution curve properties? 2) What are the ranges of the statistical descriptor values for distribution curves to be expected from actual experiments? 3) How sensitive are the resulting temperature distribution curves, and their associated descriptor values, to actual temperature gradients vs. spurious effects, i.e., curve shape contributions unrelated to temperature distribution? 4) What are the effects of water 1 H NMR line shape correction by our empirical deconvolution method on the resulting distribution curve shapes and, consequently, on the derived descriptors of temperature profiles? 5) How do temperature distribution descriptors derived from simulated spectra compare with those from measured NMR spectra? Thus, our simulations do not serve to elucidate NMR physics, but to study temperature distributions, albeit derived from water NMR resonances.

All calculations are performed employing the preprogrammed EXCEL spreadsheet provided in our previously published paper [START_REF] Lutz | Multiparametric Quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: Paradigms and Algorithms[END_REF]. We also investigate how temperature profile measurements can be used to efficiently and quantitatively study dynamic phenomena such as heat transfer and thermal conduction, as water 1 H NMR spectra can be acquired with sub-second temporal resolution. Our findings for basic hydrogels and biological samples have significant implications for studying the relationship between temperature variations and other physical, chemical and biological properties in a variety of (semi-)aqueous materials. A broad scope of research ranging from materials science and physical chemistry of soft matter to biology and medicine [START_REF] Raue | Sodium NMR relaxation: a versatile non-invasive tool for the monitoring of phase transitions and the estimation of effective pore sizes of supramolecular hydrogels. Intelligent hydrogels[END_REF][START_REF] Kopecek | Hydrogel biomaterials: a smart future?[END_REF][START_REF] Mehta | The collagen I mimetic peptide DGEA enhances an osteogenic phenotype in mesenchymal stem cells when presented from cell-encapsulating hydrogels[END_REF][START_REF] Chaudhuri | Hydrogels with tunable stress relaxation regulate stem cell fate and activity[END_REF] may benefit from our approach.

Results and Discussion

The most relevant results of our study are presented in this section. Additional data are included in the Results section of S1 Supplemental Materials 1 that also comprises Figs. S-1 through S-29, and in S2 Supplemental Materials 2 with Figures S-30 through S-45 and Tables S-1 through S-5. Taken together, the values given in our tables indicate the ranges to be expected for statistical descriptors in experimental distribution curves.

Studies of unimodal, bimodal and trimodal temperature profiles -overview

When subjected to a heating or cooling process, most materials are likely to experience irregularly shaped internal temperature distributions, due to thermal gradients that are more or less continuous and variable over time.

However, for the purpose of validating our method we initially focused on numerically simulated unimodal, bimodal and trimodal temperature distributions, because the concept of quantitative parameters of thermal heterogeneity is best tested and verified on the basis of models representing well-defined temperature distribution functions [START_REF] Lutz | Multiparametric Quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: Paradigms and Algorithms[END_REF]. In fact, most quantitative parameters suggested in this report are universally applicable to any given distribution of temperature values; their use does not depend on the existence of distinct temperature modes. Our first in-silico validations were based on temperature modes at t 1 = 32.0, t 2 = 34.5 or t 3 = 37.0 °C (Table 1, for simulated Gaussian spectral lines at 11.7 T). Different bimodal (or trimodal) temperature distributions were modeled by varying the relative intensities or amplitudes i 1 (or i 1 and i 2 ) of 32.0 °C (or 32.0 and 34.5 °C) modes, added to an i 2 (or i 3 ) 37.0 °C mode. Subsequently, statistical descriptors were derived for unimodal (Table 1, column A), bimodal (columns B through E) and trimodal (columns F through I) temperature distributions, as a function of the relative mode intensities characterizing each distribution. For the unimodal and bimodal distributions of this series, the relative intensities of the 32 °C vs. 37 °C modes were used to scale the abscissae of the diagrams in Fig. 1. Analogous validations were performed for simulated Gaussian spectral lines at 3 T (Table S-

Quantification of oligomodal temperature distributions

We have previously shown that under high-field conditions, up to three distinct regions of similar size can be identified based on differing temperature values, within a temperature range close to 5 °C (no assumptions are made with respect to the size and spectral distribution of the underlying sample volume elements) [START_REF] Lutz | Multiparametric Quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: Paradigms and Algorithms[END_REF]. This finding takes into account that any experimentally obtained NMR resonance has a finite linewidth even in the absence of any temperature gradient, essentially due to T 2 and T 2 * effects [START_REF] Lutz | Multiparametric Quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: Paradigms and Algorithms[END_REF]. In the absence of such spurious effects, our simulated uni, bi and trimodal temperature distribution curves, and the associated spectra, would have the shape of one, two or three delta functions, respectively. The influence that removing these spurious effects has on statistical descriptors of temperature distributions is discussed in S1 Supplemental Materials 1. In the following, we focus on "raw" lineshapes, i.e., before any correction for spurious effects. changes in statistical descriptors of temperature distributions with increasing mode 1 contribution (see also Table 1).

In our simulated oligomodal curves, the corresponding mode areas may strongly overlap. The implications of this situation for temperature profile quantitation will be discussed in the following. Changes of quantitative temperature distribution descriptors as a function of the shape of bimodal temperature profiles are best visualized by diagrams plotting descriptor values against increasing contribution of mode t 1 to a distribution dominated by mode t 2 , as indicated by increasing intensity ratios, i 1 :i 2 (Fig. 1, Figs. S-2 through S-4 for simulated, and Figs. S-5 and S-42 for measurement-based distributions). For distributions derived from high-field spectra, t 1 and t 2 were rather stable because peak overlap was negligible for a difference of 5 °C (Figs. 1 and S-3 for simulated, and Figs.

S-5 and S-42 for measured distributions). By contrast, the relatively broad lines derived from simulated low-field spectra resulted in high-temperature mode values decreasing, and in low-temperature mode values increasing with increasing contribution of mode t 1 (Figs. S-2 and S-4). In many cases, the temperature of the less intense mode, t 1 , could only be estimated, notably for i 1 :i 2 = 0.125 and 0.25. As was to be expected for left-skewed distributions, < for i 1 :i 2 = 0.125, 0.25 and 0.5. For i 1 :i 2 = 1.0, only simulated distributions were perfectly symmetric and, therefore, yielded identical values for, and , whereas measured distributions were slightly right-skewed due to an increased right foot (or tail) of the underlying NMR resonances (spurious effect due to magnetic-field inhomogeneity). This resulted in somewhat smaller values for vs.

(Figs.

S-5 and S-42). In summary, mode errors were small but not entirely negligible for modes of small intensity located very close to one or two modes of much larger intensity. Since our simulation predicted rather poor resolution for low-field spectra, we decided not to acquire experimental low-field NMR spectra.

The concepts of skewness, kurtosis and entropy are generally employed to characterize unimodal probability distributions [START_REF] Jensen | Operations Research Models and Methods[END_REF]. However, they can also reveal global characteristics of temperature profiles with more than one mode [START_REF] Lutz | Multiparametric Quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: Paradigms and Algorithms[END_REF]. For instance, the symmetric distributions formed by two identical, simulated high-field resonances displayed in Figs. distributions continuously increased with i 1 :i 2 increasing from 0.125 to 1.0, i.e., with decreasing asymmetry as indicated by modes of unequal peak heights. Due to the relatively broad bases of Lorentzians, the symmetry of Lorentzian distributions is less influenced by the difference in heights between the two modes, than are Gaussian distributions. This is reflected by (i) relatively low absolute skewness values, and (ii) relatively small increases in skewness with i 1 :i 2 increasing from 0.125 to 1.0 for high-field Lorentzians vs. Gaussians. However, in low-field spectra, the bases of Gaussians also become relatively broad, which approaches the skewness behavior of Gaussians to that of Lorentzians. Kurtosis continuously decreased over the same range of i 1 :i 2 values, notably for distributions composed of narrow peaks, because the overall distribution becomes less peaked with increasing t 1 contribution. More generally, distributions made up of broad and, therefore, strongly overlapping mode areas vary little in shape with increasing i 1 :i 2 ratio; therefore, quantitative descriptors of curve shape also show relatively little change, although there was a drastic difference between the kurtosis values of a single narrow Gaussian (G2 = 0, Fig. 1) and a single narrow Lorentzian (G2 = 14, Fig. S-3). In summary, the skewness, kurtosis and entropy values of our bimodal curves can be readily interpreted in accordance with the underlying line shapes [START_REF] Lutz | Multiparametric Quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: Paradigms and Algorithms[END_REF].

For trimodal distributions with modes in 2.5 °C intervals, mode temperature values (t 1 , t 2 and t 3 ) did not change with intensity ratios, i 1 :i 2 :i 3 , when simulated high-field Gaussians were compared (Table 1, columns F -I). Very small variations (by ca. 0.1 °C) with intensity ratios were observed for measured high-field (near-)Gaussians (Table S-4, columns F -I), and for simulated high-field Lorentzians (Table S-2, columns F -I). For the remaining trimodal distributions tested, the resolution was often insufficient to determine the exact value, or even the presence, of all three modes. This confirms for trimodal distributions the same conclusions presented above for bimodal distributions, namely that working under high-resolution conditions (high magnetic fields, (near-) trend toward stronger left skew and increased kurtosis from column F to column I (Table 1, simulated high-field Gaussians). This reflects the gradual transition from a rather symmetric and flat temperature profile with equal mode intensities (i 1 =i 2 =i 3 ), to a more asymmetric and peaked profile with intensities strongly increasing from lower to higher temperatures (i 1 <i 2 <i 3 ). As a consequence, the overall temperature distribution curve increasingly takes on the shape of a steep slope (see Gaussian temperature modes, but less consistently so. Obviously, the increased overlap between mode areas, caused by broader bases in all three modes, modified the characteristics of temperature distributions significantly.

As opposed to skewness and kurtosis, standard entropy only exhibited relatively small variations within any given series of mode intensity variations. However, there were consistent differences between series, i.e., between distributions based on high-field and low-field NMR spectra, and between distributions based on Gaussian and Lorentzian-like areas under the modes. This finding underlines the validity and the complementary character of the statistical descriptors of temperature distributions used in this work. These and further results for multimodal temperature distributions are presented in more detail in S1 Supplemental Materials 1.

Effects of global vs. local heat exchange in gel phantoms

The first experiment with a dedicated temperature gradient sample (gel phantom) was performed as presented in and lower them into the NMR probe (8-10 s), resulted in significant heat exchange (i) between the NMR tubes and the surrounding air, and (ii) between the two phantom compartments, even before the first spectrum was measured. For this reason, the first temperature profile exhibited two broad, partially overlapping mode areas, instead of two narrow, well-separated mode areas, thus mimicking an irregularly shaped temperature distribution which may occur in many real-life situations of rapid sample heating or cooling. [Note that particular NMR singlets with temperature-invariant chemical shifts, used as chemical-shift references, were obtained simultaneously with the water resonances, and from the same sample volumes. These singlets were constant, narrow and regularly shaped, as was to be expected. See S1 Supplemental Materials 1 for details.] The hightemperature mode was below, and the low-temperature above the starting temperature of the respective individual tube. The behavior of the temperature profile was followed over a period of ca. 3 h, with increasing time intervals between measurements. The lists of intervals for all relevant experiments are given in S1 Supplemental Lists.

Five of a total of 146 subsequent measurements performed in this series were selected to document the most significant changes in curve shape occurring over time (Fig. S-40 a -e, height-normalized profiles; see also Fig. 2). In the initial phase of these measurements, the two mode areas progressively merged to form one unimodal, phantom compartments (= local effect), as is the subsequent narrowing of the curve. By contrast, the fact that the phantom temperature progressively approached the final temperature of about 60 C is owed to heat exchange with the surrounding gas stream, influencing and, ultimately, determining the temperature of both inner and outer tubes (= global effect). Obviously, both effects were not strictly separated in time as the sample heating by the gas stream must have started as soon as the phantom was exposed to the gas stream. However, local heat exchange initially dominated and was responsible for the early, major changes in the temperature profile, while the subsequent drift of the overall sample temperature to high values was caused by the then dominating global heat exchange with the gas stream. Note that the curve shapes discussed above are based on NMR spectra processed with Lorentzian filtering (1-Hz apodization). Alternative filtering parameter values and methods can slightly modify details of curve shapes, as discussed in S1 Supplemental Materials 1. Apodization with LB = 1 Hz, followed by enhanced baseline correction.

The general behavior of our phantom design was confirmed by quantitative evaluation as described in the following paragraphs. Results from the above experiment are shown in Figure 2 (height-normalized curve shapes)

and Figure 3 (values of descriptors of thermal heterogeneity). Both figures present data for the same time points, in the same measurement series. Until 14.1 s after the first measurement, there was little change in weighted mean temperature and almost no change in median temperature (Fig. 3, top left), suggesting very little net heat exchange between the phantom as a whole and its environment (gas stream). After 32.3 s, mean and median temperatures increased in a sigmoidal fashion, approaching the final temperature (see S1 Supplemental Materials 1 for the choice of time points). At 14.1 s following the first measurement, the low and high-temperature mode areas began to merge (Fig. 2), and at 28.8 s only one mode (= global temperature maximum; Fig. 3, bottom left) was distinguishable. Interestingly, neither entropy nor kurtosis changed uniformly over time (Fig. 3, top right).

Kurtosis reached a moderately pronounced maximum when the two modes were still clearly identifiable (6.4 s), then went through a minimum at 14.1 s (before complete merger of the mode areas), followed by an increase as the temperature distribution curve became very narrow (partial Lorentzian character). This behavior is linked to the well-known fact that kurtosis can increase as a consequence of two different influences: either due to tails extending far from the mode (6.4 s, left tail, Fig. 2), or due to a peak top that is narrow compared to the tails (end of series: 2h54m). Two phases can be identified in the temporal behavior of skewness: an almost constant, moderate left skew as long as the two modes were still clearly distinguishable (6.4 s), then a slow increase up to a slightly positive value for the final narrow curve, but still close to symmetry (G1 = 0). Thus, quantifying both skewness and kurtosis permits a more detailed description of the changing temperature profile than would be possible by simply using mean, median and mode temperature values. Range and standard deviation also changed in a biphasic fashion, albeit in a more uniform way, with inflection points close to the point of complete merger of the mode areas (between 14.1 and 28.8 s). As was to be expected, peak height ratios and peak area ratios reflected the relative decrease in t 1 mode area vs. t 2 mode area (a 1 vs. a 2 ) during the early phase of the experiment, where heat exchange between phantom compartments is dominating (Fig. 3, bottom right). As a result, the combination of multiple statistical descriptors of temperature heterogeneity was able to give a detailed quantitative view of the heat exchange processes occurring over time.

We complemented this analysis with an overview "at a glance", enabling easy visualization of the evolution of all temperature profiles over the entire period of serial measurements. This was achieved by plotting all profiles together as a stacked plot (Fig. 4 a), and connecting the individual profiles as to generate a virtual threedimensional surface plot (Fig. 4 b). Top projections of this surface plot give an immediate impression of the time course of the underlying temperature distributions (Fig. 4, c -e). However, it is also very instructive to view the three-dimensional shape of the entire surface plot from different angles. This is achieved by rotating the virtual three-dimensional object on a computer screen (S1 Horizontal rotation of 3D temperature profiles, and S2 Vertical rotation of 3D temperature profiles), thus fully exploiting the fact that the surface plot is assembled of a series of complete temperature profiles. Horizontal rotation by 360° (Video S-1), and vertical rotation by 90° from front view to top projection (Video S-2) can be performed under manual control by using the video scrollbar. Close inspection reveals that, for instance, temperature modes appear as "mountain ridges" in these pseudo-3D renderings. Following the appearance and disappearance of this and other 3D structural features as a function of measurement time provides additional insight into the dynamics of thermal processes in the material of interest.

Note that in these videos, as in Fig. 4 a and b, narrow temperature profiles appear significantly higher than broad profiles, which is a consequence of using area normalization (for details see S1 Supplemental Materials 1).

Effects of heat exchange in biological tissue

Prompted by previously published findings that the temperature dependence of water chemical shift (0.01 ppm per °C, [START_REF] Lutz | Multiparametric Quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: Paradigms and Algorithms[END_REF]) is virtually identical for aqueous solutions and biological tissue, we decided to test the congruence of the results obtainable for bovine muscle tissue as a model for mammalian tissue, with those obtained above for hydrogels. Bovine muscle tissue that was initially at thermal equilibrium, was first measured being externally cooled from 45 to ca. 0 °C, then being heated from 0 °C to 80 °C by analogy with the gel experiments mentioned above. Due to the relatively small temperature gradients generated during the cooling period, changes in line shape were very minor. By contrast, line shape alterations were readily discernable for the heating period (Fig. 5)

where the initial difference between sample temperature (ca. 0 °C) and N 2 stream (80 °C) caused marked temperature gradients across the sample. These gradients even resulted in two discernible, albeit heavily However, the temperature distribution remained rather symmetric throughout the measurement series; this is reflected by the fact that mean  median  mode(s)  global maximum temperature, and that skewness values are close to zero (at start of series: small and decreasing left skew) [START_REF] Lutz | Multiparametric Quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: Paradigms and Algorithms[END_REF]. Kurtosis was highest for the first measurement, where lines were very narrow (limit of zero temperature gradient). The decreasing values for peak area and height values are consistent with a globally increasing temperature. The qualitatively similar behavior of biological tissue and gel phantoms hints at the possibility to apply our method to in-vivo measurement in future projects.

deconvolution-based curve shape correction concept previously introduced by us was validated for an in-silico model of an asymmetric temperature distribution.

Future perspectives include full implementation and evaluation of our deconvolution method for compensating spurious influences on temperature profiles, and its adjustment to temperature-dependent response functions.

Potential applications of our suggested method concern a broad range of aqueous materials used in various disciplines, such as investigations of the formation and alterations of hydrogels and similar aqueous materials currently being developed in materials science. Furthermore, research and development in food science, quality control for industrial purposes, and biomedical procedures such as in-vivo examination of tissue properties may be of particular interest for our method [START_REF] Raue | Sodium NMR relaxation: a versatile non-invasive tool for the monitoring of phase transitions and the estimation of effective pore sizes of supramolecular hydrogels. Intelligent hydrogels[END_REF][START_REF] Kopecek | Hydrogel biomaterials: a smart future?[END_REF][START_REF] Mehta | The collagen I mimetic peptide DGEA enhances an osteogenic phenotype in mesenchymal stem cells when presented from cell-encapsulating hydrogels[END_REF][START_REF] Chaudhuri | Hydrogels with tunable stress relaxation regulate stem cell fate and activity[END_REF], notably in the context of tumor treatment by hyperthermia (e.g., thermoablation) or cryotherapy, and other processes involving major perturbations of temperature regulation.

Caution should be exercised if our technique is to be applied to special materials such as hydrated solids, e.g., porous materials loaded with water. Since our method has been developed for typical aqueous materials, it does not currently account for situations characterized by a nonlinear relationship between water 1 H chemical shift and temperature. However, the latter condition may be accommodated by future developments based on our approach.

Methods

In this section, the in-silico and in-vitro procedures employed in the validation of our new approach are described; additional technical details are presented in the Methods section of S1 Supplemental Materials 1.

In-silico calculations

The purpose of our in-silico calculations was twofold: (i) to study statistical descriptors characterizing the overall temperature distributions resulting from the addition, in varying proportions, of up to three simulated temperature curves, and (ii) to enable quantitative comparisons between heterogeneity parameter values derived from models of theoretical (i.e., simulated) and measured temperature distributions. A range of well-defined temperature distributions were computer simulated to calculate at least eight descriptors for each of these distributions [START_REF] Lutz | Multiparametric Quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: Paradigms and Algorithms[END_REF].

The temperature distributions investigated were chosen as to represent different degrees of heterogeneity; they range from narrow (very homogeneous) distributions about a single mode (= dominating temperature value) to broad (very heterogeneous) distributions about up to three modes of varying intensities. The temperatures studied in our in-silico H2O and temperature, distribution patterns are invariant to shifts on the absolute temperature scale. Further specifics on the rationale [START_REF] He | Selecting the number of bins in a histogram: A decision theoretic approach[END_REF][START_REF] Hargas | Sophisticated biomedical tissue measurement using image analysis and virtual instrumentation[END_REF][START_REF] Kohler | Data analysis using stata[END_REF][START_REF] Frigg | Entropy: A Guide for the Perplexed[END_REF][START_REF] Güçlü | Maximizing the entropy of histogram bar heights to explore neural activity: a simulation study on auditory and tactile fibers[END_REF] for the choice of our in-silico parameters are given in S1 Supplemental Materials 1.

Acquisition of 1 H NMR spectra

Gel phantoms for temperature gradients

Combined test samples designed to provide genuine temperature gradients, herein called phantoms, were generated as follows. An alkaline (pH 8.2) gel sample containing 1 % agarose and 80 mM PCho was heated to about 90 °C and filled into a 5-mm NMR tube which was then put on ice. In addition, an alkaline gel sample containing 1 % agarose and 80 mM NAA was heated to about 90 °C and filled into a 10-mm NMR tube, which was then maintained at 60 °C in a water bath. Subsequently, the heated 10-mm tube was removed from the water bath, and the cooled 5-mm NMR tube was immediately inserted into the 10-mm tube, upon which the combined phantom was quickly lowered into a 10-mm NMR probe. Simultaneously, a series of 1 H NMR experiments was started by way of home-made AU programs (Bruker macros based on the C programming language; see Appendix in S1 Supplemental Materials 1). These experiments consisted of pulse-acquisition sequences [START_REF] Lutz | Inflammatory multiplesclerosis plaques generate characteristic metabolic profiles in cerebrospinal fluid[END_REF] executed in increasing intervals ranging from 0.82 s to > 1 h. In the NMR probe, the sample was actively heated

Test samples consisting of biological tissue

To demonstrate the applicability of our approach to mammalian tissue, we have designed a test sample consisting of bovine muscle tissue (lean veal). Tissue was cut into small pieces that were progressively transferred to a 10mm NMR tube partially filled with phosphate buffered saline (pH 7.2). Occasional air bubbles were removed by moving the pieces of tissue in the buffer solution. After each addition, tissue was pushed to the bottom of the tube by exerting vertical pressure to minimize any space between tissue pieces, and between tissue and the wall of the tube. Once the available space was filled with muscle tissue to a height of 1-2 cm above the sensitive volume of the probe, the supernatant buffer solution was removed with a pipette. This procedure ensured that nearly all the material in the sensitive volume consisted of muscle tissue, and that the residual buffer volume was no more than a few percent of the total volume. The filled NMR tube was brought to ca. 40 °C in a water bath. 1 H NMR spectra were acquired according to the protocols described in the preceding paragraph for temperature gradient gel phantoms, first for temp nom progressively decreasing to 278 K, then progressively increasing from 278 K to 340 K. See S1 Supplemental Materials 1 and our preceding paper [START_REF] Lutz | Multiparametric Quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: Paradigms and Algorithms[END_REF] for further details concerning acquisition, processing and post-processing of 1 H NMR spectra [START_REF] Marshall | Use of Voigt lineshape for quantification of in vivo 1H spectra[END_REF][START_REF] Soher | Noninvasive temperature mapping with MRI using chemical shift water-fat separation[END_REF].

  1 and Fig. S-2); and for simulated Lorentzian spectral lines at 11.7 T (Table S-2 and Fig. S-3) and 3 T (Table S-3 and Fig. S-4). These and similar temperature distributions based on measured 1 H NMR spectra are discussed in S1 Supplemental Materials 1.

Figure 1 .

 1 Figure 1. Temperature distributions as a function of relative contributions of mode 1 (t 1 = 32 °C) and mode 2 (t 2 = 37 °C), based on numerically simulated Gaussian water 1 H NMR resonances (line width, 7.5 Hz) for a high-field

  1 and S-3 (top right curves) result in (nearly) zero skewness (Tables 1 and S-2, columns B). Also the less ideal, but theoretically symmetric temperature profiles based on two identical resonances, simulated for low field (Figs. S-2 and S-4, top right curves) or measured (Figs. S-5 and S-42, top right curves) had very low absolute skewness values (Tables S-1, S-3 through S-5, columns B). Predictably, skewness values [1] of bimodal

Fig. S- 6 ,

 6 bottom row, for visualization of intensity ratios used in all Tables, and in Figs. S-6 and S-43). Trimodal distributions composed of Lorentzian, low-field and experimental temperature modes, exhibited skewness and kurtosis trends that corresponded to those observed for high-field

Figure S- 41

 41 Figure S-41 a. Before insertion of the narrow tube into the wide tube, the former was brought to about -15 C, and

  broad temperature distribution curve (Fig. S-40 a -c). Subsequently, the line width of the resulting curve decreased over time, and the temperature mode asymptotically approached ca. 60 C. The emergence of a uniform, stable temperature throughout the gel phantom was imposed by a temperature-regulated gas stream (Fig. S-41 a). These results can be explained by a combination of local and global heat exchange. The initially observed thermal equilibration is the result of heat exchange between the two

Figure 2 .

 2 Figure 2. Selected height-normalized temperature profiles (with delays after first measured profile) obtained from

Figure 3 .

 3 Figure 3. Descriptors of temperature heterogeneity calculated from temperature profiles presented in Fig. 2.

Figure 4 .

 4 Figure 4. Generation of surface plots by combining all 136 temperature profiles (2 h 53 min starting from the

Figure 5 .

 5 Figure 5. Selected height-normalized temperature profiles (with delays after first measured profile) obtained from a series of high-field (9.4 T) water 1 H NMR measurements on bovine muscle tissue, contained in a 10-mm

overlapping, modes for 9

 9 to 27 s after the first measurement (Fig.S-34). Note that at the end of the experiment, the muscle tissue structure degraded due to heat exposure. At this time, we decided to reduce the N 2 stream temperature to 60 °C. This resulted in an atypically broad line shape and a slightly reduced measured sample temperature (Figs. 5 and S-34, last time points). This return to somewhat lower temperatures can also be seen in the corresponding surface plot (Fig.S-35 d). It is obvious that temperature range and standard deviation exhibit a broad maximum at 18 s after the first measurement (Fig. S-34), corresponding to the largest line width (Fig.5).

  by a 60 °C, or cooled by a 0 °C nitrogen stream, respectively, during the entire series of NMR sequences (Fig. S-41 a). The experiment was repeated with both tubes being heated together to about 90 °C (Fig. S-41, top right) and adjusted to 60 °C before starting NMR acquisitions as described above. During these NMR experiments, the sample was either kept in the probe initially adjusted to room temperature (Fig. S-41 a), without sample spinning or temperature regulation; or was actively cooled by a 0 °C nitrogen stream. In a fifth series of phantom measurements, the 5-mm NMR tube was inserted into the 10-mm NMR tube shortly after the start of data acquisition (Fig. S-41 b; NMR experiments as in preceding experiments, except for minimum delay of 0.41 ms between two subsequent spectrum acquisitions). This experimental design was suitable for close inspection of the fast temperature gradient changes occurring immediately after the combination of high and low-temperature gels in the phantom.

Table 1 . Statistical parameters characterizing simulated temperature heterogeneities.

 1 Gaussian line shapes) is essential for obtaining spectra with highest sensitivity to temperature gradients, and for minimizing spurious effects on temperature profiles. Unlike our bimodal distributions, trimodal distributions have not been generated by systematically varying only one specific intensity parameter value such as i 1 :i 2 , since the focus here was on providing examples for qualitatively different curve shapes. Nonetheless, there was a global

	descriptor		A	B	C		D	E
	weighted mean,	[°C]	37.0	34.5	35.3		35.9	36.4
	weighted median,	[°C]	37.0	34.5	36.5		36.7	36.8
	mode, t 1 [°C]		n/a	32.0	32.1		32.1	32.1
	mode, t 2 [°C]		37.0	37.0	37.0		37.0	37.0
	peak area ratios (a 1 :a 2 )
	-integrated		n/a	1.00:1	0.50:1	0.25:1	0.13:1
	-theoretical		n/a	1.00:1	0.50:1	0.25:1	0.125:1
	peak height ratios (h 1 :h 2 )	n/a	1.00:1	0.50:1	0.25:1	0.125:1
	skewness, G1		0.000	0.000	-0.6		-1.314	-2.009
	kurtosis, G2		0.000	-1.795	-1.319	0.210	3.124
	standard entropy, H S		4.13	5.13	5.05		4.86	4.64
			[6.10]	[6.64]	[6.59]	[6.47]	[6.35]
	descriptor		F	G		H		I
	weighted mean,	[°C]	34.5	35.2		35.5		35.8
	weighted median,	[°C]	34.5	36.2		36.2		36.4
	mode, t 1 [°C]		32.1	32.1		32.1		32.2
	mode, t 2 [°C]		34.6	34.6		34.6		34.6
	mode, t 3 [°C]		37.0	37.0		37.0		37.0
	peak area ratios (a 1 :a 2 :a 3 )
	-integrated		1.00:0.99:1	0.51:0.24:1	0.25:0.50:1	0.12:0.50:1
	-theoretical		1.00:1.00:1	0.50:0.25:1	0.25:0.50:1	0.125:0.50:1
	peak height ratios (h 1 :h 2 :h 3 )	1.00:1.00:1	0.50:0.25:1	0.25:0.50:1	0.12:0.50:1
	skewness, G1		-0.031	-0.541		-0.767		-0.891
	kurtosis, G2		-1.266	-1.238		-0.439		0.125
	standard entropy, H S		5.63	5.45		5.43		5.30
			[6.51]	[6.55]		[6.44]		[6.37]
	Data are based on Gaussian lineshapes, corresponding to 1 H NMR resonances at 11.7 T. A, single Gaussian, 37°C; B-E: two
	Gaussians, 32°C and 37°C; B, equal intensities; C, intensity ratio i 1 :i 2 = 1:2; D, intensity ratio i 1 :i 2 = 1:4; E, intensity ratio
	i 1 :i 2 = 1:8; F-I: three Gaussians, 32°C, 34.5°C and 37°C; F: equal intensities; G, intensity ratio i 1: i 2 :i 3 = 2:1:4; H, intensity
	ratio i 1: i 2 :i 3 = 1:2:4; I, intensity ratio i 1: i 2 :i 3 = 1:4:8. Entropy values given in [] are based on peaks simulated for the same

linewidth (in ppm) as for 3 T (Table S-1). a n , h n and i n : areas, peak heights and intensities, respectively, associated with modes t n , with n = 1 to 3.

This work was supported by the Centre National de la Recherche Scientifique (CRMBM UMR 7339).

Conclusion

The first concept for measuring statistical distributions of temperature values in aqueous materials, whose theoretical basis had been provided previously, has been validated in a comprehensive study. The following quantitative descriptors of thermal heterogeneity were successfully evaluated by numerical simulation (in-silico experiments) and water proton NMR spectroscopy (in-vitro experiments): temperature means, medians, modes, ranges and standard deviations; skewness, kurtosis and entropy; as well as ratios of peak heights and areas under temperature distribution curves. The algorithms underlying our new method are now confirmed for uni-, bi-and trimodal temperature distributions. We further demonstrated that the dynamics of statistical temperature distribution during heat exchange processes in hydrogel can be characterized quantitatively by analyzing the associated temperature profiles at sub-second time resolution. We also confirmed the applicability of our approach to biological materials (muscle tissue).

Fast data acquisition has the added advantage of providing far more realistic temperature snapshots than other, notably image-based, methods based on water proton chemical shift such as MRSI-PRF (magnetic resonance spectroscopic imaging -proton resonance frequency), as the latter methods require far more time for data collection. In the presence of thermal gradients, diffusing water molecules provide a mechanism for thermal equilibration, causing changes in the temperature profile over the data acquisition time. Thus, temperature profiles are "smeared out" across the temperature ranges concerned. In our method, "smearing" takes place over a much smaller time period (a few 100 ms) than it does in the currently available image-based experiments mentioned above. For instance, in MRSI-PRF data acquisition occurs over ca. 10 s per slice, and becomes much longer for an entire volume.

To visualize and guide these analyses, we constructed virtual 3D plots complemented by examples of virtual 3D videos for close inspection of changes in temperature profiles over time. The suggested approach provided consistent and robust results for temperature ranges down to a few degrees Celsius. For very small ranges, the

Supporting