GROMACS Simulation Package Research Articles

Probing microscopic material properties inside simulated membranes through spatially resolved three-dimensional local pressure fields and surface tensions

Cellular lipid membranes are spatially inhomogeneous soft materials. Materials properties such as pressure and surface tension thus show important microscopic-scale variation that is critical to many biological functions. We present a means to calculate pressure and surface tension in a 3D-resolved manner within molecular-dynamics simulations and show how such measurements can yield important insight. We also present the first corrections to local virial and pressure fields to account for the constraints typically used in lipid simulations that otherwise cause problems in highly oriented systems such as bilayers. Based on simulations of an asymmetric bacterial ion channel in a POPC bilayer, we demonstrate how 3D-resolved pressure can probe for both short-range and long-range effects from the protein on the membrane environment. We also show how surface tension is a sensitive metric for inter-leaflet equilibrium and can be used to detect even subtle imbalances between bilayer leaflets in a membrane-protein simulation. Since surface tension is known to modulate the function of many proteins, this effect is an important consideration for predictions of ion channel function. We outline a strategy by which our local pressure measurements, which we make available within a version of the GROMACS simulation package, may be used to design optimally equilibrated membrane-protein simulations.

Study of Nanoparticle-Lipid Bilayer Interactions: Insights from Coarse-Grained Molecular Dynamics Simulations

One of the main aims in the design of engineered nanomaterials with applications in medicine, such as diagnostic and therapeutic nanoparticles (NPs), is the ability of these materials to translocate across the human cells without damaging essential tissues. The entry point of a NP to a cell is the plasma membrane. Thus, the first step into assessing the NP cytotoxicity requires a thorough understanding of the NP-membrane interaction mechanism. Extensive Molecular Dynamics (MD) simulations and free-energy calculations were employed, providing insights into the significance of NP surface chemistry and cholesterol concentration of the membrane in the NP-membrane interplay. The MARTINI coarse-grained force-field as implemented in the GROMACS simulation package was employed to model NPs and investigate their interaction with a model lipid bilayer. NPs with a diameter of 3 nm and a surface comprising different patterns of hydrophobic and hydrophilic groups were designed. MD equilibrium simulations were carried out in order to gain insight into the mechanism of NP-membrane interaction. To evaluate the molecular energetics of this interaction, Potential of Mean Force (PMF) calculations, using the Umbrella Sampling technique, were also performed. Our results demonstrate that certain surface patterns, such as an ordered distribution of hydrophilic groups, alter the interaction mechanisms and change the molecular energetics of the NP-membrane interactions. Moreover, PMF calculations showed that an increase in the cholesterol concentration of the membrane leads to a higher energy barrier of NP translocation across the membrane. The results from the present study provide valuable information on possible links between NP surface characteristics and cholesterol concentration with specific interaction types with the cell membrane, and provide insights into the design of NPs with tailored functionalities, for example direct cellular entry.

A Computational Investigation of the Effect of Membrane Curvature on G-Protein Coupled Receptor Oligomerization

Typical cellular membranes display wide structural and rheological diversity. Such properties have been suggested to impact the spatial organization of membrane proteins with consequent functional implications. Indeed, oligomers of G protein-coupled receptors (GPCRs) are actively investigated for their role in cell physiology and as promising new targets of improved clinical relevance. Inspired by recent experimental findings, we used microsecond-scale molecular dynamics (MD) simulations to investigate the effect of membrane curvature on the mobility and spatial organization of coarse-grained (CG) representations of the prototypical GPCR beta2 adrenergic receptor (B2AR) in a mixture of 80% dioleoylphosphatidylcholine (DOPC), 10% dioleoylphosphatidylglycerol (DOPG), and 10% cholesterol used as a membrane mimetic. Specifically, we explored the self-association of 16 CG B2AR molecules in either a flat lipid bilayer or a closed vesicle with a ∼32 nm diameter at a protein/lipid-cholesterol ratio of ∼1:450. All simulations were performed with the GROMACS simulation package using the 4 atoms-to-1 bead MARTINI force field. In agreement with simple hydrodynamic models, we find that receptor mobility is hindered at large curvature. A detailed analysis of receptor collisions and residence time of dimeric complexes during simulations yields results consistent with an increased, more stable self-association of receptors in the planar bilayer compared to the small vesicle system.

Combination of the CHARMM27 force field with united‐atom lipid force fields

