Mixing of nanofluids: molecular dynamics simulations and modelling

In this paper we study a mixing scheme, which has recently been proposed for microfluids, on the nanoscale. We do this by performing a series of nonequilibrium molecular dynamics simulations. On the nanoscale the chaotic mixing regime is captured. We discover a new phenomenon where the two mixing fluids exchange positions after leaving the mixing intersection. The results from the molecular dynamics simulations also reveal complex spatio-temporal stream velocity profiles generated by the mixing device. We find that these profiles can be modelled through an approximate analytical solution to the Navier–Stokes equation.     

Sensitivity of molecular dynamics simulations to the choice of the X-ray structure used to model an enzymatic reaction

A subject of great practical importance that has not received much attention is the question of the sensitivity of molecular dynamics simulations to the initial X-ray structure used to set up the calculation. We have found two cases in which seemingly similar structures lead to quite different results, and in this article we present a detailed analysis of these cases. The first case is acyl-CoA dehydrogenase, and the chief difference of the two structures is attributed to a slight shift in a backbone carbonyl that causes a key residue (the proton-abstracting base) to be in a bad conformation for reaction. The second case is xylose isomerase, and the chief difference of the two structures appears to be the ligand sphere of a Mg2+ metal cofactor that plays an active role in catalysis.     

Effect of temperature on welding of metallic nanowires investigated using molecular dynamics simulations

The mechanism of welding of Au–Au, Ag–Ag and Au–Ag nanowires (NWs) with head-to-head contact is studied using molecular dynamics simulations based on the second-moment approximation of the many-body tight-binding potential. The effect of temperature in the range of 300–900 K is investigated. Simulation results show that at the initial welding, an incomplete jointing area forms through the interactions of the van der Waals attractive force, and that the jointing area increases with increasing the extent of contact between the two NWs during the welding process and temperature. Few defects form along the (1 1 1) close-packed plane during the welding process because the acting stress is quite low. Among the three NW pairs, the Au–Au NWs have the best cold-welding quality, whereas the Au–Ag NWs have the worst cold-welding quality due to the welding of different materials. With an increase in temperature, the weld stress and the mechanical strength of the NWs significantly decrease, and the number of disordered structures increases. The welding fails when the temperature exceeds the molten temperature of the NWs.     

Trypsinogen activation as observed in accelerated molecular dynamics simulations

Serine proteases are involved in many fundamental physiological processes, and control of their activity mainly results from the fact that they are synthetized in an inactive form that becomes active upon cleavage. Three decades ago Martin Karplus's group performed the first molecular dynamics simulations of trypsin, the most studied member of the serine protease family, to address the transition from the zymogen to its active form. Based on the computational power available at the time, only high frequency fluctuations, but not the transition steps, could be observed. By performing accelerated molecular dynamics (aMD) simulations, an interesting approach that increases the configurational sampling of atomistic simulations, we were able to observe the N‐terminal tail insertion, a crucial step of the transition mechanism. Our results also support the hypothesis that the hydrophobic effect is the main force guiding the insertion step, although substantial enthalpic contributions are important in the activation mechanism. As the N‐terminal tail insertion is a conserved step in the activation of serine proteases, these results afford new perspective on the underlying thermodynamics of the transition from the zymogen to the active enzyme.     

On the behaviour of fluidic material at molecular dynamics boundary conditions used in hybrid molecular-continuum simulations

This paper presents an investigation of the behaviour of the density profile of fluidic material confined by a force field, as it occurs across molecular-continuum mechanics boundaries in multiscale, hybrid molecular-continuum simulations. A theoretical model for the density profile across the boundary is derived. Furthermore, numerical experiments to validate the density profile and thickness of the relaxation zone are performed using molecular dynamics for a Lennard-Jones fluid in gaseous, liquid and supercritical state conditions. The simulation results show excellent agreement with the theoretical model.     

Possible ligand release pathway of dipeptidyl peptidase IV investigated by molecular dynamics simulations

Dipeptidyl peptidase IV (DPP4) is an important target for the treatment of Type II diabetes mellitus. The crystal structure of DPP4 demonstrates that there are two possible pathways to the active site, a side opening and a β propeller opening. However, it still lacks quantitative evidence to illustrate which pathway is more favorable for inhibitor to enter into or release from the active site. In this study, conventional and steered molecular dynamics simulations were performed to explore the details of inhibitor Q448 release from the active site of DPP4 via the two potential pathways. The comparisons of force and work together with potentials of mean force results suggested that the side opening might be more favorable for the inhibitor to pass through. Moreover, Glu205-Glu206 and Phe357 were recognized as two "key residues" in the active site for inhibitor binding. Accordingly, suggestions for further inhibitor design were provided.     

