Harmonic analysis of large systems. II. Comparison of different protein models

A series of normal mode analyses of bovine pancreatic trypsin inhibitor (BPTI) has been performed. The results of modifying the long-range truncation of electrostatics, reducing the conformational space of the system (reduced basis normal mode analysis), and using different parameter sets and models for the potential function are reported. Both explicit (904 atoms) and polar hydrogen (580 atoms) representations of BPTI were examined and produced nearly identical normal mode vectors but slightly modified vibrational frequencies. The truncation methods—no cutoff, shift, and switch—were examined, and the use of a short switching function was found to alter harmonic motion greatly. A table relating the different cutoff methods to several previously published frequencies for BPTI indicates that the diversity of published lowest frequencies is due to the use of different electrostatic models rather than to inherent differences in the models or energy parameters. Examining reduced basis results demonstrates that a dihedral basis yields similar normal mode vectors, though the vibrational frequencies are shifted to higher values. The analysis of BPTI harmonic dynamics using a spherical harmonic reduced basis set yields significantly altered dynamics, indicating that BPTI is not well represented as a homogeneous object at low temperatures. © 1995 John Wiley & Sons, Inc.

1 This article is a U.S. Government work and, as such, is in the public domain in the United States of America.

     

Harmonic analysis of large systems. I. Methodology

Methods have been developed for the determination of vibrational frequencies and normal modes of large systems in the full conformational space (including all degrees of freedom) and in a reduced conformational space (reducing the number of degrees of freedom). The computational method, which includes Hessian generation and storage, full and iterative diagonalization techniques, and the refinement of the results, is presented. A method is given for the quasiharmonic analysis and the reduced basis quasiharmonic analysis. The underlying principle is that from the atomic fluctuations, an effective harmonic force field can be determined relative to the dynamic average structure. Normal mode analysis tools can be used to characterize quasiharmonic modes of vibration. These correspond to conventional normal modes except that anharmonic effects are included. Numerous techniques for the analyses of vibrational frequencies and normal modes are described. Criteria for the analysis of the similarity of low-frequency normal modes is presented. The approach to determining the natural frequencies and normal modes of vibration described here is general and applicable to any large system. © 1995 John Wiley & Sons, Inc.

1 This article is a U.S. Government work and, as such, is in the public domain in the United States of America.

     

Solvent interaction energy calculations on molecular dynamics trajectories: increasing the efficiency using systematic frame selection

End-point methods such as linear interaction energy (LIE) analysis, molecular mechanics generalized Born solvent-accessible surface (MM/GBSA), and solvent interaction energy (SIE) analysis have become popular techniques to calculate the free energy associated with protein-ligand binding. Such methods typically use molecular dynamics (MD) simulations to generate an ensemble of protein structures that encompasses the bound and unbound states. The energy evaluation method (LIE, MM/GBSA, or SIE) is subsequently used to calculate the energy of each member of the ensemble, thus providing an estimate of the average free energy difference between the bound and unbound states. The workflow requiring both MD simulation and energy calculation for each frame and each trajectory proves to be computationally expensive. In an attempt to reduce the high computational cost associated with end-point methods, we study several methods by which frames may be intelligently selected from the MD simulation including clustering and address the question of how the number of selected frames influences the accuracy of the SIE calculations.     

Convergence of Free Energy Profile of Coumarin in Lipid Bilayer

Atomistic molecular dynamics (MD) simulations of druglike molecules embedded in lipid bilayers are of considerable interest as models for drug penetration and positioning in biological membranes. Here we analyze partitioning of coumarin in dioleoylphosphatidylcholine (DOPC) bilayer, based on both multiple, unbiased 3 μs MD simulations (total length) and free energy profiles along the bilayer normal calculated by biased MD simulations (～7 μs in total). The convergences in time of free energy profiles calculated by both umbrella sampling and z-constraint techniques are thoroughly analyzed. Two sets of starting structures are also considered, one from unbiased MD simulation and the other from "pulling" coumarin along the bilayer normal. The structures obtained by pulling simulation contain water defects on the lipid bilayer surface, while those acquired from unbiased simulation have no membrane defects. The free energy profiles converge more rapidly when starting frames from unbiased simulations are used. In addition, z-constraint simulation leads to more rapid convergence than umbrella sampling, due to quicker relaxation of membrane defects. Furthermore, we show that the choice of RESP, PRODRG, or Mulliken charges considerably affects the resulting free energy profile of our model drug along the bilayer normal. We recommend using z-constraint biased MD simulations based on starting geometries acquired from unbiased MD simulations for efficient calculation of convergent free energy profiles of druglike molecules along bilayer normals. The calculation of free energy profile should start with an unbiased simulation, though the polar molecules might need a slow pulling afterward. Results obtained with the recommended simulation protocol agree well with available experimental data for two coumarin derivatives.     

