On the use of 3J-coupling NMR data to derive structural information on proteins

Values of 3J-couplings as obtained from NMR experiments on proteins cannot easily be used to determine protein structure due to the difficulty of accounting for the high sensitivity of intermediate 3J-coupling values (4–8 Hz) to the averaging period that must cover the conformational variability of the torsional angle related to the 3J-coupling, and due to the difficulty of handling the multiple-valued character of the inverse Karplus relation between torsional angle and 3J-coupling. Both problems can be solved by using 3J-coupling time-averaging local-elevation restraining MD simulation. Application to the protein hen egg white lysozyme using 213 backbone and side-chain 3J-coupling restraints shows that a conformational ensemble compatible with the experimental data can be obtained using this technique, and that accounting for averaging and the ability of the algorithm to escape from local minima for the torsional angle induced by the Karplus relation, are essential for a comprehensive use of 3J-coupling data in protein structure determination.


Introduction
Structural information on proteins can be derived from a variety of observable quantities, such as X-ray diffraction intensities, NMR, CD, Raman or infrared spectra to mention a few . For proteins in aqueous solution, NMR measurements provide the highest information density, i.e. ratio of the number of independent, measured values of observable quantities for a molecule and the number of independent molecular degrees of freedom.
Quantities observable using NMR techniques are 3 J-couplings, chemical shifts, nuclear Overhauser enhancement intensities (NOEs), S 2 order parameters and residual dipolar couplings (RDCs). The measured value of such an observable quantity Q is an average 〈Q〉 space,time of Q over the molecules (space) in the test tube and over a time window determined by the particular measurement technique. This means that 〈Q〉 constitutes an average over a Boltzmann-weighted set, i.e. a statistical-mechanical ensemble, of configurations ⃗ r . The weights are proportional to exp(−V(⃗ r)∕(k B T)), where V(⃗ r ) indicates the energy of a molecular configuration or structure ⃗ r, k B is Boltzmann's constant and T is the temperature.
If an observable quantity Q is dependent on the molecular configuration ⃗ r , one may try to derive an expression or function Q(⃗ r) that approximates the relation between Q and ⃗ r, which expression may then be used to derive molecular structures that are compatible with measured values of Q, i.e.
1. Q exp values are subject to uncertainty or error. 2. It is not possible to fully account for averaging over space and time inherent in the experimental measurement: inversion of the averaging operation in Eq. (1) is impossible. 3. For most bio-molecular systems the number of independent Q exp values available is much smaller than the number of degrees of freedom of the system. This means that the structure determination problem is underdetermined and can only be addressed by using an atomic model, i.e. a function V(⃗ r) specifying likely structural parameters (e.g. bond lengths and bond angles) of a system. The function V(⃗ r ) may yield low-energy values for configurations that are physically most reasonable. The fewer Q exp values are available, the larger the influence of the choice of molecular model and interaction function V(⃗ r ) and its (in)accuracy on the generated structures will be. 4. The function Q(⃗ r ) is not known or its accuracy is uncertain or low. 5. The inverse function ⃗ r(Q) of the function Q(⃗ r ) may not exist, or if it does, it may be multiple-valued, as in the case of the Karplus relation or function (Karplus 1959(Karplus , 1963 linking a 3 J-coupling (Q) to a torsional angle θ(⃗ r ) in a molecule, 6. The generation or sampling of molecular configurations ⃗ r must be biased, i.e. guided towards those that are (on average) compatible with Q exp . This is particularly challenging when the inverse function ⃗ r(Q) of the function Q(⃗ r ) is multiple-valued. The inverse function θ( 3 J) of Eq.
(2) is multiple-valued. For a range of 3 J-values there are four different θ-values satisfying Eq. (2).
Although 3 J-couplings are relatively easily measurable, their use to derive molecular structure is hampered by all six mentioned problematic aspects of procedures to derive molecular structure from measured Q exp values: (i) Small values of 3 J-couplings, such as occurring for 3 J HNHαcouplings in helical structures, are not easily precisely measured; (ii) The averaging of 3 J-couplings may cover very long time periods and is not limited by the rotational tumbling time of the molecule, as for NOEs; (iii) 3 J-couplings can only be measured for particular parts of proteins, backbone 3 J HNHα -couplings and side-chain 3 J HαHβ -couplings only address values for the φ and χ 1 torsional angles respectively; (2) 3 J( ) = a cos 2 ( ) + b cos( ) + c, Although a range of other 3 J-couplings can be measured, particularly with isotopically labelled samples, in general these show a smaller range of values and often have less well defined Karplus relation parameters. They are therefore not so useful for structure determination (Wang and Bax 1996;Perez et al. 2001); (iv) The parameters a, b and c of the Karplus relation Eq. (2) are of empirical nature and their values are commonly derived by matching X-ray diffraction derived crystal structures ( ⃗ r ) to solution NMR measured 3 J-couplings (Q) for a set of (protein) molecules (Schmidt et al. 1999;Schmidt 2007). This leads to an uncertainty of 1-2 Hz in Eq. (2) (Dolenc et al. 2010;Steiner et al. 2012); (v) The inverse function θ( 3 J) of Eq. (2) is up to fourfold multiple-valued.
A few years ago, a procedure has been proposed to effectively handle the multiple-valuedness of θ( 3 J) by searching for all possible local minima of the biasing restraining function V k J,restr (J k (θ k (⃗ r(t))), 〈J k 〉 t ; K Jr , N le , J k 0 , ΔJ fb ) that keeps the value of 〈J k 〉 of the kth 3 J-coupling close to J k 0 = J k exp ). This restraining function has a flat-bottom on the interval [J k 0 − ΔJ fb , J k 0 + ΔJ fb ] and is beyond this interval harmonic, see Fig. 1 of Smith et al. (2016) with ΔQ h = ∞. Time is denoted by t and K Jr is the overall weight or force constant of the penalty or restraining function. The local-elevation parameter N le defines the number of intervals N le of the torsional angle θ around the θ i 0 -values and their widths used in the local-elevation search and sampling algorithm (Christen et al. 2007). The restraining function, biquadratic in J k (θ k (⃗ r(t))) and 〈J k 〉 t , only yields a non-zero energy and restraining force when both, J k (θ k (⃗ r(t))) and 〈J k 〉 t , adopt values on the same side outside the flat-bottom region. When either the instantaneous value J k (θ k (⃗ r(t))) or the average 〈J k 〉 t satisfies J 0 = J exp within a given uncertainty ΔJ fb , there is no need for restraining. The local-elevation searching and sampling technique (Huber et al. 1994) is used, in which the potential energy at already visited parts of configuration space is raised in order to avoid repetitive sampling of the same parts of configuration space in a simulation. The basic idea of local-elevation structure refinement based on adaptive restraints (Christen et al. 2007) is that whenever the simulated average 〈J〉 t of the 3 J-coupling and the current value J(t) of the 3 J-coupling do not match the measured target value J 0 = J exp to within an uncertainty ΔJ fb , the force constant or weight ki (t) of the penalty or restraining function V k J,restr (J k (θ k (⃗ r(t)))), 〈J k 〉 t ; K Jr , N le , J k 0 , ΔJ fb ), acting on the current value θ k (t) of the torsional angle θ k corresponding to the 3 J-coupling, is increased. The restraining potential-energy term of a given 3 J-coupling J k is a sum of N le terms, local with respect to a particular value θ i 0 of θ k , with the weight ki (t) of the ith penalty function changing during the simulation according to for an attractive (both, J k (t) and 〈J k 〉 t larger than J k 0 ) 3 J-coupling restraint, and according to for a repulsive (both, J k (t) and 〈J k 〉 t smaller than J k 0 ) 3 J-coupling restraint. The Kronecker delta δ is defined using finite differences, and the Heaviside step function H(x;x 0 ) is defined by This means that the torsional angle θ is pushed away from the value θ(t) when both the current J(t) and the (5) averaged 〈J〉 t values are deviating more than ± ΔJ fb from the target value J 0 . In this way the whole range of θ values is searched for those values that yield 3 J-couplings close to J exp (Allison and van Gunsteren 2009). The time-averaged value 〈J〉 t of J(t) is commonly (Torda et al. 1989) exponentially damped in the restraining function V k J,restr (J k (θ k (⃗ r(t))), 〈J k 〉 t ; K Jr , N le , J k 0 , ΔJ fb ), Eqs. (6) and (7), in order to avoid that the restraining force progressively approaches zero with time.
This procedure to refine protein structure using measured 3 J-couplings was applied ) to sets of 95 backbone 3 J HNHα -couplings and 62 side-chain 3 J HαHβcouplings (58 stereo-specifically assigned 3 J HαHβ -couplings plus those of Glu 7 and Arg 45 were included in the set of restraints) measured for hen egg white lysozyme (HEWL) (Smith et al. 1991). It was concluded that (i) the weight or force constant K J,restr of the restraining function should be chosen as small as possible while large enough to bring the average 3 J-couplings close to the target values in order to approximate as good as possible the properly (using V(⃗ r )) Boltzmann-weighted conformational probability distribution, (ii) the averaging time τ J should match the experimental one but should not be larger than 1/10 of the total simulation time in order to secure sufficient statistics when averaging, and (iii) the flat-bottom Δ 3 J fb should represent the uncertainty or inaccuracy of the Karplus relation 3 J(θ).
In the present article the biquadratic time-averaging localelevation restraining (BQ-TA-LER) method ) is further investigated using NMR and X-ray data on HEWL, particularly concentrating on side chains whose 3 J HαHβ -couplings are significantly averaged.
(10) Fig. 1 Ribbon pictures of the structure of HEWL with residues for which experimentally derived backbone 3 J HNHα -coupling values (left panel, experimentally stereo-specifically assigned (set bb1) in blue, computationally stereo-specifically assigned (set bb2) in magenta) and side-chain 3 J Hα-Hβ -coupling values (right panel, experimentally stereo-specifically assigned (set sc1) in green, computationally stereospecifically assigned (set sc2) in red) are available

