Understanding and controlling chemical reactivity in biological systems requires atomic-level insight into processes that are often inaccessible to experiments. Hybrid quantum mechanics/molecular mechanics (QM/MM) simulations provide a powerful framework to describe chemical reactions in complex environments, but their practical application remains limited by fragmented software ecosystems, restricted accessibility, and methodological approximations that can compromise accuracy and reproducibility. In particular, many existing QM/MM implementations rely on ad hoc couplings, proprietary software, or truncated treatments of long-range electrostatic interactions. Here, we present a robust and fully periodic QM/MM interface between the open-source molecular dynamics engine GROMACS and the electronic structure theory code CP2K. The implementation enables efficient and reproducible QM/MM molecular dynamics and enhanced sampling simulations with a consistent treatment of long-range electrostatics under periodic boundary conditions. By combining the strengths of two widely used community codes, this interface provides a general and scalable platform for studying

Better detectors and automated data collection have generated a flood of high-resolution cryo-EM maps, which in turn has renewed interest in improving methods for determining structure models corresponding to these maps. However, automatically fitting atoms to densities becomes difficult as their resolution increases and the refinement potential has a vast number of local minima. In practice, the problem becomes even more complex when one also wants to achieve a balance between a good fit of atom positions to the map, while also establishing good stereochemistry or allowing protein secondary structure to change during fitting. Here, we present a solution to this challenge using a maximum likelihood approach by formulating the problem as identifying the structure most likely to have produced the observed density map. This allows us to derive new types of smooth refinement potential—based on relative entropy—in combination with a novel adaptive force scaling algorithm to allow balancing of force-field and density-based potentials. In a low-noise scenario, as expected from modern cryo-EM data, the relative-entropy based refinement potential outperforms alternatives, and the adaptive force scaling appears to aid all existing refinement potentials. The method is available as a component in the GROMACS molecular simulation toolkit.

The resolution revolution has increasingly enabled single-particle cryogenic electron microscopy (cryo-EM) reconstructions of previously inaccessible systems, including membrane proteins—a category that constitutes a disproportionate share of drug targets. We present a protocol for using density-guided molecular dynamics simulations to automatically refine atomistic models into membrane protein cryo-EM maps. Using adaptive force density-guided simulations as implemented in the GROMACS molecular dynamics package, we show how automated model refinement of a membrane protein is achieved without the need to manually tune the fitting force ad hoc. We also present selection criteria to choose the best-fit model that balances stereochemistry and goodness of fit. The proposed protocol was used to refine models into a new cryo-EM density of the membrane protein maltoporin, either in a lipid bilayer or detergent micelle, and we found that results do not substantially differ from fitting in solution. Fitted structures satisfied classical model-quality metrics and improved the quality and the model-to-map correlation of the x-ray starting structure. Additionally, the density-guided fitting in combination with generalized orientation-dependent all-atom potential was used to correct the pixel-size estimation of the experimental cryo-EM density map. This work demonstrates the applicability of a straightforward automated approach to fitting membrane protein cryo-EM densities. Such computational approaches promise to facilitate rapid refinement of proteins under different conditions or with various ligands present, including targets in the highly relevant superfamily of membrane proteins.

Better detectors and automated data collection have generated a flood of high-resolution cryo-EM maps, which in turn has renewed interest in improving methods for determining structure models corresponding to these maps. However, automatically fitting atoms to densities becomes difficult as their resolution increases and the refinement potential has a vast number of local minima. In practice, the problem becomes even more complex when one also wants to achieve a balance between a good fit of atom positions to the map, while also establishing good stereochemistry or allowing protein secondary structure to change during fitting. Here, we present a solution to this challenge using Bayes’ approach by formulating the problem as identifying the structure most likely to have produced the observed density map. This allows us to derive a new type of smooth refinement potential - based on relative entropy - in combination with a novel adaptive force scaling algorithm to allow balancing of force-field and density-based potentials. In a low-noise scenario, as expected from modern cryo-EM data, the Bayesian refinement potential outperforms alternatives, and the adaptive force scaling appears to also aid existing refinement potentials. The method is available as a component in the GROMACS molecular simulation toolkit.

