Molecular dynamics simulations of the enzyme Catechol-O-Methyltransferase: methodological issues1

Results from extensive 70 ns all-atom molecular dynamics simulations of catechol-O-methyltransferase (COMT) enzyme are reported. The simulations were performed with explicit TIP3P water and Mg^2?+ ions. Four different crystal structures of COMT, with and without different ligands, were used. These simulations are among the most extensive of their kind and as such served as a stability test for such simulations. On the methodological side we found that the initial energy minimization procedure may be a crucial step: particular hydrogen bonds may break, and this can initiate an irreversible loss of protein structure that becomes observable in longer time scales of the order of tens of nanoseconds. This has important implications for both molecular dynamics and quantum mechanics–molecular mechanics simulations.     

Molecular dynamics simulations of the enzyme Cu, Zn superoxide dismutase

The enzyme Cu, Zn superoxide dismutase (Cu,Zn-SOD) is a ubiquitous oxireductase, which is responsible for the cellular defense against oxidative stress caused by the high toxicity of the superoxide radical, and has been also linked to some cases of familiar amyotrophic lateral sclerosis. In the present study a set of molecular mechanics parameters for the active site of Cu,Zn-SOD has been derived. Afterward, an extensive molecular dynamics simulation has been carried out in an aqueous environment. The obtained results shed a further light on the structural flexibility of the backbone, where the active site is nested, and the solvation shell occupancy. The relatively small backbone deviation, shown by a root-mean-square deviation below 1.0 A, confirms the accuracy of the parameters. The solvent shell analysis has shown that the first solvation shell is located at about 5 A from the copper ion, generating an empty cavity with enough space to accommodate the superoxide radical. The low residence time means that a high permutation rate of water molecules in both solvation shells is consistent with the efficiency of this catalytic mechanism. Hybrid studies using ONIOM methodologies can now be done to evaluate the mechanistic implications of the explicit inclusion of the whole system.     

Molecular dynamics simulations of halide glasses

Summary Halide glasses have been extensively studied in recent years because of their potential application as infrared transmitting fibre optic materials. They are believed to be more ionic than glasses based on silica and should therefore be more amenable to molecular dynamics simulation using simple two-body potentials. Here the main features of structural models derived using such techniques are described and compared with available structural data. Possible future applications of this approach are outlined.     

Molecular dynamics simulations of peptide-surface interactions

Proteins, which are bioactive molecules, adsorb on implants placed in the body through complex and poorly understood mechanisms and directly influence biocompatibility. Molecular dynamics modeling using empirical force fields provides one of the most direct methods of theoretically analyzing the behavior of complex molecular systems and is well-suited for the simulation of protein adsorption behavior. To accurately simulate protein adsorption behavior, a force field must correctly represent the thermodynamic driving forces that govern peptide residue-surface interactions. However, since existing force fields were developed without specific consideration of protein-surface interactions, they may not accurately represent this type of molecular behavior. To address this concern, we developed a host-guest peptide adsorption model in the form of a G(4)-X-G(4) peptide (G is glycine, X is a variable residue) to enable determination of the contributions to adsorption free energy of different X residues when adsorbed to functionalized Au-alkanethiol self-assembled monolayers (SAMs). We have previously reported experimental results using surface plasmon resonance (SPR) spectroscopy to measure the free energy of peptide adsorption for this peptide model with X = G and K (lysine) on OH and COOH functionalized SAMs. The objectives of the present research were the development and assessment of methods to calculate adsorption free energy using molecular dynamics simulations with the GROMACS force field for these same peptide adsorption systems, with an oligoethylene oxide (OEG) functionalized SAM surface also being considered. By comparing simulation results to the experimental results, the accuracy of the selected force field to represent the behavior of these molecular systems can be evaluated. From our simulations, the G(4)-G-G(4) and G(4)-K-G(4) peptides showed minimal to no adsorption to the OH SAM surfaces and the G(4)-K-G(4) showed strong adsorption to the COOH SAM surface, which is in agreement with our SPR experiments. Contrary to our experimental results, however, the simulations predicted a relatively strong adsorption of G(4)-G-G(4) peptide to the COOH SAM surface. In addition, both peptides were unexpectedly predicted to adsorb to the OEG surface. These findings demonstrate the need for GROMACS force field parameters to be rebalanced for the simulation of peptide adsorption behavior on SAM surfaces. The developed methods provide a direct means of assessing, modifying, and validating force field performance for the simulation of peptide and protein adsorption to surfaces, without which little confidence can be placed in the simulation results that are generated with these types of systems.     

