PyRETIS 2: An improbability drive for rare events

The algorithmic development in the field of path sampling has made tremendous progress in recent years. Although the original transition path sampling method was mostly used as a qualitative tool to sample reaction paths, the more recent family of interface-based path sampling methods has paved the way for more quantitative rate calculation studies. Of the exact methods, the replica exchange transition interface sampling (RETIS) method is the most efficient, but rather difficult to implement. This has been the main motivation to develop the open-source Python-based computer library PyRETIS that was released in 2017. PyRETIS is designed to be easily interfaced with any molecular dynamics (MD) package using either classical or ab initio MD. In this study, we report on the principles and the software enhancements that are now included in PyRETIS 2, as well as the recent developments on the user interface, improvements of the efficiency via the implementation of new shooting moves, easier initialization procedures, analysis methods, and supported interfaced software. © 2019 The Authors. Journal of Computational Chemistry published by Wiley Periodicals, Inc.     

Structural dissimilarity sampling with dynamically self‐guiding selection

Structural dissimilarity sampling (SDS) has been proposed as an enhanced conformational sampling method for reproducing the structural transitions of a given protein. SDS consists of cycles of two steps: (1) Selections of initial structures with structural dissimilarities by referring to a measure. (2) Conformational resampling by restarting short‐time molecular dynamics (MD) simulations from the initial structures. In the present study, an efficient measure is proposed as a dynamically self‐guiding selection to accelerate the structural transitions from a reactant state to a product state as an extension to the original SDS. In the extended SDS, the inner product (IP ) between the reactant and the snapshots generated by short‐time MD simulations are evaluated and ranked according to the IP s at every cycle. Then, the snapshots with low IP s are selected as initial structures for the short‐time MD simulations. This scheme enables one to choose dissimilar and distant initial structures from the reactant, and thus the initial structures dynamically head towards the product, promoting structural transitions from the reactant. To confirm the conformational sampling efficiency, the extended SDS was applied to maltodextrin binding protein (MBP), and we successfully reproduced the structural transition from the open to closed states with submicrosecond‐order simulation times. However, a conventional long‐time MD simulation failed to reproduce the same structural transition. We also compared the performance with that obtained by the ordinary SDS and other sampling techniques that have been developed by us to characterize the possible utility of the extended SDS for actual applications. © 2017 Wiley Periodicals, Inc.     

Vase‐Kite Equilibrium of Resorcin[4]arene Cavitands Investigated Using Molecular Dynamics Simulations with Ball‐and‐Stick Local Elevation Umbrella Sampling

A key feature of resorcin[4]arene cavitands is their ability to switch between a closed/contracted (Vase ) and an open/expanded (Kite ) conformation. The mechanism and dynamics of this interconversion remains, however, elusive. In the present study, the Vase ‐Kite transitions of a quinoxaline‐based and of a dinitrobenzene‐based resorcin[4]arene are investigated using molecular dynamics (MD) simulations in three environments (vacuum, chloroform, and toluene) and at three temperatures (198.15, 248.15, and 298.15 K). The challenge of sampling the Vase ‐Kite transition, which occurs experimentally on the millisecond time scale, is overcome by calculating relative free energies using ball‐and stick local elevation umbrella sampling (B&S‐LEUS) to enhance the statistics on the relevant states and to promote interconversion transitions. Associated unbiased MD simulations also evidence for the first time a complete Vase ‐to‐Kite transition, as well as transitions between degenerate Kite 1 and Kite 2 forms and solvent‐exchange events. The calculated Vase ‐to‐Kite free‐energy changes ΔG are in qualitative agreement with the experimental magnitudes and trends. The level of quantitative agreement is, however, limited by the force‐field accuracy and, in particular, by the approximate treatment of intramolecular interactions at the classical level. The results are in line with a less stable Vase state for the dinitrobenzene compared to the quinoxaline compound, and a negative entropy change ΔS for the Vase ‐to‐Kite transition of the latter compound. Relative free energies calculated for intermediates also suggest that the Vase ‐Kite transition does not follow a concerted mechanism, but an asynchronous one with sequential opening of the flaps. In particular, the conformation involving two adjacent flaps open in a parallel direction (cis‐p) represents a likely intermediate, which has not been observed experimentally to date.     

