NMR free ligand conformations and atomic resolution dynamics

Knowledge of free ligand conformational preferences (energy minima) and conformational dynamics (rotational energy barriers) of small molecules in solution can guide drug design hypotheses and help rank ideas to bias syntheses towards more active compounds. Visualization of conformational exchange dynamics around torsion angles, by replicaexchange with solute tempering molecular dynamics (REST-MD), gives results in agreement with high resolution H NMR spectra and complements free ligand conformational analyses. Rotational energy barriers around individual bonds are 10 comparable between calculated and experimental values, making the in silico method relevant to ranking prospective design ideas in drug discovery programs, particularly across a series of analogues. Prioritizing design ideas, based on calculations and analysis of measurements across a series, efficiently guides rational discovery towards the 'right molecules' for effective medicines.

Determining the conformational profile of a free ligand in solution enhances early drug discovery programs (LaPlante et al., 2014;Blundell et al., 2016;Foloppe and Chen, 2016;Chiarparin et al., 2019). A general overview of how NMR fits into Medicinal Chemistry design cycles is illustrated in Fig. 1. An efficacious pharmaceutical that positively impacts patients' lives starts with Medicinal Chemistry teams designing the right molecule. Design teams need to understand whether a molecule readily adopts its "bioactive" conformation in solution to optimize the binding on-rate through reduced 35 conformational entropy and energetic penalty paid on conformational rearrangement to the proper binding mode. In addition, pre-organization of the free ligand in solution would indicate minimized conformational strain energy in the bound molecule.
If not, the challenge is to conceive of ideas to modify the molecule to favor this conformation. Towards this aim of optimizing the free energy of binding, it is desirable for ligands in solution to preferentially pre-organize into the bioactive binding mode (Blundell et al., 2013;Balazs et al., 2019). Molecular rigidification strategies (Fang et al., 2014;de Sena M 40 Pinheiro et al., 2019) increase pre-organization and NMR conformational analysis can deconvolute and report on the molar fraction adopting the bioactive conformation. Structure based drug design (SBDD) can be enhanced by ready access to 3D free ligand average solution conformations to complement X-ray crystallographic models of the bound ligand and proteinligand interactions (Blundell et al., 2013;Chiarparin et al., 2019;Balazs et al., 2019). Faster design cycles require quick turnover times in analyzing solution conformations of synthesized compounds. Design cycles can be accelerated through 45 faster computational schemes, efficient automation to obtain NMR spectral parameters, and recognition of conformational signatures from 1D NMR spectra (Balazs et al., 2019), also named, "SAR by 1D NMR" (Zondlo, 2019).  NMR plays a key role in synthesis support for structural identification and analysis of compound purity. NMR also can be used to enhance structure based drug design through NMR free ligand conformational analysis. This provides the relative population of the bioactive conformation in solution by determining the percentages of the minimal energy conformers. Structure based drug design is enhanced through measured extent of free-ligand pre-organization into the bioactive conformation, which lends itself to rigidification design hypotheses aimed at optimized binding to the protein target.
Herein, we demonstrate incorporation of molecular dynamics, specifically an efficient version using replica-exchange with 55 solute tempering (REST-MD) (Liu et al., 2005;Huang et al., 2007;Wang et al., 2011), into an NMR based semi-automatic drug discovery platform, to visualize rotational barriers around molecular bonds. Good agreement is demonstrated between REST-MD calculated energy barriers and NMR measurements, using a small molecule selective estrogen-receptor degrader (SERD) example from a recent Oncology R&D project (Scott et al., 2016;Scott et al., 2019;Scott et al., 2020). The theoretical and experimental data complement each other: REST-MD simplifies the interpretation of NMR conformational 60 dynamics, while the experimental NMR results can inform calculations by defining site-specific preferred torsions of the dominant conformer and experimental conformer distributions, which may influence the initial REST-MD 3D geometry and the sampling ergodicity achieved, as reflected in the resultant histograms.

