DeerLab: a comprehensive software package for analyzing dipolar electron paramagnetic resonance spectroscopy data

Dipolar electron paramagnetic resonance (EPR) spectroscopy (DEER and other techniques) enables the structural characterization of macromolecular and biological systems by measurement of distance distributions between unpaired electrons on a nanometer scale. The inference of these distributions from the measured signals is challenging due to the ill-posed nature of the inverse problem. Existing analysis tools are scattered over several applications with specialized graphical user interfaces. This renders comparison, reproducibility, and method development difficult. To remedy this situation, we present DeerLab, an open-source software package for analyzing dipolar EPR data that is modular and implements a wide range of methods. We show that DeerLab can perform one-step analysis based on separable non-linear least squares, fit dipolar multi-pathway models to multi-pulse DEER data, run global analysis with non-parametric distributions, and use a bootstrapping approach to fully quantify the uncertainty in the analysis.

1 Introduction

Dipolar electron paramagnetic resonance (EPR) spectroscopy encompasses a growing family of techniques for determining distributions of nanometer-scale distances between unpaired electrons. These distance distributions provide valuable information for the structural characterization of macromolecular or biological systems that is complementary to information obtained by other techniques. For structurally disordered or highly complex systems, where established techniques may fail, such distance distributions provide unique information. The family of dipolar EPR spectroscopy techniques includes double electron–electron resonance (DEER) (Milov et al., 1981, 1984; Pannier et al., 2000a, b), double quantum coherence (DQC) (Saxena and Freed, 1996, 1997; Borbat et al., 2013), relaxation-induced dipolar modulation enhancement (RIDME) (Kulik et al., 2001; Milikisyants et al., 2009), single-frequency technique for refocusing (SIFTER) (Jeschke et al., 2000), and several other related techniques (Borbat et al., 2013; Di Valentin et al., 2014; Hintze et al., 2016; Pribitzer et al., 2017; Borbat and Freed, 2017; Doll and Jeschke, 2017; Milikisiyants et al., 2018). All of them provide a time-domain signal that depends on the dipolar interaction between pairs of electrons. From this time-domain signal, the distance distribution is inferred.

Due to its mathematical nature, the robust inference of distance distributions from noisy dipolar EPR spectroscopy data is not straightforward. Many approaches have been proposed to tackle this problem (Pannier et al., 2000a; Jeschke et al., 2002, 2004, 2006; Bowman et al., 2004; Chiang et al., 2005a, b; Sen et al., 2007; Brandon et al., 2012; Stein et al., 2015; Dzuba, 2016; Matveeva et al., 2017; Srivastava and Freed, 2017; Edwards and Stoll, 2016; Rein et al., 2018; Timofeev et al., 2018; Worswick et al., 2018; Hustedt et al., 2018; Edwards and Stoll, 2018; Fábregas Ibáñez and Jeschke, 2019, 2020; Sweger et al., 2020), each with its pros and cons. Some of these methods have found widespread use, via software packages such as DeerAnalysis (Jeschke et al., 2006), GLADD/DD (Brandon et al., 2012), and LongDistances (Altenbach, 2020).

However, there are several major challenges with the current situation. (i) A comparative assessment of the relative merits of various methods is missing. (ii) Many methods have been argued based on anecdotal evidence from small datasets, and their performance has not been assessed comprehensively. (iii) Reproducibility of analysis results is very limited due to the lack of common platforms for data sharing and data analysis.

To remedy this situation, we introduce DeerLab, an open-source software for data analysis in dipolar EPR. It is based on the Python programming language and consists of a collection of modular functions, analogous to EasySpin (Stoll and Schweiger, 2006) and Spinach (Hogben et al., 2011). This has several distinct advantages over a graphical user interface (GUI). (i) It allows for very flexible workflow designs, easily adapting to different experimental situations. (ii) All existing methods can be directly compared on a single platform. (iii) Automation and processing of large datasets become straightforward. (iv) Scripted data analysis improves reproducibility and collaboration. (v) It provides a foundation for implementing new methodologies. (vi) It can be embedded into other software, such as tools for protein structure modeling based on distance distributions (Jeschke, 2018). The disadvantage is an accessibility barrier for potential users without programming skills. This disadvantage can be remedied by building a dedicated GUI for standard workflows as an additional software layer.

This paper is structured as follows. We start by summarizing the theoretical basics of dipolar EPR spectroscopy. Then, we illustrate the functionality of DeerLab through a series of examples. First, we show how to reproduce well-established workflows such as Tikhonov regularization and multi-Gauss fits. We then demonstrate how DeerLab can perform one-step analysis based on separable non-linear least-squares optimization, fit dipolar multi-pathway models to multi-pulse DEER data, and run global analysis with non-parametric distributions. Finally, several sections are dedicated to important aspects of uncertainty analysis, including the use of bootstrapping. The DeerLab scripts for generating all figures, as well as the corresponding distance distributions and dipolar signals, are available in the supporting information.

2 Theoretical basics

This section summarizes the central theoretical concepts of dipolar EPR spectroscopy that DeerLab is based on. For more details, see Jeschke (2012, 2016). The theory is limited to $S=\mathrm{1}/\mathrm{2}$ spins with isotropic g-values, without any orientation selection, at most two spins per protein, no exchange coupling, weak dipolar coupling, and no conformer-dependent relaxation rates.

Dipolar EPR spectroscopy techniques measure the magnetic dipole–dipole couplings between spins via the modulation of the amplitude V(t) of a spin echo as a function of the position t of one or more pump pulses. The echo amplitude is modeled as (Milov et al., 1981)

$$\begin{array}{}\text{(1)}& V\left(t\right)=V_{\mathrm{0}}\cdot V_{\mathrm{intra}}\left(t\right)\cdot V_{\mathrm{inter}}\left(t\right),\end{array}$$

where V₀ is the echo amplitude in the absence of any pump pulses, V_intra describes the pump-pulse-induced modulation due intra-molecular dipolar couplings, and V_inter describes the modulation due to intermolecular couplings. V₀ is a constant prefactor that we set to one from now on in order not to complicate the notation unnecessarily. DeerLab takes V₀ into account as a fitting parameter.

The product of intra- and inter-molecular dipolar modulations can be written in a general form as

$$\begin{array}{}\text{(2)}& V\left(t\right)=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}K(t,r)P\left(r\right)\mathrm{d}r.\end{array}$$

P(r) is the distribution of intra-molecular spin–spin distances r on the protein or other complex, normalized such that ${\int }_{\mathrm{0}}^{\mathrm{\infty }}P\left(r\right)\mathrm{d}r=\mathrm{1}$.

K(t,r) is the kernel that captures how the complete dipolar modulation is determined by the distance distribution. It includes the inter-molecular modulation (Fábregas Ibáñez and Jeschke, 2020). For standard four-pulse DEER it is

$$\begin{array}{}\text{(3)}& K(t,r)=\left[(\mathrm{1}-\mathit{\lambda })+\mathit{\lambda }{K}_{\mathrm{0}}(t,r)\right]V_{\mathrm{inter}}(t,\mathit{\lambda }).\end{array}$$

Here, λ is the modulation depth. K₀ is the elementary kernel

