Conformational Investigations in Flexible Molecules Using Orientational NMR Constraints in Combination with 3J-Couplings and NOE Distances

The downscaling of NMR tensorial interactions, such as dipolar couplings, from tens of kilohertz to a few hertz in low-order media is the result of dynamics spanning several orders of magnitudes, including vibrational modes (~ns-fs), whole-molecule reorientation (~ns) and higher barrier internal conformational exchange (<ms). In this work, we propose to employ these dynamically averaged interactions to drive an “alignment-tensor-free” molecular dynamic simulation with orientation constraints (MDOC) in order to efficiently access the conformational space sampled by flexible small molecules such as natural products. Key to this approach is the application of tensorial pseudo-force restraints which simultaneously guide the overall reorientation and conformational fluctuations based on defined memory function over the running trajectory. With the molecular mechanics force-field, which includes bond polarization theory (BPT), and complemented with other available NMR parameters such as NOEs and scalar J-couplings, MDOC efficiently arrives at dynamic ensembles that reproduce the entire NMR dataset with exquisite accuracy and theoretically reveal the systems conformational space and equilibrium. The method as well as its potential towards configurational elucidation is presented on diastereomeric pairs of flexible molecules: a small 1,4-diketone 1 with a single rotatable bond as well as a 24-ring macrolide related to the natural product mandelalide A 2.


Introduction
Tensorial NMR observables, such as chemical shift anisotropy, dipole-dipole couplings or nuclear quadrupolar splittings, give access to structurally-relevant local orientation at nuclear sites relative to the external magnetic field. This information can also be retrieved in solution at high-resolution by intentionally imposing an anisotropy to the otherwise Brownian molecular reorientation. Since their first description in solution [1], such interactions have been used to understand structure and dynamic properties of small molecules dispersed in liquid-crystalline solvents [2]. In the last two decades, a number of media have been proposed that scale these interactions down to small perturbations on the order of a few Hz which are as easy to evaluate as scalar couplings. These structure-dependent "residual" dipolar couplings (RDCs) or "residual" chemical shift anisotropy (RCSAs), when taken together, report on the geometry of the molecule as a whole, complementing the more commonly-used, exclusively short-range observables, such as chemical shifts (CS), nuclear Overhauser enhancements (NOEs) and scalar (or J-) couplings [3,4]. Whereas, for biological systems, these tensorial interactions have been

NMR Datasets
For the first assessment of the MDOC approach, a simple flexible chiral molecule was required. For this task, the 1,4-diketone 1 was synthesized as a coupling product between 2-methylnaphth-1-yl hydrazine and 1-indanone in an asymmetric Brønsted acid catalyzed dearomatizing redox cross coupling reaction [34] (Figure 1). Characteristic for this molecule are the two chiral and relatively rigid subgroups (a naphthalenone group and an indenone group) connected by a single rotatable bond. The modest stereoselectivity (dr < 2.5) during the trials of this reaction resulted in both diastereomers of the 1,4-diketone 1 being available for this study. The labels 1-SR and 1-SS are applied henceforth to represent the two diastereomers, differing only in the chirality at carbon C12, for simplicity. The diastereomers were separated by chromatography and provided two samples for dynamic investigation by NMR. structural ensemble and show encouraging potential towards guiding the correct configuration assignment.

NMR Datasets
For the first assessment of the MDOC approach, a simple flexible chiral molecule was required. For this task, the 1,4-diketone 1 was synthesized as a coupling product between 2-methylnaphth-1yl hydrazine and 1-indanone in an asymmetric Brønsted acid catalyzed dearomatizing redox cross coupling reaction [34] (Figure 1). Characteristic for this molecule are the two chiral and relatively rigid subgroups (a naphthalenone group and an indenone group) connected by a single rotatable bond. The modest stereoselectivity (dr < 2.5) during the trials of this reaction resulted in both diastereomers of the 1,4-diketone 1 being available for this study. The labels 1-SR and 1-SS are applied henceforth to represent the two diastereomers, differing only in the chirality at carbon C12, for simplicity. The diastereomers were separated by chromatography and provided two samples for dynamic investigation by NMR. The theoretical rotamers (trans, gauche (-) and gauche (+)) of 1 were first geometry-optimized and their relative energies calculated using DFT (see supporting information). Table S1 shows that the trans and the gauche (-) geometries are the most stable among the conformers of 1-SS and 1-SR, respectively, which both correspond to the configuration where the small indenone H12 is simultaneously gauche to the bulkier methyl (C11) and carbonyl (C1) groups. Nevertheless, the other The theoretical rotamers (trans, gauche (-) and gauche (+)) of 1 were first geometry-optimized and their relative energies calculated using DFT (see Supporting Information). Table S1 shows that the trans and the gauche (-) geometries are the most stable among the conformers of 1-SS and 1-SR, respectively, which both correspond to the configuration where the small indenone H12 is simultaneously gauche to the bulkier methyl (C11) and carbonyl (C1) groups. Nevertheless, the other rotamers should also be accessible, especially for 1-SR with all three relative energies within 1 kcal/mol. For 1-SS, the theoretical differences are somewhat larger (E tr − E g+ = 3.6 kcal/mol) due to the combination of the repulsive effect of the ketones and overall steric effect of the other 2 rotamers. It is expected however that these differences could be mitigated by the polar solvent. These optimized structures provide already theoretical evidence for significant rotameric exchange and will serve also as starting structures for the dynamical evaluations.
The NMR investigation provided two datasets (1A (1-SS) and 1B for (1-SR)) comprising carefully evaluated NOESY cross-peaks integrals and converted to interproton distances and a set of 3 J-coupling constants related to torsion angle distribution (for datasets, see Supplementary Information). Both the measured NOE contacts as well as the homo-and heteronuclear 3 J couplings in the initial NMR analysis indicate some level of rotameric exchange, fast on the NMR timescale. For instance, there are no distinctly large 3 J CH couplings (>9 ppm) from H12 to C1, C9 or C11 across the rotatable C10-C12 bond that would clearly hint to a static trans relation. Furthermore, the semi-quantitative evaluation of the interannular NOEs from H13 or H12 to H11 or H9 show that all pairs correlate with NOEs of similar magnitudes which also imply some form of dynamic averaging.
The datasets were extended with HH and CH RDC values (22 in the case of 1A and 24 for 1B), obtained from a set of coupled experiments in the isotropic (CDCl 3 ) and the anisotropic (PMMA/CDCl 3 ) sample for each diastereomer. An RDC-fitting procedure using MSPIN (SVD) against the corresponding three DFT rotamers identified the 1A and 1B dataset to fit best to 1-SS-tr (Q = 0.19) and 1-SR-g(-) (Q = 0.21) conformation, respectively (see Supporting Information). These rotamers match with those with the lowest calculated DFT energies. However, the quality of the best fits were less than optimal as indicated by the quality criterion (n/χ 2 ) and outlier criterion (1/χ 2 min ) far below unity (1A: n/χ 2 = 0.07, 1/χ 2 min = 0.013 and 1B: n/χ 2 = 0.05, 1/χ 2 min = 0.006), as well as only F = 6/22 (27%) and 5/24 (21%) of the computed RDCs being valid, i.e., having values within the estimated error margins. These observations indicate that overall, a majority of RDCs cannot be accurately rationalized by a lone rotamer.
The attempts to fit RDCs to a multi-rotameric distribution with a common averaged alignment tensor lead to a substantial improvement in the fitting (1A: Q = 0.16 and 1B: Q = 0.14), however the above mentioned criteria are still less than satisfactory (1A: n/χ 2 = 0.11, 1/χ 2 min = 0.01 and F = 41% (9/22) and 1B: n/χ 2 = 0.08, 1/χ 2 min = 0.015 and F = 54%(13/24)). It should be mentioned that these only modest improvements may be in part due to the assumption of an existing average alignment tensor; indeed, since the overall shape of the molecule changes substantially from one rotamer to the other, a strong dependence of the alignment to the rotameric state is to be expected. Unfortunately, the sparseness of the RDC dataset did not allow the use of multiple alignment tensors since this procedure would have required at a minimum 17 independent RDC values (3 × 5 for the Saupe (alignment) matrix + 2 for the three populations). A further extension of the dataset (e.g., with less accessible long-range HC RDC couplings) would not have improved the required low condition number since a large portion of the RDC vector orientations are confined within or close to the aromatic planes and are therefore not completely independent. The inclusion of the NOE interproton information and the 3 J-coupling (CH) in the stereofitter module of MSPIN did not give significantly better results (data not shown).
The example of the diketone 1 demonstrates, as for many molecules previously studied in past, that the linear fitting approaches like SVD constitute a powerful tool to roughly discriminate structures especially between a finite number of possible conformers; however, as a molecule exhibits increasingly more disorder, the method becomes progressively more deficient towards accurately explaining dynamic exchange. The aim of the MDOC approach introduced in the next section therefore presents a novel tool to obtain an ensemble of structures based on NMR observables that describe the conformational space for each diastereomer.

