Pulsed-field gradient nuclear magnetic resonance (PFG-NMR) is extensively used for characterization of macromolecules and macromolecular assemblies. The range of problems addressed by PFG NMR includes protein-ligand interactions (Tillett et al. 1999; Weljie et al. 2003; Lucas et al. 2004; Segev et al. 2005; Dalvit and Vulpetti 2012), protein oligomerization and aggregation (Ilyina et al. 1997; Price et al. 1999; Barhoum and Yethiraj 2010; Kheddo et al. 2016; Rabdano et al. 2017) including applications to amyloidogenesis (Tseng et al. 1999; Li et al. 2005; Narayanan and Reif 2005; Svane et al. 2008; Soong et al. 2009; Waelti et al. 2015; Hoffmann et al. 2018), protein folding (Wilkins et al. 1999; Balbach 2000; Buevich and Baum 2002; Li et al. 2007; Ramanujam et al. 2020), characterization of molecular machines and other supramolecular assemblies (Christodoulou et al. 2004; Baldwin et al. 2011; Huang et al. 2016; Leung et al. 2017; Kitevski-LeBlanc et al. 2018; Kharkov et al. 2021), probing special protein environments (Coffman et al. 1997; Wang et al. 2010; Waudby et al. 2012; Roos et al. 2016; Brady et al. 2017; Wong et al. 2020) and a host of other applications.

Oftentimes, PFG NMR measurements are faced with low signal-to-noise (S/N) ratio. This situation typically arises from low concentration of species of interest in the NMR sample. In turn, this can be a consequence of limited solubility or, otherwise, reflect the fact that sample material is difficult to prepare and/or expensive. Sometimes the situation is exacerbated by a short lifetime of the sample, which limits the experimental acquisition time. Along the same lines, when monitoring the evolution of amyloidogenic samples the data must be acquired over the short time intervals. Furthermore, such experiments often target certain specific species within the evolving mixture (e.g. monomers or low-order oligomers), which means that the actual concentration can be lower than the nominal concentration. Because of spectral overlaps, the observations may be limited to a small number of spectral lines (peaks), which further undermines the sensitivity. Alternatively, one can deliberately choose to work at low protein concentration in order to avoid the effects of obstruction (entanglement), weak dimerization effects, etc. (Wilkins et al. 1999; Liu et al. 2012; Ramanujam et al. 2020). This is particularly important to detect small changes in diffusion coefficient, e.g. due to expansion (contraction) of disordered proteins in response to experimental parameters such as temperature, pressure or pH.

The other factor responsible for low S/N ratio is associated with relaxation losses and relaxation broadening. This is particularly relevant for heavy (supra)molecular species. Bear in mind that signals from such species tend to be broadened beyond detection, such that often one has to rely on special labeling schemes in order to collect NMR data (Kitevski-LeBlanc et al. 2018). In turn, this means that only a limited number of spins (e.g. a subset of methyl sites) contribute to the observed signal, which has the effect of lowering the effective concentration of the sample. Alternatively, relaxation losses can be due to viscous environment, e.g. as a result of crowding, phase separation, etc. Since spin magnetization recovery is slowed down in the slow motion regime, the experiments have to employ long recycling delay, which further degrades S/N ratio. In some other cases, the signals come from a (moderately) flexible portions of a heavy system, e.g. histone tails in nucleosome core particle (Kato et al. 2009). In this situation, the widely used stimulated echo experiments can suffer from a certain amount of relaxation loss during the requisite long diffusion delay .

A substantial effort has been made to address the S/N problem in PFG NMR experiments by means of pulse sequence engineering. In doing so, the investigators used favorable relaxation properties of heteronuclear spin modes (Choy et al. 2002; Ferrage et al. 2003), selective excitation schemes and other designs to speed up magnetization recovery (Augustyniak et al. 2011; Augustyniak et al. 2012; Chan et al. 2015), more efficient phase coding via multiple-quantum coherences (Huang et al. 2017) and TROSY-type schemes to minimize relaxation losses (Horst et al. 2011; Didenko et al. 2011). However, to this day the majority of PFG NMR measurements are conducted on unlabeled samples using the simple and robust stimulated echo (STE) and double-stimulated echo (DSTE) sequences (Stejskal and Tanner 1965; Jerschow and Muller 1997). For these applications, poor S/N ratio often remains a pressing problem.