$$\begin{array}{}\text{(4)}& K_{\mathrm{0}}(t,r)=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\cos \left[(\mathrm{1}-\mathrm{3}{\cos }^{\mathrm{2}}\mathit{\theta })D{r}^{-\mathrm{3}}t\right]\mathrm{d}\cos \mathit{\theta },\end{array}$$

with the dipolar coupling constant

$$\begin{array}{}\text{(5)}& D={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{g_{\mathrm{e}}^{\mathrm{2}}{\mathit{\mu }}_{\mathrm{B}}^{\mathrm{2}}}{\mathrm{\hslash }}},\end{array}$$

where g_e is the g-value of the free electron, μ_B the Bohr magneton, μ₀ the magnetic constant, and ℏ the reduced Planck constant. K₀ assumes full orientation averaging and unlimited excitation bandwidth. The subscript 0 distinguishes this elementary kernel from more general kernels such as Eq. (3).

V_inter(t,λ) represents the inter-molecular modulation and is commonly called the background. It can be modeled as a stretched-exponential function

$$\begin{array}{}\text{(6)}& V_{\mathrm{inter}}(t,\mathit{\lambda })=\exp \left(-{\mathit{\kappa }}_{\mathrm{d}}\mathit{\lambda }\mathrm{|}t{\mathrm{|}}^{d/\mathrm{3}}\right),\end{array}$$

where κ_d is a decay rate constant and d is the dimensionality (Kutsovsky et al., 1990; Milov et al., 1998; Jeschke et al., 2002). Other background models are possible (Kattnig et al., 2013).

Experimentally, the echo amplitude is measured only for a discrete set of usually equally spaced time points t_i, yielding a dipolar signal vector V with n elements V_i=V(t_i). For numerical analysis, P(r) is represented as a discrete distance distribution vector P with m elements P_j=P(r_j) at equally spaced r_j. With this, Eq. (2) reads

$$\begin{array}{}\text{(7)}& \mathit{V}=\mathbf{K}\mathit{P},\end{array}$$

where K is the n×m kernel matrix with elements $(\mathbf{K}{)}_{ij}=K(t_i,r_j)\mathrm{\Delta }r$, and Δr is the increment in the distance domain.

Experimental data V_exp deviate from in Eq. (1) due to presence of noise. From experiments, it was found that the noise distribution in DEER signals is well approximated by an uncorrelated Gaussian distribution with zero mean and constant variance (Edwards and Stoll, 2016):

$$\begin{array}{}\text{(8)}& {\mathit{V}}_{\exp }=\mathit{V}+\mathcal{N}\left(\mathbf{0},{\mathit{\sigma }}^{\mathrm{2}}\mathbf{I}\right).\end{array}$$

Inferring the distance distribution from the dipolar signal formally requires inversion of the kernel matrix

$$\begin{array}{}\text{(9)}& \mathit{P}={\mathbf{K}}^{-\mathrm{1}}{\mathit{V}}_{\exp }.\end{array}$$

However, K is generally badly ill-conditioned (it has an extremely large condition number). This renders the inverse problem ill-posed, and the results obtained by Eq. (9) are highly unstable, erratic, and unreliable, especially in the presence of noise. Because of the ill-posedness, inferring distance distributions from the dipolar signals poses a major challenge in dipolar EPR spectroscopy data analysis.

3 Current approaches

Currently, two families of methods are commonly used in dipolar EPR spectroscopy data analysis. They differ in whether the distance distribution is represented as a parametric model or as a non-parametric model. Both methods stabilize the solution and are widely used since they are simple and often effective.

In the following, we discuss these two families of methods from the perspective of DeerLab.

3.1 Non-parametric distributions

If P is represented as a non-parametric vector, regularization methods are used to determine the solution. If background and modulation depth are known and fixed, the associated regularized optimization problem has the form

$$\begin{array}{}\text{(10)}& {\mathit{P}}_{\mathrm{fit}}=\underset{\mathit{P}\ge \mathrm{0}}{\mathrm{argmin}}\left({∥{\mathit{V}}_{\exp }-\mathbf{K}\mathit{P}∥}^{\mathrm{2}}+{\mathit{\alpha }}^{\mathrm{2}}\mathcal{R}\left(\mathbf{L}\mathit{P}\right)\right).\end{array}$$

The first term represents the sum of squared residuals, i.e., χ² without normalization by the noise variance, which we assume constant across the signal (see Eq. 8). It quantifies the quality of the fit of the model to the data. The second term is an additional penalty term that, together with the non-negativity constraint P≥0, stabilizes the solution. The regularization matrix L is a numerical approximation of a differential operator to impose smoothness, and α is the regularization parameter, which controls the balance between data agreement and regularization.

Regularization methods differ in the choice of the penalty norm ℛ. Tikhonov regularization (Tikhonov, 1963) was the first regularization approach introduced for dipolar data analysis (Bowman et al., 2004; Chiang et al., 2005a; Jeschke et al., 2004). Other approaches such as total variation (TV), Huber regularization, and Osher's Bregman-iterated regularization are beneficial in certain cases (Fábregas Ibáñez and Jeschke, 2019).

In the absence of the non-negativity constraint, the problem in Eq. (10) could be solved directly by ${\mathit{P}}_{\mathrm{fit}}=\stackrel{\mathrm{‾}}{\mathbf{K}}{\mathit{V}}_{\exp }$, with the α-dependent regularized pseudoinverse of the kernel matrix

$$\begin{array}{}\text{(11)}& \stackrel{\mathrm{‾}}{\mathbf{K}}={\left({\mathbf{K}}^{\mathrm{T}}\mathbf{K}+{\mathit{\alpha }}^{\mathrm{2}}{\mathbf{L}}^{\mathrm{T}}\mathbf{L}\right)}^{-\mathrm{1}}{\mathbf{K}}^{\mathrm{T}}\end{array}$$

in the case of Tikhonov regularization. However, due to the non-negativity constraint, this is not possible, and Eq. (10) is solved using non-negative least-squares optimization algorithms (Lawson and Hanson, 1974; Bro and Jong, 1997; Chen and Plemmons, 2009). $\stackrel{\mathrm{‾}}{\mathbf{K}}$ will be useful for uncertainty analysis (see Sect. 9).

The selection of α has been optimized for a large test set in previous work (Edwards and Stoll, 2018) and replicated (Fábregas Ibáñez and Jeschke, 2019), revealing that selection methods such as the Akaike information criterion (AIC) (Akaike, 1974) or general cross-validation (GCV) (Golub et al., 1979) can be superior to L-curve criteria (Hansen, 2000; Chiang et al., 2005a; Jeschke et al., 2006) if Gaussian white noise is the only source of error.

DeerLab provides a flexible regularization framework including all the aforementioned α-selection, non-iterated and iterated regularization methods, as well as a selection of solvers. In particular, implementation of Huber and TV regularization is much improved compared to our original study (Fábregas Ibáñez and Jeschke, 2019), such that these now perform similarly to Tikhonov regularization, invalidating the conclusions on the under-performance of Huber and TV regularization from that initial work.

