The role of simulation in biomathematical modeling

This paper is a general introduction to the field of biomathematical modeling. Biomathematical modeling is divided into three parts: the derivation of models, the fitting of models to data, and the simulation of data from models. This paper focuses on the simulation of data from models. The uses of simulation, the potential users of simulation, and simulation software are described.     

Molecular dynamics simulation of hydrated Nafion with a reactive force field for water

We apply a newly parameterized central force field to highlight the problem of proton transport in fuel cell membranes and show that central force fields are potential candidates to describe chemical reactions on a classical level. After a short sketch of the parameterization of the force field, we validate the obtained force field for several properties of water. The experimental and simulated radial distribution functions are reproduced very accurately as a consequence of the applied parameterization procedure. Further properties, geometry, coordination, diffusion coefficient and density, are simulated adequately for our purposes. Afterwards we use the new force field for the molecular dynamics simulation of a swollen polyelectrolyte membrane similar to the widespread Nafion 117. We investigate the equilibrated structures, proton transfer, lifetimes of hydronium ions, the diffusion coefficients, and the conductivity in dependence of water content. In a short movie we demonstrate the ability of the obtained force field to describe the bond breaking/formation, and conclude that this force field can be considered as a kind of a reactive force field. The investigations of the lifetimes of hydronium ions give us the information about the kinetics of the proton transfer in a membrane with low water content. We found the evidence for the second order reaction. Finally, we demonstrate that the model is simple enough to handle the large systems sufficient to calculate the conductivity from molecular dynamics simulations. The detailed analysis of the conductivity reveals the importance of the collective moving of hydronium ions in membrane, which might give an interesting encouragement for further development of membranes. Figure: The structure of water in one pore of the highly hydrated Nafion membranes. Figure The structure of water in one of pore of the highly hydrated Nafion membrane Electronic supplementary material The online version of this article (doi:) contains supplementary material, which is available to authorized users.     

Why the OPLS-AA force field cannot produce the β-hairpin structure of H1 peptide in solution when comparing with the GROMOS 43A1 force field?

The optimal combination of force field and water model is an essential problem that is able to increase molecular dynamics simulation quality for different types of proteins and peptides. In this work, an attempt has been made to explore the problem by studying H1 peptide using four different models based on different force fields, water models and electrostatic schemes. The driving force for H1 peptide conformation transition and the reason why the OPLS-AA force field cannot produce the β-hairpin structure of H1 peptide in solution while the GROMOS 43A1 force field can do were investigated by temperature replica exchange molecular dynamics simulation (T-REMD). The simulation using the GROMOS 43A1 force field preferred to adopt a β-hairpin structure, which was in good agreement with the several other simulations and the experimental evidences. However, the simulation using the OPLS-AA force field has a significant difference from the simulations with the GROMOS 43A1 force field simulation. The results show that the driving force in H1 peptide conformation transition is solvent exposure of its hydrophobic residues. However, the subtle balances between residue-residue interactions and residue-solvent interaction are disrupted by using the OPLS-AA force field, which induced the reduction in the number of residue-residue contact. Similar solvent exposure of the hydrophobic residues is observed for all the conformations sampled using the OPLS-AA force field. For H1 peptide which exhibits large solvent exposure of the hydrophobic residues, the GROMOS 43A1 force field with the SPC water model can provide more accurate results.     

Making optimal use of empirical energy functions: force-field parameterization in crystal space

Today's energy functions are not able yet to distinguish reliably between correct and almost correct protein models. Improving these near-native models is currently a major bottle-neck in homology modeling or experimental structure determination at low resolution. Increasingly accurate energy functions are required to complete the "last mile of the protein folding problem," for example during a molecular dynamics simulation. We present a new approach to reach this goal. For 50 high resolution X-ray structures, the complete unit cell was reconstructed, including disordered water molecules, counter ions, and hydrogen atoms. Simulations were then run at the pH at which the crystal was solved, while force-field parameters were iteratively adjusted so that the damage done to the structures was minimal. Starting with initial parameters from the AMBER force field, the optimization procedure converged at a new force field called YAMBER (Yet Another Model Building and Energy Refinement force field), which is shown to do significantly less damage to X-ray structures, often move homology models in the right direction, and occasionally make them look like experimental structures. Application of YAMBER during the CASP5 structure prediction experiment yielded a model for target 176 that was ranked first among 150 submissions. Due to its compatibility with the well-established AMBER format, YAMBER can be used by almost any molecular dynamics program. The parameters are freely available from www.yasara.org/yamber.     