Crystal water dynamics of guanosine dihydrate: analysis of atomic displacement parameters, time profile of hydrogen-bonding probability, and translocation of water by MD simulation

The dynamics of crystal water molecules of guanosine dihydrate are investigated in detail by molecular dynamics (MD) simulation. A 2 ns simulation is performed using a periodic boundary box composed of 4 x 5 x 8 crystallographic unit cells and using the particle-mesh Ewald method for calculation of electrostatic energy. The simulated average atomic positions and atomic displacement parameters are remarkably coincident with the experimental values determined by X-ray analysis, confirming the high accuracy of this simulation. The dynamics of crystal water are analyzed in terms of atomic displacement parameters, orientation vectors, order parameters, self-correlation functions of the orientation vectors, time profiles of hydrogen-bonding probability, and translocations. The simulation clarifies that the average structure is composed of various stable and transient structures of the molecules. The simulated guanosine crystal forms a layered structure, with four water sites per asymmetric unit, classified as either interlayer water or intralayer water. From a detailed analysis of the translocations of water molecules in the simulation, columns of intralayer water molecules along the c axis appear to represent a pathway for hydration and dehydration by a kind of molecular valve mechanism.     

Ab initio protein structure prediction with force field parameters derived from water-phase quantum chemical calculation

Molecular dynamics (MD) simulations are extensively used in the study of the structures and functions of proteins. Ab initio protein structure prediction is one of the most important subjects in computational biology, and many trials have been performed using MD simulation so far. Since the results of MD simulations largely depend on the force field, reliable force field parameters are indispensable for the success of MD simulation. In this work, we have modified atom charges in a standard force field on the basis of water-phase quantum chemical calculations. The modified force field turned out appropriate for ab initio protein structure prediction by the MD simulation with the generalized Born method. Detailed analysis was performed in terms of the conformational stability of amino acid residues, the stability of secondary structure of proteins, and the accuracy for prediction of protein tertiary structure, comparing the modified force field with a standard one. The energy balance between alpha-helix and beta-sheet structures was significantly improved by the modification of charge parameters.     

Simulations of vibrational spectra from classical trajectories: calibration with ab initio force fields

An algorithm allowing simulating vibrational spectra from classical time-dependent trajectories was applied for infrared absorption, vibrational circular dichroism, Raman, and Raman optical activity of model harmonic systems. The implementation of the theory within the TINKER molecular dynamics (MD) program package was tested with ab initio harmonic force fields in order to determine the feasibility for more extended MD simulations. The results suggest that sufficiently accurate frequencies can be simulated with integration time steps shorter than about 0.5 fs. For a given integration time step, lower vibrational frequencies ( approximately 0-2000 cm(-1)) could be reproduced with a higher accuracy than higher-frequency vibrational modes (e.g., O-H and C-H stretching). In principle, the algorithm also provides correct intensities for ideal systems. In applied simulations, however, the intensity profiles are affected by an unrealistic energy distribution between normal modes and a slow energy relaxation. Additionally, the energy fluctuations may cause weakening of the intensities on average. For ab initio force fields, these obstacles could be overcome by an arbitrary normal mode energy correction. For general MD simulations, averaging of many shorter MD trajectories started with randomly distributed atomic velocities provided the best spectral shapes. alpha-pinene, D-gluconic acid, formaldehyde dimer, and the acetylprolineamide molecule were used in the tests.     

Infrared absorption line shapes in the classical limit: a comparison of the classical dipole and fluctuating frequency approximations