Figure 1 shows how regularization approaches can be compared using DeerLab. Figure 1a presents an example of a low-noise dipolar signal processed via Tikhonov, TV and Osher's Bregman (Tikhonov) iterated regularization. All methods using the AIC for α-selection yield similar results. Figure 1b illustrates how regularization methods perform in the presence of strong noise. In such cases more differences arise between the different methods.

Figure 1 Analysis of background-free dipolar signals with regularization. A distance distribution is fit to (a) a low-noise signal and (b) a high-noise signal using Tikhonov regularization (blue), TV regularization (orange), and Osher's Bregman-iterated (OBIR) regularization (red). In all cases, the AIC was used for the selection of the regularization parameter (see Fig. S1). The input data (simulated) are shown as grey dots and the ground truth distance distribution is shown as grey shaded areas.

The outcomes of regularization analysis depend strongly on the choice of penalty norm, regularization operator, and α. For the remainder of this work, if not specified otherwise, we will use the Tikhonov penalty equipped with the second-order difference operator L₂ and the AIC for α-selection.

3.2 Parametric models

In an alternative representation, P is described as a parametric model P[θ] with $\left(\mathit{P}\right[\mathit{\theta }]{)}_i=P(t_i,\mathit{\theta })$, where θ is a vector of a small number of parameters. This is fit to the data using

$$\begin{array}{}\text{(12)}& {\mathit{\theta }}_{\mathrm{fit}}=\underset{\mathit{\theta }}{\mathrm{argmin}}{∥{\mathit{V}}_{\exp }-\mathbf{K}\mathit{P}\left[\mathit{\theta }\right]∥}^{\mathrm{2}}.\end{array}$$

Since solving the inverse problem by fitting a non-linear parametric model with only a few parameters is a well-conditioned problem, it can be solved without the need for regularization.

While bimodal Gaussian distributions were the first parametric models (Pannier et al., 2000a), the idea was generalized to a linear combination of N Gaussian distributions (Sen et al., 2007) (which we will refer to as a multi-Gauss model)

$$\begin{array}{}\text{(13)}& \mathit{P}\left[\mathit{\theta }\right]=\sum _{i=\mathrm{1}}^N{a}_i{\mathit{p}}_i[{\stackrel{\mathrm{‾}}{r}}_i,{\mathit{\sigma }}_i],\end{array}$$

where a_i are the amplitudes and p_i are the normalized Gaussian basis functions parameterized by their center distances ${\stackrel{\mathrm{‾}}{r}}_i$ and widths σ_i. Other parameterizations of the amplitudes can be used (Brandon et al., 2012; Stein et al., 2015).

To determine the optimal number N of Gaussians in the multi-Gauss model, Sen et al. (2007) proposed a statistical F-test, while Stein et al. (2015) and Hustedt et al. (2018) introduced the corrected Akaike information criterion (AICc) (Sugiura, 1978; Hurvich and Tsai, 1989) and the Bayesian information criterion (BIC) (Schwarz, 1978). Two more recent approaches utilize Monte Carlo simulations to determine the optimal multi-Gauss model (Dzuba, 2016; Timofeev et al., 2018). Parametric models are not limited to Gaussian basis functions. Many other basis function types (or mixtures thereof) can be employed, e.g., 3D Rice distributions (Domingo Köhler et al., 2011), spherical distributions (Ionita et al., 2008; Kattnig and Hinderberger, 2013), random-coil models (Fitzkee and Rose, 2004), or worm-like chain models (Wilhelm and Frey, 1996).

In a milestone for parametric modeling, Brandon et al. expanded the use of parametric models to include the modulation depth and a stretched-exponential background in the analysis of the signal in their software GLADD/DD (Brandon et al., 2012; Stein et al., 2015). This results in a time-domain parametric model

$$\begin{array}{}\text{(14)}& \mathit{V}\left[\mathit{\theta }\right]=\mathit{V}[\mathit{\lambda },{\mathit{\theta }}_P,{\mathit{\theta }}_B]=\mathbf{K}[\mathit{\lambda },{\mathit{\theta }}_{\mathrm{B}}]\mathit{P}\left[{\mathit{\theta }}_P\right].\end{array}$$

In this, the parameter vector θ includes not only the distance distribution parameters θ_P, but also the modulation depth λ and the background parameters θ_B. In general, a parametric time-domain model can be fit to the experimental data by solving

$$\begin{array}{}\text{(15)}& {\mathit{\theta }}_{\mathrm{fit}}=\underset{\mathit{\theta }}{\mathrm{argmin}}{∥{\mathit{V}}_{\exp }-\mathit{V}\left[\mathit{\theta }\right]∥}^{\mathrm{2}}.\end{array}$$

DeerLab expands upon this, allowing the design and fitting of any kind of time-domain or distance distribution parametric model. Automated multi-Gauss fitting and model selection using AIC, BIC, and other metrics are provided as well. In Fig. 2 we provide such an example of multi-Gauss fitting with DeerLab. The signal is fit using the time-domain model Eq. (14) and a varying number of Gaussians as the distance distribution model. Model selection based on the AIC determines the most parsimonious number of Gaussians, with a decent fit of the distance distribution.

Figure 2 Time-domain multi-Gauss fitting of a simulated four-pulse DEER signal. (a) The data are given as grey dots, the ground truth distance distribution as a shaded area, and the corresponding fits of a 3-Gauss model as blue lines. The background fit is given as dashed orange line. (b) The difference in AIC values as a function of the number of Gaussians is shown in blue, and the corresponding fits are given in the insets. The differences Δ_AIC are relative to the model with the lowest AIC value. The corresponding Akaike weights (see Sect. 10) are given next to each model. The model with the lowest AIC value (i.e., largest Akaike weight) is selected as the optimal model.

It is, however, crucially important to keep in mind that approaches based on parametric models may suffer from selection bias, i.e., the bias introduced by limiting the analysis to a specific family of models and by using a particular criterion for model selection within that family (Freedman, 1983; Lukacs et al., 2009). Parametric model fits may also be affected by confirmation bias, i.e., the tendency to process data in a way that matches one's preconceptions and avoids contradiction of prior belief (Nickerson, 1998).

4 One-step analysis

Neither parametric model fitting nor regularization is ideal. Parametric models can be fit directly to the full time-domain signal but are strongly limited by how well they can approximate the ground truth distribution. On the other hand, regularization yields non-parametric distance distributions that accommodate a much larger range of ground truths, but regularization cannot directly include other parameters such as the background and the modulation depth.

With regularization methods, it is therefore common practice to use a two-step approach: (1) fit a parametric background model to the time-domain signal and correct the signal by the fitted background (either by division or subtraction) and (2) apply regularization to the background-corrected signal. Software based on this approach includes DeerAnalysis (Jeschke et al., 2006) and LongDistances (Altenbach, 2020). Including the fitted background in the kernel for step (2), as in Eq. (3), also does not eliminate the need for the two-step approach (Fábregas Ibáñez and Jeschke, 2020).

The two-step analysis as done in DeerAnalysis is sub-optimal, because the background fit in step (1) relies on the assumption that the oscillation periods in the time-domain signal are much shorter than the overall signal length. Many experimentally observed signals do not satisfy this assumption. Therefore, the two-step analysis cannot robustly process these types of signals.

