Explicit-solvent molecular dynamics simulations of a DNA tetradecanucleotide duplex: lattice-sum versus reaction-field electrostatics

Five long-timescale (10 ns) explicit-solvent molecular dynamics simulations of a DNA tetradecanucleotide dimer are performed using the GROMOS 45A4 force field and the simple-point-charge water model, in order to investigate the effect of the treatment of long-range electrostatic interactions as well as of the box shape and size on the structure and dynamics of the molecule (starting from an idealised B-DNA conformation). Long-range electrostatic interactions are handled using either a lattice-sum (LS) method (particle–particle–particle–mesh; one simulation performed within a cubic box) or a cutoff-based reaction-field (RF) method (four simulations, with long-range cutoff distances of 1.4 or 2.0 nm and performed within cubic or truncated octahedral periodic boxes). The overall double-helical structure, including Watson–Crick (WC) base-pairing, is well conserved in the simulation employing the LS scheme. In contrast, the WC base-pairing is nearly completely disrupted in the four simulations employing the RF scheme. These four simulations result in highly distorted compact (cutoff distance of 1.4 nm) or extended (cutoff distance of 2 nm) structures, irrespective of the shape and size of the computational box. These differences observed between the two schemes seem correlated with large differences in the radial distribution function between charged entities (backbone phosphate groups and sodium counterions) within the system.     

Nanoscale, Electric Field-Driven Water Bridges in Vacuum Gaps and Lipid Bilayers

Formation of a water bridge across the lipid bilayer is the first stage of pore formation in molecular dynamic (MD) simulations of electroporation, suggesting that the intrusion of individual water molecules into the membrane interior is the initiation event in a sequence that leads to the formation of a conductive membrane pore. To delineate more clearly the role of water in membrane permeabilization, we conducted extensive MD simulations of water bridge formation, stabilization, and collapse in palmitoyloleoylphosphatidylcholine bilayers and in water–vacuum–water systems, in which two groups of water molecules are separated by a 2.8 nm vacuum gap, a simple analog of a phospholipid bilayer. Certain features, such as the exponential decrease in water bridge initiation time with increased external electric field, are similar in both systems. Other features, such as the relationship between water bridge lifetime and the diameter of the water bridge, are quite different between the two systems. Data such as these contribute to a better and more quantitative understanding of the relative roles of water and lipid in membrane electropore creation and annihilation, facilitating a mechanism-driven development of electroporation protocols. These methods can be extended to more complex, heterogeneous systems that include membrane proteins and intracellular and extracellular membrane attachments, leading to more accurate models of living cells in electric fields.     

Non-equilibrium molecular dynamics study of nanoscale thermal contact resistance

Interfaces play an important role in microscale and nanoscale heat transfer processes with molecular dynamics (MD) simulations often used to study these interfacial phenomena. In this study, two models were used to simulate thermal conduction across micro contact points and the thermal contact resistance using non-equilibrium molecular dynamics simulations with consideration of the near field radiation. When the ratio of the length of the micro contact to the length of the conduction region is less than 0.125, the influence of the near field radiation should be considered; but when the ratio is larger than 0.2, it can be neglected. When the computational domain sizes are 8.50 × 10.62 × 8.50 nm and 10.62 × 10.62 × 10.62 nm, the MD results show that the thermal contact resistance exponentially increases with decreasing area of the micro contact point and increases with increasing micro contact layer thickness. The MD thermal contact resistances in nanoscale are much larger than that of the classical thermal analysis since the material thermal conductivity reduction is ignored in the classical model. The results also show that material defects increase the thermal resistance.     

Force-field dependence on the liquid-expanded to liquid-condensed transition in DPPC monolayers

We have investigated using molecular dynamics simulations, the influence of the interaction cut-off and water model on the surface-pressure area isotherms of dipalmitoylphosphatidylcholine monolayers, where the phospholipids and the water molecules are modelled atomistically. We find that both the cut-off employed and the water model, influence the pressure area isotherms and the location of the liquid-expanded to liquid-condensed transitions. The combination of the Berger's et al. force field, with the TIP4P/2005 water model, and a long cut-off for the pair interactions ( ≥ 1.7 nm) provides a more accurate prediction of the surface pressure–area isotherm and reproduces the liquid condensed–liquid expanded transition observed in the experiments at 310 K.     

Hydration effects on the electrostatic potential around tuftsin

