A smooth particle mesh Ewald method

MDSIMAID: Automatic parameter optimization in fast electrostatic algorithms

MDSIMAID is a recommender system that optimizes parallel Particle Mesh Ewald (PME) and both sequential and parallel multigrid (MG) summation fast electrostatic solvers. MDSIMAID optimizes the running time or parallel scalability of these methods within a given error tolerance. MDSIMAID performs a run time constrained search on the parameter space of each method starting from semiempirical performance models. Recommended parameters are presented to the user. MDSIMAID's optimization of MG leads to configurations that are up to 14 times faster or 17 times more accurate than published recommendations. Optimization of PME can improve its parallel scalability, making it run twice as fast in parallel in our tests. MDSIMAID and its Python source code are accessible through a Web portal located at http://mdsimaid.cse.nd.edu.

Molecular simulations for the conformational assessment of a porphyrin–fullerene dyad in different environments

Conformational space of a porphyrin-fullerene dyad with the donor and acceptor connected by a relatively flexible linker is studied by molecular dynamics simulations in both non-polar and polar solvents, as well as in vacuum. The most probable conformations obtained from the vacuum MD simulations were optimized with semi-empirical (SE) and density functional theory (DFT) methods and the extent of the structural changes is assessed. The computational results indicate the co-existence of different conformers in both polar and nonpolar solvents showing agreement with experimental results. The most probable vacuum conformations at 300 K are similar to the ones at 0 K, while the structures most often observed in the solvents show less compact conformations. Optimization with SE and DFT calculations leads to structures, which represent relatively well the folded conformations in solvent, which validates the electronic structure calculations relevant to describing photoinduced electron-transfer in H2P-O34-C60.

Molecular dynamics simulation of truncated bovine adrenodoxin

The 128 amino acid long soluble protein adrenodoxin (Adx) is a typical member of the ferredoxin protein family that are electron carrier proteins with an iron-sulfur cofactor. Adx carries electrons from adrenodoxin reductase (AdR) to cytochrome P450s. Its binding modes to these proteins were previously characterized by site-directed mutagenesis, by X-ray crystallography for the complex Adx:AdR, and by NMR. However, no clear evidence has been provided for the driving force that promotes Adx detachment from AdR upon reduction. Here, we characterized the conformational dynamics of unbound Adx in the oxidized and reduced forms using 2-20 ns long molecular dynamics simulations. The most noticeable difference between both forms is the enhanced flexibility of the loop (47-51) surrounding the iron-sulfur cluster in the reduced form. Together with several structural displacements at the binding interface, this increased flexibility may be the key factor promoting unbinding of reduced Adx from AdR. This points to an intrinsic property of reduced Adx that drives dissociation.

GROMACS: Fast, flexible, and free

This article describes the software suite GROMACS (Groningen MAchine for Chemical Simulation) that was developed at the University of Groningen, The Netherlands, in the early 1990s. The software, written in ANSI C, originates from a parallel hardware project, and is well suited for parallelization on processor clusters. By careful optimization of neighbor searching and of inner loop performance, GROMACS is a very fast program for molecular dynamics simulation. It does not have a force field of its own, but is compatible with GROMOS, OPLS, AMBER, and ENCAD force fields. In addition, it can handle polarizable shell models and flexible constraints. The program is versatile, as force routines can be added by the user, tabulated functions can be specified, and analyses can be easily customized. Nonequilibrium dynamics and free energy determinations are incorporated. Interfaces with popular quantum-chemical packages (MOPAC, GAMES-UK, GAUSSIAN) are provided to perform mixed MM/QM simulations. The package includes about 100 utility and analysis programs. GROMACS is in the public domain and distributed (with source code and documentation) under the GNU General Public License. It is maintained by a group of developers from the Universities of Groningen, Uppsala, and Stockholm, and the Max Planck Institute for Polymer Research in Mainz. Its Web site is http://www.gromacs.org.

Peptide and protein folding and conformational equilibria: theoretical treatment of electrostatics and hydrogen bonding with implicit solvent models

Since biomolecules exist in aqueous and membrane environments, the accurate modeling of solvation, and hydrogen bonding interactions in particular, is essential for the exploration of structure and function in theoretical and computational studies. In this chapter, we focus on alternatives to explicit solvent models and discuss recent advances in generalized Born (GB) implicit solvent theories. We present a brief review of the successes and shortcomings of the application of these theories to biomolecular problems that are strongly linked to backbone H-bonding and electrostatics. This discussion naturally leads us to explore existing areas for improvement in current GB theories and our approach towards addressing a number of the key issues that remain in the refinement of these models. Specifically, the critical importance of balancing solvation forces and intramolecular forces in GB models is illustrated by examining the influence of backbone hydrogen bond strength and backbone dihedral energetics on conformational equilibria of small peptids.

Repulsive solvent-induced interaction between C60 fullerenes in water

The role of water-fullerene interactions in the behavior of C60 in aqueous solution was investigated utilizing realistic Lennard-Jones (LJ) and repulsive Weeks-Chandler-Anderson (WCA) potentials. Strong water-fullerene dispersion interactions in the LJ potential dramatically influence the hydration of the fullerene promoting the formation of a high-density hydration shell of water. In contrast to the WCA potential, the water liquid phase between fullerenes remains stable with decreasing fullerene separation, resulting in a repulsive solvent-induced contribution to the fullerene potential of mean force.

