Toward Reproducing Sequence Trends in Phosphorus Chemical Shifts for Nucleic Acids by MD/DFT Calculations

Modeling solvent effects and convergence of 31P-NMR shielding calculations with COBRAMM

Solvent effects on 31P-NMR parameters for triphenylphosphine oxide and triphenylphosphine in chloroform have been extensively investigated by testing different solvation models. The solvent is described implicitly, mixed implicitly/explicitly, and using full explicit models. Polarizable continuum model (PCM), molecular dynamic (MD) simulations, and hybrid quantum mechanics/molecular mechanics (QM/MM) calculations are used to disclose the effects of solute/solvent interactions and, more generally, the role of the embedding in NMR simulations. The results show the beneficial effect of carrying out QM/MM optimizations on top of geometries directly extracted from classical MD simulations, used to ensure representative conformational sampling. The nuclear shielding convergence has been tested against a different number of snapshots and with the inclusion of solvent shells into the QM region. An automated MD//QM/MM//GIAO protocol, implemented in the COBRAMM package, is here proposed and tested on trimethyl phosphite showing that our approach boosts the convergence of nuclear shielding satisfactorily. The present work aims to be a stepping-stone to assess proper QM/MM computational strategies in simulating chemical shifts in non-homogeneous systems like supramolecular and biological systems.

Systematic QM/MM Study for Predicting 31P NMR Chemical Shifts of Adenosine Nucleotides in Solution and Stages of ATP Hydrolysis in a Protein Environment

NMR (nuclear magnetic resonance) spectroscopy allows for important atomistic insights into the structure and dynamics of biological macromolecules; however, reliable assignments of experimental spectra are often difficult. Herein, quantum mechanical/molecular mechanical (QM/MM) calculations can provide crucial support. A major problem for the simulations is that experimental NMR signals are time-averaged over much longer time scales, and since computed chemical shifts are highly sensitive to local changes in the electronic and structural environment, sufficiently large averages over representative structural ensembles are essential. This entails high computational demands for reliable simulations. For NMR measurements in biological systems, a nucleus of major interest is 31P since it is both highly present (e.g., in nucleic acids) and easily observable. The focus of our present study is to develop a robust and computationally cost-efficient framework for simulating 31P NMR chemical shifts of nucleotides. We apply this scheme to study the different stages of the ATP hydrolysis reaction catalyzed by p97. Our methodology is based on MM molecular dynamics (MM-MD) sampling, followed by QM/MM structure optimizations and NMR calculations. Overall, our study is one of the most comprehensive QM-based 31P studies in a protein environment and the first to provide computed NMR chemical shifts for multiple nucleotide states in a protein environment. This study sheds light on a process that is challenging to probe experimentally and aims to bridge the gap between measured and calculated NMR spectroscopic properties.

Improving the accuracy of 31P NMR chemical shift calculations by use of scaling methods

Calculation of ³¹P NMR chemical shifts for a series of tri- and tetracoordinate phosphorus compounds using several basis sets and density functional theory (DFT) functionals gave a modest fit to experimental chemical shifts, but an excellent linear fit when plotted against the experimental values. The resultant scaling methods were then applied to a variety of "large" compounds previously selected by Latypov et al. and a set of stereoisomeric and unusual compounds selected here. No one method was best for all structural types. For compounds that contain P-P bonds and P-C multiple bonds, the Latypov et al. method using the PBE0 functional was best (mean absolute deviation/root mean square deviation (MAD/RMSD) = 6.9/8.5 ppm and 6.6/8.2 ppm, respectively), but for the full set of compounds gave higher deviations (MAD/RMSD = 8.2/12.3 ppm), and failed by over 60 ppm for a three-membered phosphorus heterocycle. Use of the M06-2X functional for both the structural optimization and NMR chemical shift calculation was best overall for the compounds without P-C multiple bonds (MAD/RMSD = 5.4/7.1 ppm), but failed by 30-49 ppm for compounds having any P-C multiple-bond character. Failures of these magnitudes have not been reported previously for these widely used functionals. These failures were then used to screen a variety of recommended functionals, leading to better overall methods for calculation of these chemical shifts: optimization with the M06-2X functional and NMR calculation with the PBE0 or ωB97x-D functionals gave values for MAD/RMSD = 6.9/8.5 ppm and 6.8/9.1 ppm, respectively, over an experimental chemical shift range of -181 to 356 ppm. Due to the unexplained failures observed, we recommend use of more than one method when looking at novel structures.