The electrostatic potential and component dielectric constants from molecular dynamics (MD) trajectories of tuftsin, a tetrapeptide with the amino acid sequence Thr–Lys–Pro–Arg in water and in saline solution are presented. The results obtained from the analysis of the MD trajectories for the total electrostatic potential at points on a grid using the Ewald technique are compared with the solution to the Poisson–Boltzmann (PB) equation. The latter was solved using several sets of dielectric constant parameters. The effects of structural averaging on the PB results were also considered. Solute conformational mobility in simulations gives rise to an electrostatic potential map around the solute dominated by the solute monopole (or lowest order multipole). The detailed spatial variation of the electrostatic potential on the molecular surface brought about by the compounded effects of the distribution of water and ions close to the peptide, solvent mobility, and solute conformational mobility are not qualitatively reproducible from a reparametrization of the input solute and solvent dielectric constants to the PB equation for a single structure or for structurally averaged PB calculations. Nevertheless, by fitting the PB to the MD electrostatic potential surfaces with the dielectric constants as fitting parameters, we found that the values that give the best fit are the values calculated from the MD trajectories. Implications of using such field calculations on the design of tuftsin peptide analogues are discussed. © 1999 John Wiley & Sons, Inc. Biopoly 50: 133–143, 1999     

Structure,stability and water permeation of ([D-Leu-L-Lys-(D-Gln-L-Ala)3]) cyclic peptide nanotube: a molecular dynamics study

The structural stability of 8 × ([D-Leu-L-Lys-(D-Gln-L-Ala)₃]) cyclic peptide nanotube (CPN) in water and different phospholipid bilayers were explored by 100 ns independent molecular dynamics (MD) simulations. The role of non-bonded interaction energy between the side and main chains of cyclic peptide rings in different membrane environments assessed, wherein the repulsive electrostatic interaction energy between neighbouring cyclic peptide rings was found adequate to break hydrogen bond energy thereby to crumple CPN. Further, the water permeation across the CPN channel was studied in four types of phospholipid bilayers- DMPG (1,2-Dimyristoyl-sn-glycero-3-phosphorylglycerol), DMPS (1,2-Dimyristoyl-sn-glycero-3-phosphoserine), POPC (1-Palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) and POPE (1-Palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine) from MD simulations. DMPS membrane shows higher non-bonded interaction energies (?1913.06 kJ/mol of electrostatic interaction energy and ?994.13 kJ/mol of van der Waals interaction energy) with CPN due to the presence of polar molecules in lipid structure. Thusly, the non-bonded interaction energies were essential towards the stability of CPN than hydrogen bonds between the nearby cyclic peptides. The result also reveals the role of side chains, hydrogen bonds and non-bonded interaction energies in an aqueous environment. The diffusion coefficient of water obtained from means square deviation calculation shows similar coefficients irrespective of the lipid surroundings. However, the permeation coefficients demonstrate water flow in the channel relies upon the environment.     

Analysis of dissociation process for gas hydrates by molecular dynamics simulation

The dissociation processes of methane and carbon dioxide hydrates were investigated by molecular dynamics simulation. The simulations were performed with 368 water molecules and 64 gas molecules using NPT ensembles. The TraPPE (single-site) and 5-site models were adopted for methane molecules. The EPM2 (3-site) and SPC/E models were used for carbon dioxide and water molecules, respectively. The simulations were carried out at 270 K and 5.0 MPa for hydrate stabilisation. Then, temperature was increased up to 370 K. The temperature increasing rates were 0.1–20 TK/s. The gas hydrates dissociated during increasing temperature or at 370 K. The potential models of methane molecule did not much influence the dissociation process of methane hydrate. The mechanisms of dissociation process were analysed with the coordination numbers and mean square displacements. It was found that the water cages break down first, then the gas molecules escape from the water cages. The methane hydrate was more stable than the carbon dioxide hydrate at the calculated conditions.     

Cation transport in nanopores

Using periodic boundary conditions and external electric potential field, we have simulated an ion current flow through a flexible nanopore using cations and an explicit extended simple point charge (SPC/E) water with molecular dynamics simulation. The simulation voltages range goes beyond the usual ionic channel measurements ( ± 1 V) and yields useful information about density profiles, current density distribution and current–voltage (I–V) characteristics.     

Curcumin specifically binds to the human calcium–calmodulin-dependent protein kinase IV: fluorescence and molecular dynamics simulation studies