MDOC Simulation
Two MDOC simulations were performed combining the structure configurations 1-SR and the 1-SS with their corresponding NMR constraints from datasets 1A and 1B (for the details of the simulations see Tables S10 and S11 in Supplementary Information). Datasets 1A and 1B included 1-bond CH RDCs (1A:14, 1B:14) and long-range HH RDCs (1A:8, 1B:10)), NOEs (1A: 7, 1B:5), 3J-couplings (1A: 5, 1B:5), each translated into a pseudo-force term for the MDOC simulations. The outcome was independent of the initial conformation, since the pseudo-forces calculated from the RDC have the ability to rapidly generate all possible rotamers as well as reorient the molecule as a whole, on a timescale related to the memory time τ of the mean value calculation (in this case 200 ps, Equation (3)).
In Figure 2, typical trajectories of 13 C-1 H dipolar splittings are displayed. The left panel shows the fluctuation of the dipolar coupling of the aromatic C-H group at position C6 calculated using the exponential function with a memory time τ of 200 ps, according to Equation (3). After a short simulation time of about 1 ns-approximately 5 times the exponential rise constant of the pseudo-forces-the time mean value of the dipolar coupling fluctuates mostly within a range of about ±1 Hz. The experimental error of 1.2 Hz is indicated in Figure 2 with green lines. The pseudo-force width constant ∆D (Equation (5)) was set to 0.5 Hz which is close to the error ranges of most experimental values. Figure 2 also shows the convergence of the running mean value (red) of the trajectory towards a value close to the experimental RDC and reaching a constant value after about 10 ns, indicating that this MDOC simulation period was sufficient to average the dipolar coupling through combined whole-molecule reorientation and internal dynamics. Two MDOC simulations were performed combining the structure configurations 1-SR and the 1-SS with their corresponding NMR constraints from datasets 1A and 1B (for the details of the simulations see Tables S10 and S11 in supplementary information). Datasets 1A and 1B included 1bond CH RDCs (1A:14, 1B:14) and long-range HH RDCs (1A:8, 1B:10)), NOEs (1A: 7, 1B:5), 3Jcouplings (1A: 5, 1B:5), each translated into a pseudo-force term for the MDOC simulations. The outcome was independent of the initial conformation, since the pseudo-forces calculated from the RDC have the ability to rapidly generate all possible rotamers as well as reorient the molecule as a whole, on a timescale related to the memory time τ of the mean value calculation (in this case 200 ps, Equation (3)).
In Figure 2, typical trajectories of 13 C-1 H dipolar splittings are displayed. The left panel shows the fluctuation of the dipolar coupling of the aromatic C-H group at position C6 calculated using the exponential function with a memory time τ of 200 ps, according to Equation (3). After a short simulation time of about 1 ns -approximately 5 times the exponential rise constant of the pseudoforces -the time mean value of the dipolar coupling fluctuates mostly within a range of about ±1 Hz. The experimental error of 1.2 Hz is indicated in Figure 2 with green lines. The pseudo-force width constant ΔD (Equation (5)) was set to 0.5 Hz which is close to the error ranges of most experimental values. Figure 2 also shows the convergence of the running mean value (red) of the trajectory towards a value close to the experimental RDC and reaching a constant value after about 10 ns, indicating that this MDOC simulation period was sufficient to average the dipolar coupling through combined whole-molecule reorientation and internal dynamics. In statistical analysis of traditional MD simulations, the initial time period before the system reaches a thermal equilibrium is often disregarded. The state of equilibrium in MDOC is reached when the mean values with exponential memory fluctuate only within the experimental error bounds, mostly after about 1 ns. Flowingly, the first nanosecond in the final analysis is dropped according to the decay function described in Equation (3).
The right panel of Figure 2 displays the interesting case of the CH3 group which involves one additional mode of motion -the axial reorientation of the CH3 group. Though the local order parameters for freely rotating methyl groups are typically much smaller than in the case of the aromatic CH groups, the conditions for a successful convergence are also fulfilled within these time limits. Figure 3 shows that the running ensembles of the MDOC simulation for both 1-SS and 1-SR are in good accord with the entire corresponding measured datasets, which include four categories of In statistical analysis of traditional MD simulations, the initial time period before the system reaches a thermal equilibrium is often disregarded. The state of equilibrium in MDOC is reached when the mean values with exponential memory fluctuate only within the experimental error bounds, mostly after about 1 ns. Flowingly, the first nanosecond in the final analysis is dropped according to the decay function described in Equation (3).
The right panel of Figure 2 displays the interesting case of the CH 3 group which involves one additional mode of motion-the axial reorientation of the CH 3 group. Though the local order parameters for freely rotating methyl groups are typically much smaller than in the case of the aromatic CH groups, the conditions for a successful convergence are also fulfilled within these time limits. Figure 3 shows that the running ensembles of the MDOC simulation for both 1-SS and 1-SR are in good accord with the entire corresponding measured datasets, which include four categories of NMR parameter types (one-bond and long-range RDCs, 3 J-couplings and NOEs). In stark contrast to the outcome of the linear analysis based only on RDC values fit to static rotameric DFT models, the aforementioned quality criterion (n/χ 2 ) exceeds unity for all but one parameter type, the only exception being the "long-range" ( n D, n > 1) RDCs for the 1-SS configuration with n/χ 2 = 0.76 ( Figure 3). This exception is specifically affected by the calculated H6-H8 coupling, with a deviation of 1.4 Hz from the measured RDC, whose error was estimated to be ± 0.9 Hz. In the case of the MDOC simulation of the 1-SS and 1-SR configurations with dataset 1A and 1B, there are 6 (F = 34/40 (85%)) and 4 (F = 37/41 (90%)) outliers observed, with χ −2 min values of 0.29 and 0.13 respectively (see supplementary information). Though these numbers represent a large improvement over the corresponding values in the static SVD study (black in Figure 3: 1A: #outliers = 13/22 and χ −2 min = 0.015; 1B: #outliers = 11/22 and χ −2 min = 0.010), it will be important to determine the factors contributing most to the remaining incongruity.
The analysis of the trajectories from the MDOC simulations also provided insight into the population of conformers. In Figure 4, the torsion angle ω distributions for the central C10-C12 bond are presented. Clearly, the three maxima represent the accessible rotameric states of 1 corresponding to trans, gauche(-) and gauche(+) subensembles. In the 1-SS form the trans-conformation is clearly favored whereas in the 1-SR form the gauche(-) conformation is the dominant conformation. This is in general accordance with the relative DFT energies from the geometry optimized DFT models and with the results of the SVD analysis (see supporting information). More specifically, however, a statistical analysis showed that the relative population of the individual rotameric states ratios   (7)), "smallest outlier" criterion 1/χ 2 min and fidelity F (labels) for MDOC simulations of the 1,4-diketone 1-SS and 1-SR with the experimental datasets 1A and 1B, respectively. On average, all interaction types are very well reproduced by ensembles, with n/χ 2 > 1 (overall 1A: 1.9, 1B: 1.8) and fidelity F above 85% (overall 1A: 34/40, 1B: 37/40).
The analysis of the trajectories from the MDOC simulations also provided insight into the population of conformers. In Figure 4, the torsion angle ω distributions for the central C10-C12 bond are presented. Clearly, the three maxima represent the accessible rotameric states of 1 corresponding to trans, gauche(-) and gauche(+) subensembles. In the 1-SS form the trans-conformation is clearly favored whereas in the 1-SR form the gauche(-) conformation is the dominant conformation. This is in general accordance with the relative DFT energies from the geometry optimized DFT models and with the results of the SVD analysis (see Supporting Information). More specifically, however, a statistical analysis showed that the relative population of the individual rotameric states ratios   Right inset diagram: torsion angle distribution of (C8-C9-C10-C1 in Figure 1). Left inset diagram: torsion angle distribution of (C14-C12-C13-C20 in Figure 1).
The distribution of the C8-C9-C10-C1 ( ) torsion angle shows maximum puckered forms at −40° and +40° although a continuous distribution of more planar forms are possible as well (see inset panel, Figure 5). More distinct maxima at ±35° are seen for the C14-C12-C13-C20 ( ) torsion angle within the 5-membered ring indicating a stronger preference of twisted forms. By comparison, the molecular model of 1-SR in gauche(-) conformation obtained from a DFT geometry optimization, a static torsion angle of = −6.4° and = 13.5° are observed, which are not dominant in the MDOC simulation. This Besides the rotation about the C10-C12 bond, the MDOC simulations reveal the presence of conformational isomeric exchange within the non-aromatic rings. Indeed, both rings of the molecules are not fully planar and therefore more than one twist forms of the rings are possible ( Figure 5).  f exchange within the non-aromatic rings. Indeed, both rings of the molecules are not fully planar and therefore more than one twist forms of the rings are possible ( Figure 5). Right inset diagram: torsion angle distribution of (C8-C9-C10-C1 in Figure 1). Left inset diagram: torsion angle distribution of (C14-C12-C13-C20 in Figure 1).
The distribution of the C8-C9-C10-C1 ( ) torsion angle shows maximum puckered forms at −40° and +40° although a continuous distribution of more planar forms are possible as well (see inset panel, Figure 5). More distinct maxima at ±35° are seen for the C14-C12-C13-C20 ( ) torsion angle within the 5-membered ring indicating a stronger preference of twisted forms. By comparison, the molecular model of 1-SR in gauche(-) conformation obtained from a DFT geometry optimization, a static torsion angle of = −6.4° and = 13.5° are observed, which are not dominant in the MDOC simulation. This  Figure 1). Left inset diagram: torsion angle distribution of τ 2 (C14-C12-C13-C20 in Figure 1).
The distribution of the C8-C9-C10-C1 (τ 1 ) torsion angle shows maximum puckered forms at −40 • and +40 • although a continuous distribution of more planar forms are possible as well (see inset panel, Figure 5). More distinct maxima at ±35 • are seen for the C14-C12-C13-C20 (τ 2 ) torsion angle within the 5-membered ring indicating a stronger preference of twisted forms. By comparison, the molecular model of 1-SR in gauche(-) conformation obtained from a DFT geometry optimization, a static torsion angle of τ 1 = −6.4 • and τ 2 = 13.5 • are observed, which are not dominant in the MDOC simulation. This raises the question whether this ring-twist dynamic behavior represented by this ensemble distribution is real or rather an artifact of the MDOC method. For the τ 2 dihedral angle, this question could be investigated, by considering the 3 J HH couplings of H13a and H13b to H11. Using the Altona equation [4], the 1-SR-g(-) DFT model gives a value of 3 J H13b-H12 = 8.5 Hz. The MDOC simulation which was run at a mean temperature of 313 K gave rise to a mean value of 5.8 Hz within the estimated error range (±0.75 Hz) of the experimental value of 5.1 Hz, so that one can suppose that this isomeric exchange occurs also in reality.

