Search

We present a correlation-driven molecular dynamics (CDMD) method for automated refinement of atomistic models into cryo-electron microscopy (cryo-EM) maps at resolutions ranging from near-atomic to subnanometer. It utilizes a chemically accurate force field and thermodynamic sampling to improve the real-space correlation between the modeled structure and the cryo-EM map. Our framework employs a gradual increase in resolution and map-model agreement as well as simulated annealing, and allows fully automated refinement without manual intervention or any additional rotamer- and backbone-specific restraints. Using multiple challenging systems covering a wide range of map resolutions, system sizes, starting model geometries and distances from the target state, we assess the quality of generated models in terms of both model accuracy and potential of overfitting. To provide an objective comparison, we apply several well-established methods across all examples and demonstrate that CDMD performs best in most cases.

An important and computationally demanding part of molecular dynamics simulations is the calculation of long-range electrostatic interactions. Today, the prevalent method to compute these interactions is particle mesh Ewald (PME). The PME implementation in the GROMACS molecular dynamics package is extremely fast on individual GPU nodes. However, for large scale multinode parallel simulations, PME becomes the main scaling bottleneck as it requires all-to-all communication between the nodes; as a consequence, the number of exchanged messages scales quadratically with the number of involved nodes in that communication step. To enable efficient and scalable biomolecular simulations on future exascale supercomputers, clearly a method with a better scaling property is required. The fast multipole method (FMM) is such a method. As a first step on the path to exascale, we have implemented a performance-optimized, highly efficient GPU FMM and integrated it into GROMACS as an alternative to PME. For a fair performance comparison between FMM and PME, we first assessed the accuracies of the methods for various sets of input parameters. With parameters yielding similar accuracies for both methods, we determined the performance of GROMACS with FMM and compared it to PME for exemplary benchmark systems. We found that FMM with a multipole order of 8 yields electrostatic forces that are as accurate as PME with standard parameters. Further, for typical mixed-precision simulation settings, FMM does not lead to an increased energy drift with multipole orders of 8 or larger. Whereas an ≈50 000 atom simulation system with our FMM reaches only about a third of the performance with PME, for systems with large dimensions and inhomogeneous particle distribution, e.g., aerosol systems with water droplets floating in a vacuum, FMM substantially outperforms PME already on a single node.

A Flexible, GPU - Powered Fast Multipole Method for Realistic Biomolecular Simulations in Gromacs

Ion channels play a fundamental role in maintaining vital electrochemical gradients across the cell membrane and in enabling electrical signaling in cells. Understanding their functional mechanism is crucial for facilitating drug design on this important class of membrane proteins. Key characteristics of ion channel function that are commonly quantified experimentally are ionic permeation rates and selectivities. The Computational Electrophysiology (CompEL) protocol allows the investigation of ion channels in GROMACS all-atom molecular dynamics simulations. By employing an ion/water exchange protocol in a double-membrane simulation setup, a steady state with a continuous flow of ions through the channels is achieved. The recorded ion permeation events directly reveal conductance, selectivity, and rectification behavior, which all allow a straightforward comparison to experimental quantities. In addition, the ionic pathways through the channels are readily revealed. CompEL is available in GROMACS and has a negligible impact on the simulation performance even in parallel.Recent application to the bacterial channel PorB resulted in conductance and selectivity values and revealed the pathways of traversing ions. For a mutant, pathway disruption was observed, thereby explaining its resistance to certain antibiotics. The conductance calculated from the simulations agreed very well with measured conductances. Other successful examples of CompEL include the discovery of an unexpected ionic pathway in the anti-microbial peptide dermcidin as well as the identification of a novel conduction mechanism in the potassium channel KcsA. We finally investigated the conductance and the gating mechanism of pentameric ligand-gated ion channels such as GLIC. We constructed a reaction coordinate from two crystal structures, which have been proposed to represent the open and closed states of the channel, respectively. Constraining the pore to specific positions on that reaction coordinate allowed us to sample the conductivity at these positions. CompEL confirms that the assumedly open structure is indeed conducting and cation-selective (with conductances in the pS range, comparable to published experimental values), whereas the closed structure is indeed non-conducting. Moreover, conductance values steadily increase along the closed-to-open transition and are correlated with the amount of water in the pore.

