Quasi-Hamiltonian equations of motion for internal coordinate molecular dynamics of polymers

Conventional molecular dynamics simulations of macromolecules require long computational times because the most interesting motions are very slow compared to the fast oscillations of bond lengths and bond angles that limit the integration time step. Simulation of dynamics in the space of internal coordinates, that is, with bond lengths, bond angles, and torsions as independent variables, gives a theoretical possibility of eliminating all uninteresting fast degrees of freedom from the system. This article presents a new method for internal coordinate molecular dynamics simulations of macromolecules. Equations of motion are derived that are applicable to branched chain molecules with any number of internal degrees of freedom. Equations use the canonical variables and they are much simpler than existing analogs. In the numerical tests the internal coordinate dynamics are compared with the traditional Cartesian coordinate molecular dynamics in simulations of a 56 residue globular protein. For the first time it was possible to compare the two alternative methods on identical molecular models in conventional quality tests. It is shown that the traditional and internal coordinate dynamics require the same time step size for the same accuracy and that in the standard geometry approximation of amino acids, that is, with fixed bond lengths, bond angles, and rigid aromatic groups, the characteristic step size is 4 fs, which is 2 times higher than with fixed bond lengths only. The step size can be increased up to 11 fs when rotation of hydrogen atoms is suppressed. © 1997 by John Wiley & Sons, Inc. J Comput Chem 18 : 1354–1364, 1997     

Backbone motions of free and pheromone-bound major urinary protein I studied by molecular dynamics simulation

Molecular motions of free and pheromone-bound mouse major urinary protein I, previously investigated by NMR relaxation, were simulated in 30 ns molecular dynamics (MD) runs. The backbone flexibility was described in terms of order parameters and correlation times, commonly used in the NMR relaxation analysis. Special attention was paid to the effect of conformational changes on the nanosecond time scale. Time-dependent order parameters were determined in order to separate motions occurring on different time scales. As an alternative approach, slow conformational changes were identified from the backbone torsion angle variances, and "conformationally filtered" order parameters were calculated for well-defined conformation states. A comparison of the data obtained for the free and pheromone-bound protein showed that some residues are more rigid in the bound form, but a larger portion of the protein becomes more flexible upon the pheromone binding. This finding is in general agreement with the NMR results. The higher flexibility observed on the fast (fs-ps) time scale was typically observed for the residues exhibiting higher conformational freedom on the ns time scale. An inspection of the hydrogen bond network provided a structural explanation for the flexibility differences between the free and pheromone-bound proteins in the simulations.     

MC(JBW): Simple but smart Monte Carlo algorithm for free energy simulations of multiconformational molecules

Many of the most common molecular simulation methods, including Monte Carlo (MC) and molecular or stochastic dynamics (MD or SD), have significant difficulties in sampling the space of molecular potential energy surfaces characterized by multiple conformational minima and significant energy barriers. In such cases improved sampling can be obtained by special techniques that lower such barriers or somehow direct search steps toward different low energy regions of space. We recently described a hybrid MC/SD algorithm [MC(JBW)/SD] incorporating such a technique that directed MC moves of selected torsion and bond angles toward known low energy regions of conformational space. Exploration of other degrees of freedom was left to the SD part of the hybrid algorithm. In the work described here, we develop a related but simpler simulation algorithm that uses only MC to sample all degrees of freedom (e.g., stretch, bend, and torsion). We term this algorithm MC(JBW). Using simulations on various model potential energy surfaces and on simple molecular systems (n-pentane, n-butane, and cyclohexane), MC(JBW) is shown to generate ensembles of states that are indistinguishable from the canonical ensembles generated by classical Metropolis MC in the limit of very long simulations. We further demonstrate the utility of MC(JBW) by evaluating the room temperature free energy differences between conformers of various substituted cyclohexanes and the larger ring hydrocarbons cycloheptane, cyclooctane, cyclononane, and cyclodecane. The results compare favorably with available experimental data and results from previously reported MC(JBW)/SD conformational free energy calculations. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1736–1745, 1998     

