Multiparametric quantification of thermal heterogeneity within aqueous materials by water 1H NMR spectroscopy: Paradigms and algorithms

Processes involving heat generation and dissipation play an important role in the performance of numerous materials. The behavior of (semi-)aqueous materials such as hydrogels during production and application, but also properties of biological tissue in disease and therapy (e.g., hyperthermia) critically depend on heat regulation. However, currently available thermometry methods do not provide quantitative parameters characterizing the overall temperature distribution within a volume of soft matter. To this end, we present here a new paradigm enabling accurate, contactless quantification of thermal heterogeneity based on the line shape of a water proton nuclear magnetic resonance (1H NMR) spectrum. First, the 1H NMR resonance from water serving as a "temperature probe" is transformed into a temperature curve. Then, the digital points of this temperature profile are used to construct a histogram by way of specifically developed algorithms. We demonstrate that from this histogram, at least eight quantitative parameters describing the underlying statistical temperature distribution can be computed: weighted median, weighted mean, standard deviation, range, mode(s), kurtosis, skewness, and entropy. All mathematical transformations and calculations are performed using specifically programmed EXCEL spreadsheets. Our new paradigm is helpful in detailed investigations of thermal heterogeneity, including dynamic characteristics of heat exchange at sub-second temporal resolution.

S-2 1. In-silico line shape simulation. EXCEL spreadsheets were programmed to simulate basic line shapes. For Gaussian lines, the following equation was used: For Lorentzian lines, the following equation was used: with parameter values as given above for Gaussians (unimodal distribution only).

Temperature curves as histograms; weighted mean and median temperature.
In an NMR spectrum presented as a histogram, each curve point is characterized by an abscissa value, t k , and a weight, W k , for k = 1 to k = m (Fig. 2 C). The histogram then reflects the relative prevalence of t k values in the measured sample. As indicated in the main body of this paper, the weight W = W k of a particular digital point k in a temperature curve or histogram S-3 reflects the relative frequency with which the temperature value represented by this point occurs in the measured sample ( Fig. 2 C). W k can be thought of as the sum of (real or hypothetical) individual contributions to a given digital point k of the temperature curve. It follows that is the sum of the weights of all digital points of the temperature curve [39]. By analogy to conventional histograms, n can be considered as the (real or hypothetical) total number of individual contributions to the temperature frequency distribution, and interpreted as being proportional to the total number of equal microscopic sample volume elements contributing to the weights W k of the digital points k of the temperature curve. The following equation was used for determining weighted mean temperature: where W k is the weight of digital point k of the temperature curve, t k is the discrete temperature value of digital point k, and m is the number of digital point k of the temperature curve.
As outlined in the main body, water 1   ! . Subsequently, the temperature curve point k whose cumulative sum equals the obtained half-sum is identified, i.e., the curve point k for which CSUM(k) = CSUM(h). The temperature value of curve point k = h, i.e., t k = t h , is defined as the weighted median temperature, temp.
Although this relationship may be directly applied to t k defined in the main part of this report, two adjustments are needed to obtain correct temperature values for h and, subsequently, accurate temp values 1 . First, the algorithm described above has originally been established for frequency distributions based on binned measurement of continuous parameters, as outlined in this report. In this context, the cumulative sum of each bin k includes all observations from parameter interval 1 through parameter interval k. By contrast, our temperature histograms are not based on continuous-parameter measurements that are assigned to bins because each point k represents a discrete temperature value. In fact, a t k value can be defined as a distinct data point (Fig. 2 C) located at the center of a conventional bin interval x k (Fig. 2 A). For this reason, the temperature values t k used for the identification of temp as described in the previous paragraph (Figs. S-1 A and B, black bars), need to be corrected by adding to each t k one half of i.e., the distance to the next curve point, such that t where t !" is the temperature value of curve point k h possessing the cumulative sum CSUM(k h ); t !"!! is the temperature value of curve point (k h −1) possessing the cumulative sum CSUM(k h −1); and f !"# is an interpolation factor defined as It is evident that temp = t !" in the limit where CSUM(k h ) = CSUM(h). Likewise, temp = t !"!! in the limit where CSUM(k h −1)= CSUM(h).

