Umbrella Sampling Method Research Articles

High-precision work distributions for extreme nonequilibrium processes in large systems

The distributions of work for strongly nonequilibrium processes are studied using a very general form of a large-deviation approach, which allows one to study distributions down to extremely small probabilities of almost arbitrary quantities of interest for equilibrium, nonequilibrium stationary, and even nonstationary processes. The method is applied to quickly vary the external field in a wide range B = 3 ↔ 0 for a critical (T = 2.269) two-dimensional Ising system of size L × L = 128 × 128. To obtain free-energy differences from the work distributions, they must be studied in ranges where the probabilities are as small as 10^{-240}, which is not possible using direct simulation approaches. By comparison with the exact free energies, which are available for this model for the zero-field case, one sees that the present approach allows one to obtain the free energy with a very high relative precision of 10^{-4}. This works well also for a nonzero field, i.e., for a case where standard umbrella-sampling methods are not efficient to calculate free energies. Furthermore, for the present case it is verified that the resulting distributions of work for forward and backward processes fulfill Crooks theorem with high precision. Finally, the free energy for the Ising magnet as a function of the field strength is obtained.

Molecular dynamics simulations and free energy profile of Paracetamol in DPPC and DMPC lipid bilayers

Molecular dynamics (MD) simulations and biased MD simulation were carried out for the neutral form of Paracetamol inserted in fully hydrated dipalmitoylphosphatidylcholine (DPPC) and dimyristoylphosphatidylcholine (DMPC) lipid bilayers. For comparison, fully hydrated DMPC and DPPC lipid bilayers were also simulated separately without Paracetamol. The simulation time for each system was 50 ns. At two concentrations of Paracetamol, various properties of the lipid bilayer such as area per lipid, order parameter, diffusion coefficient, radial distribution function, electrostatic potential, mass density and hydrogen bonds have been calculated. Also, the convergence in time of the free energy profile of the Paracetamol along a DPPC bilayer normal was calculated by umbrella sampling method. From the obtained results, it can be concluded that neutral form of Paracetamol shows a generally similar behaviour in DPPC and DMPC lipid bilayers. It was shown that the addition of Paracetamol causes a decrease in tail order parameter of both DPPC and DMPC lipid bilayers and the tail of Paracetamol adopts an inward orientation in the lipid bilayers. Also from the free energy profile, the high penetration barrier in the bilayer centre was determined. The behaviour of Paracetamol in different dosages in two types of lipid bilayer membranes, DPPC and DMPC has been investigated by MD simulation. For the neutral form of the drug, the best position in the membrane has been identified by free energy calculations.

Interactions of aqueous amino acids and proteins with the (110) surface of ZnS in molecular dynamics simulations

The growing usage of nanoparticles of zinc sulfide as quantum dots and biosensors calls for a theoretical assessment of interactions of ZnS with biomolecules. We employ the molecular-dynamics-based umbrella sampling method to determine potentials of mean force for 20 single amino acids near the ZnS (110) surface in aqueous solutions. We find that five amino acids do not bind at all and the binding energy of the remaining amino acids does not exceed 4.3 kJ/mol. Such energies are comparable to those found for ZnO (and to hydrogen bonds in proteins) but the nature of the specificity is different. Cysteine can bind with ZnS in a covalent way, e.g., by forming the disulfide bond with S in the solid. If this effect is included within a model incorporating the Morse potential, then the potential well becomes much deeper--the binding energy is close to 98 kJ/mol. We then consider tryptophan cage, a protein of 20 residues, and characterize its events of adsorption to ZnS. We demonstrate the relevance of interactions between the amino acids in the selection of optimal adsorbed conformations and recognize the key role of cysteine in generation of lasting adsorption. We show that ZnS is more hydrophobic than ZnO and that the density profile of water is quite different than that forming near ZnO--it has only a minor articulation into layers. Furthermore, the first layer of water is disordered and mobile.

Selectivity trend of gas separation through nanoporous graphene