Structure, dynamics, and elasticity of free 16s rRNA helix 44 studied by molecular dynamics simulations

Molecular dynamics (MD) simulations were employed to investigate the structure, dynamics, and local base-pair step deformability of the free 16S ribosomal helix 44 from Thermus thermophilus and of a canonical A-RNA double helix. While helix 44 is bent in the crystal structure of the small ribosomal subunit, the simulated helix 44 is intrinsically straight. It shows, however, substantial instantaneous bends that are isotropic. The spontaneous motions seen in simulations achieve large degrees of bending seen in the X-ray structure and would be entirely sufficient to allow the dynamics of the upper part of helix 44 evidenced by cryo-electron microscopic studies. Analysis of local base-pair step deformability reveals a patch of flexible steps in the upper part of helix 44 and in the area proximal to the bulge bases, suggesting that the upper part of helix 44 has enhanced flexibility. The simulations identify two conformational substates of the second bulge area (bottom part of the helix) with distinct base pairing. In agreement with nuclear magnetic resonance (NMR) and X-ray studies, a flipped out conformational substate of conserved 1492A is seen in the first bulge area. Molecular dynamics (MD) simulations reveal a number of reversible alpha-gamma backbone flips that correspond to transitions between two known A-RNA backbone families. The flipped substates do not cumulate along the trajectory and lead to a modest transient reduction of helical twist with no significant influence on the overall geometry of the duplexes. Despite their considerable flexibility, the simulated structures are very stable with no indication of substantial force field inaccuracies.     

Quantifying uncertainty in molecular dynamics simulations of grain boundary migration

Molecular dynamics simulations of simple bicrystal systems have been much used as a tool to explore how the migration of grain boundaries varies with their structure and with experimental conditions. In order to permit the exploration of a large parameter space, many studies are forced to rely on a small number of simulations (often a single simulation) for each configuration. The motion of a grain boundary is inherently statistical and any variability in the measured grain boundary velocity should be taken into account in subsequent analysis of trends in grain boundary mobility. Here we present the results of large numbers of simulations of equivalent boundaries, which show that this variability can be large, particularly when small systems are simulated. We show how a bootstrap resampling approach can be used to characterise the statistical uncertainty in boundary velocity using the information present in a single simulation. We show that the approach is robust across a variety of system sizes, temperatures and driving force strengths and types, and provides a good order-of-magnitude measure of the population standard deviation across multiple equivalent simulations.     

Reconstruction of the src-SH3 protein domain transition state ensemble using multiscale molecular dynamics simulations

We use an integrated computational approach to reconstruct accurately the transition state ensemble (TSE) for folding of the src-SH3 protein domain. We first identify putative TSE conformations from free energy surfaces generated by importance sampling molecular dynamics for a fully atomic, solvated model of the src-SH3 protein domain. These putative TSE conformations are then subjected to a folding analysis using a coarse-grained representation of the protein and rapid discrete molecular dynamics simulations. Those conformations that fold to the native conformation with a probability (P(fold)) of approximately 0.5, constitute the true transition state. Approximately 20% of the putative TSE structures were found to have a P(fold) near 0.5, indicating that, although correct TSE conformations are populated at the free energy barrier, there is a critical need to refine this ensemble. Our simulations indicate that the true TSE conformations are compact, with a well-defined central beta sheet, in good agreement with previous experimental and theoretical studies. A structured central beta sheet was found to be present in a number of pre-TSE conformations, however, indicating that this element, although required in the transition state, does not define it uniquely. An additional tight cluster of contacts between highly conserved residues belonging to the diverging turn and second beta-sheet of the protein emerged as being critical elements of the folding nucleus. A number of commonly used order parameters to identify the transition state for folding were investigated, with the number of native Cbeta contacts displaying the most satisfactory correlation with P(fold) values.     

Investigations of the thermostability of rubredoxin models using molecular dynamics simulations.

