Balancing the interactions of ions, water, and DNA in the Drude polarizable force field

Balancing Group I Monatomic Ion-Polar Compound Interactions for Condensed Phase Simulation in the Polarizable Drude Force Field

Molecular dynamics (MD) simulations are a commonly used method for investigating molecular behavior at the atomic level. Achieving reliable MD simulation results necessitates the use of an accurate force field. In the present work, we present a protocol to enhance the quality of group 1 monatomic ions (specifically Li+, Na+, K+, Rb+, and Cs+) with respect to their interactions with common polar model compounds in biomolecules in condensed phases in the context of the Drude polarizable force field. Instead of adjusting preexisting individual parameters for ions, model compounds, and water, we employ atom-pair specific Lennard-Jones (LJ) (known as NBFIX in CHARMM) and through-space Thole dipole screening (NBTHOLE) terms to fine-tune the balance of ion-model compound, ion-water, and model compound-water interactions. This involved establishing a protocol for the optimization of NBFIX and NBTHOLE parameters targeting the difference between molecular mechanical (MM) and quantum mechanical (QM) potential energy scans (PES). It is shown that targeting PES involving complexes that include multiple model compounds and/or ions as trimers and tetramers yields parameters that produce condensed phase properties in agreement with experimental data. Validation of this protocol involved the reproduction of experimental thermodynamic benchmarks, including solvation free energies of ions in methanol and N-methylacetamide, osmotic pressures, ionic conductivities, and diffusion coefficients within the condensed phase. These results show the importance of including more complex ion-model compound complexes beyond dimers in the QM target data to account for many-body effects during parameter fitting. The presented parameters represent a significant refinement of the Drude polarizable force field, which will lead to improved accuracy for modeling ion-biomolecular interactions.

sAMP-VGG16: Force-field assisted image-based deep neural network prediction model for short antimicrobial peptides

During the last three decades, antimicrobial peptides (AMPs) have emerged as a promising therapeutic alternative to antibiotics. The approaches for designing AMPs span from experimental trial-and-error methods to synthetic hybrid peptide libraries. To overcome the exceedingly expensive and time-consuming process of designing effective AMPs, many computational and machine-learning tools for AMP prediction have been recently developed. In general, to encode the peptide sequences, featurization relies on approaches based on (a) amino acid (AA) composition, (b) physicochemical properties, (c) sequence similarity, and (d) structural properties. In this work, we present an image-based deep neural network model to predict AMPs, where we are using feature encoding based on Drude polarizable force-field atom types, which can capture the peptide properties more efficiently compared to conventional feature vectors. The proposed prediction model identifies short AMPs (≤30 AA) with promising accuracy and efficiency and can be used as a next-generation screening method for predicting new AMPs. The source code is publicly available at the Figshare server sAMP-VGG16.

A Fast, Convenient, Polarizable Electrostatic Model for Molecular Dynamics

We present an efficient polarizable electrostatic model, utilizing typed, atom-centered polarizabilities and the fast direct approximation, designed for efficient use in molecular dynamics (MD) simulations. The model provides two convenient approaches for assigning partial charges in the context of atomic polarizabilities. One is a generalization of RESP, called RESP-dPol, and the other, AM1-BCC-dPol, is an adaptation of the widely used AM1-BCC method. Both are designed to accurately replicate gas-phase quantum mechanical electrostatic potentials. Benchmarks of this polarizable electrostatic model against gas-phase dipole moments, molecular polarizabilities, bulk liquid densities, and static dielectric constants of organic liquids show good agreement with the reference values. Of note, the model yields markedly more accurate dielectric constants of organic liquids, relative to a matched nonpolarizable force field. MD simulations with this method, which is currently parametrized for molecules containing elements C, N, O, and H, run only about 3.6-fold slower than fixed charge force fields, while simulations with the self-consistent mutual polarization average 4.5-fold slower. Our results suggest that RESP-dPol and AM1-BCC-dPol afford improved accuracy relative to fixed charge force fields and are good starting points for developing general, affordable, and transferable polarizable force fields. The software implementing these approaches has been designed to utilize the force field fitting frameworks developed and maintained by the Open Force Field Initiative, setting the stage for further exploration of this approach to polarizable force field development.

Polarizable AMOEBA Model for Simulating Mg2+·Protein·Nucleotide Complexes

Molecular mechanics (MM) simulations have the potential to provide detailed insights into the mechanisms of enzymes that utilize nucleotides as cofactors. In most cases, the activities of these enzymes also require the binding of divalent cations to catalytic sites. However, modeling divalent cations in MM simulations has been challenging. The inclusion of explicit polarization was considered promising, but despite improvements over nonpolarizable force fields and despite the inclusion of "Nonbonded-fix (NB-fix)" corrections, errors in interaction energies of divalent cations with proteins remain large. Importantly, the application of these models fails to reproduce the experimental structural data on Mg2+·Protein·ATP complexes. Focusing on these complexes, here we provide a systematic assessment of the polarizable AMOEBA model and recommend critical changes that substantially improve its predictive performance. Our key results are as follows. We first show that our recent revision of the AMOEBA protein model (AMOEBABIO18-HFC), which contains high field corrections (HFCs) to induced dipoles, dramatically improves Mg2+-protein interaction energies, reducing the mean absolute error (MAE) from 17 to 10 kcal/mol. This further supports the general applicability of AMOEBABIO18-HFC. The inclusion of many-body NB-fix corrections further reduces MAE to 6 kcal/mol, which amounts to less than 2% error. The errors are estimated with respect to vdW-inclusive density functional theory that we benchmark against CCSD(T) calculations and experiments. We also present a new model of ATP with revised polarization parameters to better capture its high field response, as well as new vdW and dihedral parameters. The ATP model accurately predicts experimental Mg2+-ATP binding free energy in the aqueous phase and provides new insights into how Mg2+ associates with ATP. Finally, we show that molecular dynamics (MD) simulations of Mg2+·Kinase·ATP complexes carried out with these improvements lead to a better agreement in global and local catalytic site structures between MD and X-ray crystallography.

NICE-FF: A non-empirical, intermolecular, consistent, and extensible force field for nucleic acids and beyond

A new non-empirical ab initio intermolecular force field (NICE-FF in buffered 14-7 potential form) has been developed for nucleic acids and beyond based on the dimer interaction energies (IEs) calculated at the spin component scaled-MI-second order Møller-Plesset perturbation theory. A fully automatic framework has been implemented for this purpose, capable of generating well-polished computational grids, performing the necessary ab initio calculations, conducting machine learning (ML) assisted force field (FF) parametrization, and extending existing FF parameters by incorporating new atom types. For the ML-assisted parametrization of NICE-FF, interaction energies of ∼18 000 dimer geometries (with IE < 0) were used, and the best fit gave a mean square deviation of about 0.46 kcal/mol. During this parametrization, atom types apparent in four deoxyribonucleic acid (DNA) bases have been first trained using the generated DNA base datasets. Both uracil and hypoxanthine, which contain the same atom types found in DNA bases, have been considered as test molecules. Three new atom types have been added to the DNA atom types by using IE datasets of both pyrazinamide and 9-methylhypoxanthine. Finally, the last test molecule, theophylline, has been selected, which contains already-fitted atom-type parameters. The performance of NICE-FF has been investigated on the S22 dataset, and it has been found that NICE-FF outperforms the well-known FFs by generating the most consistent IEs with the high-level ab initio ones. Moreover, NICE-FF has been integrated into our in-house developed crystal structure prediction (CSP) tool [called FFCASP (Fast and Flexible CrystAl Structure Predictor)], aiming to find the experimental crystal structures of all considered molecules. CSPs, which were performed up to 4 formula units (Z), resulted in NICE-FF being able to locate almost all the known experimental crystal structures with sufficiently low RMSD20 values to provide good starting points for density functional theory optimizations.