Water activity coefficients in aqueous amino acid solutions by molecular dynamics simulation: 1. Force field development

New force fields for molecular dynamics (MD) simulation of aqueous zwitterionic amino acid simulations were developed. These were especially designed to calculate activity coefficient of water in amino acid solutions with high accuracy. For example, aqueous solutions of the following amino acids were considered: glycine, alanine, α-aminobutyric acid, α-aminovalerianic acid, valine and leucine. The force fields were obtained by quantum chemical calculations using B3LYP/6-31G and MP2/6-311(d,p) model theories in combination with the Merz–Kollmann–Singh scheme. To further increase the accuracy of the force field, a polarised continuum was considered in all quantum chemical calculations. Water activity coefficients obtained from MD using different all-purpose literature force fields, namely, OPLS, AMBER ff03 and GROMOS 53A6 as well as experimental data are compared with the results utilising the new force field. The new force field is shown to give better results compared with experimental data than existing force fields.     

A Robust Water Potential Parameterisation

Abstract

We compare molecular dynamics simulation results for the properties of liquid water predicted by four novel water potential models. These models are designed as a combination of parameters taken from the dedicated but brittle TIP3P water potential, and the more flexible but less accurate parameterisations such as the Dreiding and Universal force fields. We find that a hybrid of Dreiding and TIP3P delivers the best results, yielding a density, diffusion coefficient and radial distribution function in good agreement with experiment, performing in some respects even better than the dedicated reference TIP3P model. Another Dreiding based force field predicts semi-quantitative results for the water structure and dynamics while the Universal force field based models are incapable of simulating a condensed phase of water at all, continuing to expand indefinitely. These observations are useful for selecting and designing robust water force field parameterisations that can be used for general simulation purposes.     

Plant growth and architectural modelling and its applications. Preface

Over the last decade, a growing number of scientists around the world have invested in research on plant growth and architectural modelling and applications (often abbreviated to plant modelling and applications, PMA). By combining physical and biological processes, spatially explicit models have shown their ability to help in understanding plant–environment interactions. This Special Issue on plant growth modelling presents new information within this topic, which are summarized in this preface. Research results for a variety of plant species growing in the field, in greenhouses and in natural environments are presented. Various models and simulation platforms are developed in this field of research, opening new features to a wider community of researchers and end users. New modelling technologies relating to the structure and function of plant shoots and root systems are explored from the cellular to the whole-plant and plant-community levels.     

Development of an AMBER-compatible transferable force field for poly(ethylene glycol) ethers (glymes)

An all-atom force field consistent with the general AMBER force field (GAFF) format for poly(ethylene glycol) dimethyl ether (diglyme or G2) was developed by fitting to experimental liquid densities and dielectric constants. Not surprisingly, the new force field gives excellent agreement with experimental liquid phase densities and dielectric constants over a wide temperature range. Other dynamic and thermodynamic properties of liquid G2 such as its self-diffusion coefficient, shear viscosity, and vaporization enthalpy were also calculated and compared to experimental data. For all of the properties studied, the performance of the proposed new force field is better than that of the standard GAFF force field. The force field parameters were transferred to model two other poly(ethylene glycol) ethers: monoglyme (G1) and tetraglyme (G4). The predictive ability of the modified force field for G1 and G4 was significantly better than that of the original GAFF force field. The proposed force field provides an alternative option for the simulation of mixtures containing glymes using GAFF-compatible force fields, particularly for electrochemical applications. The accuracy of a previously published force field based on the OPLS-AA format and the accuracies of two modified versions of that force field were also examined for G1, G2, and G4. It was found that the original OPLS-AA force field is superior to the modified versions of it, and that it has a similar accuracy to the proposed new GAFF-compatible force field.

