Real Cost of Speed: The Effect of a Time-Saving Multiple-Time-Stepping Algorithm on the Accuracy of Molecular Dynamics Simulations

On the Validation of Protein Force Fields Based on Structural Criteria

Molecular dynamics simulations depend critically on the quality of the force field used to describe the interatomic interactions and the extent to which it has been validated for use in a specific application. Using a curated test set of 52 high-resolution structures, 39 derived from X-ray diffraction and 13 solved using NMR, we consider the extent to which different parameter sets of the GROMOS protein force field can be distinguished based on comparing a range of structural criteria, including the number of backbone hydrogen bonds, the number of native hydrogen bonds, polar and nonpolar solvent-accessible surface area, radius of gyration, the prevalence of secondary structure elements, J-coupling constants, nuclear Overhauser effect (NOE) intensities, positional root-mean-square deviations (RMSD), and the distribution of backbone ϕ and ψ dihedral angles. It is shown that while statistically significant differences between the average values of individual metrics could be detected, these were in general small. Furthermore, improvements in agreement in one metric were often offset by loss of agreement in another. The work establishes a framework and test set against which protein force fields can be validated. It also highlights the danger of inferring the relative quality of a given force field based on a small range of structural properties or small number of proteins.

PyThinFilm: Automated Molecular Dynamics Simulation Protocols for the Generation of Thin Film Morphologies

The performance of organic optoelectronic devices, such as organic light-emitting diodes (OLEDs) and organic solar cells (OSCs), is intrinsically related to the molecular-scale morphology of the thin films from which they are composed. However, the experimental characterization of morphology at the molecular level is challenging due to the often amorphous or at best semicrystalline nature of these films. Classical molecular modeling techniques, such as molecular dynamics (MD) simulation, are increasingly used to understand the relationship between morphology and the properties of thin-film devices. PyThinFilm (github.com/ATB-UQ/PyThinFilm) is an open-source Python package which allows fully automated MD simulations of thin film growth to be performed using vacuum and/or solution deposition processes. PyThinFilm utilizes the GROMACS simulation package in combination with interaction parameters from the Automated Topology Builder (atb.uq.edu.au). Here, PyThinFilm is described along with an overview of applications in which PyThinFilm has been used to study the thin films of organic semiconductor materials typically used in OLEDs and OSCs.

On the use of multiple-time-step algorithms to save computing effort in molecular dynamics simulations of proteins

Computer simulation of proteins in aqueous solution at the atomic level of resolution is still limited in time span and system size due to limited computing power available and thus employs a variety of time-saving techniques that trade some accuracy against computational effort. Examples of such time-saving techniques are the application of constraints to particular degrees of freedom or the use of a multiple-time-step (MTS) algorithm distinguishing between particular forces when integrating Newton's equations of motion. The application of two types of MTS algorithms to bond-stretching forces versus the remaining forces in molecular dynamics (MD) simulations of a protein in aqueous solution or of liquid water is investigated and the results in terms of total energy conservation and the influence on various other properties are compared to those of MD simulations of the same systems using bond-length, and for water bond-angle, constraints. At comparable computational effort, the use of bond-length constraints in proteins leads to better energy conservation and less distorted properties than the two MTS algorithms investigated.

On the Effect of the Various Assumptions and Approximations used in Molecular Simulations on the Properties of Bio-Molecular Systems: Overview and Perspective on Issues

Computer simulations of molecular systems enable structure-energy-function relationships of molecular processes to be described at the sub-atomic, atomic, supra-atomic or supra-molecular level and plays an increasingly important role in chemistry, biology and physics. To interpret the results of such simulations appropriately, the degree of uncertainty and potential errors affecting the calculated properties must be considered. Uncertainty and errors arise from (1) assumptions underlying the molecular model, force field and simulation algorithms, (2) approximations implicit in the interatomic interaction function (force field), or when integrating the equations of motion, (3) the chosen values of the parameters that determine the accuracy of the approximations used, and (4) the nature of the system and the property of interest. In this overview, advantages and shortcomings of assumptions and approximations commonly used when simulating bio-molecular systems are considered. What the developers of bio-molecular force fields and simulation software can do to facilitate and broaden research involving bio-molecular simulations is also discussed.

The effect of different cutoff schemes in molecular simulations of proteins