Inclusion of High-Field Target Data in AMOEBA's Calibration Improves Predictions of Protein-Ion Interactions

The reliability of molecular mechanics simulations to predict effects of ion binding to proteins depends on their ability to simultaneously describe ion-protein, ion-water, and protein-water interactions. Force fields (FFs) to describe protein-water and ion-water interactions have been constructed carefully and have also been refined routinely to improve accuracy. Descriptions for ion-protein interactions have also been refined, although in an a posteriori manner through the use of "nonbonded-fix (NB-fix)" approaches in which parameters from default Lennard-Jones mixing rules are replaced with those optimized against some reference data. However, even after NB-fix corrections, there remains a significant need for improvement. This is also true for polarizable FFs that include self-consistent inducible moments. Our recent studies on the polarizable AMOEBA FF suggested that the problem associated with modeling ion-protein interactions could be alleviated by recalibrating polarization models of cation-coordinating functional groups so that they respond better to the high electric fields present near ions. Here, we present such a recalibration of carbonyls, carboxylates, and hydroxyls in the AMOEBA protein FF and report that it does improve predictions substantially─mean absolute errors in Na⁺-protein and K⁺-protein interaction energies decrease from 8.7 to 5.3 and 9.6 to 6.3 kcal/mol, respectively. Errors are computed with respect to estimates from van der Waals-inclusive density functional theory benchmarked against high-level quantum mechanical calculations and experiments. While recalibration does improve ion-protein interaction energies, they still remain underestimated, suggesting that further improvements can be made in a systematic manner through modifications in classical formalism. Nevertheless, we show that by applying our many-body NB-fix correction to Lennard-Jones components, these errors are further reduced to 2.7 and 2.6 kcal/mol, respectively, for Na⁺ and K⁺ ions. Finally, we show that the recalibrated AMOEBA protein FF retains its intrinsic reliability in predicting protein structure and dynamics in the condensed phase.

Molecular Dynamics Simulations of Protein RNA Complexes by Using an Advanced Electrostatic Model

Protein-RNA interactions are integral to the biological functions of RNA. It is well recognized that molecular dynamics (MD) simulations of protein-RNA complexes are more challenging than those of each component. The difficulty arises from the strong electrostatic interactions and the delicate balance between various types of physical forces at the interface. Previously, MD simulations of protein-RNA complexes have predominantly employed fixed-charge force fields. Although force field modifications have been developed to address problems identified in the simulations, some protein-RNA structures are still hard to reproduce by simulations. Here, we present MD simulations of two representative protein-RNA complexes using the AMOEBA polarizable force field. The van der Waals parameters were refined to reproduce accurate quantum-mechanical data of base-base and base-amino acid interactions. It was found that the refined parameters produced a more stable hydrogen-bond network in the interface. One of the complexes remained stable during the short simulations, whereas it could quickly break down in previous simulations using fixed-charge force fields. There was reversible breaking and formation of hydrogen bonds that are observed in the crystal structure, which may indicate the difference in solution and crystal structures. While further improvement and validation of the force fields are still needed, this work demonstrates that polarizable force fields are promising for the study of protein-RNA complexes.

Ion-Dependent Conformational Plasticity of Telomeric G-Hairpins and G-Quadruplexes

Telomeric DNA is guanine-rich and can adopt structures such as G-quadruplexes (GQs) and G-hairpins. Telomeric GQs influence genome stability and telomerase activity, making understanding of enzyme-GQ interactions and dynamics important for potential drug design. GQs have a characteristic tetrad core, which is connected by loop regions. Within this architecture are G-hairpins, fold-back motifs that are thought to represent the first intermediate in GQ folding. To better understand the relationship between G-hairpin motifs and GQs, we performed polarizable simulations of a two-tetrad telomeric GQ and an isolated SC11 telomeric G-hairpin. The telomeric GQ contains a G-triad, which functions as part of the tetrad core or linker regions, depending on local conformational change. This triad and another motif below the tetrad core frequently bound ions and may represent druggable sites. Further, we observed the unbiased formation of a G-triad and a G-tetrad in simulations of the SC11 G-hairpin and found that cations can be partially hydrated while facilitating the formation of these motifs. Finally, we demonstrated that K⁺ ions form specific interactions with guanine bases, while Na⁺ ions interact nonspecifically with bases in the structure. Together, these simulations provide new insights into the influence of ions on GQs, G-hairpins, and G-triad motifs.

Harnessing Deep Learning for Optimization of Lennard-Jones Parameters for the Polarizable Classical Drude Oscillator Force Field

The outcomes of computational chemistry and biology research, including drug design, are significantly influenced by the underlying force field (FF) used in molecular simulations. While improved FF accuracy may be achieved via inclusion of explicit treatment of electronic polarization, such an extension must be accompanied by optimization of van der Waals (vdW) interactions, in the context of the Lennard-Jones (LJ) formalism in the present study. This is particularly challenging due to the extensive nature of chemical space combined with the correlated nature of LJ parameters. To address this challenge, a deep learning (DL)-based parametrization framework is developed, allowing for sampling of wide ranges of LJ parameters targeting experimental condensed phase thermodynamic properties. The present work utilizes this framework to develop the LJ parameters for atoms associated with four distinct groups covering 10 different atom types. Final parameter selection was facilitated by quantum mechanical data on rare-gas interactions with the training set molecules. The chosen parameters were then validated through experimental hydration free energies and condensed phase thermodynamic properties of validation set molecules to confirm transferability. The ultimate outcome of utilizing this framework is a set of LJ parameters in the context of the polarizable Drude FF, which demonstrated improvement in the reproduction of both experimental pure solvent and crystal properties and hydration free energies of the molecules compared to the additive CHARMM General FF (CGenFF) including the ability of the Drude FF to accurately reproduce both experimental pure solvent properties and hydration free energies. The study also shows how correlations between difference in the reproduction of condensed phase data between model compounds may be used to direct the selection of new atom types and training set molecules during FF development.

High-Dimensional Parameter Search Method to Determine Force Field Mixing Terms in Molecular Simulations

Molecular dynamics (MD) force fields for lipids and ions are typically developed independently of one another. In simulations consisting of both lipids and ions, lipid-ion interaction energies are estimated using a predefined set of mixing rules for Lennard-Jones (LJ) interactions. This, however, does not guarantee their reliability. In fact, compared to the quantum mechanical reference data, Lorentz-Berthelot mixing rules substantially underestimate the binding energies of Na⁺ ions with small-molecule analogues of lipid headgroups, yielding errors on the order of 80 and 130 kJ/mol, respectively, for methyl acetate and diethyl phosphate. Previously, errors associated with mixing force fields have been reduced using approaches such as "NB-fix" in which LJ interactions are computed using explicit cross terms rather than those from mixing rules. Building on this idea, we derive explicit lipid-ion cross terms that also may implicitly include many-body cooperativity effects. Additionally, to account for the interdependency between cross terms, we optimize all cross terms simultaneously by performing high-dimensional searches using our ParOpt software. The cross terms we obtain reduce the errors due to mixing rules to below 10 kJ/mol. MD simulation of the lipid bilayer conducted using these optimized cross terms resolves the structural discrepancies between our previous simulations and small-angle X-ray and neutron scattering experiments. These results demonstrate that simulations of lipid bilayers with ions that are accurate up to structural data from scattering experiments can be performed without explicit polarization terms. However, it is worth noting that such NB-fix cross terms are not based on any physical principle; a polarizable lipid model would be more realistic and is still desired. Our approach is generic and can be applied to improve the accuracies of simulations employing mixed force fields.

