Calculating chemically accurate redox potentials for engineered flavoproteins from classical molecular dynamics free energy simulations

The tricyclic isoalloxazine nucleus of the redox cofactors flavin mononucleotide (FMN) and flavin adenine dinucleotide (FAD) acts as an electron sink in life-sustaining biological electron transfer (eT). The functional diversity of flavin-containing proteins (flavoproteins) transcends that of free flavins. A large body of experimental evidence attributes natural control of flavoprotein-mediated eT to tuning of the thermodynamic driving force by the protein environment. Understanding and engineering such modulation by the protein environment of the flavin redox potential (DeltaE(o)) is valuable in biotechnology and device design. In this study we employed classical molecular dynamics free energy simulations (MDFES), within a thermodynamic integration (TI) formalism, to calculate the change in FMN first reduction potential (DeltaDeltaE(o)(ox/sq)) imparted by 6 flavoprotein active site mutations. The combined performance of the AMBER ff03 (protein) and GAFF (cofactor) force fields was benchmarked against experimental data for mutations close to the isoalloxazine re- and si-faces that perturb the wild-type DeltaE(o)(ox/sq) value in Anabaena flavodoxin. The classical alchemical approach used in this study overestimates the magnitude of DeltaE(o) values, in common with other studies. Nevertheless, chemically accurate DeltaDeltaE(o) values--calculated to within 1 kcal mol(-1) of the experimental value--were obtained for five of the six mutations studied. We have shown that this approach is practical for quantitative in silico screening of the effect of mutations on the first reduction potential where experimental values and structural data are available for the wild-type flavoprotein. This approach promises to be useful as an integral part of future interdisciplinary strategies to engineer desired thermodynamic properties in flavoproteins of biotechnological interest.     

Unbinding of nicotine from the acetylcholine binding protein: steered molecular dynamics simulations

The ligand binding/unbinding process is critical to our understanding of the pharmacology of both the nicotinic acetylcholine receptor (nAChR) and the acetylcholine binding protein (AChBP). Steered molecular dynamics simulations were performed to learn about the unbinding process of the full agonist nicotine. Three different pulling models were designed to investigate the possible binding/unbinding pathways: radial and tangent models, and also a mixed model. Of the three, the tangent pulling model finally failed to dissociate nicotine from the ligand binding pocket. The efficiency of the pulling force profiles was superior, and the opening of the C-loop was smaller in the mixed pulling model than that in the radial model. The most favorable pathway for the cholinergic agonist nicotine to enter or leave the binding pocket is through the principal binding side, following a curvilinear track. Noticeably, it has been seen that the unbinding of the nicotine is concomitant with a global rotation of the protein-ligand complex which could be caused by the interactions of the ligand with protein at the tangent direction.     

Ranking protein–protein docking results using steered molecular dynamics and potential of mean force calculations

Crystallization of protein–protein complexes can often be problematic and therefore computational structural models are often relied on. Such models are often generated using protein–protein docking algorithms, where one of the main challenges is selecting which of several thousand potential predictions represents the most near‐native complex. We have developed a novel technique that involves the use of steered molecular dynamics (sMD) and umbrella sampling to identify near‐native complexes among protein–protein docking predictions. Using this technique, we have found a strong correlation between our predictions and the interface RMSD (iRMSD) in ten diverse test systems. On two of the systems, we investigated if the prediction results could be further improved using potential of mean force calculations. We demonstrated that a near‐native (<2.0 Å iRMSD) structure could be identified in the top‐1 ranked position for both systems. © 2016 Wiley Periodicals, Inc.     

Mechanically untying a protein slipknot: multiple pathways revealed by force spectroscopy and steered molecular dynamics simulations

Protein structure is highly diverse when considering a wide range of protein types, helping to give rise to the multitude of functions that proteins perform. In particular, certain proteins are known to adopt a knotted or slipknotted fold. How such proteins undergo mechanical unfolding was investigated utilizing a combination of single molecule atomic force microscopy (AFM), protein engineering, and steered molecular dynamics (SMD) simulations to show the mechanical unfolding mechanism of the slipknotted protein AFV3-109. Our results reveal that the mechanical unfolding of AFV3-109 can proceed via multiple parallel unfolding pathways that all cause the protein slipknot to untie and the polypeptide chain to completely extend. These distinct unfolding pathways proceed via either a two- or three-state unfolding process involving the formation of a well-defined, stable intermediate state. SMD simulations predict the same contour length increments for different unfolding pathways as single molecule AFM results, thus providing a plausible molecular mechanism for the mechanical unfolding of AFV3-109. These SMD simulations also reveal that two-state unfolding is initiated from both the N- and C-termini, while three-state unfolding is initiated only from the C-terminus. In both pathways, the protein slipknot was untied during unfolding, and no tightened slipknot conformation was observed. Detailed analysis revealed that interactions between key structural elements lock the knotting loop in place, preventing it from shrinking and the formation of a tightened slipknot conformation. Our results demonstrate the bifurcation of the mechanical unfolding pathway of AFV3-109 and point to the generality of a kinetic partitioning mechanism for protein folding/unfolding.     