The affects of differences in amino acid sequence on the temperature stability of the three-dimensional structure of the small beta-sheet protein, rubredoxin (Rd), was revealed when a set of homology models was subjected to molecular dynamics simulations at relatively high temperatures. Models of Rd from the hyperthermophile, Pyrococcus furiosus (Pf), an organism that grows optimally at 100 degrees C, were compared to three mesophilic Rds of known X-ray crystal structure. Simulations covering the limits of known Rd thermostabilities were carried out at temperatures of 300 K, 343 K, 373 K, and 413 K. They suggest that Rd stability is correlated with structural dynamics. Because the dynamic behavior of three Pf Rd models was consistently different from the dynamic behavior of the three mesophilic Rd structures, detailed analysis of the temperature-dependent dynamic behavior was carried out. The major differences between the models of the protein from the hyperthermophile and the others were: (1) an obvious temperature-dependent transition in the mesophilic structures not seen with the Pf Rd models, (2) consistent AMBER energy for the Pf Rd due to differences in nonbonded interaction terms, (3) less variation in the average conformations for the Pf Rd models with temperature, and (4) the presence of more extensive secondary structure for the Pf Rd models. These unsolvated dynamics simulations support a simple, general hypothesis to explain the hyperthermostability of Pf Rd. Its structure simplifies the conformational space to give a single minimum accessible over an extreme range of temperatures, whereas the mesophilic proteins sample a more complex conformational space with two or more minima over the same temperature range.     

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.     

Sequence dependence of amyloid fibril formation: insights from molecular dynamics simulations

The clarification of the physico-chemical determinants underlying amyloid deposition is critical for our understanding of misfolding diseases. With this purpose we have performed a systematic all-atom molecular dynamics (MD) study of a series of single point mutants of the de novo designed amyloidogenic peptide STVIIE. Sixteen different 50ns long simulations using explicit solvent have been carried out starting from four different conformations of a polymeric six-stranded beta-sheet. The simulations have provided evidence for the influence of a small number of site-specific hydrophobic interactions on the packing and stabilization of nascent aggregates, as well as the interplay between side-chain interactions and the net charge of the molecule on the strand arrangement of polymeric beta-sheets. This MD analysis has also shed light into the origin of the position dependence on mutation of beta-sheet polymerization that was found experimentally for this model system. Our results suggest that MD can be applied to detect critical positions for beta-sheet aggregation within a given amyloidogenic stretch. Studies similar to the one presented here can guide site-directed mutations or the design of drugs that specifically disrupt the key stabilizing interactions of beta-sheet aggregates.     

Improved model of hydrated calcium ion for molecular dynamics simulations using classical biomolecular force fields

Calcium ions (Ca²⁺) play key roles in various fundamental biological processes such as cell signaling and brain function. Molecular dynamics (MD) simulations have been used to study such interactions, however, the accuracy of the Ca²⁺ models provided by the standard MD force fields has not been rigorously tested. Here, we assess the performance of the Ca²⁺ models from the most popular classical force fields AMBER and CHARMM by computing the osmotic pressure of model compounds and the free energy of DNA–DNA interactions. In the simulations performed using the two standard models, Ca²⁺ ions are seen to form artificial clusters with chloride, acetate, and phosphate species; the osmotic pressure of CaAc₂ and CaCl₂ solutions is a small fraction of the experimental values for both force fields. Using the standard parameterization of Ca²⁺ ions in the simulations of Ca²⁺‐mediated DNA–DNA interactions leads to qualitatively wrong outcomes: both AMBER and CHARMM simulations suggest strong inter‐DNA attraction whereas, in experiment, DNA molecules repel one another. The artificial attraction of Ca²⁺ to DNA phosphate is strong enough to affect the direction of the electric field‐driven translocation of DNA through a solid‐state nanopore. To address these shortcomings of the standard Ca²⁺ model, we introduce a custom model of a hydrated Ca²⁺ ion and show that using our model brings the results of the above MD simulations in quantitative agreement with experiment. Our improved model of Ca²⁺ can be readily applied to MD simulations of various biomolecular systems, including nucleic acids, proteins and lipid bilayer membranes. © 2016 Wiley Periodicals, Inc. Biopolymers 105: 752–763, 2016.     

The role of disulfide bonds and N‐terminus in the structural properties of hepcidins: Insights from molecular dynamics simulations

The main purpose of this work is to analyse, by means of molecular dynamics (MD) simulations both structural and mechanical‐dynamical differences between Hepcidin‐20 and Hepcidin‐25 in both oxidized and reduced states in aqueous solution. Results indicate that the presence of disulfide bonds is essential, in both peptides, for maintaining their β‐hairpin motif. As a matter of fact, the lack of this intra‐peptide covalent interactions produces an almost immediate deviation from the oxidized, plausibly active, structure in both the systems. Interestingly, reduced Hepcidin‐25 turns out to be characterized by a highly fluctuating structure which is found to rapidly span a large number of configurations at equilibrium. On the other hand, loss of disulfide bonds in the shorter peptide, results in a more compact and relatively rigid double‐turn structure. Comparison of mechanical–dynamical properties and sidechains–sidechains interactions in oxidized Hepcidin‐20 and Hepcidin‐25 strongly suggest also the key role of N‐terminus in the aggregation tendency of Hepcidin‐25. © 2010 Wiley Periodicals, Inc. Biopolymers 93: 917–926, 2010.     