Computational and Experimental Characterization of rDNA and rRNA G-Quadruplexes

DNA G-quadruplexes in human telomeres and gene promoters are being extensively studied for their role in controlling the growth of cancer cells. G-quadruplexes have been unambiguously shown to exist both in vitro and in vivo, including in the guanine (G)-rich DNA genes encoding pre-ribosomal RNA (pre-rRNA), which is transcribed in the cell's nucleolus. Recent studies strongly suggest that these DNA sequences ("rDNA"), and the transcribed rRNA, are a potential anticancer target through the inhibition of RNA polymerase I (Pol I) in ribosome biogenesis, but the structures of ribosomal G-quadruplexes at atomic resolution are unknown and very little biophysical characterization has been performed on them to date. In the present study, circular dichroism (CD) spectroscopy is used to show that two putative rDNA G-quadruplex sequences, NUC 19P and NUC 23P and their counterpart rRNAs, predominantly adopt parallel topologies, reminiscent of the analogous telomeric quadruplex structures. Based on this information, we modeled parallel topology atomistic structures of the putative ribosomal G-quadruplexes. We then validated and refined the modeled ribosomal G-quadruplex structures using all-atom molecular dynamics (MD) simulations with the CHARMM36 force field in the presence and absence of stabilizing K⁺. Motivated by preliminary MD simulations of the telomeric parallel G-quadruplex (TEL 24P) in which the K⁺ ion is expelled, we used updated CHARMM36 force field K⁺ parameters that were optimized, targeting the data from quantum mechanical calculations and the polarizable Drude model force field. In subsequent MD simulations with optimized CHARMM36 parameters, the K⁺ ions are predominantly in the G-quadruplex channel and the rDNA G-quadruplexes have more well-defined, predominantly parallel-topology structures as compared to rRNA. In addition, NUC 19P is more structured than NUC 23P, which contains extended loops. Results from this study set the structural foundation for understanding G-quadruplex functions and the design of novel chemotherapeutics against these nucleolar targets and can be readily extended to other DNA and RNA G-quadruplexes.

CHARMM-GUI Drude prepper for molecular dynamics simulation using the classical Drude polarizable force field

Explicit treatment of electronic polarizability in empirical force fields (FFs) represents an extension over a traditional additive or pairwise FF and provides a more realistic model of the variations in electronic structure in condensed phase, macromolecular simulations. To facilitate utilization of the polarizable FF based on the classical Drude oscillator model, Drude Prepper has been developed in CHARMM-GUI. Drude Prepper ingests additive CHARMM protein structures file (PSF) and pre-equilibrated coordinates in CHARMM, PDB, or NAMD format, from which the molecular components of the system are identified. These include all residues and patches connecting those residues along with water, ions, and other solute molecules. This information is then used to construct the Drude FF-based PSF using molecular generation capabilities in CHARMM, followed by minimization and equilibration. In addition, inputs are generated for molecular dynamics (MD) simulations using CHARMM, GROMACS, NAMD, and OpenMM. Validation of the Drude Prepper protocol and inputs is performed through conversion and MD simulations of various heterogeneous systems that include proteins, nucleic acids, lipids, polysaccharides, and atomic ions using the aforementioned simulation packages. Stable simulations are obtained in all studied systems, including 5 μs simulation of ubiquitin, verifying the integrity of the generated Drude PSFs. In addition, the ability of the Drude FF to model variations in electronic structure is shown through dipole moment analysis in selected systems. The capabilities and availability of Drude Prepper in CHARMM-GUI is anticipated to greatly facilitate the application of the Drude FF to a range of condensed phase, macromolecular systems.

Molecular dynamics free energy simulations of ATP:Mg2+ and ADP:Mg2+ using the polarizable force field AMOEBA

ATPases and GTPases are two important classes of protein that play critical roles in energy transduction, cellular signaling, gene regulation and catalysis. These proteins use cofactors such as nucleoside di and tri-phosphates (NTP, NDP) and can detect the difference between NDP and NTP which then induce different protein conformations. Mechanisms that drive proteins into the NTP or NDP conformation may depend on factors such as ligand structure and how Mg²⁺ coordinates with the ligand, amino acids in the pocket and water molecules. Here, we have used the advanced electrostatic and polarizable force field AMOEBA and molecular dynamics free energy simulations (MDFE) to examine the various binding mechanisms of ATP:Mg²⁺ and ADP:Mg²⁺.We compared the ATP:Mg²⁺ binding with previous studies using non-polarizable force fields and experimental data on the binding affinity. It was found that the total free energy of binding for ATP:Mg²⁺ (-7.00 ± 2.13 kcal/mol) is in good agreement with experimental values (-8.6 ± .2 kcal/mol)¹. In addition, parameters for relevant protonation states of ATP, ADP, GTP and GDP have been derived. These parameters will allow for researchers to investigate biochemical phenomena involving NTP's and NDP's with greater accuracy than previous studies involving non-polarizable force fields.

Water-Mediated Oxidation of Guanine by a Repair Enzyme: Simulation Using the ABEEM Polarizable Force Field

The recognition mechanism of oxidative damage in organisms has long been a research hotspot. Water is an important medium in the recognition process, but its specific role remains unknown. There is a need to develop a suitable force field that can adequately describe the electrostatic, hydrogen bond, and other interactions among the molecules in the complex system of the repair enzyme and oxidized base. The developing ABEEM polarizable force field (PFF) has been used to simulate the repaired enzyme hOGG1 and oxidized DNA (PDB ID: 1EBM) in a biological environment, and the corresponding results are better than those of the fixed-charge force fields OPLS/AA and AMBER OL15. 8-Oxo-G is recognized by Gln315 of hOGG1 mainly through hydrogen bonds mediated by continuous exchange of 2 water molecules. Phe319 and Cys253 are stacked on both sides of the π planes of bases to form sandwich structures. The charge polarization effect gives an important signal to drive the exchange of water molecules and maintains the recognition of oxidation bases by enzymes. The mediated main water molecule A and mediated auxiliary water molecule B together pull Gln315 to recognize 8-oxo-G by hydrogen bond interactions. Then, the charge polarization signal of solvent water molecule C with a large absolute charge causes the absolute charge of O atoms in water molecule A or B to increase by approximately 0.2 e, and water molecule A or B leaves Gln315 and 8-oxo-G. The other water molecule and water molecule C synergistically recognize 8-oxo-G with Gln315. Even though the water molecules between Gln315 and 8-oxo-G are removed, the MD simulation results show that water molecules appear between Gln315 and 8-oxo-G in a very short time (<2 ps). The dwell time of each water molecule is approximately 60 ps. The radial distribution function and dwell time support the correctness of the above mechanism. These polarization effects and hydrogen bonding interactions cannot be simulated by a fixed-charge force field.

