Statistical analysis of ENDOR spectra.

Electron-nuclear double resonance (ENDOR) measures the hyperfine interaction of magnetic nuclei with paramagnetic centers and is hence a powerful tool for spectroscopic investigations extending from biophysics to material science. Progress in microwave technology and the recent availability of commercial electron paramagnetic resonance (EPR) spectrometers up to an electron Larmor frequency of 263 GHz now open the opportunity for a more quantitative spectral analysis. Using representative spectra of a prototype amino acid radical in a biologically relevant enzyme, the [Formula: see text] in Escherichia coli ribonucleotide reductase, we developed a statistical model for ENDOR data and conducted statistical inference on the spectra including uncertainty estimation and hypothesis testing. Our approach in conjunction with 1H/2H isotopic labeling of [Formula: see text] in the protein unambiguously established new unexpected spectral contributions. Density functional theory (DFT) calculations and ENDOR spectral simulations indicated that these features result from the beta-methylene hyperfine coupling and are caused by a distribution of molecular conformations, likely important for the biological function of this essential radical. The results demonstrate that model-based statistical analysis in combination with state-of-the-art spectroscopy accesses information hitherto beyond standard approaches.

Structural information at atomic resolution is essential for many areas of physics, chemistry, and biology. Electron spin resonance (ESR) or electron paramagnetic resonance (EPR) (1) belongs to a pool of methods, including nuclear magnetic resonance, optical spectroscopy, X-ray diffraction, and cryo-electron microscopy, which can deliver this information. For biophysical applications, the EPR-based electron–nuclear double resonance (ENDOR) technique is particularly suited for mechanistic studies involving endogenous paramagnetic centers, such as redox-active amino acid radical intermediates or metal ions and clusters, as it measures the hyperfine (hf) interaction between these centers and their surrounding magnetic nuclei (2), providing information on the protein architecture and establishing structural to functional relationships. The method has identified and characterized a wealth of paramagnetic intermediates in enzymes; for recent examples see refs. 3–7. Moreover, ENDOR spectroscopy has recently been proposed as a method of choice to measure molecular distances in the angstrom to nanometer range by 19F spin labeling (8), bridging a critical gap in pulse EPR-based distance measurements in biomolecules and complementing 19F-based NMR spectroscopy (9, 10).

ENDOR spectra of magnetic nuclei are usually recorded in the solid state (a frozen solution of a biological sample) at very low temperatures (T < 80 K) and thus display the full anisotropy of magnetic interactions, resulting in very broad lines. The signals are recorded as a change of the EPR electron-spin echo (11), which is sensitive to drifts in experimental conditions, particularly at low temperatures. Therefore, a model-free analysis of the observed signals constitutes one of the main challenges of ENDOR spectroscopy. In the last decade, progress in EPR instrumentation enabled the implementation of pulsed EPR spectrometers operating in the quasi-optical regime (12). We have recently reported a 1H ENDOR study of the essential in the wild-type (WT) subunit of Escherichia coli ribonucleotide reductase (RNR) (Fig. 1), using a commercial instrument (13) operating at 9.4 T/263 GHz. The strong orientation selection (hole burning) under high field/frequency conditions leads to sharpening of the spectra, strongly facilitating their analysis. We demonstrated the effect with the 1H spectrum of , illustrated in Fig. 1, where the sharp peaks were attributed to internal 1Hs of the radical. Additional broad resonances became visible, whose attribution presents additional challenges. From visual inspection alone, the significance of these signals could not be inferred given their shallow line shapes and the difficulty to separate them from unknown baseline distortions.

Fig. 1.

(A) Structure of and surrounding amino acids in the E. coli RNR WT- subunit from PDB 1mxr (22) in the reduced state (phenolic proton omitted). Highlighted are interactions with two protons at distances Å from the tyrosine oxygen. Protons were added with Pymol 2.2.2 in positions determined by the overall orientation of the amino acids. The O-O distance (3.8 Å) to the water molecule coordinating the proximal Fe ion is highlighted as well. (B) The 263-GHz 1H Davies ENDOR spectrum of adapted from ref. 13. Marked are resonances of the ring protons. Unknown broad features (asterisks) became observable only at 263 GHz. (Inset) Structure of a with labeling of the internal protons.

The structure and function of amino acid radicals, particularly in prototype biological machineries such as RNRs or photosystem II, have been the focus of several studies since the 1990s, due to their representative role in biological redox reactions and proton-coupled electron transfer (14, 15). In class Ia RNRs, the /diiron cofactor is essential for enzyme activity and thus for cell survival (16). In contrast to nonprotein model systems (17, 18), the conformation of s in proteins appeared much more constrained (19, 20). This finding was rationalized as due to the interaction with the protein environment that confers redox properties and thus governs biological function. Nevertheless, the mechanism of action of at the diiron cluster in class Ia RNRs remains puzzling. The current model for initiating nucleotide reduction requires that (E. coli numbering) is reduced in combination with protonation by a water molecule bound to the diferric cluster (16, 21) when the RNR subunits , are mixed with substrate and allosteric effector. In all structures of alone, however, the distance between the oxygen of the water and the oxygen of is too long for this step [3.8 Å in Protein Data Bank (PDB) 1mxr (22); Fig. 1]. Thus, the model requires a conformational change to shorten this distance, for which to date no experimental evidence exists.

