Base pair dynamics, electrostatics, and thermodynamics at the LTR-III quadruplex:duplex junction.

G-quadruplexes (GQs) play key regulatory roles within the human genome and have also been identified to play similar roles in other eukaryotes, bacteria, archaea, and viruses. Human immunodeficiency virus 1, the etiological agent of acquired immunodeficiency syndrome, can form two GQs in its long terminal repeat (LTR) promoter region, each of which act to regulate viral gene expression in opposing manners. The major LTR GQ, called LTR-III, is a distinct hybrid GQ containing a 12-nucleotide duplex loop attached to the quadruplex motif. The resulting quadruplex:duplex junction (QDJ) has been hypothesized to serve as a selective drug targeting site. To better understand the dynamics of this QDJ, we performed conventional and enhanced-sampling molecular dynamics simulations using the Drude-2017 force field. We observed unbiased and reversible formation of additional base pairs in the QDJ, between Ade4:Thy14 and Gua3:Thy14. Both base pairs were electrostatically favored, but geometric constraints within the junction may drive the formation of, and preference for, the Ade4:Thy14 base pair. Finally, we demonstrated that the base pairs are separated only by small energy barriers that may enable transitions between both base-paired states. Together, these simulations provide new insights into the dynamics, electrostatics, and thermodynamics of the LTR-III QDJ.

Structural Dynamics of the Methyl-Coenzyme M Reductase Active Site Are Influenced by Coenzyme F430 Modifications.

Methyl-coenzyme M reductase (MCR) is a central player in methane biogeochemistry, governing methanogenesis and the anaerobic oxidation of methane (AOM) in methanogens and anaerobic methanotrophs (ANME), respectively. The prosthetic group of MCR is coenzyme F430, a nickel-containing tetrahydrocorphin. Several modified versions of F430 have been discovered, including the 172-methylthio-F430 (mtF430) used by ANME-1 MCR. Here, we employ molecular dynamics (MD) simulations to investigate the active site dynamics of MCR from Methanosarcina acetivorans and ANME-1 when bound to the canonical F430 compared to 172-thioether coenzyme F430 variants and substrates (methyl-coenzyme M and coenzyme B) for methane formation. Our simulations highlight the importance of the Gln to Val substitution in accommodating the 172 methylthio modification in ANME-1 MCR. Modifications at the 172 position disrupt the canonical substrate positioning in M. acetivorans MCR. However, in some replicates, active site reorganization to maintain substrate positioning suggests that the modified F430 variants could be accommodated in a methanogenic MCR. We additionally report the first quantitative estimate of MCR intrinsic electric fields that are pivotal in driving methane formation. Our results suggest that the electric field aligned along the CH3-S-CoM thioether bond facilitates homolytic bond cleavage, coinciding with the proposed catalytic mechanism. Structural perturbations, however, weaken and misalign these electric fields, emphasizing the importance of the active site structure in maintaining their integrity. In conclusion, our results deepen the understanding of MCR active site dynamics, the enzyme's organizational role in intrinsic electric fields for catalysis, and the interplay between active site structure and electrostatics.

Pyroglutamylation modulates electronic properties and the conformational ensemble of the amyloid ß-peptide.

Alzheimer's disease (AD) is a neurodegenerative disorder that is characterized by the formation of extracellular amyloid-ß (Aß) plaques. The underlying cause of AD is unknown, however, post-translational modifications (PTMs) of Aß have been found in AD patients and are thought to play a role in protein aggregation. One such PTM is pyroglutamylation, which can occur at two sites in Aß, Glu3 and Glu11. This modification of Aß involves the truncation and charge-neutralization of N-terminal glutamate, causing Aß to become more hydrophobic and prone to aggregation. The molecular mechanism by which the introduction of pyroglutamate (pE) promotes aggregation has not been determined. To gain a greater understanding of the role that charge neutralization and truncation of the N-terminus plays on Aß conformational sampling, we used the Drude polarizable force field (FF) to perform molecular dynamics simulations on AßpE3-42 and AßpE11-42 and comparing their properties to previous simulations of Aß1-42. The Drude polarizable FF allows for a more accurate representation of electrostatic interactions, therefore providing novel insights into the role that charge plays in protein dynamics. Here, we report the parametrization of pE in the Drude polarizable FF and the effect of pyroglutamylation on Aß. We found that AßpE3-42 and AßpE11-42 alter the permanent and induced dipoles of the peptide. Specifically, we found that AßpE3-42 and AßpE11-42 have modification-specific backbone and sidechain polarization response and perturbed solvation properties that shift the Aß conformational ensemble.