NMR Datasets
Mandelalide A is a glycosylated macrolide with interesting cytotoxic properties found in a species of Lissoclinum ascidian (marine invertebrate) from South Africa [35]. This natural product has been the object of two recent total synthetic projects [36,37]. An original configuration of the large lactone ring containing 9 stereocenters was proposed based on scalar coupling constants and interproton NOEs (2p in Figure 6). The tedious syntheses ultimately led to the correction of the configuration, reallocating the entire northern part to the inverted configuration (2r in Figure 6). raises the question whether this ring-twist dynamic behavior represented by this ensemble distribution is real or rather an artifact of the MDOC method. For the dihedral angle, this question could be investigated, by considering the 3 JHH couplings of H13a and H13b to H11. Using the Altona equation [4], the 1-SR-g(-) DFT model gives a value of 3 JH13b-H12 = 8.5 Hz. The MDOC simulation which was run at a mean temperature of 313 K gave rise to a mean value of 5.8 Hz within the estimated error range (±0.75 Hz) of the experimental value of 5.1 Hz, so that one can suppose that this isomeric exchange occurs also in reality.

NMR Datasets
Mandelalide A is a glycosylated macrolide with interesting cytotoxic properties found in a species of Lissoclinum ascidian (marine invertebrate) from South Africa [35]. This natural product has been the object of two recent total synthetic projects [36,37]. An original configuration of the large lactone ring containing 9 stereocenters was proposed based on scalar coupling constants and interproton NOEs (2p in Figure 6). The tedious syntheses ultimately led to the correction of the configuration, reallocating the entire northern part to the inverted configuration (2r in Figure 6).
Since the configurations 2p and its 11-epi-2p were synthesized in larger scale in a local laboratory and characterized using NMR, this excellent material was made available to test the MDOC methodology on a challenging, larger-sized and flexible molecule. The current structural characterization of mandelalide A isomers relies on carefully evaluated RDC values, NOE distances and scalar 3 J couplings. In the case of structure 2p (11-epi-2p), a full unambiguous 1 H and 13 C assignment was obtained (including the prochiral 1 H assignment 7 CH2), and complemented with 45 (48) CH RDC values, 38 (38) 3 JHH couplings and 129 (106) NOE distances (see supplementary information). Any attempt to obtain alignment tensor evaluation using the SVD approach on calculated structural models failed (data not shown). The NMR-based structural information, however, was used as constraints in a 41-ns MDOC simulations (see supplementary information). Since the configurations 2p and its 11-epi-2p were synthesized in larger scale in a local laboratory and characterized using NMR, this excellent material was made available to test the MDOC methodology on a challenging, larger-sized and flexible molecule.
The current structural characterization of mandelalide A isomers relies on carefully evaluated RDC values, NOE distances and scalar 3 J couplings. In the case of structure 2p (11-epi-2p), a full unambiguous 1 H and 13 C assignment was obtained (including the prochiral 1 H assignment 7 CH 2 ), and complemented with 45 (48) CH RDC values, 38 (38) 3 J HH couplings and 129 (106) NOE distances (see Supplementary  Information). Any attempt to obtain alignment tensor evaluation using the SVD approach on calculated structural models failed (data not shown). The NMR-based structural information, however, was used as constraints in a 41-ns MDOC simulations (see Supplementary Information).
The MDOC simulations (see SI for details) of the madelalide A isomers are more challenging than the 1,4-diketone simulations because conformational changes of the ring system have larger energy barriers. The rotational degrees of freedom of the side chains add up to the complexity of the system. As indicated by the n/χ 2 quality criteria (Figure 7), values back-calculated from the MDOC trajectory are on average well within the experimental error bounds (n/χ 2 > 1). As for the MDOC simulation of the 1,4-diketone 1, inspection of the individual data reveals however the occasional parameter were not completely fulfilled (see labels representing the fidelity F in Figure 7). The outliers represent generally less than 10% of the overall NMR constraints. Furthermore, based on the magnitude of the outlier criterion, χ −2 min , the worst of these outliers deviate by less than one error margin. Whereas the error margins for the RDC values were determined individually, those for the NOE distances were estimated to 0.5Å. In the case of 3 J couplings, the errors were assessed to be about 1.0 Hz taking into account also the possible uncertainties of the equation of Haasnoot et al. [4]. The RMS deviation of the NOE distances was lower than 0.3 Å and that of the 3 J couplings below 0.6 Hz. The MDOC simulations (see SI for details) of the madelalide A isomers are more challenging than the 1,4-diketone simulations because conformational changes of the ring system have larger energy barriers. The rotational degrees of freedom of the side chains add up to the complexity of the system. As indicated by the n/χ 2 quality criteria (Figure 7), values back-calculated from the MDOC trajectory are on average well within the experimental error bounds (n/χ 2 > 1). As for the MDOC simulation of the 1,4-diketone 1, inspection of the individual data reveals however the occasional parameter were not completely fulfilled (see labels representing the fidelity ℱ in Figure 7). The outliers represent generally less than 10% of the overall NMR constraints. Furthermore, based on the magnitude of the outlier criterion, , the worst of these outliers deviate by less than one error margin. Whereas the error margins for the RDC values were determined individually, those for the NOE distances were estimated to 0.5Å. In the case of 3 J couplings, the errors were assessed to be about 1.0 Hz taking into account also the possible uncertainties of the equation of Haasnoot et al. [4]. The RMS deviation of the NOE distances was lower than 0.3 Å and that of the 3 J couplings below 0.6 Hz. The final mean data values of the MDOC simulations are calculated from 2000 coordinate, skipping the first nanosecond of the 41-ns trajectories. These snapshots contain the calculated NMR data (RDC, 3 J couplings and NOE distances) as averaged according to the memory function as given in Equation (3). In other words, the snapshots contain average information from every time step of the trajectory. Considering the n/χ 2 criteria (Equation (7)) or the calculated data (see supplementary information), it can be stated that the MDOC results are mostly within the estimated error bounds and no severe outliers are observed. Since the NMR data are both time-and ensemble-averaged In other words, the snapshots contain average information from every time step of the trajectory. Considering the n/χ 2 criteria (Equation (7)) or the calculated data (see Supplementary Information), it can be stated that the MDOC results are mostly within the estimated error bounds and no severe outliers are observed. Since the NMR data are both time-and ensemble-averaged values, it can also be stated that the MDOC simulation for a single molecule behaves to a good approximation in an ergodic manner. In this sense, the 2000 coordinate and data snapshots represent the final result of the simulations. Since this larger flexible molecules may undergo many conformational changes, the population of these states may be of high interest and therefore advanced methods need to be employed to analyze these populations. Figure 8 presents a representative conformer of the configuration 2p with torsion angles commonly populated within the 1000 MDOC coordinate snapshots (the Supplementary Information contains a  collection of 10 typical conformers). This conformer also exhibits the lowest total force field energy among all snapshots. Also displayed Figure 8 are the torsion distributions of some σ-bonds that display highest variability throughout the trajectory. As can be seen, many bonds sample a wide range of dihedral angles even within a closed 24-membered macrolide ring, indicating large amplitude motions and complex conformation exchange. Noteworthily, the oxane 6-ring (-C5-C6-C7-C8-C9-O-) did not change from its chair conformation throughout the simulation although the oxolane 5-ring displayed two twist states as demonstrated by σ 19 with a slight preference of g(-) conformation (Figure 9, top-right inset). values, it can also be stated that the MDOC simulation for a single molecule behaves to a good approximation in an ergodic manner. In this sense, the 2000 coordinate and data snapshots represent the final result of the simulations. Since this larger flexible molecules may undergo many conformational changes, the population of these states may be of high interest and therefore advanced methods need to be employed to analyze these populations. Figure 8 presents a representative conformer of the configuration 2p with torsion angles commonly populated within the 1000 MDOC coordinate snapshots (the supplementary information contains a collection of 10 typical conformers). This conformer also exhibits the lowest total force field energy among all snapshots. Also displayed Figure 8 are the torsion distributions of some σ-bonds that display highest variability throughout the trajectory. As can be seen, many bonds sample a wide range of dihedral angles even within a closed 24-membered macrolide ring, indicating large amplitude motions and complex conformation exchange. Noteworthily, the oxane 6-ring (-C5-C6-C7-C8-C9-O-) did not change from its chair conformation throughout the simulation although the oxolane 5-ring displayed two twist states as demonstrated by σ19 with a slight preference of g(-) conformation ( Figure 9, top-right inset). The representation of the fluctuation of the tethered rhamnose ring proved to be a special case. From its sparse and weak NOEs to the lactone ring, it can be expected to undergo conformational  (Figure 8). The dominant rhamnose ring conformation was characterized by the O-CH3 in axial position of the ring and the two OH groups and the CH3 group in equatorial positions. This conformation is strongly supported by the 3 J couplings of the protons of the rhamnose ring and rendered by the MDOC simulation, although nearly 10% rhamnose ring inversions were also generated. In Figure 9, a typical conformer of 11-epi-2p selected from a set of 1000 snapshots of an 40-ns MDOC simulation is presented, along with the distribution of key dihedral angles. It is instructive to compare the significantly different rotameric distribution of torsion angles with the values given for configuration 2p in Figure 8. Not surprisingly, the dynamics present in the smaller rings are almost identical in both diastereomers (e.g.,

Torsion Angle Distributions
). In the larger macrolide ring, however, not only are the torsion angle distributions in the direct vicinity of the epimerization site affected ( and ) but also those at moderately and very distant sites (e.g., , , and . This is indicative of long-range trans-annular steric interactions of the methyl group at position C11. However, a clearer conformational picture can only emerge from the examination of the entire structural ensemble generated by the MDOC simulation. The representation of the fluctuation of the tethered rhamnose ring proved to be a special case. From its sparse and weak NOEs to the lactone ring, it can be expected to undergo conformational changes with low correlation to the large macrolide ring. Indeed, the MDOC simulation resulted in a flexible oxygen bridge with several large scale modes of motions in order to average the parameters down to reliable mean values in the rhamnose system. The ether bond to the oxane ring (C7-O) showed nearly no preference for a torsion state while the ether bond to the rhamnose ring (O-C51) favored the trans-position to the carbon with the O-CH 3 group (Figure 8). The dominant rhamnose ring conformation was characterized by the O-CH 3 in axial position of the ring and the two OH groups and the CH 3 group in equatorial positions. This conformation is strongly supported by the 3 J couplings of the protons of the rhamnose ring and rendered by the MDOC simulation, although nearly 10% rhamnose ring inversions were also generated.

Principal Component Analysis
In Figure 9, a typical conformer of 11-epi-2p selected from a set of 1000 snapshots of an 40-ns MDOC simulation is presented, along with the distribution of key dihedral angles. It is instructive to compare the significantly different rotameric distribution of torsion angles with the values given for configuration 2p in Figure 8. Not surprisingly, the dynamics present in the smaller rings are almost identical in both diastereomers (e.g., σ 19 ). In the larger macrolide ring, however, not only are the torsion angle distributions in the direct vicinity of the epimerization site affected (σ 11 and σ 12 ) but also those at moderately and very distant sites (e.g., σ 0 , σ 4 , σ 5 and σ 17 ). This is indicative of long-range trans-annular steric interactions of the methyl group at position C11. However, a clearer conformational picture can only emerge from the examination of the entire structural ensemble generated by the MDOC simulation.

Principal Component Analysis
The representative models shown in Figure 8 were selected because their dihedral angles match a majority of the most highly populated rotameric states of the MDOC ensemble and they represent states of low force field energy. But just as for 1, it would be desirable to rationalize the dynamics of the mandelalide A isomers 2 in terms of population of major conformers with representative torsion angle combinations. Instead of only three rotameric states as in 1, there are nine "rotatable" σ-bond ( Figure 6) each with roughly three main rotamers, leading to nearly 20,000 possible combinations of angles. Though the MDOC simulations indicate that very broad distribution of angles is sampled by all of these nine key bonds, all combinations of angles are, of course, not accessible in reality. The challenge here lies in the complexity of the conformational space, and in finding the right tools to describe its essence.
A general and elegant method to evaluate the present MD data is the principal component analysis (PCA) on the basis of dihedral angles as developed by Altis et al. [38]. In its substance, PCA distills and ranks the geometrical parameters (e.g., atom coordinates) contributing most to a fluctuating system (e.g., MD trajectory) by identifying the orthogonal eigenvectors from the diagonalised co-variance matrix. In contrast to using methods based on Cartesian coordinates, PCA based on dihedral angles ("dPCA") provides a correct separation of internal and overall dynamics and places more emphasis on large structural rearrangements and less on bond length fluctuations. Since it deals with angular variables, dPCA has the added intricacy of accounting for circular statistics by transforming the angles to a complex coordinate space. A slightly adapted variant of the complex dPCA metod was used here to identify conformers with low energy.
In the case of the 2p isomer of mandelalide A, all 24 of the macrolide dihedral angles within the large ring were selected as the basis for the dPCA. Results are shown in Figure 10, where the relative contribution of each of the 24 dihedral angles to the angular variance is displayed for the first two principal components. Further components (n pc > 2) were already less significant. Out of the 24 bonds, five contribute very significantly to the overall conformational fluctuations within the MDOC ensemble, namely σ 0 , σ 1 , σ 4 , σ 11 and σ 12 , which according to Figure 8 are involved in fluctuations of large amplitude.The product of each eigenvector n with the dihedral angle matrix deliver a single angular value θ n for a given MD snapshot. In Figure 11, the (θ 1 , θ 2 ) distribution for the first two principal components obtained for the dPCA of mandelalide A 2p is presented from all 8000 data points each corresponding to an MDOC snapshot taken every 5 ps. The resulting distribution identifies three sub-regions at approximately (θ 1 , θ 2 ) = (+100, +130) [45%], (+100, −150) [40%] and (170, −150) [~5%]. In principle, these describe the subsets of large amplitude fluctuations contributing the most to the overall dynamics within the large ring system. Though the θ 1 and θ 2 angles bear no direct geometrical representation in the macrolide, they correspond to a linear combination of the angles highlighted in Figure 10. The major conformational exchanges between the three sub-regions thus involve correlated variations of these angles. Structural representation of typical conformers for the two principal subregions are also shown in Figure 11.
The question arises whether sparsely populated conformers are accurately represented in this simulation ensemble. The requirement that the ensemble meet the stringent rule of satisfying a very large number of NMR parameters within their error margins would certainly elicit emergence of low-populated states. So far, conformation subsets with populations ≥5% have proved to be essential for achieving the n/χ 2 criterion. Further investigations will however be required to establish the relationship between the abundance of short-and long-range constraints and the representation of low-energy conformations. There are several experimental and theoretical investigations that discuss precisely measured NMR parameters, especially NOE distances, as sensitive indicators on sparsely populated conformers to elucidate conformer equilibria [39,40].       . Dihedral landscape of the first two principal components of the mandelalide A isomer 2p with two typical molecules according their contribution to the principal components ( Figure 10). The left conformation is selected according to the lower θ 2 maximum of the dihedral distribution and the right molecule according to the upper θ 2 maximum. The side chains of the large ring system are displayed transparent. The arrows indicate the largest contributions to the first two principal components: both bonds to methyl-group-bearing C11 (σ 11 and σ 12 ), the bond to the ring oxygen O1 (σ 0 ) and the bond the methylene C4 (σ 4 ).
An analogous analysis was also conducted for the isomer 11-epi-2p (Supplementary Information). The corresponding dominantly contributing torsion angles are different to those of 2p: whereas σ 1 , σ 4 , σ 11 and σ 12 have similarly important co-variation, σ 0 and σ 17 are relatively more involved. Since the eigenvectors have different angular linear contributions, the dihedral landscapes for (θ 1 , θ 2 ) cannot be directly compared; still, one can notice that the 11-epi-2p conformational equilibrium is constituted of two major sub-regions ((θ 1 , θ 2 ) = (120 • ,160 • ) and (175 • ,120 • )) with well-defined distributions. The results here show that transannular steric influence reaching the opposite part of the 24-membered ring seems to favor an overall more open ring conformation for 11-epi-2p.
Intuitively, one might expect only minor differences in the conformational distribution of 2p and its isomer 11-epi-2p, whose configurations only differ in the chirality of a single methyl group, and these differences would be expected to be manifest mainly in the NMR data of nuclei in the direct vicinity of the C11-CH 3 group, but based on Figures 8 and 9 as well as in the PCA results, this is not the case -the disparities hint to divergent conformer distributions. Though it is perhaps unexpected that a single methyl group epimerization leads to considerable change in the conformer distribution of the entire ring system, the significantly different NMR spectra and datasets corroborates with this finding. By contrast, the oxolane 5-ring, the oxane 6-ring and the rhamnose ring displayed no significant difference comparing 2p and 11-epi-2p.
The resulting structural ensembles are the product of MD simulations driven by a rich set of NMR data, including NOEs, J-couplings and RDCs. Since these interactions act in reality on different time scales, no attempts were made here to extract kinetic information from these results. However, since MD simulations generally include energetic and entropic influences, these MDOC results propose that it is possible to extract a detailed picture of the conformer distribution and equilibrium in solution. The high consistency of these dynamic results with their corresponding data leads to the interrogation whether MDOC could be confidently used towards determination of configuration in flexible molecules such as natural products.

Configuration Determination with MDOC
It would be interesting to evaluate the diastereomeric discrimination potential of MDOC for flexible molecules. MDOC was already introduced for this task in relatively rigid and slightly flexible molecules [32,33], but in the present case, it is expected that orientation pseudo-forces may arrive more easily at false positives (e.g., structures with quality criterion >1) especially if the molecule's topology allows for large geometrical fluctuations due to increased number of degrees of freedom. On the other hand, the pseudo-forces from short-range NOEs and J-couplings may play an important role in hindering unrealistic ensembles. The test molecules at hand are ideally suited for this evaluation, since two diastereomers are available and can be used reciprocally for cross-validation.

Discriminating 1,4-Diketone 1-SS vs. 1-SR
As seen in Figure 12, all NMR parameters extracted from the MDOC trajectory run with dataset 1A scored better according on three criteria (n/χ 2 , 1/χ 2 min and F ) with the corresponding 1-SS configuration and, conversely, those with dataset 1B scored better (or equally well (1/χ 2 min )) with the 1-SR form. More importantly, the wrong configuration does not pass the quality criterion in both cases. The individual parameter types taken individually (1-bond RDCs, long-range RDCs, NOE distances and 3 J-couplings) also follow the trend that the criteria are consistently better for the correct configuration (Table S6 and Figure S12 in Supporting Information). The only exception regards the 3 J-couplings of dataset 1A for which the three values are satisfied almost equally for both diastereomers. Note that the three 3 J-couplings entered as constraints into the simulation influence the ω dihedral angle, and that these values are more characteristic of the conformational exchange than of the configuration. By contrast, the NOE distances seems to be highly sensitive to the configuration.  (Table S6 and Figure S12 in supporting information). The only exception regards the 3 J-couplings of dataset 1A for which the three values are satisfied almost equally for both diastereomers. Note that the three 3 J-couplings entered as constraints into the simulation influence the ω dihedral angle, and that these values are more characteristic of the conformational exchange than of the configuration. By contrast, the NOE distances seems to be highly sensitive to the configuration. Figure 12. The MDOC outcome based on datasets 1A and 1B against the correct (green) and incorrect (red) diketone 1 configuration. Displayed are the three criteria (quality criterion n/χ 2 , outlier criterion (1/χ 2 min) and data validity ℱ (labels)). For the breakdown of these parameters based on data type, see Figure S2 in the Supporting Information).