Conformational dynamics of polypeptides and proteins in the dihedral angle space and in the cartesian coordinate space: Normal mode analysis of deca-alanine

Effects of different treatments of the degrees of freedom of bond length stretching and bond angle bending in computational analysis of conformational dynamics of proteins and polypeptides are assessed. More specifically, the normal mode analysis of conformational dynamics of α-helix of deca-alanine has been carried out both in the dihedral angle space (DAS) and in the Cartesian coordinate space (CCS). Almost perfect one-to-one correspondence has been found between normal modes in the CCS with frequencies less than 128 cm^?1 and those in the DAS with frequencies less than 164 cm^?1. Patterns of atomic displacements in the corresponding modes are very similar. This indicates that the effects of fixing degrees of freedom of bond length stretching and bond angle bending on the very-low-frequency normal modes in the CCS with frequencies less than 128 cm^?1 are almost solely to increase the frequencies by about 20%. The conclusion indicates that the different treatment of these degree does not lead to qualitatively different results as long as low-frequency motions are concerned. Based on the results of calculation, mechanical property of the α-helix of deca-alanine is discussed.     

Using internal and collective variables in Monte Carlo simulations of nucleic acid structures: chain breakage/closure algorithm and associated Jacobians

This article describes a method for solving the geometric closure problem for simplified models of nucleic acid structures by using the constant bond lengths approximation. The resulting chain breakage/closure equations, formulated in the space of variable torsion and bond angles, are easy to solve, and have only two solutions. The analytical simplicity is in contrast with the high complexity of the closure problem in the torsion angle space with at most 16 solutions, which has been dealt with by several authors and was solved analytically by Wu and Deem (J. Chem. Phys. 1999, 111, 6625). The discussion on the choice of variables and associated Jacobians is focussed on the question of how conformational equilibration is affected in Monte Carlo simulations of molecular systems. In addition to the closure of the phosphate backbone, it is necessary to also solve the closure problem for the five-membered flexible furanose sugar ring. Explicit closure equations and the resulting Jacobians are given both for the complete four-variable model of the furanose ring and simulations in the phase-amplitude space of the five-membered ring, which are based on the approximate two-variable model of furanose introduced by Gabb et al. (J. Comput. Chem. 1995, 16, 667). The suggested closure algorithm can be combined with collective variables defined by translations and rotations of the monomeric nucleotide units. In comparison with simple internal coordinate moves, the resulting concerted moves describe local structural changes that have high acceptance rates and enable fast conformational equilibration. Appropriate molecular models are put forward for prospective Monte Carlo simulations of nucleic acids, but can be easily adapted to other biomolecular systems, such as proteins and lipid structures in biological membranes.     

PDB ligand conformational energies calculated quantum-mechanically

We present here a greatly updated version of an earlier study on the conformational energies of protein-ligand complexes in the Protein Data Bank (PDB) [Nicklaus et al. Bioorg. Med. Chem. 1995, 3, 411-428], with the goal of improving on all possible aspects such as number and selection of ligand instances, energy calculations performed, and additional analyses conducted. Starting from about 357,000 ligand instances deposited in the 2008 version of the Ligand Expo database of the experimental 3D coordinates of all small-molecule instances in the PDB, we created a "high-quality" subset of ligand instances by various filtering steps including application of crystallographic quality criteria and structural unambiguousness. Submission of 640 Gaussian 03 jobs yielded a set of about 415 successfully concluded runs. We used a stepwise optimization of internal degrees of freedom at the DFT level of theory with the B3LYP/6-31G(d) basis set and a single-point energy calculation at B3LYP/6-311++G(3df,2p) after each round of (partial) optimization to separate energy changes due to bond length stretches vs bond angle changes vs torsion changes. Even for the most "conservative" choice of all the possible conformational energies-the energy difference between the conformation in which all internal degrees of freedom except torsions have been optimized and the fully optimized conformer-significant energy values were found. The range of 0 to ~25 kcal/mol was populated quite evenly and independently of the crystallographic resolution. A smaller number of "outliers" of yet higher energies were seen only at resolutions above 1.3 ?. The energies showed some correlation with molecular size and flexibility but not with crystallographic quality metrics such as the Cruickshank diffraction-component precision index (DPI) and R(free)-R, or with the ligand instance-specific metrics such as occupancy-weighted B-factor (OWAB), real-space R factor (RSR), and real-space correlation coefficient (RSCC). We repeated these calculations with the solvent model IEFPCM, which yielded energy differences that were generally somewhat lower than the corresponding vacuum results but did not produce a qualitatively different picture. Torsional sampling around the crystal conformation at the molecular mechanics level using the MMFF94s force field typically led to an increase in energy.     