Differences in Conformational Sampling and Intrinsic Electric Fields Drive Ion Binding in Telomeric and TERRA G-Quadruplexes.

The formation of G-quadruplexes (GQs) occurs in guanine-rich sequences of DNA and RNA, producing highly stable and structurally diverse noncanonical nucleic acid structures. GQs play crucial roles in regulating transcription, translation, and replication and maintaining the genome, among others; thus, changes to their structures can lead to diseases such as cancer. Previous studies using polarizable molecular dynamics simulations have shown differences in ion binding properties between telomeric and telomeric repeat-containing RNA GQs despite architectural similarities. Here, we used volume-based metadynamics and repulsive potential simulations in conjunction with polarizable force fields to quantify the impact of ion binding on the GQ dynamics and ion binding free energies. Furthermore, we describe how GQs exert electric fields on their surroundings to link dynamics with variations in the electronic structure. Our findings provide new insights into the energetic, physical, and conformational properties of GQs and expose subtle but important differences between DNA and RNA GQs with the same fold.

charmm2gmx: An Automated Method to Port the CHARMM Additive Force Field to GROMACS.

CHARMM is one of the most widely used biomolecular force fields. Although developed in close connection with a dedicated molecular simulation engine of the same name, it is also usable with other codes. GROMACS is a well-established, highly optimized, and multipurpose software for molecular dynamics, versatile enough to accommodate many different force field potential functions and the associated algorithms. Due to conceptional differences related to software design and the large amount of numeric data inherent to residue topologies and parameter sets, conversion from one software format to another is not straightforward. Here, we present an automated and validated means to port the CHARMM force field to a format read by the GROMACS engine, harmonizing the different capabilities of the two codes in a self-documenting and reproducible way with a bare minimum of user interaction required. Being based entirely on the upstream data files, the presented approach does not involve any hard-coded data, in contrast with previous attempts to solve the same problem. The heuristic approach used for perceiving the local internal geometry is directly applicable for analogous transformations of other force fields.

TUPÃ: Electric field analyses for molecular simulations.

We introduce TUPÃ, a Python-based algorithm to calculate and analyze electric fields in molecular simulations. To demonstrate the features in TUPÃ, we present three test cases in which the orientation and magnitude of the electric field exerted by biomolecules help explain biological phenomena or observed kinetics. As part of TUPÃ, we also provide a PyMOL plugin to help researchers visualize how electric fields are organized within the simulation system. The code is freely available and can be obtained at https://mdpoleto.github.io/tupa/.

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.

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.

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.

Sequential Bending and Twisting around C-C Single Bonds by Mechanical Lifting of a Pre-Adsorbed Polymer.

Bending and twisting around carbon-carbon single bonds are ubiquitous in natural and synthetic polymers. Force-induced changes were so far not measured at the single-monomer level, owing to limited ways to apply local forces. We quantified down to the submolecular level the mechanical response within individual poly-pyrenylene chains upon their detachment from a gold surface with an atomic force microscope at 5 K. Computer simulations based on a dedicated force field reproduce the experimental traces and reveal symmetry-broken bent and rotated conformations of the sliding physisorbed segment besides steric hindrance of the just lifted monomer. Our study also shows that the tip-molecule bond remains intact but remarkably soft and links force variations to complex but well-defined conformational changes.

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.

HIV-1 Env gp41 Transmembrane Domain Dynamics Are Modulated by Lipid, Water, and Ion Interactions.