Statistical methods are valuable to support the interpretation of biophysical experiments (23, 24). They have already been introduced in EPR spectroscopy for the analysis of distance measurements and distance distributions (25, 26). While adding uncertainty quantification to general spectroscopic practice, these works still hinge on assumptions about functions of interest, whether explicitly made through a Bayesian prior distribution or implicitly included in a choice of a penalization functional or thresholding procedure, some of which are hard to verify in practice. Additionally, ENDOR targets a different object of interest as it uses a different physical method. Therefore, a statistical analysis of ENDOR signals requires different methodology and has been less explored, possibly also due to a high incidence of systematic errors that are hard to model. In this contribution we aim at a statistical model that does not hinge on any specific assumptions regarding the functions of interest, in particular not on any that are hard to verify. Hence, our statistical model can be generalized to other biophysical experiments. In particular, our statistical approach includes careful statistical checking of distributional assumptions and thus lays a strong foundation for our statistical hypothesis tests. Recent advances in microwave (mw) technology, which allow for long-term signal stability and thus accumulation of large datasets, now provide the opportunity to gain additional information from quantitative ENDOR analysis. In this context, development of statistical methods that disentangle signals from noise and systematic error without having to rely on any assumptions regarding the spectrum and that report uncertainty estimates becomes mandatory.

The topic of this paper is a statistical analysis of ENDOR spectra and its representative application to assess the significance of broad spectral features observed in the 263-GHz 1H spectra of the E. coli RNR . We present a statistical model for the experimental data, which is used to extract the “most likely signal.” This model takes advantage of the information hidden in each individual scan (or batch) of the spectrum that is usually lost in the process of signal averaging. The treatment presents estimation of the uncertainty in the ENDOR spectra and permits subsequent statistical tests. The mathematical approach combined with spectroscopy of various isotopically labeled mutants of ultimately uncovered a distribution of conformations of , yielding insight into the mechanism of action of this essential radical.

Statistical Modeling of ENDOR Data

The ENDOR Signal.

ENDOR spectra display resonances of magnetic nuclei, which are hf coupled to the unpaired electron spin of a paramagnetic center. The observed doublet associated with one or more equivalent nuclei directly reports the hf splitting (Fig. 1). Resonances are recorded by detecting the effect of a radio frequency (RF) pulse on the intensity of the electron spin resonance signal, here a spin echo. Here, we discuss spectra recorded with a pulse sequence called Davies ENDOR, displayed in Fig. 2 (27); however, our statistical approach is not in any way limited to this pulse sequence. If the RF pulse excites a nuclear spin transition, the intensity of the echo is modified (for details refer to ). The ENDOR spectrum is recorded by integrating the transient of the total echo intensity as a function of the RF, which is usually stochastically varied after each sequence shot to attenuate heating (28) and nuclear saturation effects (29).

Fig. 2.

(A) Pulse sequence of the Davies ENDOR experiment. Nutation angles of the individual pulses are indicated with the respective phases as subscripts. (B) Quadrature detection (QD) scheme in the 263-GHz spectrometer, leading to a real () and an imaginary () part of the signal, detected as proportional to x and y magnetization, or the expectation values and , respectively, as explained in the text. Symbols containing crosses denote mw mixers, and triangles denote video amplifiers. A circle containing an arrow indicates a phase shifter.

To analyze and distinguish possible contributions in the ENDOR spectra, we first recapitulate some general features of the EPR and ENDOR signal detection, which is illustrated in Fig. 2. In the 263-GHz spectrometer, the EPR signal is down-converted to 9 GHz, where it is subsequently split into two channels and mixed with a reference signal to permit the standard quadrature detection scheme. When the relative phase of the reference with respect to the signal channel is properly adjusted, the electron spin primary echo (two-pulse echo) appears on the real channel (11). With the phase of the mw pulses defining the x direction, the echo signal is proportional to the y magnetization or to the expectation value of the spin operator . Due to the effect of mw pulses with finite bandwidth, a dispersive echo signal appears also on the imaginary channel; however, its integration over time results in a net zero signal.

In pulse ENDOR, discussed in detail by Gemperle and Schweiger (30) using the product operator formalism, only integration over the echo in the real channel gives a detectable ENDOR effect. The physical description does not predict a signal component proportional to , which would be detectable on the orthogonal channel. Relaxation effects can lead to anomalous ENDOR effects at large thermal spin polarizations (28, 31) but, to date, out-of-phase ENDOR signal components have not been reported. Nevertheless, an ENDOR signal is usually recorded in a complex-valued form (with an in-phase and an orthogonal component corresponding to a real and an imaginary part, respectively) to allow for proper phase reconstruction. In the following we mathematically describe standard processing of ENDOR data using averaging. We then propose our drift model, which yields improved signal-to-noise ratio by taking thermal drift in the measurement apparatus into account. Enabled by this model, we here introduce statistical uncertainty quantification to the analysis of ENDOR spectra in general, yielding confidence intervals and enabling testing of a range of specific hypotheses, e.g., as to the presence of a particular baseline distortion.

ENDOR Spectra via Signal Averaging.

Measurements of an ENDOR signal are carried out covering the radio frequencies indexed by . The measurements are repeated and we refer to these repetitions as batches and number them by . Further details are given in . Thus the entire measurement data are presented as a matrix with complex-valued entries . The batches are affected by experimental noise as we expect random variation of the signal from one repetition to the next. In the standard averaging approach, only the complex measurement vector averaging over batches according to Eq. is retained:We now mathematically formulate how the obtained average ENDOR signal is related to the ENDOR spectrum proper, which we denote as , to then extract the latter by statistical methods. We first decompose the average quadrature-detected signal into one component that is independent of the RF and produces a baseline and one component that is dependent on the RF signal and represents the desired ENDOR signal: . A normalization is then applied to facilitate comparison of spectra from different experiments, e.g., from different samples: We separate a root-mean-square average ENDOR signal strength , a positive real number, from relative variations of the signal by setting . While several normalization conditions are conceivable, herein we require 1) the mean of the normalized ENDOR signal to equal zero and 2) the sum of its modulus squared to equal one:

Compared to normalization conditions based on the observed maximum and minimum of the spectrum, this approach is less sensitive to noise. Since quadrature detection is used but a real-valued spectrum is expected, a phase correction [i.e., multiplication with for suitably chosen ] is required to rotate the normalized ENDOR signal so that to finally obtain , the desired ENDOR spectrum. Ideally, for the optimal rotation angle , the imaginary component would exhibit only noise and we describe the details of how to choose in . During standard data processing, one would first subtract the offset () as a baseline from both real and imaginary input signals and then rotate the phase of the spectrum such that the ENDOR spectrum is maximized. Statistically, this is equivalent to estimating the parameters , , and in the following statistical model:We point out that while is independent of the RF, it does not simply correspond to the EPR signal in the absence of the RF pulse, but may account for an additional offset contribution from the detector. The physical interpretation of these parameters including the offset contribution is schematically illustrated in . Inspection of Eq. reveals that, even before finding is attempted ( ), there are as many observations (namely real-valued observations) as there are parameters (namely 2 for , 1 for , for , and for reduced by due to the constraints in Eq. ). Therefore, the parameter estimates , , and (given in closed form in ) yield identically zero residualsHere and throughout, we use the hat to denote a statistical estimate, i.e., a value computed from experimental data, of a parameter of interest, e.g., estimating by .

Drift Model.

The averaging procedure according to Eq. works best when all observations have the spectrum in the same phase (corresponding to constant ) and magnitude (corresponding to constant ). However, an ENDOR experiment over many hours may exhibit signal drift due to thermal and mechanical instabilities. To investigate the possible drift behavior in detail we analyzed a 1H ENDOR data matrix consisting of batches of ENDOR signals of the perdeuterated 1,3-bis(diphenylene)-2-phenylallyl (BDPA) radical. Its 1H ENDOR spectrum consists of a single matrix line only and therefore facilitates straightforward decomposition of the ENDOR signal into the frequency range that contains an ENDOR effect and the frequency range that contains only baseline. We apply principal component analysis (PCA) separately for each batch and frequency range (see for details). The results indicate that the amplitudes averaged over either or , respectively, drift in a parallel fashion across batches. Furthermore, the phases of the leading principal components of for the two frequency ranges also drift across the batches. Additionally, the phases differ between the two frequency ranges by an offset of approximately , which stays constant across batches (). Overall, it becomes evident that there are baseline distortions and drift across the whole measurement for the parameters , , and in Eq. . Therefore, the standard averaging procedure is not suitable and can cause severe distortions. Moreover, from the PCA in the frequency region of the baseline we observe an additional signal orthogonal to the spectrum which is discussed in .

Based on these observations, we propose a refined model in Eq. for extracting the ENDOR spectrum directly from the complex data matrix :Here, is a Gaussian noise model with real and imaginary components with standard deviation and , respectively. Since the real and imaginary components of result from mixing a single measurement signal with a reference in a quadrature detector, we expect correlation between the noise in the two channels and explicitly allow for (partial) correlation expressed by the correlation coefficient :As evident from Eq. , there is a close connection to the averaging model: The scalar parameters and are replaced by complex vectors with components and (), respectively, to capture variation of the average signal and the amplitude and phase of RF effects across batches enumerated by . To ensure identifiability of the parameters , and , we retain the normalization conditions of Eq. .

We employ maximum-likelihood estimation to obtain complex-valued estimates ; our iterative algorithm for finding the maximum-likelihood estimator (MLE) is described in and can be viewed as a generalization of established iterative methods for finding the singular-value decomposition (SVD) of the matrix . Note that simultaneous access to all batches is required to obtain maximum-likelihood estimates as the full data matrix must be operated on.

Uncertainty estimation by bias-corrected bootstrap.

To estimate the uncertainty of the ENDOR spectra we employ the parametric bootstrap method (see as well as ref. 32 for a technical overview).

We choose a number of bootstrap samples (and hence ) for , which we use to perform bias correction () of the MLE of and . Bootstrap samples for (and hence ) obtained using the bias-corrected MLE are then employed to obtain confidence intervals for and for statistical hypothesis testing.

Shape of the Orthogonal Component .

Closer inspection of our previously reported ENDOR data at 263 GHz revealed that a systematic distortion is present in the orthogonal component . This systematic distortion is best illustrated in the previously mentioned BDPA 1H ENDOR dataset reported in . The approximate cosine shape of this disturbance—referred to as the wave—suggests that it results from a resonance phenomenon of some specific spectrometer component, yet its origin remains unknown. To the best of our knowledge, this distortion is present in all ENDOR spectra recorded at our 263-GHz spectrometer, both at room temperature and at low temperatures, including examples using other ENDOR pulse sequences. Based on these observations, we do not assume that this phenomenon arises as part of the ENDOR effect.

Specifically for the 263-GHz setup, we therefore adopt the following procedure to choose to best separate the desired ENDOR spectrum from the wave component. After computing the parameter estimate , the ENDOR spectrum of interest () is obtained by solving the following nested optimization problem: Find the optimal angle that rotates in such a way that the norm (i.e., the square-rooted sum of squared deviations from zero) of the resulting residual imaginary part is minimal. Here, the residual imaginary part refers to what is left of the imaginary part after subtraction of a cosine function whose amplitude, phase, and frequency are chosen to minimize the norm of the residual imaginary part. Following the optimization we obtain the estimated ENDOR spectrum and the estimated wave distortion as real and imaginary parts of the rotated , respectively. See for more details.

Comparison of Model Fit.

An algorithm for carrying out the estimation given the complex data matrix has been implemented in Python. This applies either the averaging procedure in Eqs. – or the drift model from Eqs. and subject to the constraints in Eq. and follows up by solving the optimization problem described in .

Averaging procedure.

Results of applying the averaging procedure to an example dataset (1H ENDOR of in E. coli RNR WT-) are summarized in Fig. 3. Fig. 3 presents the EPR spectrum, Fig. 3 provides the resulting ENDOR spectrum , and the resulting wave is shown in Fig. 3. Fig. 3 is used to critically assess the averaging model’s goodness of fit and is referred to in greater detail below. Indeed, as the complex data matrix is available, we can calculate and inspect the residualswhich, contrary to those of Eq. , are not identically zero and are plotted in Fig. 3. It is apparent from Fig. 3 that the residuals do not follow any of the well-known distributions and a plot of the frequency-averaged residuals in Fig. 3 reveals the origin of this problem: The residuals are not identically distributed across batches as their mean value changes over batches. This behavior is consistent with the experimental observation that phase and intensity of the echo signal drift with time.