Assessment of the DNA partial specific volume and hydration layer properties from CHARMM Drude polarizable and additive MD simulations

In this study we report on the accurate computation of the biomolecular partial specific volume (PSV) from explicit-solvent molecular dynamics (MD) simulations. The case of DNA is considered, and the predictions from two state-of-the-art biomolecular force fields, the CHARMM36 additive (C36) and Drude polarizable models, are presented. Unlike most of the existing approaches to assess the biomolecular PSV, our proposed method bypasses the need for the arbitrarily defined volume partitioning scheme into the intrinsic solute and solvent contributions. At the same time, to assess the density of the hydration layer water, we combine our simulation analysis approach with some of the existing fixed-size methods to determine the solute's intrinsic volume, and also propose our own approach to compute all required quantities exclusively from MD simulations. Our findings provide useful insights into the properties of the hydration layer, specifically its size and density, parameters of great importance to the variety of techniques used to model hydrodynamic and structural properties of biological molecules. The computed PSV values are found to be in close agreement with the values obtained from analytical ultracentrifugation (AUC) experiments performed on canonical B-form duplex DNAs and single-stranded DNAs forming G-quadruplex structures. Since the biomolecular PSV represents an important quantitative measure of solute-solvent interactions, near quantitative agreement with AUC measurements is indicative of the quality of the all-atom models used in the MD simulations, particularly the reliability of the CHARMM force-field parameters for nucleic acids, water, mobile ions, and interactions among these entities.

An electrostatic energy-based charge model for molecular dynamics simulation

The interactions of the polar chemical bonds such as C=O and N-H with an external electric field were investigated, and a linear relationship between the QM/MM interaction energies and the electric field along the chemical bond is established in the range of weak to intermediate electrical fields. The linear relationship indicates that the electrostatic interactions of a polar group with its surroundings can be described by a simple model of a dipole with constant moment under the action of an electric field. This relationship is employed to develop a general approach to generating an electrostatic energy-based charge (EEC) model for molecules containing single or multiple polar chemical bonds. Benchmark test studies of this model were carried out for (CH₃)₂-CO and N-methyl acetamide in explicit water, and the result shows that the EEC model gives more accurate electrostatic energies than those given by the widely used charge model based on fitting to the electrostatic potential (ESP) in direct comparison to the energies computed by the QM/MM method. The MD simulations of the electric field at the active site of ketosteroid isomerase based on EEC demonstrated that EEC gave a better representation of the electrostatic interaction in the hydrogen-bonding environment than the Amber14SB force field by comparison with experiment. The current study suggests that EEC should be better suited for molecular dynamics study of molecular systems with polar chemical bonds such as biomolecules than the widely used ESP or RESP (restrained ESP) charge models.

Cation competition and recruitment around the c-kit1 G-quadruplex using polarizable simulations

Nucleic acid-ion interactions are fundamentally important to the physical, energetic, and conformational properties of DNA and RNA. These interactions help fold and stabilize highly ordered secondary and tertiary structures, such as G-quadruplexes (GQs), which are functionally relevant in telomeres, replication initiation sites, and promoter sequences. The c-kit proto-oncogene encodes for a receptor tyrosine kinase and is linked to gastrointestinal stromal tumors, mast cell disease, and leukemia. This gene contains three unique GQ-forming sequences that have proposed antagonistic effects on gene expression. The dominant GQ, denoted c-kit1, has been shown to decrease expression of c-kit transcripts, making the c-kit1 GQ a promising drug target. Toward disease intervention, more information is needed regarding its conformational dynamics and ion binding properties. Therefore, we performed molecular dynamics simulations of the c-kit1 GQ with K⁺, Na⁺, Li⁺, and mixed salt solutions using the Drude-2017 polarizable force field. We evaluated GQ structure, ion sampling, core energetics, ion dehydration and binding, and ion competition and found that each analysis supported the known GQ-ion specificity trend (K⁺ > Na⁺ > Li⁺). We also found that K⁺ ions coordinate in the tetrad core antiprismatically, whereas Na⁺ and Li⁺ align coplanar to guanine tetrads, partially because of their attraction to surrounding water. Further, we showed that K⁺ occupancy is higher around the c-kit1 GQ and its nucleobases than Na⁺ and Li⁺, which tend to interact with backbone and sugar moieties. Finally, we showed that K⁺ binding to the c-kit1 GQ is faster and more frequent than Na⁺ and Li⁺. Such descriptions of GQ-ion dynamics suggest the rate of dehydration as the dominant factor for preference of K⁺ by DNA GQs and provide insight into noncanonical nucleic acids for which little experimental data exist.

Ion Binding Properties and Dynamics of the bcl-2 G-Quadruplex Using a Polarizable Force Field

G-quadruplexes (GQs) are topologically diverse, highly thermostable noncanonical nucleic acid structures that form in guanine-rich sequences in DNA and RNA. GQs are implicated in transcriptional and translational regulation and genome maintenance, and deleterious alterations to their structures contribute to diseases such as cancer. The expression of the B-cell lymphoma 2 (Bcl-2) antiapoptotic protein, for example, is under transcriptional control of a GQ in the promoter of the bcl-2 gene. Modulation of the bcl-2 GQ by small molecules is of interest for chemotherapeutic development but doing so requires knowledge of the factors driving GQ folding and stabilization. To develop a greater understanding of the electrostatic properties of the bcl-2 promoter GQ, we performed molecular dynamics simulations using the Drude-2017 polarizable force field and compared relevant outcomes to the nonpolarizable CHARMM36 force field. Our simulation outcomes highlight the importance of dipole-dipole interactions in the bcl-2 GQ, particularly during the recruitment of a bulk K⁺ ion to the solvent-exposed face of the tetrad stem. We also predict and characterize an "electronegative pocket" at the tetrad-long loop junction that induces local backbone conformational change and may induce local conformational changes at cellular concentrations of K⁺. These outcomes suggest that moieties within the bcl-2 GQ can be targeted by small molecules to modulate bcl-2 GQ stability.

Contributions and competition of Mg2+ and K+ in folding and stabilization of the Twister ribozyme

Native folded and compact intermediate states of RNA typically involve tertiary structures in the presence of divalent ions such as Mg²⁺ in a background of monovalent ions. In a recent study, we have shown how the presence of Mg²⁺ impacts the transition from partially unfolded to folded states through a "push-pull" mechanism where the ion both favors and disfavors the sampling of specific phosphate-phosphate interactions. To further understand the ion atmosphere of RNA in folded and partially folded states results from atomistic umbrella sampling and oscillating chemical potential grand canonical Monte Carlo/molecular dynamics (GCMC/MD) simulations are used to obtain atomic-level details of the distributions of Mg²⁺ and K⁺ ions around Twister RNA. Results show the presence of 100 mM Mg²⁺ to lead to increased charge neutralization over that predicted by counterion condensation theory. Upon going from partially unfolded to folded states, overall charge neutralization increases at all studied ion concentrations that, while associated with an increase in the number of direct ion-phosphate interactions, is fully accounted for by the monovalent K⁺ ions. Furthermore, K⁺ preferentially interacts with purine N7 atoms of helical regions in partially unfolded states, thereby potentially stabilizing the helical regions. Thus, both secondary helical structures and formation of tertiary structures leads to increased counterion condensation, thereby stabilizing those structural features of Twister. Notably, it is shown that K⁺ can act as a surrogate for Mg²⁺ by participating in specific interactions with nonsequential phosphate pairs that occur in the folded state, explaining the ability of Twister to self-cleave at submillimolar Mg²⁺ concentrations.