Brownian Dynamics, Molecular Dynamics, and Monte Carlo modeling of colloidal systems

This paper serves as an introductory review of Brownian Dynamics (BD), Molecular Dynamics (MD), and Monte Carlo (MC) modeling techniques. These three simulation methods have proven to be exceptional investigative solutions for probing discrete molecular, ionic, and colloidal motions at their basic microscopic levels. The review offers a general study of the classical theories and algorithms that are foundational to Brownian Dynamics, Molecular Dynamics, and Monte Carlo simulations. Important topics of interest include fundamental theories that govern Brownian motion, the Langevin equation, the Verlet algorithm, and the Metropolis method. Brownian Dynamics demonstrates advantages over Molecular Dynamics as pertaining to the issue of time-scale separation. Monte Carlo methods exhibit strengths in terms of ease of implementation. Hybrid techniques that combine these methods and draw from these efficacies are also presented. With their rigorous microscopic approach, Brownian Dynamics, Molecular Dynamics, and Monte Carlo methods prove to be especially viable modeling methods for problems with challenging complexities such as high-level particle concentration and multiple particle interactions. These methods hold promising potential for effective modeling of transport in colloidal systems.

A simple topological representation of protein structure: implications for new, fast, and robust structural classification

A topological representation of proteins is developed that makes use of two metrics: the Euclidean metric for identifying natural nearest neighboring residues via the Delaunay tessellation in Cartesian space and the distance between residues in sequence space. Using this representation, we introduce a quantitative and computationally inexpensive method for the comparison of protein structural topology. The method ultimately results in a numerical score quantifying the distance between proteins in a heuristically defined topological space. The properties of this scoring scheme are investigated and correlated with the standard Calpha distance root-mean-square deviation measure of protein similarity calculated by rigid body structural alignment. The topological comparison method is shown to have a characteristic dependence on protein conformational differences and secondary structure. This distinctive behavior is also observed in the comparison of proteins within families of structural relatives. The ability of the comparison method to successfully classify proteins into classes, superfamilies, folds, and families that are consistent with standard classification methods, both automated and human-driven, is demonstrated. Furthermore, it is shown that the scoring method allows for a fine-grained classification on the family, protein, and species level that agrees very well with currently established phylogenetic hierarchies. This fine classification is achieved without requiring visual inspection of proteins, sequence analysis, or the use of structural superimposition methods. Implications of the method for a fast, automated, topological hierarchical classification of proteins are discussed.

Long range interactions on wires: A reciprocal space based formalism

There are many atomic scale systems in materials, chemistry, and biology that can be effectively modeled as finite in two of the physical spatial dimensions and periodically replicated in the third including nanoscale metallic and semiconducting wires, carbon nanotubes, and DNA. However, it is difficult to design techniques to treat long range forces in these systems without truncation or recourse to slowly convergent supercells or computationally inefficient Poisson solvers. In this paper, a rigorous reciprocal space based formalism which permits long range forces on wires to be evaluated simply and easily via a small modification of existing methods for three dimensional periodicity is derived. The formalism is applied to determine long range interactions both between point particles using an Ewald-like approach and the continuous charge distributions that appear in electronic structure calculations. In this way, both empirical force field calculations and, for example, plane-wave based density functional theory computations on wires can be performed easily. The methodology is tested on model and realistic systems including a lithium doped carbon nanotube.

Computer simulation of solid and liquid benzene with an atomistic interaction potential derived from ab initio calculations

Molecular dynamics atomistic simulations of solid and liquid benzene have been performed, employing a model intermolecular potential derived from quantum mechanical calculations. The ab initio database includes approximately 200 geometries of the benzene dimer with interaction energies computed at the MP2 level of theory. The accuracy of the modeled force field results is satisfactory. The thermodynamic and structural properties, calculated in the condensed phases, are compared with experimental data and previous simulation results. Single particle and collective dynamical properties are also investigated through the calculation of translational and rotational diffusion coefficients, reorientational dynamics, and viscosities. The agreement of these data with experimental measurements confirms the reliability of the proposed force field.

An efficient hybrid explicit/implicit solvent method for biomolecular simulations

We present a new hybrid explicit/implicit solvent method for dynamics simulations of macromolecular systems. The method models explicitly the hydration of the solute by either a layer or sphere of water molecules, and the generalized Born (GB) theory is used to treat the bulk continuum solvent outside the explicit simulation volume. To reduce the computational cost, we implemented a multigrid method for evaluating the pairwise electrostatic and GB terms. It is shown that for typical ion and protein simulations our method achieves similar equilibrium and dynamical observables as the conventional particle mesh Ewald (PME) method. Simulation timings are reported, which indicate that the hybrid method is much faster than PME, primarily due to a significant reduction in the number of explicit water molecules required to model hydration effects.

Reversible temperature and pressure denaturation of a protein fragment: a replica exchange molecular dynamics simulation study