Improving IDP theoretical chemical shift accuracy and efficiency through a combined MD/ADMA/DFT and machine learning approach

This work extends the multi-scale computational scheme for the quantum mechanics (QM) calculations of Nuclear Magnetic Resonance (NMR) chemical shifts (CSs) in proteins that lack a well-defined 3D structure. The scheme couples the sampling of an intrinsically disordered protein (IDP) by classical molecular dynamics (MD) with protein fragmentation using the adjustable density matrix assembler (ADMA) and density functional theory (DFT) calculations. In contrast to our early investigation on IDPs (Pavlíková Přecechtělová et al., J. Chem. Theory Comput., 2019, 15, 5642-5658) and the state-of-the art NMR calculations for structured proteins, a partial re-optimization was implemented on the raw MD geometries in vibrational normal mode coordinates to enhance the accuracy of the MD/ADMA/DFT computational scheme. In addition, machine-learning based cluster analysis was performed on the scheme to explore its potential in producing protein structure ensembles (CLUSTER ensembles) that yield accurate CSs at a reduced computational cost. The performance of the cluster-based calculations is validated against results obtained with conventional structural ensembles consisting of MD snapshots extracted from the MD trajectory at regular time intervals (REGULAR ensembles). CS calculations performed with the refined MD/ADMA/DFT framework employed the 6-311++G(d,p) basis set that outperformed IGLO-III calculations with the same density functional approximation (B3LYP) and both explicit and implicit solvation. The partial geometry optimization did not universally improve the agreement of computed CSs with the experiment but substantially decreased errors associated with the ensemble averaging. A CLUSTER ensemble with 50 structures yielded ensemble averages close to those obtained with a REGULAR ensemble consisting of 500 MD frames. The cluster based calculations thus required only a fraction of the computational time.

Probabilistic Interpretation of NMR J-Couplings Determines BI-BII State Equilibria in DNA

Interpretation of ³J_P,H3' NMR scalar spin-spin coupling constants in DNA becomes more reliable by including distinct structural states such as BI and BII, using the weighted-static or, better still, the recently implemented adiabatic-MD (Ad-MD) method. The calculation method employs an adiabatic ("Ad") dependence of ³J_P,H3' coupling on NMR-assigned torsion angle, ε, weighted by P(ε) probability distribution calculated by molecular dynamics (MD). Ad-MD calculations enable cross-validation of the bsc1, OL15, and OL21 force fields and various parametrizations of the Karplus equation describing the dependence of ³J_P,H3' coupling on ε torsion (KE). The mean absolute deviation of Ad-MD ³J_P,H3' couplings from the experimental values in Dickerson-Drew DNA is comparable to the scatter of ³J_P,H3' couplings among four separate NMR experiments. A commonly accepted assumption of homogeneity of one kind of structure-dynamic state within DNA (BI or BII) is questionable because the principal characteristics of relevant P(ε) probabilities (shapes and positioning) vary with DNA sequence. The theory outlined in the present work sets limits to future reparameterization of MD force fields, as relevant to NMR data.

The Ad-MD method to calculate NMR shift including effects due to conformational dynamics: The 31 P NMR shift in DNA

A method for averaging of NMR parameters by molecular dynamics (MD) has been derived from the method of statistical averaging in MD snapshots, benchmarked and applied to structurally dynamic interpretation of the ³¹ P NMR shift (δ_31P ) in DNA phosphates. The method employs adiabatic dependence of an NMR parameter on selected geometric parameter(s) that is weighted by MD-calculated probability distribution(s) for the geometric parameter(s) (Ad-MD method). The usage of Ad-MD for polymers is computationally convenient when one pre-calculated structural dependence of an NMR parameter is employed for all chemically equivalent units differing only in dynamic behavior. The Ad-MD method is benchmarked against the statistical averaging method for δ_31P in the model phosphates featuring distinctively different structures and dynamic behavior. The applicability of Ad-MD is illustrated by calculating ³¹ P NMR spectra in the Dickerson-Drew DNA dodecamer. δ_31P was calculated with the B3LYP/IGLO-III/PCM(water) and the probability distributions for the torsion angles adjacent to the phosphorus atoms in the DNA phosphates were calculated using the OL15 force field.