Recent developments in empirical atomistic force fields for nucleic acids and applications to studies of folding and dynamics

Nucleic acids play critical roles in carrying genetic information, participating in catalysis, and preserving chromosomal structure. Despite over a century of study, efforts to understand the dynamics and structure-function relationships of DNA and RNA at the atomic level are still ongoing. Molecular dynamics (MD) simulations augment experiments by providing atomistic resolution and quantitative relationships between structure and conformational energy. Steady advancements in computer hardware, software, and atomistic force fields (FFs) over 40 years have facilitated new discoveries. Here, we review nucleic acid FF development with emphasis on recent refinements that have improved descriptions of important nucleic acid properties. We then discuss several key examples of successes and challenges in modeling nucleic acid structure and dynamics using the latest FFs.

Transferable interactions of Li+ and Mg2+ ions in polarizable models

Therapeutic implications of Li⁺, in many cases, stem from its ability to inhibit certain Mg²⁺-dependent enzymes, where it interacts with or substitutes for Mg²⁺. The underlying details of its action are, however, unknown. Molecular simulations can provide insights, but their reliability depends on how well they describe relative interactions of Li⁺ and Mg²⁺ with water and other biochemical groups. Here, we explore, benchmark, and recommend improvements to two simulation approaches: the one that employs an all-atom polarizable molecular mechanics (MM) model and the other that uses a hybrid quantum and MM implementation of the quasi-chemical theory (QCT). The strength of the former is that it describes thermal motions explicitly and that of the latter is that it derives local contributions from electron densities. Reference data are taken from the experiment, and also obtained systematically from CCSD(T) theory, followed by a benchmarked vdW-inclusive density functional theory. We find that the QCT model predicts relative hydration energies and structures in agreement with the experiment and without the need for additional parameterization. This implies that accurate descriptions of local interactions are essential. Consistent with this observation, recalibration of local interactions in the MM model, which reduces errors from 10.0 kcal/mol to 1.4 kcal/mol, also fixes aqueous phase properties. Finally, we show that ion-ligand transferability errors in the MM model can be reduced significantly from 10.3 kcal/mol to 1.2 kcal/mol by correcting the ligand's polarization term and by introducing Lennard-Jones cross-terms. In general, this work sets up systematic approaches to evaluate and improve molecular models of ions binding to proteins.

Improved description of ligand polarization enhances transferability of ion-ligand interactions

The reliability of molecular mechanics (MM) simulations in describing biomolecular ion-driven processes depends on their ability to accurately model interactions of ions simultaneously with water and other biochemical groups. In these models, ion descriptors are calibrated against reference data on ion-water interactions, and it is then assumed that these descriptors will also satisfactorily describe interactions of ions with other biochemical ligands. The comparison against the experiment and high-level quantum mechanical data show that this transferability assumption can break down severely. One approach to improve transferability is to assign cross terms or separate sets of non-bonded descriptors for every distinct pair of ion type and its coordinating ligand. Here, we propose an alternative solution that targets an error-source directly and corrects misrepresented physics. In standard model development, ligand descriptors are never calibrated or benchmarked in the high electric fields present near ions. We demonstrate for a representative MM model that when the polarization descriptors of its ligands are improved to respond to both low and high fields, ligand interactions with ions also improve, and transferability errors reduce substantially. In our case, the overall transferability error reduces from 3.3 kcal/mol to 1.8 kcal/mol. These improvements are observed without compromising on the accuracy of low-field interactions of ligands in gas and condensed phases. Reference data for calibration and performance evaluation are taken from the experiment and also obtained systematically from "gold-standard" CCSD(T) in the complete basis set limit, followed by benchmarked vdW-inclusive density functional theory.

Same fold, different properties: polarizable molecular dynamics simulations of telomeric and TERRA G-quadruplexes

DNA and RNA sequences rich in guanine can fold into noncanonical structures called G-quadruplexes (GQs), which exhibit a common stem structure of Hoogsteen hydrogen-bonded guanine tetrads and diverse loop structures. GQ sequence motifs are overrepresented in promoters, origins of replication, telomeres, and untranslated regions in mRNA, suggesting roles in modulating gene expression and preserving genomic integrity. Given these roles and unique aspects of different structures, GQs are attractive targets for drug design, but greater insight into GQ folding pathways and the interactions stabilizing them is required. Here, we performed molecular dynamics simulations to study two bimolecular GQs, a telomeric DNA GQ and the analogous telomeric repeat-containing RNA (TERRA) GQ. We applied the Drude polarizable force field, which we show outperforms the additive CHARMM36 force field in both ion retention and maintenance of the GQ folds. The polarizable simulations reveal that the GQs bind bulk K+ ions differently, and that the TERRA GQ accumulates more K+ ions, suggesting different ion interactions stabilize these structures. Nucleobase dipole moments vary as a function of position and also contribute to ion binding. Finally, we show that the TERRA GQ is more sensitive than the telomeric DNA GQ to water-mediated modulation of ion-induced dipole-dipole interactions.

Molecular dynamics simulations of alkaline earth metal ions binding to DNA reveal ion size and hydration effects

Classical molecular dynamics simulations reveal size-dependent trends of alkaline earth metal ions binding to DNA are due to ion size and hydration behavior.

Molecular Dynamics Simulations of Ionic Liquids and Electrolytes Using Polarizable Force Fields

Many applications in chemistry, biology, and energy storage/conversion research rely on molecular simulations to provide fundamental insight into structural and transport properties of materials with high ionic concentrations. Whether the system is comprised entirely of ions, like ionic liquids, or is a mixture of a polar solvent with a salt, e.g., liquid electrolytes for battery applications, the presence of ions in these materials results in strong local electric fields polarizing solvent molecules and large ions. To predict properties of such systems from molecular simulations often requires either explicit or mean-field inclusion of the influence of polarization on electrostatic interactions. In this manuscript, we review the pros and cons of different treatments of polarization ranging from the mean-field approaches to the most popular explicit polarization models in molecular dynamics simulations of ionic materials. For each method, we discuss their advantages and disadvantages and emphasize key assumptions as well as their adjustable parameters. Strategies for the development of polarizable models are presented with a specific focus on extracting atomic polarizabilities. Finally, we compare simulations using polarizable and nonpolarizable models for several classes of ionic systems, discussing the underlying physics that each approach includes or ignores, implications for implementation and computational efficiency, and the accuracy of properties predicted by these methods compared to experiments.

Ion-Hydroxyl Interactions: From High-Level Quantum Benchmarks to Transferable Polarizable Force Fields

