RNA folding using quantum computers.

The 3-dimensional fold of an RNA molecule is largely determined by patterns of intramolecular hydrogen bonds between bases. Predicting the base pairing network from the sequence, also referred to as RNA secondary structure prediction or RNA folding, is a nondeterministic polynomial-time (NP)-complete computational problem. The structure of the molecule is strongly predictive of its functions and biochemical properties, and therefore the ability to accurately predict the structure is a crucial tool for biochemists. Many methods have been proposed to efficiently sample possible secondary structure patterns. Classic approaches employ dynamic programming, and recent studies have explored approaches inspired by evolutionary and machine learning algorithms. This work demonstrates leveraging quantum computing hardware to predict the secondary structure of RNA. A Hamiltonian written in the form of a Binary Quadratic Model (BQM) is derived to drive the system toward maximizing the number of consecutive base pairs while jointly maximizing the average length of the stems. A Quantum Annealer (QA) is compared to a Replica Exchange Monte Carlo (REMC) algorithm programmed with the same objective function, with the QA being shown to be highly competitive at rapidly identifying low energy solutions. The method proposed in this study was compared to three algorithms from literature and, despite its simplicity, was found to be competitive on a test set containing known structures with pseudoknots.

Lipid21: Complex Lipid Membrane Simulations with AMBER.

We extend the modular AMBER lipid force field to include anionic lipids, polyunsaturated fatty acid (PUFA) lipids, and sphingomyelin, allowing the simulation of realistic cell membrane lipid compositions, including raft-like domains. Head group torsion parameters are revised, resulting in improved agreement with NMR order parameters, and hydrocarbon chain parameters are updated, providing a better match with phase transition temperature. Extensive validation runs (0.9 µs per lipid type) show good agreement with experimental measurements. Furthermore, the simulation of raft-like bilayers demonstrates the perturbing effect of increasing PUFA concentrations on cholesterol molecules. The force field derivation is consistent with the AMBER philosophy, meaning it can be easily mixed with protein, small molecule, nucleic acid, and carbohydrate force fields.

mRNA codon optimization with quantum computers.

Reverse translation of polypeptide sequences to expressible mRNA constructs is a NP-hard combinatorial optimization problem. Each amino acid in the protein sequence can be represented by as many as six codons, and the process of selecting the combination that maximizes probability of expression is termed codon optimization. This work investigates the potential impact of leveraging quantum computing technology for codon optimization. A Quantum Annealer (QA) is compared to a standard genetic algorithm (GA) programmed with the same objective function. The QA is found to be competitive in identifying optimal solutions. The utility of gate-based systems is also evaluated using a simulator resulting in the finding that while current generations of devices lack the hardware requirements, in terms of both qubit count and connectivity, to solve realistic problems, future generation devices may be highly efficient.

X-ray Crystallographic and Molecular Dynamic Analyses of Drosophila melanogaster Embryonic Muscle Myosin Define Domains Responsible for Isoform-Specific Properties.

Drosophila melanogaster is a powerful system for characterizing alternative myosin isoforms and modeling muscle diseases, but high-resolution structures of fruit fly contractile proteins have not been determined. Here we report the first x-ray crystal structure of an insect myosin: the D melanogaster skeletal muscle myosin II embryonic isoform (EMB). Using our system for recombinant expression of myosin heavy chain (MHC) proteins in whole transgenic flies, we prepared and crystallized stable proteolytic S1-like fragments containing the entire EMB motor domain bound to an essential light chain. We solved the x-ray crystal structure by molecular replacement and refined the resulting model against diffraction data to 2.2 Å resolution. The protein is captured in two slightly different renditions of the rigor-like conformation with a citrate of crystallization at the nucleotide binding site and exhibits structural features common to myosins of diverse classes from all kingdoms of life. All atom molecular dynamics simulations on EMB in its nucleotide-free state and a derivative homology model containing 61 amino acid substitutions unique to the indirect flight muscle isoform (IFI) suggest that differences in the identity of residues within the relay and the converter that are encoded for by MHC alternative exons 9 and 11, respectively, directly contribute to increased mobility of these regions in IFI relative to EMB. This suggests the possibility that alternative folding or conformational stability within these regions contribute to the observed functional differences in Drosophila EMB and IFI myosins.