An extensible interface for QM/MM molecular dynamics simulations with AMBER

We present an extensible interface between the AMBER molecular dynamics (MD) software package and electronic structure software packages for quantum mechanical (QM) and mixed QM and classical molecular mechanical (MM) MD simulations within both mechanical and electronic embedding schemes. With this interface, ab initio wave function theory and density functional theory methods, as available in the supported electronic structure software packages, become available for QM/MM MD simulations with AMBER. The interface has been written in a modular fashion that allows straight forward extensions to support additional QM software packages and can easily be ported to other MD software. Data exchange between the MD and QM software is implemented by means of files and system calls or the message passing interface standard. Based on extensive tests, default settings for the supported QM packages are provided such that energy is conserved for typical QM/MM MD simulations in the microcanonical ensemble. Results for the free energy of binding of calcium ions to aspartate in aqueous solution comparing semiempirical and density functional Hamiltonians are shown to demonstrate features of this interface. © 2013 Wiley Periodicals, Inc.     

Binding structures of tri‐N‐acetyl‐β‐glucosamine in hen egg white lysozyme using molecular dynamics with a polarizable force field

Lysozyme is a well‐studied enzyme that hydrolyzes the β‐(1,4)‐glycosidic linkage of N‐acetyl‐β‐glucosamine (NAG)_n oligomers. The active site of hen egg‐white lysozyme (HEWL) is believed to consist of six subsites, A‐F that can accommodate six sugar residues. We present studies exploring the use of polarizable force fields in conjunction with all‐atom molecular dynamics (MD) simulations to analyze binding structures of complexes of lysozyme and NAG trisaccharide, (NAG)₃. MD trajectories are applied to analyze structures and conformation of the complex as well as protein–ligand interactions, including the hydrogen‐bonding network in the binding pocket. Two binding modes (ABC and BCD) of (NAG)₃ are investigated independently based on a fixed‐charge model and a polarizable model. We also apply molecular mechanics with generalized born and surface area (MM‐GBSA) methods based on MD using both nonpolarizable and polarizable force fields to compute binding free energies. We also study the correlation between root‐mean‐squared deviation and binding free energies of the wildtype and W62Y mutant; we find that for this prototypical system, approaches using the MD trajectories coupled with implicit solvent models are equivalent for polarizable and fixed‐charge models. © 2012 Wiley Periodicals, Inc.     

Increasing the time step with mass scaling in Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations

Born‐Oppenheimer ab initio QM/MM molecular dynamics simulation with umbrella sampling is a state‐of‐the‐art approach to calculate free energy profiles of chemical reactions in complex systems. To further improve its computational efficiency, a mass‐scaling method with the increased time step in MD simulations has been explored and tested. It is found that by increasing the hydrogen mass to 10 amu, a time step of 3 fs can be employed in ab initio QM/MM MD simulations. In all our three test cases, including two solution reactions and one enzyme reaction, the resulted reaction free energy profiles with 3 fs time step and mass scaling are found to be in excellent agreement with the corresponding simulation results using 1 fs time step and the normal mass. These results indicate that for Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations with umbrella sampling, the mass‐scaling method can significantly reduce its computational cost while has little effect on the calculated free energy profiles. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

Conformational Selection in Glycomimetics: Human Galectin‐1 Only Recognizes syn‐Ψ‐Type Conformations of β‐1,3‐Linked Lactose and Its C‐Glycosyl Derivative