We determine the reversible folding-unfolding of the C-terminal (41-56) fragment of protein G as a function of density and temperature using replica exchange molecular dynamics simulations. We employ a total of 253 replicas, covering the temperature range between 320 and 515 K and the density range between 0.96 and 1.16 g cm(-3). Using the root mean square deviation from the folded structure as a quantitative measure, we are able to obtain the fraction of folded states, and can thus establish the free energy difference between the folded and the unfolded states of the protein fragment as a function of temperature and pressure. For the pressure denaturation the weakening of the hydrophobic interaction between the bulky side chains is found to be crucial at lower temperatures, leading to an apparent destabilization of the folded backbone structure at elevated pressures.

Molecular dynamics simulations of the 136 unique tetranucleotide sequences of DNA oligonucleotides. I. Research design and results on d(CpG) steps

We describe herein a computationally intensive project aimed at carrying out molecular dynamics (MD) simulations including water and counterions on B-DNA oligomers containing all 136 unique tetranucleotide base sequences. This initiative was undertaken by an international collaborative effort involving nine research groups, the "Ascona B-DNA Consortium" (ABC). Calculations were carried out on the 136 cases imbedded in 39 DNA oligomers with repeating tetranucleotide sequences, capped on both ends by GC pairs and each having a total length of 15 nucleotide pairs. All MD simulations were carried out using a well-defined protocol, the AMBER suite of programs, and the parm94 force field. Phase I of the ABC project involves a total of approximately 0.6 mus of simulation for systems containing approximately 24,000 atoms. The resulting trajectories involve 600,000 coordinate sets and represent approximately 400 gigabytes of data. In this article, the research design, details of the simulation protocol, informatics issues, and the organization of the results into a web-accessible database are described. Preliminary results from 15-ns MD trajectories are presented for the d(CpG) step in its 10 unique sequence contexts, and issues of stability and convergence, the extent of quasiergodic problems, and the possibility of long-lived conformational substates are discussed.

Stochastic Liouville equation simulation of multidimensional vibrational line shapes of trialanine

The line shapes detected in coherent femtosecond vibrational spectroscopies contain direct signatures of peptide conformational fluctuations through their effect on vibrational frequencies and intermode couplings. These effects are simulated in trialanine using a Green's function solution of a stochastic Liouville equation constructed for four collective bath coordinates (two Ramachandran angles affecting the mode couplings and two diagonal energies). We find that fluctuations of the Ramachandran angles which hardly affect the linear absorption can be effectively probed by two-dimensional spectra. The signal generated at k(1)+k(2)-k(3) is particularly sensitive to such fluctuations.

Collective diffusion model for water permeation through microscopic channels

Water permeation through nanometric channels, generally a coupled many-body process, is described in this study by a single collective coordinate. A collective diffusion model is proposed in which water movement at equilibrium is characterized as an unbiased diffusion along this coordinate and water transport in the presence of a chemical potential difference is described as one-dimensional diffusion in a linear potential. The model allows one to determine the osmotic permeability of a water channel from equilibrium simulations.

A double-deletion method to quantifying incremental binding energies in proteins from experiment: example of a destabilizing hydrogen bonding pair

The contribution of a specific hydrogen bond in apoflavodoxin to protein stability is investigated by combining theory, experiment and simulation. Although hydrogen bonds are major determinants of protein structure and function, their contribution to protein stability is still unclear and widely debated. The best method so far devised to estimate the contribution of side-chain interactions to protein stability is double mutant cycle analysis, but the interaction energies so derived are not identical to incremental binding energies (the energies quantifying net contributions of two interacting groups to protein stability). Here we introduce double-deletion analysis of 'isolated' residue pairs as a means to precisely quantify incremental binding. The method is exemplified by studying a surface-exposed hydrogen bond in a model protein (Asp96/Asn128 in apoflavodoxin). Combined substitution of these residues by alanines slightly destabilizes the protein due to a decrease in hydrophobic surface burial. Subtraction of this effect, however, clearly indicates that the hydrogen-bonded groups in fact destabilize the native conformation. In addition, molecular dynamics simulations and classic double mutant cycle analysis explain quantitatively that, due to frustration, the hydrogen bond must form in the native structure because when the two groups get approximated upon folding their binding becomes favorable. We would like to remark that 1), this is the first time the contribution of a specific hydrogen bond to protein stability has been measured by experiment; and 2), more hydrogen bonds need to be analyzed to draw general conclusions on protein hydrogen bond energetics. To that end, the double-deletion method should be of help.

Atomic mean-square displacements in proteins by molecular dynamics: a case for analysis of variance

Information on protein internal motions is usually obtained through the analysis of atomic mean-square displacements, which are a measure of variability of the atomic positions distribution functions. We report a statistical approach to analyze molecular dynamics data on these displacements that is based on probability distribution functions. Using a technique inspired by the analysis of variance, we compute unbiased, reliable mean-square displacements of the atoms and analyze them statistically. We applied this procedure to characterize protein thermostability by comparing the results for a thermophilic enzyme and a mesophilic homolog. In agreement with previous experimental observations, our analysis suggests that the proteins surface regions can play a role in the different thermal behavior.

Misfolding pathways of the prion protein probed by molecular dynamics simulations