ICFF: a new method to incorporate implicit flexibility into an internal coordinate force field

We introduce a new method to accurately "project" a Cartesian force field onto an internal coordinate molecular model with fixed-bond geometry. The algorithm automatically generates the Internal Coordinate Force Field (ICFF), which is a close approximation of the "source" Cartesian force field. The ICFF method reduces the number of free variables in a model by at least 10-fold and facilitates the fast convergence of geometry optimizations, an advantage that is critical for many applications such as the docking of flexible ligands or conformational modeling of macromolecules. Although covalent geometry is fixed in an ICFF model, implicit flexibility is incorporated into the force field parameters in the following two ways. First, we formulate an empirical torsion energy term in ICFF as a sixfold Fourier series and develop a procedure to calculate the Fourier coefficients from the conformational energy profiles of the fully flexible Cartesian model. The ICFF torsion parameters thus represent not only torsion component of the source force field, but also bond bending, bond stretching, and "1-4" van der Waals interactions. Second, we use a soft polynomial repulsion function for "1-5" and "1-6" interactions to mimic the flexibility of bonds, connecting these atoms. Also, we suggest a way to use a local part of the Cartesian force field to automatically generate fixed covalent geometries, compatible with the ICFF energy function. Here, we present an implementation of the ICFF algorithm, which employs the MMFF94s Cartesian force field as a "source." Extensive benchmarking of ICFF with a representative set of organic molecules demonstrates that the implicit flexibility model accurately reproduces MMFF94s equilibrium conformational energy differences (RMSD approximately 0.64 kcal) and, most importantly, detailed torsion energy profiles (RMSD approximately 0.37 kcal). This accuracy is characteristic of the method, because all the ICFF parameters (except one scaling factor in the "1-5,1-6" repulsion term) are derived directly from the source Cartesian force field and do not depend on any particular molecular set. In contrast, the rigid geometry model with the MMFF94s energy function yields highly biased estimations in this test with the RMSD exceeding 1.2 kcal for the equilibrium energy comparisons and approximately 3.4 kcal for the torsion energy profiles.     

Molecular dynamics coupled with a virtual system for effective conformational sampling

An enhanced conformational sampling method is proposed: virtual‐system coupled canonical molecular dynamics (VcMD). Although VcMD enhances sampling along a reaction coordinate, this method is free from estimation of a canonical distribution function along the reaction coordinate. This method introduces a virtual system that does not necessarily obey a physical law. To enhance sampling the virtual system couples with a molecular system to be studied. Resultant snapshots produce a canonical ensemble. This method was applied to a system consisting of two short peptides in an explicit solvent. Conventional molecular dynamics simulation, which is ten times longer than VcMD, was performed along with adaptive umbrella sampling. Free‐energy landscapes computed from the three simulations mutually converged well. The VcMD provided quicker association/dissociation motions of peptides than the conventional molecular dynamics did. The VcMD method is applicable to various complicated systems because of its methodological simplicity. © 2018 Wiley Periodicals, Inc.     

Rotamer decomposition and protein dynamics: Efficiently analyzing dihedral populations from molecular dynamics