Simulating Water Exchange to Buried Binding Sites.

Traditional molecular dynamics (MD) simulations of proteins, which relies on integration of Newton's equations of motion, cannot efficiently equilibrate water occupancy for buried cavities in proteins. This leads to slow convergence of thermodynamic averages for such systems. We have addressed this challenge by efficiently integrating standard Metropolis Monte Carlo (MC) translational water moves with MD in the AMBER simulation package. The translational moves allow water to easily enter or exit buried sites in a thermodynamically correct way during a simulation. To maximize efficiency, the algorithm avoids moves that only interchange waters within the bulk around the protein instead focusing on moves that can transfer water between bulk and the protein interior. In addition, a steric grid allows avoidance of moves that would lead to obvious steric clashes, and a fast grid-based energy evaluation is used to reduce the number of expensive full energy calculations. The potential energy distribution produced using MC/MD was found to be statistically indistinguishable from that of control simulations using only MD, and the algorithm effectively equilibrated water across steric barriers and into binding pockets that are not accessible with pure MD. The MC/MD method introduced here should be of increasing utility for applications spanning protein folding, the elucidation of protein mechanisms, and free energy calculations for computer-aided drug design. It is available in version 18 release of the widely disseminated AMBER simulation package.

pH-dependent conformational dynamics of beta-secretase 1: A molecular dynamics study.

Beta-secretase 1 (BACE-1) is an aspartyl protease implicated in the overproduction of ß-amyloid fibrils responsible for Alzheimer disease. The process of ß-amyloid genesis is known to be pH dependent, with an activity peak between solution pH of 3.5 and 5.5. We have studied the pH-dependent dynamics of BACE-1 to better understand the pH dependent mechanism. We have implemented support for graphics processor unit (GPU) accelerated constant pH molecular dynamics within the AMBER molecular dynamics software package and employed this to determine the relative population of different aspartyl dyad protonation states in the pH range of greatest ß-amyloid production, followed by conventional molecular dynamics to explore the differences among the various aspartyl dyad protonation states. We observed a difference in dynamics between double-protonated, mono-protonated, and double-deprotonated states over the known pH range of higher activity. These differences include Tyr 71-aspartyl dyad proximity and active water lifetime. This work indicates that Tyr 71 stabilizes catalytic water in the aspartyl dyad active site, enabling BACE-1 activity.

GPU-Accelerated Molecular Dynamics and Free Energy Methods in Amber18: Performance Enhancements and New Features.

We report progress in graphics processing unit (GPU)-accelerated molecular dynamics and free energy methods in Amber18. Of particular interest is the development of alchemical free energy algorithms, including free energy perturbation and thermodynamic integration methods with support for nonlinear soft-core potential and parameter interpolation transformation pathways. These methods can be used in conjunction with enhanced sampling techniques such as replica exchange, constant-pH molecular dynamics, and new 12-6-4 potentials for metal ions. Additional performance enhancements have been made that enable appreciable speed-up on GPUs relative to the previous software release.

Fast and flexible gpu accelerated binding free energy calculations within the amber molecular dynamics package.

Alchemical free energy (AFE) calculations based on molecular dynamics (MD) simulations are key tools in both improving our understanding of a wide variety of biological processes and accelerating the design and optimization of therapeutics for numerous diseases. Computing power and theory have, however, long been insufficient to enable AFE calculations to be routinely applied in early stage drug discovery. One of the major difficulties in performing AFE calculations is the length of time required for calculations to converge to an ensemble average. CPU implementations of MD-based free energy algorithms can effectively only reach tens of nanoseconds per day for systems on the order of 50,000 atoms, even running on massively parallel supercomputers. Therefore, converged free energy calculations on large numbers of potential lead compounds are often untenable, preventing researchers from gaining crucial insight into molecular recognition, potential druggability and other crucial areas of interest. Graphics Processing Units (GPUs) can help address this. We present here a seamless GPU implementation, within the PMEMD module of the AMBER molecular dynamics package, of thermodynamic integration (TI) capable of reaching speeds of >140 ns/day for a 44,907-atom system, with accuracy equivalent to the existing CPU implementation in AMBER. The implementation described here is currently part of the AMBER 18 beta code and will be an integral part of the upcoming version 18 release of AMBER. © 2018 Wiley Periodicals, Inc.