Fig. 3.

Results of applying the averaging model to an example ENDOR measurement. (A) Echo-detected 263-GHz EPR spectrum of in E. coli RNR WT-. (B) The 1H ENDOR spectrum of obtained as in Eq. . (C) The approximately cosine-shaped distortion (the wave). (D) Kernel density estimate of residuals from Eq. . (E) Real component of frequency-averaged residuals plotted versus batch number. (F) Imaginary component of frequency-averaged residuals plotted versus batch number. Experimental settings for ENDOR are given in . The data were recorded in batches.

Drift model fit.

Results of fitting the drift model to the example dataset ( in E. coli RNR WT- from Fig. 3) are summarized in Fig. 4. Analyses of drift model fit for a representative subset of all datasets in the study are available in . Comparison with the results from the averaging procedure in Fig. 3 shows that both the ENDOR spectrum (in Fig. 3) and the wave (in Fig. 3) resulting from fitting the two models are similar except for improved signal-to-noise ratio (S/N). This is more clearly visible in . To see how this improved S/N results from taking the drift into account, consider the following hypothetical example: Assume the phase of the detector drifts at a constant rate covering one full period () over the course of the experiment; mathematically this corresponds to letting . This would lead to the signal covering the entire range of directions in the complex plane so that contributions pointing in opposite directions will cancel. Mathematically, this is seen by plugging Eq. into Eq. . Thus, such a drift would lead to complete cancellation and the spectrum obtained through the averaging procedure would be nothing but noise. In real data cases, phase drift is likely to lead to partial cancellation and thus accounts for varying degrees of improvement in signal-to-noise ratio (), depending on the amount of drift present in the data.

Fig. 4.

Results of applying the drift model to an example ENDOR measurement. (A) The 1H ENDOR spectrum of from in Eq. . (B) The distortion (the wave). (C and D) Real (black) and imaginary (red) components of and , respectively. (E) Kernel density estimate of residuals . (F) Dependence of the real (black) and imaginary (red) components of the averaged residuals on the radio frequency. (G) Standard deviation of the real (black) and imaginary (red) components of the residuals over radio frequencies versus batch number . (H and I) Empirical quantile—theoretical quantile plots for the real and imaginary components of , respectively. Shaded gray and red regions indicate pointwise 95% confidence intervals in A–D, very narrow in A and even more so in C.

The previously reported broad features in the spectrum (Fig. 1, asterisks) are reproduced using the drift model (Fig. 4) and visual comparison suggests that they are indeed part of the ENDOR spectrum and not a result of either random variations due to noise or the wave. Such an assessment is crucially supported by the availability of confidence intervals (displayed in Fig. 4) from the drift model, which are absent in the averaging procedure. Formal statistical hypothesis tests to further support our observation of broad features are carried out in Statistical Analysis of the Spectra. The eponymous drift of and is clearly observed in Fig. 4 , respectively. Detailed analysis of the residuals is provided in Fig. 4 . Correlation of the real and imaginary components of the residuals is clearly visible in Fig. 4 and persists after averaging over batches as is demonstrated in Fig. 4. Both components of the residuals are remarkably well modeled by normal distributions with mean zero. We visualize this using a normal quantile–quantile plot that plots the sorted sample values of the residuals versus the corresponding quantiles under the normal distribution: Deviations from the diagonal indicate differences between empirical and theoretical distributions. As Fig. 4 shows, these deviations appear very small. Indeed, Kolmogorov–Smirnov tests for goodness of fit of the variance-standardized residuals for normality yield P values of and for the real and imaginary components, respectively, each of a sample size of 108,000. Therefore, the null hypothesis that the residuals follow a normal distribution cannot be rejected. This is a remarkable finding given the sample size and demonstrates the superior fit of the drift model both compared to the averaging model and judged against commonly encountered goodness of fit in statistical practice.

Application: Testing for Broad Resonances Reveals Conformational Distributions

The drift model enables extraction of spectra with an estimate of their associated uncertainty. These estimated uncertainties can, in turn, be used to decide whether broad features in ENDOR spectra, as discussed in the Introduction in the case of , are statistically significant. To this end, we propose two different statistical tests.

Flatness Testing.

The flatness test tests the null hypothesis that the spectrum is flat over a specified frequency region of interest , mathematically speaking that it equals an unknown constant, quantifying the evidence against flatness by the P value . To accomplish this, we define a test statistic based on the summed squares of the differences of from its average value over the region of interest according to Eq. :where denotes the number of radio frequency indexes in . We recall that whenever the P value is below the chosen significance level, the null hypothesis of flatness is rejected. Here, is the observed value of the test statistic whereas denotes the random variation of the test statistic under the null hypothesis, i.e., assuming that the deviation from flatness arises only from the modeled noise. We then use a bias-corrected bootstrap sample from the distribution of to estimate the P value resulting from the test as well as a -based asymptotic approximation . We inspect the standard significance levels of , and Benjamini–Hochberg corrected for multiple testing (a total of 56 flatness tests are conducted herein). Details are given in and all values are reported in .

Difference Testing.

To assess the effect of partial deuteration when comparing the ENDOR spectra of the and of the -2- samples, a test ascertaining whether an observed difference between and is statistically significant is desirable. To avoid comparison of P values from two individual flatness tests for the individual spectra, we propose a difference test considering both measurements and their uncertainties simultaneously. The method’s null hypothesis is that a straight line least-squares fitted to the difference between and in the frequency region of interest of the spectrum has slope zero; i.e., it is horizontal. We again employ bias-corrected bootstrap samples of both spectra. Details on how the test computes the resulting P value as well as its Gaussian-based asymptotic approximation are given in . Since we will perform a total of 16 difference tests, we employ a Benjamini–Hochberg correction for multiple testing.