Calculating potentials of mean force and diffusion coefficients from nonequilibrium processes without Jarzynski's equality

In general, the direct application of the Jarzynski equality (JE) to reconstruct potentials of mean force (PMFs) from a small number of nonequilibrium unidirectional steered molecular-dynamics (SMD) paths is hindered by the lack of sampling of extremely rare paths with negative dissipative work. Such trajectories that transiently violate the second law of thermodynamics are crucial for the validity of JE. As a solution to this daunting problem, we propose a simple and efficient method, referred to as the FR method, for calculating simultaneously both the PMF U(z) and the corresponding diffusion coefficient D(z) along a reaction coordinate z for a classical many-particle system by employing a small number of fast SMD pullings in both forward (F) and time reverse (R) directions, without invoking JE. By employing Crooks [Phys. Rev. E 61, 2361 (2000)] transient fluctuation theorem (that is more general than JE) and the stiff-spring approximation, we show that (i) the mean dissipative work W(d) in the F and R pullings is the same, (ii) both U(z) and W(d) can be expressed in terms of the easily calculable mean work of the F and R processes, and (iii) D(z) can be expressed in terms of the slope of W(d). To test its viability, the FR method is applied to determine U(z) and D(z) of single-file water molecules in single-walled carbon nanotubes (SWNTs). The obtained U(z) is found to be in very good agreement with the results from other PMF calculation methods, e.g., umbrella sampling. Finally, U(z) and D(z) are used as input in a stochastic model, based on the Fokker-Planck equation, for describing water transport through SWNTs on a mesoscopic time scale that in general is inaccessible to MD simulations.     

Atomistic molecular dynamics simulations of chemical force microscopy

Chemical force microscopy and related force measurement techniques have emerged as powerful tools for studying fundamental interactions central to understanding adhesion and tribology at the molecular scale. However, detailed interpretation of these interactions requires knowledge of chemical and physical processes occurring in the region of the tip-sample junction that experiments cannot provide, such as atomic-scale motions and distribution of forces. In an effort to address some of these open issues, atomistic molecular dynamics simulations were performed modeling a chemical force microscope stylus covered with a planar C12 alkylthiolate self-assembled monolayer (SAM) interacting with a solid wall. A complete loading-unloading sequence was simulated under conditions of near-constant equilibrium, approximating the case of infinitely slow tip motion. In the absence of the solid wall, the stylus film existed in a fluid state with structural and dynamic properties similar to those of the analogous planar SAM at an elevated temperature. When the wall was brought into contact with the stylus and pressed against it, a series of reversible changes occurred culminating with solidification of the SAM film at the largest compressive force. During loading, the chemical composition of the contact changed, as much of the film's interior was exposed to the wall. At all tip heights, the distribution of forces within the contact zone was uneven and subject to large local fluctuations. Analysis using the Johnson-Kendall-Roberts, Derjaguin-Muller-Toporov, and Hertz contacts mechanics models revealed significant deviations from the simulation results, with the JKR model providing best overall agreement. Some of the discrepancies found would be overlooked in an actual experiment, where, unlike the simulations, contact area is not separately known, possibly producing a misleading or incorrect interpretation of experimental results. These shortcomings may be improved upon by using a model that correctly accounts for the finite thickness of the compliant components and nonlinear elastic effects.     

Introduction of steered molecular dynamics into UNRES coarse‐grained simulations package

In this article, an implementation of steered molecular dynamics (SMD) in coarse‐grain UNited RESidue (UNRES) simulations package is presented. Two variants of SMD have been implemented: with a constant force and a constant velocity. The huge advantage of SMD implementation in the UNRES force field is that it allows to pull with the speed significantly lower than the accessible pulling speed in simulations with all‐atom representation of a system, with respect to a reasonable computational time. Therefore, obtaining pulling speed closer to those which appear in the atomic force spectroscopy is possible. The newly implemented method has been tested for behavior in a microcanonical run to verify the influence of introduction of artificial constrains on keeping total energy of the system. Moreover, as time dependent artificial force was introduced, the thermostat behavior was tested. The new method was also tested via unfolding of the Fn3 domain of human contactin 1 protein and the I27 titin domain. Obtained results were compared with Gø‐like force field, all‐atom force field, and experimental results. © 2017 Wiley Periodicals, Inc.     