Discriminating between 4 Stereoisomers of Mandelalide A 2
Next, cross-validations of the datasets for the mandelalide A isomers 2 were tested in MDOC simulation against the incorrect configurations to address the diastereomeric discriminating ability of MDOC in more complex flexible molecules with various chiral centers. To this end, MDOC simulations with identical parameters (see SI) were applied to four configurations 2p, 11-epi-2p, 2r and 11-epi-2r, and the accuracy of the simulated NMR parameters were examined based again on the three aforementioned criteria. The results are displayed in Figure 13. Although the quality criterion is met (>1) for all configurations, the isomers 2r and 11-epi-2r with their inverted northern part show a substantial drop in the quality of fit relative to the correct ones (2p and 11-epi-2p, respectively). However, the inversion of the methyl group at position 11 does not lead to a large deterioration of the quality -especially with NMR dataset 2p (4.64 to 4.20). Closer inspection however showed that the small drop in quality is caused mainly by a number of long-range NOE distances, whereas the Figure 12. The MDOC outcome based on datasets 1A and 1B against the correct (green) and incorrect (red) diketone 1 configuration. Displayed are the three criteria (quality criterion n/χ 2 , outlier criterion (1/χ 2 min ) and data validity F (labels)). For the breakdown of these parameters based on data type, see Figure S2 in the Supporting Information).