Ion channels are of universal importance for all cell types and play key roles in cellular physiology and pathology. Increased insight into their functional mechanisms is crucial to enable drug design on this important class of membrane proteins, and to enhance our understanding of some of the fundamental features of cells. This review presents the concepts behind the recently developed simulation protocol Computational Electrophysiology (CompEL), which facilitates the atomistic simulation of ion channels in action. In addition, the review provides guidelines for its application in conjunction with the molecular dynamics software package GROMACS. We first lay out the rationale for designing CompEL as a method that models the driving force for ion permeation through channels the way it is established in cells, i.e., by electrochemical ion gradients across the membrane. This is followed by an outline of its implementation and a description of key settings and parameters helpful to users wishing to set up and conduct such simulations. In recent years, key mechanistic and biophysical insights have been obtained by employing the CompEL protocol to address a wide range of questions on ion channels and permeation. We summarize these recent findings on membrane proteins, which span a spectrum from highly ion-selective, narrow channels to wide diffusion pores. Finally we discuss the future potential of CompEL in light of its limitations and strengths. This article is part of a Special Issue entitled: Membrane Proteins edited by J.C. Gumbart and Sergei Noskov.

Molecular electrostatics, notably in proteins and other macromolecules, is complicated by the presence of titratable sites which occur in different forms whose charge distribution differs, e.g., due to protonation or reduction in response to changes in their environment. This variability is often crucial for molecular function and interaction properties and thus has to be included for a realistic description of electrostatics. The computation of electrostatic interactions is also the computationally most demanding part of a molecular dynamics (MD) simulation. To address these issues, we combine a fast multipole method (FMM) with λ-dynamics for the open source molecular dynamics package GROMACS.λ-dynamics bridges discrete, physical site forms via continuous λ-variables to allow the variable sites to interconvert between their forms, thus adding the variability of the charge distribution as physical detail missing in standard MD. The FMM ideally complements by enabling efficient, non-redundant computation of the electrostatic forces needed also in standard MD and of the additional electrostatic interaction energies required to compute the forces that drive the interconversion between the different site forms. In addition, the FMM allows us to take full advantage of highly parallel supercomputers beyond 10,000s of CPU cores, as predicted for the upcoming exascale machines, by avoiding time consuming communication between compute nodes.Our λ-dynamics method extends existing constant-pH λ-dynamics beyond protonation by allowing for titratable sites with any number of forms of potentially entirely different chemical structure, including binding of different ligand types. Allowing for many different site forms, and thus λ-variables, increases the size of the configuration space that needs be sampled. While interconversion between the site forms via continuous λ-variables increases sampling efficiency, it also leads to the occurrence of unphysical intermediate states. A new bias potential that explicitly controls the occurrence of unphysical intermediate states avoids their overrepresentation and ensures sampling of physically relevant configurations.

Most plasmalemmal proteins organize in submicrometer-sized clusters whose architecture and dynamics are still enigmatic. With syntaxin 1 as an example, we applied a combination of far-field optical nanoscopy, biochemistry, fluorescence recovery after photobleaching (FRAP) analysis, and simulations to show that clustering can be explained by self-organization based on simple physical principles. On average, the syntaxin clusters exhibit a diameter of 50 to 60 nanometers and contain 75 densely crowded syntaxins that dynamically exchange with freely diffusing molecules. Self-association depends on weak homophilic protein-protein interactions. Simulations suggest that clustering immobilizes and conformationally constrains the molecules. Moreover, a balance between self-association and crowding-induced steric repulsions is sufficient to explain both the size and dynamics of syntaxin clusters and likely of many oligomerizing membrane proteins that form supramolecular structures.

Hybrid state formation is a crucial step of the tRNA translocation process. After peptide bond formation, a deacylated tRNA in the P site (here tRNA-fMet) moves to a hybrid P/E configuration while the A-site tRNA (here tRNA-Val) still occupies a classical A/A state. In this process the CCA tail of the tRNA-fMet must move by about 4.3 nm and the elbow region moves by ∼2.5 nm forming contact with the L1 stalk. The anticodon stem-loop remains almost fixed at its 30S position, the 30S body and head undergo a slight (∼5°) rotation relative to their classical state. The atomic detailed pathway of this transition is still unresolved. Here, starting from structures obtained from previous equilibrium simulations, we performed a series of 200ns all-atom MD simulations in which we drive the system using a biasing potential that maximizes correlation to the target state ensemble. Analysis of the trajectories in terms of the applied biasing potential and the dwelling time in configurational space appears to indicate the existence of several intermediate states separated by barriers along the hybrid-state formation pathway. Interestingly, most of the bottlenecks along the trajectories are found to involve interactions between the CCA tail of the tRNA-fMet and specific 50S residues. The CCA tail interacts initially with the P-loop of helix 80 of the 23S RNA, it then contacts helix 93 and U2076 of helix 75, finally reaching the E site through a different set of interactions which include His19 of ribosomal protein L28. In contrast, the other large scale motions of the L1 stalk and 30S subunit appear to be less restricted along the pathway.