The most desirable approach is to simultaneously fit both the time-domain parameters θ and a non-parametric distance distribution P to the time-domain signal V_exp in one step, i.e.,

$$\begin{array}{}\text{(16)}& \left({\mathit{\theta }}_{\mathrm{fit}},{\mathit{P}}_{\mathrm{fit}}\right)=\underset{\mathit{\theta },\mathit{P}\ge \mathrm{0}}{\mathrm{argmin}}F(\mathit{\theta },\mathit{P}),\end{array}$$

with the regularized objective function

$$\begin{array}{}\text{(17)}& F(\mathit{\theta },\mathit{P})={∥{\mathit{V}}_{\exp }-\mathbf{K}\left[\mathit{\theta }\right]\mathit{P}∥}^{\mathrm{2}}+{\mathit{\alpha }}^{\mathrm{2}}\mathcal{R}\left(\mathbf{L}\mathit{P}\right).\end{array}$$

Equation (16) can be solved using a variety of constrainable non-linear optimization algorithms by combining θ and P into a single parameter vector and applying all the constraints for θ and P (Altenbach, 2020). However, this does not take full advantage of the special structure of the problem, i.e., that it is a penalized least-squares problem with a model for V that is linear in P and non-linear in θ.

To take advantage of this structure, DeerLab implements a nested subspace optimization algorithm based on separable non-linear least-squares (Budil et al., 1996; Golub and Pereyra, 2003; Sima and Van Huffel, 2007). It separates the θ and P spaces and uses a non-linear least-squares algorithm to solve

$$\begin{array}{}\text{(18)}& {\mathit{\theta }}_{\mathrm{fit}}=\underset{\mathit{\theta }}{\mathrm{argmin}}F\left(\mathit{\theta },\mathit{P}\left[\mathit{\theta }\right]\right),\end{array}$$

where P[θ] is the optimal non-parametric distance distribution for a given θ, determined via regularization

$$\begin{array}{}\text{(19)}& \mathit{P}\left[\mathit{\theta }\right]=\underset{{\mathit{P}}^{\prime }\ge \mathrm{0}}{\mathrm{argmin}}F\left(\mathit{\theta },{\mathit{P}}^{\prime }\right)\end{array}$$

as in Eq. (10), using a dedicated non-negative linear least-squares algorithm. (Note that P is now a parametric model since it depends on θ, although this includes only modulation depth and background parameters.) Once θ_fit is obtained, P_fit is calculated as P[θ_fit]. The algorithm is iterative and is illustrated in Fig. 3. It starts with an initial guess for θ and determines the associated P from Eq. (19). This P is then used by the algorithm of Eq. (18) to determine the next θ, which is then used again by the algorithm in Eq. (19) to get the next P[θ], and so on until convergence is reached. The optimization in Eq. (19) can utilize any form of regularization, and it can be run with a fixed regularization parameter α or optimize it each time.

Figure 3 One-step analysis of dipolar signals using separable non-linear least-squares optimization. A set of parameters θ_n (blue) is computed by optimization of Eq. (18) for some given experimental signal V_exp (grey dots). For each θ_n, a corresponding distance distribution P_n (orange) is computed by optimization of Eq. (19). This procedure is repeated until a minimum of the objective function is found for an optimal parameter set θ_fit, leading to the corresponding optimal distribution P_fit.

Figure 4 Comparison of one-step and two-step analysis. A stretched-exponential background and a non-parametric distance distribution were fit to a simulated four-pulse DEER signal. The analysis was done by either fitting both simultaneously (blue) or by fitting the background followed by the distribution (two-step analysis, orange). The analysis was repeated on the same signal with increasing truncation (a–d). The two-step analysis fails in case (d). The data are given as grey dots, and the fitted signal and background are given as solid and dashed colored lines, respectively. The true background is given as a grey dashed line for reference. The non-parametric distance distributions obtained via Tikhonov regularization are given as, respectively, colored lines and the ground truth as a shaded area.

Figure 4 shows an example that compares this one-step approach to the traditional two-step analysis, using progressively more truncated dipolar signals. Analogous comparisons between two-step analysis and fully parametric models have been previously reported (Brandon et al., 2012; Stein et al., 2015). For a sufficiently long signal (Fig. 4a) both approaches yield similar results, as the signal is long enough for the oscillations to decay, facilitating the separate fit of the background in the two-step analysis. If the signal is truncated as in Fig. 4b–d, the two-step analysis cannot properly fit the background anymore. In contrast, one-step analysis correctly identifies the background in all cases and recovers the underlying non-parametric distance distribution with reasonable fidelity. It is important to keep in mind that highly truncated signals such as the ones shown in Fig. 4d can fail to provide correct results if the measurement noise is too large (see Fig. S2 in the Supplement).

This shows that DeerLab's ability to simultaneously fit the background and a non-parametric distribution opens up the possibility of fitting non-parametric distance distributions to signals which might have been deemed unanalyzable in the past.

5 Multi-pathway models

Dipolar EPR spectroscopy found widespread use with four-pulse DEER. Since then, experimental dipolar EPR spectroscopy has developed into a set of diverse techniques, with signals exhibiting a variety of features. However, the data processing has not evolved much from its four-pulse DEER origins.

The key operator in the analysis is the dipolar kernel K(t,r) (see Eq. 2). It represents the experiment; i.e., it describes how the time-domain signal is obtained from a given distance distribution. If the kernel cannot account for a feature in the signal, it is because the kernel model is incomplete. In these cases, it is preferable to improve the kernel rather than to tweak the signal into an artificial four-pulse DEER signal.

One example of this is RIDME, where it is known that the measured signals in $S>\mathrm{1}/\mathrm{2}$ systems contain overtones not present in four-pulse DEER (Razzaghi et al., 2014). If disregarded, these overtones cause distortions in the distance distribution if the four-pulse DEER kernel is used to analyze the data. This can be avoided by including the overtones in the model (Keller et al., 2017). In DeerLab, one can include a set of overtones with a background in the kernel model, which can be used to directly fit primary RIDME data via, e.g., regularization.

Other examples are multi-pulse DEER sequences, which generally are not fully modeled. All multi-pulse DEER experiments feature modulations in addition to the basic modulation of Eq. (3). Despite their clear dipolar origin, these are regarded as undesirable “artifacts”: the “2+1 artifact” in four-pulse DEER (Jeschke, 2012; Teucher and Bordignon, 2018) and “artifacts” or “residues” in five-pulse and seven-pulse DEER (Borbat et al., 2013; Spindler et al., 2015; Breitgoff et al., 2017b, a; Doll and Jeschke, 2017; Milikisiyants et al., 2018). Several experimental and processing approaches have been published that aim to remove these contributions from the total signal to recover the idealized dipolar evolution function (Borbat et al., 2013; Spindler et al., 2015; Teucher and Bordignon, 2018; Breitgoff et al., 2017b, a). These approaches introduce further experimental or theoretical complexity.

However, these additional contributions are actual dipolar signals. They may even provide strong oscillations at times when the oscillations from the main signal have decayed, thus increasing the signal-to-noise ratio of the experiment. Instead of removing these contributions due to the lack of a proper model, it is advantageous to extend the model to explicitly include these contributions.