Better detectors and automated data collection has created a flood of high-resolution cryo-EM density maps, which in turn has renewed interest in improving methods for determining atomic structure models corresponding to these maps. Automatically fitting atoms to a density is increasingly difficult as the resolution increases and the refinement potential has a large number of local minima. In practice the problem is even more complex when one also wants to balance fit to maps with maintaining good stereochemistry and allowing protein secondary structure to change as part of the fitting.

Ligand-gated ion channels are critical mediators of electrochemical signal transduction across evolution. Biophysical and pharmacological development in this family relies on high-quality structural data in multiple, subtly distinct functional states; however, structural data remain limited, particularly for the native closed state. Here, we report cryo-electron microscopy structures of the proton-gated Gloeobacter violaceus ligand-gated ion channel (GLIC) under resting and activating conditions (pH 7, 5, and 3). Predominant classes in all three conditions featured constricted pores, consistent with a low maximal open probability; nonetheless, densities in predicted gating regions were better defined at low pH, implicating key electrostatic networks in channel activation. Molecular dynamics simulations further reflected domain stabilization at lower pH, and greater alignment with previous X-ray structures. Further analysis of datasets collected in each condition enabled the reconstruction of minority classes, including an alternative low-pH state with an expanded pore. In addition to providing new structures of closed GLIC in multiple conditions without crystallization, these results offer dynamic insight into an ion channel's heterogeneous resting state, with activating conditions condensing the energy landscape on a pathway towards gating.

The resolution revolution has increasingly enabled single-particle cryo-electron microscopy (cryo-EM) reconstructions of previously inaccessible systems, including membrane proteins that constitute a disproportionate share of drug targets. However, manual model-building remains a barrier to efficient structure determination, often involving considerable time, expertise, and potential user bias. We present a protocol for using density-guided molecular dynamics simulations to automatically refine atomistic models into membrane-protein cryo-EM maps. Previous simulations-based refinement methods have often required additional user input, such as secondary-structure restraints, to avoid model distortion, because they weigh refinement forces against stereochemical features ad hoc. Using adaptive-force density-guided simulations as implemented in the GROMACS molecular dynamics package, we show how automated model refinement of membrane proteins is achieved while dynamically balancing stereochemistry and goodness-of-fit. We applied this protocol to refine models of multiple <350-kD membrane proteins, both in lipid bilayers and detergent micelles. Fitted structures satisfied classical model-quality metrics and in many cases rivalled manually built models, while taking hours instead of days to produce. Furthermore, comparison of fits using alternative parameters demonstrated the influence of solubilization conditions and model-map similarity metrics. This work demonstrates the applicability of a straightforward automated approach to fitting membrane-protein cryo-EM densities. Such computational approaches promise to facilitate determination of complex or challenging structures, including targets in the highly relevant superfamily of membrane proteins.

Data accompanying prepared manuscript to describe novel density-guided simulation algorithms.

Pentameric ligand-gated ion channels (pLGICs) play central roles for signal conduction in the nervous system, where closely related anionic or cationic channels exhibit a broad range responses to various neurotransmitters. In addition to the primary agonist, pLGICs are also highly sensitive to allosteric modulation, which makes them highly interesting model systems to understand conformational transitions and state-specific stabilisation - in particular bacterial homologs such as the pH-gated GLIC channel from Gloeobacter violaceus. It has however been surprisingly difficult to characterise both the various states as well as the conformational transitions between them. Here, we present a series of new high-resolution cryo-EM structures of the GLIC channel obtained at a range of different pH values, with local resolutions of 3-5Å. These structures show interesting differences to previous X-ray structures, including a significant bias towards non-conducting transmembrane pores at all pH values, while the overall orientation of the pore-lining helices and ECD vary more between structures. To improve the structural modeling, we combined the 3D density reconstruction with a newly developed framework for Bayesian density fitting of models in molecular dynamics simulations, which enables us to rapidly achieve high-resolution structural models, to determine pH-dependent structural differences, and not least to assess the flexibility of different parts of the structure. This provides new insight into the structure of ligand gated ion channels under different conditions, it raises interesting questions about the exact local pH environment on cryo-EM grids, and it provides promising results about our ability to automatically fit channel models to electron densities of different resolution.