Temperature skewness, kurtosis and entropy.
Our equations for skewness and kurtosis calculation were adopted from the statistics module of the EXCEL spreadsheet and adapted to temperature distributions: where s is the nominal standard deviation, a parameter analogous to the standard deviation of the mean based on individual observations; t l is the temperature value of the l-th (real or hypothetical) individual contribution to the measured temperature distribution; l is a parameter analogous to the index that counts individual observations in conventional skewness and S-7 kurtosis calculation, from l = 1 to l = n. As pointed out elsewhere in more detail, skewness and kurtosis become independent of n for n greater than several times the number of digital points m 1 . Therefore, to obtain consistent results it is sufficient to scale the weights W ! such that n is much greater than m (W ! is typically on the order of several thousand or higher).
The interpretation of temperature profiles in terms of skewness and kurtosis is straightforward, based on the definitions presented above and in the main body. Skewness > 0 indicates a higher prevalence of high-temperature vs. low-temperature sample volumes. This implies that the low-temperature (< mean temperature) sample volume fraction consists of regions whose temperature values vary little and are close to the mean, while there is a relatively large temperature variability within the high-temperature sample volume fraction. The opposite applies to distributions with skewness < 0. For kurtosis > 0, a considerable proportion of the sample volume has temperature values that fall within a very narrow range around the center of the temperature distribution, while the residual sample volume has temperature values that are spread over a large range beyond the center (tail-heavy temperature distribution). By contrast, for kurtosis < 0 the bulk of the sample volume has temperature values distributed over a significant range around the mean (flat temperature distribution), but very little sample volume has extreme temperature values (shoulder-heavy temperature distribution). Accordingly, both skewness and kurtosis are parameters useful for characterization of temperature heterogeneity.
When the concept of entropy is applied to the analysis of the weights W k of a temperature distribution, S-8 where H n is the normalized entropy.

Temperature modes, ranges and volume regions.
In the presence of major temperature differences, more than two or three separate modes can potentially be detected. In the opposite case, i.e., when the temperature difference between two volume regions is on the order of 2.5 °C or smaller, it will be difficult to clearly separate modes in temperature curves derived from standard water 1 H NMR spectra. Instead, temperature "maxima" may appear as shoulders on broad lines, but can still be used to characterize temperature heterogeneity.
Relatively flat temperature distributions that do not possess well-defined modes may be frequently found in practical situations, but cannot be mimicked easily in silico or in vitro.
However, well-defined unimodal or oligomodal distributions are amenable to illustrating our paradigm and demonstrating the basic proof of principle.
Quantitative characterization of temperature distributions with multiple distinguishable temperature ranges can be achieved by calculating or estimating the overall physical sizes of these regions. Since the temperature curves are based on NMR spectra, the area under a temperature curve is, under appropriate measurement conditions, proportional to the amount of the compound measured. For this reason, we suggest to use 1 H NMR-derived temperature curves to quantify overall volume regions defined by specific temperature values. Relative sizes of individual areas under a particular temperature curve can be obtained as a result of (i) fitting of multiple analytical curves in some favorable cases, or (ii) integration over defined temperature ranges. The relative areas (e.g., area ratios) calculated in this way directly yield the corresponding relative volumes, provided that the free-water content does not vary significantly across the regions studied. To obtain accurate area ratios by integration (ii), intensities of curve points are summed over chosen temperature ranges. Fitting of multiple analytical functions with Gaussian/Lorentzian line shapes can be achieved with conventional NMR spectrometer software. Although the latter procedure is frequently referred to as "deconvolution", it should not be confounded with the homonymous mathematical operation also known as "faltung" (folding) which is used by us in the context of "deconvolution" for water signal correction (see below). Here, we introduce concepts designed to cancel out these spurious effects. The feasibility of our approach is demonstrated below; its full implementation, as well as the evaluation of its potential and limitations are an ongoing effort.