Computer simulations offer a valuable way to study membrane systems, from simple lipid bilayers to large transmembrane protein complexes and lipid-nucleic acid complexes for drug delivery. Their accuracy depends on the quality of the force field parameters used to describe the components of a particular system. We have implemented the widely used CHARMM22 and CHARMM27 force fields in the GROMACS simulation package to (i) combine the CHARMM22 protein force field with two sets of united-atom lipids parameters; (ii) allow comparisons of the lipid CHARMM27 force field with other lipid force fields or lipid-protein force field combinations. Our tests do not show any particular issue with the combination of the all-atom CHARMM22 force field with united-atoms lipid parameters, although pertinent experimental data are lacking to assess the quality of the lipid-protein interactions. The conversion utilities allow automatic generation of GROMACS simulation files with CHARMM nucleic acids and protein parameters and topologies, starting from pdb files using the standard GROMACS pdb2gmx method. CMAP is currently not implemented.

Molecular simulation of structure, thermodynamic and transport properties of polyacrylonitrile, polystyrene and their alternating copolymers in high temperatures

This manuscript focuses on the molecular simulation of an industrially important copolymer, namely poly(styrene-co-acrylonitrile) (usually known with the acronym SAN). This polymer is glassy at ambient temperatures with glass transition temperature around 100–110°C. In our case, the copolymer is studied at rather high temperature, in the range 150–300°C, corresponding to industrial processing conditions.In the initial stage of the work, emphasis was put on the molecular dynamics (MD) code development for the simulation of SAN and on the development of an optimum atomistic force field. Polar interactions were accounted implicitly through appropriate adjustment of the Lennard–Jones parameters. All simulations and optimizations were performed using the GROMACS simulation package that is used widely by the polymer and biopolymer research community. Both homo-polymers (polyacrylonitrile and polystyrene) and the copolymer at different compositions were simulated. The force field was tuned to high temperature experimental data at 1bar and used subsequently for simulations at higher pressures. Volumetric and structural properties of the homo-polymers and copolymers were calculated. Simulation results are in very good agreement with limited experimental data available.In the second stage of the work, the solubility of five gases (namely He, Ne, Ar, Kr and Xe), four n-alkanes (namely CH4, C3H8, n-C4H10 and n-C5H12) and toluene in the homo-polymers and copolymers was calculated using the Widom particle insertion methodology in the temperature range 450–550K. Finally, the diffusion coefficient of acrylonitrile and styrene in the same polymers and temperatures was calculated.Due to lack of available experimental data for the mixtures, the accuracy of the reported calculations cannot be quantitatively assessed. Nevertheless, the trends observed with respect to the effect of temperature and polymer chemical structure on the physical properties evaluated should be useful for future polymer selection for specific applications. Clearly, additional experimental data are needed in order to refine the models developed here and validate the predictions reported.

An Improved United Atom Force Field for Simulation of Mixed Lipid Bilayers

We introduce a new force field (43A1-S3) for simulation of membranes by the Gromacs simulation package. Construction of the force fields is by standard methods of electronic structure computations for bond parameters and charge distribution and specific volumes and heats of vaporization for small-molecule components of the larger lipid molecules for van der Waals parameters. Some parameters from the earlier 43A1 force field are found to be correct in the context of these calculations, while others are modified. The validity of the force fields is demonstrated by correct replication of X-ray form factors and NMR order parameters over a wide range of membrane compositions in semi-isotropic NTP 1 atm simulations. 43-A1-S3 compares favorably with other force fields used in conjunction with the Gromacs simulation package with respect to the breadth of phenomena that it accurately reproduces.

Conformational Transitions and Proton Conduction in the Multidrug Efflux Pump AcrB

The increasing problem of multi-drug resistance (MDR) in cancer therapy or bacterial infections is to a large degree caused by multi-drug efflux pumps in pathogenic cells. In Escherichia coli, a major resistance mechanism against antibiotics is based on a tripartite multi-drug export complex comprising the inner membrane translocase AcrB, the membrane-fusion protein AcrA and the outer-membrane channel TolC. AcrB functions as the engine of this complex, using proton motive force to expel a wide variety of unrelated toxic compounds such as antibiotics, disinfectants or detergents. The molecular mechanisms of how proton conduction through AcrB is coupled to drug expulsion are not fully understood yet. Here we report a combination of normal mode analysis (NMA) and molecular dynamics (MD) simulation to investigate conformational transitions occurring in the AcrB reaction cycle and to identify residues crucial for proton conduction. In the crystallographic structure of AcrB each monomer is trapped in a different conformation, representing consecutive states in the transport mechanism. Applying the elastic network NMA variant of minimum action pathway (Kim et al. 2002), we computed transitions between these states. The resulting c-alpha trajectories were then converted back to all atom in an approach of steered energy minimization. We also performed multi-copy MD simulations of AcrB embedded in a phospholipid/water environment using the GroMACS simulation package. Mapping the proton conduction pathway was done on the basis of protein-internal water dynamics and monitoring their frequency of forming hydrogen bonds to adjacent residues.