DeerLab includes such an extended model for multi-pulse DEER, derivable from spin density matrix dynamics. In this model, which we call the dipolar multi-pathway model, the overall signal is a combination of several dipolar signals arising from dipolar pathways of varying amplitudes and refocusing times (Salikhov et al., 1981; Kutsovsky et al., 1990; Kattnig et al., 2013; Borbat et al., 2013; Spindler et al., 2015). The total signal is given by Eq. (2) with the general kernel

$$\begin{array}{}\text{(20)}& \begin{array}{rl}K(t,r)& =\left[{\mathrm{\Lambda }}_{\mathrm{0}}+\sum _{p=\mathrm{1}}^N{\mathit{\lambda }}_p{K}_{\mathrm{0}}\left(n_p\right(t-T_p),r)\right]\\ & \cdot \prod _{p=\mathrm{1}}^N{V}_{\mathrm{inter}}\left(n_p\left(t-T_p\right),{\mathit{\lambda }}_p\right),\end{array}\end{array}$$

where Λ₀ is the total contribution of the unmodulated dipolar pathways, the index p runs over all N modulated dipolar pathways, T_p are the refocusing times of the individual modulated dipolar pathways, λ_p are the amplitudes of the modulated dipolar pathways, n_p are the harmonics of the individual modulated pathways, and the background function is as in Eq. (6). A schematic representation is shown in Fig. 5. The kernel for standard four-pulse DEER from Eq. (3) is a special case of Eq. (20), with N=1, T₁=0, n₁=1 ${\mathrm{\Lambda }}_{\mathrm{0}}=\mathrm{1}-\mathit{\lambda }$, and λ₁=λ.

Figure 5 Schematic representation of a dipolar multi-pathway DEER signal. (a) Background decays (orange dashed lines) and dipolar evolution functions (blue solid lines) are shown for two different dipolar pathways. The unmodulated component Λ₀ is given as a solid grey line. (b) The overall signal (black) is given by the sum over all pathway dipolar evolution functions (blue) times the product over all pathway background decays (orange).

Discretization of the kernel in Eq. (20) gives the same expression as in Eq. (7),

$$\begin{array}{}\text{(21)}& \mathit{V}=\mathbf{K}\left[\mathit{\theta }\right]\mathit{P},\end{array}$$

with the parameter set θ comprising Λ₀, all λ_p, n_p, T_p, as well as κ_d and d.

Figure 6 shows examples of how DeerLab can analyze multi-pulse DEER data in terms of a multi-pathway model, thus alleviating the need for experimental correction protocols or signal pre-processing (beyond phase correction). Note that experimental schemes which generate multiple datasets with some pathways shifted with respect to each other (Breitgoff et al., 2017a) can profit from global analysis (vide infra) to stabilize the accurate estimation of pathway parameters.

Figure 6 One-step analysis of multi-pulse DEER signals with the dipolar multi-pathway model for (a) four-pulse DEER with the “2+1 artifact” and (b) five-pulse DEER with its four-pulse “artifact”. All data were fit with two modulated dipolar pathways. The data are given as grey dots, the signal and distribution fits are given as solid blue lines, the background fit is given as an orange dashed line, and the ground truth is given as a shaded area.

This analysis of multi-pulse DEER signals shows the benefits of using full models for dipolar signals instead of removing or avoiding “artifacts” to try to match partial models. DeerLab provides a compact framework for testing, developing, and applying more complete models.

6 Global analysis

Global analysis denotes the situation where a single model is fit simultaneously to multiple datasets. In dipolar EPR, this was first used for fitting a multi-Gauss distribution model to several DEER datasets (Brandon et al., 2012). Later, Tikhonov regularization was employed to fit short and long DEER traces simultaneously (Rein et al., 2018).

Global analysis is implemented in DeerLab for both parametric and non-parametric distance distribution models and for an arbitrary number of dipolar signals. The global optimization problem of a set of M dipolar signals V_exp,i using a model depending on parameters θ and on N non-parametric distance distributions P_j is

$$\begin{array}{}\text{(22)}& \left({\mathit{\theta }}_{\mathrm{fit}},\mathit{\left\{}{\mathit{P}}_{\mathrm{fit}}\mathit{\right\}}\right)=\underset{\mathit{\theta },\mathit{\left\{}\mathit{P}\mathit{\right\}}\ge \mathrm{0}}{\mathrm{argmin}}\left[F(\mathit{\theta },\mathit{\{}\mathit{P}\mathit{\left\}}\right)+G\left(\mathit{\right\{}\mathit{P}\mathit{\left\}}\right)\right],\end{array}$$

with $\mathit{\left\{}\mathit{P}\mathit{\right\}}=\mathit{\{}{\mathit{P}}_{\mathrm{1}},\mathrm{\dots },{\mathit{P}}_N\mathit{\}}$ and

$$\begin{array}{}\text{(23)}& {\displaystyle }F(\mathit{\theta },\mathit{\{}\mathit{P}\mathit{\left\}}\right)=\sum _{i=\mathrm{1}}^M{w}_i{\displaystyle \frac{\Vert {\mathit{V}}_{\exp ,i}-{\mathit{V}}_i\left[\mathit{\theta },\mathit{\left\{}\mathit{P}\mathit{\right\}}\right]{\Vert }^{\mathrm{2}}}{{\mathit{\sigma }}_i^{\mathrm{2}}}},\text{(24)}& {\displaystyle }G\left(\mathit{\right\{}\mathit{P}\mathit{\left\}}\right)={\mathit{\alpha }}^{\mathrm{2}}\sum _{j=\mathrm{1}}^N\mathcal{R}\left(\mathbf{L}{\mathit{P}}_j\right),\end{array}$$

with a similar expression without {P} and G({P}) for a fully parametric model. σ_i are the noise levels of the individual signals. The parameter vector θ includes the parameters needed to generate all signals V_i, where each V_i typically depends only on a subset of the parameters.

The quantities w_i in Eq. (23) are the weights that determine the contribution of each signal to the objective function. The default is w_i=1, meaning that each data point from each signal contributes equally, given its noise level, to the objective function. Different values of w_i can be used to indicate preferential weighing.

If all signals V_i derive from a single distance distribution, then we can use Eq. (22) together with

$${\mathit{V}}_i[\mathit{\theta },\mathit{P}]={\mathbf{K}}_i\left[\mathit{\theta }\right]\mathit{P}$$

or

$$\begin{array}{}\text{(25)}& {\mathit{V}}_i\left[\mathit{\theta }\right]={\mathbf{K}}_i\left[\mathit{\theta }\right]\mathit{P}\left[\mathit{\theta }\right].\end{array}$$

Here, each K_i describes a different experiment (different pulse sequence, different trace length, etc.) and depends on a subset of the parameters in θ. The distance distribution P can be either parametric (in which case θ includes the distribution parameters) or it can be non-parametric. In the latter case, Eq. (22) is solved using separable non-linear least squares. As an example, in Fig. 7 we simultaneously fit a four-pulse DEER signal and a five-pulse DEER signal with its secondary four-pulse pathway contribution, using a model with a single non-parametric distribution but separate backgrounds and modulation depths for the two signals. The distance distribution underlying both signals is nicely recovered.