CHARMM force field and molecular dynamics simulations of protonated polyethylenimine

As a gene delivery vector, polyethylenimine (PEI) shows one of the highest transfection efficiencies, while effectively protecting DNA from enzyme degradation. The distinctive charge pattern of protonated PEI is widely considered responsible for fundamental process such as DNA condensation into PEI/DNA polyplexes (which are able to enter cells via endocytosis), proton sponge effect (which triggers the release of polyplexes from endosome), and release of DNA from polyplexes (to be further processed inside the nucleus). Our investigations are largely motivated by the crucial need for a realistic molecular mechanics force field (FF) for PEI, and, accordingly, we focus on two major issues: (1) development of a new atomistic (CHARMM) FF for PEI in different protonation states, rigorously derived from high‐quality ab initio calculations performed on model polymers, and (2) molecular dynamics investigations of solvated PEI, providing a detailed picture of the dynamic structuring thereof in dependence on their size and protonation state. The modeled PEI chains are essentially described in terms of gyration radius, end‐to‐end distance, persistence length, radial distribution functions, coordination numbers, and diffusion coefficients. They turn out to be more rigid than in other computational studies and we find diffusion coefficients in fair agreement with experimental data. The developed atomistic FF proves adequate for the realistic modeling of the size and protonation behavior of linear PEI, either as individual chains or composing polyplexes. © 2017 Wiley Periodicals, Inc.     

Vesicle shapes from molecular dynamics simulations

Lipid bilayer membranes are known to form various structures such as large sheets or vesicles. When the two leaflets of the bilayer have an equal composition, the membrane preferentially forms a flat sheet or a spherical vesicle. However, a difference in the composition of the two leaflets may result in a curved bilayer or in a wide variety of vesicle shapes. Vesicles with different shapes have already been shown in experiments and diverse vesicle shapes have been predicted theoretically from energy minimization of continuous curves. Here we present a molecular dynamics study of the effect of small changes in the phospholipid headgroups on the spontaneous curvature of the bilayer and on the resulting vesicle shape transformations. Small asymmetries in the bilayers already result in high spontaneous curvature and large vesicle deformations. Vesicle shapes that are formed include ellipsoids, discoids, pear-shaped vesicles, cup-shaped vesicles, as well as budded vesicles. Comparison of these vesicles with theoretically derived vesicle shapes shows both resemblances and differences.     

Modern protein force fields behave comparably in molecular dynamics simulations

Several molecular dynamics simulations were performed on three proteins--bovine apo-calbindin D9K, human interleukin-4 R88Q mutant, and domain IIA of bacillus subtilis glucose permease--with each of the AMBER94, CHARMM22, and OPLS-AA force fields as implemented in CHARMM. Structural and dynamic properties such as solvent-accessible surface area, radius of gyration, deviation from their respective experimental structures, secondary structure, and backbone order parameters are obtained from each of the 2-ns simulations for the purpose of comparing the protein portions of these force fields. For one of the proteins, the interleukin-4 mutant, two independent simulations were performed using the CHARMM22 force field to gauge the sensitivity of some of these properties to the specific trajectory. In general, the force fields tested performed remarkably similarly with differences on the order of those found for the two independent trajectories of interleukin-4 with CHARMM22. When all three proteins are considered together, no force field showed any consistent trend in variations for most of the properties monitored in the study.     

Liquid properties of dimethyl ether from molecular dynamics simulations using Ab Initio force fields

We have used molecular dynamic simulations to study the structural and dynamical properties of liquid dimethyl ether (DME) with a newly constructed ab initio force field in this article. The ab initio potential energy data were calculated at the second order Møller‐Plesset (MP2) perturbation theory with Dunning's correlation consistent basis sets (up to aug‐cc‐pVQZ). We considered 17 configurations of the DME dime for the orientation sampling. The calculated MP2 potential data were used to construct a 3‐site united atom force field model. The simulation results are compared with those using the empirical force field of Jorgensen and Ibrahim (Jorgensen and Ibrahim, J Am Chem Soc 1981, 103, 3976) and with available experimental measurements. We obtain quantitative agreements for the atom‐wise radial distribution functions, the self‐diffusion coefficients, and the shear viscosities over a wide range of experimental conditions. This force field thus provides a suitable starting point to predict liquid properties of DME from first principles intermolecular interactions with no empirical data input a priori. © 2012 Wiley Periodicals, Inc.     

Direct calculation of 1-octanol-water partition coefficients from adaptive biasing force molecular dynamics simulations