In light of recent high-resolution cryo-electron-microscopy (cryo-EM) structures of large complexes, modeling atom coordinates into cryo-EM densities faces new challenges. Though high-resolution structures provide more data, finding the best-fitting structure with current refinement methods is more difficult with increased resolution, because the used refinement potentials become more rugged the higher the resolution. Currently used refinement potentials are defined by empirically chosen measures of similarity between a calculated cryo-EM density and the given experimental map, e.g. cross-correlation or absolute distance.Here, we present a new refinement potential that is based on a statistical physics model of the cryo-EM measuring and reconstruction process using Bayesian statistics. Our method contains previously previously developed algorithms as limiting cases.The minima of the refinement potential and its shape both influence the efficiency of the refinement algorithm; a smoother energy landscape allows a more efficient exploration of the minima, i.e. fitting structures, in this landscape. Compared to earlier methods, our refinement energy landscape is smoother, allowing more efficient sampling of the energy landscape. Further, our refinement protocol provides an appropriate refinement force constant and takes into account the thermal fluctuation of the atoms. Additionally, our algorithm allows us to generate molecular dynamics ensembles that represent the simultaneous input from multiple cryo-EM maps.

Near-atom resolution in 3D cryo-electron microscopy brings two new challenges for refining atom coordinates to densities. Because their underlying potentials become more rugged with increased resolution, established refinement methods are easily trapped in local minima. Second, single structures no longer well describe the inherent structural diversity that can be resolved with cryo-EM.Here, we developed a method to derive cryo-EM refinement potentials from a Bayesian approach. Our refinement potential takes a physical model of the cryo-EM measuring and reconstruction process into account using Bayesian statistics. The result is a potential that statistically correctly reflects the given EM-data and is smooth even at high resolution. Previously developed algorithms are contained as limiting cases.With our method we are able to represent the configurational dynamics that is captured in cryo-EM density maps through a series of features. The smoothness of our refinement energy landscape allows efficient sampling and refinement while also taking into account thermal fluctuations. We provide an refinement force constant and potential from the Bayes approach that truthfully represents the complete distribution of atom configurations underlying the given EM maps. Thus, in combination with traditional molecular dynamics simulation we create refined ensembles that - as a whole - represent a given cryo-EM map.We further use the advantage of the Bayes approach to generate molecular dynamics ensembles that represent the simultaneous input from multiple cryo-EM maps together to capture the physically relevant transitions between different states as resolved by cryo-EM. Overall, our methods enables complete use of the data in cryo-EM maps and provides structural interpretation from ensembles that will aid understanding ever more complex cryo-EM data.

The movement of tRNAs during translocation is accompanied by large conformational changes of the ribosome such as intersubunit rotations. Here, we present a method to cover conformational changes of the ribosome that occur on timescales not accessible to equilibrium Molecular Dynamics (MD) simulations through the combination of X-ray crystallography, cryo-EM data and MD simulations. Cryo-electron microscopy (Cryo-EM) provides medium/low-resolution density maps for many intermediate states of large molecular complexes. In contrast, X-ray crystallography provides high-resolution structures, usually limited to few states. To obtain pathways connecting the intermediate states of tRNA translocation, we start MD simulations from crystal structures with an additional biasing potential. This biasing potential, which maximizes the correlation between atomic model and cryo-EM map (Tama 2008), allows us to drive the ribosome from one intermediate state to another, covering the full translocation pathway. This method of cryo-EM driven MD was implemented in the high-throughput and highly parallel MD simulation package GROMACS.

Presently, simulations of ion channel permeation usually rely on nonatomistic Brownian dynamics calculations, indirect interpretation of energy maps or external electric fields. In biological cells, however, electric fields across bilayers are generally established by electrochemical ionic gradients. Small charge imbalances across the bilayers also evoke electric fields applied in electrophysiological experiments on a microscopic scale.We present a computational method which enables the direct simulation of ion flux through membrane channels driven by biologically realistic electrochemical gradients. This also makes the simulation of reversal potential experiments possible. As it is implemented in a highly efficient molecular dynamics package, simulation timescales relevant for physiology and experiment, for instance single-channel electrophysiology, are achieved.We illustrate the use of our method by applying it to the bacterial channel PorB from pathogenic Neisseria meningitidis. PorB inserts into the inner mitochondrial membrane of target cells during Neisserial infection and triggers their apoptosis by dissipating the potential across the membrane. These channels must also be passable for antibiotics during treatment. We show that the method accurately predicts ion conductance and selectivity and elucidates detailed ion conduction mechanisms. PorB mutants resistant to antibiotics display a markedly altered ion pathway and selectivity.