Infrared spectroscopy is a versatile technique for probing the structure and dynamics of condensed-phase systems. Simulating infrared absorption spectra with molecular dynamics (MD) offers a powerful means to establish a molecular-level interpretation of experimental results, as well as a basis for the parametrization of more accurate simulation force-fields. Two distinct methods for the calculation of infrared absorption line shapes of high-frequency (Planck's omega/k(B)T>1) vibrational probes from MD simulations are examined: The classical dipole approximation (CDA) and the fluctuating frequency approximation (FFA). Although these two formalisms result in expressions for the infrared absorption line shape that appear very different, both approximations are shown to yield identical results for the infrared line shape of a harmonic system in the condensed-phase. The equivalence of the FFA and CDA is also demonstrated in the case where the transition dipole of the oscillator fluctuates in response to the environment (i.e., where the Condon approximation has been relaxed). Finally we examine the effects of solute anharmonicity and demonstrate that the CDA and FFA are not equivalent in general, and the magnitude of the deviations increases with anharmonicity. We conclude that the calculation of infrared absorption line shapes via the CDA is a promising alternative to the FFA approach in cases where it may be difficult or undesirable to employ the latter, particularly when the effects of anharmonicity are small.     

On the calculation of velocity-dependent properties in molecular dynamics simulations using the leapfrog integration algorithm

Widely used programs for molecular dynamics simulation of (bio)molecular systems are the Verlet and leapfrog algorithms. In these algorithms, the particle velocities are less accurately propagated than the positions. Important quantities for the simulation such as the temperature and the pressure involve the squared velocities at full time steps. Here, we derive an expression for the squared particle velocity at full time step in the leapfrog scheme, which is more accurate than the standardly used one. In particular, this allows us to show that the full time step kinetic energy of a particle is more accurately computed as the average of the kinetic energies at previous and following half steps than as the square of the average velocity as implemented in various molecular dynamics codes. Use of the square of the average velocity introduces a systematic bias in the calculation of the instantaneous temperature and pressure of a molecular dynamics system. We show the consequences when the system is coupled to a thermostat and a barostat.     

Mechanism of electron and hole localization in poly(dimethylsilane) radical ions

The mechanism of excess electron and hole localizations in radical ions of poly(dimethylsilane) (PDMS) has been investigated by means of molecular dynamics (MD) and extended Hückel methods. Oligo(dimethylsilane) composed of 100 monomer units of dimethylsilane, CH3(Si(CH3)2)(n)CH3 (n = 100), were used as a model of PDMS. Both wings of the oligomer were capped by a methyl group. First, the geometry of PDMS with a regular all-trans form was fully optimized by MM2+ energy gradient method. Next, the MD calculation was carried out for PDMS at 300 K. The structure of PDMS was gradually deformed as a function of simulation time, especially the dihedral angle of Si-Si-Si-Si backbone that was randomized. At time zero when the structure has the regular all-trans form, both the excess electron and hole were completely delocalized on the Si backbone of PDMS. After thermal activation, the localization of the electron and hole was found. The mechanism of the localization was discussed on the basis of theoretical results.     

Dcdftbmd: Divide-and-Conquer Density Functional Tight-Binding Program for Huge-System Quantum Mechanical Molecular Dynamics Simulations

Dcdftbmd is a Fortran 90/95 program that enables efficient quantum mechanical molecular dynamics (MD) simulations using divide-and-conquer density functional tight-binding (DC-DFTB) method. Based on the remarkable performance of previous massively parallel DC-DFTB energy and gradient calculations for huge systems, the code has been specialized to MD simulations. Recent implementations and modifications including DFTB extensions, improved computational speed in the DC-DFTB computational steps, algorithms for efficient initial guess charge prediction, and free energy calculations via metadynamics technique have enhanced the capability to obtain atomistic insights in novel applications to nanomaterials and biomolecules. The energy, structure, and other molecular properties are also accessible through the single-point calculation, geometry optimization, and vibrational frequency analysis. The available functionalities are outlined together with efficiency tests and simulation examples. © 2019 Wiley Periodicals, Inc.     

Study of protein fluctuation with an effective inter-Cα atomic potential derived from average distances between amino acids in proteins

The conformational dynamics around the native structure of bovine pancreatic trypsin inhibitor (BPTI) in both the oxidized and reduced forms was analyzed by a Monte Carlo method using an approximate residue–residue potential derived from the statistics of average distances between C^α atoms of residues as proposed by the present author (T. Kikuchi, J Comput Chem 1996, 17, 226–237). The results from using this effective potential are similar to those from molecular dynamics simulations, taking all atoms into account, and are consistent with temperature factors from an X-ray analysis and disulfide formation from a kinetic experiment. This agreement suggests that the essential nature of the potential energy surface formed by the potential around the native structure closely mimics the actual energy landscape within the resolution of C^α atomic fluctuation. Furthermore, it is expected that the potential we found can describe the basic properties of folding kinetics. Examination of the fluctuation property of the native structure of BPTI threaded by a sequence from cytochrome b562 reveals differences specific to the sequence and this result also shows that the dynamical properties obtained in our calculations are not only ascribed to the geometrical constraints of the initial conformation but also the force field specifically produced by a sequence. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 713–719, 1999     

离子液体的量化计算及分子动力学模拟研究进展

离子液体是由有机阳离子和无机/有机阴离子构成的盐类,一般在室温或接近于室温下呈液态,因此常被称为室温离子液体(RTIL).依据不同的划分标准,离子液体有多种分类方式:根据年代的不同可将离子液体分为第一代、第二代及第三代离子液体,例如:烷基咪唑和烷基吡啶的金属卤化物盐等[1];根据阳离子的不同可将离子液体分为季鏻     

Molecular modeling and dynamics of neuropeptide Y

Summary A combination of molecular modeling and molecular dynamics (MD) is used to determine a theoretical structure for neuropeptide Y (NPY). Starting with the X-ray structure for avian pancreatic polypeptide (APP), the substituted amino acids were mutated, the side chains oriented to local potential energy minima, and the entire structure minimized and subjected to an MD simulation. Comparison of the resulting NPY structure with APP X-ray and MD results showed secondary structural elements to be maintained and RMS fluctuations to be similar, although differences in both were observed. The approach presented offers a means to study the structure-function relationships of NPY and other similar polypeptides when combined with pharmacological measurements.Abbreviations NPY Neuropeptide Y - APP Avian pancreatic polypeptide - ABNR Adopted-basis Newton Raphson - MD Molecular dynamics     

Langevin model of the temperature and hydration dependence of protein vibrational dynamics

The modification of internal vibrational modes in a protein due to intraprotein anharmonicity and solvation effects is determined by performing molecular dynamics (MD) simulations of myoglobin, analyzing them using a Langevin model of the vibrational dynamics and comparing the Langevin results to a harmonic, normal mode model of the protein in vacuum. The diagonal and off-diagonal Langevin friction matrix elements, which model the roughness of the vibrational potential energy surfaces, are determined together with the vibrational potentials of mean force from the MD trajectories at 120 K and 300 K in vacuum and in solution. The frictional properties are found to be describable using simple phenomenological functions of the mode frequency, the accessible surface area, and the intraprotein interaction (the displacement vector overlap of any given mode with the other modes in the protein). The frictional damping of a vibrational mode in vacuum is found to be directly proportional to the intraprotein interaction of the mode, whereas in solution, the friction is proportional to the accessible surface area of the mode. In vacuum, the MD frequencies are lower than those of the normal modes, indicating intramolecular anharmonic broadening of the associated potential energy surfaces. Solvation has the opposite effect, increasing the large-amplitude vibrational frequencies relative to in vacuum and thus vibrationally confining the protein atoms. Frictional damping of the low-frequency modes is highly frequency dependent. In contrast to the damping effect of the solvent, the vibrational frequency increase due to solvation is relatively temperature independent, indicating that it is primarily a structural effect. The MD-derived vibrational dynamic structure factor and density of states are well reproduced by a model in which the Langevin friction and potential of mean force parameters are applied to the harmonic normal modes.     

An efficient method for computing excess free energy of liquid

We present a new theoretical method for efficient calculation of free energy of liquid. This interaction entropy method allows one to compute entropy and free energy of liquid from standard single step MD (molecular dynamics) simulation directly in liquid state without the need to perform MD simulations at many intermediate states as required in thermodynamic integration or free energy perturbation methods. In this new approach, one only needs to evaluate the interaction energy of a single (fixed) liquid molecule with the rest of liquid molecules as a function of time from a standard MD simulation of liquid and the fluctuation of distribution of this interaction energy is then used to calculate the interaction entropy of the liquid. Explicit theoretical derivation of this interaction entropy approach is provided and numerical calculations for the benchmark liquid water system were carried out using three different water models. Numerical analysis of the result was performed and comparison of the computational result with experimental data and other theoretical results were provided. Excellent agreement of calculated free energies with the experimental data using TIP4P model is obtained for liquid water.     

Molecular Dynamics Simulation of Glass Transition Behavior of Polyimide Ensemble

The prediction of glass transition temperature from chemical structure has a great significance to select and design new high-properties materials. However, for the estimation and correlation methods, the deficiency of parameters for newer groups will lead to invalidity of Tg prediction or greater deviation from experiment. In the present work, we predicted Tg for a polyimide (PI) ensemble with rigid moieties, and analyzed structural factor that regards to the rotation barrier of the bridging…     

Solvophobic and steric effects of side groups on polymer folding: molecular modeling studies of amine-functionalized m-poly(phenyleneethynylene) foldamers in aqueous solution

The folding behavior of five different amine-functionalized m-poly(phenyleneethynylene) (m-PPE) oligomers containing 24 phenyl rings (12 residues, where a residue includes 2 phenyl rings) in water was examined by using a combination of molecular dynamics (MD) and replica exchange molecular dynamics (REMD) simulation techniques. The REMD method employed the highly parallelized GROMACS MD software and a modified OPLS-AA force field to simulate 44 replicas of each solvated system in parallel, with temperatures ranging from 300 to 577 K. Our results showed that the REMD method was more effective in predicting the helical conformation of the m-PPE in water, from an extended structure, than canonical MD methods in the same simulation time. Furthermore, we observed from canonical MD simulations of the explicitly solvated helical m-PPEs at 300 K that the radius of gyration, average helix inner diameter, and average helix pitch of the helical structure all pass through a minima when the side group is R = OC(2)H(5) as R is changed from R = H through OC(4)H(9).