Discriminating between 4 Stereoisomers of Mandelalide A 2
Next, cross-validations of the datasets for the mandelalide A isomers 2 were tested in MDOC simulation against the incorrect configurations to address the diastereomeric discriminating ability of MDOC in more complex flexible molecules with various chiral centers. To this end, MDOC simulations with identical parameters (see SI) were applied to four configurations 2p, 11-epi-2p, 2r and 11-epi-2r, and the accuracy of the simulated NMR parameters were examined based again on the three aforementioned criteria. The results are displayed in Figure 13. Although the quality criterion is met (>1) for all configurations, the isomers 2r and 11-epi-2r with their inverted northern part show a substantial drop in the quality of fit relative to the correct ones (2p and 11-epi-2p, respectively). However, the inversion of the methyl group at position 11 does not lead to a large deterioration of the quality -especially with NMR dataset 2p (4.64 to 4.20). Closer inspection however showed that the small drop in quality is caused mainly by a number of long-range NOE distances, whereas the simulated RDCs and J-coupling could not readily distinguish either configurations. In contrast, a more distinct drop in quality is observed for the 11-epi-2p dataset applied to configuration 2p. simulated RDCs and J-coupling could not readily distinguish either configurations. In contrast, a more distinct drop in quality is observed for the 11-epi-2p dataset applied to configuration 2p. Figure 13. Quality criteria (n/χ 2 , 1/χ 2 min and ℱ) for MDOC simulations of four configurations with two different NMR data sets. The configuration that belongs to the data is indicated in green.
From these results, it can be stated that in complex flexible molecules with various rotatable bonds and chiral centers, the MDOC-based discrimination of configurations on the basis of one-bond RDCs alone will be impracticable in a majority of cases; however, with the inclusion of additional high quality parameters such as NOE distances, some level of success should be expected. The example of mandelalide A already shows potential in this regard. With the refinement of the method and with the use of advanced statistical tools, it is conceivable that MDOC could provide a robust and general approach to assist in the difficult task of determining stereochemistry in larger natural products.