Ligand Binding Pathways and Conformational Transitions of the HIV Protease.

It is important to determine the binding pathways and mechanisms of ligand molecules to target proteins to effectively design therapeutic drugs. Molecular dynamics (MD) is a promising computational tool that allows us to simulate protein-drug binding at an atomistic level. However, the gap between the time scales of current simulations and those of many drug binding processes has limited the usage of conventional MD, which has been reflected in studies of the HIV protease. Here, we have applied a robust enhanced simulation method, Gaussian accelerated molecular dynamics (GaMD), to sample binding pathways of the XK263 ligand and associated protein conformational changes in the HIV protease. During two of 10 independent GaMD simulations performed over 500-2500 ns, the ligand was observed to successfully bind to the protein active site. Although GaMD-derived free energy profiles were not fully converged because of insufficient sampling of the complex system, the simulations still allowed us to identify relatively low-energy intermediate conformational states during binding of the ligand to the HIV protease. Relative to the X-ray crystal structure, the XK263 ligand reached a minimum root-mean-square deviation (RMSD) of 2.26 Å during 2.5 µs of GaMD simulation. In comparison, the ligand RMSD reached a minimum of only ∼5.73 Å during an earlier 14 µs conventional MD simulation. This work highlights the enhanced sampling power of the GaMD approach and demonstrates its wide applicability to studies of drug-receptor interactions for the HIV protease and by extension many other target proteins.

A Kepler Workflow Tool for Reproducible AMBER GPU Molecular Dynamics.

With the drive toward high throughput molecular dynamics (MD) simulations involving ever-greater numbers of simulation replicates run for longer, biologically relevant timescales (microseconds), the need for improved computational methods that facilitate fully automated MD workflows gains more importance. Here we report the development of an automated workflow tool to perform AMBER GPU MD simulations. Our workflow tool capitalizes on the capabilities of the Kepler platform to deliver a flexible, intuitive, and user-friendly environment and the AMBER GPU code for a robust and high-performance simulation engine. Additionally, the workflow tool reduces user input time by automating repetitive processes and facilitates access to GPU clusters, whose high-performance processing power makes simulations of large numerical scale possible. The presented workflow tool facilitates the management and deployment of large sets of MD simulations on heterogeneous computing resources. The workflow tool also performs systematic analysis on the simulation outputs and enhances simulation reproducibility, execution scalability, and MD method development including benchmarking and validation.

CRISPR-Cas9 conformational activation as elucidated from enhanced molecular simulations.

CRISPR-Cas9 has become a facile genome editing technology, yet the structural and mechanistic features underlying its function are unclear. Here, we perform extensive molecular simulations in an enhanced sampling regime, using a Gaussian-accelerated molecular dynamics (GaMD) methodology, which probes displacements over hundreds of microseconds to milliseconds, to reveal the conformational dynamics of the endonuclease Cas9 during its activation toward catalysis. We disclose the conformational transition of Cas9 from its apo form to the RNA-bound form, suggesting a mechanism for RNA recruitment in which the domain relocations cause the formation of a positively charged cavity for nucleic acid binding. GaMD also reveals the conformation of a catalytically competent Cas9, which is prone for catalysis and whose experimental characterization is still limited. We show that, upon DNA binding, the conformational dynamics of the HNH domain triggers the formation of the active state, explaining how the HNH domain exerts a conformational control domain over DNA cleavage [Sternberg SH et al. (2015) Nature, 527, 110-113]. These results provide atomic-level information on the molecular mechanism of CRISPR-Cas9 that will inspire future experimental investigations aimed at fully clarifying the biophysics of this unique genome editing machinery and at developing new tools for nucleic acid manipulation based on CRISPR-Cas9.

Combined quantum-mechanical molecular mechanics calculations with NWChem and AMBER: Excited state properties of green fluorescent protein chromophore analogue in aqueous solution.