We note explicitly that our slope-based approach has high power for detecting monotonous differences covering a broad spectral range and is therefore well suited to test for the linear component of the suspected broad features in the spectra. The region of interest has to be chosen with care, as the test inherently has low power if positive and negative slopes of the difference between the spectra are canceled out on average within .

Statistical Analysis of the Spectra.

To analyze the origin of the broad features in the spectrum (Fig. 1) with the aid of statistical methods, three 2H-isotope–labeled protein samples were prepared and investigated. All protein mutants were catalytically active and contained about 1 per (Incorporation of D2-Y, D4-Y, and D7-Y in the E. coli RNR Protein). Echo-detected EPR spectra at 263 GHz of the four samples, i.e., (WT), -2- (D2), -2- (D4), and -2- (D7), are reported in . The spectra are clearly dominated by -factor anisotropy and their comparison nicely demonstrates the expected effect of partial deuteration: The large hf interaction with one of the -methylene protons, which produced the doublet structure visible in Fig. 3, in the EPR spectrum of , is not observed in the -deuterated analogs (D2 and D7).

The 1H Davies ENDOR spectra of the four samples were typically measured at five different resonance positions in the 263-GHz EPR line, of which only the position denoted is shown in Fig. 5—all other positions for WT and D2 data are provided in . Note that the spectra in Fig. 5 and in are obtained from the drift model as opposed to the spectrum in Fig. 1, where the averaging model was employed. The five resonance positions, at which the ENDOR spectra were typically recorded, are denoted by , mT, , mT, and and they are illustrated in . All displayed spectra were analyzed with the drift model and frequency-wise 95% confidence intervals were computed using bootstrap samples to obtain error bars.

Fig. 5.

The 263-GHz 1H Davies ENDOR spectra (recorded at EPR resonance position ) of the four protein samples with different 2H content of and results of flatness tests—the deuteration scheme is given in Left Insets. (A) Wild-type data, (B) D2 data, (C) D4 data, (D) D7 data. The -approximated P values obtained from the flatness test are included and significance (after Benjamini–Hochberg [BH] correction across all flatness tests in the paper) is indicated by *** for P values achieving the 0.001 BH-corrected significance level and by * for P values achieving the 0.05 BH-corrected significance level. The orange center line indicates the 1H Larmor frequency.

The Davies ENDOR spectra show the expected disappearance of the ring 1 resonances in the D4 and D7 samples (Fig. 5 at MHz). The broad features (Fig. 1) seem absent in spectra in Fig. 5 ; however, ascertaining the presence of the broad features in the spectrum in Fig. 5 by simple inspection is difficult and requires statistical analysis. Additionally, we note that surprisingly little variation is observed among isotopologues in the central area of the spectra (white central areas in Fig. 5). Thus, the lines attributed to ring 1 resonances (±1 MHz) are only reduced in intensity. The large number of resonance lines in the spectrum of the D7 sample is the first clear indication of multiple contributions due to the protein environment.

We apply the flatness test to all ENDOR spectra in those spectral ranges lacking resonances due to known hf couplings of the ring protons. The results are depicted in Fig. 5. For all spectra shown, two or four regions (I to IV) of interest were examined: the left and right marginal regions (regions I and IV in Fig. 5) as well as the intermediate regions, i.e., the parts where broad features are suspected (regions II and III for spectra in Fig. 5 ). For all spectra in Fig. 5 , these regions were defined using the same offsets relative to the proton Larmor frequency. Regions III and IV were chosen to be mirror images of regions II and I, respectively, following the approximate mirror symmetry of the ENDOR spectra with respect to the proton Larmor frequency. The resulting P values and an indication of Benjamini–Hochberg-corrected significance (by asterisks) are also included in Fig. 5. We observe that in all marginal regions, the P values are above the standard 0.05 significance level (even before correction for multiple testing), consistent with our null hypothesis of flatness. Only in the intermediate regions of spectra in Fig. 5 do we observe P values that are surprisingly small. Since we have already confirmed goodness of fit of the drift model (), we can rely on the P values and cannot attribute the surprisingly small P values to model misfit. Therefore, either an extremely low probability event yields the surprisingly small P values or the null hypothesis (flatness) is wrong in the regions concerned.

We thus conclude that statistically significant deviations of the spectrum from flatness are found in the intermediate regions for the WT and the D4 isotopologue. Similar results are obtained in an analogous fashion for the remaining spectra at other orientations as summarized in .

We apply the difference test to the WT-D2 pair of ENDOR spectra and the small P values indicate statistically significant differences in all intermediate regions (Fig. 6). In contrast to this, the difference test does not reject the null hypothesis (i.e., no difference in the spectra) in any marginal regions. Since the test has sufficient statistical power to detect differences (as seen in the intermediate regions), it is reasonable to assume that the spectra are in fact not different in the marginal regions.

Fig. 6.

Result of the difference test for the ENDOR spectra of the wild-type (black) and D2 (red) samples recorded at the five different EPR resonance positions illustrated in . For each resonance position (, mT, , mT, and ), the respective plot (A, B, C, D, and E) shows the ENDOR spectra aligned horizontally and vertically for best agreement. The regions are indicated by colored shading and the approximate P values are provided along with a Benjamini–Hochberg-corrected indication of significance, taking into account all difference test P values (*** corresponding to the significance level and * corresponding to the significance level). P values are based on pairs of simulated spectra. The null hypothesis of identity of the spectra is rejected in the regions of interest II and III for resonance positions , mT and .

Thus, while the flatness test confirms that the WT spectrum is not flat in the intermediate regions, the difference test confirms that the WT and the D2 spectra are different from one another in these regions. In the marginal regions, on the other hand, D2 and WT spectra seem to both be flat and show no differences from one another. This supports the attribution of the differences in the intermediate regions to deuteration because other potential variations between measurements (e.g., in detector offsets, RF heating, etc.) would likely have caused statistically significant differences or nonflatness in other regions, too.

