Molecular Dynamics as a Mathematical Mapping. II. Partial Derivatives in the Microcanonical Ensemble

Abstract

Performing molecular dynamics in a fully continuous and differentiable framework can be viewed as a deterministic mathematical mapping between, on one side, the force field parameters that describe the potential energy interactions and input macroscopic conditions, and, on the other, the calculated corresponding macroscopic properties of the bulk molecular system.

Within this framework, it is possible to apply standard methods of variational calculus for the computation of the partial derivatives of the molecular dynamics mapping based on the integration of either the adjoint equations or the sensitivity equations of the classical Newtonian equations of motion. We present procedures for these computations in the standard microcanonical (N, V, E) ensemble, and compare the computational efficiency of the two approaches. The general formulations developed are applied to the specific example of bulk ethane fluid.

With these procedures in place, it is now possible to compute the partial derivatives of any property determined by molecular dynamics with respect to any input property and any potential parameter. Moreover, these derivatives are computed to essentially the same level of numerical accuracy as the output properties themselves.     

A Theoretical Study of the Aqueous Hydration of Canonical B d(CGCGAATTCGCG): Monte Carlo Simulation and Comparison with Crystallographic Ordered Water Sites

Abstract

Monte Carlo computer simulation is described for the dodecamer d(CGCGAATTCGCG) together with 1777 water molecules at an environmental density of 1 gm/cc in a cubic cell under periodic boundary conditions. Water-water interactions were treated using the TIP4P potential and the solute water interactions by TIP4P spliced with the non-bonded interactions from the AMBER 3.0 force field. The simulation was subjected to proximity analysis to obtain solute coordination numbers and pair interaction energies for each solute atom. Hydration density distributions partitioned into contributions from the major groove side, the minor groove side and the sugar-phosphate backbone were examined, and the probabilities of occurence for one- and two-water bridges in the simulation were enumerated. The results were compared with observations of crystallographic ordered water sites from x-ray diffraction studies on the native dodecamer by Dickerson and coworkers.     

Combined use of periodic reaction field and coarse-grained molecular dynamics simulations. I. phospholipid monolayer systems

Abstract

A periodic reaction field based on a linear-combination-based isotropic periodic sum (LIPS) method was applied for coarse-grained molecular dynamics simulations of zwitterionic lipid systems. In phospholipid monolayer systems with various number of lipid molecules, the density profile, lipid orientation and surface tension were mainly calculated using the periodic reaction field and Ewald sum. The results from the periodic reaction field were almost equal to that from the Ewald sum. It is concluded that the periodic reaction field method has a great possibility to provide a high accuracy in determining coarse-grained zwitterionic lipid systems.     

Molecular Dynamics as a Mathematical Mapping. III. Efficient Evaluation of the Differentiable Force Functions and Their Derivatives

Abstract

The fully continuous and differentiable framework for performing molecular dynamics calculations introduced in parts I and II of this paper [1,2] requires the evaluation of rather complex force functions and their spatial partial derivatives. This paper presents an efficient interpolation scheme for the evaluation of these quantities over a finite spatial domain.

The modified force function is approximated by a linear combination of Hermite cubic basis functions such that both the interpolant of the force and its spatial derivatives are continuous across the grid boundaries. In order to achieve better accuracy for a given grid size, a nonuniform rectilinear grid is constructed via iterative refinement procedure. The latter guarantees the accuracy of the force computed by interpolation within any specified tolerance > ε O.

For many potential functions of practical interest, it is possible for polynomial interpolants to be constructed for parts of the force functions which are independent of the potential parameters and system density (the so-called “separable force functions”). In such cases, a single interpolation grid which is applicable for a wide range of potential parameters and system densities can be constructed a priori.     

The Calculation of 3D Solubility Parameters Using Molecular Models

Abstract

