Modeling solvent effects in molecular dynamics simulations of proteins

A series of computer simulations has been carried out on bovine pancreatic trypsin inhibitor using various models to mimic the effects of explicit bulk solvent on the structure of the protein. The solvent properties included are the polarization of the solute by the polar bulk solvent and the restraining effect on the motional freedom of the solute due to frictional drag at the solvent–protein surface interface. The former has been included by using a distance–dependent dielectric permittivity to screen the electrostatic interactions, whereas the latter is simulated by adding a limited number of solvent molecules near the protein surface. To achieve the proper mobility of the water molecules, their motion was restrained by adding a harmonic restraining force. It was found that a very small force constant was sufficient to model the static and dynamical behavior of the fully solvated solute, but that it was necessary to include enough explicit waters to occupy the first solvation shell. © 1992 John Wiley & Sons, Inc.     

Calculations of solute and solvent entropies from molecular dynamics simulations

The translational, rotational and conformational (vibrational) entropy contributions to ligand-receptor binding free energies are analyzed within the standard formulation of statistical thermodynamics. It is shown that the partitioning of the binding entropy into different components is to some extent arbitrary, but an appropriate method to calculate both translational and rotational entropy contributions to noncovalent association is by estimating the configurational volumes of the ligand in the bound and free states. Different approaches to calculating solute entropies using free energy perturbation calculations, configurational volumes based on root-mean-square fluctuations and covariance matrix based quasiharmonic analysis are illustrated for some simple molecular systems. Numerical examples for the different contributions demonstrate that theoretically derived results are well reproduced by the approximations. Calculation of solvent entropies, either using total potential energy averages or van't Hoff plots, are carried out for the case of ion solvation in water. Although convergence problems will persist for large and complex simulation systems, good agreement with experiment is obtained here for relative and absolute ion hydration entropies. We also outline how solvent and solute entropic contributions are taken into account in empirical binding free energy calculations using the linear interaction energy method. In particular it is shown that empirical scaling of the nonpolar intermolecular ligand interaction energy effectively takes into account size dependent contributions to the binding free energy.     

Comparison of a quantum random number generator with pseudorandom number generators for their use in molecular Monte Carlo simulations

Four pseudorandom number generators were compared with a physical, quantum‐based random number generator using the NIST suite of statistical tests, which only the quantum‐based random number generator could successfully pass. We then measured the effect of the five random number generators on various calculated properties in different Markov‐chain Monte Carlo simulations. Two types of systems were tested: conformational sampling of a small molecule in aqueous solution and liquid methanol under constant temperature and pressure. The results show that poor quality pseudorandom number generators produce results that deviate significantly from those obtained with the quantum‐based random number generator, particularly in the case of the small molecule in aqueous solution setup. In contrast, the widely used Mersenne Twister pseudorandom generator and a 64‐bit Linear Congruential Generator with a scrambler produce results that are statistically indistinguishable from those obtained with the quantum‐based random number generator. © 2017 Wiley Periodicals, Inc.     

Mixed monte carlo/molecular dynamics simulations in explicit solvent

A mixed Monte Carlo/Molecular Dynamics method using the trial moves for peptide backbone sampling known as Concerted Rotations with Angles was implemented. The algorithm was used to study polyalanine systems. Equivalent results to conventional Molecular Dynamics were obtained for simulations of Ala₆ in implicit solvent. To test the efficiency of the implemented method, several 150 ns simulations of Ala₁₂ in explicit water were performed. The results show that the present method yields significantly faster formation of secondary structure than the conventional Molecular Dynamics simulations. This opens the possibility to selectively sample alanine‐rich regions of larger peptides or proteins. It remains to be established whether hydrophilic amino acid residues can be successfully treated with the present methodology. © 2012 Wiley Periodicals, Inc.     

Simulation of DNA electrophoresis in systems of large number of solvent particles by coarse-grained hybrid molecular dynamics approach

Simulation of DNA electrophoresis facilitates the design of DNA separation devices. Various methods have been explored for simulating DNA electrophoresis and other processes using implicit and explicit solvent models. Explicit solvent models are highly desired but their applications may be limited by high computing cost in simulating large number of solvent particles. In this work, a coarse-grained hybrid molecular dynamics (CGH-MD) approach was introduced for simulating DNA electrophoresis in explicit solvent of large number of solvent particles. CGH-MD was tested in the simulation of a polymer solution and computation of nonuniform charge distribution in a cylindrical nanotube, which shows good agreement with observations and those of more rigorous computational methods at a significantly lower computing cost than other explicit-solvent methods. CGH-MD was further applied to the simulation of DNA electrophoresis in a polymer solution and in a well-studied nanofluidic device. Simulation results are consistent with observations and reported simulation results, suggesting that CGH-MD is potentially useful for studying electrophoresis of macromolecules and assemblies in nanofluidic, microfluidic, and microstructure array systems that involve extremely large number of solvent particles, nonuniformly distributed electrostatic interactions, bound and sequestered water molecules.     