Single-Base Lesions and Mismatches Alter the Backbone Conformational Dynamics in DNA

DNA damage has been implicated in numerous human diseases, particularly cancer, and the aging process. Single-base lesions and mismatches in DNA can be cytotoxic or mutagenic and are recognized by a DNA glycosylase during the process of base excision repair. Altered local dynamics and conformational properties in damaged DNAs have previously been suggested to assist in recognition and specificity. Herein, we use solution nuclear magnetic resonance to quantify changes in BI-BII backbone conformational dynamics due to the presence of single-base lesions in DNA, including uracil, dihydrouracil, 1,N6-ethenoadenine, and T:G mismatches. Stepwise changes to the %BII and ΔG of the BI-BII dynamic equilibrium compared to those of unmodified sequences were observed. Additionally, the equilibrium skews toward endothermicity for the phosphates nearest the lesion/mismatched base pair. Finally, the phosphates with the greatest alterations correlate with those most relevant to the repair of enzyme binding. All of these results suggest local conformational rearrangement of the DNA backbone may play a role in lesion recognition by repair enzymes.

Deconvolution of conformational exchange from Raman spectra of aqueous RNA nucleosides

Ribonucleic acids (RNAs) are key to the central dogma of molecular biology. While Raman spectroscopy holds great potential for studying RNA conformational dynamics, current computational Raman prediction and assignment methods are limited in terms of system size and inclusion of conformational exchange. Here, a framework is presented that predicts Raman spectra using mixtures of sub-spectra corresponding to major conformers calculated using classical and ab initio molecular dynamics. Experimental optimization allowed purines and pyrimidines to be characterized as predominantly syn and anti, respectively, and ribose into exchange between equivalent south and north populations. These measurements are in excellent agreement with Raman spectroscopy of ribonucleosides, and previous experimental and computational results. This framework provides a measure of ribonucleoside solution populations and conformational exchange in RNA subunits. It complements other experimental techniques and could be extended to other molecules, such as proteins and carbohydrates, enabling biological insights and providing a new analytical tool.

Structural interpretation of the 31P NMR chemical shifts in thiophosphate and phosphate: key effects due to spin–orbit and explicit solvent

Structural interpretation of the ³¹P NMR shifts measured in different molecules including thiophosphate or phosphate group was obtained by means of theoretical calculations including the effects of geometry, molecular dynamics, solvent, relativistic effects and the effect of NMR reference.

The benchmark of 31P NMR parameters in phosphate: a case study on structurally constrained and flexible phosphate

A benchmark for structural interpretation of the ³¹P NMR shift and the ²J_P,C spin–spin coupling in the phosphate group was obtained by means of theoretical calculations and measurements in diethylphosphate and 5,5-dimethyl-2-hydroxy-1,3,2-dioxaphosphinane 2-oxide.

The influence of Mg2+ coordination on 13C and 15N chemical shifts in CKI1RD protein domain from experiment and molecular dynamics/density functional theory calculations

Sequence dependence of (13) C and (15) N chemical shifts in the receiver domain of CKI1 protein from Arabidopsis thaliana, CKI1RD , and its complexed form, CKI1RD •Mg(2+), was studied by means of MD/DFT calculations. MD simulations of a 20-ns production run length were performed. Nine explicitly hydrated structures of increasing complexity were explored, up to a 40-amino-acid structure. The size of the model necessary depended on the type of nucleus, the type of amino acid and its sequence neighbors, other spatially close amino acids, and the orientation of amino acid NH groups and their surface/interior position. Using models covering a 10 and a 15 Å environment of Mg(2+), a semi-quantitative agreement has been obtained between experiment and theory for the V67-I73 sequence. The influence of Mg(2+) binding was described better by the 15 Å as compared to the 10 Å model. Thirteen chemical shifts were analyzed in terms of the effect of Mg(2+) insertion and geometry preparation. The effect of geometry was significant and opposite in sign to the effect of Mg(2+) binding. The strongest individual effects were found for (15) N of D70, S74, and V68, where the electrostatics dominated; for (13) Cβ of D69 and (15) N of K76, where the influences were equal, and for (13) Cα of F72 and (13) Cβ of K76, where the geometry adjustment dominated. A partial correlation between dominant geometry influence and torsion angle shifts upon the coordination has been observed.