Molecular mechanics methods have matured into powerful methods to understand the dynamics and flexibility of macromolecules and especially proteins. As multinanosecond to microsecond length molecular dynamics (MD) simulations become commonplace, advanced analysis tools are required to generate scientifically useful information from large amounts of data. Some of the key degrees of freedom to understand protein flexibility and dynamics are the amino acid residue side chain dihedral angles. In this work, we present an easily automated way to summarize and understand the relevant dihedral populations. A tremendous reduction in complexity is achieved by describing dihedral timeseries in terms of histograms decomposed into Gaussians. Using the familiar and widely studied protein lysozyme, it is demonstrated that our approach captures essential properties of protein structure and dynamics. A simple classification scheme is proposed that indicates the rotational state population for each dihedral angle of interest and allows a decision if a given side chain or peptide backbone fragment remains rigid during the course of an MD simulation, adopts a converged distribution between conformational substates or has not reached convergence yet. © 2012 Wiley Periodicals, Inc.     

Conformational flexibility of the pentasaccharide LNF-2 deduced from NMR spectroscopy and molecular dynamics simulations

Human milk oligosaccharides (HMOs) are important as prebiotics since they stimulate the growth of beneficial bacteria in the intestine and act as receptor analogues that can inhibit the binding of pathogens. The conformation and dynamics of the HMO Lacto-N-fucopentaose 2 (LNF-2), α-L-Fucp-(1 → 4)[β-D-Galp-(1 → 3)]-β-D-GlcpNAc-(1 → 3)-β-D-Galp-(1 → 4)-D-Glcp, having a Lewis A epitope, has been investigated employing NMR spectroscopy and molecular dynamics (MD) computer simulations. 1D (1)H,(1)H-NOESY experiments were used to obtain proton-proton cross-relaxation rates from which effective distances were deduced and 2D J-HMBC and 1D long-range experiments were utilized to measure trans-glycosidic (3)J(CH) coupling constants. The MD simulations using the PARM22/SU01 force field for carbohydrates were carried out for 600 ns with explicit water as solvent which resulted in excellent sampling for flexible glycosidic torsion angles. In addition, in vacuo MD simulations were performed using an MM3-2000 force field, but the agreement was less satisfactory based on an analysis of heteronuclear trans-glycosidic coupling constants. LNF-2 has a conformationally well-defined region consisting of the terminal branched part of the pentasaccharide, i.e., the Lewis A epitope, and a flexible β-D-GlcpNAc-(1 → 3)-β-D-Galp-linkage towards the lactose unit, which is situated at the reducing end. For this β-(1 → 3)-linkage a negative ψ torsion angle is favored, when experimental NMR data is combined with the MD simulation in the analysis. In addition, flexibility on a similar time scale, i.e., on the order of the global overall molecular reorientation, may also be present for the ? torsion angle of the β-D-Galp-(1 → 4)-D-Glcp-linkage as suggested by the simulation. It was further observed from a temperature variation study that some (1)H NMR chemical shifts of LNF-2 were highly sensitive and this study indicates that Δδ/ΔT may be an additional tool for revealing conformational dynamics of oligosaccharides.     

Statistical thermodynamics of internal rotation in a hindering potential of mean force obtained from computer simulations

A method of statistical estimation is applied to the problem of one-dimensional internal rotation in a hindering potential of mean force. The hindering potential, which may have a completely general shape, is expanded in a Fourier series, the coefficients of which are estimated by fitting an appropriate statistical-mechanical distribution to the random variable of internal rotation angle. The function of reduced moment of inertia of an internal rotation is averaged over the thermodynamic ensemble of atomic configurations of the molecule obtained in stochastic simulations. When quantum effects are not important, an accurate estimate of the absolute internal rotation entropy of a molecule with a single rotatable bond is obtained. When there is more than one rotatable bond, the "marginal" statistical-mechanical properties corresponding to a given internal rotational degree of freedom are reduced. The method is illustrated using Monte Carlo simulations of two public health relevant halocarbon molecules, each having a single internal-rotation degree of freedom, and a molecular dynamics simulation of an immunologically relevant polypeptide, in which several dihedral angles are analyzed.     