Pentameric ligand gated ion channels (pLGICs) are important in neuronal communication between synapses. Endogenous ligands are known to mediate both the opening and closing of pLGICs and these proteins are also the target for a number of neurotoxins, animal venoms, as well as pharmacological agents. In the past decade a number of crystal structures have elucidated gating mechanisms as well as open, closed, and desensitized states, often in the presence of agonists or antagonists. More recently, cryo-electron microscopy (cryo-EM) has also produced high resolution structures of channels in the pLGIC family and has shed further light on ligand mediated activation. Given the importance of these proteins to pharmaceutical intervention, a method to understand how ligands bind to pLGICs could have important clinical impacts. We have recently implemented a tool that allows inclusion of cryo-EM densities as an input to the popular molecular dynamics package GROMACS. Using a fitting procedure that uses successively higher resolution density maps derived from experimental cryo-EM data, we demonstrate the ability to drive a simulation to a conformation corresponding to EM density. Here we use recently solved cryo-EM structures of members of the pLGIC family to demonstrate the use of this procedure to give insight into modes of ligand binding in this important class of molecules. The atomic structures derived from this method are then compared to the models built directly from the EM density as well as to available crystal structures. We believe this novel technique has much promise to help elucidate mechanisms of drug binding in the increasing number of systems that have been characterized in a ligand bound state via cryo-EM.