Although the cellular monomeric form of the benign prion protein is now well characterized, a model for the monomer of the misfolded conformation (PrP(Sc)) remains elusive. PrP(Sc) quickly aggregates into highly insoluble fibrils making experimental structural characterization very difficult. The tendency to aggregation of PrP(Sc) in aqueous solution implies that the monomer fold must be hydrophobic. Here, by using molecular dynamics simulations, we have studied the cellular mouse prion protein and its D178N pathogenic mutant immersed in a hydrophobic environment (solution of CCl4), to reveal conformational changes and/or local structural weaknesses of the prion protein fold in unfavorable structural and thermodynamic conditions. Simulations in water have been also performed. Although observing in general a rather limited conformation activity in the nanosecond timescale, we have detected a significant weakening of the antiparallel beta-sheet of the D178N mutant in CCl4 and to a less extent in water. No weakening is observed for the native prion protein. The increase of beta-structure in the monomer, recently claimed as evidence for misfolding to PrP(Sc), has been also observed in this study irrespective of the thermodynamic or structural conditions, showing that this behavior is very likely an intrinsic characteristic of the prion protein fold.

Phosphatidylethanolamine-phosphatidylglycerol bilayer as a model of the inner bacterial membrane

Phosphatidylethanolamine (PE) and phosphatidylglycerol (PG) are the main lipid components of the inner bacterial membrane. A computer model for such a membrane was built of palmitoyloleoyl PE (POPE) and palmitoyloleoyl PG (POPG) in the proportion 3:1, and sodium ions (Na+) to neutralize the net negative charge on each POPG (POPE-POPG bilayer). The bilayer was simulated for 25 ns. A final 10-ns trajectory fragment was used for analyses. In the bilayer interfacial region, POPEs and POPGs interact readily with one another via intermolecular hydrogen (H) bonds and water bridges. POPE is the main H-bond donor in either PEPE or PEPG H-bonds; PGPG H-bonds are rarely formed. Almost all POPEs are H-bonded and/or water bridged to either POPE or POPG but PE-PG links are favored. In effect, the atom packing in the near-the-interface regions of the bilayer core is tight. Na+ does not bind readily to lipids, and interlipid links via Na+ are not numerous. Although POPG and POPE comprise one bilayer, their bilayer properties differ. The average surface area per POPG is larger and the average vertical location of the POPG phosphate group is lower than those of POPE. Also, the alkyl chains of POPG are more ordered and less densely packed than the POPE chains. The main conclusion of this study is that in the PE-PG bilayer PE interacts more strongly with PG than with PE. This is a likely molecular-level event behind a regulating mechanism developed by the bacteria to control its membrane permeability and stability consisting in changes of the relative PG/PE concentration in the membrane.

Penetratin-membrane association: W48/R52/W56 shield the peptide from the aqueous phase

Using molecular dynamics simulations, we studied the mode of association of the cell-penetrating peptide penetratin with both a neutral and a charged bilayer. The results show that the initial peptide-lipid association is a fast process driven by electrostatic interactions. The homogeneous distribution of positively charged residues along the axis of the helical peptide, and especially residues K46, R53, and K57, contribute to the association of the peptide with lipids. The bilayer enhances the stability of the penetratin helix. Oriented parallel to the lipid-water interface, the subsequent insertion of the peptide through the bilayer headgroups is significantly slower. The presence of negatively charged lipids considerably enhances peptide binding. Lateral side-chain motion creates an opening for the helix into the hydrophobic core of the membrane. The peptide aromatic residues form a pi-stacking cluster through W48/R52/W56 and F49/R53, protecting the peptide from the water phase. Interaction with the penetratin peptide has only limited effect on the overall membrane structure, as it affects mainly the conformation of the lipids which interact directly with the peptide. Charge matching locally increases the concentration of negatively charged lipids, lateral lipid diffusion locally decreases. Lipid disorder increases, through decreased order parameters of the lipids interacting with the penetratin side chains. Penetratin molecules at the membrane surface do not seem to aggregate.

Dynamic water networks in cytochrome C oxidase from Paracoccus denitrificans investigated by molecular dynamics simulations

We present a molecular dynamics study of cytochrome c oxidase from Paracoccus denitrificans in the fully oxidized state, embedded in a fully hydrated dimyristoylphosphatidylcholine lipid bilayer membrane. Parallel simulations with different levels of protein hydration, 1.125 ns each in length, were carried out under conditions of constant temperature and pressure using three-dimensional periodic boundary conditions and full electrostatics to investigate the distribution and dynamics of water molecules and their corresponding hydrogen-bonded networks inside cytochrome c oxidase. The majority of the water molecules had residence times shorter than 100 ps, but a few water molecules are fixed inside the protein for up to 1.125 ns. The hydrogen-bonded network in cytochrome c oxidase is not uniformly distributed, and the degree of water arrangement is variable. The average number of solvent sites in the proton-conducting K- and D-pathways was determined. In contrast to single water files in narrow geometries we observe significant diffusion of individual water molecules along these pathways. The highly fluctuating hydrogen-bonded networks, combined with the significant diffusion of individual water molecules, provide a basis for the transfer of protons in cytochrome c oxidase, therefore leading to a better understanding of the mechanism of proton pumping.

Experimental validation of molecular dynamics simulations of lipid bilayers: a new approach