Conclusions
In this report, the MDOC method is proposed as a tensor-free approach for the conformational analysis of flexible molecules based on RDCs measured at high resolution in orienting media. One advantage of the MDOC method is that a vast range of interconversions are efficiently driven by the action of orientational pseudo-forces derived from tensorial properties, so that large compatible ensemble of conformers are generated. Since the expression of each orientational pseudo-energy possesses multiple minima, it is expected that the simulated ensemble may include non-native conformations, leading to inaccurate population equilibria. It is however another advantage of the MDOC method that it can be combined with traditional short-range scalar constraints like NOE distances or 3 J couplings, which are easily measured in solution and which should limit the generation of unrealistic rotamers. The entire NMR dataset can be brought together in a single simulation to generate a consistent picture regarding the structure of molecules in solution and their Figure 13. Quality criteria (n/χ 2 , 1/χ 2 min and F ) for MDOC simulations of four configurations with two different NMR data sets. The configuration that belongs to the data is indicated in green.
From these results, it can be stated that in complex flexible molecules with various rotatable bonds and chiral centers, the MDOC-based discrimination of configurations on the basis of one-bond RDCs alone will be impracticable in a majority of cases; however, with the inclusion of additional high quality parameters such as NOE distances, some level of success should be expected. The example of mandelalide A already shows potential in this regard. With the refinement of the method and with the use of advanced statistical tools, it is conceivable that MDOC could provide a robust and general approach to assist in the difficult task of determining stereochemistry in larger natural products.