Materials and methods
Energy minimisations and molecular dynamics simulations were performed using the GROMOS bio-molecular simulation software (Schmid et al. 2011a(Schmid et al. , 2012van Gunsteren et al. 2019).

Molecular model
The protein was modelled using the GROMOS bio-molecular force field 54A7 (Poger et al. 2010;Schmid et al. 2011b).
In view of the pH used in the experimental NMR measurements, pH = 3.5, only Glu 35 was protonated and His was doubly protonated (Bartik et al. 1994). The simple point charge (SPC) model (Berendsen et al. 1981) was used to describe the solvent molecules in the rectangular periodic box. To compensate for the overall positive charge of the protein, 10 Cl − ions were included in the solution. All bond lengths and the bond angle of the water molecules were kept rigid with a relative geometric precision of 10 -4 using the SHAKE algorithm (Ryckaert et al. 1977), allowing for a 2 fs MD time step in the leap-frog algorithm (Hockney and Eastwood 1981) used to integrate the equations of motion. For the non-bonded interactions a triple-range method (van Gunsteren et al. 1986) with cut-off radii of 0.8/1.4 nm was used. Short-range (within 0.8 nm) van der Waals and electrostatic interactions were evaluated every time step based on a charge-group pair list (van Gunsteren et al. 2019). Medium-range van der Waals and electrostatic interactions, between pairs at a distance larger than 0.8 nm and shorter than 1.4 nm, were evaluated every fifth time step (10 fs), at which time point the pair list was updated, and kept constant between updates. Outside the larger cut-off radius (1.4 nm) a reaction-field approximation (Barker and Watts 1973;Tironi et al. 1995) with a relative dielectric permittivity of 61 (Heinz et al. 2001) was used. Minimum-image periodic boundary conditions were applied.

Simulation set-up
Four X-ray crystal structures were used as initial structures for the energy minimisations followed by MD simulations.
1. Structure "2VB1" of the Protein Data Bank (PDB) (Berman et al. 2000), derived from a triclinic unit cell at 0.065 nm resolution at T = 100 K. It contains multiple side-chain conformations for 46 residues. 2. Structure "4LZT" of the PDB, derived from a triclinic unit cell at 0.095 nm resolution at T = 295 K. It contains multiple side-chain conformations for 8 residues.
3. Structure "1IEE" of the PDB, derived from a tetragonal unit cell at 0.094 nm resolution at T = 110 K. It contains multiple side-chain conformations for 33 residues. 4. Structure "1AKI" of the PDB, derived from an orthorhombic unit cell at 0.15 nm resolution at T = 298 K. It contains no multiple side-chain conformations.
An initial structure was first energy minimised in vacuo to release possible strain induced by small differences in bond lengths, bond angles, improper dihedral angles, and short distance non-bonded contacts between the force-field parameters and the X-ray structure. Subsequently, the protein was put into a rectangular box filled with water molecules, such that the minimum solute-wall distance was 1.0 nm, and water molecules closer than 0.23 nm from the solute were removed. This resulted in boxes with 12,157 water molecules for the initial protein structures. In order to relax unfavourable contacts between atoms of the solute and the solvent, a second energy minimisation was performed for the protein in the periodic box with water while keeping the atoms of the solute harmonically position-restrained (van Gunsteren et al. 2019) with a force constant of 25,000 kJ mol −1 nm −2 .
The resulting protein-water configuration was used as initial configuration for the MD simulation. In order to avoid artificial deformations in the protein structure due to relatively high-energy atomic interactions still present in the system, the MD simulation was started at T = 60 K and then the temperature was slowly raised to T = 308 K. Initial atomic velocities were sampled from a Maxwell distribution at T = 60 K. The equilibration scheme consisted of five short 20 ps simulations at temperatures 60, 120, 180, 240 and 308 K at constant volume. During the first four of the equilibration periods, the solute atoms were harmonically restrained to their positions in the initial structures with force constants of 25,000, 2500, 250, and 25 kJ mol −1 nm −2 . The temperature was kept constant using the weak-coupling algorithm (Berendsen et al. 1984) with a relaxation or coupling time τ Τ = 0.1 ps. Solute and solvent were separately coupled to the heat bath. Following this equilibration procedure, the simulations were performed at a reference temperature of 308 K and a reference pressure of 1 atm. The pressure was kept constant using the weak-coupling algorithm (Berendsen et al. 1984) with a coupling time τ p = 0.5 ps and an isothermal compressibility κ T = 4.575 × 10 -4 (kJ mol −1 nm −3 ) −1 . The centre of mass motion of the system was removed every 1000 time steps (2 ps).

J-coupling restraining
Two sets of backbone 3 J HN-Hα couplings and two sets of sidechain 3 J Hα-Hβ couplings for restraining  were used, see Tables 1, 2, 3 and 4. 1. A set (bb1) of 95 backbone 3 J HN-Hα -coupling values, see Table II of Smith et al. (1991) from which the values for 11 glycine residues were omitted, because these had not been stereo-specifically assigned. 2. A set (bb2) of 22 experimentally stereo-specifically unassigned backbone 3 J HN-Hα -coupling values for the 11 glycine residues, see Table II of Smith et al. (1991). 10 of these were stereo-specifically assigned based on the 3 J HN-Hα -coupling values calculated from the four unrestrained MD simulations starting from the four mentioned X-ray structures (see below) as either 4 or 3 of the unrestrained MD simulations suggested the same stereo-specific assignment. Gly 104 could not be stereo-specifically assigned using this criterion. Instead it was stereo-specifically assigned based on the 3 J HN-Hαcoupling values calculated for the four X-ray structures. 3. A set (sc1) of 58 3 J Hα-Hβ -coupling values, see Tables  III and IV of Smith et al. (1991), which were stereospecifically assigned using experimental data. 4. A set (sc2) of 38 out of 40 experimentally stereo-specifically unassigned 3 J Hα-Hβ -coupling values, see Table  III of Smith et al. (1991), which were stereo-specifically assigned based on the 3 J Hα-Hβ -coupling values calculated from the four unrestrained MD simulations starting from the four mentioned X-ray structures (see below) as either 4 or 3 of the unrestrained MD simulations suggested the same stereo-specific assignment. Glu 7 could not be stereo-specifically assigned using this criterion. It could also not be stereo-specifically assigned using the 3 J Hα-Hβcoupling values calculated for the four X-ray structures.
The experimentally derived 3 J HN-Hα -coupling values for Val 2, Thr 51, Asp 66, Cys 115, Thr 118 and Ile 124 lie outside the Karplus curve, so were set to 9.7 Hz, which is the maximum of the Karplus curve used (Pardi et al. 1984). None of the experimentally derived 3 J Hα-Hβ -coupling values lie outside the Karplus curve used (deMarco et al. 1978). The nomenclature for the H α2 and H α3 atoms in Gly residues and the H β , H β2 and H β3 atoms in the side chains was defined as in Markley et al. (1998), see Fig. 3.
The initial implementation of the algorithm proposed in Smith et al. (2016) in the GROMOS MD program was incorrect. The factors t −1 and dt in the expressions of Eqs. (6) and (7) (Eqs. (21) and (22) in Smith et al. 2016) were not implemented, which meant that the local-elevation weight factors ω φ or ω θ were not time-averaged and of different magnitude than intended. They would thus only stay constant instead of decrease when the restraints were satisfied. This implementation error has been corrected. As a consequence, the values of the restraining force constants K Jr used in the previous work  are much smaller than the value of the force constant K Jr used in the present work. In Smith et al. (2016), K Jr was varied between 5 and 50⋅10 -4 kJ mol −1 Hz −4 . Here, K Jr is set to 50 kJ mol −1 Hz −4 , because the weights ω θ are much smaller. In Smith et al. (2016), the memory relaxation time was varied between 5 and 50 ps, because of the length of 2 ns of the many test simulations. Here it could be extended to τ Q ≡ τ J = 500 ps, in view of the 20 ns length of the MD simulations. In Smith et al. (2016), the flat-bottom parameter of the restraining potential energy term was varied between 0.5 and 1 Hz. Here the value Δ 3 J fb = 1.0 Hz is used, which means a flat bottom of 2 Hz width , representing the uncertainty of the Karplus relation used.

MD simulations performed
Four unrestrained MD simulations, starting from the four mentioned X-ray crystal structures, were performed: each 20 ns long. The average solute temperatures were 311 K and the solvent temperatures 312 K.

Analysis of atomic trajectories and X-ray structures
Trajectory energies and atomic coordinates were stored at 5 ps intervals and used for analysis (Eichenberger et al. 2011).
In view of the various factors contributing to an uncertainty of about 2 Hz inherent to the Karplus relation linking structure and 3 J-couplings, as discussed in the Introduction, a deviation of less than 2 Hz between 3 J-coupling values calculated from X-ray or MD trajectory structures and 3 J-coupling values derived from experiment is considered insignificant.
The GROMOS force fields treat aliphatic carbons as united CH, CH 2 and CH 3 atoms. So inter-hydrogen distances involving the aliphatic hydrogen atoms were calculated using virtual atomic positions for CH and pro-chiral CH 2 (van Gunsteren et al. 1985) and pseudo-atomic positions for CH 3 (Wüthrich et al. 1983) for those hydrogen atoms (van Gunsteren et al. 2019). The pseudo-atom NOE distance bound corrections of Wüthrich et al. (1983) were used ). The set of NOE distance bounds (Smith et al. 1993;Schwalbe et al. 2001) can be found in Table S1 of Supporting Information, together with values of the five simulations starting from the 2VB1 X-ray crystal structure. The NOE between Trp 28 HZ3 and Leu 56 HG was reassigned as between Trp 28 HZ3 and Leu 56 HD* following reassessment of the experimental spectra. Interhydrogen distances were calculated as 〈r −3 〉 −1/3 , i.e. using r −3 averaging over the trajectory structures, where r indicates the actual hydrogen-hydrogen distance. In view of the uncertainty inherent to the calculation of NOE bounds and r −3 averaged distances, deviations from experiment of less than 0.1 nm are considered insignificant.
S 2 order parameters were calculated using the ensemble averaging expression (Henry and Szabo 1985) where τ indicates the time-averaging window, here 1 ns, shorter than the rotational correlation time of 5.7 ns of HEWL in solution (Smith et al. 1993), are the normalised x-, y-, and z-components of the vector r XY ≡ r X − r Y connecting atoms X and Y, and r XY ≡ |r XY | its length (Hansen et al. 2014). To obtain a dimensionless quantity the term in curly brackets is multiplied with the 6th power of the effective length (r eff XY ) of the vector r XY . Because in the present work bond length constraints are used, the length (11) of r XY is essentially constant over time and thus equal to its effective value r eff XY . The set of 79 experimentally derived S 2 values (Buck et al. 1995;Moorman et al. 2012) together with the S 2 values calculated from the five MD simulations started from the 2VB1 X-ray crystal structure (Smith et al. 2020) can be found in Supporting Information, Table S2.
In view of the uncertainty inherent to the derivation of S 2 XY (exp)-values from relaxation experiments and inherent to the calculation of S 2 XY (MD)-values from MD simulation, a deviation of less than 0.2 between simulation and experiment is considered insignificant.
Atom-positional root-mean-square fluctuations, i.e. around their average positions, in the MD trajectories were calculated using the above-mentioned superposition too.
The secondary structure assignment was done with the program DSSP, based on the Kabsch-Sander rules (Kabsch and Sander 1983).
Hydrogen bonds were identified according to a geometric criterion: a hydrogen bond was assumed to exist if the hydrogen-acceptor distance was smaller than 0.25 nm and the donor-hydrogen-acceptor angle was larger than 135°. Table 1 lists 95 backbone 3 J HNHα -coupling values derived and stereo-specifically assigned based on NMR measurements (set bb1), as calculated from the four unrestrained MD simulations starting from four X-ray crystal structures, and as calculated from the four X-ray structures. The mean of the four MD or X-ray values and the root-mean-square deviation (RMSD) from it are also presented. Deviations from the experimental values of more than 2 Hz are denoted in italics. The mean values of the MD simulations show 14 deviations larger than 2 Hz, while for the X-ray structures there is only one such case (Trp 111). The great majority of the large differences between MD simulation and experiment are found for residues (11), for which large 3 J HNHα -coupling values (> 9 Hz) have been observed (Ser 36, Gln 41, Asn 65, Arg 68, Thr 69, Leu 84, Trp 108, Arg 114, Cys 115, Thr 118 and Ile 124). This can be explained by considering the Karplus curve in Fig. 2 (left panel). Its maximum of 9.7 Hz lies at φ ≈ 240° or − 120° and 3 J HNHα -coupling values larger than 8 Hz are only found for φ-angle values between 210° and 270°. Any motion around φ ≈ 240° will lower the calculated 3 J HNHα -coupling value. We note that other parametrisations of the Karplus curve show maxima of 10 Hz (Wang and Bax 1996) or even 11 Hz (Brüschweiler and Case 1994), see Fig. 1 in Dolenc et al. (2010). The remaining three residues with 3 J HNHα -coupling value differences with experiment larger than 2 Hz are found for Ser 50, Ile 78 and Asn 103, whose experimental values are 7.8, 8.0 and 8.2 Hz, respectively. For the X-ray structures, the only 3 J HNHα -coupling value differing more than 2 Hz from experiment (7.1 Hz) is for Trp 111, where all four values are smaller than 4.9 Hz. Use of the Karplus relation of Brüschweiler and Case (1994) (Fig. 2, left panel, blue lines) slightly lowers the 3 J HNHα -coupling values smaller than 6.3 Hz and increases values larger than 6.3 Hz, the increase being about 1 Hz for 3 J HNHα -coupling values larger than 8.5 Hz. It does not improve significantly the agreement of the 3 J HNHα -coupling values obtained from the MD simulations and the X-ray structures with the experimental ones. Table 2 lists 22 experimentally stereo-specifically unassigned backbone 3 J HNHα -coupling values derived from NMR measurements (set bb2), as calculated from the four unrestrained MD simulations starting from four X-ray crystal structures, and as calculated from the four X-ray structures. The mean of the four MD or X-ray values and the root-meansquare deviation (RMSD) from it are also presented. The stereo-specific assignment of the experimental values for the α2 Re and the α3 Si hydrogens in glycine residues is based on the criterion that the 3 J HN-Hα -coupling values calculated from the four unrestrained MD simulations starting from the four X-ray structures do suggest in 4 or 3 of the unrestrained MD simulations the same assignment. Only Gly 104 could not be stereo-specifically assigned using this criterion. For this residue the stereo-specific assignment was based on the 3 J HN-Hα -coupling values calculated from the four X-ray structures using a corresponding criterion. For only two residues the X-ray structures would suggest stereo-specific assignments different from the ones based on the MD trajectories, for Gly 67, which shows a rather small difference of 0.7 Hz between the two experimental 3 J HN-Hα -coupling values, and for Gly 126, with a somewhat larger difference of 1.2 Hz between the two experimental 3 J HN-Hα -coupling values. Deviations of the MD or X-ray calculated 3 J HN-Hα -coupling values from the experimental ones of more than 2 Hz are denoted in italics. The MD trajectories show four of such deviations, for Gly 49 α2 Re in the MD_2VB1 simulation, for Gly 104 α2 Re and α3 Si in the MD_4LZT simulation and α3 Si in the MD_1AKI simulation. The X-ray structures show seven of such deviations, for Gly 102 α2 Re and α3 Si in the 1IEE and 1AKI structures, for Gly 104 α3 Si in the 2VB1 structure and for Gly 126 α3 Si in the 4LZT and 1AKI structures. Table 3 lists 58 side-chain 3 J HαHβ -coupling values derived and stereo-specifically assigned based on NMR measurements (set sc1), as calculated from the four unrestrained MD simulations starting from four X-ray crystal structures, and as calculated from the four X-ray structures. The mean of the four MD or X-ray values and the root-mean-square deviation (RMSD) from it are also presented. Deviations from the experimental values larger than 2 Hz are denoted in italics. Out of 58*4 = 232 3 J HαHβ -coupling values, the MD trajectories yield 58 values 2 Hz larger than experiment, and the X-ray structures 55 values. The MD trajectories show a variation in 3 J HαHβ -coupling values larger than 2 Hz for seven hydrogens in five residues, Tyr 23 β 2 , Asn 27 β 2 , Phe 34 β 2 and β 3 , Asn 46 β 2 and β 3 , and Cys 127 β 2 . The X-ray structures show two such cases, for Val 99 and Val 109. In all these cases the mean of the four MD or X-ray values lies within 2 Hz from experiment, but the four values contributing to the mean show a large variation of 2.2-3.7 Hz for the MD values and 3.4-4.6 for the X-ray values. Apparently, MD trajectories and different X-ray structures contain different side-chain χ 1 -angle conformers for these residues. Table 4 lists 40 experimentally stereo-specifically unassigned side-chain 3 J HαHβ -coupling values derived from NMR measurements (set sc2), as calculated from the four unrestrained MD simulations starting from four X-ray crystal structures, and as calculated from the four X-ray structures. The mean of the four MD or X-ray values and the rootmean-square deviation (RMSD) from it are also presented. The stereo-specific assignment of the experimental values for the β 2 and β 3 hydrogens is based on the criterion that the 3 J HαHβ -coupling values calculated from the four unrestrained MD simulations starting from the four X-ray structures do suggest in 4 or 3 of the unrestrained MD simulations the same stereo-specific assignment. Only Glu 7 could not be stereo-specifically assigned using this criterion. For four residues the X-ray structures would suggest stereo-specific assignments different from the ones based on the MD trajectories, for Asn 19, which shows a rather small difference of 0.9 Hz between the two experimental 3 J HαHβ -coupling values, for Asn 37, with a larger difference of 3.9 Hz between the two experimental 3 J HαHβ -coupling values, for Asn 74, with a difference of 6.6 Hz between the two 3 J HαHβ -coupling values, and for Asn 77, with a difference of 3.4. Hz between the two 3 J HαHβ -coupling values. Deviations of the MD or X-ray calculated 3 J HαHβ -coupling values from the experimental ones of more than 2 Hz are denoted in italics.   Smith et al. (1991). The values within brackets in the six "MD simulation (X-ray crystal)" columns represent the 3 J HNHα -couplings calculated from the four X-ray crystal structures, their mean and RMSD values. The stereo-specific assignment of the experimental values for the α2 Re and the α3 Si hydrogens in glycine residues is based on the criterion that the 3 J HN-Hα -coupling values calculated from the four unrestrained MD simulations starting from the four X-ray structures do suggest in 4 or 3 of the unrestrained MD simulations the same stereo-specific assignment. Only Gly 104 could not be stereo-specifically assigned using this criterion. For this residue the stereo-specific assignment was based on the 3 J HN-Hα -coupling values calculated from the four X-ray structures using a corresponding criterion. MD or X-ray values differing more than 2 Hz from the experimental value and the maximum of the Karplus curve are denoted using italics