Here we address this problem by presenting a simple data processing scheme, which allows one to improve the precision of the standard NMR diffusion measurements. In what follows, we will briefly explain the idea behind the proposed scheme (a more detailed description is presented in Materials & Methods section; see also (Kharkov et al. 2021)). For the sake of argument, let us consider an NMR spectrum consisting of a single Lorentzian line that has been recorded with a very poor S/N ratio. In dealing with such data, the optimal processing strategy is to fit the spectrum with a Lorentzian function, integrate the fitted contour and thus obtain the (best attainable) estimate for signal intensity. The benefit to this approach is that we know the shape of the spectral line and utilize this prior knowledge in order to access the information of interest (i.e. to find the signal intensity). Clearly, this is a well-known strategy, which is broadly used in the experimental NMR practice. In particular, this strategy has been implemented in the programs NMRPipe (Delaglio et al. 1995) and FuDA (Hansen et al. 2007), where it is adapted for 2D and 3D spectra. Of note, the algorithm can deal not only with isolated peaks, but also with clusters of two or more peaks (which are treated as a superposition of several multidimensional Lorentzian contours or, depending on data acquisition scheme and window function applied, other types of contours). However, what happens if the spectral shape is complex and cannot easily be represented as a sum of several Lorentzians?

The case in question is a set of 1D spectra from STE or DSTE measurements on a biomacromolecular sample (here subscript enumerates the spectra collected with different gradient strengths ). In interpreting these data, one usually selects a certain region within the spectrum and determines the integral intensity of the signal over this region, . Typically, the selected region contains multiple (from tens to hundreds) heavily overlapped spectral lines, giving rise to a rather unique spectral envelope. This envelope cannot be readily or stably fitted using some sort of a generic template function, such as a combination of several Lorentizans. In lieu of such generic template, we seek to build a case-specific function that simply reproduces the relevant portion of the spectral envelope.

A simple way to obtain reasonable is to add all spectra in the STE (DSTE) series, . Obviously, this should reduce the amount of noise compared to each individual spectrum . However, on second thought it becomes clear that this strategy is sub-optimal. Indeed, the spectra corresponding to the strongest gradients usually contain precious little signal. Including these spectra into the sum does not improve , but rather contaminates it with noise. A better strategy is to form the weighted sum of the spectra, , where is the zeroth-order approximation for the integral of interest (obtained e.g. by straightforward integration of the -th spectrum over the selected region). This ansatz is rooted directly in the optimum filtration theory (Ernst et al. 1987), yielding the best possible S/N ratio for . Using the resulting “low-noise” template , we can further fit the individual spectra – and thus obtain more accurate intensities, (similar to the familiar fitting strategy using Lorentizan templates).