Conclusions
In this report, the MDOC method is proposed as a tensor-free approach for the conformational analysis of flexible molecules based on RDCs measured at high resolution in orienting media. One advantage of the MDOC method is that a vast range of interconversions are efficiently driven by the action of orientational pseudo-forces derived from tensorial properties, so that large compatible ensemble of conformers are generated. Since the expression of each orientational pseudo-energy possesses multiple minima, it is expected that the simulated ensemble may include non-native conformations, leading to inaccurate population equilibria. It is however another advantage of the MDOC method that it can be combined with traditional short-range scalar constraints like NOE distances or 3 J couplings, which are easily measured in solution and which should limit the generation of unrealistic rotamers. The entire NMR dataset can be brought together in a single simulation to generate a consistent picture regarding the structure of molecules in solution and their dynamic equilibrium. For the 2p isomer of mandelalide A, for instance, it can at least be stated that the resulting ensemble are in accordance with no less than 212 NMR constraints. Further investigations will show to what extent this method can be generalized.
The cross-validation applied on different configurations of the same molecules offered the opportunity to show that MDOC has the ability to discriminate between diastereomers although the level of contrast depends on the quality of the NMR data and on how strongly the stereogenic centers interact with the rest of the structure and induce overall changes in the conformer distribution. In the case of the relatively small diketone pair, 1-SS and 1-SR, differing at a single chiral position, the MDOC results fit the incorrect diastereomer quite well by conventional standard testing (e.g., RMS or Q), but according to the n/χ 2 and 1/χ 2 criteria, the outcome is much clearer. In the case of the much larger mandelalide A, configuration 2p can readily be distinguished from 2r but in the case of the epimers 2p and 11-epi-2p the difference is barely significant. Nevertheless, all 4 configurations pass the n/χ 2 test and have similar in 1/χ 2 values. This is not surprising since a majority of the short-range constraints (NOEs and J) are only weakly affected by the chiral difference at distant positions. It would thus be desirable to take these aspects into consideration and develop more advanced statistical tools in order to confidently use MDOC for configurational evaluations in flexible molecules.
As a consequence, MDOC may become a method of choice for the NMR investigation of flexible molecules where the structure and distribution of conformers are of interest. To solve these questions can be of high interest in all areas of chemistry where conformational selection is of interest such as in the study of interactions of biological active molecules with their targets or in the mechanistic understanding of stereoselective catalytic reactions.

Materials
The monomers, methyl methacrylate (MMA, 99%, Sigma-Aldrich, St. Louis, MO, USA) and ethylene glycoldimethacrylate (EGDMA, 98%, Sigma-Aldrich) were purified prior to the experiments by passing the neat liquids through a short column filled with basic alumina in order to remove the polymerization inhibitor. The radical initiator 2,2'-azobis(2,4-dimethyl-4-methoxyvaleronitrile (V-70) was purchased from Wako (Neuss, Germany), and acetone-d 6 and chloroform-d 1 (99.9% of D atoms) were purchased from Cambridge Isotope Laboratories (Tewksbury, MA, USA). Twenty five mg of a single diastereomer of the 1,4 diketone was kindly donated by the List group, and was a product of an asymmetric cross coupling catalysis reaction [34]. About 2 mg of "pseudo"-mandelalide and 11-epi-"pseudo"-mandelalide (recognized as isomers of mandelalide A as a result of its synthesis, following the inversion of the northern part) were obtained from the Fürstner group as the end product of a multistep synthesis by Willwacher et al. [41] 4.1.2. Sample Alignment Alignment of the diketone and of the mandelalide was achieved using reversible compression of polymethylmetacrylate (PMMA) gels prepared as described by Gayathri et al. [42]. PMMA gel sticks of about 2 mm in diameter and 25 mm long and crosslink density of 0.3 mol% were pre-swollen in chloroform-d 1 , inserted into a 5-mm NMR tube. Residual monomers were washed out of the polymer stick gel by applying several compession-and-release cycles in the deuterated solvent (chloroform-d 1 ) using a home-made Teflon plunger. This was repeated until no monomer could be observed in the 1D 1 H NMR spectrum. The gels exhibited quadrupolar splitting (∆ν Q ) of the solvent signal of form 40 to 47 Hz when maximally compressed, and no quadrupolar splitting (∆ν Q = 0) when fully relaxed as assessed by 1D 2H NMR. The homogeneity of the alignment was also assessed using 2 H-mapping approaches [43]. The dissolved analytes (diketones or mandelalides) were dispersed into the gel using a pumping actions as described above. Samples could be measured immediately.  [44] were recorded to measure the CH splittings corresponding to scalar couplings ( 1 J CH ) or to scalar plus residual dipolar couplings ( 1 T CH = 1 J CH + 1 D CH ) in the isotropic and anisotropic samples. This experiment was typically parametrized with 16k x 512 points in the acquisition matrix, 8 scans per increment and 3-s relaxation delay. A homo-decoupled CLIP-RESET-HSQC [45] was also measured in the case of the mandelalide sample as a 3D matrix of 512 × 512 × 12 (chunks) points, with chunk duration of 16.255 ms, 8 scans per increment and relaxation delay of 1 s. All raw, processed and analysed data are made available as NMReDATA records [46].

Conformational Analysis
The initial input structures for the two possible diastereomers of diketone 1 (1-SS and 1-SR) in three different rotamers (trans, gauche+ and gauche-) were generated with the hybrid density function B3LYP, cc-pVTZ as a basis set in Gaussian-09. For the conformation analysis of these structures against the measured RDCs based on the classical alignment method using SVD, the program Mspin (Mestrelabs, Santiago de Compostela, Spain) was employed. All inputs are given in the Supporting Informations.