Neurotransmitter release at the synapse requires fusion of synaptic vesicles with the presynaptic plasma membrane. SNAREs are the core constituents of the protein machinery responsible for this membrane fusion, but the actual fusion mechanism remains unclear. Here, we have simulated neuronal SNARE-mediated membrane fusion in molecular detail. In our simulations, membrane fusion progresses through an inverted micelle fusion intermediate before reaching the hemifused state. We show that at least one single SNARE complex is required for fusion, as has also been confirmed in a recent in vitro single-molecule fluoresence study. Further, the transmembrane regions of the SNAREs were found to play a vital role in the initiation of fusion by causing distortions of the lipid packing of the outer membrane leaflets, and the C termini of the transmembrane regions are associated with the formation of the fusion pores. The inherent mechanical stress in the linker region of the SNARE complex was found to drive both the subsequent formation and expansion of fusion pores. Our simulations also revealed that the presence of homodimerizations between the transmembrane regions leads to the formation of unstable fusion intermediates that are under high curvature stress. We show that multiple SNARE complexes mediate membrane fusion in a cooperative and synchronized process. Finally, we show that after fusion, the zipping of the SNAREs extends into the membrane region, in agreement with the recently resolved X-ray structure of the fully assembled state.

Presently, most simulations of ion channel function rely upon nonatomistic Brownian dynamics calculations, indirect interpretation of energy maps, or application of external electric fields. We present a computational method to directly simulate ion flux through membrane channels based on biologically realistic electrochemical gradients. In close analogy to single-channel electrophysiology, physiologically and experimentally relevant timescales are achieved. We apply our method to the bacterial channel PorB from pathogenic Neisseria meningitidis, which, during Neisserial infection, inserts into the mitochondrial membrane of target cells and elicits apoptosis by dissipating the membrane potential. We show that our method accurately predicts ion conductance and selectivity and elucidates ion conduction mechanisms in great detail. Handles for overcoming channel-related antibiotic resistance are identified.

We introduce a computational toolset, named GROmaρs, to obtain and compare time-averaged density maps from molecular dynamics simulations. GROmaρs efficiently computes density maps by fast multi-Gaussian spreading of atomic densities onto a three-dimensional grid. It complements existing map-based tools by enabling spatial inspection of atomic average localization during the simulations. Most importantly, it allows the comparison between computed and reference maps (e.g., experimental) through calculation of difference maps and local and time-resolved global correlation. These comparison operations proved useful to quantitatively contrast perturbed and control simulation data sets and to examine how much biomolecular systems resemble both synthetic and experimental density maps. This was especially advantageous for multimolecule systems in which standard comparisons like RMSDs are difficult to compute. In addition, GROmaρs incorporates absolute and relative spatial free-energy estimates to provide an energetic picture of atomistic localization. This is an open-source GROMACS-based toolset, thus allowing for static or dynamic selection of atoms or even coarse-grained beads for the density calculation. Furthermore, masking of regions was implemented to speed up calculations and to facilitate the comparison with experimental maps. Beyond map comparison, GROmaρs provides a straightforward method to detect solvent cavities and average charge distribution in biomolecular systems. We employed all these functionalities to inspect the localization of lipid and water molecules in aquaporin systems, the binding of cholesterol to the G protein coupled chemokine receptor type 4, and the identification of permeation pathways through the dermicidin antimicrobial channel. Based on these examples, we anticipate a high applicability of GROmaρs for the analysis of molecular dynamics simulations and their comparison with experimentally determined densities.

We describe a versatile method to enforce the rotation of subsets of atoms, e.g., a protein subunit, in molecular dynamics (MD) simulations. In particular, we introduce a "flexible axis" technique that allows realistic flexible adaptions of both the rotary subunit as well as the local rotation axis during the simulation. A variety of useful rotation potentials were implemented for the GROMACS 4.5 MD package. Application to the molecular motor F(1)-ATP synthase demonstrates the advantages of the flexible axis approach over the established fixed axis rotation technique.