Molecular dynamics simulations of lipid bilayers

Computer simulation methods are becoming increasingly widespread as tools for studying the structure and dynamics of lipid bilayer membranes. The length scale and time scale accessible to atomic-level molecular dynamics simulations are rapidly increasing, providing insight into the relatively slow motions of molecular reorientation and translation and demonstrating that effects due to the finite size of the simulation cell can influence simulation results. Additionally, significant advances have been made in the complexity of membrane systems studied, including bilayers with cholesterol, small solute molecules, and lipid-protein and lipid-DNA complexes. Especially promising is the progress that continues to be made in the comparison of simulation results with experiment, both to validate the simulation algorithms and to aid in the interpretation of existing experimental data.     

Molecular dynamics simulations of polyelectrolyte adsorption

We have performed molecular dynamics simulations of polyelectrolyte adsorption at oppositely charged surfaces from dilute polyelectrolyte solutions. In our simulations, polyelectrolytes were modeled by chains of charged Lennard-Jones particles with explicit counterions. We have studied the effects of the surface charge density, surface charge distribution, solvent quality for the polymer backbone, strength of the short-range interactions between polymers and substrates on the polymer surface coverage, and the thickness of the adsorbed layer. The polymer surface coverage monotonically increases with increasing surface charge density for almost all studied systems except for the system of hydrophilic polyelectrolytes adsorbing at hydrophilic surfaces. In this case the polymer surface coverage saturates at high surface charge densities. This is due to additional monomer-monomer repulsion between adsorbed polymer chains, which becomes important in dense polymeric layers. These interactions also preclude surface overcharging by hydrophilic polyelectrolytes at high surface charge densities. The thickness of the adsorbed layer shows monotonic dependence on the surface charge density for the systems of hydrophobic polyelectrolytes for both hydrophobic and hydrophilic surfaces. Thickness is a decreasing function of the surface charge density in the case of hydrophilic surfaces while it increases with the surface charge density for hydrophobic substrates. Qualitatively different behavior is observed for the thickness of the adsorbed layer of hydrophilic polyelectrolytes at hydrophilic surfaces. In this case, thickness first decreases with increasing surface charge density, then it begins to increase.     

Molecular dynamics simulations of biotin carboxylase

Biotin carboxylase catalyzes the ATP-dependent carboxylation of biotin and is one component of the multienzyme complex acetyl-CoA carboxylase that catalyzes the first committed step in fatty acid synthesis in all organisms. Biotin carboxylase from Escherichia coli, whose crystal structures with and without ATP bound have been determined, has served as a model system for this component of the acetyl-CoA carboxylase complex. The two crystal structures revealed a large conformational change of one domain relative to the other domains when ATP is bound. Unfortunately, the crystal structure with ATP bound was obtained with an inactive site-directed mutant of the enzyme. As a consequence the structure with ATP bound lacked key structural information such as for the Mg2+ ions and contained altered conformations of key active-site residues. Therefore, nanosecond molecular dynamics studies of the wild-type biotin carboxylase were undertaken to supplant and amend the results of the crystal structures. Specifically, the protein-metal interactions of the two catalytically critical Mg2+ ions bound in the active site are presented along with a reevaluation of the conformations of active-site residues bound to ATP. In addition, the regions of the polypeptide chain that serve as hinges for the large conformational change were identified. The results of the hinge analysis complemented a covariance analysis that identified the individual structural elements of biotin carboxylase that change their conformation in response to ATP binding.     

Molecular dynamics simulations of cardiolipin bilayers