G_permute: Permutation-reduced phase space density compaction

Biomolecular processes are governed by free energy changes and thus depend on a fine-tuned interplay between entropy and enthalpy. To calculate accurate values for entropies from simulations is particularly challenging for the solvation shell of proteins, which contributes crucially to the total entropy of solvated proteins, due to the diffusive motion of the solvent molecules. Accordingly, for each frame of a Molecular dynamics (MD) trajectory, our software relabels the solvent molecules, such that the resulting configuration space volume is reduced by a factor of N! with N being the number of solvent molecules. The combinatorial explosion of a naive implementation is here overcome by transforming the task into a linear assignment problem, for which algorithms with complexity O ( N 3 ) exist. We have shown in previous research that the solvent entropy can be estimated from such a compacted trajectory by established entropy estimation methods. In this paper, we describe the software implementation which also allows applications beyond entropy estimation, such as the permutation of lipids in membrane bilayers. Program summary Program title: g_permute Catalogue identifier: AECJ_v1_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AECJ_v1_0.html Program obtainable from: CPC Program Library, Queen's University, Belfast, N. Ireland Licensing provisions: GPL No. of lines in distributed program, including test data, etc.: 45 173 No. of bytes in distributed program, including test data, etc.: 2 730 678 Distribution format: tar.gz Programming language: C Computer: PC-compatible running Linux Operating system: Linux RAM: Dependent on the number of solvent molecules, min 12 582 912 bytes Classification: 3, 4.8, 4.9 External routines: liblap (included); From GROMACS-3.3.1: libgmx (not included) Nature of problem: Estimating the entropy of solvent molecules from a molecular dynamics simulation trajectory cannot be performed on ordinary trajectories. Solution method: Compacting the configuration space of molecules by exploiting their permutation symmetry. Applies to trajectories either compatible to those obtained with the GROMACS simulation package [1] or multi-model pdb (Protein Data Bank) files. Restrictions: In rare cases the time to find a solution for the linear assignment problem can be very long. Running time: Dependent on trajectory length and number of molecules to be permuted.

Physicochemical properties of blue fluorescent protein determined via molecular dynamics simulation

Blue fluorescent protein (BFP) is a mutant of green fluorescent protein (GFP), where the chromophore has been modified to shift the emitted fluorescence into the blue spectral region. In this study, MD calculations were performed with the GROMACS simulation package and AMBER force field to investigate the dependence of BFPs physicochemical properties on temperature and applied pressure. The MD approach enabled us to calculate the compressibility of protein itself, separately from the nontrivial contribution of the hydration shell, which is difficult to achieve experimentally. The computed compressibility of BFP (3.94 x10(-5) MPa(-1)) is in agreement with experimental values of globular proteins. The center-of-mass diffusion coefficient of BFP and its dependence on temperature and pressure, which plays an important role in its application as a probe for intracellular liquid viscosity measurement, was calculated and found to be in good agreement with photobleaching recovery experimental data. We have shown that decreased temperature as well as applied pressure increases the water viscosity, but the concomitant decrease of the BFP diffusion coefficient behaves differently from Stokes-Einstein formula. It is shown that the number of hydrogen bonds around the protein grows with pressure, which explains the aforementioned deviation. Pressure also reduces root mean square (RMS) fluctuations, especially those of the most flexible residues situated in the loops. The analysis of the RMS fluctuations of the backbone Calpha atoms also reveals that the most rigid part of BFP is the center of the beta-barrel, in accord with temperature B factors obtained from the Protein Data Bank.

Functional and Computational Analysis of the Role of the Adjacent to Metal Ion-Dependent Adhesion Site (ADMIDAS) of Integrin αIIbβ3.