Open image in new window

Graphical abstract Transferability of an AMBER-compatible force field parameterized for G2 to other glymes

     

Analysis of dissociation process for gas hydrates by molecular dynamics simulation

The dissociation processes of methane and carbon dioxide hydrates were investigated by molecular dynamics simulation. The simulations were performed with 368 water molecules and 64 gas molecules using NPT ensembles. The TraPPE (single-site) and 5-site models were adopted for methane molecules. The EPM2 (3-site) and SPC/E models were used for carbon dioxide and water molecules, respectively. The simulations were carried out at 270 K and 5.0 MPa for hydrate stabilisation. Then, temperature was increased up to 370 K. The temperature increasing rates were 0.1–20 TK/s. The gas hydrates dissociated during increasing temperature or at 370 K. The potential models of methane molecule did not much influence the dissociation process of methane hydrate. The mechanisms of dissociation process were analysed with the coordination numbers and mean square displacements. It was found that the water cages break down first, then the gas molecules escape from the water cages. The methane hydrate was more stable than the carbon dioxide hydrate at the calculated conditions.     

Recent progress in general force fields of small molecules

Recent advances in computational hardware and free energy algorithms enable a broader application of molecular simulation of binding interactions between receptors and small-molecule ligands. The underlying molecular mechanics force fields (FFs) for small molecules have also achieved advancements in accuracy, user-friendliness, and speed during the past several years (2018–2020). Besides the expansion of chemical space coverage of ligand-like molecules among major popular classical additive FFs and polarizable FFs, new charge models have been proposed for better accuracy and transferability, new chemical perception of avoiding predefined atom types have been applied, and new automated parameterization toolkits, including machine learning approaches, have been developed for users’ convenience.     

Polarizable Water Model for the Coarse-Grained MARTINI Force Field

Coarse-grained (CG) simulations have become an essential tool to study a large variety of biomolecular processes, exploring temporal and spatial scales inaccessible to traditional models of atomistic resolution. One of the major simplifications of CG models is the representation of the solvent, which is either implicit or modeled explicitly as a van der Waals particle. The effect of polarization, and thus a proper screening of interactions depending on the local environment, is absent. Given the important role of water as a ubiquitous solvent in biological systems, its treatment is crucial to the properties derived from simulation studies. Here, we parameterize a polarizable coarse-grained water model to be used in combination with the CG MARTINI force field. Using a three-bead model to represent four water molecules, we show that the orientational polarizability of real water can be effectively accounted for. This has the consequence that the dielectric screening of bulk water is reproduced. At the same time, we parameterized our new water model such that bulk water density and oil/water partitioning data remain at the same level of accuracy as for the standard MARTINI force field. We apply the new model to two cases for which current CG force fields are inadequate. First, we address the transport of ions across a lipid membrane. The computed potential of mean force shows that the ions now naturally feel the change in dielectric medium when moving from the high dielectric aqueous phase toward the low dielectric membrane interior. In the second application we consider the electroporation process of both an oil slab and a lipid bilayer. The electrostatic field drives the formation of water filled pores in both cases, following a similar mechanism as seen with atomistically detailed models.     

Mesoscopic lateral diffusion in lipid bilayers

The lateral diffusion in bilayers is modeled with a multiscale mesoscopic simulation. The methodology consists of two simulations, where the first employs atomistic models to obtain exact results for the mesoscopic model. The second simulation takes the results obtained from the first to parameterize an effective force field that is employed in a new coarse-grained model. The multiscale aspect of this scheme occurs at the point where the microscopic time-averaged results of the first simulation are employed to parameterize the second simulation that operates in a higher spatial and temporal domain. The results of both simulation schemes give quantitative information on the details associated with lipid lateral diffusion.     

All-atom and coarse-grained force fields for polydimethylsiloxane