The gp41 transmembrane domain (TMD) of the envelope glycoprotein of the human immunodeficiency virus modulates the conformation of the viral envelope spike, the only druggable target on the surface of the virion. Targeting the envelope glycoprotein with small-molecule and antibody therapies requires an understanding of gp41 TMD dynamics, which is often challenging given the difficulties in describing native membrane properties. Here, atomistic molecular dynamics simulations of a trimeric, prefusion gp41 TMD in a model, asymmetric viral membrane that mimics the native viral envelope were performed. Water and chloride ions were observed to permeate the membrane and interact with the highly conserved arginine bundle, (R696)3, at the center of the membrane and influenced TMD stability by creating a network of hydrogen bonds and electrostatic interactions. We propose that this (R696)3 - water - anion network plays an important role in viral fusion with the host cell by modulating protein conformational changes within the membrane. Additionally, R683 and R707 at the exofacial and cytofacial membrane-water interfaces, respectively, are anchored in the lipid headgroup region and serve as a junction point for stabilization of the termini. The membrane thins as a result of the tilting of the gp41 trimer with nearby lipids increasing in volume, leading to an entropic driving force for TMD conformational change. These results provide additional detail and perspective on the influence of certain lipid types on TMD dynamics and a rationale for targeting key residues of the TMD for therapeutic design. These insights into the molecular details of TMD membrane anchoring will build toward a greater understanding of the dynamics that lead to viral fusion with the host cell.

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.

Molecular dynamics simulations using the drude polarizable force field on GPUs with OpenMM: Implementation, validation, and benchmarks.

Presented is the implementation of the Drude force field in the open-source OpenMM simulation package allowing for access to graphical processing unit (GPU) hardware. In the Drude model, electronic degrees of freedom are represented by negatively charged particles attached to their parent atoms via harmonic springs, such that extra computational overhead comes from these additional particles and virtual sites representing lone pairs on electronegative atoms, as well as the associated thermostat and integration algorithms. This leads to an approximately fourfold increase in computational demand over additive force fields. However, by making the Drude model accessible to consumer-grade desktop GPU hardware it will be possible to perform simulations of one microsecond or more in less than a month, indicating that the barrier to employ polarizable models has largely been removed such that polarizable simulations with the classical Drude model are readily accessible and practical.

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.

Small Molecule Thermochemistry: A Tool for Empirical Force Field Development.

Spectroscopic analysis of compounds is typically combined with density functional theory, for instance, for assigning vibrational frequencies, limiting application to relatively small compounds. Accurate classical force fields could, in principle, complement these quantum-chemical tools. A relatively simple way to validate vibrational frequencies is by computing thermochemical properties. We present such a validation for over 1800 small molecules using the harmonic approximation. Two popular empirical force fields (GAFF and CGenFF) are compared to experimental data and results from Gaussian-4 quantum-chemical calculations. Frequency scaling factors of 1.035 (CGenFF) and 1.018 (GAFF) are derived from the zero-point energies. The force field calculations have larger deviation from experiment than the G4 method for standard entropy, but for heat capacity the results are comparable. For internal thermal energy and zero-point energy the deviations from G4 are relatively small. The work suggests that with some tuning force fields could indeed complement DFT in spectroscopical applications.

Structure and Dynamics of FosA-Mediated Fosfomycin Resistance in Klebsiella pneumoniae and Escherichia coli.

Fosfomycin exhibits broad-spectrum antibacterial activity and is being reevaluated for the treatment of extensively drug-resistant pathogens. Its activity in Gram-negative organisms, however, can be compromised by expression of FosA, a metal-dependent transferase that catalyzes the conjugation of glutathione to fosfomycin, rendering the antibiotic inactive. In this study, we solved the crystal structures of two of the most clinically relevant FosA enzymes: plasmid-encoded FosA3 from Escherichia coli and chromosomally encoded FosA from Klebsiella pneumoniae (FosAKP). The structure, molecular dynamics, catalytic activity, and fosfomycin resistance of FosA3 and FosAKP were also compared to those of FosA from Pseudomonas aeruginosa (FosAPA), for which prior crystal structures exist. E. coli TOP10 transformants expressing FosA3 and FosAKP conferred significantly greater fosfomycin resistance (MIC, >1,024 µg/ml) than those expressing FosAPA (MIC, 16 µg/ml), which could be explained in part by the higher catalytic efficiencies of the FosA3 and FosAKP enzymes. Interestingly, these differences in enzyme activity could not be attributed to structural differences at their active sites. Instead, molecular dynamics simulations and hydrogen-deuterium exchange experiments with FosAKP revealed dynamic interconnectivity between its active sites and a loop structure that extends from the active site of each monomer and traverses the dimer interface. This dimer interface loop is longer and more extended in FosAKP and FosA3 than in FosAPA, and kinetic analyses of FosAKP and FosAPA loop-swapped chimeric enzymes highlighted its importance in FosA activity. Collectively, these data yield novel insights into fosfomycin resistance that could be leveraged to develop new strategies to inhibit FosA and potentiate fosfomycin activity.