Molecular simulations of nanoscale systems invariably involve assumptions and approximations to describe the electrostatic interactions, which are long-ranged in nature. One approach is the use of cutoff schemes with a reaction-field contribution to account for the medium outside the cutoff scheme. Recent reports show that macroscopic properties may depend on the exact choice of cutoff schemes in modern day simulations. In this work, a systematic analysis of the effects of different cutoff schemes was performed using a set of 52 proteins. We find no statistically significant differences between using a twin-range or a single-range cutoff scheme. Applying the cutoff based on charge groups or based on atomic positions, does lead to significant differences, which is traced to the cutoff noise for energies and forces. While group-based cutoff schemes show increased cutoff noise in the potential energy, applying an atomistic cutoff leads to artificial structure in the solvent at the cutoff distance. Carefully setting the temperature control, or using an atomistic cutoff for the solute and a group-based cutoff for the solvent significantly reduces the effects of the cutoff noise, without introducing structure in the solvent. This study aims to deepen the understanding of the implications different cutoffs have on molecular dynamics simulations.

Understanding the Activated Form of a Class-I Fusion Protein: Modeling the Interaction of the Ebola Virus Glycoprotein 2 with a Lipid Bilayer

The fusion of the viral and target cell membranes is a key step in the life cycle of all enveloped viruses. Here, a range of structural data is used to generate an evidence-based model of the active conformation of an archetypical type-I fusion protein, the Ebola glycoprotein 2 (GP2). The stability of the trimeric complex is demonstrated using molecular dynamics and validated by simulating the interaction of the complex with a lipid bilayer. In this model, the fusion peptides project away from the central helix bundle parallel to the target membrane. This maximizes contact with the host membrane, enhances lateral stability, and would explain why, when activated, viral fusion proteins are trimeric.

The Effect of Using a Twin-Range Cutoff Scheme for Nonbonded Interactions: Implications for Force-Field Parametrization?

Recently, concerns have been voiced regarding the validity of the GROMOS force fields, being parametrized using a twin-range cutoff scheme, in which longer ranged nonbonded forces and energies are updated less frequently than shorter ranged ones. Here we demonstrate that the influence of such a scheme on the thermodynamic, structural, and dynamic properties used in the parametrization of the GROMOS force fields is minor. We find root-mean-square differences of maximally 0.5 kJ/mol for the solvation free energy and heat of vaporization and of maximally 0.4% for the density. Slightly larger differences are observed when switching from a group-based to an atom-based cutoff scheme. In cases where the twin-range cutoff scheme does result in minor differences compared to a single-range cutoff these are well within the deviation from the experimentally measured values.

Evolution and Morphology of Thin Films Formed by Solvent Evaporation: An Organic Semiconductor Case Study

The crucial role played by the solution-vapor interface in determining the growth and morphology of an organic semiconductor thin film formed by solvent evaporation has been examined in atomic detail. Specifically, how the loss of individual solvent molecules from the surface of the solution induces solute assembly has been studied using molecular dynamics simulations. The system consisted of bis(2-phenylpyridine) (acetylacetonate)iridium(III) [Ir(ppy)₂(acac)] and 4,4'-bis(N-carbazolyl)-1,1'-biphenyl (CBP) in chloroform at 310 K. The simulations clearly indicate that (a) the system does not undergo uniform phase separation (spinodal decomposition), (b) solute aggregation initiates at the solution-vapor interface, (c) the distribution of solvent in the film is nonhomogeneous, (d) this nonhomogeneous distribution can induce preferential alignment of host molecules, and (e) a portion of the solvent likely remains trapped within the film. The work not only demonstrates the ability to directly model evaporation in atomic detail on the relevant length scales but also shows that atomistic simulations have the potential to shed new light on morphological properties of a wide range of organic semiconductor devices manufactured using solution-processing methods.

A Computational and Modeling Study of the Reaction Mechanism of Staphylococcus aureus Monoglycosyltransferase Reveals New Insights on the GT51 Family of Enzymes