Comparison of 3 J-coupling values calculated from 3 J-coupling time-averaging local-elevation restraining MD trajectories with NMR measured values
Val 2 10.0 (9.7) 7.9 (1.6) 9.0 (0.7) 6.7 (1.6) 7.5 (1.7) 9.0 (0.7) Phe 3 7.4 7.5 (1.7) 8.0 (1.7) 8.9 (1.2) 6.0 (1.9) 7.8 (1.5) Cys 6 5.     Experimental values from Table II of Smith et al. (1991). The value within brackets in the column "Experimental value" represents the maximum in the Karplus relation (Pardi et al. 1984) used for the calculation of the 3 J HNHα -couplings. The root-mean-square fluctuations (RMSF) of the 3 J-couplings in the MD simulations are given within parentheses. MD values differing more than 2 Hz from the experimental value and the maximum of the Karplus curve are denoted using italics experiment for the 58 side-chain 3 J HαHβ -couplings. No deviations larger than 2 Hz are observed. This is also the case when restraining towards all 213 experimental backbone and side-chain 3 J-coupling values (sets bb1, bb2, sc1 and sc2). Table 8 lists 40 experimentally stereo-specifically unassigned side-chain 3 J HαHβ -coupling values derived from NMR measurements (set sc2 plus Glu 7) and values from the unrestrained and 3 J-coupling time-averaging local-elevation restrained MD simulations starting from the 2VB1 X-ray crystal structure using different sets of backbone and side-chain restraints. The unrestrained MD simulation shows 13 3 J HαHβ -coupling values (in italics) that deviate more than 2 Hz from the experimental values (residues 68, 72(2), 77(2), 85, 86, 101(2), 106,  (22) in Hz derived from NMR measurements (set bb2) and values from the unrestrained and 3 J-coupling time-averaging local-elevation restrained MD simulations starting from the 2VB1 X-ray crystal structure and using different sets of backbone and side-chain restraints Experimental values from Table II of Smith et al. (1991). The root-mean-square fluctuations (RMSF) of the 3 J-couplings in the MD simulations are given within parentheses. Stereo-specific assignments for restraining were based on the 3 J HN-Hα -coupling values calculated from the four unrestrained MD simulations starting from the four X-ray structures in case 4 or 3 of the unrestrained MD simulations suggested the same stereo-specific assignment. Only Gly 104 could not be stereo-specifically assigned using this criterion. For this residue the stereo-specific assignment was based on the 3 J HN-Hα -coupling values calculated from the four X-ray structures using a corresponding criterion. MD values differing more than 2 Hz from the experimental value and the maximum of the Karplus curve are denoted using italics   (58) in Hz derived and stereo-specifically assigned based on NMR measurements (set sc1) and from the unrestrained and 3 J-coupling time-averaging local-elevation restrained MD simulations starting from the 2VB1 X-ray crystal structure and using different sets of backbone and side-chain restraints

125
(2), 128). Restraining towards the sets bb1 and bb2 of 95 and 22 backbone 3 J HNHα -coupling values yields 13 deviations larger than 2 Hz, and so no improvement of the agreement between simulated and experimental 3 J HαHβ -coupling values, as one would expect. Backbone 3 J HNHα -coupling restraining does not improve the agreement between simulation and experiment for the sidechain 3 J HαHβ -couplings, apart from a few cases. This is also observed using the sc1 set of 58 restraints (14 deviations larger than 2 Hz). 3 J HαHβ -coupling time-averaging local-elevation restraining towards the sets sc1 and sc2 of 58 and 38 target side-chain 3 J HαHβ -coupling values leads, as expected, to good agreement between simulation and experiment for the 40 side-chain 3 J HαHβ -couplings. Only two deviations (Glu 7 β 2 and β 3 ) larger than 2 Hz are observed. This is also the case when restraining towards all 213 experimental backbone and side-chain 3 J-coupling values (sets bb1, bb2, sc1 and sc2), only for Glu 7 β 2 the deviation is larger than 2 Hz. Note that Glu 7 is not part of the set sc2 of restraints. Tables 9, 10, 11 and 12 summarise the agreement between experimental 3 J-coupling values and those calculated from the four X-ray structures, from the four unrestrained MD simulation trajectories and from the 3 J-coupling time-averaging local-elevation restrained MD simulation trajectories. Using backbone 3 J HNHα -coupling restraints (sets bb1 and bb2) only, the experimental 3 J HNHα -coupling values are reproduced and using 3 J HαHβ -coupling restraints (sc1 or sc1 and sc2) only, the experimental 3 J HαHβ -coupling values are reproduced. Not surprisingly, there appears to be insignificant mutual influence between the different sets of restraints. Using all 213 experimental backbone and sidechain 3 J-coupling restraints (sets bb1, bb2, sc1 and sc2) no significant deviations from experimental 3 J-coupling values are observed.

