Improving the Accuracy of Model-based Quantitative NMR

Abstract. Low spectral resolution and extensive peak overlap are the common challenges that preclude quantitative analysis of NMR data with the established peak integration method. While numerous model-based approaches overcome these obstacles and enable quantification, they intrinsically rely on rigid assumptions about functional forms for peaks, which are often insufficient to account for all unforeseen imperfections in experimental data. Indeed, even in spectra with well separated peaks whose integration is possible, model-based methods often achieve suboptimal results, which in turn raises the question of their 5 validity for more challenging datasets. We address this problem with a simple model adjustment procedure, which draws its inspiration directly from the peak integration approach that is almost invariant to lineshape deviations. Specifically, we aim to recover all useful signal left in the residual after model fitting and use it to adjust the intensity estimates of modelled peaks. We propose an alternative objective function, which we found particularly useful for correcting imperfect phasing of the data – a critical step in the processing pipeline. Application of our method to the analysis of experimental data shows the accuracy 10 improvement of 20-40% compared to the simple least squares model fitting.


Introduction
Proposed thirty years ago (Miller and Greene, 1989;Bretthorst, 1990;Chylla and Markley, 1995), model-based approaches for quantitative NMR data analysis (qNMR) are getting wider acceptance as an effective alternative to the established peak integration (Kriesten et al., 2008;Krishnamurthy et al., 2017;Kern et al., 2018). Based on the idea that an experimental 15 spectrum can be represented as a collection of parametric lineshapes, e.g. Lorentzians with certain positions, widths, and heights, they offer a principled mechanism to resolve overlapping peaks and are less susceptible to noise (Matviychuk et al., 2017). By adjusting its parameters, the model is fitted to match experimental data, which eventually determines the sought concentrations of chemical species in the analysed sample. The reduction of a spectrum to a frequency-intensity table of peaks (Krishnamurthy, 2013) allows for easier automation of post-processing tasks and simplifies the analysis of large arrayed 20 datasets (Kriesten et al., 2008;Alsmeyer et al., 2004). Finally, quantum mechanical formulations minimize the number of free model parameters and are inherently invariant with respect to the spectrometer field strength (Kuprov et al., 2007;Tiainen et al., 2014;Dashti et al., 2017); they enable the analysis of highly complex low-resolution spectra acquired on medium-field benchtop instruments and are found successful in modern practical applications (Matviychuk et al., 2019).
There have been proposed numerous model-based approaches to the problem of qNMR formulated either in the time ( Vanhamme et al., 2001;Krishnamurthy, 2013;Rubtsov and Griffin, 2007) or the frequency domain (Mierisová and Ala-Korpela, 2001;Poullet et al., 2008). Notably, the latter typically depend upon phase and baseline correction of the spectra before fitting signal models to them (Cobas et al., 2008;Kriesten et al., 2008). In contrast, time-domain methods that work with the FID signal are often regarded as being able to lift this requirement (Krishnamurthy, 2013). However, we note that when the model fitting is performed in the least squares sense -as done most commonly -both variants of the problem formulations are 30 equivalent and result in the same solution. Thus, even though explicit data preprocessing steps can be often obviated by timedomain methods, they inevitably include the phasing parameters in a certain form, either as angles of complex-valued intensity estimates for separate resonances (Krishnamurthy, 2013;Rubtsov and Griffin, 2007;Kung et al., 1983) or as independently optimized parameters of a linear phasing model (Matviychuk et al., 2017). On the other hand, the phasing parameters can also be estimated from the complex-valued frequency domain data (Sokolenko et al., 2019). Similarly, the baseline effects that are 35 often observed over wide spectral ranges, appear in the leading time points of the original FID signal. Hence, these distortions also need to be taken into account in the time domain analysis, either by masking or weighting the early time samples.
Despite numerous advantages, model-based qNMR is often found suboptimal in seemingly easy cases: when peaks in the spectrum are well resolved, and the signal-to-noise ratio (SNR) is sufficently high, the peak integration after careful phase and baseline correction typically achieves higher quantification accuracy, as we observe later in Sec. 4.1. This can be explained 40 by the high sensitivity of most model-based qNMR algorithms to any unforeseen distortions in the experimental data, such as imperfections of peak shapes and their deviations from the assumed ideal Lorentzians. Indeed, model misspecification and its inability to faithfully represent the data biases the estimates of concentrations along with the associated uncertainties (White, 1981(White, , 1982Grünwald and van Ommen, 2017); this produces misleading results becoming one of the major points of criticism of model-based qNMR. To overcome this obstacle, several generalizations of the peak lineshape function have been proposed 45 over time, most notably the Voigt lineshape (Humlíček, 1982;Marshall et al., 1997;Bruce et al., 2000) and other combinations of Lorentzian and Gaussian terms (Kriesten et al., 2008;Schoenberger et al., 2016). Nevertheless, peak shape deviations in experimental spectra can often be very hard to model explicitly within the parametric framework, as they typically reflect multiple independent physical processes, such as diffusion, magnetic field inhomogeneity, higher-order coupling effects, etc. Reference deconvolution methods (Morris et al., 1997;Metz et al., 2000;Osorio-Garcia et al., 2011) offer an effective mechanism to 50 eliminate complex distortion patterns common for all peaks in the spectrum, e.g. arising due to the lack of shimming. However, they can not easily address possible differences in shapes of separate peaks, for example as a result of small long-range couplings, whose effects become even more noticeable at lower magnetic field strengths.
Since model-based qNMR is the only viable option for the analysis of complicated spectra with multiple overlapping peaks, it is of utmost importance to develop accurate algorithms that are robust to possible model misspecifications. The main goal of 55 this work is to bridge the performance gap between the peak integration and model-based qNMR by combining the strength of both approaches. Specifically, after fitting a model to the data, we observe that the residual -instead of being purely noise -often contains non-stochastic elements pertinent to the useful NMR signals. We propose to explicitly incorporate this unaccounted remainder into the model-based analysis, as would have been done with peak integration. On the other hand, Preprint. Discussion started: 28 January 2020 c Author(s) 2020. CC BY 4.0 License. neither of the existing phase and baseline correction methods take into account prior information about the studied system, which is conveniently employed in our approach in the form of an adjustable model. As a result, our alternative optimization procedure achieves better baseline and phase correction than the usual least-squares model fitting and improves model-based quantification of both well-resolved and overlapped data.
In the next section, we briefly review the main idea of model-based qNMR and introduce the notation for the problem of estimating the concentrations of components in a mixture. We then proceed by studying the weaknesses of the traditional least-65 squares model fitting and propose our alternative optimization criterion. Section 3 describes the setup for our simulations and laboratory experiments; their results are presented in Section 4.

Theory
This section provides a theoretical background for our method. First we review the general principle of qNMR: given a mixture of known chemical species, we are set to estimate their unknown relative concentrations (mole fractions) using the NMR data.

70
In model-based qNMR, an ideal model that represents the studied mixture is fitted to the experimental data, and its found optimal parameters -specifically the intensities of model components -determine the estimates of concentrations of chemical species. Here we discuss the consequences of model misspecification, and propose our model adjustment method to improve the accuracy of quantification.
2.1 Overview and the main idea of model-based qNMR 75 We choose to formulate our method in the frequency domain using real-valued spectra. Even though discarding the imaginary counterpart of complex-valued data entails reduction of SNR by a factor of 2, this will allow us to develop an adjustment algorithm for our model-fitting method inspired by the peak integration, which traditionally operates with real-valued spectra. Thus, an experimental spectrum is obtained using the discrete Fourier transform of the acquired FID followed by the usual first-order phase correction with parameters ϕ 0 and ϕ 1 . It is formally represented as an n×1 column vector y = Re F (y T ) e −i(ϕ0+ϕ1f ) , 80 where f = − 1 2 ∆t ≤ f ≤ 1 2 ∆t is the vector of frequency values corresponding to the particular sampling (dwell) time of the FID, ∆t.
Next we define a model matrix Z whose columns contain signature spectra for all K analyzed chemical species evaluated on the same frequency grid f . A typical signature model is a combination of P elemental peaks with different chemical shifts, widths, and intensities b p : (1) Here u p (f |f p , α p ) defines an ideal Lorentzian peak with central frequency f p and full width at half maximum αp π (both expressed in Hz) evaluated at the frequency f , We note that f p = B 0 δ p − f 0 , where B 0 and f 0 are the operating frequency of spectrometer in MHz and the spectral offset 90 respectively, which are used to convert the frequency units of the chemical shift δ p from ppm to Hz; α p is the decay rate of the corresponding FID signal in the time domain. Chemical shifts and widths, at least for certain peaks, can vary independently and usually reflect the specifics of experimental conditions. For example, the chemical shift of the proton in a hydroxyl group is famously related to the pH value of the sample. On the other hand, relative intensities b p of peaks pertaining to the same chemical necessarily remain constant, as they are defined by the atomic composition of the molecule. In the present work, we 95 use the quantum mechanical approach for modelling the signatures of chemical species (Matviychuk et al., 2019). It allows us to minimize the number of free parameters and produce relevant model spectra at any field strength of the spectrometer.
Finally, to account for a possibly imperfect baseline in the experimental data, we augment the model matrix Z with several basis vectors of the form f l−1 for l = 1, . . . , L, which serves to model any polynomial baseline of degree up to L (we use L = 1 in all our experiments in Section 4 to correct for the constant offset in the spectra).

100
With this notation, the complete model spectrum is expressed as x = Zc, where c is a vector of component intensities.
The main idea of the model-based quantification is to find a model x that is as close to the measured data as possible; the corresponding vector of intensities c is used to estimate the concentrations. To formalize this idea, we define the residual spectrum r = y − x and note that r implicitly depends on the set of model parameters -chemical shifts, peak widths, as well as the phasing values -which we denote collectively as θ. The model fitting is typically done in the least-squares sense by 105 minimizing the Euclidean norm of the residual: It is well known that, given the model matrix associated with the optimal set of model parameters Z, the vector of intensities can be estimated in closed form,

110
where Z is obtained as a result of unconstrained minimization of the non-linear variable projection functional L = I − Z T Z y 2 .
It can be shown that the criterion of Eq. 3 stems from the assumption that the measured signal is generated as an instance of the model affected by isotropic Gaussian noise, y = Zc + n. This plausible assumption is supported by the principle of maximum entropy and the central limit theorem, which made the least squares minimization -along with the existence of a simple solution -the most popular setting for the model fitting problem. However, as any mathematical model of the physical 115 world, this approach has certain limitations. We discuss them in more detail in the following subsection.

Model misspecification
The optimality conditions of the least squares fit only hold if the assumed signal model is capable of describing the experimental data given some set of parameters. Unfortunately, the most common assumption that underlies model-based qNMR -that an FID decays mono-exponentially producing spectral peaks with simple shapes -often does not hold in practice. Such effects 120 as diffusion and the magnetic field inhomogeneity cause the resulting peak shapes to deviate from the ideal Lorentzians, Preprint. Discussion started: 28 January 2020 c Author(s) 2020. CC BY 4.0 License. to which a model of Eq. 1 can no longer be perfectly fit. In turn, this leads to incorrect estimates of the intensities of the components c and erroneous (biased) quantification results (White, 1981;Deegan Jr., 1976). This stimulated the development of more complex signal models that account for second and higher order effects in the FID, such as Voigt (Marshall et al., 1997), generalized Lorentzian-Gaussian (Kriesten et al., 2008;Alsmeyer et al., 2004;Schoenberger et al., 2016), or flexible custom lineshapes in numerous spectral deconvolution methods (Cobas and Sýkora, 2009). These approaches were found to be very successful in cases when different peaks in the spectrum, even if they overlap, can be attributed to separate resonances with similar distortions, as often seen in high-resolution data acquired with a high-field instrument. Unfortunately, at the medium field strengths of benchtop instruments, these approaches become less effective. Higher-order coupling between neighboring and distant protons often cause different 1 H peaks to show different asymmetric distortions due to separation of transition 130 resonances (Kuprov et al., 2007). Quantum mechanical models were found useful for describing such data but also can not guarantee the perfect fit of complex spectra (Tiainen et al., 2014;Matviychuk et al., 2019).
In this work, instead of trying to refine the peak shape model, which can complicate the analysis and often bears the risk of overfitting, we propose to alter the optimization criterion in order to completely remove any unaccounted signal from the residual. As an illustrative example, in Fig. 1, we consider a spectrum of thiamine in D 2 O acquired on a high-field spectrometer.

135
The high spectral resolution and low level of noise in this dataset make it possible to achieve very accurate quantification results with conventional peak integration. Surprisingly, this appears to be a difficult case for a simple model-based method. The top panel in Fig. 2 demonstrates the least squares fit of Lorentzian peaks to the measured data obtained by minimizing Eq. 3 with respect to the positions and widths of all peaks and the phasing parameters. Close examination of the fit reveals significant deviations between the experimental and fitted peak shapes. To compensate for the model misspecification, the least-squares 140 fitting distorts the phasing of the spectrum and introduces a notable offset in the baseline. Even though these imperfections are relatively small, less than 1% of the average peak height, they are comparable to the level of random noise and can affect quantification. Furthermore, the mismatch between the model and the data can be easily observed in the residual spectrum: instead of being purely random Gaussian, as postulated in the model assumptions, it is dominated by large spikes where the model peaks do not fit the data perfectly. The resulting magnitude range of the residual is approximately 100 times higher 145 than the actual noise floor and is at least 20% of the average peak height. Thus, even though the found model admits to the requirements of the optimization criterion, it can not completely explain and account for the measured spectrum. Finally, we note that in this, and many other examples, more flexible peak models (e.g. Lorentzian-Gaussian) still do not eliminate the misspecification error completely.
This observation motivates our proposed approach and distinguishes it from other model-based quantification algorithms: 150 instead of relying on the top-down fitting of a supposedly ideal model, we employ a bottom-up view and aim to find a model spectrum, which after subtraction from the experimental data would lead to exclusively noise in the residual. We present our solution in the following subsection. Overlapping peaks Figure 1. Examples of spectra of thiamine acquired with a high-field (top) and a medium-field (bottom) spectrometers.

Outline of the proposed adjustment algorithm
The above example demonstrates that the conventional least squares minimization criterion, while being convenient to use, To develop our solution, we start with the least-squares model fit as described above and represent the residual spectrum r as a sum of three distinct components: a signal remaining solely due to the imperfect model fit that could potentially be accounted 160 for with more flexible signal models, a slow changing residual baseline that arises to compensate for the imperfect fit of the Preprint. Discussion started: 28 January 2020 c Author(s) 2020. CC BY 4.0 License.
peaks, and the random noise: Our strategy is to isolate the first term in the above decomposition, r m , and incorporate it directly into the fitted model, adjusting the corresponding component intensities c. 165 We draw the inspiration for our method from the conventional peak integration procedure and note that if the spectrum y is perfectly phased, its total area under the curve can be found as the sum of all fitted models and the misfit term of the residual, accordingly. Specifically, we define the resulting component intensities after adjustment as: for each k = 1, . . . , K, where W is an n × K matrix of row-normalized non-negative weights, k W i,k = 1, that determine the allocation rule of the residual among the K components at each point in the spectrum i = 1, . . . , n. Note that if the model is fitted perfectly, and r m = 0, the normalization in the adjustment rule of Eq. 6 preserves the original intensities, 175 c = c. In our experiments in Section 4, we found it particularly effective to assume that the misfit error of each component is proportional to its value at frequency i; then the allocation matrix is defined as: With the assumption that Eq. 6 is capable of recovering the true model intensities, the model adjustment problem reduces to the isolation of the misfit term r m in Eq. 5. For this, we start by explicitly removing the random noise from the residual 180 spectrum, which can be accomplished with any suitable 1D denoising algorithm. We found that soft thresholding of wavelet coefficients is particularly effective for this purpose (Donoho, 1995): it removes the stochastic deviations but preserves the spiky features of the residual that are due to model misspecification. In our examples in Section 4, we use symlets with eight vanishing moments and set universal thresholds proportional to the level-dependent estimates of noise on each wavelet decomposition level. The resulting signal after denoising, r = D (r), is assumed to be purely deterministic.

185
Next, we proceed by smoothing the denoised residual to extract its slowly-changing component, r b = S (D (r)), which encompasses the error introduced by an incomplete baseline correction. While any type of low-pass filtering can be used for this purpose, it is known that median filters -i.e. replacing each point in a signal with the sample median of its w m neighbours -are especially suitable for removing sharp spike artifacts (Tukey, 1977;Mallows, 1979). To summarize, we define r m = D (r) − S (D (r)), however, we note that the representation of a residual according to Eq. 5 is inherently an ill- Preprint. Discussion started: 28 January 2020 c Author(s) 2020. CC BY 4.0 License.
As in peak integration, it is paramount in our method that the analysed spectrum y is perfectly phased before adjusting the models. Fortunately, Eq. 5 provides a convenient way for accurate phase correction, which -unlike most other methodstakes into account information supplied in the form of signature models. We note that the recovered baseline r b = S (D (r)) 195 implicitly depends on all model parameters and observe that it is particularly sensitive to the linear phasing values, ϕ 0 and ϕ 1 .
We demonstrate this with a series of simulations. For this, we generate a spectrum of a single Lorentzian peak with slightly disturbed phase (see Fig. 3); this can, for example, correspond to a residual error after the usual phase correction procedure.
We fit this peak with a zero-phase signature model by minimizing Eq. 3 only with respect to the chemical shift and peak width thus intentionally keeping the phase error in the fit. As shown in Fig. 4, the least squares fit of an imperfectly phased peak leads 200 to erroneous estimates of its position and width. Consequently, phase error of approximately 0.4 rad is able to cause an error in the peak intensity estimate of about 5%. Thus, it is important to eliminate any phasing imperfections if the desired accuracy of quantification lies below 5%.
Although there is no random noise added to the spectrum in the above example, the phasing error manifests itself in the residual in Fig. 3. Notably, the asymmetry of an imperfectly phased peak induces positive and negative tails in the residual; 205 this effect is greatly emphasized by median filtering, which creates a sharp transition in r b . Naturally, we desire to recover the residual baseline as smooth as possible, and thus penalising such sharp edges is an especially effective strategy for fine-tuning of the phasing parameters. In this work, we found that the same criterion of Eq. 3 but now applied only to the extracted residual baseline works most effectively for this purpose, i.e. to adjust the phasing, we minimize min ϕ0,ϕ1

210
Furthermore, multistage median filtering with increasing window widths w m tends to improve the smoothing results, which agrees with recent works (Arias-Castro and Donoho, 2009). Intuitively, w m = 0 corresponds to no smoothing, and the problem of Eq. 8 reduces to the original least squares formulation of model fitting (Eq. 3), except for the removed noise. On the other hand, very broad filters on the order of the total spectral width tend to produce smooth results; in turn, this makes them ineffective for the above optimization, whose goal is to remove sharp edges by adjusting the phase. Notably, in our simulations, 215 the minimum is attained at the true value of the phasing parameter ϕ 0 using median filters at least eight times wider than the width of the peak (see Fig. 5). Thus, we propose to solve the problem of Eq. 8 iteratively: given an average peak width in the spectrum, FWHM, we start with a median filter of size at least w m > 10 · FWHM, minimize Eq. 8 with respect to the phasing parameters, and then increase w m to recover a smooth residual baseline r b at the final optimization stage.
We apply the same method to adjust both phasing parameters in the experimental data of thiamine in D 2 O. Fig. 6 displays 220 the values of the cost function in Eq. 8 plotted with respect to the deviation in the phasing parameters ϕ 0 and ϕ 1 from their current values. We note that the phasing of the spectrum after the usual least squares fit with Eq. 3 is suboptimal in the sense of our criterion based on the smoothness of the baseline. The proposed adjustment reduces the fitting cost (note the lower minima after adjustment, especially for ϕ 0 ) and significantly improves the phasing of the resulting spectrum as shown in   phasing. The proposed phase adjustment reduces peak-to-peak deviations in the residual baseline r b by more than four times, and further filtering with a wider window produces an almost flat baseline, which does not exceed the natural level of noise, as desired. This recovered baseline is now safely removed from the phased spectrum without affecting the areas under the peaks.
Unlike the least-squares criterion of Eq. 3, smoothing the residual baseline with Eq. 8 does not directly penalize the misfit 230 between the model and the data; in general it leads to slightly higher mean-squared error between x and y. However, it allows us to isolate the remaining unaccounted signal r m in the residual from the baseline and noise, which would have been naturally  Aldrich; the purity of both species is 99.8% by weight. 245 We use a Mettler Toledo AX205 balance with instrument accuracy of 0.1 mg (provided in the calibration protocol of the manufacturer). By means of the accuracy of the laboratory balance and error propagation, the uncertainty of the gravimetrically determined mole fraction was estimated to be 1.29 · 10 −5 mol/mol.
In this experiment, the data were acquired on a high-field NMR spectrometer with a 9.4 T vertical superconducting magnet (Ascend 400, console: Avance 3 HD 400, Bruker Biospin, Rheinstetten, Germany), which correspond to a proton Larmor points by zero-filling. The SNR in these datasets was estimated to be 100 − 10 4 depending on the specific peak considered.

Data Processing and Quantification
Peak integration and qGSD analysis were carried out with the software Mnova (version 14.0.1, Mestrelab Research, Santiago de Compostela, Spain). In each case, automatic phase (global, whitening) and baseline (Whittaker smoother or polynomial fit of the third degree) corrections were applied followed by visual inspection and manual adjustment where necessary. Integration boundaries for each peak are chosen based on their full width at half-maximum (FWHM) and are set at least 50 × FWHM.

265
Quantitative GSD was run with manual range selection and five improvement cycles.
The least-squares fitting and the proposed adjustment algorithm were implemented in a custom software written in Python 3.5.
For each sample s, we report the results of quantification with all methods in terms of the root mean square error (RMSE s ) in mole fractions computed with respect to the values obtained gravimetrically, x grav s,k : where x est s,k is the mole fraction of the k th species estimated in the s th sample and expressed in mol/mol. The average RMSE is computed over all S samples, RMSE avg = 1 S S s=1 RMSE s .
In this section, we apply the proposed adjustment procedure for model-based quantification to two sets of samples. First, we study the performance of our algorithm using a sample of thiamine in D 2 O and compare the relative ratios of its peaks with 275 the known ground truth. In the second example, we analyse a set of organic mixtures prepared gravimetrically. In both cases, we look at data acquired with a high-field spectrometer as well as a medium-field benchtop instrument.

A Sample of Thiamine in D 2 O
For the first series of experiments, we prepare a sample of 0.5 M thiamine dissolved in D 2 O, which we referred to previously in Section 2. We choose this compound for its very characteristic 1 H NMR spectrum: it exhibits a pair of well-separated peaks 280 at 5.45 and 7.9 ppm, a pair of partially overlapping peaks at 2.42 and 2.52 ppm, and a pair of triplets at 3.07 and 3.77 ppm (see Fig. 1). This allows us to test various aspects of our model, including the quantum mechanical formulation of coupled spin systems (see Matviychuk et al. (2019) for more detail). Since the true ratios of peak intensities are known and fixed, by estimating them separately -as if the peaks belonged to several chemical species in unknown concentrations -we can unequivocally compare different qNMR methods in terms of their accuracy. In these examples, we refer to and estimate the 285 intensity of the three peaks that comprise each triplet collectively.
We consider a spectrum acquired with a high field spectrometer, and also analyse the same sample with a benchtop system.
In the former case, the conventional peak integration readily achieves errors in relative mole fractions as low as 0.001 (see Fig. 9). Furthermore, quantitative global spectrum deconvolution (qGSD) performs excellently in this example and shows significant improvement compared to the standard GSD algorithm (Cobas et al., 2008;Cobas and Sýkora, 2009). The ratios of 290 peak intensities estimated with the least squares model fitting (LS) are also more accurate than the GSD results. However, due to inevitable lineshape misspecifications, the LS method can not outperform the peak integration of well separated peaks in a low noise spectrum as well as qGSD, which relies on a sophisticated deconvolution of each peak individually. On the other hand, the proposed adjustment algorithm significantly improves the LS results and brings the root mean-squared-error (RMSE) to the level achieved with peak integration (see Fig. 10).

295
It is instructive to look at the possible improvement of the least squares fitting results with more representative lineshape models. Specifically, we include second and third order terms in the real and complex-valued decay model of FID as done in (Matviychuk et al., 2017), which contribute additional weighting parameters to be fitted; we observe that the second order FIDs correspond to linear combinations of Lorentzian and Gaussian lines in the spectrum, while complex-valued polynomial decay models allow us to address peaks asymmetry. Furthermore, we consider a custom-written version of the reference 300 deconvolution method (Metz et al., 2000), in which we estimate the convolution kernel given a Lorentzian model fitted to the experimental data as described above. This method has the highest potential to represent various lineshape deviations but nevertheless is restricted by using the same convolution kernel for all peaks. The quantification results in these cases are summarized in Fig. 11 in terms of the RMS errors in peak ratios as well as the values of fitting objective of Eq. 3. As expected, more complex signal models allow us to fit the experimental data better and reduce the norm of the residual r. However, better Note that the adjustment method improves the accuracy of model-based quantification for all peaks in both cases.
LS fit does not always entail lower quantification errors, which signifies possible overfitting. On the other hand, the proposed phase adjustment by minimizing r b leads to slight increase of the total residual norm r , but the following distribution of r m among the reference signals according to Eq. 5 significantly reduces the quantification error.  Figure 10. Root mean square errors (RMSE) of estimating peak ratios in spectra of thiamine computed by averaging the results in Fig. 9.
The peak integration and the qGSD algorithm are very effective for the analysis of well resolved high-field data. The proposed adjustment method improves the LS quantification results even when peak overlap is present, such as in the medium-field spectra.
Analysis of the same sample acquired on a benchtop spectrometer is a more challenging task for the established methods.
Partially overlapping methyl peaks at 2.5 ppm make it difficult to define ranges for their integration. On the other hand, GSD 310 methods, that lack information about the underlying molecular system, often tend to define a third broad peak that overlaps with these two main resonances to compensate for lineshape broadening near the baseline. Although this produces a plausible fit overall, the spurious extra peak does not have any physical meaning and is difficult to account for in the final quantification results (we attribute its area to the closest resonance in our analysis). The distortions in lineshape of the triplets at 3-4 ppm are caused by higher-order coupling effects rather then a magnetic field inhomogeneity, and thus modelling them with a quantum 315 mechanical approach is especially effective here. Even though the RMSE of the LS model fit is more than two times lower than that of peak integration and qGSD, the proposed adjustment method allows us to further reduce it by almost 40%. The remaining error is mostly due to one of the overlapping methyl peaks, whose intensity estimate is particularly sensitive to small deviations in phasing parameters. As in the high-field example, the bottom panel in Fig. 11 shows the quantification accuracy achieved with alternative lineshape models. The proposed phase adjustment method with the distribution of the residual has 320 lower quantification error that the reference deconvolution approach and does not require fitting any additional lineshape parameters, unlike the the higher-order peak models.

A set of Organic Mixtures
Next, we study a set of 22 mixtures of methanol, ethanol, methyl acetate, and ethyl acetate prepared gravimetrically in varying relative concentrations ranging from 0.02 to 0.95 mol/mol for each component; as before, we measure their spectra with a   Figure 11. Results of LS fitting of high field (top) and medium field (bottom) data using models that take into account higher order lineshape distortions. From left to right: first order Lorenzian, second order symmetric model (weighted combination of Lorentzian and Gaussian lineshapes), third order symmetric and assymetric peak models. In the reference deconvolution approach, a single kernel is estimated for all peaks based on the difference between the experimental spectrum and a fitted Lorentzian model. For the proposed method, the results are reported after the phase adjustment (Eq. 8) and the final distribution of the residual (Eq. 6). Note that lower norms of the residual achieved by fitting more flexible models to the data do not always entail reduced quantification errors, whereas the proposed adjustment method achieves this goal more efficiently. Preprint. Discussion started: 28 January 2020 c Author(s) 2020. CC BY 4.0 License.
high-field and a medium-field benchtop spectrometer (see Fig. 12). Using these datasets, we estimate the mole fractions of each chemical species in the mixtures and compare them with the gravimetric values. To define the models for all chemical species, we use their complete quantum mechanical formulations and find the corresponding chemical shifts and J-couplings using the high-field data (Matviychuk et al., 2019). We use these parameters to initialise corresponding models at the lower field strength of the benchtop instrument and then refine them by minimizing Eq. 3. For comparison, we apply the qGSD algorithm to the high-field datasets: specifically, we include all non-overlapping peaks in our analysis and also peaks of proton in the acetate groups (see two overlapping peaks at 2.0 ppm in Fig. 12). Where possible, we also take into account the peaks of protons in the hydroxyl groups. On the other hand, severe peak overlap in the benchtop spectra precludes their accurate assignments and deconvolution, which makes this dataset extremely challenging to analyse with the traditional methods; therefore, methods based on fitting of quantum mechanical models are especially useful Preprint. Discussion started: 28 January 2020 c Author(s) 2020. CC BY 4.0 License. algorithm, and the proposed adjustment procedure reduces the average error in mole fractions by almost 50%. As expected, the analysis of the benchtop data results in slightly higher errors in mole fractions -with RMSE of up to 0.04 mol/mol for 340 certain samples, which nevertheless is acceptable in many practical applications. However, the average error across all samples, RMSE avg , is comparable to that achieved by the model-based methods (qGSD and LS) with the high field dataset, and the proposed model adjustment further reduces it by 25% on average. Finally, we note that occasionally -especially with the benchtop data -the adjustment of the LS fit results in slightly higher quantification errors (e.g. see the results for the seventh sample in Fig. 13). The increase in the error is likely due to an imperfect mechanism of distributing the residual among the 345 overlapping signature models. The assumption that the error is proportional to the component intensity as postulated in Eq. 7 may appear insufficient in these cases; the development of more accurate allocation rules is a topic of ongoing work. However, the increase is usually less than 0.005 parts in mole fractions, and if the entire data set is considered the average error is reduced, as already noted.

350
We proposed an effective and computationally simple mechanism to improve the accuracy of model-based quantification in NMR data analysis. The proposed adjustment procedure aims to account for all useful signal left in the residual after the usual least squares fit, which can signify a case of model misspecification -a problem notoriously difficult to avoid in most modelbased qNMR methods. Our alternative optimization criterion explicitly relies on the denoising of residual and smoothing the