Figure 7 Global analysis of a four-pulse DEER signal and a five-pulse DEER signal with its secondary four-pulse pathway contribution, both derived from the same distance distribution. The simulated data are given as grey dots, the ground truth distribution is given as a shaded area, and the signal and distribution fits are given as solid blue lines.

Note, however, that the analysis of dipolar signals of different length or obtained under different dynamical decoupling conditions may be inconsistent if different conformers have different dephasing rates (Baber et al., 2015).

Another common global analysis situation is when the measured signals stem from several samples containing a mixture of chemical or structural components, each with its own distinct distribution, and related to each other via additional conditions such as a chemical equilibrium. The simplest case is a system equilibrated between two forms A and B (A⇌B), e.g., a protein–ligand binding equilibrium or a monomer–dimer equilibrium. In this case, the dipolar signals are described as

$$\begin{array}{}\text{(26)}& \begin{array}{rl}{\mathit{V}}_i[\mathit{\theta },\mathit{\{}\mathit{P}\mathit{\left\}}\right]& ={\mathit{V}}_i\left[\mathit{\right\{}{\mathit{\theta }}_K,x_A\mathit{\}},\mathit{\{}{\mathit{P}}_A,{\mathit{P}}_B\mathit{\left\}}\right]\\ & ={\mathbf{K}}_i\left[{\mathit{\theta }}_K\right]\left[x_{A,i}{\mathit{P}}_A+\left(\mathrm{1}-x_{A,i}\right){\mathit{P}}_B\right],\end{array}\end{array}$$

with both component distributions P_A and P_B being fitted along the parameter set $\mathit{\theta }=\mathit{\{}{\mathit{\theta }}_K,\mathit{\{}{x}_{A,i}\mathit{\left\}}\mathit{\right\}}$ via Eq. (22). The mole fractions x_A,i depend on the location of the equilibrium, which might vary among the samples via ligand concentration, matrix composition, and other factors. Either the mole fractions or the underlying equilibrium constant can be included among the fitting parameters.

Such titration or dose–response datasets have been analyzed using multi-Gaussian distribution models for the component distributions (Stein et al., 2015; Martens et al., 2016; Collauto et al., 2017; Barth et al., 2018; Jagessar et al., 2020). As discussed above, however, non-parametric distributions may often be a preferable choice. DeerLab enables global fitting of an arbitrary number of datasets with regularization approaches and, thus, analysis of titration datasets in terms of non-parametric distributions. As an example, Fig. 8 shows such an analysis of a protein–ligand binding assay, using signals with different noise levels, trace lengths, backgrounds, and modulation depths (Fig. 8a) with a model that includes the bound-protein mole fractions among the parameters. The global analysis gives good fits to the time-domain data and results in two non-parametric distributions that capture the underlying ground truth well (Fig. 8b). The extracted mole fractions are consistent with the underlying binding equilibrium (Fig. 8c). Alternatively, one can skip the separate determination of the mole fractions and include the dissociation constant directly as a parameter in the global analysis of the dipolar data. Even in cases where one or both extremes of the binding curve are experimentally unavailable, such an analysis is still feasible albeit at the cost of larger uncertainty in the fitted molar fraction or dissociation constants. In summary, this example shows that DeerLab allows global analysis of titration datasets using non-parametric distributions.

Figure 8 Global fitting of titration data of a protein–ligand binding equilibrium with non-parametric distributions. In panel (a) four-pulse DEER traces with different trace lengths, modulation depths, backgrounds, and noise levels were simulated for different ligand concentrations added to a protein concentration of 1 µM. These parameters as well as the mole fractions of bound/unbound protein for each trace and two non-parametric distance distributions for the bound and unbound states (via Tikhonov regularization) were fit simultaneously. The simulated data are given as grey dots and the fitted signals and backgrounds are given as solid and dashed blue lines, respectively. In panel (b) the distance distribution fits for the unbound (red) and bound (orange) states are given as well as the combined fitted distribution (blue) for the different ligand concentrations. The ground truth sum distributions are given as grey shaded areas. In panel (c) the fitted mole fractions of the unbound (red) and bound (orange) states are given as colored dots. The solid grey lines represent the ground truth for a dissociation constant of K_d=5.65 µM.

7 Global and local minima

Even when a fitted model agrees with the experimental data and a minimum of the objective function has been located, it does not mean that the only or the best fit has been found. The objective function can have multiple minima, meaning that the located minimum is not necessarily the global minimum (see Fig. 9). Which minimum is found depends mainly on where the search is started. The boundaries set on the parameter space can also influence the outcome if a minimum is located outside these boundaries. Other factors, e.g., the numerical algorithms and convergence parameters used for the optimization, can affect this as well.

Figure 9 Global vs. local minima. During optimization of a parameter θ, by minimization of an objective function f(θ) (black line), several local minima (blue) might be found instead of the global minimum (orange). The global minimum can be found by varying the starting point, the lower/upper boundaries θ_upper and θ_lower (to find minima which might be outside the bounds, e.g., θ_local,3), or the convergence parameters ε of the numerical solver such that some local minima might be ignored.

While there are dedicated global optimization algorithms, the simplest approach to find a global minimum is to repeat the optimization process with different starting values in order to explore the parameter space more fully (Sen et al., 2007). After sampling enough starting points, a set of minima is found and the one with the lowest objective function value is taken as the global minimum (see Fig. 9).

While these procedures can be costly, it is recommended to routinely check that variation of algorithmic parameters does not yield a lower minimum.

8 Goodness of fit

After obtaining a fitted model, an important step is to assess the goodness of the fit. If the fit is not good, either the optimization failed to locate an appropriate minimum (see previous section) or the model is inapplicable (e.g., oversimplified) for the given experimental data, and a better model needs to be used.

Figure 10 Uncertainty analysis. In panel (a) a simulated noisy four-pulse DEER signal with exponential background (grey dots) is fitted with a non-parametric distance distribution (red line). In panel (b), the red bars show the histogram of normalized residual values and the reference standard Gaussian distribution as a grey shaded area. In panels (c) and (d) the fitted distance distribution is given as blue/orange lines as well as the 50 % and 95 % confidence intervals (shaded areas) obtained via the standard method (orange, top) and via bootstrapping (blue, bottom). The ground truth distribution is shown as a shaded grey area for reference. The estimated distribution of the fitted modulation depth λ (e) and background decay rate constant κ (f) obtained via the standard method are given as orange lines. The kernel density estimations obtained via bootstrapping are given as blue lines. The true values of λ and κ are given as a dashed grey line for reference. All bootstrap results were obtained from 1000 bootstrap samples and the regularization parameter was optimized using AIC for the original signal and fixed throughout the analysis.

To assess the goodness of fit, several procedures can be utilized. A direct test is to compare the histogram of the normalized residuals $|V_{\exp ,i}-V_i[\mathit{\theta }\left]\right|/\mathit{\sigma }$ to the standard normal distribution 𝒩(0,1). Here, σ is the noise standard deviation estimated either from a dataset with several, individually stored scans (Edwards and Stoll, 2016), from the standard deviation of the (flat) imaginary part of the signal, or from the residuals of the signal minus the fitted model. Accuracy of the estimate of σ decreases in the sequence indicated, due to a possible imbalance in quadrature channels and the slight model inadequacy that is unavoidable for ill-posed problems. An example of this is shown in Fig. 10b. If the histogram deviates strongly from a Gaussian distribution, the fit is considered inadequate. Alternatively, the comparison can also be based on a statistical test.