By means of molecular dynamics (MD) simulations, we demonstrate that porous graphene can efficiently separate gases according to their molecular sizes. The flux sequence from the classical MD simulation is H2>CO2≫N2>Ar>CH4, which generally follows the trend in the kinetic diameters. This trend is also confirmed from the fluxes based on the computed free energy barriers for gas permeation using the umbrella sampling method and kinetic theory of gases. Both brute-force MD simulations and free-energy calcualtions lead to the flux trend consistent with experiments. Case studies of two compositions of CO2/N2 mixtures further demonstrate the separation capability of nanoporous graphene.

A practical perspective on the implementation of hyperdynamics for accelerated simulation.

Consideration is given to several practical issues arising during the implementation of hyperdynamics, a methodology that extends the time scale of the conventional molecular dynamics simulation potentially by orders of magnitude. First, the methodology is reformulated in terms of the transition rate based on the buffer region approach (buffer rate), which can describe transitions in more general contexts than the transition state theory (TST). It will be shown that hyperdynamics can exactly preserve the buffer rate as well as the TST rate, which broadens the scope of the method. Next, the originally proposed scheme to compute the boost factor on-the-fly is reviewed and some alternative methods, one of which uses the umbrella sampling method, are presented. Finally, the methodology is validated in the context of a 1-dimensional example potential and a 3-dimensional simulation of the motion of an atomic force microscope tip moving along a surface.

DNA i-motif provides steel-like tough ends to chromosomes

ABSTRACTWe studied the structure and mechanical properties of DNA i-motif nanowires by means of molecular dynamics computer simulations. We built up to 230 nm-long nanowires, based on a repeated TC5 sequence from NMR crystallographic data, fully relaxed and equilibrated in water. The unusual C●C+ stacked structure, formed by four ssDNA strands arranged in an intercalated tetramer, is here fully characterized both statically and dynamically. By applying stretching, compression and bending deformations with the steered molecular dynamics and umbrella sampling methods, we extract the apparent Young’s and bending moduli of the nanowire, as well as estimates for the tensile strength and persistence length. According to our results, i-motif nanowires share similarities with structural proteins, as far as their tensile stiffness, but are closer to nucleic acids and flexible proteins, as far as their bending rigidity is concerned. Curiously enough, their tensile strength makes such DNA fragments tough as mild steel or a nickel alloy. Besides their yet to be clarified biological significance, i-motif nanowires may qualify as interesting candidates for nanotechnology templates, due to such outstanding mechanical properties.

Characterizing Conformational Ensemble and Free Energy Landscape of ABC Exporters using a Novel System-Specific Sampling Approach

During their transport cycle, ATP-binding cassette (ABC) exporters undergo large-scale conformational changes between inward-facing (IF) and outward-facing (OF) states. Using a novel approach based on designing system-specific reaction coordinates and employing nonequilibrium simulations, we have recently characterized the IF-OF conformational transition pathway of MsbA, a bacterial ABC exporter whose structure has been solved in multiple functional states including an OF conformation and two IF conformations (termed IF-closed and IF-open based on the extent to which the cytoplasmic side is open). Due to low resolution of MsbA crystal structures in IF conformation, only Cα atoms have been reported; however our unique system-specific sampling strategy results in a stable all-atom model of MsbA in the IF conformation in the presence of explicit membrane and water. In order to more accurately characterize the conformational ensemble of nucleotide-free apo MsbA, we also employ a free energy calculation technique based on a combination of umbrella sampling and replica exchange methods, using a reaction coordinate defined based on the orientation of transmembrane helices. The results indicate that the IF-closed conformation is associated with a relatively deep minimum and the IF-open-like conformations are also thermally accessible. However, the deepest free energy basin is associated with an IF conformation, more open than IF-closed and less open than IF-open, resembling the crystal structures of P-glycoprotein (a MsbA homolog) which are obtained at higher resolutions. The overall picture emerging from the reconstructed free energies is that MsbA is fairly flexible in its resting state in the absence of nucleotides and substrates. The approach proposed here provides a framework to study large-scale conformational changes of membrane transporters whose computational investigation at an atomic resolution may not be currently feasible using conventional methods.

Energetic View on Membrane Pore Formation

Energetic View on Membrane Pore Formation

Interactions of Fullerene (C60) and its Hydroxyl Derivatives with Lipid Bilayer: A Coarse-Grained Molecular Dynamics Simulation