Other quantities and agreement with NOE atomatom distance bounds and S 2 order-parameter values
The application of backbone 3 J HNHα -coupling time-averaging local-elevation restraining in MD simulation does not significantly influence the secondary structure of the protein, with the four main α-helices, the two 3 10 -helices and the triple-stranded anti-parallel β-sheet in the protein all being maintained (see Fig. S1 and S2 of Supporting Information). However, there are some subtle differences in loop regions, particularly around residues Gly 102-Gly 104 and around Gly 117-Thr 118. For example, there are increases in the populations of the hydrogen bonds 107 NH-104 O and 118 NH-115 O and a decrease in the population of the hydrogen bond 104 NH-101 O compared to the unrestrained simulations. These regions are known to be mobile in solution, Gly 102, Asn 103 and Thr 118 all having lower backbone NH order parameters (0.72, 0.52 and 0.72, respectively; Buck et al. 1995) and these regions were less well defined in the NMR structure of HEWL (Schwalbe et al. 2001) having few longer range NOE identified for them. They provide an example of the extra conformational insights that could be obtained with time-averaging local-elevation restraining to glycine 3 J HNHα -couplings although further experimental data would be needed to confirm the details of the hydrogen bond population changes observed here. As expected, the time-averaging local-elevation restraining does slightly enhance the atomic mobility, as observed from the backbone atom-positional fluctuations, see Figure S3 in Supporting Information.