By combining the bottom-up and top-down approaches, we have developed a new all-atom (AA) force field from quantum mechanics and experimental data and a new coarse grained (CG) force field from AA simulation and experimental data, for polydimethylsiloxane (PDMS). The AA force field is developed based on the TEAM force field database. The CG force field uses a mapping rule that splits the connecting oxygen into neighbouring CG beads to maintain the charge neutrality of the beads, analytical functional forms including anharmonic terms in the valence terms, and the temperature-dependent free-energy functional form to describe the inter-bead interactions. Broad range of thermodynamic properties of PDMS including density, surface tension, solubility parameter, radius of gyration and glass transition temperature are calculated to validate the force fields, and good agreements with the experimental data are obtained.     

Space as the final frontier in stochastic simulations of biological systems

Recent technological and theoretical advances are only now allowing the simulation of detailed kinetic models of biological systems that reflect the stochastic movement and reactivity of individual molecules within cellular compartments. The behavior of many systems could not be properly understood without this level of resolution, opening up new perspectives of using computer simulations to accelerate biological research. We review the modeling methodology applied to stochastic spatial models, also to the attention of non-expert potential users. Modeling choices, current limitations and perspectives of improvement of current general-purpose modeling/simulation platforms for biological systems are discussed.     

Comparison of a QM/MM force field and molecular mechanics force fields in simulations of alanine and glycine "dipeptides" (Ace-Ala-Nme and Ace-Gly-Nme) in water in relation to the problem of modeling the unfolded peptide backbone in solution

We compare the conformational distributions of Ace-Ala-Nme and Ace-Gly-Nme sampled in long simulations with several molecular mechanics (MM) force fields and with a fast combined quantum mechanics/molecular mechanics (QM/MM) force field, in which the solute's intramolecular energy and forces are calculated with the self-consistent charge density functional tight binding method (SCCDFTB), and the solvent is represented by either one of the well-known SPC and TIP3P models. All MM force fields give two main states for Ace-Ala-Nme, beta and alpha separated by free energy barriers, but the ratio in which these are sampled varies by a factor of 30, from a high in favor of beta of 6 to a low of 1/5. The frequency of transitions between states is particularly low with the amber and charmm force fields, for which the distributions are noticeably narrower, and the energy barriers between states higher. The lower of the two barriers lies between alpha and beta at values of psi near 0 for all MM simulations except for charmm22. The results of the QM/MM simulations vary less with the choice of MM force field; the ratio beta/alpha varies between 1.5 and 2.2, the easy pass lies at psi near 0, and transitions between states are more frequent than for amber and charmm, but less frequent than for cedar. For Ace-Gly-Nme, all force fields locate a diffuse stable region around phi = pi and psi = pi, whereas the amber force field gives two additional densely sampled states near phi = +/-100 degrees and psi = 0, which are also found with the QM/MM force field. For both solutes, the distribution from the QM/MM simulation shows greater similarity with the distribution in high-resolution protein structures than is the case for any of the MM simulations.     

The bridge between first-principles calculations and grand canonical Monte Carlo simulations: Morse and Lennard-Jones force fields

We use the modelled adsorbate (methane) and adsorbent (carbon nanotube) to explore the choice of site-to-site force fields in order to bridge the gap between first-principles calculation and grand canonical Monte Carlo. By using the mature classical Lennard-Jones potential-based results as benchmark, it is found if the first-principles calculation is performed at the higher MP2/6-311G** level, Morse potential model should be adopted to achieve the force field parameters. It is expected that this work provides useful information for the building of the multiscale simulation method and the selection of accurate force field models.     

A gradual update method for simulating the steady-state solution of stiff differential equations in metabolic circuits

Numerical simulation of differential equation systems plays a major role in the understanding of how metabolic network models generate particular cellular functions. On the other hand, the classical and technical problems for stiff differential equations still remain to be solved, while many elegant algorithms have been presented. To relax the stiffness problem, we propose new practical methods: the gradual update of differential-algebraic equations based on gradual application of the steady-state approximation to stiff differential equations, and the gradual update of the initial values in differential-algebraic equations. These empirical methods show a high efficiency for simulating the steady-state solutions for the stiff differential equations that existing solvers alone cannot solve. They are effective in extending the applicability of dynamic simulation to biochemical network models. Electronic supplementary material  The online version of this article (doi:) contains supplementary material, which is available to authorized users.