In this paper we describe the use of molecular mechanics models to examine detailed intermolecular interactions within the liquid state of a common nonionic surfactant system, nonyl phenol ethoxylate (NPE). Using constant energy molecular dynamics simulations we have studied the relative strengths of dispersive interactions versus polar interactions and have estimated three dimensional solubility parameters for NPE systems as a function of temperature and ethylene oxide content. The predictions at 300 K are in good agreement with three dimensional solubility parameters predicted using group contribution tables. Models of the amorphous liquid state were represented by single molecular structures of NPE in a periodic cell. The solubility parameters predicted with these models were in good agreement with those values derived from models having eight NPE molecules packed into a cell with the exception of the electrostatic interactions, which are the most sensitive to system size effects.     

New nonbonded interactions calculation strategy for rectangular systems. I. Preliminary molecular dynamics study of solvated Na(+) ion

A new molecular nonbonded interactions treatment strategy is proposed in the context of rectangular periodic boundary conditions simulations. Several molecular dynamics simulations are performed on a sodium ion in aqueous solution. Box sizes are modified from a cubic to a rectangular shape. The results are compared with those found using a classical spherical cutoff. This new method yields ion-oxygen radial distribution functions in good agreement with experimental results, thus showing its reliability. Severe perturbations in the structural orientation of water molecules in the first shell with the increase of the box length are observed under the classical cutoff method. However, these distorting effects are reduced with the present nonbonded interactions treatment.     

Charge grouping approaches to calculation of electrostatic forces in molecular dynamics of macromolecules.

Calculation of long-range electrostatic interactions is the most time-consuming step in theoretical simulation of the structure and dynamics of macromolecules. In practice very short cutoff distances are used, which may distort the behavior of the model system. We describe two accurate approaches to calculation of electrostatic forces based on hierarchical grouping of charges into cubes. The first is similar to the O(NlogN) algorithm developed by Barnes, J. and Hut, P., Nature (London) 324, 446-449 (1986), for simulation of a gravitational motion of N bodies. The second approach we formulate for a system with periodic boundary conditions in the nearest image approximation. The calculation of electrostatic interactions and a charge grouping procedure are faster than O(N2). The average inaccuracy in the force introduced by the grouping does not exceed 1%. We describe a small modification of the same approach which makes it suitable for long strongly charged polymers as well. This accurate approach to calculation of electrostatic interactions is illustrated with an example of the dynamics of ions near DNA. Quick equilibration of the ionic distribution is observed during molecular dynamics simulation if electrostatic forces are properly calculated, while the behavior and distribution of ions are less realistic when the conventional cutoff distances are used.     

Conformational Flexibility of Two RNA Trimers Explored by Computational Tools and Database Search

Abstract

Two RNA sequences, AAA and AUG, were studied by the conformational search program CICADA and by molecular dynamics (MD) in the framework of the AMBER force field, and also via thorough PDB database search. CICADA was used to provide detailed information about conformers and conformational interconversions on the energy surfaces of the above molecules. Several conformational families were found for both sequences. Analysis of the results shows differences, especially between the energy of the single families, and also in flexibility and concerted conformational movement. Therefore, several MD trajectories (altogether 16 ns) were run to obtain more details about both the stability of conformers belonging to different conformational families and about the dynamics of the two systems. Results show that the trajectories strongly depend on the starting structure. When the MD start from the global minimum found by CICADA, they provide a stable run, while MD starting from another conformational family generates a trajectory where several different conformational families are visited. The results obtained by theoretical methods are compared with the thorough database search data. It is concluded that all except for the highest energy conformational families found in theoretical result also appear in experimental data.

Registry numbers:

adenylyl-(3′ →5′)-adenylyl-(3′ →5′)-adenosine [917-44-2]

adenylyl-(3′ →5′)-uridylyl-(3′ →5′)-guanosine [3494-35-7]     

Vorinostat,a possible alternative to metronidazole for the treatment of amebiasis caused by Entamoeba histolytica

Abbreviations SAHA suberoylanilide hydroxamic acid