The intrinsic mechanics of B-DNA in solution characterized by NMR

Experimental characterization of the structural couplings in free B-DNA in solution has been elusive, because of subtle effects that are challenging to tackle. Here, the exploitation of the NMR measurements collected on four dodecamers containing a substantial set of dinucleotide sequences provides new, consistent correlations revealing the DNA intrinsic mechanics. The difference between two successive residual dipolar couplings (ΔRDCs) involving C6/8-H6/8, C3′-H3′ and C4′-H4′ vectors are correlated to the ³¹P chemical shifts (δP), which reflect the populations of the BI and BII backbone states. The δPs are also correlated to the internucleotide distances (D_inter) involving H6/8, H2′ and H2″ protons. Calculations of NMR quantities on high resolution X-ray structures and controlled models of DNA enable to interpret these couplings: the studied ΔRDCs depend mostly on roll, while D_inter are mainly sensitive to twist or slide. Overall, these relations demonstrate how δP measurements inform on key inter base parameters, in addition to probe the BI↔BII backbone equilibrium, and shed new light into coordinated motions of phosphate groups and bases in free B-DNA in solution. Inspection of the 5′ and 3′ ends of the dodecamers also supplies new information on the fraying events, otherwise neglected.

Simulations Meet Experiment to Reveal New Insights into DNA Intrinsic Mechanics

The accurate prediction of the structure and dynamics of DNA remains a major challenge in computational biology due to the dearth of precise experimental information on DNA free in solution and limitations in the DNA force-fields underpinning the simulations. A new generation of force-fields has been developed to better represent the sequence-dependent B-DNA intrinsic mechanics, in particular with respect to the BI ↔ BII backbone equilibrium, which is essential to understand the B-DNA properties. Here, the performance of MD simulations with the newly updated force-fields Parmbsc0εζOLI and CHARMM36 was tested against a large ensemble of recent NMR data collected on four DNA dodecamers involved in nucleosome positioning. We find impressive progress towards a coherent, realistic representation of B-DNA in solution, despite residual shortcomings. This improved representation allows new and deeper interpretation of the experimental observables, including regarding the behavior of facing phosphate groups in complementary dinucleotides, and their modulation by the sequence. It also provides the opportunity to extensively revisit and refine the coupling between backbone states and inter base pair parameters, which emerges as a common theme across all the complementary dinucleotides. In sum, the global agreement between simulations and experiment reveals new aspects of intrinsic DNA mechanics, a key component of DNA-protein recognition.

Refinement of the Sugar-Phosphate Backbone Torsion Beta for AMBER Force Fields Improves the Description of Z- and B-DNA

Z-DNA duplexes are a particularly complicated test case for current force fields. We performed a set of explicit solvent molecular dynamics (MD) simulations with various AMBER force field parametrizations including our recent refinements of the ε/ζ and glycosidic torsions. None of these force fields described the ZI/ZII and other backbone substates correctly, and all of them underpredicted the population of the important ZI substate. We show that this underprediction can be attributed to an inaccurate potential for the sugar-phosphate backbone torsion angle β. We suggest a refinement of this potential, β(OL1), which was derived using our recently introduced methodology that includes conformation-dependent solvation effects. The new potential significantly increases the stability of the dominant ZI backbone substate and improves the overall description of the Z-DNA backbone. It also has a positive (albeit small) impact on another important DNA form, the antiparallel guanine quadruplex (G-DNA), and improves the description of the canonical B-DNA backbone by increasing the population of BII backbone substates, providing a better agreement with experiment. We recommend using β(OL1) in combination with our previously introduced corrections, εζ(OL1) and χ(OL4), (the combination being named OL15) as a possible alternative to the current β torsion potential for more accurate modeling of nucleic acids.

Searching the Sequence Space for Potent Aptamers Using SELEX in Silico