A novel protocol has been developed for comparing the structural properties of lipid bilayers determined by simulation with those determined by diffraction experiments, which makes it possible to test critically the ability of molecular dynamics simulations to reproduce experimental data. This model-independent method consists of analyzing data from molecular dynamics bilayer simulations in the same way as experimental data by determining the structure factors of the system and, via Fourier reconstruction, the overall transbilayer scattering-density profiles. Multi-nanosecond molecular dynamics simulations of a dioleoylphosphatidylcholine bilayer at 66% RH (5.4 waters/lipid) were performed in the constant pressure and temperature ensemble using the united-atom GROMACS and the all-atom CHARMM22/27 force fields with the GROMACS and NAMD software packages, respectively. The quality of the simulated bilayer structures was evaluated by comparing simulation with experimental results for bilayer thickness, area/lipid, individual molecular-component distributions, continuous and discrete structure factors, and overall scattering-density profiles. Neither the GROMACS nor the CHARMM22/27 simulations reproduced experimental data within experimental error. The widths of the simulated terminal methyl distributions showed a particularly strong disagreement with the experimentally observed distributions. A comparison of the older CHARMM22 with the newer CHARMM27 force fields shows that significant progress is being made in the development of atomic force fields for describing lipid bilayer systems empirically.

Simulation of the early stages of nano-domain formation in mixed bilayers of sphingomyelin, cholesterol, and dioleylphosphatidylcholine

It is known from experimental studies that lipid bilayers composed of unsaturated phospholipids, sphingomyelin, and cholesterol contain microdomains rich in sphingomyelin and cholesterol. These domains are similar to "rafts" isolated from cell membranes, although the latter are much smaller in lateral size. Such domain formation can be a result of very specific and subtle lipid-lipid interactions. To identify and study these interactions, we have performed two molecular dynamics simulations, of 200-ns duration, of dioleylphosphatidylcholine (DOPC), sphingomyelin (SM), and cholesterol (Chol) systems, a 1:1:1 mixture of DOPC/SM/Chol, and a 1:1 mixture of DOPC/SM. The simulations show initial stages of the onset of spontaneous phase-separated domains in the systems. On the simulation timescale cholesterol favors a position at the interface between the ordered SM region and the disordered DOPC region in the ternary system and accelerates the process of domain formation. We find that the smooth alpha-face of Chol preferentially packs next to SM molecules. Based on a comparative analysis of interaction energies, we find that Chol molecules do not show a preference for SM or DOPC. We conclude that Chol molecules assist in the process of domain formation and the process is driven by entropic factors rather than differences in interaction energies.

Toc12, a novel subunit of the intermembrane space preprotein translocon of chloroplasts

Translocation of proteins across membranes is essential for the biogenesis of each cell and is achieved by proteinaceous complexes. We analyzed the translocation complex of the intermembrane space from chloroplasts and identified a 12-kDa protein associated with the Toc machinery. Toc12 is an outer envelope protein exposing a soluble domain into the intermembrane space. Toc12 contains a J-domain and stimulates the ATPase activity of DnaK. The conformational stability and the ability to stimulate Hsp70 are dependent on a disulfide bridge within the loop region of the J-domain, suggesting a redox-regulated activation of the chaperone. Toc12 is associated with Toc64 and Tic22. Its J-domain recruits the Hsp70 of outer envelope membrane to the intermembrane space translocon and facilitates its interaction to the preprotein.

Model for the dynamics of a water cluster in an x-ray free electron laser beam

A microscopic sample placed into a focused x-ray free electron laser beam will explode due to strong ionization on a femtosecond time scale. The dynamics of this Coulomb explosion has been modeled by Neutze et al. [Nature (London) 406, 752 (2000)] for a protein, using computer simulations. The results suggest that by using ultrashort exposures, structural information may be collected before the sample is destroyed due to radiation damage. In this paper a method is presented to include the effect of screening by free electrons in the sample in a molecular dynamics simulation. The electrons are approximated by a classical gas, and the electron distribution is calculated iteratively from the Poisson-Boltzmann equation. Test simulations of water clusters reveal the details of the explosion dynamics, as well as the evolution of the free electron gas during the beam exposure. We find that inclusion of the electron gas in the model slows down the Coulomb explosion. The hydrogen atoms leave the sample faster than the oxygen atoms, leading to a double layer of positive ions. A considerable electron density is located between these two layers. The fact that the hydrogens are found to explode much faster than the oxygens means that the diffracting part of the sample stays intact somewhat longer than the sample as a whole.

HIV-1 protease molecular dynamics of a wild-type and of the V82F/I84V mutant: possible contributions to drug resistance and a potential new target site for drugs

The protease from type 1 human immunodeficiency virus (HIV-1) is a critical drug target against which many therapeutically useful inhibitors have been developed; however, the set of viral strains in the population has been shifting to become more drug-resistant. Because indirect effects are contributing to drug resistance, an examination of the dynamic structures of a wild-type and a mutant could be insightful. Consequently, this study examined structural properties sampled during 22 nsec, all atom molecular dynamics (MD) simulations (in explicit water) of both a wild-type and the drug-resistant V82F/I84V mutant of HIV-1 protease. The V82F/I84V mutation significantly decreases the binding affinity of all HIV-1 protease inhibitors currently used clinically. Simulations have shown that the curling of the tips of the active site flaps immediately results in flap opening. In the 22-nsec MD simulations presented here, more frequent and more rapid curling of the mutant's active site flap tips was observed. The mutant protease's flaps also opened farther than the wild-type's flaps did and displayed more flexibility. This suggests that the effect of the mutations on the equilibrium between the semiopen and closed conformations could be one aspect of the mechanism of drug resistance for this mutant. In addition, correlated fluctuations in the active site and periphery were noted that point to a possible binding site for allosteric inhibitors.