NAMFIS plus NMR line shape analysis 65
The ability of NMR to provide information on conformational dynamics, in addition to giving information on conformational preference, is useful in small molecule drug discovery. In Fig. 2a, each peak is doubled for compound 2 (1:1 ratio), instantly recognizable to NMR users as slow exchange of rotamers due to hindered bond rotation (measured half-life ~ 0.5 sec). As a guide to the eye, the signal(s) for the benzylic CH proton at ~ 5 ppm is/are highlighted in Fig. 2. The NMR spectrum reports two dominant conformers, equally populated, for the free ligand in solution for compound 2. The bioactive 70 conformation is one of a family of conformers (shown in green) with some flexibility around the pendant base. The alternative conformation (shown in orange) is the other, giving ~ 50% of the compound locked in a non-bioactive conformation. In contrast, Fig. 2b shows that compound 3 has a single set of sharp signals due to fast exchange (corresponds to a typical half-life of ~ nanoseconds, Δv½ = 2.8 Hz). The 1 H NMR spectrum has a single isotropic chemical shift for both of the protons within one CH2 functional group (these are not diastereotopic), an indication of local flexibility quickly picked 75 up by an edited 13 C HSQC spectrum. Appreciation of the temporal dependence of free ligand exchange dynamics on NMR spectra can quickly inform medicinal chemists on local flexibility around bonds of newly synthesized molecules. This analysis, combined with potency data and matched molecular pair analysis, or a full comparison across a congeneric series, can provide critical insights into structure-activity relationships (Balazs et al., 2019).
Together with information regarding relative populations of conformations in solution, information about the magnitude of 80 the rotational energy barrier between conformations, i.e. between one rotational isomer and another, is relevant information. The challenge has been to get quick, easy, and comprehensive, yet accurate torsional profiles. Building a practicable and prospective visualization of conformational exchange dynamics around torsion angles into an NMR conformational analysis platform increases the potential to impact design, by highlighting the potential energy penalty of restricting torsions. To evaluate the dynamics, incorporating REST-MD into the workflow met the goal of expanding the current free ligand 85 conformational analysis platform to make use of kinetic parameters from NMR spectra, e.g. signal line widths, while keeping within practical time limit constraints for Medicinal Chemistry design cycles.   function of the simulation time (radial plots, from the center at the start and spiraling outward) or as histograms superimposed on the torsion energy profiles (kcal/mol vs bond angle across each color coded bond in the molecule). Profiles are readily compared between molecular bonds that are pre-organized (light blue with y-axis maximum in the plot at 25 kcal/mol, with a bimodal radial plot and two energy minima), and flexible (pink with 3 kcal/mol maximum y-axis value, three energy minima, equally populated, and a randomly populated radial plot). (b) The barrier to rotation of the dimethyl is calculated to be, based on the lower of the two 115 barriers, ~20 kcal/mol. Whereas experimentally both energy minima are equally populated as seen by the 1:1 ratio by NMR, the sampling conditions of the rigid "blue" bond were insufficient in the simulation to equally populate both wells. The NMR data in such cases clearly informs on the calculated predictions. A separate REST-MD simulation for the des-methyl compound, 3, was < 6 kcal/mol calculated rotational energy barrier for the same "blue" bond, consistent with sharp lines and ready conversions between the two conformers, with a broadened range of torsions, albeit still bimodal.

NMR measured rotational energy barrier
Methylation is a familiar and fundamental structural rigidification tool in a Medicinal Chemist's toolbox. In bond rotation on the NMR timescale (~ nanoseconds). Typically such barriers to rotation at room temperature correspond to ~ 5 kcal/mol (LaPlante et al., 2011a;LaPlante et al., 2011b;Wipf et al., 2015). The NMR spectrum of the ensemble of rapidly exchanging conformations reflects the Boltzmann weighted average of the chemical shifts, J-couplings, and interproton distances, with a single set of sharp peaks.
To locate the bond with the hindered rotation, chemical intuition is often sufficient. Using variable temperature NMR and/or 130 exchange spectroscopy, the rotational energy barrier and the torsional rotation half-life of exchange can be determined. Fig.   4 shows 1 H NMR spectra as a function of eight different temperatures. The spectrum near room temperature has two equal rotameric populations undergoing slow exchange on the NMR timescale, and highlighted in the figure. With increasing temperature the peaks coalesce and then begin to narrow. Increasing the temperature not only increases the rotation rate of the aromatic ring, it also increases the rate of fluctuations in the pendant base and the CH2CF3 groups and between axial vs. 135 equatorial methyl orientation in ring C. Overall, this drives a shift to higher ppm for the exchange averaged signal with increasing temperature (Fig. 4). In order to determine the barrier to rotation around the aromatic ring, it was important to collect data within a temperature range where exchange rates were dominated by the dynamics of interest in order to follow a simple two-state model for analysis. Therefore, three temperatures at 300, 305, and 310 K were chosen and exchange spectroscopy was performed with 145 selective inversion on the peak near 5.3 ppm. At each of the three temperatures 8 mixing times (100, 200, 300, 400, 500, 700, 1000, and 2000 milliseconds) were used to determine the first order rate kinetics, with the fitted value for k given in tabular form in Fig. 4. The half-life was derived from ln(2)/k. The fitted plot of ln(k) vs. 1/RT is shown, revealing a value of 19.9 kcal/mol for the barrier to rotation.