To isolate functional nucleic acids that bind to defined targets with high affinity and specificity, which are known as aptamers, the systematic evolution of ligands by exponential enrichment (SELEX) methodology has emerged as the preferred approach. Here, we propose a computational approach, SELEX in silico, that allows the sequence space to be more thoroughly explored regarding binding of a certain target. Our approach consists of two steps: (i) secondary structure-based sequence screening, which aims to collect the sequences that can form a desired RNA motif as an enhanced initial library, followed by (ii) sequence enrichment regarding target binding by molecular dynamics simulation-based virtual screening. Our SELEX in silico method provided a practical computational solution to three key problems in aptamer sequence searching: design of nucleic acid libraries, knowledge of sequence enrichment, and identification of potent aptamers. Six potent theophylline-binding aptamers, which were isolated by SELEX in silico from a sequence space containing 4(13) sequences, were experimentally verified to bind theophylline with high affinity: Kd ranging from 0.16 to 0.52 μM, compared with the dissociation constant of the original aptamer-theophylline, 0.32 μM. These results demonstrate the significant potential of SELEX in silico as a new method for aptamer discovery and optimization.

Vibrational Properties of the Phosphate Group Investigated by Molecular Dynamics and Density Functional Theory

The phosphate group (PO2(-)) is an important building block occurring in many components of living matter including nucleic acids. It provides distinct features in vibrational spectra and is useful as a local probe of NA conformation and interactions with the environment. For this purpose, it is desirable to explore in detail various factors influencing spectral shapes of characteristic phosphate vibrations. In the present study, effects of the solvent and conformational averaging are analyzed for simple model molecules, dimethylphosphate, ethylmethylphosphate, and ethylmethylthiophosphate. Infrared absorption (IR) and Raman spectra were measured and calculated using a combination of molecular dynamics (MD) and density functional theory (DFT). To fully understand the link between the structure and the spectra, the solvent has to be explicitly included in the computational modeling. The results indicate that vibrational properties of the phosphate moiety are very sensitive to its conformation and interactions with the aqueous environment indeed. Polarizable continuum solvent models without explicit water molecules provided significantly worse agreement with the experiment. The combined MD/DFT approach captures well spectral characteristics for the model systems and constitutes the most reliable basis for exploration of phosphate vibrational properties in biomolecular structural studies.

Dynamic hydrogen bonding and DNA flexibility in minor groove binders: molecular dynamics simulation of the polyamide f-ImPyIm bound to the Mlu1 (MCB) sequence 5'-ACGCGT-3' in 2:1 motif

Molecular dynamics simulations of the DNA 10-mer 5'-CCACGCGTGG-3' alone and complexed with the formamido-imidazole-pyrrole-imidazole (f-ImPyIm) polyamide minor groove binder in a 2:1 fashion were conducted for 50 ns using the pbsc0 parameters within the AMBER 12 software package. The change in DNA structure upon binding of f-ImPyIm was evaluated via minor groove width and depth, base pair parameters of Slide, Twist, Roll, Stretch, Stagger, Opening, Propeller, and x-displacement, dihedral angle distributions of ζ, ε, α, and γ determined using the Curves+ software program, and hydrogen bond formation. The dynamic hydrogen bonding between the f-ImPyIm and its cognate DNA sequence was compared to the static image used to predict sequence recognition by polyamide minor groove binders. Many of the predicted hydrogen bonds were present in less than 50% of the simulation; however, persistent hydrogen bonds between G5/15 and the formamido group of f-ImPyIm were observed. It was determined that the DNA is wider in the Complex than without the polyamide binder; however, there is flexibility in this particular sequence, even in the presence of the f-ImPyIm as evidenced by the range of minor groove widths the DNA exhibits and the dynamics of the hydrogen bonding that binds the two f-ImPyIm ions to the minor groove. The Complex consisting of the DNA and the 2 f-ImPyIm binders shows slight fraying of the 5' end of the 10-mer at the end of the simulation, but the portion of the oligomer responsible for recognition and binding is stable throughout the simulation. Several structural changes in the Complex indicate that minor groove binders may have a more active role in inhibiting transcription than just preventing binding of important transcription factors.

Solid-state NMR studies of nucleic acid components

Recent applications of solid-state NMR spectroscopy to studies of nucleic acids and their components.

Accurate ab initio prediction of NMR chemical shifts of nucleic acids and nucleic acids/protein complexes