Bacterial glycosyltransferases of the GT51 family are key enzymes in bacterial cell wall synthesis. Inhibiting cell wall synthesis is a very effective approach for development of antibiotics, as this can lead to either bacteriostatic or bactericidal effects. Even though the existence of this family has been known for over 50 years, only one potent inhibitor exists, which is an analog of the lipid IV product and derived from a natural product. Drug development focused on bacterial transglycosylase has been hampered due to little being know about its structure and reaction mechanism. In this study, Staphylococcus aureus monoglycosyltransferase was investigated at an atomistic level using computational methods. Classical molecular dynamics simulations were used to reveal information about the large-scale dynamics of the enzyme-substrate complex and the importance of magnesium in structure and function of the protein, while mixed mode quantum mechanics/molecular mechanics calculations unveiled a novel hypothesis for the reaction mechanism. From these results, we present a new model for the binding mode of lipid II and the reaction mechanism of the GT51 glycosyltransferases. A metal-bound hydroxide catalyzed reaction mechanism yields an estimated free energy barrier of 16.1 ± 1.0 kcal/mol, which is in line with experimental values. The importance of divalent cations is also further discussed. These findings could significantly aid targeted drug design, particularly the efficient development of transition state analogues as potential inhibitors for the GT51 glycosyltransferases.

Rings, Hexagons, Hetals, and Dipolar Moment Sink-Sources: The Fanciful Behavior of Water around Cyclodextrin Complexes

The basket-like geometry of cyclodextrins (CDs), with a cavity able to host hydrophobic groups, makes these molecules well suited for a large number of fundamental and industrial applications. Most of the established CD-based applications rely on trial and error studies, often ignoring key information at the atomic level that could be employed to design new products and to optimize their use. Computational simulations are well suited to fill this gap, especially in the case of CD systems due to their low number of degrees of freedom compared with typical macromolecular systems. Thus, the design and validation of solid and efficient methods to simulate and analyze CD-based systems is key to contribute to this field. The behavior of supramolecular complexes critically depends on the media where they are embedded, so the detailed characterization of the solvent is required to fully understand these systems. In the present work, we use the inclusion complex formed by two α-CDs and one sodium dodecyl sulfate molecule to test eight different parameterizations of the GROMOS and AMBER force fields, including several methods aimed to increase the conformational sampling in computational molecular dynamics simulation trajectories. The system proved to be extremely sensitive to the employed force field, as well as to the presence of a water/air interface. In agreement with previous experiments and in contrast to the results obtained with AMBER, the analysis of the simulations using GROMOS showed a quick adsorption of the complex to the interface as well as an extremely exotic behavior of the water molecules surrounding the structure both in the bulk aqueous solution and at the water surface. The chirality of the CD molecule seems to play an important role in this behavior. All together, these results are expected to be useful to better understand the behavior of CD-based supramolecular complexes such as adsorption or aggregation driving forces, as well as to introduce new methods able to speed up general MD simulations.

Molecular Dynamics Simulation of Proteins

Molecular dynamics simulations allow the conformational motion of a molecule such as a protein to be followed over time at atomic-level detail. Several choices need to be made prior to running a simulation, including the software, which molecules to include in the simulation, and the force field used to describe their behavior. Guidance on making these choices and other important aspects of running MD simulations is outlined here.

Reaction-field electrostatics in molecular dynamics simulations: development of a conservative scheme compatible with an atomic cutoff

Shifting and switching schemes are developed to enable strict energy conservation in molecular dynamics simulations relying on reaction-field electrostatic (as well as Lennard-Jones) interactions with an atom-based cutoff truncation.

Spontaneous Curvature, Differential Stress, and Bending Modulus of Asymmetric Lipid Membranes

Lipid bilayers can exhibit asymmetric states, in which the physical characteristics of one leaflet differ from those of the other. This most visibly manifests in a different lipid composition, but it can also involve opposing lateral stresses in each leaflet that combine to an overall vanishing membrane tension. Here, we use theoretical modeling and coarse-grained simulation to explore the interplay between a compositional asymmetry and a nonvanishing differential stress. Minimizing the total elastic energy leads to a preferred spontaneous curvature that balances torques due to both bending moments and differential stress, with sometimes unexpected consequences. For instance, asymmetric flat bilayers, whose specific areas in each leaflet are matched to those of corresponding tensionless symmetric flat membranes, still exhibit a residual differential stress because the conditions of vanishing area strain and vanishing bending moment differ. We also measure the curvature rigidity of asymmetric bilayers and find that a sufficiently strong differential stress, but not compositional asymmetry alone, can increase the bending modulus. The likely cause is a stiffening of the compressed leaflet, which appears to be related to its gel transition but not identical with it. We finally show that the impact of cholesterol on differential stress depends on the relative strength of elastic and thermodynamic driving forces: if cholesterol solvates equally well in both leaflets, it will redistribute to cancel both leaflet tensions almost completely, but if its partitioning free energy prefers one leaflet over the other, the resulting distribution bias may even create differential stress. Because cells keep most of their lipid bilayers in an asymmetric nonequilibrium steady state, our findings suggest that biomembranes are elastically more complex than previously thought: besides a spontaneous curvature, they might also exhibit significant differential stress, which could strongly affect their curvature energetics.