Cardiolipin is a key lipid component in the inner mitochondrial membrane, where the lipid is involved in energy production, cristae structure, and mechanisms in the apoptotic pathway. In this article we used molecular dynamics computer simulations to investigate cardiolipin and its effect on the structure of lipid bilayers. Three cardiolipin/POPC bilayers with different lipid compositions were simulated: 100, 9.2, and 0% cardiolipin. We found strong association of sodium counterions to the carbonyl groups of both lipid types, leaving in the case of 9.2% cardiolipin virtually no ions in the aqueous compartment. Although binding occurred primarily at the carbonyl position, there was a preference to bind to the carbonyl groups of cardiolipin. Ion binding and the small headgroup of cardiolipin gave a strong ordering of the hydrocarbon chains. We found significant effects in the water dipole orientation and water dipole potential which can compensate for the electrostatic repulsion that otherwise should force charged lipids apart. Several parameters relevant for the molecular structure of cardiolipin were calculated and compared with results from analyses of coarse-grained simulations and available X-ray structural data.     

Molecular dynamics simulations and free energy calculations on the enzyme 4‐hydroxyphenylpyruvate dioxygenase

4‐Hydroxyphenylpyruvate dioxygenase is a relevant target in both pharmaceutical and agricultural research. We report on molecular dynamics simulations and free energy calculations on this enzyme, in complex with 12 inhibitors for which experimental affinities were determined. We applied the thermodynamic integration approach and the more efficient one‐step perturbation. Even though simulations seem well converged and both methods show excellent agreement between them, the correlation with the experimental values remains poor. We investigate the effect of slight modifications on the charge distribution of these highly conjugated systems and find that accurate models can be obtained when using improved force field parameters. This study gives insight into the applicability of free energy methods and current limitations in force field parameterization. © 2011 Wiley Periodicals, Inc. J Comput Chem 2011     

Molecular dynamics simulations of lithium alkyl carbonates

A quantum chemistry-based many-body polarizable force field has been developed for two model solid-electrolyte interphase (SEI) components: dilithium ethylene dicarbonate (Li(2)EDC) and lithium methyl carbonate (LiMC). Molecular dynamics (MD) simulations of amorphous Li(2)EDC and LiMC were performed at temperatures from 393 to 600 K. Simulations reveal that Li(+) is coordinated by approximately 4.6 oxygen atoms from -COO(-) groups coming from different alkyl carbonate molecules. Charge transport in Li(2)EDC was found to be almost entirely due to Li(+). The temperature dependence of the ionic conductivity of the SEI model compounds Li(2)EDC and LiMC was found to be significantly stronger than that of liquid electrolytes (e.g., ethylene carbonate + LiTFSI), yielding extrapolated conductivities of the Li(2)EDC on order of 10(-10) S/cm at -30 degrees C.     

Molecular dynamics simulations of peptide carboxylate hydration

Aqueous solvation of carboxylate groups, as present in the glycine zwitterion and the dipeptide aspartylalanine, is studied employing a force-field that includes distributed multipole electrostatics and induction contributions (Amoebapro: P. Ren and J. W. Ponder, J. Comput. Chem., 2002, 23, 1497; P. Ren and J. W. Ponder, J. Phys. Chem. B, 2003, 107, 5933; J. W. Ponder and D. A. Case, Adv. Protein Chem., 2003, 66, 27). Radial and orientation distribution functions, as well as hydration numbers, are calculated and compared with existing simulation data derived from Car-Parrinello molecular dynamics (CPMD), and also distributed-charge force-fields. Connections are also made with experimental data for solvation of carboxylates in water. Our findings show that Amoebapro yields carboxylate solvation properties in very good agreement with CPMD results, significantly closer agreement than can be obtained from traditional force-fields. We also demonstrate that the influence of solvation on the conformation of the dipeptide is markedly different using Amoebapro compared with the other force-fields.     

Molecular dynamics simulations of local field factors

In the present work, the authors evaluate a scheme based on molecular dynamics to derive local field factors. These are given without any assumption of a cavity by fitting the Langevin functions to the order parameters obtained from the molecular dynamics simulations. The local field factors so obtained, with the detailed chromophore-solvent interactions and solvent structures taken into account, are much smaller than those calculated from the conventional Onsager and Lorentz models. A numerical demonstration is given for two typical organic chromophore molecules, p-nitroaniline and p-nitro-N,N-dimethylaniline dissolved in chloroform.     

Molecular dynamics simulations of supramolecular polymer rheology