Peptide loop-closure kinetics from microsecond molecular dynamics simulations in explicit solvent

End-to-end contact formation rates of several peptides were recently measured by tryptophan triplet quenching (Lapidus et al. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 7220). Motivated by these experiments, we study loop-closure kinetics for two peptides of different lengths, Cys-(Ala-Gly-Gln)n-Trp (n = 1, 2), in multiple all-atom explicit-solvent molecular dynamics simulations with different initial conditions and force fields. In 150 simulations of approximately 20 ns each, we collect data covering 1.0 and 0.8 micros for the penta-peptide simulated with the AMBER and CHARMM force fields, respectively, and about 0.5 micros each with the two force fields for the octa-peptide. These extensive simulations allow us to analyze the dynamics of peptides in the unfolded state with atomic resolution, thus probing early events in protein folding, and to compare molecular dynamics simulations directly with experiment. The calculated lifetimes of the tryptophan triplet state are in the range of 50-100 ns, in agreement with experimental measurements. However, end-to-end contacts form more rapidly, with characteristic times less than 10 ns. The contact formation rates for the two force fields are similar despite differences in the respective ensembles of peptide conformations.     

Reduced matrix of topological distances with a minimum number of independent parameters: distance vectors and molecular codes

Entries of the topological distance matrix are shown to be functions of entries of the reduced distance matrix, which has a smaller size. The latter entries are expressed through a minimum number m of independent parameters (m ≤ 2n − 3). The expanded distance matrices, whose sum constitutes the reduced matrix, are defined. The reduced vectors have lower degeneracy than the corresponding vectors proposed by Randić as molecular codes. Nondegenerate sets of reduced-matrix entries are proposed as molecular codes.     

Trees with the minimum Wiener number

The Wiener number (𝒲) of a connected graph is the sum of distances for all pairs of vertices. As a graphical invariant, it has been found extensive application in chemistry. Considering the family of trees with n vertices and a fixed maximum vertex degree, we derive some methods that can strictly reduce 𝒲 by shifting leaves. And then, by a process, we prove that the dendrimer on n vertices is the unique graph reaching the minimum Wiener number. © 2000 John Wiley & Sons, Inc. Int J Quant Chem 78: 331–340, 2000     

Energy conservation in molecular dynamics simulations of classical systems

Classical Newtonian dynamics is analytic and the energy of an isolated system is conserved. The energy of such a system, obtained by the discrete "Verlet" algorithm commonly used in molecular dynamics simulations, fluctuates but is conserved in the mean. This is explained by the existence of a "shadow Hamiltonian" H [S. Toxvaerd, Phys. Rev. E 50, 2271 (1994)], i.e., a Hamiltonian close to the original H with the property that the discrete positions of the Verlet algorithm for H lie on the analytic trajectories of H. The shadow Hamiltonian can be obtained from H by an asymptotic expansion in the time step length. Here we use the first non-trivial term in this expansion to obtain an improved estimate of the discrete values of the energy. The investigation is performed for a representative system with Lennard-Jones pair interactions. The simulations show that inclusion of this term reduces the standard deviation of the energy fluctuations by a factor of 100 for typical values of the time step length. Simulations further show that the energy is conserved for at least one hundred million time steps provided the potential and its first four derivatives are continuous at the cutoff. Finally, we show analytically as well as numerically that energy conservation is not sensitive to round-off errors.     

The solute-solvent system: solvent constraints on the conformational dynamics of acetylcholine

The objective of this study was to determine if and how a solvent influences internal motions in a solute molecule. Acetylcholine was chosen as the object of study given its interesting molecular structure and major biological significance. Molecular dynamics simulations were carried out in the vacuum (10 ns), water (5 ns), methanol (5 ns), and octanol (1.5 ns). Seven clusters of conformers were identified, namely, +g+g, -g-g, +gt, -gt, t+g, t-g, and tt, where the gauche and trans labels refer to the dihedral angles tau(2) and tau(3), respectively. As expected, the relative proportion of these conformational clusters was highly solvent-dependent and corresponded to a progressive loss of conformational freedom with increasing molecular weight of the solvent. More importantly, the conformational clusters were used to calculate instantaneous and median angular velocity (omega and omega(M), respectively) and instantaneous and median angular acceleration (alpha and alpha(M), respectively). Angular velocity and angular acceleration were both found to decrease markedly with increasing molecular weight of the solvent, i.e., vacuum (epsilon = 1) > water > methanol > octanol. The decrease from the vacuum to octanol was approximately 40% for tau(2) and approximately 60% for tau(3). Such solvent-dependent constraints on a solute's internal motions may be biologically and pharmacologically relevant.     

Treatment of multibody interactions in molecular simulations of systems with general bond networks