Combined quantum mechanical molecular mechanics (QM/MM) calculations have become a popular methodology for efficient and accurate description of large molecular systems. In this work we introduce our development of a QM/MM framework based on two well-known codes-NWChem and AMBER. As an initial application area we are focused on excited state properties of small molecules in an aqueous phase using an analogue of the green fluorescent protein (GFP) chromophore as a particular test case. Our approach incorporates high level coupled cluster theory for the analysis of excited states providing a reliable theoretical analysis of effects of an aqueous solvation environment on the photochemical properties of the GFP chromophore. Using a systematic approach, which involves comparison of gas phase and aqueous phase results for different protonation states and conformations, we resolve existing uncertainties regarding the theoretical interpretation of experimental data. We observe that the impact of aqueous environment on charged states generally results in blue shifts of the absorption spectra, but the magnitude of the effect is sensitive to both protonation state and conformation and can be rationalized based on charge movement into the area of higher/lower external electrostatic potentials. At neutral pH levels the experimentally observed absorption signal is most likely coming from the phenol protonated form. Our results also show that the high level electron correlated method is essential for a proper description of excited states of GFP. © 2017 Wiley Periodicals, Inc.

Differential flap dynamics in l,d-transpeptidase2 from mycobacterium tuberculosis revealed by molecular dynamics.

Despite the advances in tuberculosis treatment, TB is still one the most deadly infectious diseases and remains a major global health quandary. Mycobacterium tuberculosis (Mtb) is the only known mycobacterium with a high content of 3→3 crosslinks in the biosynthesis of peptidoglycan, which is negligible in most bacterial species. An Mtb lacking LdtMt2 leads to alteration of the colony morphology and loss of virulence which makes this enzyme an attractive target. Regardless of the vital role of LdtMt2 for cell wall survival, the impact of ligand binding on the dynamics of the ß-hairpin flap is still unknown. Understanding the structural and dynamical behaviour of the flap regions provides clear insight into the design of the effective inhibitors against LdtMt2. Carbapenems, an specific class of ß-lactam family, have been shown to inactivate this enzyme. Herein a comprehensive investigation of the flap dynamics of LdtMt2 complex with substrate and three carbapenems namely, ertapenem, imipenem and meropenem is discussed and analyzed for the first account using 140 ns molecular dynamics simulations. The structural features (RMSD, RMSF and Rg) derived by MD trajectories were analyzed. Distance analysis, particularly tip-tip SER135-ASN167 index, identified conformational changes in terms of flap opening and closure within binding process. Principal component analysis (PCA) was employed to qualitatively understand the divergent effects of different inhibitors on the dominant motion of each residue. To probe different internal dynamics induced by ligand binding, dynamic cross-correlation marix (DCCM) analysis was used. The binding free energies of the selected complexes were assessed using MM-GBSA method and per residue free energy decomposition analysis were performed to characterize the contribution of the key residues to the total binding free energies.

Striking Plasticity of CRISPR-Cas9 and Key Role of Non-target DNA, as Revealed by Molecular Simulations.

The CRISPR (clustered regularly interspaced short palindromic repeats)-Cas9 system recently emerged as a transformative genome-editing technology that is innovating basic bioscience and applied medicine and biotechnology. The endonuclease Cas9 associates with a guide RNA to match and cleave complementary sequences in double stranded DNA, forming an RNA:DNA hybrid and a displaced non-target DNA strand. Although extensive structural studies are ongoing, the conformational dynamics of Cas9 and its interplay with the nucleic acids during association and DNA cleavage are largely unclear. Here, by employing multi-microsecond time scale molecular dynamics, we reveal the conformational plasticity of Cas9 and identify key determinants that allow its large-scale conformational changes during nucleic acid binding and processing. We show how the "closure" of the protein, which accompanies nucleic acid binding, fundamentally relies on highly coupled and specific motions of the protein domains, collectively initiating the prominent conformational changes needed for nucleic acid association. We further reveal a key role of the non-target DNA during the process of activation of the nuclease HNH domain, showing how the nontarget DNA positioning triggers local conformational changes that favor the formation of a catalytically competent Cas9. Finally, a remarkable conformational plasticity is identified as an intrinsic property of the HNH domain, constituting a necessary element that allows for the HNH repositioning. These novel findings constitute a reference for future experimental studies aimed at a full characterization of the dynamic features of the CRISPR-Cas9 system, and-more importantly-call for novel structure engineering efforts that are of fundamental importance for the rational design of new genome-engineering applications.