Another method is to examine the reduced χ²-value:

$$\begin{array}{}\text{(27)}& {\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{dof}}}}{\displaystyle \frac{{∥V_{\exp }-V\left[{\mathit{\theta }}_{\mathrm{fit}}\right]∥}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}.\end{array}$$

Here, N_dof is the number of degrees of freedom, which can be taken as approximately equal to the number of data points minus the number of model parameters. A good fit is characterized by ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}\approx \mathrm{1}$. Note that the use of ${\mathit{\chi }}_{\mathrm{red}}^{\mathrm{2}}$ for non-linear models is not rigorous. Also, the notion of (effective) number of parameters for a non-parametric distribution model, estimated as $\mathrm{tr}\left(\mathbf{K}\stackrel{\mathrm{‾}}{\mathbf{K}}\right)$ with $\stackrel{\mathrm{‾}}{\mathbf{K}}$ defined in Eq. (11), is not straightforward (Edwards and Stoll, 2016).

Additional methods for assessing goodness of fit are discussed in Budil et al. (1996) and Hansen et al. (2012).

9 Uncertainty analysis

Up to this point, we took advantage of knowing the ground truth when we assessed quality of the solutions. However, in experimental work the ground truth is unknown. The presence of noise in experimental signals introduces uncertainty about the underlying noise-free dipolar signal and results in uncertainty about the model parameters. This, in turn, affects the strength of the conclusions that can safely be drawn. For example, $r_{\mathrm{0}}=(\mathrm{3.2}\pm \mathrm{0.1})$ nm supports much more confident conclusions about r₀ than $r_{\mathrm{0}}=(\mathrm{3.2}\pm \mathrm{0.9})$ nm. It is therefore crucial to always report uncertainty estimates for all fitted parameters and distance distributions extracted from experimental data, for model parameters θ_i as well as vector elements P_j of non-parametric distance distributions. Reporting fitted values without accompanying uncertainty estimates is meaningless. Several approaches have been proposed for uncertainty estimation, including validation of the regularization model (Jeschke et al., 2006; Altenbach, 2020), iterative scanning of the χ²-surface (Brandon et al., 2012; Stein et al., 2015), covariance matrices (Stein et al., 2015; Hustedt et al., 2018), and Bayesian inference (Edwards and Stoll, 2016; Sweger et al., 2020).

DeerLab provides two separate methods for uncertainty estimation for both model parameters θ_i and for vector elements P_j of non-parametric distance distributions.

The first method estimates parameter uncertainties from the covariance matrix (Budil et al., 1996). For a fully parametric model with parameter vector θ, this is well established (Hustedt et al., 2018). The covariance matrix for θ is

$$\begin{array}{}\text{(28)}& {\mathbf{\Sigma }}_{\mathit{\theta }}={\mathit{\sigma }}^{\mathrm{2}}{\left({\mathbf{J}}^{\mathrm{T}}\mathbf{J}\right)}^{-\mathrm{1}},\end{array}$$

where J is the Jacobian matrix of derivatives with elements $J_{ij}=\partial V_j\left[\mathit{\theta }\right]/\partial {\mathit{\theta }}_i$ evaluated at θ=θ_fit and is calculated using numerical differentiation.

There is also a simple way to obtain the covariance matrix Σ_P for the elements of a non-parametric P determined via regularization, but only if the non-negativity constraint for P is disregarded. It is obtained by propagating the time-domain noise covariance matrix σ²I (see Eq. 8) to the distance domain by (Weese, 1992; Kasper et al., 2002)

$$\begin{array}{}\text{(29)}& {\mathbf{\Sigma }}_{\mathit{P}}={\mathit{\sigma }}^{\mathrm{2}}\stackrel{\mathrm{‾}}{\mathbf{K}}{\stackrel{\mathrm{‾}}{\mathbf{K}}}^{\mathrm{T}},\end{array}$$

with $\stackrel{\mathrm{‾}}{\mathbf{K}}$ defined in Eq. (11). Note that this can be utilized even for models that depend on other parameters, since V always depends linearly on P, as shown above.

From the covariance matrices in Eqs. (28) and (29), the standard errors of a parameter θ_i or a distribution vector element P_i are obtained as

$${\mathit{\sigma }}_{{\mathit{\theta }}_i}=\sqrt{({\mathbf{\Sigma }}_{\mathit{\theta }}{)}_{ii}}$$

and

$$\begin{array}{}\text{(30)}& {\mathit{\sigma }}_{P_i}=\sqrt{({\mathbf{\Sigma }}_{\mathit{P}}{)}_{ii}}.\end{array}$$

These are used to estimate symmetric confidence intervals (CIs) around the fitted parameter for a confidence level γ with boundaries

$${\mathit{\theta }}_{\mathrm{fit},i}\pm z_{\mathit{\gamma }}{\mathit{\sigma }}_{{\mathit{\theta }}_i}$$

and

$$\begin{array}{}\text{(31)}& P_{\mathrm{fit},i}\pm z_{\mathit{\gamma }}{\mathit{\sigma }}_{P_i},\end{array}$$

where z_γ is the γ-quantile of a standard normal or Student's t-distribution.

This method for estimating CIs is simple and stands on a sound theoretical basis. However, it has several limitations. (1) It approximates the parameter likelihood by a Gaussian distribution centered around the fitted θ_fit and P_fit. While this is often a reasonable approximation, it can still fail in many cases. (2) It assumes that the parameters are unbounded. This is not fulfilled in the analysis of dipolar signals, as most parameters are constrained to a certain range, e.g., $\mathrm{0}\le \mathit{\lambda }\le \mathrm{1}$ and P_i≥0. Although the boundary conditions can be imposed by cropping the CIs in Eq. (31) to the parameter range (Wang, 2008), the calculation of the CIs is still based on an unbounded assumption. Hence, the CIs do not include the additional information provided by the constraints. (3) In the case of a model that depends on both parameters θ and a non-parametric distribution P, the covariances between θ_i and P_j are not accounted for, potentially leading to underestimation of uncertainties.

A more general and accurate method for the estimation of parameter uncertainty involves bootstrapping (Efron and Tibshirani, 1986; Banks et al., 2010). This is a Monte Carlo resampling method based on generated synthetic signals with different noise realizations. The variant implemented in DeerLab first generates N synthetic traces V_k by adding N different noise realizations to the fitted model V[θ_fit,P_fit]. The noise is drawn from a Gaussian distribution with standard deviation σ estimated from the experimental data. Then, the N bootstrap traces V_k are analyzed in the same fashion as the original dataset V_exp, resulting in N+1 fitted parameter vectors θ_fit,k and distance distributions P_fit,k. The distributions of the parameter values and the distance distribution vector elements are then taken as approximations of the underlying parameter uncertainty distributions.