Together, both statistical tests therefore unequivocally confirm the presence of the broad features in the WT and D4 spectra suspected upon visual inspection and support their attribution to the -methylene protons.

Evidence for a Distribution of Hyperfine Couplings.

Given the association of the newly discovered broad features in the ENDOR spectra of with the -methylene protons ( and ), the question arises: What is the origin of their peculiar line shape? Considering that one of these two protons () produces the dominant, almost isotropic hf splitting (ca. 60 MHz) observed in the EPR spectrum (Fig. 3), the smaller interaction to is likely responsible for the broad features. In our previous analyses (13, 33) we attributed the coupling to the sharp inner peaks at (0.5 MHz; Fig. 5). However, from the D7 spectrum in Fig. 5 it becomes evident that most sharp peaks in the central part of the 1H ENDOR spectrum (within 2 MHz) are conserved upon deuteration and must arise from the protein environment. This finding together with the establishment of the new broad features caused us to reexamine the 1H assignments.

Previous density functional theory (DFT) studies (18, 34) have predicted that the isotropic part of the and couplings strongly depends on the ring dihedral angle (defined by the , , , and atoms; Fig. 7, Inset) of the tyrosyl radical. To illustrate this dependence, we have plotted the hf tensor components of and as a function of this dihedral angle (Fig. 7) as calculated by DFT on a small tyrosyl radical model system (). As expected, the hf couplings are nearly isotropic. We note that for orientations of close to the ring plane (corresponding to a dihedral of ) its hf coupling is nearly zero. Nevertheless, very subtle reorientation of the radical can strongly affect the coupling. Thus, we presume that the broad features in the ENDOR spectrum arise from a distribution of ring dihedrals in . A distribution of couplings due to structural heterogeneities has been previously reported from 2H-electron spin echo envelope modulation studies of tyrosyl radicals in frozen solutions produced by light irradiation (17). However, to our knowledge, distributions of couplings specifically have not been reported for protein .

Fig. 7.

Distribution of ring dihedral angles and simulations of 263-GHz 1H ENDOR spectra of wild-type . (A) Principal components of the hyperfine tensors for the s (green and blue) as a function of the ring dihedral, as obtained from DFT calculation (, triangles; , circles; , crosses). The range included in the simulation is highlighted in gray. (Inset) Overlay of model structures with three different ring dihedral angles, corresponding to the highest-probability structure (, solid) and those corresponding to 0.5% probability (shaded). (B) Normalized probabilities of the ring dihedral angles of the simulated distribution. (C–E) Representative 263-GHz ENDOR spectra of (black) at those EPR resonance positions for which the broad features are clearly visible, i.e., (C), mT (D), and (E), and spectral simulations (red) with values of Table 1. The distribution of the resonances is highlighted in green.

Table 1.

Hyperfine coupling parameters used for the EPR and ENDOR spectral simulations

Nucleus	A_1/MHz	A_2/MHz	A_3/MHz	α/°	β/°	γ/°
H^3	−27.2	−7.9	−19.8	25	0	0
H5	−27.2	−7.9	−19.8	155	0	0
H^2	4.9	7.7	1.7	−25	0	0
H^6	4.9	7.7	1.7	−155	0	0
Hα	1.4 (1.2 to 1.6)	−2.1 (−2.9 to −1.8)	−1.8 (−2.2 to 1.6)	83 (73 to 88)	106 (95 to 82)	−40 (−35 to −43)
	1.2	−1.0	−4.1	0	20	−60
H^β_1	60.2 (54.7 to 61.6)	53.7 (48.3 to 54.9)	55.5 (49.9 to 56.9)	3 (43 to −19)	33 (38 to 36)	−1 (−35 to 18)
	59	55	56	0	0	0
H^β_2	12.8 (6.9 to 21.9)	7.0 (1.3 to 15.9)	7.5 (1.8 to 16.5)	−81 (−94 to −77)	135 (172 to 99)	66 (80 to 62)
	8	−1	−1	7	−4	−70

As values for , , and depend on the ring dihedral angle, values for the highest-probability conformation () are given along with themaximum and minimum values (in parentheses). For comparison, values for , , and from our previous report (13) are included in italic font. For the ring protons previously reported values (13) were used. Euler angles follow the z, convention and describe the rotation of the A tensor to the tensor. Positive values correspond to clockwise rotation when viewed along the rotation axis. The values were 2.00915, 2.00460, 2.00225 (13).

To validate the assignment of the broad features to , we performed spectral simulations combined with DFT predictions for the couplings. For simplicity, we describe the distribution of ring dihedral angles using a normal distribution with only two parameters (mean and half-width ) (Fig. 7). These two parameters were adjusted to reproduce the new broad features in the ENDOR spectra as well as the resonances of . All ring 1H couplings were adopted from the previous reports in the literature (13, 20, 33, 35). Representative ENDOR simulations for values of and and including only the internal H hf couplings are depicted along with the experimental spectra in Fig. 7 . The relatively narrow distribution (width of ) is sufficient to spread the hf values of over a range up to 22 MHz and to qualitatively reproduce the line shape of the broad absorptions. Interestingly, the coupling covers a range of ca. 10 MHz and is less sensitive to the conformational distribution. This is due to the dependence of the hf coupling on the dihedral angle combined with the particular mean value of the dihedral angle in (Fig. 7 ). Nevertheless, the distribution significantly improves the simulation of the ENDOR resonances (). A more elaborate mathematical procedure to analyze the distribution in terms of possible radical conformations, as for instance proposed in ref. 17, is in progress but goes beyond the scope of this work. It will require a global statistical analysis of the ENDOR spectra including the wider range of the resonances.