Calcium–calmodulin-dependent protein kinase IV (CAMK4) plays significant role in the regulation of calcium-dependent gene expression, and thus, it is involved in varieties of cellular functions such as cell signaling and neuronal survival. On the other hand, curcumin, a naturally occurring yellow bioactive component of turmeric possesses wide spectrum of biological actions, and it is widely used to treat atherosclerosis, diabetes, cancer, and inflammation. It also acts as an antioxidant. Here, we studied the interaction of curcumin with human CAMK4 at pH 7.4 using molecular docking, molecular dynamics (MD) simulations, fluorescence binding, and surface plasmon resonance (SPR) methods. We performed MD simulations for both neutral and anionic forms of CAMK4-curcumin complexes for a reasonably long time (150 ns) to see the overall stability of the protein–ligand complex. Molecular docking studies revealed that the curcumin binds in the large hydrophobic cavity of kinase domain of CAMK4 through several hydrophobic and hydrogen-bonded interactions. Additionally, MD simulations studies contributed in understanding the stability of protein–ligand complex system in aqueous solution and conformational changes in the CAMK4 upon binding of curcumin. A significant increase in the fluorescence intensity at 495 nm was observed (λ_exc = 425 nm), suggesting a strong interaction of curcumin to the CAMK4. A high binding affinity (K_D = 3.7 × 10^?8 ± .03 M) of curcumin for the CAMK4 was measured by SPR further indicating curcumin as a potential ligand for the CAMK4. This study will provide insights into designing a new inspired curcumin derivatives as therapeutic agents against many life-threatening diseases.     

Analysis of crystal growth of methane hydrate using molecular dynamics simulation

Methane hydrate is a crystalline compound with methane molecules enclosed in cages formed by hydrogen-bonded water molecules. Understanding the mechanism of nucleation and crystal growth from methane vapour and liquid water is important for all hydrate applications. However, processes near the water/methane interface are still unclear. In this work, we focused on the crystal growth of methane hydrate seeds located near the water/methane interface. We performed molecular dynamics (MD) simulation and analysed the crystal growth of the hydrate seed at the interface. New cages formed in the liquid water phase were stabilised when they shared faces with the hydrate seed. We also investigated the crystal growth rate as the time development of the number of methane molecules trapped in hydrate cages, based on the trajectory of the MD simulation. The calculated growth rate in the direction that covers the interface was 1.38 times that in the direction towards the inside of the water phase.     

Hybrid MD-Nernst Planck model of α-hemolysin conductance properties

Motivated by experiments in which an applied electric field translocates polynucleotides through an α-hemolysin protein channel causing ionic current transient blockade, a hybrid simulation model is proposed to predict the conductance properties of the open channel. Time scales corresponding to ion permeation processes are reached using the Poisson–Nernst–Planck (PNP) electro-diffusion model in which both solvent and local ion concentrations are represented as a continuum. The diffusion coefficients of the ions (K⁺ and Cl^?) input in the PNP model are, however, calculated from all-atom molecular dynamics (MD). In the MD simulations, a reduced representation of the channel is used. The channel is solvated in a 1?M KCl solution, and an external electric field is applied. The pore specific diffusion coefficients for both ionic species are reduced 5–7 times in comparison to bulk values. Significant statistical variations (17–45%) of the pore-ions diffusivities are observed. Within the statistics, the ionic diffusivities remain invariable for a range of external applied voltages between 30 and 240?mV. In the 2D-PNP calculations, the pore stem is approximated by a smooth cylinder of radius ～9?Å with two constriction blocks where the radius is reduced to ～6?Å. The electrostatic potential includes the contribution from the atomistic charges. The MD-PNP model shows that the atomic charges are responsible for the rectifying behaviour and for the slight anion selectivity of the α-hemolysin pore. Independent of the hierarchy between the anion and cation diffusivities, the anionic contribution to the total ionic current will dominate. The predictions of the MD-PNP model are in good agreement with experimental data and give confidence in the present approach of bridging time scales by combining a microscopic and macroscopic model.     

Molecular dynamics simulations of palmitoyloleoylphosphatidylglycerol bilayers

Phosphatidylglycerol (PG) is one of the important components of biological membranes, but there is a paucity of experimental data to test the accuracy of molecular dynamics (MD) simulations. This work consists of testing the accuracy of the CHARMM36 (C36) lipid force field on 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoglycerol (POPG) lipid bilayers. MD simulations of POPG lipid bilayers are compared to recently available X-ray and neutron scattering and deuterium NMR measurements. Overall, the C36 lipid force field accurately represents the X-ray and neutron form factors, bilayer and hydrocarbon thicknesses and chain deuterium order parameters. The surface area per lipid from MD simulations with C36 (67.7 ± 0.2 Å²) is in excellent agreement with the experimentally determined value of 66.0 ± 1.3 Å². C36 outperforms the lipid force field developed by Berger et al. [15] and suggests that past studies with this force field may result in lateral areas that are too small. Moreover, our studies give some insight into the structural model used in experiments and suggest that the functional form for the head group may not be Gaussian-like. Based on our simulations, the POPG lipid in the C36 lipid force field is well parameterised and can be used for other PG lipids and membrane models with mixed lipids.     