Using equilibrium and nonequilibrium molecular dynamics simulations, we studied the equilibrium and rheological properties of dilute and semidilute solutions of head-to-tail associating polymers. In our simulation model, a spontaneous complementary reversible association between the donor and the acceptor groups at the ends of oligomers was achieved by introducing a combination of truncated pseudo-Coulombic attractive potential and Lennard Jones repulsive potential between donor, acceptor, and neighboring groups. We have calculated the equilibrium properties of supramolecular polymers, such as the ring/chain equilibrium, average molecular weight, and molecular weight distribution of self-assembled chains and rings, which all agree well with previous analytical and computer modeling results. We have investigated shear thinning of solutions of 8- and 20-bead associating oligomers with different association energies at different temperatures and oligomer volume fractions. All reduced viscosity data for a given oligomer length can be collapsed into one master curve, exhibiting two power-law regions of shear-thinning behavior with an exponent of -0.55 at intermediate ranges of the reduced shear rate β and -0.8 (or -0.9) at larger shear rates. The equilibrium viscosity of supramolecular solutions with different oligomer lengths and associating energies is found to obey a power-law scaling dependence on oligomer volume fraction with an exponent of 1.5, in agreement with the experimental observations for several dilute or semidilute solutions of supramolecular polymers. This implies that dilute and semidilute supramolecular polymer solutions exhibit high polydispersity but may not be sufficiently entangled to follow the reptation mechanism of relaxation.     

Molecular dynamics simulations of the oxidation of aluminum nanoparticles

The oxidation of aluminum nanoparticles is studied with classical molecular dynamics and the Streitz-Mintmire (Streitz, F. H.; Mintmire, J. W. Phys. Rev. B 1994, 50, 11996) electrostatic plus (ES+) potential that allows for the variation of electrostatic charge on all atoms in the simulation. The structure and charge distributions of bulk crystalline alpha-Al(2)O(3), a surface slab of alpha-Al(2)O(3) with an exposed (0001) basal plane, and an isolated Al(2)O(3) nanoparticle are studied. Constant NVT simulations of the oxidation of aluminum nanoparticles are also performed with different oxygen exposures. The calculations simulate a thermostated one-time exposure of an aluminum nanoparticle to different numbers of surface oxygen atoms. In the first set of oxidation studies, the overall approximate ratios of Al to O in the nanoparticle are 1:1 and 2:1. The nanoparticles are annealed to 3000 K and are then cooled to 500, 1000, or 1500 K. The atomic kinetic energy is scaled during the simulation to maintain the desired temperature. The structure and charge distributions in the oxidized nanoparticles differ from each other and from those of the bulk Al(2)O(3) phases. In the Al(1)O(1) simulation, an oxide shell forms that stabilizes the shape of the particle, and thus the original structure of the nanoparticle is approximately retained. In the case of Al(1)O(0.5), there is insufficient oxygen to form a complete oxide shell, and the oxidation results in particles of irregular shapes and rough surfaces. The particle surface is rough, and the nanoparticle is deformed.     

Molecular dynamics simulations of the melting of 1,3,3-trinitroazetidine

Physical properties of condensed-phase 1,3,3-trinitroazetidine (TNAZ) have been computed with molecular dynamics (MD) and a nonreactive, fully flexible force field formulated by combining the intramolecular interactions obtained from the Generalized AMBER Force Field and the rigid-molecule force field developed by Sorescu-Rice-Thompson [J. Phys. Chem. B 1997, 101, 798] (AMBER-SRT). The results are compared with MD calculations, using the AMBER force field. The predicted densities of crystalline TNAZ from both force fields are about 10% lower than the experimental value. The calculated thermodynamic melting point at 1 atm from the AMBER-SRT force field is 390 K, in good agreement with the measured value of 374 K, while the AMBER force field predicts a thermodynamic melting point of 462 K. The lattice parameters and the molecular and crystal structures calculated with the AMBER-SRT force field are in excellent agreement with experiment. Simulations with the AMBER-SRT force field were also used to generate the isotherm of TNAZ up to 4 GPa and the bulk modulus and its pressure derivative.     

Molecular dynamics simulations of the melting of aluminum nanoparticles