EhHDAC Histone Deacetylase from Entamoeba histolytica

Rg Radius of gyration

RMSD root-mean-square deviation

RMSF root-mean-square fluctuation

MDS molecular dynamics simulation

VMD Visual Molecular Dynamics

NAMD Nanoscale Molecular Dynamics

PBC periodic boundary conditions

PME Particle Mesh Ewald

3D three-dimensional

Cα alpha carbon

FDA Food and Drug Administration

ns nanoseconds

GPU CUDA Graphics Processing Unit Compute Unified Device Architecture

Communicated by Ramaswamy H. Sarma     

Mean field stochastic boundary molecular dynamics simulation of a phospholipid in a membrane

Computer simulations of phospholipid membranes have been carried out by using a combined approach of molecular and stochastic dynamics and a mean field based on the Marcelja model. First, the single-chain mean field simulations of Pastor et al. [(1988) J. Chem. Phys. 89, 1112-1127] were extended to a complete dipalmitoylphosphatidylcholine molecule; a 102-ns Langevin dynamics simulation is presented and compared with experiment. Subsequently, a hexagonally packed seven-lipid array was simulated with Langevin dynamics and a mean field at the boundary and with molecular dynamics (and no mean field) in the center. This hybrid method, mean field stochastic boundary molecular dynamics, reduces bias introduced by the mean field and eliminates the need for periodic boundary conditions. As a result, simulations extending to tens of nanoseconds may be carried out by using a relatively small number of molecules to model the membrane environment. Preliminary results of a 20-ns simulation are reported here. A wide range of motions, including overall reorientation with a nanosecond decay time, is observed in both simulations, and good agreement with NMR, IR, and neutron diffraction data is found.     

An Efficient,Box Shape Independent Non-Bonded Force and Virial Algorithm for Molecular Dynamics

Abstract

A notation is introduced and used to transform a conventional specification of the non-bonded force and virial algorithm in the case of periodic boundary conditions into an alternative specification. The implementation of the transformed specification is simpler and typically a factor of 1.5 faster than a conventional implementation. Moreover, it is generic with respect to the shape of the simulated system, i.e. the same routines can be used to handle triclinic boxes, truncated octahedron boxes etc. An implementation of this method is presented, and the speed achieved on various machines is given. Essence of the new method is that the number of calculations of image particle positions is strongly reduced during non-bonded force calculations.     

Towards molecular dynamics simulation of large proteins with a hydration shell at constant pressure.

Molecular dynamics simulation of a large protein in explicit water with periodic boundary conditions is extremely demanding in terms of computation time. Consequently, we have sought approximations of the solvent environment that model its important features. Here, we describe our SAPHYR (Shell Approximation for Protein HYdRation) model in which the protein is surrounded by a shell of water molecules maintained at constant pressure. In addition to the usual pairwise interatomic interactions, these water molecules are subjected to forces approximating van der Waals and dipole-dipole interactions with the implicit surrounding bulk solvent. The SAPHYR model is tested for a system of one argon atom in water and for the protein ubiquitin, and then applied to cytochrome P450cam, a protein with over 400 residues. The results demonstrate that structural and dynamic properties of the simulated systems are improved by use of the SAPHYR model, and that this model provides a significant computational saving over simulations with periodic boundary conditions.     

Molecular Modelling of The Mechanism of Action of Organic Clay-Swelling Inhibitors

Abstract

It is well known that the sodium smectite class of clays swells macroscopically in contact with water, whereas under normal conditions the potassium form does not. In recent work using molecular simulation methods, we have provided a quantitative explanation both for the swelling behaviour of sodium smectite clays and the lack of swelling of potassium smectites [1]. In the present paper, we apply similar modelling methods to study the mechanism of inhibition of clay-swelling by a range of organic molecules.