A multistate empirical valence bond description of protonatable amino acids

The multistate empirical valence bond (MS-EVB) model, which was developed for molecular dynamics simulations of proton transport in water and biomolecular systems, is extended for the modeling of protonatable amino acid residues in aqueous environments, specifically histidine and glutamic acid. The parameters of the MS-EVB force field are first determined to reproduce the geometries and energetics of the gas phase amino acid-water clusters. These parameters are then optimized to reproduce experimental pK(a) values. The free energy profiles for acid ionization and the corresponding pK(a) values are calculated by MS-EVB molecular dynamics simulations utilizing the umbrella sampling technique, with the center of excess charge coordinate chosen as the dissociation reaction coordinate. A general procedure for fitting the MS-EVB parameters is formulated, which allows for the parametrization of other amino acid residues with protonatable groups and the subsequent use of the MS-EVB approach for molecular dynamics simulations of proton transfer processes in proteins involving protonation/deprotonation of the protonatable amino acid groups.     

A Method to Explore Protein Side Chain Conformational Variability Using Experimental Data

Experimentally measured values of molecular properties or observables of biomolecules such as proteins are generally averages over time and space, which do not contain su?cient information to determine the underlying conformational distribution of the molecules in solution. The relationship between experimentally measured NMR ³J‐coupling values and the corresponding dihedral angle values is a particularly complicated case due to its nonlinear, multiple‐valued nature. Molecular dynamics (MD) simulations at constant temperature can generate Boltzmann ensembles of molecular structures that are free from a priori assumptions about the nature of the underlying conformational distribution. They suffer, however, from limited sampling with respect to time and conformational space. Moreover, the quality of the obtained structures is dependent on the choice of force ?eld and solvation model. A recently proposed method that uses time‐averaging with local‐elevation (LE) biasing of the conformational search provides an elegant means of overcoming these three problems. Using a set of side chain ³J‐coupling values for the FK506 binding protein (FKBP), we ?rst investigate the uncertainty in the angle values predicted theoretically. We then propose a simple MD‐based technique to detect inconsistencies within an experimental data set and identify degrees of freedom for which conformational averaging takes place or for which force ?eld parameters may be de?cient. Finally, we show that LE MD is the best method for producing ensembles of structures that, on average, ?t the experimental data.     

Force‐field parameters of the Ψ and Φ around glycosidic bonds to oxygen and sulfur atoms

The Ψ and Φ torsion angles around glycosidic bonds in a glycoside chain are the most important determinants of the conformation of a glycoside chain. We determined force‐field parameters for Ψ and Φ torsion angles around a glycosidic bond bridged by a sulfur atom, as well as a bond bridged by an oxygen atom as a preparation for the next study, i.e., molecular dynamics free energy calculations for protein‐sugar and protein‐inhibitor complexes. First, we extracted the Ψ or Φ torsion energy component from a quantum mechanics (QM) total energy by subtracting all the molecular mechanics (MM) force‐field components except for the Ψ or Φ torsion angle. The Ψ and Φ energy components extracted (hereafter called “the remaining energy components”) were calculated for simple sugar models and plotted as functions of the Ψ and Φ angles. The remaining energy component curves of Ψ and Φ were well represented by the torsion force‐field functions consisting of four and three cosine functions, respectively. To confirm the reliability of the force‐field parameters and to confirm its compatibility with other force‐fields, we calculated adiabatic potential curves as functions of Ψ and Φ for the model glycosides by adopting the Ψ and Φ force‐field parameters obtained and by energetically optimizing other degrees of freedom. The MM potential energy curves obtained for Ψ and Φ well represented the QM adiabatic curves and also these curves' differences with regard to the glycosidic oxygen and sulfur atoms. Our Ψ and Φ force‐fields of glycosidic oxygen gave MM potential energy curves that more closely represented the respective QM curves than did those of the recently developed GLYCAM force‐field. © 2009 Wiley Periodicals, Inc., J Comput Chem, 2009     

Using operators to expand the block matrices forming the Hessian of a molecular potential