Coarse-grained molecular dynamic simulations were employed to study the interactions of fullerene (C60) and its hydroxyl derivatives (C60(OH)n, n = 4, 5, 6, 8, 12, and 16) with a lipid bilayer composed of dipalmitoylphosphatidylcholine molecules. It was found that the C60 moves towards the center of the bilayer and laid between central and peripheral regions of the bilayer. The potential mean force was calculated to estimate free energy profile when pulling the fullerene from its initial position to the center of the bilayer using an umbrella sampling method. Results showed that the hydrophobic region of the membrane acts as a barrier to transport a nonpolar C60 molecule through the bilayer. This makes a deep minimum in the free energy profile between the center and head regions of membrane. Various numbers of polar functional groups (–OH) were then used to make derivatives of fullerene and change the hydrophilic of the molecule. It was found that optimal number of hydroxyl groups to facilitate the transportation of C60(OH)n through the bilayer is 4.

Structure and Mechanical Characterization of DNA i-Motif Nanowires by Molecular Dynamics Simulation

We studied the structure and mechanical properties of DNA i-motif nanowires by means of molecular dynamics computer simulations. We built up to 230 nm-long nanowires, based on a repeated TC5 sequence from crystallographic data, fully relaxed and equilibrated in water. The unusual C⋅C+ stacked structure, formed by four ssDNA strands arranged in an intercalated tetramer, is here fully characterized both statically and dynamically. By applying stretching, compression, and bending deformations with the steered molecular dynamics and umbrella sampling methods, we extract the apparent Young’s and bending moduli of the nanowire, as well as estimates for the tensile strength and persistence length. According to our results, the i-motif nanowire shares similarities with structural proteins, as far as its tensile stiffness, but is closer to nucleic acids and flexible proteins, as far as its bending rigidity is concerned. Furthermore, thanks to its very thin cross section, the apparent tensile toughness is close to that of a metal. Besides their yet to be clarified biological significance, i-motif nanowires may qualify as interesting candidates for nanotechnology templates, due to such outstanding mechanical properties.

Dynamics of ions in the selectivity filter of the KcsA channel

The statistical and dynamical properties of ions in the selectivity filter of the KcsA ion channel are considered on the basis of molecular dynamics (MD) simulations of the KcsA protein embedded in a lipid membrane surrounded by an ionic solution. A new approach to the derivation of a Brownian dynamics (BD) model of ion permeation through the filter is discussed, based on unbiased MD simulations. It is shown that depending on additional assumptions, ion's dynamics can be described either by under-damped Langevin equation with constant damping and white noise or by Langevin equation with a fractional memory kernel. A comparison of the potential of the mean force derived from unbiased MD simulations with the potential produced by the umbrella sampling method demonstrates significant differences in these potentials. The origin of these differences is an open question that requires further clarifications.

QuanPol: A full spectrum and seamless QM/MM program

The quantum chemistry polarizable force field program (QuanPol) is implemented to perform combined quantum mechanical and molecular mechanical (QM/MM) calculations with induced dipole polarizable force fields and induced surface charge continuum solvation models. The QM methods include Hartree-Fock method, density functional theory method (DFT), generalized valence bond theory method, multiconfiguration self-consistent field method, Møller-Plesset perturbation theory method, and time-dependent DFT method. The induced dipoles of the MM atoms and the induced surface charges of the continuum solvation model are self-consistently and variationally determined together with the QM wavefunction. The MM force field methods can be user specified, or a standard force field such as MMFF94, Chemistry at Harvard Molecular Mechanics (CHARMM), Assisted Model Building with Energy Refinement (AMBER), and Optimized Potentials for Liquid Simulations-All Atom (OPLS-AA). Analytic gradients for all of these methods are implemented so geometry optimization and molecular dynamics (MD) simulation can be performed. MD free energy perturbation and umbrella sampling methods are also implemented.

Two‐dimensional replica‐exchange method for predicting protein–ligand binding structures

We have developed a two-dimensional replica-exchange method for the prediction of protein-ligand binding structures. The first dimension is the umbrella sampling along the reaction coordinate, which is the distance between a protein binding pocket and a ligand. The second dimension is the solute tempering, in which the interaction between a ligand and a protein and water is weakened. The second dimension is introduced to make a ligand follow the umbrella potential more easily and enhance the binding events, which should improve the sampling efficiency. As test cases, we applied our method to two protein-ligand complex systems (MDM2 and HSP 90-alpha). Starting from the configuration in which the protein and the ligand are far away from each other in each system, our method predicted the ligand binding structures in excellent agreement with the experimental data from Protein Data Bank much faster with the improved sampling efficiency than the replica-exchange umbrella sampling method that we have previously developed.