Importance of time-averaging
Values of measured 3 J-couplings may the result from considerable conformational averaging. This is illustrated in Figs. 4, 5, 6, 7 and 8 for one backbone φ-angle and four side-chain χ 1 -angles. Figure 4 shows for Asn 103 the variation of the backbone φ-angle determining the 3 J HNHα -coupling, of this 3 J HNHαcoupling with an experimental value of 8.2 Hz (green line), and the presence of the Asn 103 NH-Asp 101 OD1 hydrogen bond as function of time for the unrestrained (left panels) and the 3 J-coupling time-averaging local-elevation restraining to all 213 experimental 3 J-coupling values (right panels) MD simulations starting from the 2VB1 X-ray structure. A variation of the φ-angle of about 30° around its average value of − 69° leads to the 3 J HNHα -coupling covering the range 2-10 Hz (RMSF 1.4 Hz) and an average value of 5.4 Hz (red line) in the unrestrained MD simulation, and in the 3 J-coupling time-averaging local-elevation restraining simulation to a distribution of 3 J HNHα -couplings shifted to larger values (RMSF 1.5 Hz) with an average of 8.1 Hz. The restraining shifts the average 3 J HNHα -coupling more than 2 Hz towards the experimental value and slightly increases its fluctuation. The Asn 103 NH-Asp 101 OD1 hydrogen bond regularly populated in the unrestrained simulation Table 11 Number of deviations, | 3 J HαHβ (exp) − 3 J HαHβ (MD or X-ray)|, for the 58 side-chain 3 J HαHβ -coupling values derived and stereo-specifically assigned based on NMR measurements (set sc1), in four X-ray crystal structures, in the four unrestrained MD simulations starting from these, and in the 3 J-coupling time-averaging localelevation restrained MD simulations starting from the 2VB1 X-ray crystal structure and using different sets of backbone and side-chain restraints Crystal structure or simulation Size of 3 J HαHβ deviation (in Hz) 1-2 2-3 3-4 4-5 > 5  Table 12 Number of deviations, | 3 J HαHβ (exp) − 3 J HαHβ (MD or X-ray)|, for the 38 side-chain 3 J HαHβ -coupling values derived but stereo-specifically unassigned from NMR measurements (set sc2), in four X-ray crystal structures, in the four unrestrained MD simulations starting from these, and in the 3 J-coupling time-averaging localelevation restrained MD simulations starting from the 2VB1 X-ray crystal structure and using different sets of backbone and side-chain restraints Stereo-specific assignments for restraining were based on the 3 J Hα-Hβcoupling values calculated from the four unrestrained MD simulations starting from the four X-ray structures in case 4 or 3 of the unrestrained MD simulations suggested the same stereo-specific assignment. Only Glu 7 could not be stereo-specifically assigned using this criterion. Its 3 J Hα-Hβ -coupling values were not included in the calculation Crystal structure or simulation Size of 3 J HαHβ deviation (in Hz) 1-2 2-3 3-4 4-5 > 5    (Smith et al. 2020), in the four unrestrained MD simulations starting from four different X-ray crystal structures and in the 3 J-coupling time-averaging local-elevation restrained MD simulations starting from the 2VB1 X-ray crystal structure and using different sets of backbone and side-chain restraints. Order parameter values are reported in Table S2