The human lectin galectin‐1 (hGal‐1) translates sugar signals, that is, β‐galactosides, into effects on the level of cells, for example, growth regulation, and has become a model for studying binding of biopharmaceutically relevant derivatives. Bound‐state conformations of Galβ‐C‐(1→3)‐Glcβ‐OMe ( 1 ) and its βGal‐(1→3)‐βGlc‐OMe disaccharide parent compound were studied by using NMR spectroscopy (transferred (TR)‐NOESY data), assisted by docking experiments and molecular dynamics (MD) simulations. The molecular recognition process involves a conformational selection event. Although free C‐glycoside access four distinct conformers in solution, hGal‐1 recognizes shape of a local minimum of compound 1 , the syn‐Φ/syn‐Ψ conformer, not the structure at global minimum. MD simulations were run to explain, in structural terms, the observed geometry of the complex.     

All‐atom Molecular Dynamic Simulations Combined with the Chemical Shifts Study on the Weak Interactions of Ethanol‐water System

All‐atom molecular dynamics (MD) simulation combined with chemical shifts was performed to investigate the interactions over the entire concentration range of the ethanol (EtOH)‐water system. The results of the simulation were adopted to explain the NMR experiments by hydrogen bonding analysis. The strong hydrogen bonds and weak C–H···O contacts coexist in the mixtures through the analysis of the radial distribution functions. And the liquid structures in the whole concentration of EtOH‐water mixtures can be classified into three regions by the statistic analysis of the hydrogen‐bonding network in the MD simulations. Moreover, the chemical shifts of the hydrogen atom are in agreement with the statistical results of the average number hydrogen bonds in the MD simulations. Interestingly, the excess relative extent of η_rel^E calculated by the MD simulations and chemical shifts in the EtOH aqueous solutions shows the largest deviation at x_EtOH≈0.18. The excess properties present good agreement with the excess enthalpy in the concentration dependence.     

Comparison of radii sets,entropy, QM methods,and sampling on MM‐PBSA,MM‐GBSA,and QM/MM‐GBSA ligand binding energies of F. tularensis enoyl‐ACP reductase (FabI)

To validate a method for predicting the binding affinities of FabI inhibitors, three implicit solvent methods, MM‐PBSA, MM‐GBSA, and QM/MM‐GBSA were carefully compared using 16 benzimidazole inhibitors in complex with Francisella tularensis FabI. The data suggests that the prediction results are sensitive to radii sets, GB methods, QM Hamiltonians, sampling protocols, and simulation length, if only one simulation trajectory is used for each ligand. In this case, QM/MM‐GBSA using 6 ns MD simulation trajectories together with GB^neck2, PM3, and the mbondi2 radii set, generate the closest agreement with experimental values (r² = 0.88). However, if the three implicit solvent methods are averaged from six 1 ns MD simulations for each ligand (called “multiple independent sampling”), the prediction results are relatively insensitive to all the tested parameters. Moreover, MM/GBSA together with GB^HCT and mbondi, using 600 frames extracted evenly from six 0.25 ns MD simulations, can also provide accurate prediction to experimental values (r² = 0.84). Therefore, the multiple independent sampling method can be more efficient than a single, long simulation method. Since future scaffold expansions may significantly change the benzimidazole's physiochemical properties (charges, etc.) and possibly binding modes, which may affect the sensitivities of various parameters, the relatively insensitive “multiple independent sampling method” may avoid the need of an entirely new validation study. Moreover, due to large fluctuating entropy values, (QM/)MM‐P(G)BSA were limited to inhibitors’ relative affinity prediction, but not the absolute affinity. The developed protocol will support an ongoing benzimidazole lead optimization program. © 2015 Wiley Periodicals, Inc.     

Enhanced sampling method in molecular simulations using genetic algorithm for biomolecular systems