The bootstrap method, while costly, has several important advantages over the method based on covariance matrices. All information provided by parameter constraints is included in the estimation. Additionally, model-free estimations of parameter distributions are obtained, without the need to assume a Gaussian distribution. These can be statistically analyzed in multiple ways to quantify the parameter uncertainty. For example, in analogy to above, one can define the confidence interval for a parameter θ with confidence level $\mathit{\gamma }=\mathrm{1}-\mathit{\alpha }$ as

$$\begin{array}{}\text{(32)}& \left({\mathit{\theta }}_{\mathrm{1}-\mathit{\alpha }/\mathrm{2}},{\mathit{\theta }}_{\mathit{\alpha }/\mathrm{2}}\right),\end{array}$$

where ${\mathit{\theta }}_{\mathrm{1}-\mathit{\alpha }/\mathrm{2}}$ and θ_α∕2 are the $(\mathrm{1}-\mathit{\alpha }/\mathrm{2})$th and (α∕2)th percentiles of the bootstrapped θ distribution.

Figure 10 shows an example of parameter uncertainty estimation using both the standard method and bootstrap. Although both methods lead to similar confidence intervals, the bootstrapped solution provides narrower intervals thanks to the use of the additional information provided by parameter boundaries and the non-negativity constraint P≥0. While the standard method provides an easily accessible uncertainty estimation, the use of bootstrapping is recommended for producing final results.

While the distance distribution and its uncertainty analysis provide the complete information of the data, in application work it is often of interest to determine summarizing quantities such as the distance mode, the mean distance, or the standard deviation of the distance distribution. It is important to report uncertainties for these quantities as well. If the uncertainty analysis is based on covariance matrices, these uncertainties can be calculated from the covariance matrices Σ_P and Σ_θ using error propagation via the Jacobian. For instance, the variance of the mean distance $\stackrel{\mathrm{‾}}{r}$ is

$$\begin{array}{}\text{(33)}& {\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{r}}^{\mathrm{2}}=\mathbf{J}\left(\begin{array}{cc}{\mathbf{\Sigma }}_{\mathit{P}}& \mathrm{0}\\ \mathrm{0}& {\mathbf{\Sigma }}_{\mathit{\theta }}\end{array}\right){\mathbf{J}}^{\mathrm{T}},\end{array}$$

where J is the Jacobian containing the derivatives of $\stackrel{\mathrm{‾}}{r}$ with respect to all model parameters (θ and P). Note that this method does not work for the distance mode. In general, it is preferable to use the more accurate bootstrap method to determine the distance mode, mean, median, and similar quantities. The means and uncertainties of these quantities are easily calculated from their histograms obtained from the ensembles of fitted θ and P.

10 Model comparison

The above uncertainty analysis is performed under the assumption of a specific model for the distance distribution (parametric or non-parametric), the background, and the experiment. All estimated uncertainties capture variability only within this assumed model. The possibility that the data could be explained equally well, or better, by other models is not incorporated.

Model selection approaches, as outlined above for parametric and non-parametric distance distribution (selection of number of Gaussians, regularization parameter selection), provide convenient quantitative decision criteria for picking one distribution model over another in a principled fashion. On the other hand, background models are often chosen ad hoc. Finally, the choice of the standard DEER experiment model in Eq. (3) is based on physical assumptions (no orientation selection, no bandwidth limitations, no conformer-dependent phase memory times, no exchange couplings, etc.) that might not all be fully valid for a given experimental situation. Whether principled or ad hoc, any model selection eliminates model uncertainty from the analysis and leads to bias.

It is therefore preferable to compare and report the relative performance of a series of plausible models, without picking a winner. This can be accomplished using Akaike weights, defined as (Burnham and Anderson, 2003)

$$\begin{array}{}\text{(34)}& w_{\mathrm{AIC},i}={\displaystyle \frac{{\mathrm{e}}^{-{\mathrm{\Delta }}_{\mathrm{AIC},i}/\mathrm{2}}}{{\sum }_{k=\mathrm{1}}^M{\mathrm{e}}^{-{\mathrm{\Delta }}_{\mathrm{AIC},k}/\mathrm{2}}}},\end{array}$$

where Δ_AIC,i is the difference between the AIC value of model i and the lowest AIC value within the set of models. The Akaike weights give the probability that model i is the best among the set of M candidate models, given the data. This can be used to compare a set of different parametric models or to compare a series of non-parametric distribution models differing in the regularization parameter α.

Figure 11a illustrates this for a set of parametric models differing in complexity (number of components) and type of basis function. The analysis finds Akaike weights of about 25 % for the 2-Gauss, 3-Gauss and 1-skewed Gauss models, indicating that the data do not provide enough evidence to clearly identify a best model. Figure 11b shows a comparison for non-parametric distribution models obtained with different α-values, where the AIC is calculated using $\mathrm{tr}\left(\mathbf{K}\stackrel{\mathrm{‾}}{\mathbf{K}}\right)$ (Edwards and Stoll, 2018). It is apparent that models over a range of α-values are similarly likely. Therefore, in neither case is there a clear “winner” model. Note that uncertainty in the regularization parameter α could also be propagated to the resulting distance distribution and thereby included in the confidence bands. In cases where this is not applicable, and depending on the conclusions that one wants to draw, it may be prudent to list all models that fit the data reasonably well as opposed to picking the model with the highest weight.

Figure 11 Model comparison using Akaike weights. The Akaike weights (in percentages) are given for (a) a set of parametric distance distribution models with varying number n of ordinary (G-n), skewed (sG-n) and generalized (gG-n) Gaussians and (b) for non-parametric distance distributions determined over a set of values for the regularization parameter α. For each model in panels (a) and (b), the corresponding distance distribution is shown above the bar. Panels (a) and (b) are based on different simulated DEER signals of different distributions.

Despite these simple numerical procedures for model comparison, researchers need to use careful judgment in which models are included in a comparison and explicitly disclose the reasoning behind all model choices. Also, if the correct model is not included, then even such model comparison will lead to biased or incorrect conclusions.

11 Concluding remarks

Dipolar EPR spectroscopy requires reliable and robust data processing tools. The associated software should also be flexible and adaptable to quickly incorporate new developments in the field. DeerLab collects many existing as well as several novel data analysis methods in a flexible, robust, and reliable manner.

By allowing analysis workflows to be shared, DeerLab provides a significant step towards solving reproducibility issues plaguing dipolar EPR spectroscopy analysis. With DeerLab, the analysis scripts can be provided along with published data, leading to improved reproducibility. Through an online repository (see “Code availability” below), all DeerLab versions remain available, further enhancing reproducibility.

DeerLab can serve as a powerful method development tool. It allows the discovery and testing of new processing techniques. We illustrated this by introducing several extensions: one-step analysis with non-parametric distributions, the dipolar multi-pathway model, global analysis with non-parametric distribution models, and uncertainty estimation using bootstrap. The discussion of new methods or comparisons between established ones based on statistical arguments is also largely facilitated by such a scriptable tool.

Although DeerLab does not come with a GUI, it can serve as the data processing engine for intuitive and dedicated GUI-based data analysis tools. These are essential for the robust and successful application of routine dipolar EPR spectroscopy in fields such as structural biology or materials science.

In conclusion, DeerLab provides a unified and extensible open-source platform for data analysis in dipolar EPR spectroscopy, opening up a new world of data processing workflows.