Molecular Mechanisms, Thermodynamics, and Dissociation Kinetics of Knob-Hole Interactions in Fibrin

Polymerization of fibrin, the primary structural protein of blood clots and thrombi, occurs through binding of knobs 'A' and 'B' in the central nodule of fibrin monomer to complementary holes 'a' and 'b' in the γ- and β-nodules, respectively, of another monomer. We characterized the A:a and B:b knob-hole interactions under varying solution conditions using molecular dynamics simulations of the structural models of fibrin(ogen) fragment D complexed with synthetic peptides GPRP (knob 'A' mimetic) and GHRP (knob 'B' mimetic). The strength of A:a and B:b knob-hole complexes was roughly equal, decreasing with pulling force; however, the dissociation kinetics were sensitive to variations in acidity (pH 5-7) and temperature (T = 25-37 °C). There were similar structural changes in holes 'a' and 'b' during forced dissociation of the knob-hole complexes: elongation of loop I, stretching of the interior region, and translocation of the moveable flap. The disruption of the knob-hole interactions was not an "all-or-none" transition as it occurred through distinct two-step or single step pathways with or without intermediate states. The knob-hole bonds were stronger, tighter, and more brittle at pH 7 than at pH 5. The B:b knob-hole bonds were weaker, looser, and more compliant than the A:a knob-hole bonds at pH 7 but stronger, tighter, and less compliant at pH 5. Surprisingly, the knob-hole bonds were stronger, not weaker, at elevated temperature (T = 37 °C) compared with T = 25 °C due to the helix-to-coil transition in loop I that helps stabilize the bonds. These results provide detailed qualitative and quantitative characteristics underlying the most significant non-covalent interactions involved in fibrin polymerization.

Investigation and optimisation of the drying of reduced glutathione by steered molecular dynamics simulation

The drying of reduced glutathione from a series of aqueous–ethanol binary solutions at 300 K (below human body temperature) and 330 K (above human body temperature) was investigated in detail by steered molecular simulation and an umbrella sampling method with the Gromacs software package and Gromos96(53a6) united atomic force field. The results show that electrostatic interactions between glutathione and solvent represent the main resistance to drying. When the aqueous solution was gradually changed to pure ethanol, the energy of electrostatic interaction between glutathione and solvent molecules increased by 445.088 kJ/mol, and the drying potential of mean force (PMF) free energy also fell by 253.040 kJ/mol. However, an increase in temperature from 300 to 330 K in the aqueous solution only results in an increase of 23.013 kJ/mol in electrostatic interaction energy and a decrease of 34.956 kJ/mol in drying PMF free energy. Furthermore, we show that hydrogen bonding is the major form of electrostatic interaction involved, and directly affects the drying of glutathione. Therefore, choosing water-miscible solvents that minimise hydrogen-bond formation with glutathione will enhance its drying rate, and this is likely to be more efficient than increasing the temperature of the process. Thus, a power-saving technology can be used to produce the high bioactivity medicines.

Insight into the Oseltamivir Resistance R292K Mutation in H5N1 Influenza Virus: A Molecular Docking and Molecular Dynamics Approach

H5N1 is a subtype of the influenza A virus that can cause disease in humans and many other animal species. Oseltamivir (Tamiflu) is a potent and selective antiviral drug employed to fight the flu virus in infected individuals by inhibiting neuraminidase (NA), a flu protein responsible for the release and spread of the progeny virions. However, oseltamivir resistance has become a critical problem. In particular, influenza strains with a R292K NA mutation are highly resistant to the oseltamivir. Though the biological functions of the mutations have previously been characterized, the structural basis behind the reduced catalytic activity and reduced protein level is not clear. In this study, molecular docking and molecular dynamics (MD) approach were employed to investigate the structural and dynamical effects throughout the protein structure and specifically, at the drug-binding pocket. Furthermore, potential of mean force was analyzed using explicit solvent MD simulations with the umbrella sampling method to explore the free energy of binding. It is believed that this study provides valuable guidance for the resistance management of oseltamivir and designing of more potent antiviral inhibitor.