We propose a molecular simulation method using genetic algorithm (GA) for biomolecular systems to obtain ensemble averages efficiently. In this method, we incorporate the genetic crossover, which is one of the operations of GA, to any simulation method such as conventional molecular dynamics (MD), Monte Carlo, and other simulation methods. The genetic crossover proposes candidate conformations by exchanging parts of conformations of a target molecule between a pair of conformations during the simulation. If the candidate conformations are accepted, the simulation resumes from the accepted ones. While conventional simulations are based on local update of conformations, the genetic crossover introduces global update of conformations. As an example of the present approach, we incorporated genetic crossover to MD simulations. We tested the validity of the method by calculating ensemble averages and the sampling efficiency by using two kinds of peptides, ALA3 and (AAQAA)₃. The results show that for ALA3 system, the distribution probabilities of backbone dihedral angles are in good agreement with those of the conventional MD and replica-exchange MD simulations. In the case of (AAQAA)₃ system, our method showed lower structural correlation of α-helix structures than the other two methods and more flexibility in the backbone ψ angles than the conventional MD simulation. These results suggest that our method gives more efficient conformational sampling than conventional simulation methods based on local update of conformations. © 2018 Wiley Periodicals, Inc.     

Comprehensive Characterization of the Self-Folding Cavitand Dynamics

The conformational equilibria and guest exchange process of a resorcin[4]arene derived self-folding cavitand receptor have been characterized in detail by molecular dynamics simulations (MD) and ¹H EXSY NMR experiments. A multi-timescale strategy for exploring the fluxional behaviour of this system has been constructed, exploiting conventional MD and accelerated MD (aMD) techniques. The use of aMD allows the reconstruction of the folding/unfolding process of the receptor by sampling high-energy barrier processes unattainable by conventional MD simulations. We obtained MD trajectories sampling events occurring at different timescales from ns to s: 1) rearrangement of the directional hydrogen bond seam stabilizing the receptor, 2) folding/unfolding of the structure transiting partially open intermediates, and 3) guest departure from different folding stages. Most remarkably, reweighing of the biased aMD simulations provided kinetic barriers that are in very good agreement with those determined experimentally by ¹H NMR. These results constitute the first comprehensive characterization of the complex dynamic features of cavitand receptors. Our approach emerges as a valuable rational design tool for synthetic host-guest systems     

Estimating kinetic rates from accelerated molecular dynamics simulations: alanine dipeptide in explicit solvent as a case study

Molecular dynamics (MD) simulation is the standard computational technique used to obtain information on the time evolution of the conformations of proteins and many other molecular systems. However, for most biological systems of interest, the time scale for slow conformational transitions is still inaccessible to standard MD simulations. Several sampling methods have been proposed to address this issue, including the accelerated molecular dynamics method. In this work, we study the extent of sampling of the phi/psi space of alanine dipeptide in explicit water using accelerated molecular dynamics and present a framework to recover the correct kinetic rate constant for the helix to beta-strand transition. We show that the accelerated MD can drastically enhance the sampling of the phi/psi conformational phase space when compared to normal MD. In addition, the free energy density plots of the phi/psi space show that all minima regions are accurately sampled and the canonical distribution is recovered. Moreover, the kinetic rate constant for the helix to beta-strand transition is accurately estimated from these simulations by relating the diffusion coefficient to the local energetic roughness of the energy landscape. Surprisingly, even for such a low barrier transition, it is difficult to obtain enough transitions to accurately estimate the rate constant when one uses normal MD.     

On the “Tertiary Structure” of Poly‐Carbenes; Self‐Assembly of sp3‐Carbon‐Based Polymers into Liquid‐Crystalline Aggregates