Data for molecular dynamics simulations of B-type cytochrome c oxidase with the Amber force field.

Cytochrome c oxidase (CcO) is a vital enzyme that catalyzes the reduction of molecular oxygen to water and pumps protons across mitochondrial and bacterial membranes. This article presents parameters for the cofactors of ba3-type CcO that are compatible with the all-atom Amber ff12SB and ff14SB force fields. Specifically, parameters were developed for the CuA pair, heme b, and the dinuclear center that consists of heme a3 and CuB bridged by a hydroperoxo group. The data includes geometries in XYZ coordinate format for cluster models that were employed to compute proton transfer energies and derive bond parameters and point charges for the force field using density functional theory. Also included are the final parameter files that can be employed with the Amber leap program to generate input files for molecular dynamics simulations with the Amber software package. Based on the high resolution (1.8 Å) X-ray crystal structure of the ba3-type CcO from Thermus thermophilus (Protein Data Bank ID number PDB: 3S8F), we built a model that is embedded in a POPC lipid bilayer membrane and solvated with TIP3P water molecules and counterions. We provide PDB data files of the initial model and the equilibrated model that can be used for further studies.

Water exit pathways and proton pumping mechanism in B-type cytochrome c oxidase from molecular dynamics simulations.

Cytochrome c oxidase (CcO) is a vital enzyme that catalyzes the reduction of molecular oxygen to water and pumps protons across mitochondrial and bacterial membranes. While proton uptake channels as well as water exit channels have been identified for A-type CcOs, the means by which water and protons exit B-type CcOs remain unclear. In this work, we investigate potential mechanisms for proton transport above the dinuclear center (DNC) in ba3-type CcO of Thermus thermophilus. Using long-time scale, all-atom molecular dynamics (MD) simulations for several relevant protonation states, we identify a potential mechanism for proton transport that involves propionate A of the active site heme a3 and residues Asp372, His376 and Glu126(II), with residue His376 acting as the proton-loading site. The proposed proton transport process involves a rotation of residue His376 and is in line with experimental findings. We also demonstrate how the strength of the salt bridge between residues Arg225 and Asp287 depends on the protonation state and that this salt bridge is unlikely to act as a simple electrostatic gate that prevents proton backflow. We identify two water exit pathways that connect the water pool above the DNC to the outer P-side of the membrane, which can potentially also act as proton exit transport pathways. Importantly, these water exit pathways can be blocked by narrowing the entrance channel between residues Gln151(II) and Arg449/Arg450 or by obstructing the entrance through a conformational change of residue Tyr136, respectively, both of which seem to be affected by protonation of residue His376.

Solvation thermodynamic mapping of molecular surfaces in AmberTools: GIST.

The expulsion of water from surfaces upon molecular recognition and nonspecific association makes a major contribution to the free energy changes of these processes. In order to facilitate the characterization of water structure and thermodynamics on surfaces, we have incorporated Grid Inhomogeneous Solvation Theory (GIST) into the CPPTRAJ toolset of AmberTools. GIST is a grid-based implementation of Inhomogeneous Fluid Solvation Theory, which analyzes the output from molecular dynamics simulations to map out solvation thermodynamic and structural properties on a high-resolution, three-dimensional grid. The CPPTRAJ implementation, called GIST-cpptraj, has a simple, easy-to-use command line interface, and is open source and freely distributed. We have also developed a set of open-source tools, called GISTPP, which facilitate the analysis of GIST output grids. Tutorials for both GIST-cpptraj and GISTPP can be found at ambermd.org. © 2016 Wiley Periodicals, Inc.

A broken-symmetry density functional study of structures, energies, and protonation states along the catalytic O-O bond cleavage pathway in ba3 cytochrome c oxidase from Thermus thermophilus.