Given the need for modern researchers to produce open, reproducible scientific output, the lack of standards and best practices for sharing data and workflows used to produce and analyze molecular dynamics (MD) simulations have become an important issue in the field. There are now multiple well-established packages to perform molecular dynamics simulations, often highly tuned for exploiting specific classes of hardware, and each with strong communities surrounding them, but with very limited interoperability/transferability options. Thus, the choice of the software package often dictates the workflow for both simulation production and analysis. The level of detail in documenting the workflows and analysis code varies greatly in published work, hindering reproducibility of the reported results and the ability for other researchers to build on these studies. An increasing number of researchers are motivated to make their data available, but many challenges remain in order to effectively share and reuse simulation data. To discuss these and other issues related to best practices in the field in general, we organized a workshop in November 2018 ( https://bioexcel.eu/events/workshop-on-sharing-data-from-molecular-simulations/). Here, we present a brief overview of this workshop and topics discussed. We hope this effort will spark further conversation in the MD community to pave the way towards more open, interoperable and reproducible outputs coming from research studies using MD simulations.

A free energy landscape estimation method based on the well-known Gaussian mixture model (GMM) is used to compare the efficiencies of thermally enhanced sampling methods with respect to regular molecular dynamics. The simulations are carried out on two binding states of calmodulin, and the free energy estimation method is compared with other estimators using a toy model. We show that GMM with cross-validation provides a robust estimate that is not subject to overfitting. The continuous nature of Gaussians provides better estimates on sparse data than canonical histogramming. We find that diffusion properties determine the sampling method effectiveness, such that diffusion-dominated apo calmodulin is most efficiently sampled by regular molecular dynamics, while holo calmodulin, with its rugged free energy landscape, is better sampled by enhanced sampling methods.

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.

After the resolution revolution, cryo-EM moves from sets of individual structures, to understanding properties of biological samples as a whole, like pathways between resolved structures and binding affinities. Extracting these molecule mechanics as probability distribution of arbitrary observables, like open and close configuation, or small molecule binding probabilities however, is still challenging despite increasing resolution and numbers of resolved structures of a single complex, due to the non-trivial connection three dimensional structures and all the possible molecule conformations they represent. Here, we present how to resolve all molecule conformations that are represented by three dimensional structures and pathways between multiple reconstructed densities. The opening and closing dynamics of adenylate kinase are used to verify our the method, using ensemble-averaged density maps of open and closed state. In the ribosome-selB complex, our method further refines structural transitions for GTPase activation. In our work, the free energy for each molecule conformation under a cryo-EM experiment is derived from a probabilistic formulation of the complete cryo-EM measurement process. It is shown why a set of reconstructed densities is a good approximation to the rigorous but computationally challenging calculation of conformation weights from all cryo-EM images. It is found that morphing densities and refining structures into cryo-EM maps is tightly coupled through the measure of goodness-of-fit. With these insights, we asess how much of all underlying configurations are represented by single densities and sample pathways between densities, that represent physical intermediate states of the ribosome-selB complex upon GTPase activation.

A free energy landscape estimation-method based on Bayesian inference is presented and used for comparing the efficiency of thermally enhanced sampling methods with respect to regular molecular dynamics, where the simulations are carried out on two binding states of calmodulin. The proposed free energy estimation method (the GM method) is compared to other estimators using a toy model showing that the GM method provides a robust estimate not subject to overfitting. The continuous nature of the GM method, as well as predictive inference on the number of basis functions, provide better estimates on sparse data. We find that the free energy diffusion proper- ties determine sampling method effectiveness, such that the diffusion dominated apo-calmodulin is most efficiently sampled by regular molecular dynamics, while the holo with its rugged free energy landscape is better sampled by enhanced methods.

Beneath any density map resolved by cryo electron microscopy (cryo-EM) lies an ensemble of all-atom structures. However, established refinement methods are designed with a single best-fitting structures in mind. Further, established refinement protocols require intricate sampling schemes and/or additional constraints on, e.g., secondary structure due to the ruggedness of their refinement potentials. Here, we introduce a goodness-of-fit measure of all-atom model to density that reflects the specific physics of cryo-EM in contrast to general measures like cross-correlation or density sum at atom positions. Its smoothness allows for all-atom refinement against cryo-EM densities without any further constraints on, e.g., secondary structure, and capturing large conformational transitions between initial model and target density. Bayesian cryo-EM ensemble refinement was tested for refinement of closed conformation adenylate-kinase (AKE) against a target density created from an ensemble of structures in the open conformation. Our ensemble is closer to the underlying ensemble as measured by Kullback-Leibler divergence after dimensionality reduction by 1.1 bits compared to a cross-correlation based potential, and 1.5 bits compared to inverted-density potential. Bayesian refinement on novel cryo-EM data yielded excellent agreement in integrated Fourier shell correlation curves (>0.5 for a 6 Å density map). Obliterating further constraints, our refinement protocol enables secondary structure prediction while providing an ensemble of high-quality structures. Overall, our method yields reliable refinement results for challenging cases using a Bayesian potential and structure ensembles to describe cryo-EM densities.

During protein synthesis, tRNA molecules move from the ribosome's aminoacyl to peptidyl to exit sites, with the two ribosomal subunits remaining associated through intersubunit bridges, despite rapid large-scale intersubunit rotation. Using molecular dynamics simulations, we here investigate conformational motions during spontaneous translocation, as well as the underlying energetics and kinetics. Resolving fast transitions between states, we find that tRNA motions govern the transition rates within the pre- and post-translocation states. The L1 stalk drives the tRNA from the peptidyl site and links intersubunit rotation to translocation. Displacement of tRNAs is controlled by ‘sliding’ and ‘stepping’ mechanisms involving conserved L6, L5, and L1 residues, thus ensuring binding to the ribosome despite large-scale tRNA movement through maintaining constant binding affinity. Intersubunit rotations exhibit remarkably fast intrinsic submicrosecond dynamics, which requires a fine-tuned flat free energy landscape, as any larger barrier would slow down the conformational motions. Maintaining such subtle balance between the many interactions involved is remarkable, in particular considering the large shifts the many intersubunit bridges undergo. Based on the observed occupancies of intersubunit contacts during our simulations, peripheral clusters were found to maintain strong steady interactions by changing contacts in the course of rotation. The peripheral B1 bridges are stabilized by a changing contact pattern of charged residues that adapts to the rotational state. In contrast, steady strong interactions of the B4 bridge are ensured by the flexible helix H34 following the movement of protein S15. The total contribution of the tRNAs -- which contact both subunits -- to the intersubunit binding enthalpy is almost constant, despite their different positions in the ribosome. These mechanisms keep the intersubunit interaction strong and steady during rotation, thereby preventing dissociation and enabling rapid rotation.

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.