The self‐assembly of poly(ethylidene acetate) (st‐PEA) into van der Waals‐stabilized liquid‐crystalline (LC) aggregates is reported. The LC behavior of these materials is unexpected, and unusual for flexible sp³‐carbon backbone polymers. Although the dense packing of polar ester functionalities along the carbon backbone of st‐PEA could perhaps be expected to lead directly to rigid‐rod behavior, molecular modeling reveals that individual st‐PEA chains are actually highly flexible and should not reveal rigid‐rod induced LC behavior. Nonetheless, st‐PEA clearly reveals LC behavior, both in solution and in the melt over a broad elevated temperature range. A combined set of experimental measurements, supported by MM/MD studies, suggests that the observed LC behavior is due to self‐aggregation of st‐PEA into higher‐order aggregates. According to MM/MD modeling st‐PEA single helices adopt a flexible helical structure with a preferred trans‐gauche syn‐syn‐anti‐anti orientation. Unexpectedly, similar modeling experiments suggest that three of these helices can self‐assemble into triple‐helical aggregates. Higher‐order assemblies were not observed in the MM/MD simulations, suggesting that the triple helix is the most stable aggregate configuration. DLS data confirmed the aggregation of st‐PEA into higher‐order structures, and suggest the formation of rod‐like particles. The dimensions derived from these light‐scattering experiments correspond with st‐PEA triple‐helix formation. Langmuir–Blodgett surface pressure–area isotherms also point to the formation of rod‐like st‐PEA aggregates with similar dimensions as st‐PEA triple helixes. Upon increasing the st‐PEA concentration, the viscosity of the polymer solution increases strongly, and at concentrations above 20 wt % st‐PEA forms an organogel. STM on this gel reveals the formation of helical aggregates on the graphite surface–solution interface with shapes and dimensions matching st‐PEA triple helices, in good agreement with the structures proposed by molecular modeling. X‐ray diffraction, WAXS, SAXS and solid state NMR spectroscopy studies suggest that st‐PEA triple helices are also present in the solid state, up to temperatures well above the melting point of st‐PEA. Formation of higher‐order aggregates explains the observed LC behavior of st‐PEA, emphasizing the importance of the “tertiary structure” of synthetic polymers on their material properties.     

Molecular Dynamics Simulation of Sub‐Transition for Polyethersulfone

Molecular dynamics (MD) simulations of a polyethersulfone (PES) chain are carried out in the amorphous state by using the Dreiding 2.21 force field at four temperatures. Two types of molecular motion, i.e. rotations of phenylene rings and torsions of large segments containing two oxygen atoms, two sulfur atoms, and five phenylene rings on the backbone, are simulated. The modeling results show that the successive phenylene rings should be in‐phase cooperative rotations, whereas the successive large segments should be out‐of‐phase cooperative torsions. By calculating the diffusion coefficient for the phenylene ring rotations, it is found that this rotation contributes to the β‐transition of PES.     

Validation of the GROMOS 54A7 Force Field Regarding Mixed α/β‐Peptide Molecules

A molecular‐dynamics (MD) simulation study of two heptapeptides containing α‐ and β‐amino acid residues is presented. According to NMR experiments, the two peptides differ in dominant fold when solvated in MeOH: peptide 3 adopts predominantly β‐hairpin‐like conformations, while peptide 8 adopts a 14/15‐helical fold. The MD simulations largely reproduce the experimental data. Application of NOE atom? atom distance restraining improves the agreement with experimental data, but reduces the conformational sampling. Peptide 3 shows a variety of conformations, while still agreeing with the NOE and ³J‐coupling data, whereas the conformational ensemble of peptide 8 is dominated by one helical conformation. The results confirm the suitability of the GROMOS 54A7 force field for simulation or structure refinement of mixed α/β‐peptides in MeOH.     

Accelerated Molecular Dynamics Study of Z/E Isomerization of Azobenzene: Kramers’ Theory Validation

The reaction rate in a condensed phase is retarded by solvent friction and deviates from the transition state theory (TST) in a highly viscous solvent. The Z/E isomerization kinetics of 4-dimethylamino-4′-nitroazobenzene in solution have been experimentally analyzed in a wide range of pressure in our previous studies. Herein the non-TST behaviors of the rate constant were theoretically modelled by means of Kramers’ theory and its extensions and by molecular dynamics (MD) simulations. The Z/E isomerization takes place in milliseconds or longer and so is not reachable by normal MD simulations. The accelerated MD method was employed to simulate the rare events using a biased isomerization potential barrier.     

Interaction potential models for bulk Zns,Zns nanoparticle,and Zns nanoparticle‐PMMA from first‐principles