This new data processing scheme has been implemented as a part of the web server DDfit (Diffusion Data fit, https://ddfit.bio-nmr.spbu.ru/). The server allows one to load the data files in a standard Bruker format, indicate the boundaries of the integration region and in a matter of several seconds obtain the table of . If desired, the server can also fit the results using the famous Stejskal-Tanner equation for the STE experiment (Stejskal and Tanner 1965) or its suitably modified version for the DSTE experiment (Jerschow and Muller 1997), produce the graph with the fitted data and report the determined diffusion coefficient. Simple interfaces are available to handle the results of the experiments recorded with the standard Bruker pulse sequences; otherwise, more flexible interfaces are offered to deal with user-designed sequences.

The performance of DDfit has been tested on both simulated and experimental diffusion data. The new scheme is shown to be both accurate and precise, achieving considerably better precision than the standard processing schemes such as offered, for example, in the Topspin package (Zick 2016) or state-of-the-art MestreNova program (Mestrelab 2024).

The weighting scheme is based on the idea of matched filtration, which maximizes the signal-to-noise ratio of the resultant sum spectrum (Peebles 2000; Ernst et al. 1987).

(iv) All spectra are approximated using the ansatz . The term is intended to refine the baseline correction performed in step (i). The term represents the spectrum per se, assuming that each spectrum can be expressed as a scaled version of the model spectrum. The coefficients , , and are fitted by minimizing the following target function:

In doing so, standard linear least-squares procedures, such as e.g. numpy.linalg.lstsq, can be used. The coefficients are taken to be the experimental values of , i.e. the refined signal intensities that are subsequently used to determine the diffusion coefficient.

The schematic layout of the DDfit server (https://ddfit.bio-nmr.spbu.ru/) is shown in Fig. 1. Currently, it is assumed that diffusion data are acquired on a Bruker spectrometer and stored in Bruker format (in principle, routines exist to convert Varian or JEOL data to Bruker format (Bruker 2020; Helmus and Jaroniec 2013) before processing them with DDfit, but this would likely require some manual editing of the spectrometer files). As a first step, a series of FIDs from the diffusion experiment are pre-processed using Bruker’s TopSpin software, viz. multiplied by a window function, Fourier-transformed and phased. The resulting spectral data are transferred to DDfit in a form of standard 2rr file. Also transmitted to DDfit are procs file (contains several parameters of interest, such as spectrometer operating frequency and spectral width), proc2s file (supports reading of the 2rr matrix by the nmrglue accessory (Helmus and Jaroniec 2013)), and difflist file (contains a list of gradient strengths expressed in the units of G/cm). The gradient strengths in difflist are for rectangular gradient shapes (in the case of shaped gradients, the shape factor is factored in ).

* For DSTE experiments, the results are interpreted using Stejskal-Tanner equation with Jerschow-Muller modifications.

Figure 1. Block diagram illustrating the processing of NMR diffusion data by TopSpin and DDfit.

The main page of the DDfit website offers user a choice between STE and DSTE experiments; besides, an option to determine signal intensities irrespective of the experimental details is also available. For both STE and DSTE, the user can select one of the standard Bruker experiments: stegp1s for STE using simple monopolar gradients, stebpgp1s for STE using bipolar gradients, dstegp3s for DSTE using monopolar gradients, or dstebpgp3s for DSTE using bipolar gradients. Otherwise, non-standard (in-house) variants of the same pulse sequences can also be processed.

For any specific choice of the pulse sequence, the user is presented with a small web form where he/she needs to enter several parameters pertaining to this particular experimental setup. For example, in the case of dstegp3s experiment one has to indicate the observation spin (selected from a dropdown menu, typically ¹H), the length of the gradient pulse p30, the duration of the diffusion delay d20, the length of the hard 90° pulse p1 and the duration of the gradient recovery delay d16. Likewise, for the other three Bruker sequences the user must enter several relevant parameters specific to the chosen experiment.

In the case of non-standard STE or DSTE sequences, the user is presented with a generalized scheme of the sequence which contains generic names for pulses and delays. The idea is that the experimentalist should make a connection between the generic variables and the actual parameters of the experiment. The so-defined variables are then entered into the web form associated with this particular type of pulse sequence.

Finally, the user also needs to indicate the boundaries of the spectral region to be used for signal intensity determination, and . The user can choose to enter these values in the units of ppm or otherwise in points. The criteria to choose a suitable spectral region are described in the previous section.

Once all input parameters are entered, the job is submitted for execution on a dedicated DDFit server. As a first step, the series of spectra are fitted using the optimal model as described in the previous section. If user request is to evaluate signal intensities (“integrals only” panel), the program generates a table of intensities, , and stops there. The table can be displayed for visual inspection or downloaded as an ascii file. As a part of the output, the server also generates an image of the spectrum (the one acquired under the lowest gradient strength), where the integration region is highlighted in red. The server also produces another image, showing a stack of experimental spectra (following baseline correction, restricted to the integration region). In addition, it offers a basic interactive plotting facility, allowing one to zoom into selected regions of the spectrum, adjust the plotting scale, etc. All images can be conveniently saved in various graphics formats. The user can generate a unique link to the results page and store it for future reference.

Otherwise, if user request is to analyze the diffusion data (“stimulated echo” or “double-stimulated echo” panels), the program additionally fits the obtained data to Stejskal-Tanner equation or Stejskal-Tanner equation with Jerschow-Muller modifications. In this case, DDfit additionally displays a graph illustrating the data fitting and reports the determined diffusion coefficient together with the standard error derived from the covariance matrix. This part of DDfit service is trivial and uses a standard minimization routine. All calculations, including spectral intensity determination and Stejskal-Tanner fitting, take no more than several seconds.

DDfit also offers a demo page, demonstrating the analyses of DSTE data acquired in-house for a sample of hen egg-white lysozyme (see next section). All spectrometer files in the demo package are downloadable. Finally, there is also a help page, providing the description of the algorithm, usage instructions, and a review of stimulated echo sequences.

The DDfit server has been built using Python Flask framework, with Nginx as a frontend and uWSGI as a gateway interface. The application makes use of several libraries: nmrglue (reading of spectrometer data), scipy (minimization), matplotlib (graphics), Bokeh (interactive graphics) and uuid module in Python (generation of unique webpage identifiers). The server is hosted on a desktop PC equipped with Intel i5-4440 CPU and 16 GB of RAM running Ubuntu 24.04.3.

DDfit algorithm was first tested on synthetic diffusion data. Synthetic FIDs have been generated for ubiquitin using the simulation platform Spinach 2.6.5625 (Edwards et al. 2014). The list of proton chemical shifts was from the Spinach demo file 1d3z.bmrb. The spectrometer operating frequency was assumed to be 500 MHz and the spectral width was taken to be 15 ppm. For simplicity, we turned off spin relaxation in Spinach simulations; instead, a uniform line broadening was introduced by applying the exponential multiply window function to the simulated FID (LB of 15 Hz). To simulate the effect of translational diffusion, 15 copies of ubiquitin FID have been multiplied by Stejskal-Tanner factors:

where is a proton gyromagnetic ratio, is the diffusion delay (100 ms), is the duration of rectangular gradient pulse (5.4 ms), is a series of fifteen gradient strengths uniformly distributed between 0.95 and 35.1 G/cm and is the translational coefficient of ubiquitin taken to be 11.2×10^− 7 cm² s^− 1 (Ramanujam et al. 2020).

Separately, we have generated an FID of solvent water resonating at 4.7 ppm. The broadening LB of the water signal was set to 2 Hz and the diffusion coefficient was taken to be 181×10^− 7 cm² s^− 1 corresponding to HDO diffusion in D₂O (Stefaniuk et al. 2022), see below. The protein and water FIDs were then combined assuming that the sample was prepared with 100 µM ubiquitin in D₂O solvent with 2% residual H₂O content. In doing so, we assumed that highly labile protons in ubiquitin were fully displaced with deuterons while amide protons were retained (the latter assumption is inconsequential since only the aliphatic portion of the simulated spectrum is further used for signal evaluation). The synthetic diffusion data were thus prepared in a form of 2D complex matrix with dimensions 15×4096.

As a next step, we have generated =16,383 complex noise matrices of the same size, 15×4096. The matrices contained uncorrelated Gaussian noise with the zero mean; this is a sufficiently accurate representation of spectrometer noise associated with quadrature detection (Grage and Akke 2003). The magnitude of the noise was set to a small fraction of the protein FID amplitude, leading to the average S/N of 65 in the noised 1D spectra in the absence of gradient or otherwise 63 under the weakest gradient (where S/N is defined as the ratio of the maximum amplitude of the protein signal in the spectrum and the root-mean-square amplitude of the noise). The so-obtained realizations of the noised gradient-encoded diffusion data were then processed using the DDfit algorithm as further detailed in the Results section.

DDfit algorithm has also been tested on experimental diffusion data acquired in-house from a sample of 0.5 mM hen egg-white lysozyme (HEWL, Sigma-Aldrich). The protein was dissolved in D₂O-based solvent (0.01% NaN₃, 5% dioxane, pH 2.0) and the data were collected at 15°C; under these conditions the sample remains folded and stable (Mariño et al. 2015). The measurements were conducted on Bruker Avance III 500 MHz spectrometer equipped with BBI probe using double-stimulated echo pulse sequence (dstebpgp3s). The diffusion delay was 200 ms and the length of bipolar gradient pulse was 5 ms (two pulses with an opposite polarity of 2.5 ms each, smoothed square shape SMSQ10.100, integral shape factor 0.9). The data were recorded with twenty gradient strengths , incremented linearly from 0.95 to 46.53 G/cm; only twelve points, from 0.95 to 28.01 G/cm have been used in data analyses. When recording the data, the gradients were arranged in pseudo-random order. The number of scans was 16 (a minimum number to accommodate the phase cycle). In total, 77 repeat diffusion experiments have been recorded back-to-back in 64 hours.

We begin by exploring different data processing schemes using the simulated STE dataset as described in Materials and Methods. In brief, this is the series of gradient-encoded 1D proton spectra, representing a sample of ubiquitin in D₂O with residual content of H₂O. The methyl portion of the spectra, from = 1.33 ppm to = 0.30 ppm, is used for signal evaluation. This region offers a highest amount of protein signal while it is also maximally distant from the water resonance. The analyses are conducted for = 16,383 replicas of the simulated STE dataset containing a moderate amount of random Gaussian noise. The data were pre-processed using EM window function with line broadening of 15 Hz (matched to the simulated homogeneous line broadening). The average S/N ratio of the simulated spectra within the chosen region () under the weakest gradient is 63.

As a first step, we test the simplified processing scheme where we do not perform any baseline correction, but rather directly proceed to spectra integration. The treatment, which involves numeric integration of the spectra within () interval followed by Stejskal-Tanner fitting of the so-obtained integrals, has been repeated times in a Monte-Carlo style simulation. The fitted values of diffusion coefficient for ubiquitin, , are summarized in Fig. 2a. We further selected the dataset corresponding to the median of the obtained distribution; for this particular dataset we present the stack of the simulated STE spectra plotted in the interval from to , Fig. 2b, and the results of Stejskal-Tanner fitting, Fig. 2c.

The inspection of Fig. 2a immediately indicates that values obtained from this procedure are systematically overestimated compared to the target value encoded in the simulated data (the histogram is shifted to the right from the vertical line, which indicates the target). The reason for this apparent bias is quite obvious since the integrated signal contains a contribution from the (extending far upfield) wing of the residual water signal. As a consequence, the diffusion decay profile in Fig. 2c consists of the two components: the dominant contribution from the slowly diffusing protein and a minor contribution from the rapidly diffusing water (affecting mainly the initial points of the diffusion decay profile). Predictably, fitting this combination with the Stejskal-Tanner equation equipped with the single parameter results in overestimation of the diffusion coefficient.

Specifically, the simulations yield (11.73 ± 0.17)×10^− 7 cm² s^− 1 (here and in what follows we report the average value together with the standard deviation from histogram such as illustrated in Fig. 2a). This is 5% overestimated relative to the target value of 11.20×10^− 7 cm² s^− 1. Given that diffusion coefficient is inversely proportional to the cubic root of molecular mass, , this effect corresponds to a roughly 15% loss in the estimated mass of the protein. Under certain circumstances the error on this scale can lead to erroneous conclusions. For example, it may be misinterpreted as the evidence of protein compaction, viz. transition from a strongly disordered state to a molten globule.

In summary, the most basic data processing strategy, which focuses on a methyl region far away from the water signal and which neglects to perform any baseline correction, yields the results that are highly reproducible (i.e. precise), but at the same time are biased (i.e. inaccurate).

As a next step, we have implemented the protocol involving the simplest baseline correction scheme – namely, 2-point baseline correction. Specifically, the baseline is approximated by a straight line connecting a pair of points, and . After the base is subtracted, the signal is numerically integrated same as before.

The outcome of the Monte-Carlo analyses for this protocol is illustrated in Figs. 2d-f, showing the histogram, the stack of simulated STE spectra corresponding to the median value, and the Stejskal-Tanner fit for this particular dataset. Obviously, application of the 2-point baseline correction scheme eliminates bias in determination of the diffusion coefficient – the histogram is now centered at the vertical line which marks the target value of . At the same time the width of the distribution is significantly increased, indicating greater uncertainty. The mean value of now amounts to (11.19 ± 0.57)×10^− 7 cm² s^− 1, which is exactly on target, but with three-fold increased uncertainty range.

These observations can be readily rationalized. Indeed, 2-point baseline correction takes out the solvent contribution to the spectrum (along with a portion of protein signal), thus eliminating the bias toward fast diffusion. It is easy to see that this scheme works well only when the solvent spectrum within the region of interest, from to , can be accurately approximated by a straight line. This is indeed the case for the wing of the water signal since the selected region (from 1.33 to 0.30 ppm) corresponds to a relatively narrow band far removed from the water resonance (4.7 ppm). Generally, this strategy should work best when choosing a narrow integration region far away from strong solvent signals.

At the same time, 2-point baseline correction scheme suffers a setback in that the magnitude of correction strongly depends on random noise contained in the two spectral points, and . In other words, there is no mechanism to average out noise in this simple procedure (more on this in what follows). In turn, this leads to significant random variations in obtained values.

As a result, the improved data processing strategy, employing 2-point baseline correction scheme, produces the results that are accurate, but not very precise.

Finally, we illustrate the performance of the scheme at the core of the DDfit algorithm, Figs. 2g-i. Same as before, the treatment begins with a simple 2-point baseline correction followed by the calculation of numeric integrals, . After that the procedure takes a different course. The integrals are used as weights in calculating of the model spectrum with optimal signal-to-noise properties, see Eq. (1). The construct is then used to fit the individual spectra, , see Eq. (2). In this construct, the term absorbs the noise-induced error from the 2-point baseline correction, while the coefficient captures the overall intensity of the protein signal, equivalent to the refined integral .

As can be seen from Fig. 2g, this method is bias-free and at the same time relatively insensitive to random noise in the simulated data. The mean value of , (11.19 ± 0.27)×10^− 7 cm² s^− 1, perfectly well reproduces the target, 11.20×10^− 7 cm² s^− 1, with two-fold decrease in uncertainty compared to the previously discussed scheme, cf. Figure 2d.

In Fig. 2h, we focus on the simulated dataset corresponding to the median value. Shown in the panel is the model spectrum (thick red line), as well as the series of the fitted spectra (thin colored lines). The Stejskal-Tanner fit of the so-obtained integrals, Fig. 2i, shows only a very limited amount of scatter, significantly less than previously observed in Fig. 2f.

In conclusion, in our simulated experiments the DDfit algorithm proves to be highly accurate and, at the same time, reasonably precise.

At this point one may ask if a more sophisticated baseline correction scheme could improve the situation and, perhaps, obviate the need in the DDfit scheme. As it appears, the short answer to this question is “no”. In fact, some highly effective baseline correction schemes have been developed for NMR studies of small molecules and their mixtures – particularly, in the context of ¹³C spectroscopy (Carlos Cobas et al. 2006; Xi and Rocke 2008). In this latter case, however, the spectrum consists of a limited number of sharp signals with little or no interference from solvent; the baseline is well sampled throughout the spectrum and typically has a simple shape (characteristic baseline roll) which can be readily captured as a part of the baseline removal procedure. In contrast, ¹H spectra of proteins are cluttered with many (significantly broadened) spectral lines. The baseline is dominated by a residual water signal, which extends far up- and downfield from the H₂O resonance frequency, has rather irregular shape and is gradient-dependent. Other major components of the buffer solution can also make a contribution to the baseline. The baseline is well sampled only at the edges of the spectrum, but cannot be directly traced between 0 and 9 ppm. All of this makes it essentially impossible to develop an effective baseline removal scheme for use in proton-based PFG NMR experiments on protein samples.

To demonstrate this point, we have tested the spline baseline correction scheme, which is one of the few more advanced options in the Bruker TopSpin repertoire. On some occasions this scheme has been actually used to process PFG NMR data from protein samples (Augé et al. 2009; Demers and Mittermaier 2009). The results of our tests are summarized in Fig. S1. In the situation where the baseline cannot be sampled well (discussed above), splines prove to be highly sensitive to noise, i.e. unstable. As a result, spline-based procedure has an effect of amplifying noise, yielding mean of (11.46 ± 2.70)×10^− 7 cm² s^− 1. This is a staggering ten-fold jump in uncertainty compared to the DDfit treatment. We have also observed similar unsatisfactory performance when using other advanced baseline correction methods, such as Whittaker smoother algorithm (Cobas 2018) available in MestreNova (Willcott 2009) (results not shown). The fundamental reason for this lack of success is the nature of protein ¹H PFG NMR spectra, as described above.

To test the performance of DDfit relative to other data processing schemes on experimental data we have recorded 77 back-to-back DSTE experiments on a sample of 0.5 mM hen egg-white lysozyme (see Materials and Methods). Two spectral regions have been selected for signal integration: region I, extending from 2.23 to -0.05 ppm, which is the portion of the spectrum with the highest signal intensity (S/N ratio of 109), and region II, extending from − 0.05 to -0.34 ppm, featuring four partially overlapped upfield-shifted methyl resonances belonging to residues Ile-98, Met-105, Leu-8 and Leu-17 (S/N ratio of 34) (Dobson et al. 1984; Redfield and Dobson 1988; Maeno et al. 2009).

Figure 3 illustrates three different signal-processing schemes as applied to the experimental datasets, spectral region I. The treatment neglecting baseline correction is clearly a failure, resulting in distorted profile, see Fig. 3c. As it happens, the data acquired at lower gradient strengths, i.e. in the presence of sizeable residual water signal, suffer from a persistent artefact in a form of lowered baseline (noticeable toward the right edge of Fig. 3b). This leads to underestimation of the respective integrals, causing characteristic distortions in the diffusion profile. The effect is responsible for the strong bias in the obtained values, see Fig. 3a, as well as their poor reproducibility, = (3.49 ± 0.47)×10^− 7 cm² s^− 1. While in the case of (idealized) simulated data one can forego the baseline correction step with relatively little consequence, cf. Figure 2a-c, this omission becomes a major issue when dealing with the real-life data.

The situation is rescued by using a simple 2-point baseline correction scheme. The baseline-corrected spectra no longer suffer from downward (or upward) shifts, see Fig. 3e, and the profiles, Fig. 3f, are free from the previously observed systematic distortions. The values obtained from 77 back-to-back experiments produce the expected Gaussian-shaped distribution, Fig. 3d, with the mean value of (7.64 ± 0.21)×10^− 7 cm² s^− 1. This result is compatible with the available literature-based estimates. Specifically, STE results from Arata laboratory (Price et al. 1999) (288 K, pH 3, 1.5 mM lysozyme in 90% H₂O – 10% D₂O solvent) after correction for solvent viscosity and protein concentration yield ≈ 6.7×10⁻⁷ cm² s⁻¹, whereas similar STE results from Byrd laboratory (Altieri et al. 1995) (298 K, pH 2.3, 2 mM lysozyme in D₂O-based solvent) after correction for temperature and protein concentration yield ≈ 8.3×10⁻⁷ cm² s⁻¹. The two estimates, therefore, straddle the result obtained in this study.

The situation is further improved by using our DDfit processing scheme, Fig. 3g-i. In this case, the values determined from seventy-seven consecutive DSTE experiments prove to be highly reproducible, with the average of (7.65 ± 0.05)×10^− 7 cm² s^− 1. Thus, the use of DDfit scheme leads to a 4-fold reduction in error margin compared to the conventional processing strategy, i.e. 2-point baseline correction followed by signal integration. Reiterating our conclusion from the analyses of the simulated data, we observe that DDfit algorithm leads to both accurate and precise results also in the case of experimentally measured data.

As a sidenote, we also tested an alternative baseline correction scheme. Same as in the previous section, we used splines for baseline correction followed by direct integration of the signal in the () interval. While spline-approximated baseline for a given (high-intensity) spectrum is visually reasonable, see Fig. S2, the scheme is unstable and, as already pointed out, leads to noise amplification. The mean value from this method is (7.72 ± 0.60)×10^− 7 cm² s^− 1, signifying 12-fold increase in the magnitude of error compared to the DDfit scheme.

Finally, we have performed the same series of tests focusing on region II in the experimental spectra, where the signal is 25-fold weaker than in region I (referring to the integral intensities after 2-point baseline correction). The results are summarized in Fig. S3. In brief, the treatment neglecting baseline correction is a complete failure (compromised by variable baseline shifts as discussed above). The problem goes away when baseline correction is applied. The scheme using simple 2-point baseline correction followed by direct integration yields = (7.52 ± 0.69)×10^− 7 cm² s^− 1, whereas the DDfit treatment using 2-point baseline correction followed by model-based integration yields = (7.43 ± 0.29)×10^− 7 cm² s^− 1. The results are consistent with each other and also with those previously obtained from region I (see above). Once again, the precision of the DDfit algorithm is considerably better than that of the standard data processing scheme, with 2.4-fold improvement in the confidence interval.

Model-based integration of spectral signals is generally a powerful approach, utilizing certain prior information about the characteristics of these signals. For example, the popular NMRPipe package (Delaglio et al. 1995) includes options to treat a stack of 2D spectra, assuming that a number of parameters for a selected spectral peak are a given and are reproduced across the series of spectra (e.g. ¹H^N chemical shift, ¹⁵N chemical shift, and the peak shape). Furthermore, certain assumptions can usually be made about the peak shape, e.g. in solution-state spectra the shapes are typically Lorentzian or, otherwise, can be described by a convolution of Lorentzian contour with the Fourier transform of the (known) window function (Van Horn et al. 2010). NMRPipe routine that utilizes such prior information has been widely used to process data from relaxation experiments, J-coupling and residual dipolar coupling measurements, exchange measurements and other types of experiments (Xue et al. 2014; Wang et al. 2021; Yuwen et al. 2018; Maltsev et al. 2012; Ward and Skrynnikov 2012; Chevelkov et al. 2010; Wylie et al. 2011).

The distinction of the protein diffusion experiments, as discussed in this paper, is that there is no a priori model for the spectral signal of interest (such as Lorentzian function for an isolated spectral peak). Instead, we deal with a sample-specific spectral pattern , which consists of many overlapped spectral lines that fall within the region of interest, (). To address this problem, we have developed a recipe to calculate as an intensity-weighted sum of the gradient-encoded spectra. This recipe is rooted in the optimal filtration theory and ensures the best quality of the model. Armed with , the algorithm performs a model-based integration of spectra, arriving at uniquely precise results.

The obvious assumption here is that the spectral shape remains unchanged, up to a scaling factor, across the series of diffusion spectra . This is certainly justified for samples that contain only one sort of protein species. This assumption, however, becomes questionable for mixtures (e.g. fully developed monomer-dimer equilibrium) and ought to be tested before our protocol is used.

One should also bear in mind that any processing scheme for PFG NMR diffusion data should necessarily include the baseline correction step. We have investigated this aspect in some detail and found that the simple 2-point baseline correction works best for both standard processing schemes and our new algorithm. In contrast, more sophisticated variants, such as spline-based corrections, lead to significant noise amplification.

With these caveats in mind, our algorithm offers an efficient tool to analyze diffusion data from proteins and other large biomolecules. The method has been implemented in a form of web server DDFit, accessible at https://ddfit.bio-nmr.spbu.ru/. The server is easy to use for standard Bruker-based experiments as well as custom-built sequences, with calculations taking no more than several seconds. We anticipate that this user-friendly web application, which determines protein diffusion coefficients with much better precision than other existing programmatic solutions, will become a useful addition to the ever-expanding bio-NMR toolchest.