All simulation parameters and comparison with previous assignments are reported in Table 1. We note that the distribution of dihedral angles also affects the coupling of , albeit to a lesser extent in comparison with and . This coupling is included in the simulation and gives rise to signals close to the matrix line. Simulations were also performed assuming only one conformation of the , i.e., the most probable one with dihedral angle . This yielded significantly sharper and more structured features not consistent with the spectra (). Notably, the EPR spectrum is also consistent with the observed distribution, which produces a line broadening of less than 0.3 mT. The broadening is within the intrinsic EPR line width at 263 GHz (13) and therefore not visible in the EPR spectrum, as confirmed by the observable splitting (Fig. 3) as well as by the EPR simulation performed with the same parameters as the ENDOR simulations ().

Inspection of various structures of the E. coli RNR determined by X-ray diffraction (PDBs 1rib, 1pfr, 1mxr, 5ci4) (22, 36–38), in which is protonated and reduced, reveals that in the crystals the ring dihedral angle of ranges from to . The g tensor was found to be reoriented by to compared to the side-chain orientation in the reduced form (22, 39). From initial ENDOR studies at low EPR frequencies, the ring dihedral in frozen solution was estimated to be around (20, 35). All these values are in principle close to or within the range of the distribution identified here, but taking them individually, they would not be consistent with the spectra. We thus conclude that the dihedral angle approximately follows a normal distribution with mean value and strictly positive standard deviation . Our results give evidence that in its protein pocket has some conformational flexibility that has never been observed before, despite substantial spectroscopic efforts since the 1980s.

Conclusion and Outlook

We have demonstrated that statistical analysis is able to model ENDOR spectra, taking into account the information hidden in multiple individual experiments (batches) that is lost in the usual process of signal averaging. This permits us to reconstruct the ENDOR signal with confidence intervals, valuable for statistical hypothesis testing, without having to rely on assumptions on the unknown spectrum. This approach can distinguish experimental baseline distortions from subtle absorption features, such as less populated conformations or distributions that are important to understand biological function. Moreover, the approach is general as it can in principle not only be adapted to ENDOR spectra at any microwave frequency but also expanded for other types of spectroscopic methods. Specifically, for ENDOR the availability of a statistical model that fits the measured data well allows uncertainty quantification and opens the way for more sophisticated data modeling and interpretation, for instance to systematically achieve better agreement of simulated with measured spectra.

Applying this approach to data collected for the essential tyrosyl radical in the wild-type -subunit of E. coli RNR as well as three additional mutants, in which the radical was partially deuterated, allowed the unambiguous distinction of broad spectral features from the baseline. We verify the presence of previously unrecognized hf couplings using two statistical tests, a flatness test and a difference test. The result establishes that the coupling is contributed by the and is consistent with a distribution of conformations. Such a distribution has not been observed before due to the difficulties of analyzing broad ENDOR resonances in the powder pattern of multiple contributing nuclei.

A wealth of studies have been reported to characterize the detailed electronic structure and the role of tyrosyl radicals in proteins, with the of E. coli RNR being a prototype. The observation of a conformational distribution of strongly suggests that very small conformational changes might be sufficient to couple and uncouple this radical from the adjacent diiron cofactor and from its hydrogen bond donor, likely a water molecule at this cluster. Thus, our finding draws attention to the role of less populated states in such a highly regulated mechanism, that is, the reduction of an unusually stable protein radical required for enzyme activity and cell survival. Future global analysis of the 1H ENDOR spectra of including the spectral range using our statistical approach will potentially allow extraction of detailed information on the radical conformation and on the protein architecture. Such experiments will require trapping the catalytically active complex of the enzyme using known constructs of E. coli RNR subunits , substrate, and effector (4, 40). The experiments are ongoing but challenging and beyond the scope of this paper. Nevertheless, the approach presented here paves the way for this kind of spectroscopic analysis in RNR but also in many other biological systems.

Materials and Methods

The 4-(2-hydroxyethyl)-1-piperazineethanesulfonic acid (HEPES), ampicillin (Amp), M9 minimal salts, L-arabinose, and streptomycin sulphate were purchased from Sigma-Aldrich. The deuterated tyrosines (-Y, 2,3,5,6--Y and 2,3,5,6,-Y) were purchased from Cambridge Isotope Laboratories with a degree of deuteration . Buffer A comprises 50 mM tris(hydroxymethyl)-aminomethane (Tris), 5% glycerol (vol/vol), and 0.5 mM phenylmethylsulfonyl fluoride (PMSF) at pH 7.6. Assay buffer contains 50 mM HEPES, 15 mM , and 1 mM ethylenediaminetetraacetic acid (EDTA) at pH 8. Wild-type E. coli RNR was prepared as previously reported (41).

Incorporation of -Y, -Y, and -Y in the E. coli RNR Protein.