Broken-symmetry density functional calculations have been performed on the [Fea3, CuB] dinuclear center (DNC) of ba3 cytochrome c oxidase from Thermus thermophilus in the states of [Fea3(3+)-(HO2)(-)-CuB(2+), Tyr237(-)] and [Fea3(4+)[double bond, length as m-dash]O(2-), OH(-)-CuB(2+), Tyr237˙], using both PW91-D3 and OLYP-D3 functionals. Tyr237 is a special tyrosine cross-linked to His233, a ligand of CuB. The calculations have shown that the DNC in these states strongly favors the protonation of His376, which is above propionate-A, but not of the carboxylate group of propionate-A. The energies of the structures obtained by constrained geometry optimizations along the O-O bond cleavage pathway between [Fea3(3+)-(O-OH)(-)-CuB(2+), Tyr237(-)] and [Fea3(4+)[double bond, length as m-dash]O(2-)HO(-)-CuB(2+), Tyr237˙] have also been calculated. The transition of [Fea3(3+)-(O-OH)(-)-CuB(2+), Tyr237(-)] → [Fea3(4+)[double bond, length as m-dash]O(2-)HO(-)-CuB(2+), Tyr237˙] shows a very small barrier, which is less than 3.0/2.0 kcal mol(-1) in PW91-D3/OLYP-D3 calculations. The protonation state of His376 does not affect this O-O cleavage barrier. The rate limiting step of the transition from state A (in which O2 binds to Fea3(2+)) to state PM ([Fea3(4+)[double bond, length as m-dash]O(2-), OH(-)-CuB(2+), Tyr237˙], where the O-O bond is cleaved) in the catalytic cycle is, therefore, the proton transfer originating from Tyr237 to O-O to form the hydroperoxo [Fea3(3+)-(O-OH)(-)-CuB(2+), Tyr237(-)] state. The importance of His376 in proton uptake and the function of propionate-A/neutral-Asp372 as a gate to prevent the proton from back-flowing to the DNC are also shown.

Simulation of lipid bilayer self-assembly using all-atom lipid force fields.

In this manuscript we expand significantly on our earlier communication by investigating the bilayer self-assembly of eight different types of phospholipids in unbiased molecular dynamics (MD) simulations using three widely used all-atom lipid force fields. Irrespective of the underlying force field, the lipids are shown to spontaneously form stable lamellar bilayer structures within 1 microsecond, the majority of which display properties in satisfactory agreement with the experimental data. The lipids self-assemble via the same general mechanism, though at formation rates that differ both between lipid types, force fields and even repeats on the same lipid/force field combination. In addition to zwitterionic phosphatidylcholine (PC) and phosphatidylethanolamine (PE) lipids, anionic phosphatidylserine (PS) and phosphatidylglycerol (PG) lipids are represented. To our knowledge this is the first time bilayer self-assembly of phospholipids with negatively charged head groups is demonstrated in all-atom MD simulations.

Single Active Site Mutation Causes Serious Resistance of HIV Reverse Transcriptase to Lamivudine: Insight from Multiple Molecular Dynamics Simulations.

Molecular dynamics simulations, binding free energy calculations, principle component analysis (PCA), and residue interaction network analysis were employed in order to investigate the molecular mechanism of M184I single mutation which played pivotal role in making the HIV-1 reverse transcriptase (RT) totally resistant to lamivudine. Results showed that single mutations at residue 184 of RT caused (1) distortion of the orientation of lamivudine in the active site due to the steric conflict between the oxathiolane ring of lamivudine and the side chain of beta-branched amino acids Ile at position 184 which, in turn, perturbs inhibitor binding, (2) decrease in the binding affinity by (~8 kcal/mol) when compared to the wild-type, (3) variation in the overall enzyme motion as evident from the PCA for both systems, and (4) distortion of the hydrogen bonding network and atomic interactions with the inhibitor. The comprehensive analysis presented in this report can provide useful information for understanding the drug resistance mechanism against lamivudine. The results can also provide some potential clues for further design of novel inhibitors that are less susceptible to drug resistance.