The 1-octanol-water partition coefficient log K(ow) of a solute is a key parameter used in the prediction of a wide variety of complex phenomena such as drug availability and bioaccumulation potential of trace contaminants. In this work, adaptive biasing force molecular dynamics simulations are used to determine absolute free energies of hydration, solvation, and 1-octanol-water partition coefficients for n-alkanes from methane to octane. Two approaches are evaluated; the direct transfer of the solute from 1-octanol to water phase, and separate transfers of the solute from the water or 1-octanol phase to vacuum, with both methods yielding statistically indistinguishable results. Calculations performed with the TIP4P and SPC∕E water models and the TraPPE united-atom force field for n-alkanes show that the choice of water model has a negligible effect on predicted free energies of transfer and partition coefficients for n-alkanes. A comparison of calculations using wet and dry octanol phases shows that the predictions for log K(ow) using wet octanol are 0.2-0.4 log units lower than for dry octanol, although this is within the statistical uncertainty of the calculation.     

Estimating entropies from molecular dynamics simulations

While the determination of free-energy differences by MD simulation has become a standard procedure for which many techniques have been developed, total entropies and entropy differences are still hardly ever computed. An overview of techniques to determine entropy differences is given, and the accuracy and convergence behavior of five methods based on thermodynamic integration and perturbation techniques was evaluated using liquid water as a test system. Reasonably accurate entropy differences are obtained through thermodynamic integration in which many copies of a solute are desolvated. When only one solute molecule is involved, only two methods seem to yield useful results, the calculation of solute-solvent entropy through thermodynamic integration, and the calculation of solvation entropy through the temperature derivative of the corresponding free-energy difference. One-step perturbation methods seem unsuitable to obtain entropy estimates.     

ReaxFF reactive force field for molecular dynamics simulations of hydrocarbon oxidation

To investigate the initial chemical events associated with high-temperature gas-phase oxidation of hydrocarbons, we have expanded the ReaxFF reactive force field training set to include additional transition states and chemical reactivity of systems relevant to these reactions and optimized the force field parameters against a quantum mechanics (QM)-based training set. To validate the ReaxFF potential obtained after parameter optimization, we performed a range of NVT-MD simulations on various hydrocarbon/O2 systems. From simulations on methane/O2, o-xylene/O2, propene/O2, and benzene/O2 mixtures, we found that ReaxFF obtains the correct reactivity trend (propene > o-xylene > methane > benzene), following the trend in the C-H bond strength in these hydrocarbons. We also tracked in detail the reactions during a complete oxidation of isolated methane, propene, and o-xylene to a CO/CO2/H2O mixture and found that the pathways predicted by ReaxFF are in agreement with chemical intuition and our QM results. We observed that the predominant initiation reaction for oxidation of methane, propene, and o-xylene under fuel lean conditions involved hydrogen abstraction of the methyl hydrogen by molecular oxygen forming hydroperoxyl and hydrocarbon radical species. While under fuel rich conditions with a mixture of these hydrocarbons, we observed different chemistry compared with the oxidation of isolated hydrocarbons including a change in the type of initiation reactions, which involved both decomposition of the hydrocarbon or attack by other radicals in the system. Since ReaxFF is capable of simulating complicated reaction pathways without any preconditioning, we believe that atomistic modeling with ReaxFF provides a useful method for determining the initial events of oxidation of hydrocarbons under extreme conditions and can enhance existing combustion models.     

The distributed diagonal force decomposition method for parallelizing molecular dynamics simulations

Parallelization is an effective way to reduce the computational time needed for molecular dynamics simulations. We describe a new parallelization method, the distributed-diagonal force decomposition method, with which we extend and improve the existing force decomposition methods. Our new method requires less data communication during molecular dynamics simulations than replicated data and current force decomposition methods, increasing the parallel efficiency. It also dynamically load-balances the processors' computational load throughout the simulation. The method is readily implemented in existing molecular dynamics codes and it has been incorporated into the CHARMM program, allowing its immediate use in conjunction with the many molecular dynamics simulation techniques that are already present in the program. We also present the design of the Force Decomposition Machine, a cluster of personal computers and networks that is tailored to running molecular dynamics simulations using the distributed diagonal force decomposition method. The design is expandable and provides various degrees of fault resilience. This approach is easily adaptable to computers with Graphics Processing Units because it is independent of the processor type being used.     

G protein- and agonist-bound serotonin 5-HT2A receptor model activated by steered molecular dynamics simulations