Aggregation Thermodynamics of Sodium Octanoate Micelles Studied by Means of Molecular Dynamics Simulations

The present work is aimed at studying the computation of the thermodynamic potentials that describe the stability of anionic surfactant molecules in micellar aggregates. We report a set of molecular dynamics simulations of a sodium octanoate micelle in aqueous solution using the umbrella sampling method along with the Jarzynski equality in order to compute the potential of mean force for the dissociation process of one surfactant molecule from a previously assembled micellar aggregate. The Jarzynski average was computed at several different temperatures in order to estimate the Gibbs energy of association for the octanoate anion, which was split into its enthalpic and entropic contributions. We also estimated the contributions arising from the polar head and the apolar tail for each thermodynamic potential, and a detailed picture emerged from these simulations. The aggregation is driven mostly by the Gibbs energy contribution arising from the hydrophobic tail, which was large enough to cancel out the unfavorable contribution from the polar head. Although the association process may be ascribed mostly to the transfer of the apolar tail to the micellar core, it should be noted that the polar head also contributed with a favorable entropic term to the overall Gibbs energy. These findings were rationalized by comparing the energetic and structural patterns of the hydration process of a free monomer in solution to an aggregated molecule. The interaction energy distributions presented at least two discernible populations and each population was related to a different structural pattern, as characterized by the radial distribution functions. Altogether, the changes in both the energy and structure of the hydration layer are consistent with the entropy-driven association of the surfactant into the micellar aggregate.

Hydration Structure at the α-Al2O3 (0001) Surface: Insights from Experimental Atomic Force Spectroscopic Data and Atomistic Molecular Dynamics Simulations

Atomic force spectroscopic data obtained in water on the (0001) face of fully hydroxylated α-Al2O3 substrate using a silicon tip are presented. The data are obtained by implementing the Brownian force profile reconstruction method (BFPR), originally proposed by Ashby and Lieber [J. Am. Chem. Soc. 2004, 126, 16973]. The method allows for an accurate reconstruction of the high stiffness force–distance curve that highlights the hydration structure. The experimental data are interpreted with the aid of massive atomistic molecular dynamics simulations in which one silicon dioxide disc of ∼2 nm diameter represents the silicon atomic force microscopy (AFM) tip. The umbrella sampling method is employed to obtain the water-mediated surface–disc force profile. It is possible to distinguish two hydration layers confined between tip and surface in both simulations and experiment. Small variations of the disc features yield some differences in the simulated force–distance curve, and the small disc size is responsible ...

Computing the phase diagram of binary mixtures: A patchy particle case study

We investigate the phase behaviour of 2D mixtures of bi-functional and three-functional patchy particles and 3D mixtures of bi-functional and tetra-functional patchy particles by means of Monte Carlo simulations and Wertheim theory. We start by computing the critical points of the pure systems and then we investigate how the critical parameters change upon lowering the temperature. We extend the successive umbrella sampling method to mixtures to make it possible to extract information about the phase behaviour of the system at a fixed temperature for the whole range of densities and compositions of interest.

Self-Learning Adaptive Umbrella Sampling Method for the Determination of Free Energy Landscapes in Multiple Dimensions.

The potential of mean force describing conformational changes of biomolecules is a central quantity that determines the function of biomolecular systems. Calculating an energy landscape of a process that depends on three or more reaction coordinates might require a lot of computational power, making some of multidimensional calculations practically impossible. Here, we present an efficient automatized umbrella sampling strategy for calculating multidimensional potential of mean force. The method progressively learns by itself, through a feedback mechanism, which regions of a multidimensional space are worth exploring and automatically generates a set of umbrella sampling windows that is adapted to the system. The self-learning adaptive umbrella sampling method is first explained with illustrative examples based on simplified reduced model systems, and then applied to two non-trivial situations: the conformational equilibrium of the pentapeptide Met-enkephalin in solution and ion permeation in the KcsA potassium channel. With this method, it is demonstrated that a significant smaller number of umbrella windows needs to be employed to characterize the free energy landscape over the most relevant regions without any loss in accuracy.