Ion descriptors in molecular mechanics models are calibrated against reference data on ion-water interactions. It is then typically assumed that these descriptors will also satisfactorily describe interactions of ions with other functional groups, such as those present in biomolecules. However, several studies now demonstrate that this transferability assumption produces, in many different cases, large errors. Here we address this issue in a representative polarizable model and focus on transferability of cationic interactions from water to a series of alcohols. Both water and alcohols use hydroxyls for ion-coordination, and, therefore, this set of molecules constitutes the simplest possible case of transferability. We obtain gas phase reference data systematically from "gold-standard" quantum Monte Carlo and CCSD(T) methods, followed by benchmarked vdW-corrected DFT. We learn that the original polarizable model yields large gas phase water → alcohol transferability errors - the RMS and maximum errors are 2.3 and 5.1 kcal/mol, respectively. These errors are, nevertheless, systematic in that ion-alcohol interactions are overstabilized, and systematic errors typically imply that some essential physics is either missing or misrepresented. A comprehensive analysis shows that when both low- and high-field responses of ligand dipole polarization are described accurately, then transferability improves significantly - the RMS and maximum errors in the gas phase reduce, respectively, to 0.9 and 2.5 kcal/mol. Additionally, predictions of condensed phase transfer free energies also improve. Nevertheless, within the limits of the extrathermodynamic assumptions necessary to separate experimental estimates of salt dissolution into constituent cationic and anionic contributions, we note that the error in the condensed phase is systematic, which we attribute, at least, partially to the parametrization in long-range electrostatics. Overall, this work demonstrates a rational approach to boosting transferability of ionic interactions that will be applicable broadly to improving other polarizable and nonpolarizable models.

Force Fields for Small Molecules

Molecular dynamics (MD) simulations have been widely applied to computer-aided drug design (CADD). While MD has been used in a variety of applications such as free energy perturbation and long-time simulations, the accuracy of the results from those methods depends strongly on the force field used. Force fields for small molecules are crucial, as they not only serve as building blocks for developing force fields for larger biomolecules but also act as model compounds that will be transferred to ligands used in CADD. Currently, a wide range of small molecule force fields based on additive or nonpolarizable models have been developed. While these nonpolarizable force fields can produce reasonable estimations of physical properties and have shown success in a variety of systems, there is still room for improvements due to inherent limitations in these models including the lack of an electronic polarization response. For this reason, incorporating polarization effects into the energy function underlying a force field is believed to be an important step forward, giving rise to the development of polarizable force fields. Recent simulations of biological systems have indicated that polarizable force fields are able to provide a better physical representation of intermolecular interactions and, in many cases, better agreement with experimental properties than nonpolarizable, additive force fields. Therefore, this chapter focuses on the development of small molecule force fields with emphasis on polarizable models. It begins with a brief introduction on the importance of small molecule force fields and their evolution from additive to polarizable force fields. Emphasis is placed on the additive CHARMM General Force Field and the polarizable force field based on the classical Drude oscillator. The theory for the Drude polarizable force field and results for small molecules are presented showing their improvements over the additive model. The potential importance of polarization for their application in a wide range of biological systems including CADD is then discussed.

Polarizable force field for RNA based on the classical drude oscillator

RNA molecules are highly dynamic and capable of adopting a wide range of complex, folded structures. The factors driving the folding and dynamics of these structures are dependent on a balance of base pairing, hydration, base stacking, ion interactions, and the conformational sampling of the 2'-hydroxyl group in the ribose sugar. The representation of these features is a challenge for empirical force fields used in molecular dynamics simulations. Toward meeting this challenge, the inclusion of explicit electronic polarization is important in accurately modeling RNA structure. In this work, we present a polarizable force field for RNA based on the classical Drude oscillator model, which represents electronic degrees of freedom via negatively charged particles attached to their parent atoms by harmonic springs. Beginning with parametrization against quantum mechanical base stacking interaction energy and conformational energy data, we have extended the Drude-2017 nucleic acid force field to include RNA. The conformational sampling of a range of RNA sequences were used to validate the force field, including canonical A-form RNA duplexes, stem-loops, and complex tertiary folds that bind multiple Mg2+ ions. Overall, the Drude-2017 RNA force field reproduces important properties of these structures, including the conformational sampling of the 2'-hydroxyl and key interactions with Mg2+ ions. © 2018 Wiley Periodicals, Inc.

Mapping the Drude polarizable force field onto a multipole and induced dipole model

The induced dipole and the classical Drude oscillator represent two major approaches for the explicit inclusion of electronic polarizability into force field-based molecular modeling and simulations. In this work, we explore the equivalency of these two models by comparing condensed phase properties computed using the Drude force field and a multipole and induced dipole (MPID) model. Presented is an approach to map the electrostatic model optimized in the context of the Drude force field onto the MPID model. Condensed phase simulations on water and 15 small model compounds show that without any reparametrization, the MPID model yields properties similar to the Drude force field with both models yielding satisfactory reproduction of a range of experimental values and quantum mechanical data. Our results illustrate that the Drude oscillator model and the point induced dipole model are different representations of essentially the same physical model. However, results indicate the presence of small differences between the use of atomic multipoles and off-center charge sites. Additionally, results on the use of dispersion particle mesh Ewald further support its utility for treating long-range Lennard Jones dispersion contributions in the context of polarizable force fields. The main motivation in demonstrating the transferability of parameters between the Drude and MPID models is that the more than 15 years of development of the Drude polarizable force field can now be used with MPID formalism without the need for dual-thermostat integrators nor self-consistent iterations. This opens up a wide range of new methodological opportunities for polarizable models.

Efficient minimization of multipole electrostatic potentials in torsion space

The development of models of macromolecular electrostatics capable of delivering improved fidelity to quantum mechanical calculations is an active field of research in computational chemistry. Most molecular force field development takes place in the context of models with full Cartesian coordinate degrees of freedom. Nevertheless, a number of macromolecular modeling programs use a reduced set of conformational variables limited to rotatable bonds. Efficient algorithms for minimizing the energies of macromolecular systems with torsional degrees of freedom have been developed with the assumption that all atom-atom interaction potentials are isotropic. We describe novel modifications to address the anisotropy of higher order multipole terms while retaining the efficiency of these approaches. In addition, we present a treatment for obtaining derivatives of atom-centered tensors with respect to torsional degrees of freedom. We apply these results to enable minimization of the Amoeba multipole electrostatics potential in a system with torsional degrees of freedom, and validate the correctness of the gradients by comparison to finite difference approximations. In the interest of enabling a complete model of electrostatics with implicit treatment of solvent-mediated effects, we also derive expressions for the derivative of solvent accessible surface area with respect to torsional degrees of freedom.

AMOEBA Polarizable Atomic Multipole Force Field for Nucleic Acids

The AMOEBA polarizable atomic multipole force field for nucleic acids is presented. Valence and electrostatic parameters were determined from high-level quantum mechanical data, including structures, conformational energy, and electrostatic potentials, of nucleotide model compounds. Previously derived parameters for the phosphate group and nucleobases were incorporated. A total of over 35 μs of condensed-phase molecular dynamics simulations of DNA and RNA molecules in aqueous solution and crystal lattice were performed to validate and refine the force field. The solution and/or crystal structures of DNA B-form duplexes, RNA duplexes, and hairpins were captured with an average root-mean-squared deviation from NMR structures below or around 2.0 Å. Structural details, such as base pairing and stacking, sugar puckering, backbone and χ-torsion angles, groove geometries, and crystal packing interfaces, agreed well with NMR and/or X-ray. The interconversion between A- and B-form DNAs was observed in ethanol-water mixtures at 328 K. Crystal lattices of B- and Z-form DNA and A-form RNA were examined with simulations. For the RNA tetraloop, single strand tetramers, and HIV TAR with 29 residues, the simulated conformational states, 3 J-coupling, nuclear Overhauser effect, and residual dipolar coupling data were compared with NMR results. Starting from a totally unstacked/unfolding state, the rCAAU tetranucleotide was folded into A-form-like structures during ∼1 μs molecular dynamics simulations.