Surface Shear Viscosity and Interleaflet Friction from Nonequilibrium Simulations of Lipid Bilayers

Nonequilibrium simulation protocols based on shear deformations are applied to determine the surface viscosity and interleaflet friction of lipid bilayers. At high shear rates, a non-Newtonian shear thinning regime is observed, but lower shear rates yield a Newtonian plateau and results that are consistent with equilibrium measurements based on fluctuation-dissipation theorems. Application to all-atom bilayers modeled with the CHARMM36 parameter set yields values for the surface viscosity that are consistent with microscopic measurements based on membrane protein diffusion but are approximately 10 times lower than more macroscopic experimental measurements. The interleaflet friction is about 10 times lower than experimental measurements. Trends across different lipids, temperatures, and ternary liquid-disordered phase mixtures produce results that are consistent with experimental diffusion constants. Application of the protocol to the liquid-ordered phase fails to yield a Newtonian plateau, suggesting more complex rheology.

Lipid Head Group Parameterization for GROMOS 54A8: A Consistent Approach with Protein Force Field Description

Membranes are a crucial component of both bacterial and mammalian cells, being involved in signaling, transport, and compartmentalization. This versatility requires a variety of lipid species to tailor the membrane’s behavior as needed, increasing the complexity of the system. Molecular dynamics simulations have been successfully applied to study model membranes and their interactions with proteins, elucidating some crucial mechanisms at the atomistic detail and thus complementing experimental techniques. An accurate description of the functional interplay of the diverse membrane components crucially depends on the selected parameters that define the adopted force field. A coherent parameterization for lipids and proteins is therefore needed. In this work, we propose and validate new lipid head group parameters for the GROMOS 54A8 force field, making use of recently published parametrizations for key chemical moieties present in lipids. We make use additionally of a new canonical set of partial charges for lipids, chosen to be consistent with the parameterization of soluble molecules such as proteins. We test the derived parameters on five phosphocholine model bilayers, composed of lipid patches four times larger than the ones used in previous studies, and run 500 ns long simulations of each system. Reproduction of experimental data like area per lipid and deuterium order parameters is good and comparable with previous parameterizations, as well as the description of liquid crystal to gel-phase transition. On the other hand, the orientational behavior of the head groups is more realistic for this new parameter set, and this can be crucial in the description of interactions with other polar molecules. For that reason, we tested the interaction of the antimicrobial peptide lactoferricin with two model membranes showing that the new parameters lead to a weaker peptide–membrane binding and give a more realistic outcome in comparing binding to antimicrobial versus mammal membranes.

Effect of Triclosan and Chloroxylenol on Bacterial Membranes

Triclosan and chloroxylenol are broad-spectrum biocides used extensively in healthcare and consumer products. They have been suggested to perturb the structure of bacterial membranes, but studies so far have not considered that most bacterial membranes contain large amounts of branched-chain lipids. Here, molecular dynamics simulation is used to examine the effect of the two biocides on membranes consisting of lipids with methyl-branched chains, cyclopropanated chains, and nonbranched chains. It is shown that triclosan and chloroxylenol induced a phase transition in membranes from a liquid-crystalline to a liquid-ordered phase irrespective of the presence and nature of branching groups. At high concentration, chloroxylenol promoted chain interdigitation. Our results suggest that triclosan and chloroxylenol decrease the degree of fluidity of membranes and that this effect is more pronounced in bacterial membranes. As a result, their biocidal activity could be associated with a change in the function of membrane proteins.

Computational Modeling of Realistic Cell Membranes

Cell membranes contain a large variety of lipid types and are crowded with proteins, endowing them with the plasticity needed to fulfill their key roles in cell functioning. The compositional complexity of cellular membranes gives rise to a heterogeneous lateral organization, which is still poorly understood. Computational models, in particular molecular dynamics simulations and related techniques, have provided important insight into the organizational principles of cell membranes over the past decades. Now, we are witnessing a transition from simulations of simpler membrane models to multicomponent systems, culminating in realistic models of an increasing variety of cell types and organelles. Here, we review the state of the art in the field of realistic membrane simulations and discuss the current limitations and challenges ahead.