DIRECT-ID: An automated method to identify and quantify conformational variations--application to ß2 -adrenergic GPCR.

The conformational dynamics of a macromolecule can be modulated by a number of factors, including changes in environment, ligand binding, and interactions with other macromolecules, among others. We present a method that quantifies the differences in macromolecular conformational dynamics and automatically extracts the structural features responsible for these changes. Given a set of molecular dynamics (MD) simulations of a macromolecule, the norms of the differences in covariance matrices are calculated for each pair of trajectories. A matrix of these norms thus quantifies the differences in conformational dynamics across the set of simulations. For each pair of trajectories, covariance difference matrices are parsed to extract structural elements that undergo changes in conformational properties. As a demonstration of its applicability to biomacromolecular systems, the method, referred to as DIRECT-ID, was used to identify relevant ligand-modulated structural variations in the ß2 -adrenergic (ß2 AR) G-protein coupled receptor. Micro-second MD simulations of the ß2 AR in an explicit lipid bilayer were run in the apo state and complexed with the ligands: BI-167107 (agonist), epinephrine (agonist), salbutamol (long-acting partial agonist), or carazolol (inverse agonist). Each ligand modulated the conformational dynamics of ß2 AR differently and DIRECT-ID analysis of the inverse-agonist vs. agonist-modulated ß2 AR identified residues known through previous studies to selectively propagate deactivation/activation information, along with some previously unidentified ligand-specific microswitches across the GPCR. This study demonstrates the utility of DIRECT-ID to rapidly extract functionally relevant conformational dynamics information from extended MD simulations of large and complex macromolecular systems.

Parametrization of halogen bonds in the CHARMM general force field: Improved treatment of ligand-protein interactions.

A halogen bond is a highly directional, non-covalent interaction between a halogen atom and another electronegative atom. It arises due to the formation of a small region of positive electrostatic potential opposite the covalent bond to the halogen, called the 'sigma hole.' Empirical force fields in which the electrostatic interactions are represented by atom-centered point charges cannot capture this effect because halogen atoms usually carry a negative charge and therefore interact unfavorably with other electronegative atoms. A strategy to overcome this problem is to attach a positively charged virtual particle to the halogen. In this work, we extend the additive CHARMM General Force Field (CGenFF) to include such interactions in model systems of phenyl-X, with X being Cl, Br or I including di- and trihalogenated species. The charges, Lennard-Jones parameters, and halogen-virtual particle distances were optimized to reproduce the orientation dependence of quantum mechanical interaction energies with water, acetone, and N-methylacetamide as well as experimental pure liquid properties and relative hydration free energies with respect to benzene. The resulting parameters were validated in molecular dynamics simulations on small-molecule crystals and on solvated protein-ligand complexes containing halogenated compounds. The inclusion of positive virtual sites leads to better agreement across experimental observables, including preservation of ligand binding poses as a direct result of the improved representation of halogen bonding.

Implementation of extended Lagrangian dynamics in GROMACS for polarizable simulations using the classical Drude oscillator model.

Explicit treatment of electronic polarization in empirical force fields used for molecular dynamics simulations represents an important advancement in simulation methodology. A straightforward means of treating electronic polarization in these simulations is the inclusion of Drude oscillators, which are auxiliary, charge-carrying particles bonded to the cores of atoms in the system. The additional degrees of freedom make these simulations more computationally expensive relative to simulations using traditional fixed-charge (additive) force fields. Thus, efficient tools are needed for conducting these simulations. Here, we present the implementation of highly scalable algorithms in the GROMACS simulation package that allow for the simulation of polarizable systems using extended Lagrangian dynamics with a dual Nosé-Hoover thermostat as well as simulations using a full self-consistent field treatment of polarization. The performance of systems of varying size is evaluated, showing that the present code parallelizes efficiently and is the fastest implementation of the extended Lagrangian methods currently available for simulations using the Drude polarizable force field.