A 5-HT(2A) receptor model was constructed by homology modeling based on the β(2)-adrenergic receptor and the G protein-bound opsin crystal structures. The 5-HT(2A) receptor model was transferred into an active conformation by an agonist ligand and a G(αq) peptide in four subsequent steered molecular dynamics (MD) simulations. The driving force for the transformation was the addition of several known intermolecular and receptor interhelical hydrogen bonds enforcing the necessary helical and rotameric movements. Subsquent MD simulations without constraints confirmed the stability of the activated receptor model as well as revealed new information about stabilizing residues and bonds. The active 5-HT(2A) receptor model was further validated by retrospective ligand screening of more than 9400 compounds, whereof 182 were known ligands. The results show that the model can be used in drug discovery for virtual screening and structure-based ligand design as well as in GPCR activation studies.     

Enthalpies of mixing predicted using molecular dynamics simulations and OPLS force field

Enthalpy of mixing (EOM) is one of the most basic thermodynamic properties of mixtures. To assess feasibility of predicting EOM using force field simulation methods, fifteen (15) representative binary mixtures were investigated using MD simulations based on OPLS and TIP4P force fields. The simulation conditions and errors were carefully examined. The precision level of 0.04 kJ/mol was obtained for calculated EOM data. However, the predictions, measured by deviations from experimental data, were only qualitatively correct. Among various factors influencing the accuracy of predictions, force field quality representing interactions among different molecules plays the most significant role. Using methanol/benzene and ethanol/benzene as examples, we demonstrated that non-additive interaction terms between polarizable atoms can be used to significantly improve the quality of predictions. In addition, it appears that charge-dependent LJ parameters are required in order to represent the polarization effects accurately.     

Paramfit: Automated optimization of force field parameters for molecular dynamics simulations

The generation of bond, angle, and torsion parameters for classical molecular dynamics force fields typically requires fitting parameters such that classical properties such as energies and gradients match precalculated quantum data for structures that scan the value of interest. We present a program, Paramfit, distributed as part of the AmberTools software package that automates and extends this fitting process, allowing for simplified parameter generation for applications ranging from single molecules to entire force fields. Paramfit implements a novel combination of a genetic and simplex algorithm to find the optimal set of parameters that replicate either quantum energy or force data. The program allows for the derivation of multiple parameters simultaneously using significantly fewer quantum calculations than previous methods, and can also fit parameters across multiple molecules with applications to force field development. Paramfit has been applied successfully to systems with a sparse number of structures, and has already proven crucial in the development of the Assisted Model Building with Energy Refinement Lipid14 force field. © 2014 Wiley Periodicals, Inc.     

Ab initio based force field and molecular dynamics simulations of crystalline TATB

An all-atom force field for 1,3,5-triamino-2,4,6-trinitrobenzene (TATB) is presented. The classical intermolecular interaction potential for TATB is based on single-point energies determined from high-level ab initio calculations of TATB dimers. The newly developed potential function is used to examine bulk crystalline TATB via molecular dynamics simulations. The isobaric thermal expansion and isothermal compression under hydrostatic pressures obtained from the molecular dynamics simulations are in good agreement with experiment. The calculated volume-temperature expansion is almost one dimensional along the c crystallographic axis, whereas under compression, all three unit cell axes participate, albeit unequally.     

Comparative analysis of nucleotide translocation through protein nanopores using steered molecular dynamics and an adaptive biasing force

The translocation of nucleotide molecules across biological and synthetic nanopores has attracted attention as a next generation technique for sequencing DNA. Computer simulations have the ability to provide atomistic‐level insight into important states and processes, delivering a means to develop a fundamental understanding of the translocation event, for example, by extracting the free energy of the process. Even with current supercomputing facilities, the simulation of many‐atom systems in fine detail is limited to shorter timescales than the real events they attempt to recreate. This imposes the need for enhanced simulation techniques that expand the scope of investigation in a given timeframe. There are numerous free energy calculation and translocation methodologies available, and it is by no means clear which method is best applied to a particular problem. This article explores the use of two popular free energy calculation methodologies in a nucleotide‐nanopore translocation system, using the α‐hemolysin nanopore. The first uses constant velocity‐steered molecular dynamics (cv‐SMD) in conjunction with Jarzynski's equality. The second applies an adaptive biasing force (ABF), which has not previously been applied to the nucleotide‐nanpore system. The purpose of this study is to provide a comprehensive comparison of these methodologies, allowing for a detailed comparative assessment of the scientific merits, the computational cost, and the statistical quality of the data obtained from each technique. We find that the ABF method produces results that are closer to experimental measurements than those from cv‐SMD, whereas the net errors are smaller for the same computational cost. © 2014 The Authors Journal of Computational Chemistry Published by Wiley Periodicals, Inc.