Structural and barrier properties of the skin ceramide lipid bilayer: a molecular dynamics simulation study

Skin provides excellent protection against the harsh external environment and foreign substances. The lipid matrix of the stratum corneum, which contains various kinds of ceramides, plays a major role in the barrier function of the skin. Here we report a study of the effects of ceramide type on the structural and transport properties of ceramide bilayers using molecular dynamics (MD) simulations. Specifically, the effects of headgroup chemistry (number and positions of hydroxyl groups) and tail structure (unsaturation of the sphingoid moiety) on the structural and transport properties of various ceramide bilayers at 310 K were analyzed. Theoretical results for structural properties such as area per lipid, bilayer thickness, lateral arrangement, order parameter, and hydrogen bonding are reported here and compared with corresponding experimental data. Our study revealed that the presence of a double bond disrupts the bilayer packing, which leads to a low area compressibility modulus, a large area per lipid, and low bilayer thickness. Furthermore, the effect of structural changes on water permeation was studied using steered MD simulations. Water permeation was found to be influenced by headgroup polarity, chain packing, and the ability of the water to hydrogen bond with the ceramides. The molecular-level information obtained from the current study should aid the design of mixed bilayer systems with desired properties and provide the basis for the development of higher order coarse-grained models.

Molecular dynamics simulation study of the effect of halothane on mixed DPPC/DPPE phospholipid membranes

We report results of a molecular dynamics simulation study of the effect of one general anesthetic, halothane, on some properties of mixed DPPC/DPPE phospholipid membranes. This is a suitable model for the study of simple, two-phospholipid membrane systems. From the simulation runs, we determined several membrane properties for five different molecular proportions of DPPC/DPPE. The effect of halothane on the studied membrane properties (area per lipid molecule, density of membrane, order parameter, etc.) was rather small. The distribution of halothane is not uniform through the bilayer thickness. Instead, there is a maximum of anesthetic concentration around 1.2 nm from the center of the membrane. The anesthetic molecule is located close to the phospholipid headgroups. The position of the halothane density maximum depends slightly on the DPPC/DPPE molar proportion. Snapshots taken over the plane of the membrane, as well as calculated two-dimensional radial distribution functions show that the anesthetic has no preference for either phospholipid (DPPC or DPPE). Our results indicate that this anesthetic molecule has only small effects on DPPC/DPPE mixed membranes. In addition, halothane displays no preferential location around DPPC or DPPE. This is probably due to the hydrophobic nature of halothane and to the fact that the chosen phospholipids have the same hydrophobic tails.

Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations

Background: Molecular dynamics (MD) simulations have become an important tool to provide insight into molecular processes involving biomolecules such as proteins, DNA, carbohydrates and membranes. As these processes cover a wide range of time scales, multiple time-step integration methods are often employed to increase the speed of MD simulations. For example, in the twin-range (TR) scheme, the nonbonded forces within the long-range cutoff are split into a short-range contribution updated every time step (inner time step) and a less frequently updated mid-range contribution (outer time step). The presence of different time steps can, however, cause numerical artefacts. Methods: The effects of multiple time-step algorithms at interfaces between polar and apolar media are investigated with MD simulations. Such interfaces occur with biological membranes or proteins in solution. Results: In this work, it is shown that the TR splitting of the nonbonded forces leads to artificial density increases at interfaces. The presence of the observed artefacts was found to be independent of the interface shape and the thermostatting method used. It is further shown that integration with an impulse-wise reversible reference system propagation algorithm (RESPA) only shifts the occurrence of density artefacts towards larger outer time steps. Using a single-range (SR) treatment of the nonbonded interactions, on the other hand, resolves the density issue for pairlist-update periods of up to 40 fs. Conclusion: A SR scheme avoids numerical artefacts and offers an interesting alternative to TR RESPA with respect to performance optimization.

Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations

Background: Molecular dynamics (MD) simulations have become an important tool to provide insight into molecular processes involving biomolecules such as proteins, DNA, carbohydrates and membranes. As these processes cover a wide range of time scales, multiple time-step integration methods are often employed to increase the speed of MD simulations. For example, in the twin-range (TR) scheme, the nonbonded forces within the long-range cutoff are split into a short-range contribution updated every time step (inner time step) and a less frequently updated mid-range contribution (outer time step). The presence of different time steps can, however, cause numerical artefacts. Methods: The effects of multiple time-step algorithms at interfaces between polar and apolar media are investigated with MD simulations. Such interfaces occur with biological membranes or proteins in solution. Results: In this work, it is shown that the TR splitting of the nonbonded forces leads to artificial density increases at interfaces for weak coupling and Nosé-Hoover (chain) thermostats. It is further shown that integration with an impulse-wise reversible reference system propagation algorithm (RESPA) only shifts the occurrence of density artefacts towards larger outer time steps. Using a single-range (SR) treatment of the nonbonded interactions or a stochastic dynamics thermostat, on the other hand, resolves the density issue for pairlist-update periods of up to 40 fs. Conclusion: TR schemes are not advisable to use in combination with weak coupling or Nosé-Hoover (chain) thermostats due to the occurrence of significant numerical artifacts at interfaces.

Density artefacts at interfaces caused by multiple time-step effects in molecular dynamics simulations

Background: Molecular dynamics (MD) simulations have become an important tool to provide insight into molecular processes involving biomolecules such as proteins, DNA, carbohydrates and membranes. As these processes cover a wide range of time scales, multiple time-step integration methods are often employed to increase the speed of MD simulations. For example, in the twin-range (TR) scheme, the nonbonded forces within the long-range cutoff are split into a short-range contribution updated every time step (inner time step) and a less frequently updated mid-range contribution (outer time step). The presence of different time steps can, however, cause numerical artefacts. Methods: The effects of multiple time-step algorithms at interfaces between polar and apolar media are investigated with MD simulations. Such interfaces occur with biological membranes or proteins in solution. Results: In this work, it is shown that the TR splitting of the nonbonded forces leads to artificial density increases at interfaces for weak coupling and Nosé-Hoover (chain) thermostats. It is further shown that integration with an impulse-wise reversible reference system propagation algorithm (RESPA) only shifts the occurrence of density artefacts towards larger outer time steps. Using a single-range (SR) treatment of the nonbonded interactions or a stochastic dynamics thermostat, on the other hand, resolves the density issue for pairlist-update periods of up to 40 fs. Conclusion: TR schemes are not advisable to use in combination with weak coupling or Nosé-Hoover (chain) thermostats due to the occurrence of significant numerical artifacts at interfaces.

Rational Design of Thermostable Carbonic Anhydrase Mutants Using Molecular Dynamics Simulations

The stability of enzymes is critical for their application in industrial processes, which generally require different conditions from the natural enzyme environment. Both rational and random protein engineering approaches have been used to increase stability, with the latter requiring extensive experimental effort for the screening of variants. Moreover, some general rules addressing the molecular origin of protein thermostability have been established. Herein, we demonstrate the use of molecular dynamics simulations to gain molecular level understanding of protein thermostability and to engineer stabilizing mutations. Carbonic anhydrase (CA) is an enzyme with a high potential for biotechnological carbon capture applications, provided it can be engineered to withstand the high temperature process environments, inevitable in most gas treatment units. In this study, we used molecular dynamics simulations at 343, 353, and 363 K to study the relationship between structure flexibility and thermostability in bacterial α-CAs and applied this knowledge to the design of mutants with increased stability. The most thermostable α-CA known, TaCA from Thermovibrio ammonificans, had the most rigid structure during molecular dynamics simulations, but also showed regions with high flexibility. The most flexible amino acids in these regions were identified from root mean square fluctuation (RMSF) studies, and stabilizing point mutations were predicted based on their capacity to improve the calculated free energy of unfolding. Disulfide bonds were also designed at sites with suitable geometries and selected based on their location at flexible sites, assessed by B-factor calculation. Molecular dynamics simulations allowed the identification of five mutants with lower RMSF of the overall structure at 400 K, compared to wild-type TaCA. Comparison of free-energy landscapes between wild-type TaCA and the most promising mutants, Pro165Cys-Gln170Cys and Asn140Gly, showed an increased conformational stability of the mutants at 400 K.

Verifiability in computer-aided research: the role of digital scientific notations at the human-computer interface

Most of today's scientific research relies on computers and software for processing scientific information. Examples of such computer-aided research are the analysis of experimental data or the simulation of phenomena based on theoretical models. With the rapid increase of computational power, scientific software has integrated more and more complex scientific knowledge in a black-box fashion. As a consequence, its users do not know, and do not even have a chance of finding out, which assumptions and approximations their computations are based on. This black-box nature of scientific software has made the verification of much computer-aided research close to impossible. The present work starts with an analysis of this situation from the point of view of human-computer interaction in scientific research. It identifies the key role of digital scientific notations at the human-computer interface, reviews the most popular ones in use today, and describes a proof-of-concept implementation of Leibniz, a language designed as a verifiable digital scientific notation for models formulated as mathematical equations.