Polarizable Empirical Force Field for Halogen-Containing Compounds Based on the Classical Drude Oscillator

The quality of the force field is crucial to ensure the accuracy of simulations used in molecular modeling, including computer-aided drug design (CADD). To perform more accurate modeling and simulations of halogenated molecules, in this study the polarizable force field based on the classical Drude oscillator model was extended to both aliphatic and aromatic systems using halogenated ethane and benzene model compounds for the halogens F, Cl, Br, and I. The force field parameters were optimized targeting quantum mechanical dipole moments, water interactions, and molecular polarizabilities as well as experimental observables, including enthalpies of vaporization, molecular volumes, hydration free energies, and dielectric constants. The developed halogenated polarizable force field is capable of reproducing QM relative energies and geometries of both halogen bonds and halogen-hydrogen bond donor interactions at an unprecedented level due to the inclusion of a virtual particle and anisotropic atomic polarizability on the halogen and, notably, the inclusion of Lennard-Jones parameters on the halogen Drude particle. The model was validated on the basis of its ability to accurately reproduce pure solvent properties for halogenated naphthalenes and alkanes, including species analogous to those used as refrigerants. Accordingly, it is anticipated that the model will be applicable for the study of halogenated derivatives in CADD as well as in other chemical and biophysical studies.

Effect of Monovalent Ion Parameters on Molecular Dynamics Simulations of G-Quadruplexes

G-quadruplexes (GQs) are key noncanonical DNA and RNA architectures stabilized by desolvated monovalent cations present in their central channels. We analyze extended atomistic molecular dynamics simulations (∼580 μs in total) of GQs with 11 monovalent cation parametrizations, assessing GQ overall structural stability, dynamics of internal cations, and distortions of the G-tetrad geometries. Majority of simulations were executed with the SPC/E water model; however, test simulations with TIP3P and OPC water models are also reported. The identity and parametrization of ions strongly affect behavior of a tetramolecular d[GGG]₄ GQ, which is unstable with several ion parametrizations. The remaining studied RNA and DNA GQs are structurally stable, though the G-tetrad geometries are always deformed by bifurcated H-bonding in a parametrization-specific manner. Thus, basic 10-μs-scale simulations of fully folded GQs can be safely done with a number of cation parametrizations. However, there are parametrization-specific differences and basic force-field errors affecting the quantitative description of ion-tetrad interactions, which may significantly affect studies of the ion-binding processes and description of the GQ folding landscape. Our d[GGG]₄ simulations indirectly suggest that such studies will also be sensitive to the water models. During exchanges with bulk water, the Na⁺ ions move inside the GQs in a concerted manner, while larger relocations of the K⁺ ions are typically separated. We suggest that the Joung-Cheatham SPC/E K⁺ parameters represent a safe choice in simulation studies of GQs, though variation of ion parameters can be used for specific simulation goals.

Polarizable Force Field for DNA Based on the Classical Drude Oscillator: II. Microsecond Molecular Dynamics Simulations of Duplex DNA

The structure and dynamics of DNA are governed by a sensitive balance between base stacking and pairing, hydration, and interactions with ions. Force-field models that include explicit representations of electronic polarization are capable of more accurately modeling the subtle details of these interactions versus commonly used additive force fields. In this work, we validate our recently refined polarizable force field for DNA based on the classical Drude oscillator model, in which electronic degrees of freedom are represented as negatively charged particles attached to their parent atoms via harmonic springs. The previous version of the force field, called Drude-2013, produced stable A- and B-DNA trajectories on the order of hundreds of nanoseconds, but deficiencies were identified that included weak base stacking ultimately leading to distortion of B-DNA duplexes and unstable Z-DNA. As a result of extensive refinement of base nonbonded terms and bonded parameters in the deoxyribofuranose sugar and phosphodiester backbone, we demonstrate that the new version of the Drude DNA force field is capable of simulating A- and B-forms of DNA on the microsecond time scale and the resulting conformational ensembles agree well with a broad set of experimental properties, including solution X-ray scattering profiles. In addition, simulations of Z-form duplex DNA in its crystal environment are stable on the order of 100 ns. The revised force field is to be called Drude-2017.

Polarizable Force Field for DNA Based on the Classical Drude Oscillator: I. Refinement Using Quantum Mechanical Base Stacking and Conformational Energetics

Empirical force fields seek to relate the configuration of a set of atoms to its energy, thus yielding the forces governing its dynamics, using classical physics rather than more expensive quantum mechanical calculations that are computationally intractable for large systems. Most force fields used to simulate biomolecular systems use fixed atomic partial charges, neglecting the influence of electronic polarization, instead making use of a mean-field approximation that may not be transferable across environments. Recent hardware and software developments make polarizable simulations feasible, and to this end, polarizable force fields represent the next generation of molecular dynamics simulation technology. In this work, we describe the refinement of a polarizable force field for DNA based on the classical Drude oscillator model by targeting quantum mechanical interaction energies and conformational energy profiles of model compounds necessary to build a complete DNA force field. The parametrization strategy employed in the present work seeks to correct weak base stacking in A- and B-DNA and the unwinding of Z-DNA observed in the previous version of the force field, called Drude-2013. Refinement of base nonbonded terms and reparametrization of dihedral terms in the glycosidic linkage, deoxyribofuranose rings, and important backbone torsions resulted in improved agreement with quantum mechanical potential energy surfaces. Notably, we expand on previous efforts by explicitly including Z-DNA conformational energetics in the refinement.

DNA dielectrophoresis: Theory and applications a review

Dielectrophoresis is the migration of an electrically polarizable particle in an inhomogeneous electric field. This migration can be exploited for several applications with (bio)molecules or cells. Dielectrophoresis is a noninvasive technique; therefore, it is very convenient for (selective) manipulation of (bio)molecules or cells. In this review, we will focus on DNA dielectrophoresis as this technique offers several advantages in trapping and immobilization, separation and purification, and analysis of DNA molecules. We present and discuss the underlying theory of the most important forces that have to be considered for applications with dielectrophoresis. Moreover, a review of DNA dielectrophoresis applications is provided to present the state-of-the-art and to offer the reader a perspective of the advances and current limitations of DNA dielectrophoresis.

Noncanonical α/γ Backbone Conformations in RNA and the Accuracy of Their Description by the AMBER Force Field

The sugar-phosphate backbone of RNA can exist in diverse rotameric substates, giving RNA molecules enormous conformational variability. The most frequent noncanonical backbone conformation in RNA is α/γ = t/t, which is derived from the canonical backbone by a crankshaft motion and largely preserves the standard geometry of the RNA duplex. A similar conformation also exists in DNA, where it has been extensively studied and shown to be involved in DNA-protein interactions. However, the function of the α/γ = t/t conformation in RNA is poorly understood. Here, we present molecular dynamics simulations of several prototypical RNA structures obtained from X-ray and NMR experiments, including canonical and mismatched RNA duplexes, UUCG and GAGA tetraloops, Loop E, the sarcin-ricin loop, a parallel guanine quadruplex, and a viral pseudoknot. The stability of various noncanonical α/γ backbone conformations was analyzed with two AMBER force fields, ff99bsc0χ_OL3 and ff99bsc0χ_OL3 with the recent εζ_OL1 and β_OL1 corrections for DNA. Although some α/γ substates were stable with seemingly well-described equilibria, many were unstable in our simulations. Notably, the most frequent noncanonical conformer α/γ = t/t was unstable in both tested force fields. Possible reasons for this instability are discussed. Our work reveals a potentially important artifact in RNA force fields and highlights a need for further force field refinement.