Compensation of spurious line shape modifications by using the convolution theorem.
As outlined in the main body of this manuscript, temperature curves derived from water proton NMR spectra may be significantly influenced by effects of T 2 , T 2 *, magnetic-field inhomogeneity and other factors. All of these define NMR line shape and line width in the absence of temperature gradients. In fact, measured "raw" temperature curves can be thought of as "true" temperature distribution curves (defined solely by the temperature dependence of δ H20 ) combined with concomitant line shape contributions caused by said spurious effects. It is therefore possible to correct experimentally obtained raw temperature distribution curves for modifications caused by spurious effects. As a first approximation, we propose to use a S-10 reference water NMR resonance from the sample in question, at thermal equilibrium. During the acquisition of this spectrum, the sample temperature should be close to the average temperature of the subsequent measurement(s) to be performed in the presence of temperature gradients.
The temperature curve(s) derived from the latter experiment(s) should then be corrected by taking into account the line shape and line width of the reference spectrum obtained at thermal equilibrium.
The separation of two concomitant contributions to NMR line shapes can be obtained by a process known as deconvolution. In physical and mathematical terms, deconvolution of two and, inversely, the Fourier transform of the product of two functions is equal to the convolution of their individual Fourier transforms where ! = signal function, ! = response function, and ℱ = Fourier transform.
Both convolution and deconvolution are widely applied in high-resolution NMR spectroscopy, in particular for spectral filtering (apodization) and resolution enhancement by Lorentzian-Gaussian line shape transformation. We suggest to apply an equivalent approach to the deconvolution of the raw temperature distribution curve (convolved function, ! ! " ) with the reference curve (response function, ! ). This is to be achieved by dividing the FID acquired in the presence of temperature gradients by the inverse Fourier transform (iFT) of the reference spectrum. The The suggested deconvolution procedure is demonstrated in Fig. S-2 for a numerically simulated model. In our example, both signal and response functions are Lorentzian. FIDs were modeled by exponentially decaying sinusoids, and the corresponding Lorentzians were obtained by

S-12
Fourier transformation. Note that the FID corresponding to the response function has the form of an exponential curve (equivalent to an on-resonance FID). The development of appropriate algorithms for these and further sophisticated corrections will be included in future projects centered about quantitative characterization of temperature heterogeneity by water 1 H NMR spectroscopy.

Principle of intravoxel vs. intervoxel thermal heterogeneity.
In the context of this study, it is important to distinguish between 'temperature regions' (corresponding to temperature ranges) that can be derived from temperature curves, on one hand, and physical regions within a sample on the other. The temperature distributions obtained from 1 H NMR spectra approximate histogram-like temperature profiles. Therefore, they do not contain information on the spatial organization of the underlying sample. In particular, no assumptions are made as to the physical shapes, sizes or locations of the physical macro or microregions responsible for the temperature regions, beyond the trivial fact that, e.g., all physical regions measured in one single-voxel 1 H NMR experiment are contained within the same macroscopic voxel (voxel = volume element). For instance, the appearance of two distinct temperature regions in a particular temperature profile may be generated by two large volume regions, each one characterized by a distinct characteristic temperature range. However, the same temperature distribution pattern may be produced, in principle, by a mixture of a large number of small volume regions whose temperature values fall into one of the distinct temperature ranges

S-13
represented by the temperature profile. In either case, the ratios of volumes characterized by different temperature values can be estimated as described above.
Our technique takes into account thermal heterogeneity independently of its spatial distribution, as opposed to heterogeneities based on temperature maps (such as temperature MR imaging [3][4][5][6][7][8] ). One of these currently existing temperature mapping techniques (MRSI-PRF, magnetic resonance spectroscopic imaging -proton resonance frequency 8  usually represent an average temperature value for the entire voxel ( Fig. S-4, right panels).
Although intravoxel temperature gradients may be negligible for (i) small thermal gradients throughout the sample and (ii) very small imaging voxels, this would no longer be true for larger thermal gradients within the sample and/or moderate spatial resolution as obtained in many MRSI-PRF experiments.
Thermal macroheterogeneity refers to temperature heterogeneity that is detectable in temperature maps because it is based on temperature differences between ( above). In the new statistical temperature profiling method presented in this report, the measured voxel comprises the entire volume of interest, unless multivoxel 1 H NMR spectroscopy is performed which produces one spectrum per voxel (briefly discussed in the following subsection). In fact, the issue of macroheterogeneity vs. microheterogeneity is not specific to temperature maps, but is common to all voxel-based imaging 1,9 .
Bottom right: Image-based histograms representing the temperature for the green and red-framed segments in our model, based on apparent voxel temperature values (see green and red arrows). The hue of each histogram bar represents its apparent temperature, which is identical to that of the corresponding voxels and segments. Note the presence of only one bar in each histogram; this artifact is a result of the partial volume effect in this example of microscopic thermal heterogeneity, and may erroneously suggest thermal homogeneity within the corresponding voxels and segments of the body volume. The height n rel of each histogram bar is normalized to the corresponding histograms at the bottom left.

Comparison of spectrum-based vs. image-based thermal heterogeneity
analysis. In our approach to quantifying thermal heterogeneity, all descriptors are based on digital points of temperature curves derived from water 1 H NMR spectra. A crucial aspect of our new paradigm is that these quantitative statistical parameters describe global features of thermal heterogeneity within a selected volume. To the best of our knowledge, none of the currently available temperature measurement methods provides such descriptors. As opposed to voxelbased imaging, our spectroscopic approach indiscriminately takes into account both macroscopic and microscopic heterogeneity, and allows for very fast measurements. Changes in temperature profiles over time can be followed at a rate of several measurements per second, but the temporal resolution is ultimately dependent on physicochemical properties of the aqueous material to be studied. However, the proposed method can be easily combined with one such method that is also based on water 1 H NMR (MRSI-PRF). In this way, both statistical and spatial information on temperature distribution can be gained from the same sample, provided the available NMR instrument is equipped for imaging. In fact, our new method could be used to directly analyze water 1 H NMR spectra obtained through MRSI-PRF, albeit at the expense of temporal resolution.