NMR chemical shift predictions based on empirical methods are nowadays indispensable tools during resonance assignment and 3D structure calculation of proteins. However, owing to the very limited statistical data basis, such methods are still in their infancy in the field of nucleic acids, especially when non-canonical structures and nucleic acid complexes are considered. Here, we present an ab initio approach for predicting proton chemical shifts of arbitrary nucleic acid structures based on state-of-the-art fragment-based quantum chemical calculations. We tested our prediction method on a diverse set of nucleic acid structures including double-stranded DNA, hairpins, DNA/protein complexes and chemically-modified DNA. Overall, our quantum chemical calculations yield highly/very accurate predictions with mean absolute deviations of 0.3–0.6 ppm and correlation coefficients (r²) usually above 0.9. This will allow for identifying misassignments and validating 3D structures. Furthermore, our calculations reveal that chemical shifts of protons involved in hydrogen bonding are predicted significantly less accurately. This is in part caused by insufficient inclusion of solvation effects. However, it also points toward shortcomings of current force fields used for structure determination of nucleic acids. Our quantum chemical calculations could therefore provide input for force field optimization.

NMR studies of DNA support the role of pre-existing minor groove variations in nucleosome indirect readout

We investigated how the intrinsic sequence-dependent properties probed via the phosphate linkages (BI ↔ BII equilibrium) influence the preferred shape of free DNA, and how this affects the nucleosome formation. First, this exploits NMR solution studies of four B-DNA dodecamers that together cover 39 base pairs of the 5' half of the sequence 601, of special interest for nucleosome formation. The results validate our previous prediction of a systematic, general sequence effect on the intrinsic backbone BII propensities. NMR provides new evidence that the backbone behavior is intimately coupled to the minor groove width. Second, application of the backbone behavior predictions to the full sequence 601 and other relevant sequences demonstrates that alternation of intrinsic low and high BII propensities, coupled to intrinsic narrow and wide minor grooves, largely coincides with the sinusoidal variations of the DNA minor groove width observed in crystallographic structures of the nucleosome. This correspondence is much poorer with low affinity sequences. Overall, the results indicate that nucleosome formation involves an indirect readout process implicating pre-existing DNA minor groove conformations. It also illustrates how the prediction of the intrinsic structural DNA behavior offers a powerful framework to gain explanatory insight on how proteins read DNA.

Conformational Sampling by Ab Initio Molecular Dynamics Simulations Improves NMR Chemical Shift Predictions

Car-Parrinello molecular dynamics simulations were performed for N-methyl acetamide as a small test system for amide groups in protein backbones, and NMR chemical shifts were calculated based on the generated ensemble. If conformational sampling and explicit solvent molecules are taken into account, excellent agreement between the calculated and experimental chemical shifts is obtained. These results represent a landmark improvement over calculations based on classical molecular dynamics (MD) simulations especially for amide protons, which are predicted too high-field shifted based on the latter ensembles. We were able to show that the better results are caused by the solute-solvents interactions forming shorter hydrogen bonds as well as by the internal degrees of freedom of the solute. Inspired by these results, we propose our approach as a new tool for the validation of force fields due to its power of identifying the structural reasons for discrepancies between the experimental and calculated data.

Toward Improved Description of DNA Backbone: Revisiting Epsilon and Zeta Torsion Force Field Parameters

We present a refinement of the backbone torsion parameters ε and ζ of the Cornell et al. AMBER force field for DNA simulations. The new parameters, denoted as εζOL1, were derived from quantum-mechanical calculations with inclusion of conformation-dependent solvation effects according to the recently reported methodology (J. Chem. Theory Comput. 2012, 7(9), 2886-2902). The performance of the refined parameters was analyzed by means of extended molecular dynamics (MD) simulations for several representative systems. The results showed that the εζOL1 refinement improves the backbone description of B-DNA double helices and G-DNA stem. In B-DNA simulations, we observed an average increase of the helical twist and narrowing of the major groove, thus achieving better agreement with X-ray and solution NMR data. The balance between populations of BI and BII backbone substates was shifted towards the BII state, in better agreement with ensemble-refined solution experimental results. Furthermore, the refined parameters decreased the backbone RMS deviations in B-DNA MD simulations. In the antiparallel guanine quadruplex (G-DNA) the εζOL1 modification improved the description of non-canonical α/γ backbone substates, which were shown to be coupled to the ε/ζ torsion potential. Thus, the refinement is suggested as a possible alternative to the current ε/ζ torsion potential, which may enable more accurate modeling of nucleic acids. However, long-term testing is recommended before its routine application in DNA simulations.