We derive compact expressions of the second‐order derivatives of bond length, bond angle, and proper and improper torsion angle potentials, in terms of operators represented in two orthonormal bases. Hereby, simple rules to generate the Hessian of an internal coordinate or a molecular potential can be formulated. The algorithms we provide can be implemented efficiently in high‐level programming languages using vectorization. Finally, the method leads to compact expressions for a second‐order expansion of an internal coordinate or a molecular potential. © 2014 Wiley Periodicals, Inc.     

Virtual‐system‐coupled adaptive umbrella sampling to compute free‐energy landscape for flexible molecular docking

A novel enhanced conformational sampling method, virtual‐system‐coupled adaptive umbrella sampling (V‐AUS), was proposed to compute 300‐K free‐energy landscape for flexible molecular docking, where a virtual degrees of freedom was introduced to control the sampling. This degree of freedom interacts with the biomolecular system. V‐AUS was applied to complex formation of two disordered amyloid‐β (Aβ_30–35) peptides in a periodic box filled by an explicit solvent. An interpeptide distance was defined as the reaction coordinate, along which sampling was enhanced. A uniform conformational distribution was obtained covering a wide interpeptide distance ranging from the bound to unbound states. The 300‐K free‐energy landscape was characterized by thermodynamically stable basins of antiparallel and parallel β‐sheet complexes and some other complex forms. Helices were frequently observed, when the two peptides contacted loosely or fluctuated freely without interpeptide contacts. We observed that V‐AUS converged to uniform distribution more effectively than conventional AUS sampling did. © 2015 Wiley Periodicals, Inc.     

1,2-Dibromoethyl-trichlorosilane (CH2BrCHBrSiCl3): conformational structure and vibrational properties by gas-phase electron diffraction, infrared and Raman spectroscopy, and ab initio molecular orbital and density functional theory calculations

The molecular structure and conformational properties of 1,2-dibromoethyl-trichlorosilane (CH2BrCHBrSiCl3) have been investigated using gas-phase electron diffraction (GED) data recorded at a temperature of 100 degrees C, together with ab initio molecular orbital (MO) and density functional theory (DFT) calculations, infrared (IR) and Raman spectroscopy in the liquid and solid phases, and normal coordinate analysis (NCA). The molecule exists in the gas- and liquid phases as a mixture of three conformers, gauche(-) [G(-)], with a refined torsion angle phi(BrCCBr)=-71(6) degrees, anti [A], with a torsion angle phi(BrCCBr) approximately -170 degrees , and gauche(+) [G(+)], with a torsion angle phi(BrCCBr) approximately +70 degrees . The second torsion angle of importance, the rotation about the CSi bond, has been refined to a value of +175(13) degrees . Torsion angles were only refined for the more abundant G(-) conformer. In the solid phase, only the G(-) conformer was observed. The temperature-dependent Raman spectra have provided an estimate of the relative conformational entropies, DeltaS. The obtained composition from GED refinements was (%) G(-)/A/G(+)=64(27)/23(13)/13(18) (values with estimated 2sigma uncertainties), giving a conformational stability order in agreement with both the Raman enthalpy measurements and the ab initio MO and DFT calculations using the 6-311G(d) basis set and scaled zero-point energies. Relevant structural parameter values obtained from the GED refinements (with the ab initio HF values used as constraints) were as follows (G(-) values with estimated 2sigma uncertainties): bond lengths (r(g)):r(C-C)=1.501(18)A, r(SiC)=1.865(15)A, r(CBr)=1.965(8)A (average), r(SiCl)=2.028(3)A (average). Bond angles ( anglealpha):angleCCSi=114.1(33) degrees , angleC1C2Br=114.0(21) degrees , angleCSiCl=109.6(7) degrees (average). Experimental IR/Raman and obtained vibrational wavenumbers based on both the unscaled, fixed-scaled as well as the scale-refined quantum-mechanical force fields [HF/6-311G(d)] are presented. The results are discussed and compared with some similar molecules from the literature.