Theory: RDC-Based Orientation Constraints
The general methodology of MDOC has already been outlined [30,31] and described more specifically for RDCs [32,33]. Since new methodological features are introduced in this paper some essential points of the method are discussed in this section.
Prerequisite of MDOC is a molecular mechanics force field that is flexible enough to calculate the relative energies of most organic molecules and provide structures that compare well to diffraction experiments or more elaborated ab initio or DFT calculations. The COSMOS-NMR [47,48] force field that was used in this case has one distinctive advantage over most other force fields -it uses partial atomic charges from a quantum chemical method [48][49][50] (Bond Polarization Theory -BPT) to calculate the electrostatic energy. Since these charges can be recalculated in the course of an MD simulation, all mutual polarization can be included into the electrostatic energy. This particular performance of the force field is an issue especially if orientational constraints are applied, since their ambiguities/degeneracies have a tendency to drive the MD towards unrealistic structures. By contrast, MD with pseudo-forces based only on scalar constraints are more self-contained.
The pseudo-energy in the case of tensorial properties like dipolar couplings has a special form: The first double sum (denoted with αβ -the indices α and β are used to denote the coordinates x, y and z) in Equation (1) runs over all nine components of the dipolar tensor D and all tensor components are regarded as constraints. Other pseudo-energies for scalar constraints like NOE distances or scalar J couplings have only the second summation over the experimental values (denoted with i). The constant k is used to convert the expression into kJ/mol and to adjust the strength of the pseudo-forces. Measured RDC values represent the z-components of diagonal traceless tensors whereas the off-diagonal elements are averaged to zero by the rapid reorientation of the molecule around the director. The measured tensors D exp are given in the laboratory coordinate system, oriented with its z-axis relative to the direction of the external magnetic field. The calculation of dipolar coupling tensor D theo can easily be performed in the coordinate system that is oriented parallel to the vector that connects the two coupled nuclei. In this case the dipolar tensor is diagonal and its principal values can be calculated from the product of the gyromagnetic ratios of the two nuclei divided by the third power of their distance (also see Supplementary Information). To calculate the differences between the calculated and experimental tensor components, the tensor D theo has to be transformed from its principal axis system to the laboratory coordinate system using proper transformation matrices T i : In Equation (2) the dipolar coupling of 1 H and a 13 C nucleus in a C-H bond is chosen as example. This coupling is nearly the same for all C-H bonds since the bond distance is nearly constant (in the present simulations, different values only for C(sp 3 )-H and C(sp 2 )-H couplings are used-see Supplementary Information). The double sum in Equation (2) runs over all components of the transformation matrices T i being mostly different for all sites i because of the different orientation of the bonds with respect to the laboratory system.
The dipolar tensor components obtained from Equation (2) cannot be directly compared with experimental tensor components since they vary within several tens of kHz whereas the experimental RDC have mean values that span only several Hz. Therefore a theoretical mean value is calculated during the time of the MD simulation. For this mean value calculation, a special form is chosen that implies a memory time τ: The exponential decaying memory function was first introduced by Torda et al. [51] and leads to a time-dependent property on the time scale of the decay or memory time τ. In the context of MDOC, this time τ can be considered as life time of an orientation of a molecule or molecular part. Since the pseudo-forces depend on this lifetime, τ should be considered as the time scale for orientational and conformational changes. After 5τ, less than 1% of the original orientation as expressed by the instantaneous value of D is left and the pseudo-forces change on the same time scale τ because of the mismatch of calculated and experimental data. In this respect, MDOC can be regarded as an accelerated MD simulation and orientational averages can be reached in moderate simulation times. A value of τ of 200 picoseconds was generally used whereas the simulations were performed over durations one or two order of magnitudes longer than τ − i.e., over several nanoseconds.
Since in MD simulations the integration of the equations of motion is performed in time steps, the integration in Equation (3) can be transformed into a summation. To speed up the calculation a recursion formula for this integral was developed (see Supplementary Information). The formalism was used for all types of tensorial and scalar constraints as for instance for NOE distances and scalar couplings. In MD simulations, the energies and pseudo-energies (Equation (1)) are not directly used but the forces that are calculated from these energies as first derivatives with respect to the coordinates of the nuclei A i and B i : From Equation (4), the magnitude of the forces on the coupling atoms grows linearly with the difference between experimental and calculated RDC. This dependence occurs for harmonic potentials like Equation (1) but for other molecular potentials like those of bonds harmonic behavior occurs only for deviations not far from the minimum. For large deviations as in the case of pseudo-forces harmonic potentials lead to very high pseudo-forces and consequently to unrealistic deformations of the structure. Therefore, the difference D Equation (5) represents a hyperbolic tangent dependence that becomes constant if the difference between calculated and experimental RDC gets much larger than a width parameter ∆D. This width parameter can be evaluated as the experimental error of the RDC.
Inserting Equation (2) into Equation (4) reveals that orientational pseudo-forces imply derivatives of the transformation matrices T i in Equation (2) with respect to the coordinates coupling nuclei. The transformation matrices can be obtained from sets of orthogonal unit vectors that span the local bond (or interconnection) coordinates systems and therefore the derivatives can be obtained from these unit vectors (see Sternberg et al. [31]). In contrast to 1-bond RDCs( 1 D), long-range RDCs ( n D) can also be used as distance constraints in analogy to NOE distances, providing the mean distribution of orientations of nuclear connection vectors (or alignment tensor) is known. However, in this investigation, only the orientational derivatives of the RDC pseudo-energies are used.
It should be remarked that the magnitude of the pseudo-forces is controlled by the differences between experimental values and mean values calculated according to Equation (3) but the direction of the orientational pseudo-forces is calculated from derivatives of the actual transformation matrices (see Equation (2)). In this way a mismatch in orientation acts immediately. The problem is that at start of the MDOC, there are no mean values and the strong pseudo-forces would completely distort the molecule. For this reason, the pseudo-forces are introduced in exponential fashion using 1 − e −t/ρ with ρ the time constant for the rise of the pseudo-forces. This value was set in the simulation to the same value as the memory time τ (see Equation (3)).
In former investigations, the orientation and motions of molecules in lipid membranes were studied. These liquid crystalline systems display a high degree of order that is expressed using the so called Saupe order parameter S (for review see Limmer [52]). The molecules in these investigations reached order parameters up to S~0.8. In the present investigation weakly oriented media such as stretched gels were used allowing order parameters in the range of 10 −2 to 10 −3 . Therefore an order parameter S am is introduced into the calculations that accounts for the reduction of the dipolar splitting caused by the alignment medium (in analogy to the order parameter S bilayer introduced by Marsan et al. [53]). S am is multiplied to the calculated dipolar couplings to prevents too high pseudo-forces by shifting the dipolar splittings from the kHz range to several Hertz. S am should be chosen large enough that its product with the maximum splitting is in any case much larger than any observed splitting. If, for example, the static C-H splitting of 47.96 kHz is multiplied with S am (0.004, see Table S10 in Supplementary Information) a maximum attainable CH splitting of 191.8 Hz is obtained, that is much larger than any observed splitting.

Theory: Scalar Constraints
In addition to the RDC, NOE distances and 3 J HH indirect scalar couplings are used as constraints. For the calculation of the time mean values in both cases the average with an exponential memory as given in Equation (3) was used. The pseudo-forces are augmented with hyperbolic tangent weight function as given in Equation (5). For the calculation of the distances the following average is applied [54]: (the mean value calculation is indicated using the symbol < >). Equation (6) is valid for rigid molecules that undergo isotropic reorientations [49] In MDOC simulations, it is possible to directly use NOE distances of protons of CH 3 groups as constraints. If for the CH 3 protons also one bond RDC are used as constraints these groups rotate fast in the MDOC simulations and therefore the same distance constraint can be applied to all three protons of a CH 3 group.
For the calculation of the 3 J HH couplings the equation according to Haasnoot et al. [4] was used. In this Karplus-like equation a correction term is added that accounts for the influence of the electronegativity of neighbor substituents on the coupling protons. Using this equation an RMS (root mean square) deviation between calculation and experiment in the range of ca. 0.6 Hz should be possible.

Evaluation of Quality
In many investigations, the RMS deviation between calculated and measured values is used to estimate the performance of theoretical or computational methods. The disadvantage of the RMS criterion is that (i) the valuable information about the error of the measurement is disregarded and (ii) since the RMS value has the unit of the measurement different data sets as for instance distances and frequencies cannot be compared properly. Therefore the n/χ 2 criterion is used, as introduced by Intelmann et al. [55]: The quality of a computation is expressed as n/χ 2 (with n the number of measured values of the property P) giving values larger than 1.0 when the calculations are on average within the experimental error bounds.
As introduced by Tzvetkova et al. [32], a second tighter indicator for the quality of a computed dataset focuses on the outliers. This outlier criterion is formulated as 1/χ 2 min (i.e., the smallest value among all (1/χ 2 i ). A dataset is perfectly reproduced by a computed ensembles within error margins, when this the condition 1/χ 2 min > 1 is met. Finally, to convey the overall fidelity of a computed dataset, the ratio of valid data F = (n − n outliers )/n will also be employed.