Effect O6‐guanine alkylation on DNA flexibility studied by comparative molecular dynamics simulations

Alkylation of guanine at the O6 atom is a highly mutagenic DNA lesion because it alters the coding specificity of the base causing G:C to A:T transversion mutations. Specific DNA repair enzymes, e.g. O⁶‐alkylguanin‐DNA‐Transferases (AGT), recognize and repair such damage after looping out the damaged base to transfer it into the enzyme active site. The exact mechanism how the repair enzyme identifies a damaged site within a large surplus of undamaged DNA is not fully understood. The O⁶‐alkylation of guanine may change the deformability of DNA which may facilitate the initial binding of a repair enzyme at the damaged site. In order to characterize the effect of O⁶‐methyl‐guanine (O⁶‐MeG) containing base pairs on the DNA deformability extensive comparative molecular dynamics (MD) simulations on duplex DNA with central G:C, O⁶‐MeG:C or O⁶‐MeG:T base pairs were performed. The simulations indicate significant differences in the helical deformability due to the presence of O⁶‐MeG compared to regular undamaged DNA. This includes enhanced base pair opening, shear and stagger motions and alterations in the backbone fine structure caused in part by transient rupture of the base pairing at the damaged site and transient insertion of water molecules. It is likely that the increased opening motions of O⁶‐MeG:C or O⁶‐MeG:T base pairs play a decisive role for the induced fit recognition or for the looping out of the damaged base by repair enzymes. © 2014 Wiley Periodicals, Inc. Biopolymers 103: 23–32, 2015.     

Immunogenicity of peptide-vaccine candidates predicted by molecular dynamics simulations

We present an in silico, structure-based approach for design and evaluation of conformationally restricted peptide-vaccines. In particular, we designed four cyclic peptides of ten or 11 residues mimicking the crystallographically observed beta-turn conformation of a predicted immunodominant loop of PorA from Neisseria meningitidis. Conformational correctness and stability of the peptide designs, as evaluated by molecular dynamics simulations, correctly predicted the immunogenicity of the peptides. We observed a peptide-induced functional antibody response that, remarkably, exceeded the response induced by the native protein in outer membrane vesicles, without losing specificity for related strains. The presented approach offers tools for a priori design and selection of peptide-vaccine candidates with full biological activity. This approach could be widely applicable: to outer membrane proteins of Gram-negative bacteria, and to other epitopes in a large range of pathogens.     

Convergence of molecular dynamics simulations of membrane proteins

The central question in evaluating almost any result from a molecular dynamics simulation is whether the calculation has converged. Unfortunately, assessing the ergodicity of a single trajectory is very difficult to do. In this work, we assess the sampling of molecular dynamics simulations of the membrane protein rhodopsin by comparing the results from 26 independent trajectories, each run for 100 ns. By examining principal components and cluster populations, we show that rhodopsin's fluctuations are not well described by 100 ns of dynamics, and that the sampling is not fully converged even for individual loops. The results serve as a reminder of the caution required when interpreting molecular dynamics simulations of macromolecules.     

Ligand unbinding pathways from the vitamin D receptor studied by molecular dynamics simulations

Molecular dynamics simulation techniques have been used to study the unbinding pathways of 1α,25-dihydroxyvitamin D₃ from the ligand-binding pocket of the vitamin D receptor (VDR). The pathways observed in a large number of relatively short (<200 ps) random acceleration molecular dynamics (RAMD) trajectories were found to be in fair agreement, both in terms of pathway locations and deduced relative preferences, compared to targeted molecular dynamics (TMD) and streered molecular dynamics simulations (SMD). However, the high-velocity ligand expulsions of RAMD tend to favor straight expulsion trajectories and the observed relative frequencies of different pathways were biased towards the probability of entering a particular exit channel. Simulations indicated that for VDR the unbinding pathway between the H1–H2 loop and the β-sheet between H5 and H6 is more favorable than the pathway located between the H1–H2 loop and H3. The latter pathway has been suggested to be the most likely unbinding path for thyroid hormone receptors (TRs) and a likely path for retinoic acid receptor. Ligand entry/exit through these two pathways would not require displacement of H12 from its agonistic position. Differences in the packing of the H1, H2, H3 and β-sheet region explain the changed relative preference of the two unbinding pathways in VDR and TRs. Based on the crystal structures of the ligand binding domains of class 2 nuclear receptors, whose members are VDR and TRs, this receptor class can be divided in two groups according to the packing of the H1, H2, H3 and β-sheet region. Electronic supplementary material  The online version of this article (doi:) contains supplementary material, which is available to authorized users.