Extracting water and ion distributions from solution x-ray scattering experiments

Small-angle X-ray scattering measurements can provide valuable information about the solvent environment around biomolecules, but it can be difficult to extract solvent-specific information from observed intensity profiles. Intensities are proportional to the square of scattering amplitudes, which are complex quantities. Amplitudes in the forward direction are real, and the contribution from a solute of known structure (and from the waters it excludes) can be estimated from theory; hence, the amplitude arising from the solvent environment can be computed by difference. We have found that this "square root subtraction scheme" can be extended to non-zero q values, out to 0.1 Å(-1) for the systems considered here, since the phases arising from the solute and from the water environment are nearly identical in this angle range. This allows us to extract aspects of the water and ion distributions (beyond their total numbers), by combining experimental data for the complete system with calculations for the solutes. We use this approach to test molecular dynamics and integral-equation (3D-RISM (three-dimensional reference interaction site model)) models for solvent structure around myoglobin, lysozyme, and a 25 base-pair duplex DNA. Comparisons can be made both in Fourier space and in terms of the distribution of interatomic distances in real space. Generally, computed solvent distributions arising from the MD simulations fit experimental data better than those from 3D-RISM, even though the total small-angle X-ray scattering patterns are very similar; this illustrates the potential power of this sort of analysis to guide the development of computational models.

An Estimation of Hybrid Quantum Mechanical Molecular Mechanical Polarization Energies for Small Molecules Using Polarizable Force-Field Approaches

In this work, we report two polarizable molecular mechanics (polMM) force field models for estimating the polarization energy in hybrid quantum mechanical molecular mechanical (QM/MM) calculations. These two models, named the potential of atomic charges (PAC) and potential of atomic dipoles (PAD), are formulated from the ab initio quantum mechanical (QM) response kernels for the prediction of the QM density response to an external molecular mechanical (MM) environment (as described by external point charges). The PAC model is similar to fluctuating charge (FQ) models because the energy depends on external electrostatic potential values at QM atomic sites; the PAD energy depends on external electrostatic field values at QM atomic sites, resembling induced dipole (ID) models. To demonstrate their uses, we apply the PAC and PAD models to 12 small molecules, which are solvated by TIP3P water. The PAC model reproduces the QM/MM polarization energy with a R² value of 0.71 for aniline (in 10,000 TIP3P water configurations) and 0.87 or higher for other 11 solute molecules, while the PAD model has a much better performance with R² values of 0.98 or higher. The PAC model reproduces reference QM/MM hydration free energies for 12 solute molecules with a RMSD of 0.59 kcal/mol. The PAD model is even more accurate, with a much smaller RMSD of 0.12 kcal/mol, with respect to the reference. This suggests that polarization effects, including both local charge distortion and intramolecular charge transfer, can be well captured by induced dipole type models with proper parametrization.

Balancing the Interactions of Mg2+ in Aqueous Solution and with Nucleic Acid Moieties For a Polarizable Force Field Based on the Classical Drude Oscillator Model

Mg²⁺ ions are important in biological systems, particularly in stabilizing compact RNA folds. Mg²⁺ is strongly polarizing, and representing its interactions in heterogeneous environments is a challenge for empirical force field development. To date, the most commonly used force fields in molecular dynamics simulations utilize a pairwise-additive approximation for electrostatic interactions, which cannot account for the significant polarization response in systems containing Mg²⁺. In the present work, we refine the interactions of Mg²⁺ with water, Cl^- ions, and nucleic acid moieties using a polarizable force field based on the classical Drude oscillator model. By targeting gas-phase quantum mechanical interaction energies and geometries of hydrated complexes, as well as condensed-phase osmotic pressure calculations, we present a model for Mg²⁺ that yields quantitative agreement with experimental measurements of water dissociation free energy and osmotic pressure across a broad range of concentrations. Notable is the direct modeling of steric repulsion between the water Drude oscillators and Mg²⁺ to treat the Pauli exclusion effects associated with overlap of the electron clouds of water molecules in the first hydration shell around Mg²⁺. Combined with the refined interactions with nucleic acid moieties, the present model represents a significant advancement in simulating nucleic acid systems containing Mg²⁺.

Divalent Ion Dependent Conformational Changes in an RNA Stem-Loop Observed by Molecular Dynamics

We compare the performance of five magnesium (Mg²⁺) ion models in simulations of an RNA stem loop which has an experimentally determined divalent ion dependent conformational shift. We show that despite their differences in parametrization and resulting van der Waals terms, including differences in the functional form of the nonbonded potential, when the RNA adopts its folded conformation, all models behave similarly across ten independent microsecond length simulations with each ion model. However, when the entire structure ensemble is accounted for, chelation of Mg²⁺ to RNA is seen in three of the five models, most egregiously and likely artifactual in simulations using a 12-6-4 model for the Lennard-Jones potential. Despite the simple nature of the fixed point-charge and van der Waals sphere models employed, and with the exception of the likely oversampled directed chelation of the 12-6-4 potential models, RNA–Mg²⁺ interactions via first shell water molecules are surprisingly well described by modern parameters, allowing us to observe the spontaneous conformational shift from Mg²⁺ free RNA to Mg²⁺ associated RNA structure in unrestrained molecular dynamics simulations.

Nucleic acid polymeric properties and electrostatics: Directly comparing theory and simulation with experiment

Nucleic acids are biopolymers that carry genetic information and are also involved in various gene regulation functions such as gene silencing and protein translation. Because of their negatively charged backbones, nucleic acids are polyelectrolytes. To adequately understand nucleic acid folding and function, we need to properly describe its i) polymer/polyelectrolyte properties and ii) associating ion atmosphere. While various theories and simulation models have been developed to describe nucleic acids and the ions around them, many of these theories/simulations have not been well evaluated due to complexities in comparison with experiment. In this review, I discuss some recent experiments that have been strategically designed for straightforward comparison with theories and simulation models. Such data serve as excellent benchmarks to identify limitations in prevailing theories and simulation parameters.

An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications

Molecular mechanics force fields that explicitly account for induced polarization represent the next generation of physical models for molecular dynamics simulations. Several methods exist for modeling induced polarization, and here we review the classical Drude oscillator model, in which electronic degrees of freedom are modeled by charged particles attached to the nuclei of their core atoms by harmonic springs. We describe the latest developments in Drude force field parametrization and application, primarily in the last 15 years. Emphasis is placed on the Drude-2013 polarizable force field for proteins, DNA, lipids, and carbohydrates. We discuss its parametrization protocol, development history, and recent simulations of biologically interesting systems, highlighting specific studies in which induced polarization plays a critical role in reproducing experimental observables and understanding physical behavior. As the Drude oscillator model is computationally tractable and available in a wide range of simulation packages, it is anticipated that use of these more complex physical models will lead to new and important discoveries of the physical forces driving a range of chemical and biological phenomena.