Fig. 5
Variation of the sidechain χ 1 -angle (degree) determining the 3 J Hα-Hβ -coupling, of this 3 J Hα-Hβ -coupling (Hz) with an experimental value of 9.5 Hz (green line), and the presence of the Thr 89 OG1-HG1-Asp 87 OD1 hydrogen bond for residue Thr 89 as function of time from the unrestrained MD simulation MD_2VB1 (left panels) and from the 3 J-coupling time-averaging local-elevation restraining (to all 213 experimental 3 J-coupling values) MD simulation MD_VB1_bb1 + bb2 + sc1 + sc2 (right panels), each starting from the 2VB1 X-ray structure. Red lines: average 3 J-coupling value in the MD simulations 87 OD1 hydrogen bond is somewhat reduced by the restraining. The 2VB1 X-ray χ 1 -angle values of − 67° or − 69° result in 3 J Hα-Hβ -couplings of 12.8 Hz and 12.7 Hz respectively, deviations of more than 3 Hz from the experimental value of 9.5 Hz. Figure 6 shows for Val 99 the variation of the sidechain χ 1 -angle determining the 3 J Hα-Hβ -coupling and of this 3 J Hα-Hβ -coupling with an experimental value of 6.3 Hz (green line), as function of time from the unrestrained (left panels) and 3 J-coupling time-averaging local-elevation restraining to all 213 experimental 3 J-coupling values (right panels) MD simulations starting from the 2VB1 X-ray structure. The unrestrained simulation shows a stable χ 1 -angle value of − 62° with little variation (RMSF 10°) resulting in an average 3 J Hα-Hβ -coupling of 3.0 Hz (red line) (RMSF 1.6 Hz). It does not reproduce the experimental value of 6.3 Hz. The four X-ray structures deviate with 3 J Hα-Hβ -coupling values of 12.8 Hz (2VB1 and 4LZT), 2.2 Hz (1IEE) and 2.4 Hz (1AKI) Fig. 6 Variation of the side-chain χ 1 -angle (degree) determining the 3 J Hα-Hβ -coupling, and of this 3 J Hα-Hβ -coupling (Hz) with an experimental value of 6.3 Hz (green line), for residue Val 99 as function of time from the unrestrained MD simulation MD_VB1 (left panels) and from the 3 J-coupling time-averaging local-elevation restrain-ing (to all 213 experimental 3 J-coupling values) MD simulation MD_VB1_bb1 + bb2 + sc1 + sc2 (right panels), each starting from the 2VB1 X-ray structure. Red lines: average 3 J-coupling value in the MD simulations Fig. 7 Variation of the side-chain χ 1 -angle (degree) determining the two, β 2 and β 3 , 3 J Hα-Hβ -couplings, of these two 3 J Hα-Hβ -couplings (Hz) with experimental values of 5.6 Hz and 6.6 Hz respectively (green lines), and the presence of the Asn 103 NH-Asp 101 OD1 hydrogen bond for residue Asp 101 as function of time from the unrestrained MD simulation MD_2VB1 (left panels) and from the 3 J-coupling time-averaging local-elevation restraining (to all 213 experimental 3 J-coupling values) MD simulation MD_2VB1_bb1 + bb2 + sc1 + sc2 (right panels), each starting from the 2VB1 X-ray structure. Red lines: average 3 J-coupling value in the MD simulations even more from the experimental value. Using 3 J-coupling local-elevation restraining other conformations of the χ 1angle are accessed, resulting in a larger variation (RMSF 2.1 Hz) and raising the average to 5.3 Hz (red line). The 2VB1 X-ray structure contains two conformations, with χ 1angles of 176° (main conformation) and − 53° (alternative conformation) and 3 J Hα-Hβ -couplings of 12.8 Hz and 4.2 Hz respectively. Averaging over different conformations seems to occur in aqueous solution. Figure 7 shows for Asp 101 the variation of the side-chain χ 1 -angle determining the two, β 2 and β 3 , 3 J Hα-Hβ -couplings, of these two 3 J Hα-Hβ -couplings (Hz) with experimental values of 5.6 Hz and 6.6 Hz respectively (green lines), and the presence of the Asn 103 NH-Asp 101 OD1 hydrogen bond as function of time from the unrestrained (left panels) and 3 J-coupling local-elevation restraining to all 213 experimental 3 J-coupling values (right panels) MD simulations starting from the 2VB1 X-ray structure. The unrestrained simulation shows a stable χ 1 -angle value of − 169° with little variation (RMSF 11°) resulting in average 3 J Hα-Hβ -couplings of 2.5 Hz (RMSF 0.8 Hz) and 12.2 Hz (RMSF 0.9 Hz) (red lines). These values deviate significantly from the experimental values of 5.6 Hz and 6.6 Hz. The four X-ray structures deviate with χ 1 -angles of − 89° and 3 J Hα-Hβ -coupling values of 10.5 Hz (β 2 ) and 1.8 Hz (β 3 ) (2VB1 and 4LZT), a χ 1 -angle of − 167° and 3 J Hα-Hβ -couplings of 2.2 Hz (β 2 ) and 12.4 Hz (β 3 ) (1IEE) and a χ 1 -angle of − 139° and 3 J Hα-Hβ -couplings of 2.4 Hz (β 2 ) and 8.4 Hz (β 3 ) (1AKI) even more from the experimental values. Using 3 J-coupling time-averaging local-elevation restraining other conformations of the χ 1angle are accessed, resulting in a larger variations (RMSF 2.8 Hz (β 2 ) and 2.2 Hz (β 3 )) of the 3 J Hα-Hβ -couplings and raising the β 2 average to 5.6 Hz and lowering the β 3 average to 5.9 Hz (red line). Figure 8 shows for Asn 106 the variation of the side-chain χ 1 -angle determining the two, β 2 and β 3 , 3 J Hα-Hβ -couplings, of these two 3 J Hα-Hβ -couplings with experimental values of 10.5 Hz and 3.6 Hz respectively (green lines), and the presence of the Ala 107 NH-Asn 106 OD1 hydrogen bond as function of time from the unrestrained (left panels) and 3 J-coupling local-elevation restraining to all 213 experimental 3 J-coupling values (right panels) MD simulations starting from the 2VB1 X-ray structure. The unrestrained simulation shows after 4 ns a stable χ 1 -angle value of about 60° with little variation (RMSF 9°) resulting in average 3 J Hα-Hβcouplings of 4.7 Hz (RMSF 2.8 Hz) and 3.9 Hz (RMSF 2.7 Hz) (red lines). The β 2 value deviates significantly from the experimental value of 10.5 Hz, while the β 3 value is close to the experimental value of 3.6 Hz. The 2VB1 X-ray structure contains two χ 1 -angle conformations, − 96° (main conformation) and − 169° (alternative conformation), with 3 J Hα-Hβ -coupling values of 9.3 Hz (β 2 ) and 2.1 Hz (β 3 ). The other three X-ray structures show only one conformation, with a χ 1 -angle of − 70° and 3 J Hα-Hβ -couplings of 12.6 Hz (β 2 ) and 2.3 Hz (β 3 ) (4LZT), a χ 1 -angle of − 64° and 3 J Hα-Hβcouplings of 12.9 Hz (β 2 ) and 3.0 Hz (β 3 ) (1IEE) and a χ 1angle of − 70° and 3 J Hα-Hβ -couplings of 12.6 Hz (β 2 ) and 2.4 Hz (β 3 ) (1AKI). Using 3 J-coupling time-averaging Fig. 8 Variation of the side-chain χ 1 -angle (degree) determining the two, β 2 and β 3 , 3 J Hα-Hβ -couplings, of these two 3 J Hα-Hβ -couplings (Hz) with experimental values of 10.5 Hz and 3.6 Hz respectively (green lines), and the presence of the Ala 107 NH-Asn 106 OD1 hydrogen bond for residue Asn 106 as function of time from the unrestrained local-elevation restraining other conformations of the χ 1angle are accessed, resulting in smaller variations (RMSF 2.3 Hz (β 2 ) and 1.9 Hz (β 3 )) of the 3 J Hα-Hβ -couplings, raising the β 2 average to 11.1 Hz and lowering the β 3 average to 3.1 Hz (red lines). The Ala 107 NH-Asn 106 OD1 hydrogen bond disappears when applying 3 J-coupling restraining.