Molecular dynamics simulations are performed to determine the melting points of aluminum nanoparticles of 55-1000 atoms with the Streitz-Mintmire [Phys. Rev. B 1994, 50, 11996] variable-charge electrostatic plus potential. The melting of the nanoparticles is characterized by studying the temperature dependence of the potential energy and Lindemann index. Nanoparticles with less than 850 atoms show bistability between the solid and liquid phases over temperature ranges below the point of complete melting. The potential energy of a nanoparticle in the bistable region alternates between values corresponding to the solid and liquid phases. This bistability is characteristic of dynamic coexistence melting. At higher temperatures, only the liquid state is stable. Nanoparticles with more than 850 atoms undergo a sharp solid-liquid-phase transition characteristic of the bulk solid phase. The variation of the melting point with the effective nanoparticle radius is also determined.     

Molecular dynamics simulations of charged dendrimers: low-to-intermediate half-generation PAMAMs

We have performed atomistic molecular dynamics simulations of PAMAM dendrimers of generations 0.5, 1.5, 2.5, 3.5, and 4.5. The simulated systems comprise the charged dendrimer and its counterions embedded in a dielectric continuum (i.e., without explicit solvent). Structural properties of these dendrimers, like the radius of gyration, the principal moments of inertia, and the segment density profiles, were evaluated from the simulations. The average radius of gyration obtained for the intermediate half-generations 2.5, 3.5, and 4.5 follows the same scaling law that was previously inferred from simulations of full-generation PAMAMs, Rg approximately M1/3, and is characteristic of space-filling objects. The low half-generations 0.5 and 1.5 deviate, however, to greater Rg values. The shape of the smaller dendrimers is approximately that of a prolate ellipsoid, which becomes more spherical for higher generations. The segment density profiles show features identical to those obtained in other simulations of flexible-chain dendrimers, like dendron-backfolding. Two slightly different configurations, in terms of size and shape, were identified for generation 2.5. The radial distributions of counterions extracted from the simulations compare well with the solutions of Poisson-Boltzmann cell model, and the dendrimer's effective charge was estimated using the Bjerrum criterion. The influence of electrostatic interactions in the dendrimer's conformation due to repulsion between the charged end-groups and its relation to counterion effects is discussed for the several generations simulated. The form factors calculated from the simulations are compared with the model of a homogeneous ellipsoid of revolution. The overall results are in agreement with the previously established morphological transition of PAMAM dendrimers toward a more spherical and compact conformation above generations 3 or 4.     

Molecular dynamics simulations of aqueous solutions of ethanolamines

We report on molecular dynamics simulations performed at constant temperature and pressure to study ethanolamines as pure components and in aqueous solutions. A new geometric integration algorithm that preserves the correct phase space volume is employed to study molecules having up to three ethanol chains. The most stable geometry, rotational barriers, and atomic charges were obtained by ab initio calculations in the gas phase. The calculated dipole moments agree well with available experimental data. The most stable conformation, due to intramolecular hydrogen bonding interactions, has a ringlike structure in one of the ethanol chains, leading to high molecular stability. All molecular dynamics simulations were performed in the liquid phase. The interaction parameters are the same for the atoms in the ethanol chains, reducing the number of variables in the potential model. Intermolecular hydrogen bonding is also analyzed, and it is shown that water associates at low water mole fractions. The force field reproduced (within 1%) the experimental liquid densities at different temperatures of pure components and aqueous solutions at 313 K. The excess and partial molar volumes are analyzed as a function of ethanolamine concentration.     

Molecular dynamics simulations of surface-initiated melting of nitromethane

The melting of nitromethane initiated at solid-vacuum interfaces has been investigated using molecular dynamics nvt simulations with a realistic force field [D. C. Sorescu et al., J. Phys. Chem. B 104, 8406 (2000)]. The calculated melting point (251+/-5 K) is in good agreement with experiment (244.73 K) and values obtained previously (approximately 255.5 and 266.5+/-8 K) using other simulation methods [P. M. Agrawal et al., J. Chem. Phys. 119, 9617 (2003)]. Analyses of the molecular orientations and diffusion during the simulations as functions of the distance from the exposed surfaces show that the melting is a direct crystal-to-liquid transition, in which the molecules first gain rotational freedom, then mobility. There is a slight dependence of the melting temperature on the exposed crystallographic face.