Justin A. Lemkul

CHARMM-GUI Input Generator for NAMD, Gromacs, Amber, Openmm, and CHARMM/OpenMM Simulations using the CHARMM36 Additive Force Field

Proper treatment of nonbonded interactions is essential for the accuracy of molecular dynamics (MD) simulations, especially in studies of lipid bilayers. The use of the CHARMM36 force field (C36 FF) in different MD simulation programs can result in disagreements with published simulations performed with CHARMM due to differences in the protocols used to treat the long-range and 1-4 nonbonded interactions. In this study, we systematically test the use of the C36 lipid FF in NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM. A wide range of Lennard-Jones (LJ) cutoff schemes and integrator algorithms were tested to find the optimal simulation protocol to best match bilayer properties of six lipids with varying acyl chain saturation and head groups. MD simulations of a 1,2-dipalmitoyl-sn-phosphatidylcholine (DPPC) bilayer were used to obtain the optimal protocol for each program. MD simulations with all programs were found to reasonably match the DPPC bilayer properties (surface area per lipid, chain order parameters, and area compressibility modulus) obtained using the standard protocol used in CHARMM as well as from experiments. The optimal simulation protocol was then applied to the other five lipid simulations and resulted in excellent agreement between results from most simulation programs as well as with experimental data. AMBER compared least favorably with the expected membrane properties, which appears to be due to its use of the hard-truncation in the LJ potential versus a force-based switching function used to smooth the LJ potential as it approaches the cutoff distance. The optimal simulation protocol for each program has been implemented in CHARMM-GUI. This protocol is expected to be applicable to the remainder of the additive C36 FF including the proteins, nucleic acids, carbohydrates, and small molecules.

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.

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.

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.

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.

Characterization of Mg2+ Distributions around RNA in Solution

Binding of metal ions is an important factor governing the folding and dynamics of RNA. Shielding of charges in the polyanionic backbone allows RNA to adopt a diverse range of folded structures that give rise to their many functions within the cell. Some RNA sequences fold only in the presence of Mg2+, which may be bound via direct interactions or occupy the more diffuse “ion atmosphere” around the RNA. To understand the driving forces for RNA folding, it is important to be able to fully characterize the distribution of metal ions around the RNA. In this work, a combined Grand Canonical Monte Carlo-Molecular Dynamics (GCMC-MD) method is applied to characterize Mg2+ distributions around folded RNA structures. The GCMC-MD approach identifies known inner- and outer-shell Mg2+ coordination, while also predicting new regions occupied by Mg2+ that are not observed in crystal structures but that may be relevant in solution, including the case of the Mg2+ riboswitch, for which alternate Mg2+ binding sites may have implications for its function. This work represents a significant step forward in establishing a structural and thermodynamic description of RNA–Mg2+ interactions and their role in RNA structure and function.

Induced Dipole–Dipole Interactions Influence the Unfolding Pathways of Wild-Type and Mutant Amyloid β-Peptides

Amyloid-forming proteins undergo a structural transition from α-helical to disordered conformations and, ultimately, cross-β fibrils. The unfolding and aggregation of the amyloid β-peptide (Aβ) have been implicated in the development and progression of Alzheimer’s disease (AD) and cerebral amyloid angiopathy (CAA). However, the events underlying the initial structural transition leading to the disease state remain unclear. Although most cases are sporadic, several genetic variants exist that alter the electrostatic properties of Aβ and lead to more rapid unfolding and more severe phenotypes. In the present study, the enhanced unfolding is shown to be due to the mutated side chains altering the local peptide-bond dipole moments leading to local destabilization of the α-helix, as determined from polarizable molecular dynamics (MD) simulations of wild-type (WT) Aβ fragments and several common mutations. The local perturbation of the helix then leads to progressive unwinding of the α-helix in a cooperative fashion due to decreases in adjacent (i ± 1) and hydrogen-bonded (i + 4) peptide-bond dipole moments. Side-chain dynamics, subsequent variations in dipole moments, and ultimately the response in the peptide-bond dipole moments are all modulated by solvent dielectric properties based on simulations in water versus ethanol. The polarizable simulation results, along with simulations using the additive CHARMM36 force field, further indicate that cooperativity due to the alignment of peptide bonds leading to enhanced dipole moments is a fundamental force in stabilizing α-helices.

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.

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 (β2AR) G‐protein coupled receptor. Micro‐second MD simulations of the β2AR 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 β2AR differently and DIRECT‐ID analysis of the inverse‐agonist vs. agonist‐modulated β2AR 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.

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

ABSTRACT 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.