Complexation of phosphatidylcholine lipids with cholesterol

It is postulated that the specific interactions between cholesterol and lipids in biological membranes are crucial in the formation of complexes leading subsequently to membrane domains (so-called rafts). These interactions are studied in molecular dynamics simulations performed on a dipalmitoylphosphatidylcholine (DPPC)-cholesterol bilayer mixture and a dilauroylphosphatidylcholine (DLPC)-cholesterol bilayer mixture, both having a cholesterol concentration of 40 mol %. Complexation of the simulated phospholipids with cholesterol is observed and visualized, exhibiting 2:1 and 1:1 stoichiometries. The most popular complex is found to be 1:1 in the case of DLPC, whereas the DPPC system carries a larger population of 2:1 complexes. This difference in the observed populations of complexes is shown to be a result of differences in packing geometry and phospholipid conformation due to the differing tail length of the two phosphatidylcholine lipids. Furthermore, aggregation of these complexes appears to form hydrogen-bonded networks in the system containing a mixture of cholesterol and DPPC. The CH...O hydrogen bond plays a crucial role in the formation of these complexes as well as the hydrogen bonded aggregates. The aggregation and extension of such a network implies a possible means by which phospholipid:cholesterol domains form.

Molecular dynamics simulation of a palmitoyl-oleoyl phosphatidylserine bilayer with Na+ counterions and NaCl

Two 40 ns molecular dynamics simulations of a palmitoyl-oleoyl phosphatidylserine (POPS) lipid bilayer in the liquid crystalline phase with Na(+) counterions and NaCl were carried out to investigate the structure of the negatively charged lipid bilayer and the effect of salt (NaCl) on the lipid bilayer structure. Na(+) ions were found to penetrate deep into the ester region of the water/lipid interface of the bilayer. Interaction of the Na(+) ions with the lipid bilayer is accompanied by a loss of water molecules around the ion and a simultaneous increase in the number of ester carbonyl oxygen atoms binding the ion, which define an octahedral and square pyramidal geometry. The amine group of the lipid molecule is involved in the formation of inter- and intramolecular hydrogen bonds with the carboxylate and the phosphodiester groups of the lipid molecule. The area per lipid of the POPS bilayer is unaffected by the presence of 0.15M NaCl. There is a small increase in the order parameter of carbon atoms in the beginning of the alkyl chain in the presence of NaCl. This is due to a greater number of Na(+) ions being coordinated by the ester carbonyl oxygen atoms in the water/lipid interface region of the POPS bilayer.

Implementation and testing of stable, fast implicit solvation in molecular dynamics using the smooth-permittivity finite difference Poisson-Boltzmann method

A fast stable finite difference Poisson-Boltzmann (FDPB) model for implicit solvation in molecular dynamics simulations was developed using the smooth permittivity FDPB method implemented in the OpenEye ZAP libraries. This was interfaced with two widely used molecular dynamics packages, AMBER and CHARMM. Using the CHARMM-ZAP software combination, the implicit solvent model was tested on eight proteins differing in size, structure, and cofactors: calmodulin, horseradish peroxidase (with and without substrate analogue bound), lipid carrier protein, flavodoxin, ubiquitin, cytochrome c, and a de novo designed 3-helix bundle. The stability and accuracy of the implicit solvent simulations was assessed by examining root-mean-squared deviations from crystal structure. This measure was compared with that of a standard explicit water solvent model. In addition we compared experimental and calculated NMR order parameters to obtain a residue level assessment of the accuracy of MD-ZAP for simulating dynamic quantities. Overall, the agreement of the implicit solvent model with experiment was as good as that of explicit water simulations. The implicit solvent method was up to eight times faster than the explicit water simulations, and approximately four times slower than a vacuum simulation (i.e., with no solvent treatment).

Plasticity of influenza haemagglutinin fusion peptides and their interaction with lipid bilayers

A detailed molecular dynamics study of the haemagglutinin fusion peptide (N-terminal 20 residues of the HA2 subunits) in a model bilayer has yielded useful information about the molecular interactions leading to insertion into the lipids. Simulations were performed on the native sequence, as well as a number of mutant sequences, which are either fusogenic or nonfusogenic. For the native sequence and fusogenic mutants, the N-terminal 11 residues of the fusion peptides are helical and insert with a tilt angle of approximately 30 degrees with respect to the membrane normal, in very good agreement with experimental data. The tilted insertion of the native sequence peptide leads to membrane bilayer thinning and the calculated order parameters show larger disorder of the alkyl chains. These results indicate that the lipid packing is perturbed by the fusion peptide and could be used to explain membrane fusion. For the nonfusogenic sequences investigated, it was found that most of them equilibrate parallel to the interface plane and do not adopt a tilted conformation. The presence of a charged residue at the beginning of the sequence (G1E mutant) resulted in a more difficult case, and the outcomes do not fall straightforwardly into the general picture. Sequence searches have revealed similarities of the fusion peptide of influenza haemagglutinin with peptide sequences such as segments of porin, amyloid alpha eta peptide, and a peptide from the prion sequence. These results confirm that the sequence can adopt different folds in different environments. The plasticity and the conformational dependence on the local environment could be used to better understand the function of fusion peptides.