Effect of membrane tension on the electric field and dipole potential of lipid bilayer membrane

The dipole potential of lipid bilayer membrane controls the difference in permeability of the membrane to oppositely charged ions. We have combined molecular dynamics (MD) simulations and experimental studies to determine changes in electric field and electrostatic potential of 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) lipid bilayer in response to applied membrane tension. MD simulations based on CHARMM36 force field showed that electrostatic potential of DOPC bilayer decreases by ~45mV in the physiologically relevant range of membrane tension values (0 to 15dyn/cm). The electrostatic field exhibits a peak (~0.8×10(9)V/m) near the water/lipid interface which shifts by 0.9? towards the bilayer center at 15dyn/cm. Maximum membrane tension of 15dyn/cm caused 6.4% increase in area per lipid, 4.7% decrease in bilayer thickness and 1.4% increase in the volume of the bilayer. Dipole-potential sensitive fluorescent probes were used to detect membrane tension induced changes in DOPC vesicles exposed to osmotic stress. Experiments confirmed that dipole potential of DOPC bilayer decreases at higher membrane tensions. These results are suggestive of a potentially new mechanosensing mechanism by which mechanically induced structural changes in the lipid bilayer membrane could modulate the function of membrane proteins by altering electrostatic interactions and energetics of protein conformational states.     

Molecular dynamics and quantum chemical studies of solvent effects on cyclo glycylglycine and glycylalanine dipeptides

Hydrogen bond (H-bond) interactions between the two cyclo dipeptides, cyclo(glycyl-glycine) (CGG) and cyclo(glycyl-alanine) (CGA), and water have been studied using molecular dynamics (MD) and quantum chemical methods. The MD studies have been carried out on CGG and CGA in water using fixed charge force field (AMBER ff03) for over 10 ns with a MD time step of 2 fs. The results of this study show that the solvation pattern influences the conformations of the cyclo dipeptides. Following molecular simulations, post Hartree–Fock and density functional theory methods have been used to explore the molecular properties of the cyclo dipeptides in gaseous and aqueous phase environments. The self-consistent reaction field theory has been used to optimise the cyclopeptides in diethyl ether (? = 4.3) and water (? = 78.5), and the solvent effects have been analysed. A cluster of eight water molecules leads to the formation of first solvation shell of CGG and CGA and the strong H-bonding mainly contributes to the interaction energies. The H-bond interactions have been analysed by the calculation of electron density ρ(r) and its Laplacian ▽²ρ(r) at bond critical points using atoms in molecules theory. The natural bond orbital analysis was carried out to reveal the nature of H-bond interactions. In the solvated complexes, the keto carbons registered the maximum NMR chemical shifts.     

Binding affinity of hydroxamate inhibitors of matrix metalloproteinase-2

Here we report molecular dynamics (MD) and free energy perturbation (FEP) simulations applied to hydroxamate-matrix metalloproteinase-2 (MMP-2) complex systems. We have developed some new force field parameters for the hydroxamate functional group that were not included in the AMBER94 force field but were necessary in our simulations. For the representation of the active zinc center, a bonded model was adopted in which restrained electrostatic potential fitting (RESP) charges were used as the electrostatic representation of this model. Using the resulted bonded model, FEP simulations predict the relative binding free energy in good agreement with the experimental value. By analyzing the molecular dynamics (MD) trajectories of the two complex systems, we can provide an explanation of why one of the two inhibitors is favored over the other. The results provide a chemical insight into the interactions between inhibitor and enzyme, and can indicate changes in the inhibitor that would enhance inhibitor–enzyme interactions.Figure The scheme of the binding site     

Fate of methane bubbles released by pockmarks in Lake Constance