Importance of escaping from torsional-angle energy minima
In Figs. 9, 10, 11, 12 and 13 the local-elevation potential energies after the time-averaging local-elevation restraining (to all 213 experimental 3 J-coupling values) MD simulation MD_2VB1_bb1 + bb2 + sc1 + sc2 for the φ-angle of Asn 103 and the χ 1 -angles of Thr 89, Val 99, Asp 101 and Asn 106, are shown. The time series of these angles, the corresponding 3 J-couplings and some hydrogen bonds in this simulation were shown in Figs. 4, 5, 6, 7 and 8 (right panels).
The  (Fig. 5,   Fig. 9 Local-elevation 3 J HNHα -coupling restraining potential energy as function of the backbone φ-angle for residue Asn 103, built-up during the 3 J-coupling time-averaging local-elevation restraining (to all 213 experimental 3 J-coupling values) MD simulation MD_2VB1_ bb1 + bb2 + sc1 + sc2 starting from the 2VB1 X-ray structure Fig. 10 Local-elevation 3 J HαHβ -coupling restraining potential energy as function of the side-chain χ 1 -angle for residue Thr 89, built-up during the 3 J-coupling time-averaging local-elevation restraining (to all 213 experimental 3 J-coupling values) MD simulation MD_2VB1_ bb1 + bb2 + sc1 + sc2 starting from the 2VB1 X-ray structure Fig. 11 Local-elevation 3 J HαHβ -coupling restraining potential energy as function of the side-chain χ 1 -angle for residue Val 99, built-up during the 3 J-coupling time-averaging local-elevation restraining (to all 213 experimental 3 J-coupling values) MD simulation MD_2VB1_ bb1 + bb2 + sc1 + sc2 starting from the 2VB1 X-ray structure upper right panel), shifting the 3 J Hα-Hβ -coupling from 4.8 Hz in the unrestrained simulation to 10.4 Hz, near the experimental value of 9.5 Hz. Figure 10 shows the local-elevation potential-energy term built-up in the 3 J-coupling restraining simulation for χ 1 -angle values in the range [− 10°, + 80°] and around − 55°, that causes this shift. The build-up of local-elevation potential energy around − 170° is caused by the χ 1 -angle occasionally visiting this region in the 3 J-coupling restraining simulation (Fig. 5, upper right panel). We note that the solvent accessibility of the side chain of Thr 89 in the 2VB1 X-ray structure is 76%.
The side-chain χ 1 -angle of Val 99 in the unrestrained simulation initially covers values around 70° and 170° (Fig. 6, upper left panel) yielding 3 J Hα-Hβ -couplings between 2 and 13 Hz. After about 1 ns, the χ 1 -angle stabilizes around − 62° with values in the range [− 50°, − 85°] yielding a 3 J Hα-Hβcoupling of 3.0 Hz, much lower than the experimental value of 6.3 Hz. In the Karplus curve for this 3 J-coupling (dashed line in the right panel of Fig. 2) this value corresponds to four χ 1 -angle values of about − 129°, − 37°, + 37° and + 129°. The Karplus curve indicates that a slight shift of the χ 1 -angle distribution towards less negative values, for example − 40°, would yield a 3 J Hα-Hβ -coupling of about 6 Hz. Time-averaging local-elevation 3 J-coupling restraining indeed induces this slight shift towards χ 1 -angle values in the range [− 40°, − 65°] (Fig. 6, upper right panel), but also makes the χ 1 -angle repeatedly move over all angle values but those between − 40° and + 40°, thereby reaching very large 3 J Hα-Hβ -coupling values, in order to raise the average 3 J Hα-Hβ -coupling value from 3.0 Hz in the unrestrained simulation towards the experimental value of 6.3 Hz. Figure 11 shows the local-elevation potential-energy term built-up in the 3 J-coupling restraining simulation for χ 1 -angle values around − 80°, + 80° and 180°. It shows minima for − 126°, + 130° and for [− 35°, + 35°]. The χ 1 -angle range [− 35°, + 35°] is disfavoured by the χ 1 dihedral-angle potential-energy term and the non-bonded van der Waals interaction of the force field. We note that the solvent accessibility of the side chain of Val 99 in the 2VB1 X-ray structure is only 7%. This side chain is surrounded by the side chains of Tyr 20, Trp 28, Ile 98 and Tyr 108.
The last two examples involve longer side chains for which the C β -atom is connected to two hydrogens, H β2 and H β3 . If their 3 J Hα-Hβ -couplings have been stereo-specifically assigned, the corresponding χ 1 -angle can be restrained using two local-elevation 3 J Hα-Hβ -coupling restraining potentialenergy terms. This imposes more restriction on the χ 1 -angle motion than in the previously discussed cases.
The side-chain χ 1 -angle of Asp 101 covers in the unrestrained simulation values around − 170° (Fig. 7, upper left panel) yielding average 3 J Hα-Hβ -couplings of 2.5 Hz for β 2 and 12.2 Hz for β 3 , to be compared to experimental values of 5.6 Hz and 6.6 Hz respectively. In the Karplus curve for the 3 J Hα-Hβ2 -coupling (solid line in the right panel of Fig. 2) the experimental value of 5.6 Hz corresponds to four χ 1angle values of about − 116°, − 4°, + 76° and + 164°. In the Karplus curve for the 3 J Hα-Hβ3 -coupling (dashed line in the right panel of Fig. 2) the experimental value of 6.6 Hz corresponds to four χ 1 -angle values of about − 129°, − 37°, + 37° and + 129°. The Karplus curves suggest a χ 1 -angle distribution dominantly around − 125°. χ 1 -angle values around + 55° would bring both 3 J Hα-Hβ -couplings near their experimental values, but would require averaging of 3 J Hα-Hβ -couplings. Fig. 12 Local-elevation 3 J HαHβ2 -coupling (dashed line) and 3 J HαHβ3coupling (dotted line) restraining potential energies and their sum (solid line) as function of the side-chain χ 1 -angle for residue Asp 101, built-up during the 3 J-coupling time-averaging local-elevation restraining (to all 213 experimental 3 J-coupling values) MD simulation MD_2VB1_bb1 + bb2 + sc1 + sc2 starting from the 2VB1 X-ray structure Fig. 13 Local-elevation 3 J HαHβ2 -coupling (dashed line) and 3 J HαHβ3coupling (dotted line) restraining potential energies and their sum (solid line) as function of the side-chain χ 1 -angle for residue Asn 106, built-up during the 3 J-coupling time-averaging local-elevation restraining (to all 213 experimental 3 J-coupling values) MD simulation MD_2VB1_bb1 + bb2 + sc1 + sc2 starting from the 2VB1 X-ray structure Time-averaging local-elevation 3 J-coupling restraining indeed shifts the χ 1 -angle distribution towards values around − 125° (Fig. 7, upper right panel), but also makes the χ 1angle repeatedly move over all angle values. Figure 12 shows the β 2 and β 3 local-elevation potential-energy terms (dashed and dotted lines, respectively) and their sum (solid line) builtup in the 3 J-coupling restraining simulation. Major buildup is observed for χ 1 -angle values around − 165°, − 70° and + 65°. The first two keep the χ 1 -angle around − 125°. Yet, excursions to other χ 1 -angle values seem required to push the average 3 J Hα-Hβ2 -and 3 J Hα-Hβ3 -couplings towards the corresponding experimental values of 5.6 Hz and 6.6 Hz, respectively. We note that the solvent accessibility of the side chain of Asp 101 in the 2VB1 X-ray structure is 42%.
The side-chain χ 1 -angle of Asn 106 in the unrestrained simulation initially covers values around − 70° and − 170° (Fig. 8, upper left panel) yielding 3 J Hα-Hβ2 -and 3 J Hα-Hβ3couplings both fluctuating between 2 and 13 Hz. After about 4 ns, the χ 1 -angle stabilizes around + 60° yielding an average 3 J Hα-Hβ2 -coupling of 4.7 Hz and an average 3 J Hα-Hβ3 -coupling of 3.9 Hz, the former much lower than the experimental value of 10.5 Hz, the latter close to the experimental value of 3.6 Hz. In the Karplus curve for the 3 J Hα-Hβ2 -coupling (solid line in the right panel of Fig. 2) the experimental value of 10.5 Hz corresponds to two χ 1 -angle values of about − 89° and − 31°. In the Karplus curve for the 3 J Hα-Hβ3 -coupling (dashed line in the right panel of Fig. 2) the experimental value of 3.6 Hz corresponds to four χ 1 -angle values of about − 111°, − 58°, + 58° and + 111°. The Karplus curves suggest a χ 1 -angle distribution dominantly around − 75° with much averaging in the range [− 100°, − 45°]. Time-averaging local-elevation 3 J-coupling restraining indeed shifts the χ 1 -angle distribution towards values around between − 100° and − 45° (Fig. 8, upper right panel). Figure 13 shows the β 2 and β 3 local-elevation potentialenergy terms (dashed and dotted lines, respectively) and their sum (solid line) built-up in the 3 J-coupling restraining simulation. Major build-up is observed for χ 1 -angle values between 0° and 90° and around − 170°. The small build-up at − 65° is due to the β 2 restraint, keeping the 3 J Hα-Hβ2 -coupling below 12 Hz. The averaging between the two local-elevation energy minima in the range [− 100°, − 25°] results in average 3 J Hα-Hβ2 -and 3 J Hα-Hβ3 -couplings of 11.1 Hz and 3.1 Hz, close to the experimental values of 10.5 Hz and 3.6 Hz, respectively. We note that the solvent accessibility of the side chain of Asn 106 in the 2VB1 X-ray structure is 76%. Figures 4,5,6,7 and 8 show examples of dihedral angles for which the 3 J-couplings in the unrestrained simulations do not agree with the experimental values. The local-elevation potential-energy term used in the restraining simulations changes the dihedral-angle distribution such that the average 3 J-coupling matches experiment. Thus such a potentialenergy term contains information on a possible modification of the dihedral-angle term of the force field used. If a consistent picture of possible modifications would emerge from simulations of a collection of proteins, such a modification could be incorporated into the force field.