Long time molecular dynamics for enhanced conformational sampling in biomolecular systems

Although molecular dynamics methods are commonly used to drive biomolecular simulations, the technique provides insufficient sampling to impact studies of the 200-300 residue proteins of greatest interest. One severe limitation of molecular dynamics is that the integrators are restricted by resonance phenomena to small time steps (Delta t<8 fs) much slower then the time scales of important structural and solvent rearrangements. Here, a novel set of equations of motion and a reversible, resonance-free, integrator are designed which permit step sizes on the order of 100 fs to be used.

Lipid-protein interactions of integral membrane proteins: a comparative simulation study

The interactions between membrane proteins and their lipid bilayer environment play important roles in the stability and function of such proteins. Extended (15-20 ns) molecular dynamics simulations have been used to explore the interactions of two membrane proteins with phosphatidylcholine bilayers. One protein (KcsA) is an alpha-helix bundle and embedded in a palmitoyl oleoyl phosphatidylcholine bilayer; the other (OmpA) is a beta-barrel outer-membrane protein and is in a dimyristoyl phosphatidylcholine bilayer. The simulations enable analysis in detail of a number of aspects of lipid-protein interactions. In particular, the interactions of aromatic amphipathic side chains (i.e., Trp, Tyr) with lipid headgroups, and "snorkeling" interactions of basic side chains (i.e., Lys, Arg) with phosphate groups are explored. Analysis of the number of contacts and of H-bonds reveal fluctuations on an approximately 1- to 5-ns timescale. There are two clear bands of interacting residues on the surface of KcsA, whereas there are three such bands on OmpA. A large number of Arg-phosphate interactions are seen for KcsA; for OmpA, the number of basic-phosphate interactions is smaller and shows more marked fluctuations with respect to time. Both classes of interaction occur in clearly defined interfacial regions of width approximately 1 nm. Analysis of lateral diffusion of lipid molecules reveals that "boundary" lipid molecules diffuse at about half the rate of bulk lipid. Overall, these simulations present a dynamic picture of lipid-protein interactions: there are a number of more specific interactions but even these fluctuate on an approximately 1- to 5-ns timescale.

The range and shielding of dipole-dipole interactions in phospholipid bilayers

In molecular dynamics simulations of lipid bilayers, the structure is sensitive to the precise treatment of electrostatics. The dipole-dipole interactions between headgroup dipoles are not long-ranged, but the area per lipid and, through it, other properties of the bilayer are very sensitive to the detailed balance between the perpendicular and in-plane components of the headgroup dipoles. This is affected by the detailed properties of the cutoff scheme or if long-range interactions are included by Ewald or particle-mesh Ewald techniques. Interaction between the in-plane components of the headgroup dipoles is attractive and decays as the inverse sixth power of distance. The interaction is screened by the square of a dielectric permittivity close to the value for water. Interaction between the components perpendicular to the membrane plane is repulsive and decays as the inverse third power of distance. These interactions are screened by a dielectric permittivity of the order 10. Thus, despite the perpendicular components being much smaller in magnitude than the in-plane components, they will dominate the interaction energies at large distances.

Atomic force algorithms in density functional theory electronic-structure techniques based on local orbitals

Electronic structure methods based on density-functional theory, pseudopotentials, and local-orbital basis sets offer a hierarchy of techniques for modeling complex condensed-matter systems with a wide range of precisions and computational speeds. We analyze the relationships between the algorithms for atomic forces in this hierarchy of techniques, going from empirical tight-binding through ab initio tight-binding to full ab initio. The analysis gives a unified overview of the force algorithms as applied within techniques based either on diagonalization or on linear-scaling approaches. The use of these force algorithms is illustrated by practical calculations with the CONQUEST code, in which different techniques in the hierarchy are applied in a concerted manner.

Diffusion and electrophoretic mobility of single-stranded RNA from molecular dynamics simulations

Hydrodynamic properties of small single-stranded RNA homopolymers with three and six nucleotides in free solution are determined from molecular dynamics simulations in explicit solvent. We find that the electrophoretic mobility increases with increasing RNA length, consistent with experiment. Diffusion coefficients of RNA, corrected for finite-size effects and solvent viscosity, agree well with those estimated from experiments and hydrodynamic calculations. The diffusion coefficients and electrophoretic mobilities satisfy a Nernst-Einstein relation in which the effective charge of RNA is reduced by the charge of transiently bound counterions. Fluctuations in the counterion atmosphere are shown to enhance the diffusive spread of RNA molecules drifting along the direction of the external electric field. As a consequence, apparent diffusion coefficients measured by capillary zone electrophoresis can be significantly larger than the actual values at certain experimental conditions.

Phase behavior of polarizable spherocylinders in external fields

Applied electric fields are known to induce significant changes in the properties of systems of polarizable molecules or particles. For rod-shaped molecules, the field-induced behavior can be rather surprising, as in the case of the negative electric birefringence of concentrated solutions of rodlike polyelectrolytes. We have investigated the interplay of shape anisotropy and field-induced anisotropy in molecular dynamics simulations of systems of polarizable soft spherocylinders in an electric field, in the limit of infinitely anisotropic polarizability, taking full account of mutual induction effects. We find a novel crystalline structure (K(2)) in the high-field limit, whose formation is driven by interactions between induced dipoles. For high pressures, the phase diagram exhibits a polar nematic phase between the hexagonal close-packed crystal phase and the K(2) phase. We also compare this system with an analogous system of spherocylinders with permanent electric dipoles and find that qualitatively similar behavior is obtained in the limit of strong coupling of the permanent dipoles to the external field.