One of the most formidable difficulties in the computer programming of molecular simulations is the sometimes complicated bookkeeping required for keeping track of internal coordinates and their derivatives. A completely general method for keeping track of stretch (two-body), bend (three-body), and torsion, wag, and other four-body interactions for ANY bond network is presented. Computer code using this method for calculating internal coordinates and their derivatives can be used for completely different types of bond networks, no matter how complex, with little or no modification. The method is designed to incorporate recent improved formulas for calculating internal coordinates and their derivatives to ensure the most optimal calculation sequence. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 : 1513–1522, 1997     

Polymer brushes in solvents of variable quality: molecular dynamics simulations using explicit solvent

The structure and thermodynamic properties of a system of end-grafted flexible polymer chains grafted to a flat substrate and exposed to a solvent of variable quality are studied by molecular dynamics methods. The macromolecules are described by a coarse-grained bead-spring model, and the solvent molecules by pointlike particles, assuming Lennard-Jones-type interactions between pairs of monomers (epsilon(pp)), solvent molecules (epsilon(ss)), and solvent monomer (epsilon(ps)), respectively. Varying the grafting density sigma(g) and some of these energy parameters, we obtain density profiles of solvent particles and monomers, study structural properties of the chain (gyration radius components, bond orientational parameters, etc.), and examine also the profile of the lateral pressure P( parallel)(z), keeping in the simulation the normal pressure P( perpendicular) constant. From these data, the reduction of the surface tension between solvent and wall as a function of the grafting density of the brush has been obtained. Further results include the stretching force on the monomer adjacent to the grafting site and its variation with solvent quality and grafting density, and dynamic characteristics such as mobility profiles and chain relaxation times. Possible phase transitions (vertical phase separation of the solvent versus lateral segregation of the polymers into "clusters," etc.) are discussed, and a comparison to previous work using implicit solvent models is made. The variation of the brush height and the interfacial width of the transition zone between the pure solvent and the brush agrees qualitatively very well with corresponding experiments.     

Elongation of ordered peptide aggregate of an amyloidogenic hexapeptide NFGAIL observed in molecular dynamics simulations with explicit solvent

The mechanisms by which amyloidogenic peptides and proteins form soluble toxic oligomers remain elusive. We have studied the formation of partially ordered tetramers and well-ordered octamers of an amyloidogenic hexapeptide NFGAIL (residues 22-27 of the human islet amyloid polypeptide) in our previous work. Continuing the effort, we here probe the beta-sheet elongation process by a combined total of 2.0 micros molecular dynamics simulations with explicit solvent. In a set of 10 simulations with the peptides restrained to the extended conformation, we observed that the main growth mode was elongation along the beta-sheet hydrogen bonds through primarily a two-stage process. Driven by hydrophobic forces, the peptides initially attached to the surface of the ordered oligomer, moved quickly to the beta-sheet edges, and formed stable beta-sheet hydrogen bonds. Addition of peptides to the existing oligomer notably improved the order of the peptide aggregate in which labile outer layer beta-sheets were stabilized, which provides good templates for further elongation. These simulations suggested that elongation along the beta-sheet hydrogen bonds occurs at the intermediate stage when low-weight oligomers start to form. We did not observe significant preference toward either parallel or antiparallel beta-sheets at the elongation stage for this peptide. In another set of 10 unrestrained simulations, the dominant growth mode was disordered aggregation. Taken together, these results offered a glimpse at the molecular events leading to the formation of ordered and disordered low-weight oligomers.     

Improving efficiency of large time-scale molecular dynamics simulations of hydrogen-rich systems

A systematic analysis is performed on the effectiveness of removing degrees of freedom from hydrogen atoms and/or increasing hydrogen masses to increase the efficiency of molecular dynamics simulations of hydrogen-rich systems such as proteins in water. In proteins, high-frequency bond-angle vibrations involving hydrogen atoms limit the time step to 3 fs, which is already a factor of 1.5 beyond the commonly used time step of 2 fs. Removing these degrees of freedom from the system by constructing hydrogen atoms as dummy atoms, allows the time step to be increased to 7 fs, a factor of 3.5 compared with 2 fs. Additionally, a gain in simulation stability can be achieved by increasing the masses of hydrogen atoms with remaining degrees of freedom from 1 to 4 u. Increasing hydrogen mass without removing the high-frequency degrees of freedom allows the time step to be increased only to 4 fs, a factor of two, compared with 2 fs. The net gain in efficiency of sampling configurational space may be up to 15% lower than expected from the increase in time step due to the increase in viscosity and decrease in diffusion constant. In principle, introducing dummy atoms and increasing hydrogen mass do not influence thermodynamical properties of the system and dynamical properties are shown to be influenced only to a moderate degree. Comparing the maximum time step attainable with these methods (7 fs) to the time step of 2 fs that is routinely used in simulation, and taking into account the increase in viscosity and decrease in diffusion constant, we can say that a net gain in simulation efficiency of a factor of 3 to 3.5 can be achieved. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 786–798, 1999