An ab initio derived transferable polarizable force‐field has been developed for Zinc sulphide (ZnS) nanoparticle (NP) and ZnS NP‐PMMA nanocomposite. The structure and elastic constants of bulk ZnS using the new force‐field are within a few percent of experimental observables. The new force‐field show remarkable ability to reproduce structures and nucleation energies of nanoclusters (Zn₁S₁‐Zn₁₂S₁₂) as validated with that of the density functional theory calculations. A qualitative agreement of the radial distribution functions of Zn? O, in a ZnS nanocluster‐PMMA system, obtained using molecular mechanics molecular dynamics (MD) and ab initio MD (AIMD) simulations indicates that the ZnS–PMMA interaction through Zn? O bonding is explained satisfactorily by our force‐field. © 2015 Wiley Periodicals, Inc.     

Conformational Dynamics of Minimal Elastin‐Like Polypeptides: The Role of Proline Revealed by Molecular Dynamics and Nuclear Magnetic Resonance

Previous molecular dynamics studies of the elastin‐like peptide (ELP) GVG(VPGVG) predict that this ELP undergoes a conformational transition from an open to a more compact closed state upon an increase in temperature. These structural changes occurring in this minimal elastin model at the so‐called inverse temperature transition (ITT), which takes place when elastin is heated to temperatures of about 20–40 ^oC, are investigated further in this work by means of a combined theoretical and experimental approach. To do this, additional extensive classical molecular dynamics (MD) simulations of the capped octapeptide are carried out, analyzed, and compared to data obtained from homonuclear magnetic resonance (NMR) spectroscopy of the same octapeptide. Moreover, in the previous simulations, the proline residue in the ELP is found to act as a hinge, thereby allowing for the large‐amplitude opening and closing conformational motion of the ITT. To explore the role of proline in such elastin repeating units, a point mutant (P5I), which replaces the proline residue with an isoleucine residue, is also investigated using the aforementioned theoretical and experimental techniques. The results show that the site‐directed mutation completely alters the properties of this ELP, thus confirming the importance of the highly conserved proline residue in the ITT. Furthermore, a correlation between the two different methods employed is seen. Both methods predict the mutant ELP to be present in an unstructured form and the wild type ELP to have a β‐turn‐like structure. Finally, the role of the peptidyl cis to trans isomerization of the proline hinge is assessed in detail.     

Fast,metadynamics‐based method for prediction of the stereochemistry‐dependent relative free energies of ligand–receptor interactions

The computational approach applicable for the molecular dynamics (MD)‐based techniques is proposed to predict the ligand–protein binding affinities dependent on the ligand stereochemistry. All possible stereoconfigurations are expressed in terms of one set of force‐field parameters [stereoconfiguration‐independent potential (SIP)], which allows for calculating all relative free energies by only single simulation. SIP can be used for studying diverse, stereoconfiguration‐dependent phenomena by means of various computational techniques of enhanced sampling. The method has been successfully tested on the β₂‐adrenergic receptor (β₂‐AR) binding the four fenoterol stereoisomers by both metadynamics simulations and replica‐exchange MD. Both the methods gave very similar results, fully confirming the presence of stereoselective effects in the fenoterol‐β₂‐AR interactions. However, the metadynamics‐based approach offered much better efficiency of sampling which allows for significant reduction of the unphysical region in SIP. © 2014 Wiley Periodicals, Inc.     

Selective sampling of transition paths

In this short paper, we introduce an approximate method for the quick estimate of rate constants based on a simple sampling method of reactive transition paths over high energy barriers. It makes use of the previously introduced accelerated molecular dynamics (MD) simulation method to generate initial points for trajectory shooting. The accelerated MD simulations, although with the loss of real dynamics, lead to a quick calculation of thermodynamic properties and at the same time produce an ensemble of configurations with an enhanced sampling over the phase space that is more "reactive." The forward/backward trajectory shooting as that used in the transition path sampling method is then initiated from the configurations obtained from accelerated MD simulations to generate transition paths on the original unbiased potential. This method selectively enhances sampling of successful trajectories and at the same time accelerates significantly the calculation of rate constants.