In the eastern part of Lake Constance, the second largest pre-alpine lake in Europe, about five hundred pockmarks (morphological depressions on the lake floor) were recently discovered of which ~40 % release methane bubbles. The carbon isotopic composition of the escaping gas indicated that the methane is of biogenic origin. In our study, we investigated the fate of the released methane bubbles, i.e., the dissolution, oxidation or transport of the bubbles to the surface. At a littoral pockmark site (PM12, 12 m water depth) and a profundal pockmark (PM80, 80 m water depth), we analysed the dissolved methane concentrations and the methane isotopic carbon signature in the water column. At PM80, higher methane concentrations (up to 1,523 nM), compared to the control site and the surface waters (225 ± 72 nM), were recorded only on some occasions and only in the bottom water, despite the fact that the released bubbles were dissolving within the hypolimnion based on bubble modeling. The isotope data suggest that most of the dissolved methane is oxidized below 40 m water depth. The isotopic signature of the methane in the surface water at PM80, however, differed from that of the methane in the hypolimnion; therefore, the surface methane at this profundal site is most likely an export product from the littoral zone. Assuming an initial bubble diameter of 5 mm, we calculated that these small bubbles would reach the surface, but approximately 96 % of the methane would have dissolved from the bubble into the hypolimnion. At PM12, we observed higher concentrations of dissolved methane (312 ± 52 nM) with no significant differences between seasons or between control sites versus pockmark site. In the shallow water, divers estimated the bubble size to be 10–15 mm, which from a release depth of 12 m would barely dissolved into the water column. The isotopic signature also indicated that there had been almost no methane oxidation in the shallow water column. Thus, the water depth of bubble release as well as the initial bubble size determine whether the methane enters the atmosphere largely unhindered (shallow site) or if the released methane is incorporated into the profundal water column.     

Importance of polarization effect in the study of metalloproteins: Application of polarized protein specific charge scheme in predicting the reduction potential of azurin

Molecular dynamics (MD) simulation is commonly used in the study of protein dynamics, and in recent years, the extension of MD simulation to the study of metalloproteins is gaining much interest. Choice of force field is crucial in MD studies, and the inclusion of metal centers complicates the process of accurately describing the electrostatic environment that surrounds the redox centre. Herein, we would like to explore the importance of including electrostatic contribution from both protein and solvent in the study of metalloproteins. MD simulations with the implementation of thermodynamic integration will be conducted to model the reduction process of azurin from Pseudomonas aeruginosa. Three charge schemes will be used to derive the partial charges of azurin. These charge schemes differ in terms of the amount of immediate environment, respective to copper, considered during charge fitting, which ranges from the inclusion of copper and residues in the first coordination sphere during density functional theory charge fitting to the comprehensive inclusion of protein and solvent effect surrounding the metal centre using polarized protein‐specific charge scheme. From the simulations conducted, the relative reduction potential of the mutated azurins respective to that of wild‐type azurin (ΔE_cal) were calculated and compared with experimental values. The ΔE_cal approached experimental value with increasing consideration of environmental effect hence substantiating the importance of polarization effect in the study of metalloproteins. This study also attests the practicality of polarized protein‐specific charge as a computational tool capable of incorporating both protein environment and solvent effect into MD simulations. Proteins 2014; 82:2209–2219. © 2014 Wiley Periodicals, Inc.     

Affinity enhancement of nanobody binding to EGFR: in silico site-directed mutagenesis and molecular dynamics simulation approaches

Epidermal growth factor receptor (EGFR), a transmembrane glycoprotein, is overexpressed in many cancers such as head-neck, breast, prostate, and skin cancers for this reason it is a good target in cancer therapy and diagnosis. In nanobody-based cancer diagnosis and treatment, nanobodies with high affinity toward receptor (e.g. EGFR) results in effective treatment or diagnosis of cancer. In this regard, the main aim of this study is to develop a method based on molecular dynamic (MD) simulations for designing of 7D12 based nanobody with high affinity compared with wild-type nanobody. By surveying electrostatic and desolvation interactions between different residues of 7D12 and EGFR, the critical residues of 7D12 that play the main role in the binding of 7D12 to EGFR were elucidated and based on these residues, five logical variants were designed. Following the 50 ns MD simulations, pull and umbrella sampling simulation were performed for 7D12 and all its variants in complex with EGFR. Binding free energy of 7D12 (and all its variants) with EGFR was obtained by weighted histogram analysis method. According to binding free energy results, GLY101 to GLU mutation showed the highest binding affinity but this variant is unstable after 50 ns MD simulations. ALA100 to GLU mutation shows suitable binding enhancement with acceptable structural stability. Suitable position and orientation of GLU in residue 100 of 7D12 against related amino acids of EGFR formed some extra hydrogen and electrostatic interactions which resulted in binding enhancement.