The pBAD-nrdB plasmids were transformed into Top10 competent cells and grown overnight at 37 C on lysogeny broth agar plates supplemented with 100 g/mL Amp. A 10-mL starter culture containing the Amp was inoculated with a single colony and grown until saturation (ca. 12 h at 37 C) in glycerol minimal medium (0.05% [wt/vol] glucose, 0.2 g/L of all amino acids except Tyr and Cys), 1 mM of the respective, partially deuterated tyrosine, 100 g/mL Amp, and 1× heavy metals. The starter culture (10 mL) was used to inoculate the same medium supplemented with glyphosate (1 g/L). The culture was grown at 200 rpm, 37 C until the optical density at 600 nm reached 0.6 and was induced with 0.05% arabinose. The culture was grown for an additional 4 h and harvested by centrifugation (4,000 × g, 10 min, 4 C). Total yields of 7.6 g (-Y), 6.2 g (-Y), and 5.8 g (-Y) of wet cell paste were obtained. The proteins were purified using diaminoethyl cellulose and Q-Sepharose chromatography methods as previously described (41), providing the protein containing the variously deuterated tyrosines (34 mg/g of cell paste for -Y, 43 mg/g of cell paste for -Y, and 45 mg/g of cell paste for -Y). Incorporation of the partially deuterated amino acids was verified by electron spray ionization mass spectrometry. Expected mass gains for 16 tyrosines in each were 32, 64, and 112 unified atomic mass units (u) for -Y, -Y, and -Y, respectively. Observed were 32.5, 63.8, and 106.5 u, the latter likely reflecting some washout of the -proton during incorporation. Protein was diluted by addition of a portion of glycerol-free deuterated buffer A (0.5 mL) and the centrifugation was repeated. After six centrifugations in total, the protein was taken up in 0.4 mL of glycerol-free deuterated buffer A. Concentrations were determined by ultraviolet/visible (UV/Vis) ( = 131,000 L) and the samples were diluted further to the final measurement concentration of 300 M. Radical content was determined by UV/Vis spectrometry using the dropline method (42) to be ca. 1 per dimer. Protein activity was measured for all samples using the spectrophotometric assay (42). The assay contained, in a volume of 300 L, 0.75 M 2, 0.15 M , 1 mM cytidine diphosphate (CDP), 3 mM adenosine triphosphate (ATP), 30 M thioredoxin (TR), 0.5 M thioredoxin reductase (TRR), and 0.2 M NADPH in assay buffer. The reaction at 25 C was initiated with addition of and the absorbance monitored at 340 nm ( = 6.22 ) on a Varian Cary 3 spectrophotometer. The observed absorbance change at 340 nm was corrected for the background change prior to addition of . The average of three independent measurements for each sample provided the following activities: 7,330 90 nmolmin−1 for -Y-, 6,840 160 nmolmin−1 for -Y-, and 7,100 90 nmol-Y-.

The 263-GHz ENDOR Spectroscopy.

Protein solutions for ENDOR were prepared at 0 C, loaded into Suprasil quartz capillaries (outer diameter, 0.33 mm; inner diameter, 0.2 mm), and frozen in liquid . For sample insertion the ENDOR resonator was precooled in a liquid bath and after mounting of the sample was transferred into the precooled cryostat (Oxford Instruments) at 80 K. Temperature of the cryostat was adjusted using a low-temperature helium gas flow and the ITC-5035 temperature controller (Oxford Instruments).

The 263-GHz spectrometer (Bruker ELEXSYS E780) and the ENDOR resonator (Bruker E9501510) have previously been described in detail (13). RF pulses are generated by a two-channel RF pulse-forming unit (Bruker DICE-II) and amplified by a 125-W RF amplifier (Amplifier Research; model 125W1000). Amplified RF power is conducted by a 50-âĎę semirigid coaxial cable (type RG 214U) to the ENDOR coil and terminated after the coil by a 50- load. To mitigate possible heating effects (43) and ENDOR signal saturation (29), a “stochastic” acquisition mode was used in contrast to a linear sweep of the RF range. Prior to each ENDOR measurement, an echo-detected EPR spectrum was recorded to determine the excitation position for ENDOR. Comparison of echo-detected EPR spectra before and after the ENDOR measurements indicated negligible field drifts of about 0.2 mT during the measurement time (ca. 22 h). Davies ENDOR spectra of the protein samples were measured at 5 K using microwave pulses of equal amplitude. If not otherwise indicated in the figure legends, the preparation pulse length was set to 200 ns, while for spin echo detection a corresponding sequence with a value of 350 ns was used. RF pulse length was 60 s, and the repetition time of the sequence was 12 ms at one shot per point. Acquisition times were typically between 10 and 20 h.

Acquisition of ENDOR Spectra.

For statistical purposes, the availability of every single measurement is desirable but this exceeds the on-board storage capacity of the detection unit in the case of several thousand repetitions. As a compromise, several experiments are accumulated (i.e., averaged) in what we herein term a batch to balance averaging time per batch and memory requirements. In all experiments described herein, 50 individual experiments were accumulated in less than 5 min per batch and we typically used batches. A data matrix was recorded at five representative excitation positions in the EPR lines to cover most molecular orientations according to their g values ().

DFT Calculations.

The dipeptide Ac-Tyr-Me was chosen as a model system as it allows meaningful use of the dihedrals (defined by the peptide backbone atoms , N, , and C) and (defined by N, , C, and ) to constrain the overall amino acid orientation to the values observed in the solid-state structure determined by X-ray diffraction (PDB 1mxr) (22). Following the results of a single-crystal EPR study of the system, the backbone dihedral (defined by atoms C, , , and ) was constrained to 169. and the ring dihedral was varied in steps from to . The radical structure was then optimized (Orca 4.0.1.2) (44) without any further constraints using the B3LYP functional in conjunction with the def2-TZVPP basis set (45). Normal self consistent field convergence criteria and a dispersion correction (D3BJ) were used (46). Subsequently, hf couplings were calculated using the same parameter combination.

Spectral Simulations.

The EPR and ENDOR spectra were simulated simultaneously using a home-written algorithm in MATLAB that computes EPR and ENDOR transitions in the high-field approximation (previously described in refs. 13 and 33). The algorithm was adapted for time efficiency. The probability of the various conformations was calculated from the values and describing the normal distribution. For each conformation with a significant probability (0.5%), the EPR spectrum and ENDOR spectra at the five positions in the EPR line were simulated. To increase resolution, dihedral angle steps of 0. were chosen and the associated hf coupling parameters of the -protons and the -proton were obtained from spline interpolation of the DFT data. For the four ring protons, possible inequivalence and dependence on the ring dihedral were disregarded. For the hyperfine couplings of and , an empirical scaling factor of 0.955 was applied to the principal values of the DFT-derived hf tensors to improve the overall agreement with the experimental ENDOR line positions. The deviation of this factor from unity likely represents the accuracy limits of the DFT method, as the hf coupling size depends strongly on the degree of spin delocalization. The corresponding spectra were then averaged over all ring dihedral values, weighted by their corresponding probability. The simulated EPR spectra were convolved with a Gaussian line shape of mT width, and the ENDOR spectra are simulated with Lorentzian line shapes with a width of approximately 20 kHz for a single transition.