NMR informs calculations 150
REST-MD generates a large ensemble of (~1000) conformations, in explicit solvent. REST-MD was run with Desmond (Schrödinger, 2020a) with the pendant base initially oriented either forward or backward relative to the tricyclic core for 1.
The resultant calculated energy barrier of ~ 20 kcal/mol (Fig. 3b) is in agreement with the NMR determined value of 19.9 kcal/mol. The REST-MD visualization confirmed chemical intuition that the source of the rotameric species is the bond between the tricyclic core and the aromatic ring. Advantageously, the full torsional profile report from the REST-MD 155 simulation can be run prospectively to rank design ideas, for instance to test a hypothesis around rigidification and the degree of bioactive pre-organization induced. The ability of REST-MD to evaluate torsion angles prior to synthesis can also flag chemists to check the NMR for site-specific dynamics information at the time of structural verification. Such information could alert the team to add a diagnostic selective-NOE measurement to the standard acquisition suite, to test a free ligand conformational hypothesis post synthesis, while the solution sample is in the spectrometer for structural identification. 160 Conversely, the NMR can supply experimental details inaccessible to the calculations, particularly helpful within a chemical series, as the lessons are generally translatable across the structural analogs. For instance, the 1 H NMR spectrum of the dimethyl compound 2 showed a 1:1 ratio between the two exchanging conformations (Fig. 2). Whereas REST-MD conformational sampling shows only one of the two minima populated (Fig. 5). The fragment based energy calculation, shown as a solid line in the REST-MD torsional profiles, is consistently the same, even using a very short simulation time 165 (e.g. picoseconds). The histograms vary based on initial conformation and number of replicates run. Starting with an initial conformation with the pendant base facing "forward" relative to the tricyclic core, the radial plot of torsion angle representation as a function of time starting at the center and spiraling outward, only populates the ~ +90º bond torsion angle during the 50 ns simulation that has a temperature range of 300 -1263 K (12 replicates, 50 nanoseconds). Analogously, the histogram has one energy minimum populated and the number of times the ~ +90º torsion was observed is fairly narrowly 170 distributed (Fig. 5, top). With the same initial conformation, increasing the temperature range to a high of 3302 K (20 replicates, 50 nanoseconds), showed some evidence of sampling of the opposite conformation in the radial plot (Fig. 5,   middle). Starting the simulation with the 3D conformation switched to put the pendant base towards the back instead and running 20 replicates for a higher sampling temperature shows both minima were sampled (Fig. 5, bottom). Once this compound has been synthesized, it then becomes experimentally apparent from the 1H NMR spectrum at 300 K in DMSO-d6 175 that both minima are equally populated (Fig. 2a). In this manner the experimental results can be fed-back to the calculations to refine details and gauge areas of caution during interpretation.

REST-MD/NMR synergy in drug discovery
What REST-MD adds to the existing NMR platform is visualization of conformational dynamics, by providing calculated rotational energy barriers across all bonds. This complements NMR spectral data to give insight into flexibility / rigidity at atomic resolution. Together, REST-MD and NMR conformational analysis allows us to utilize all the spectral information, thermodynamic and kinetic, gathered from 1 H NMR spectra: chemical shifts, J-couplings, NOEs and linewidths, to energy penalty paid if the bioactive conformation is not highly populated in solution. This can help rationalize the cost to benefit ratio of effort invested in designing an increase in the percent bioactive conformation by restricting rotation. As drug discovery requires optimization of several parameters, knowing when binding has been optimized can shift design focus and resources towards improving physicochemical properties.
The added benefit of REST-MD is its ability to deliver prospective information regarding structural ideas of compounds, not 200 yet synthesized. Accurate predictions of free ligand solution conformational dynamics can help rank compounds to focus synthesis prioritization and flag supplemental experiments, such as selective NOEs on atom pairs to quickly ascertain expected conformations.
While the full torsional profile is powerful on its own, the REST-MD results also provide an extensively sampled conformational ensemble in explicit solvent that can be clustered and fed forward for use in NMR conformational analysis. 205 Taking the selected conformer set forward for QM geometry refinements and calculations of NMR chemical shift and coupling constants provides the modeled parameter set used for further NAMFIS based analysis. Using the conformer set generated by REST-MD is particularly helpful for higher molecular weight small molecules which begin to self-associated during low mode MD conformational searches using a polarizable continuum model to emulate solvent, resulting in a set highly biased towards collapsed conformations. 210