Experimentally, it is known that polyalkylene glycols (polyethers) of intermediate to high relative molecular mass are effective inhibitors of smectite clay swelling. We use a range of atomistic simulation techniques, including Monte Carlo and molecular dynamics, to investigate the interactions between a selection of these compounds, water, and a model smectite clay mineral. These interactions occur by means of organised intercalation of water and organic molecules within the galleries between individual clay layers.

The atomic interaction potentials deployed in this work are not as highly optimised as those used in our clay-cation-water work [1]. Nevertheless, our simulations yield trends and results that are in qualitative and sometimes semi-quantitative agreement with experimental findings on similiar (but not identical) systems. The internal energy of adsorption of simple polyethers per unit mass on the model clay is not significantly different from that for water adsorption; our Monte Carlo studies indicate that entropy is the driving force for the sorption of the simpler organic molecules inside the clay layers: a single long chain polyethylene glycol can displace a large number of water molecules, each of whose translational entropy is greatly enhanced when outside the clay. Hydrophobically modified polyalkylene glycols also enjoy significant van der Waals interactions within the layers which they form within the clay galleries.

In conjunction with experimental studies, our work furnishes valuable insights into the relative effectiveness of the compounds considered and reveals the generic features that high performance clay-swelling inhibitors should possess. For optimal inhibitory activity, these compounds should be reasonably long chain linear organic molecules with localised hydrophobic and hydrophilic regions along the chain. On intercalation of these molecules within the clay layers, the hydrophobic regions provide an effective seal against ingress of water, while the hydrophilic ones enhance the binding of the sodium cations to the clay surface, preventing their hydration and the ensuing clay swelling.     

On the Approximation of Solvent Effects on the Conformation and Dynamics of Cyclosporin A by Stochastic Dynamics Simulation Techniques

Abstract

The molecular simulation technique of stochastic dynamics (SD) is tested by application to the immunosuppressive drug cyclosporin A (CPA). Two stochastic dynamics simulations are performed, one (SD_CCl4 ) with atomic friction coefficients proportional to the viscosity of the nonpolar solvent CCl₄, and one (SD_H2O) with atomic friction coefficients corresponding to an aqueous solution. The atomic friction coefficients are also taken proportional to an approximate expression for the atomic accessible surface area. The properties of both stochastic dynamics simulations are compared to those of two full molecular dynamics (MD) simulations of cyclosporin A, one in a box with 591 CCl₄ molecules, and one in a box with 632 H₂O molecules.

The properties of cyclosporin A as found in the molecular dynamics simulation in CCl₄ are well reproduced by the SD_CCl4 simulation. This indicates that the neglect of a mean force reresenting the average solvent effects on the solute is justified in the case of nonpolar solvents. For polar solvents, like water, this mean force may not be neglected. The SD_H2O simulation of cyclosporin A clearly fails to reproduce the amount of hydrogen bonding found in the molecular dynamics stimulation of cyclosporin A in water.

A comparison with a molecular dynamics simulation of cyclosporin A in vacuo shows that both the SD_CCl4 and the SD_H2O simulation come closer to the properties of the molecular dynamics simulations in CCl₄ and in H₂O than a molecular dynamics simulation in vacuo.     

Long-range electrostatic effects on peptide folding.

The reversible folding/unfolding of a short peptide in solution is studied by molecular dynamics simulations. The effects of long-range electrostatic interactions are examined and found to be important both for the equilibrium between folded and unfolded states and the dynamics of the folding process. The neglect of long-range electrostatics leads to an increased population of unfolded states and increased structural fluctuations. When such interactions are taken into account, the peptide unfolds and folds to the experimentally determined structure several times during a 25 ns trajectory, with approximately equal populations of folded and unfolded states in the neighborhood of its proposed melting temperature. The effect of using spherical boundary conditions rather than periodic ones does not appear to have any major effect on the folding dynamics.     

Investigations of Nematic-Isotropic Transition by Means of Constant Pressure Molecular Dynamics Simulations

Abstract