Conclusions
Although 3 J-coupling constants are relatively easy to obtain from NMR experiments, their use in structure determination of proteins has been rather limited due to different aspects of the measurement and the relation between a 3 J-coupling and molecular structure . The Karplus relation between structure and 3 J-coupling value is multi-valued, with up to four torsional angle values mapping to a single 3 J-coupling value. In addition, intermediate 3 J-coupling values (4-8 Hz) are sensitive to the experimental averaging period, which is rather long. This is in contrast to some other NMR measurable quantities such as NOE intensities. The difficulty of accounting for conformational averaging and for the multi-valued function of a torsional angle in terms of a 3 J-coupling has severely hampered the use of 3 J-couplings in protein structure determination or refinement. However, with the advent of time-averaging local-elevation restraining MD simulation (Christen et al. 2007;Smith et al. 2016) both problems could be solved. Time-averaging can be taken into account in a simulation (Torda et al. 1989) and a molecular conformation can be induced to escape from an (incorrect) local minima due to restraining based on the Karplus relation by use of the local-elevation algorithm (Huber et al. 1994). This makes a comprehensive use of 3 J-coupling data in structure determination possible .
The application of 3 J-coupling time-averaging localelevation restraining to the protein HEWL shows that this technique is able to produce a conformational ensemble compatible with the experimental 3 J-coupling data. Analysis of the conformations underlying the 3 J-couplings shows that conformational averaging plays an essential role in a number of cases and that finding an alternative minimum energy conformation for backbone φ or side-chain χ 1 angles is also of importance.
The 3 J-coupling time-averaging local-elevation restraining does improve the agreement with 1630 NOE atom-atom distance bounds for HEWL, but only if all 213 backbone and side-chain 3 J-coupling restraints are applied. It has no significant effect upon the agreement of the conformational ensemble with values of 121 backbone S 2 NH and 79 sidechain S 2 CH and S 2 NH order parameters for HEWL. This one would more or less expect, considering the different degrees of freedom involved.
The results for the backbone 3 J HNHα -couplings based on the parametrisation of the Karplus relation due to Pardi et al. (1984) show that this parametrisation is not capable of reproducing large 3 J HNHα -coupling values, due to a maximum of 9.7 Hz of this Karplus relation. One could use alternative parametrisations, such as the one of Wang and Bax (1996) with a maximum of 10 Hz or the one of Brüschweiler and Case (1994) with a maximum of 11 Hz. However, in the present case of HEWL, no significant change of the agreement with experiment is observed by using the latter Karplus relation in the analysis.
Of particular interest from the work reported here, is the description of the behaviour of the side-chains with a high level of conformational mobility. Examples such as the data for the side-chains of Asp 101 and Asn 106 shown in Figs. 7 and 8 indicate that any hydrogen bonds involving flexible side-chains, which are present in crystal structures of the protein, may be very fluctuating or completely absent in solution. Indeed, the simulation results show just how conformationally disordered the surface of the protein really is. This needs to be kept in mind when X-ray structures are being used in areas such as drug design or to help with the interpretation of data from receptor binding or mutational studies.