NMR Spectroscopy
1 H NMR spectra were recorded at 300 K on a 500 MHz NEO or a 600 MHz AVIII Bruker spectrometer with TCI cryoprobes. Solutions were made from 1-4 mg solid freshly dissolved in DMSO-d6. Spectra were acquired with a 30 degree hard pulse, a 1 sec delay, 2 dummy shots, and signal averaged over 16 transients. A spectral width of ~ 20 ppm with 64k 215 points was used. Spectral analysis was performed using Advanced Chemistry Development, Inc. (

Exchange Spectroscopy
1 H NMR spectra were recorded on a 500 MHz NEO at 300,305,310,340,345,350,355,360,365,and 373 K. For the 1D selective exchange spectroscopy at 300, 300, and 310 K, the mixing times used were 100, 200, 300, 400, 500, 700, 1000, and 2000 milliseconds. The spectra were integrated with consistent integral ranges and calibrating the integral of the inverted 245 peak to 100 to consistently normalize the data (Hu and Krishnamurthy, 2006). An excel spreadsheet was used to calculate the fractional intensity increase as a function of mixing time to fit exchange rate and half-life (Bovey, 1988;Li, et al., 2007).

REST-MD
Two different initial molecular conformations were run where the variation was place on the relative position of the pendant base to the tricyclic ring: (i) "forward" or (ii) "backward", using the same atom numbering for all conformations sampled for 250 the same compound, to simplify later steps in the workflow. Molecular protonation states at pH 7.0 ± 0.0 were used for the MD simulations. The force field builder in Maestro (Schrödinger, 2020b) was used to customize the OPLS3e force field for the system builder where a NaCl salt concentration of 0.15 M was used and the base was neutralized by addition of 1 Na+ ion during creation of the explicit water shell for solvation using the predefined SPC model and an orthorhombic box shape of 10 Å x 10 Å x 10 Å using the "buffer" box size calculation method. Desmond (Schrödinger, 2020a) replica exchange with 255 solute tempering molecular dynamics was run with 12 replicas giving a temperature range of 300 K to typically ~1300 K, for a total of 50 ns for extensive sampling of conformational space during the simulation.

Simulation Interaction Diagram
The plots automatically generated in Maestro (Schrödinger, 2020b) provide solid lines tracing out the barrier to rotation in kcal/mol as a function of the torsion angle. The histograms provide the resulting distribution of 1002 conformers under our 260 routine sampling conditions. Radial plots show the evolution of the simulation time from the start, at the center.

Clustering of Conformers
Ligands, without the solvent shell, were extracted from the REST-MD trajectories for both sets of initial conformers (forward and backward). To aid a quick visual inspection of the results, conformers was superimposed using the SMARTS method and the substructure smiles string of c12c(C)c(CN)[nH]c1cccc2 to align the conformers relative to the rigid tricyclic 265 ring. Conformers were clustered in Maestro (Schrödinger, 2020b) by atomic RMSD, discarding mirror-image conformers, selecting the option of one structure (nearest to centroid) per cluster, thus reducing the full set down to representative diverse conformers, typically 15 -40.

QM Calculations
Chemical Computing Group's (Molecular Operating Environment (MOE), 2019.01) conformational search GUI was 270 employed to generate input files for Gaussian 16 (Revision C.01) after importing conformers into a molecular database and using the wash function to neutralize charged species not observed by NMR in DMSO-d6 solutions. QM geometry refinement, chemical shift calculations and coupling constant calculations were carried out with the GIAO DFT method at the B3LYP/6-31G* level with PCM solvent modeling using a dielectric constant of 78.4. Geometry optimization keywords were set with opt=(tight,RecalcFC=5,MaxCycles=5000) and Int=SuperFineGrid. 275

Conformer Distribution
MOE's Spectral Analysis GUI was employed for least squares fits of chemical shifts to determine conformer distributions; the option for couplings and NOE's was used selectively.

Conclusions
The REST-MD protocols described above provide rapid and prospective access to torsional energy barriers and 280 conformational states for drug-like molecules. The REST-MD calculations accurately reproduce and visualize NMR dynamics which synergistically work with experimental conformational exchange dynamics obtained from 1D 1 H NMR spectra. Integration of REST-MD into our NMR conformational analysis platform has enabled visualization of atomic level information by all medicinal chemists and can be readily used to guide design hypotheses toward molecules with improved potency and or physicochemical properties. 285 This new methodology has been applied across more than 10 early oncology projects in 2020, both small molecules and PROTACs, to answer questions around conformational preference (populations) and dynamics (rotational barriers).
In addition, NMR provides design teams with information on the presence of intramolecular hydrogen bonds (IMHB), and the combined influences on properties such as potency, permeability and oral bioavailability. Diverse applications have enabled refinement of the approach, and represents a step towards the goal of routine use for prospective design and 290 determination of experimentally based conformation-activity relationships.

Author contribution
AB and EC drove the earliest drafts of this manuscript, to which all authors have contributed. MP developed and optimized computational workflows. DL and ND determined the energy barrier to rotation by NMR. AB ran REST-MD simulations in Schrödinger and NAMFIS-based analyses in MOE. All authors contributed valuable discussions to the preparation of this 295 manuscript.

Competing interests
All authors are shareholders in AstraZeneca PLC.