Constant pressure molecular dynamics simulations, which secure the system to be under hydrostatic pressure, are used to simulate the behavior of liquid crystals consisting of anisotropic molecules with both translational and orientational freedom. In order to investigate to what extent can the properties known to real liquid crystalline phases be explained by the anisotropy of the shape of the molecules alone, the molecular dynamic (MD) simulation uses purely repulsive short-range pair potentials representing soft spherocylinders. A clear change in the microscopic as well as the macroscopic physical properties are observed near the phase transition from nematic liquid crystal to isotropic liquid.     

Theoretical probes of conformational fluctuations in S-peptide and RNase A/3′–UMP enzyme product complex

The dynamic properties of the RNase A/3′–UMP enzyme/product complex and the S-peptide of RNase A have been investigated by molecular dynamics simulations using suitable generalization of ideas introduced to probe the energy landscape in structural glasses. We introduce two measures, namely, the kinetic energy fluctuation metric and the force metric, both of which are used to calculate the time needed for sampling the conformation space of the molecules. The calculation of the fluctuation metric requires a single trajectory whereas the force metric is computed using two independent trajectories. The vacuum MD simulations show that for both systems the time required for kinetic energy equipartitioning is surprisingly long even at high temperatures. We show that the force metric is a powerful means of probing the nature and relative importance of conformational substates which determine the dynamics at low temperatures. In particular the time dependence of the non-bonded force metric is used to demonstrate that at low temperatures the system is predominantly localized hi a single cluster of conformational substates. The force metric is used to show that relaxation of long range (in sequence space) interactions must be mediated by a sequence of local dihedral angle transitions. We also argue that the time needed for compact structure formation is intimately related to the time needed for the relaxation of the dihedral angle degrees of freedom. The tame for non-bonded interactions, which drive protein molecules to fold under appropriate conditions, to relax becomes extremely long as the temperature is lowered suggesting that the formation of maximally compact structure hi proteins must be a very slow process. © 1993 Wiley-Liss, Inc.     

Activated Complex Theory of Barite Scale Control Processes

A theoretical framework useful for the estimation of scale dissolution rate constants is introduced. The model consists of (a) a dissolution mechanism, (b) a quantum mechanical force field capable of describing the bulk and surface properties of barite, (c) a complete study of the structure of the dissolver molecule complex in solution, (d) a three dimensional periodic system useful for mapping the barite-dissolver interactions, including the localization of the activated complex and (e) a rate expression to estimate the dissolution rate constants from the properties of the activated complex. Our results show that molecular modeling, through a combination of molecular mechanics and high level quantum mechanical calculations, provide a new and insightful information about scale control processes. This revised version was published online in June 2006 with corrections to the Cover Date.     

REACH coarse-grained biomolecular simulation: transferability between different protein structural classes

Coarse graining of protein interactions provides a means of simulating large biological systems. The REACH (Realistic Extension Algorithm via Covariance Hessian) coarse-graining method, in which the force constants of a residue-scale elastic network model are calculated from the variance-covariance matrix obtained from atomistic molecular dynamics (MD) simulation, involves direct mapping between scales without the need for iterative optimization. Here, the transferability of the REACH force field is examined between protein molecules of different structural classes. As test cases, myoglobin (all α), plastocyanin (all β), and dihydrofolate reductase (α/β) are taken. The force constants derived are found to be closely similar in all three proteins. An MD version of REACH is presented, and low-temperature coarse-grained (CG) REACH MD simulations of the three proteins are compared with atomistic MD results. The mean-square fluctuations of the atomistic MD are well reproduced by the CGMD. Model functions for the CG interactions, derived by averaging over the three proteins, are also shown to produce fluctuations in good agreement with the atomistic MD. The results indicate that, similarly to the use of atomistic force fields, it is now possible to use a single, generic REACH force field for all protein studies, without having first to derive parameters from atomistic MD simulation for each individual system studied. The REACH method is thus likely to be a reliable way of determining spatiotemporal motion of a variety of proteins without the need for expensive computation of long atomistic MD simulations.