The role of the β3 MIDAS in αIIbβ3 ligand binding is well established, but the role of the nearby ADMIDAS is less well defined. Thus, we studied HEK293 cells expressing normal αIIbβ3 (normal cells) or the ADMIDAS mutants β3 D126A and D127A (mutant cells). Both mutant cells adhered as well or better than normal cells to immobilized fibrinogen under static conditions in the presence of either Ca2+/Mg2+ or Mn2+. Under low shear flow conditions (0.15 dyne/cm2), adhesion of normal cells and D126A mutant cells to fibrinogen was similar in the presence of either Ca2+/Mg2+ or Mn2+. Adherent D126A mutant cells, however, demonstrated greater resistance to detachment at increasing shear rates in the presence of Ca2+/Mg2+ (e.g., at 20.4 dynes/cm2, only 40 ± 10% of normal cells remained vs 85 ± 8% of D126A mutant cells; mean ± SD; p<0.001). Substituting Mn2+ for Ca2+/Mg2+ increased the resistance to detachment of the normal cells (60 ± 20% remaining at 20.4 dynes/cm2; p=0.01), but the value was still less than the mutant cells in the presence of either Ca2+/Mg2+ (see above; p<0.01) or Mn2+ (84 ± 4%; p<0.01). The increased strength of adhesion we observed in the αIIbβ3 ADMIDAS mutant cells is similar to that found in α4β7 ADMIDAS mutant cells (Chen et al, JBC 2004) and is consistent with the findings in isolated β3 βA (I-like) domains (Pesho et al. JBC 2006). The binding of 7E3, whose epitope is near the ADMIDAS, to the D126 mutant cells was similar to its binding to the normal cells, but 7E3 binding to the D127A mutant cells was reduced by 89 ± 7% (n = 4; p<0.001). 7E3 decreased adhesion of normal cells to fibrinogen by 88 ± 4%, but it only decreased adhesion of D126A mutant cells by 3 ± 9%, and it did not inhibit adhesion of D127A cells at all. To provide a structural context for the role of the ADMIDAS in ligand binding to αIIbβ3, we compared results from nanosecond time-scale molecular dynamics (MD) simulations of the cyclic peptide ligand eptifibatide in complex with either the fully hydrated normal αIIbβ3 or the D126A mutant in the presence of Ca2+/Mg2+. Calculations were carried out using the OPLS all-atom force-field of the GROMACS simulation package. With respect to normal αIIbβ3, the mutant receptor demonstrated reduced fluctuations in the β3 207–210 and 335 regions and increased fluctuations in the β3 282–284 region. In addition, the ADMIDAS metal ion moved ~3 Å away from the MIDAS and became more solvent exposed. Rearrangements of the coordination of the ADMIDAS involving S123, D251, and D127 were also observed in the D126A mutant compared to normal αIIbβ3. Steered MD simulations were used to investigate the unbinding of eptifibatide from its binding site. The unbinding force for the D126A mutant was similar to that for the normal αIIbβ3. Quantitative estimations of the binding energies of eptifibatide to normal and D126A mutant αIIbβ3 from Molecular Mechanics/Poisson Boltzman Surface Area analysis of the MD trajectories also yielded similar results. Thus, the much greater resistance of D126A mutant cells to detachment from fibrinogen at increasing shear rates does not appear to be explained by differences in fibrinogen-αIIbβ3 interactions at the sites involved in the binding of eptifibatide. Potential alternative mechanisms involve differences in fibrinogen's access to the binding site, interactions with other sites, or changes in fibrinogen avidity due to receptor clustering.

Long-range interactions and parallel scalability in molecular simulations

Typical biomolecular systems such as cellular membranes, DNA, and protein complexes are highly charged. Thus, efficient and accurate treatment of electrostatic interactions is of great importance in computational modeling of such systems. We have employed the GROMACS simulation package to perform extensive benchmarking of different commonly used electrostatic schemes on a range of computer architectures (Pentium-4, IBM Power 4, and Apple/IBM G5) for single processor and parallel performance up to 8 nodes—we have also tested the scalability on four different networks, namely Infiniband, GigaBit Ethernet, Fast Ethernet, and nearly uniform memory architecture, i.e. communication between CPUs is possible by directly reading from or writing to other CPUs' local memory. It turns out that the particle–mesh Ewald method (PME) performs surprisingly well and offers competitive performance unless parallel runs on PC hardware with older network infrastructure are needed. Lipid bilayers of sizes 128, 512 and 2048 lipid molecules were used as the test systems representing typical cases encountered in biomolecular simulations. Our results enable an accurate prediction of computational speed on most current computing systems, both for serial and parallel runs. These results should be helpful in, for example, choosing the most suitable configuration for a small departmental computer cluster.

Molecular dynamics simulation of barnacle cement

Barnacle cement is an underwater adhesive that is used for permanent settlement. Its main components are insoluble protein complexes that have not been fully studied. In present article, we chose two proteins of barnacle cement for study, 36-KD protein and Mrcp-100K protein. In order to investigate the characteristic of above two proteins, we introduced the method of molecular modeling. And the simulation package GROMACS was used to simulate the behavior of these proteins. In this article, before the simulations, we introduce some theories to predict the time scale for polymer relaxation. During the simulation, we mainly focus on two properties of these two proteins: structural stability and adhesive force to substrate. First, we simulate the structural stability of two proteins in water, and then the stability of 36-KD protein in seawater environment is investigated. We find that the stability varies in the different environments. Next, to study adhesive ability of two proteins, we simulate the process of peeling the two proteins from the substrate (graphite). Then, we analyze the main reasons of these results. We find that hydrogen bonds in proteins play an important role in the protein stability. In the process of the peeling, we use Lennard–Jones 12-6 potential to calculate the van der Waals interactions between proteins and substrate.