Balancing Force Field Protein-Lipid Interactions To Capture Transmembrane Helix-Helix Association

Atomistic simulations have recently been shown to be sufficiently accurate to reversibly fold globular proteins and have provided insights into folding mechanisms. Gaining similar understanding from simulations of membrane protein folding and association would be of great medical interest. All-atom simulations of the folding and assembly of transmembrane protein domains are much more challenging, not least due to very slow diffusion within the lipid bilayer membrane. Here, we focus on a simple and well-characterized prototype of membrane protein folding and assembly, namely the dimerization of glycophorin A, a homodimer of single transmembrane helices. We have determined the free energy landscape for association of the dimer using the CHARMM36 force field. We find that the native structure is a metastable state, but not stable as expected from experimental estimates of the dissociation constant and numerous experimental structures obtained under a variety of conditions. We explore two straightforward approaches to address this problem and demonstrate that they result in stable dimers with dissociation constants consistent with experimental data.

Consensus Conformations of Dinucleoside Monophosphates Described with Well-Converged Molecular Dynamics Simulations

Dinucleoside monophosphates (DNMPs) have been described using various experimental approaches as flexible molecules which generate ensembles populating at least a small set of different conformations in solution. However, due to limitations of each approach in its ability to delineate the ensemble of conformations, an accurate and quantitative description of certain conformational features has not been performed for all DNMPs. Here, we apply a temperature replica-exchange molecular dynamics approach to fully and quickly converge conformational distributions of all RNA DNMPs immersed in the TIP3P water model using the AMBER ff14 force field. For a selection of DNMPs, the conformational ensembles were also generated when immersed in the OPC water model using alternative AMBER and CHARMM force fields. The OPC water model and other force field choices did not introduce new conformational classes but shifted the populations among existing conformations. Except for pyrimidine-pyrimidine dinucleosides, all other DNMPs populated four major conformations (which are defined in the main text and labeled A-form, Ladder, Inverted, and Sheared), in addition to an Extended form. Pyrimidine-pyrimidines did not generate the Sheared conformation. Distinguishing features and stabilizing factors of each conformation were identified and assessed based on the known experimental interpretations. The configuration of the glycosidic bond and the nonbonding interactions of hydrogen bond acceptors with the 2'-hydroxyl group were found to play determining roles in stabilizing particular conformations which could serve as a guide for potential force field modifications to improve the accuracy. Additionally, we computed stacking free energies based on the DNMP conformational distributions and found significant discrepancies with a previous study. Our investigation determined that the AMBER force field was incorrectly implemented in the previous study. In the future, this simulation approach can be used to quickly analyze the effects of new force field modifications in shifting the conformational populations of DNMPs, and can can be further applied to foresee such effects in larger RNA motifs including tetranucleotides and tetraloops.

Validation and Comparison of Force Fields for Native Cyclodextrins in Aqueous Solution

Molecular dynamics simulations of native α-, β-, and γ-cyclodextrin in aqueous solution have been conducted with the goal to investigate the performance of the CHARMM36 force field, the AMBER-compatible q4md-CD force field, and five variants of the GROMOS force field. The properties analyzed are structural parameters derived from X-ray diffraction and NMR experiments as well as hydrogen bonds and hydration patterns, including hydration free enthalpies. Recent revisions of the torsional-angle parameters for carbohydrate systems within the GROMOS family of force fields lead to a significant improvement of the agreement between simulated and experimental NMR data. Therefore, we recommend using the variant 53A6_GLYC instead of 53A6 and 56A6_{CARBO_R} or 2016H66 instead of 56A6_CARBO to simulate cyclodextrins in solution. The CHARMM36 and q4md-CD force fields show a similar performance as the three recommended GROMOS parameter sets. A significant difference is the more flexible nature of the cyclodextrins modeled with the CHARMM36 and q4md-CD force fields compared to the three recommended GROMOS parameter sets.

The effect of hydronium ions on the structure of phospholipid membranes

This work studies the mechanisms by which hydronium ions modulate the structure of phospholipid bilayers.