DL_POLY_3: the CCP5 national UK code for molecular-dynamics simulations

DL_POLY_3 is a general-purpose molecular-dynamics simulation package embedding a highly efficient domain decomposition (DD) parallelization strategy. It was developed at Daresbury Laboratory under the auspices of the Engineering and Physical Sciences Research Council. Written to support academic research, it has a wide range of applications and will run on a wide range of computers; from single-processor workstations to multi-processor computers, with accent on the efficient use of multi-processor power. A new DD adaptation of the smoothed particle mesh Ewald method for calculating long-range forces in molecular simulations, incorporating a novel three-dimensional fast Fourier transform (the Daresbury Advanced Fourier Transform), makes it possible to simulate systems of the order of one million particles and beyond. DL_POLY_3 structure, functionality, performance and availability are described in this feature paper.

Competition of hydrophobic and Coulombic interactions between nanosized solutes

The solvation of charged, nanometer-sized spherical solutes in water, and the effective, solvent-induced force between two such solutes are investigated by constant temperature and pressure molecular dynamics simulations of model solutes carrying various charge patterns. The results for neutral solutes agree well with earlier findings, and with predictions of simple macroscopic considerations: substantial hydrophobic attraction may be traced back to strong depletion ("drying") of the solvent between the solutes. This hydrophobic attraction is strongly reduced when the solutes are uniformly charged, and the total force becomes repulsive at sufficiently high charge; there is a significant asymmetry between anionic and cationic solute pairs, the latter experiencing a lesser hydrophobic attraction. The situation becomes more complex when the solutes carry discrete (rather than uniform) charge patterns. Due to antagonistic effects of the resulting hydrophilic and hydrophobic "patches" on the solvent molecules, water is once more significantly depleted around the solutes, and the effective interaction reverts to being mainly attractive, despite the direct electrostatic repulsion between solutes. Examination of a highly coarse-grained configurational probability density shows that the relative orientation of the two solutes is very different in explicit solvent, compared to the prediction of the crude implicit solvent representation. The present study strongly suggests that a realistic modeling of the charge distribution on the surface of globular proteins, as well as the molecular treatment of water, are essential prerequisites for any reliable study of protein aggregation.

Direct simulation of transmembrane helix association: role of asparagines

The forces contributing to the association of transmembrane helices in folded membrane proteins have received considerable attention recently. In this study we investigate the importance of hydrogen bonding by studying the effect of a single Asn residue in the center of an otherwise hydrophobic transmembrane peptide using computer simulations. We use the model peptide MS1 which has been derived from the leucine zipper coiled-coil dimer of the transcription factor peptide GCN4-P1. We follow the trajectory of 36 initially monomeric MS1 transmembrane helical peptides in a membrane-mimicking octane layer as they associate into larger structures. These peptides predominately form dimers. The interaction between the polar asparagine residues, capable of simultaneously being a hydrogen-bond donor and acceptor, contributes strongly to the stability of associated helices. Only dimers with interhelical hydrogen bonds form stable structures, whereas aggregates without any hydrogen-bonding interactions form very transient structures. We examine the hydrogen-bonding patterns and find that there are two forms of dimer, one with symmetric hydrogen bonds and one with asymmetric hydrogen bonds. Based on the structures in our simulation we propose a model with a monomer <--> symmetric dimer <--> asymmetric dimer <--> trimer equilibrium.

Molecular dynamics simulation study of the interaction of trehalose with lipid membranes

The interactions of the cryoprotective agent trehalose with a lipid membrane made of 1,2-dipalmitoyl-sn-glycero-3-phosphocholine at 323 K were studied by means of molecular dynamics simulations. It was observed that trehalose binds to the phospholipid headgroups with its main axis parallel to the membrane normal. Trehalose establishes hydrogen bonds with the carbonyl and phosphate groups and replaces water molecules from the lipid headgroup. Notably, the number of hydrogen bonds (HBs) that the membrane made with its environment was conserved after trehalose binding. The HBs between lipid and trehalose have a longer lifetime than those established between lipid and water. The binding of the sugar does not produce changes either in the lipid area or in the lipid order parameter. The effect of trehalose on the dipole potential is in agreement with experimental results. The contribution of the different components to the membrane dipole potential was analyzed. It was observed that the binding of trehalose produces changes in the different components and the sugar itself contributes to the surface potential due to the polarization of its hydroxyl in the interface.

Mesoscopic lateral diffusion in lipid bilayers

The lateral diffusion in bilayers is modeled with a multiscale mesoscopic simulation. The methodology consists of two simulations, where the first employs atomistic models to obtain exact results for the mesoscopic model. The second simulation takes the results obtained from the first to parameterize an effective force field that is employed in a new coarse-grained model. The multiscale aspect of this scheme occurs at the point where the microscopic time-averaged results of the first simulation are employed to parameterize the second simulation that operates in a higher spatial and temporal domain. The results of both simulation schemes give quantitative information on the details associated with lipid lateral diffusion.