Free Energy Path Research Articles

Mechanisms of Directed Self-Assembly in Cylindrical Hole Confinements

The directed self-assembly of block copolymers in cylindrical holes is a promising technology for lithographic patterning, particularly in the context of vertical interconnect accesses. While the hole-shrink process for single cylinders has been extensively explored, the proliferation of morphological defects remains a significant technological barrier. We use a coarse-grained model to explore morphologies that form within cylindrical confinements for combinations of template surface energies. We identify metastable defect morphologies, in addition to the desired cylindrical morphology, in majority-wetting sidewall templates. We use our coarse-grained model and the string method to identify transition pathways between defective morphologies and the cylindrical morphology to elucidate the mechanism of defect annihilation within the confinements; the transition pathway from a disordered state is also identified. This work demonstrates that the minimum free energy path for the formation of a cylinder goes th...

SSAGES: Software Suite for Advanced General Ensemble Simulations.

Molecular simulation has emerged as an essential tool for modern-day research, but obtaining proper results and making reliable conclusions from simulations requires adequate sampling of the system under consideration. To this end, a variety of methods exist in the literature that can enhance sampling considerably, and increasingly sophisticated, effective algorithms continue to be developed at a rapid pace. Implementation of these techniques, however, can be challenging for experts and non-experts alike. There is a clear need for software that provides rapid, reliable, and easy access to a wide range of advanced sampling methods and that facilitates implementation of new techniques as they emerge. Here we present SSAGES, a publicly available Software Suite for Advanced General Ensemble Simulations designed to interface with multiple widely used molecular dynamics simulations packages. SSAGES allows facile application of a variety of enhanced sampling techniques-including adaptive biasing force, string methods, and forward flux sampling-that extract meaningful free energy and transition path data from all-atom and coarse-grained simulations. A noteworthy feature of SSAGES is a user-friendly framework that facilitates further development and implementation of new methods and collective variables. In this work, the use of SSAGES is illustrated in the context of simple representative applications involving distinct methods and different collective variables that are available in the current release of the suite. The code may be found at: https://github.com/MICCoM/SSAGES-public.

First-principles theoretical assessment of catalysis by confinement: NO-O2 reactions within voids of molecular dimensions in siliceous crystalline frameworks.

Density functional theory methods that include dispersive forces are used to show how voids of molecular dimensions enhance reaction rates by the mere confinement of transition states analogous to those involved in homogeneous routes and without requiring specific binding sites or structural defects within confining voids. These van der Waals interactions account for the observed large rate enhancements for NO oxidation in the presence of purely siliceous crystalline frameworks. The minimum free energy paths for NO oxidation within chabazite (CHA) and silicalite (SIL) frameworks involve intermediates similar in stoichiometry, geometry, and kinetic relevance to those involved in the homogeneous route. The termolecular transition state for the kinetically-relevant cis-NOO2NO isomerization to trans-NOO2NO is strongly stabilized by confinement within CHA (by 36.3 kJ mol-1 in enthalpy) and SIL (by 39.2 kJ mol-1); such enthalpic stabilization is compensated, in part, by concomitant entropy losses brought forth by confinement (CHA: 44.9; SIL: 45.3, J mol-1 K-1 at 298 K). These enthalpy and entropy changes upon confinement agree well with those measured and combine to significantly decrease activation free energies and are consistent with the rate enhancements that become larger as temperature decreases because of the more negative apparent activation energies in confined systems compared with homogeneous routes. Calculated free energies of confinement are in quantitative agreement with measured rate enhancements and with their temperature sensitivity. Such quantitative agreements reflect preeminent effects of geometry in determining the van der Waals contributions from contacts between the transition states (TS) and the confining walls and the weak effects of the level of theory on TS geometries. NO oxidation reactions are chosen here to illustrate these remarkable effects of confinement because detailed kinetic analysis of rate data are available, but also because of their critical role in the treatment of combustion effluents and in the synthesis of nitric acid and nitrates. Similar effects are evident from rate enhancements by confinement observed for Diels-Alder and alkyne oligomerization reactions. These reactions also occur in gaseous media at near ambient temperatures, for which enthalpic stabilization upon confinement of their homogeneous transition states becomes the preeminent component of activation free energies.

Adaptive Finite Temperature String Method in Collective Variables.

Here we present a modified version of the on-the-fly string method for the localization of the minimum free energy path in a space of arbitrary collective variables. In the proposed approach the shape of the biasing potential is controlled by only two force constants, defining the width of the potential along the string and orthogonal to it. The force constants and the distribution of the string nodes are optimized during the simulation, improving the convergence. The optimized parameters can be used for umbrella sampling with a path CV along the converged string as the reaction coordinate. We test the new method with three fundamentally different processes: chloride attack to chloromethane in bulk water, alanine dipeptide isomerization, and the enzymatic conversion of isochorismate to piruvate. In each case the same set of parameters resulted in a rapidly converging simulation and a precise estimation of the potential of mean force. Therefore, the default settings can be used for a wide range of processes, making the method essentially parameter free and more user-friendly.

Reaction coordinates and transition states in enzymatic catalysis

Enzymatic reactions are complex chemical processes taking place in complex dynamic environments. Theoretical characterization of these reactions requires the determination of the reaction coordinate and the transition state ensemble. This is not an easy task because many degrees of freedom may be involved in principle. We present recent efforts to find good enzymatic reaction coordinates and the implications of these findings in the interpretation of enzymatic efficiency. In particular, we analyze different strategies based on the use of minimum free energy paths and direct localization of the dividing surface on multidimensional free energy surfaces. Another strategy is based on the generation of reactive trajectories, using the transition path sampling method, from which transition state configurations can be harvested. Most of the applications carried out until now coincide to stress the change in the nature of the reaction coordinate, in terms of the participation of the chemical and environmental degrees of freedom, as the reaction advances. The degrees of freedom of the chemical system are dominant at the transition state while environmental participation can be more important at early or late stages of the process. WIREs Comput Mol Sci 2018, 8:e1329. doi: 10.1002/wcms.1329This article is categorized under: Structure and Mechanism > Computational Biochemistry and Biophysics Molecular and Statistical Mechanics > Free Energy Methods Software > Simulation Methods

Rock climbing: A local-global algorithm to compute minimum energy and minimum free energy pathways.

The calculation of minimum energy or minimum free energy paths is an important step in the quantitative and qualitative studies of chemical and physical processes. The computations of these coordinates present a significant challenge and have attracted considerable theoretical and computational interest. Here we present a new local-global approach to study reaction coordinates, based on a gradual optimization of an action. Like other global algorithms, it provides a path between known reactants and products, but it uses a local algorithm to extend the current path in small steps. The local-global approach does not require an initial guess to the path, a major challenge for global pathway finders. Finally, it provides an exact answer (the steepest descent path) at the end of the calculations. Numerical examples are provided for the Mueller potential and for a conformational transition in a solvated ring system.

Ca2+ -Cl- Association in Water Revisited: the Role of Cation Hydration.

We investigate the dissociation of a Ca2+ -Cl- pair in water using classical molecular dynamics simulations with a polarizable interaction potential, parameterized from ab initio calculations. By computing the potential of mean force as a function not only of the interionic distance but also of the coordination numbers by water molecules, we show that it is necessary to use a collective variable describing the cation hydration in order to capture the dissociation mechanism. In the contact ion pair, the Ca2+ cation has a first coordination sphere containing 5 or 6 water molecules. The minimum free-energy path for dissociation involves a two-step process: First one or two additional water molecules enter the cation coordination shell, increasing the coordination number up to 7 with an almost fixed interionic distance. Then the dissociation of the ionic pair occurs at this fixed coordination number.

Equilibrium Ensembles for Insulin Folding from Bias-Exchange Metadynamics

Earliest events in the aggregation process, such as single molecule reconfiguration, are extremely important and the most difficult to characterize in experiments. To this end, we have used well-tempered bias exchange metadynamics simulations to determine the equilibrium ensembles of an insulin molecule under amyloidogenic conditions of low pH and high temperature. A bin-based clustering method that uses statistics accumulated in bias exchange metadynamics trajectories was employed to construct a detailed thermodynamic and kinetic model of insulin folding. The highest lifetime, lowest free-energy ensemble identified consisted of native conformations adopted by a folded insulin monomer in solution, namely, the R-, the Rf-, and the T-states of insulin. The lowest free-energy structure had a root mean square deviation of only 0.15 nm from native x-ray structure. The second longest-lived metastable state was an unfolded, compact monomer with little similarity to the native structure. We have identified three additional long-lived, metastable states from the bin-based model. We then carried out an exhaustive structural characterization of metastable states on the basis of tertiary contact maps and per-residue accessible surface areas. We have also determined the lowest free-energy path between two longest-lived metastable states and confirm earlier findings of non-two-state folding for insulin through a folding intermediate. The ensemble containing the monomeric intermediate retained 58% of native hydrophobic contacts, however, accompanied by a complete loss of native secondary structure. We have discussed the relative importance of nativelike versus nonnative tertiary contacts for the folding transition. We also provide a simple measure to determine the importance of an individual residue for folding transition. Finally, we have compared and contrasted this intermediate with experimental data obtained in spectroscopic, crystallographic, and calorimetric measurements during early stages of insulin aggregation. We have also determined stability of monomeric insulin by incubation at a very low concentration to isolate protein-protein interaction effects.

Minimum free-energy paths for the self-organization of polymer brushes.

A methodology to calculate minimum free-energy paths based on the combination of a molecular theory and the improved string method is introduced and applied to study the self-organization of polymer brushes under poor solvent conditions. Polymer brushes in a poor solvent cannot undergo macroscopic phase separation due to the physical constraint imposed by the grafting points; therefore, they microphase separate forming aggregates. Under some conditions, the theory predicts that the homogeneous brush and the aggregates can exist as two different minima of the free energy. The theoretical methodology introduced in this work allows us to predict the minimum free-energy path connecting these two minima as well as the morphology of the system along the path. It is shown that the transition between the homogeneous brush and the aggregates may involve a free-energy barrier or be barrierless depending on the relative stability of the two morphologies and the chain length and grafting density of the polymer. In the case where a free-energy barrier exists, one of the morphologies is a metastable structure and, therefore, the properties of the brush as the quality of the solvent is cycled are expected to display hysteresis. The theory is also applied to study the adhesion/deadhesion transition between two opposing surfaces modified by identical polymer brushes and it is shown that this process may also require surpassing a free-energy barrier.

Effective Riemannian Diffusion Model for Conformational Dynamics of Biomolecular Systems.

We present a Riemannian formalism for effective diffusion of biomolecules in collective variable spaces that provides a robust framework for conformational free energy calculation methods. Unlike their Euclidean counterparts, the Riemannian potential of mean force (PMF) and minimum free energy path (MFEP) are invariant under coordinate transformations. The presented formalism can be readily employed to modify the collective variable based enhanced sampling techniques, such as umbrella sampling (US) commonly used in biomolecular simulations, to take into account the role of intrinsic geometry of collective variable space. Although our model is mathematically equivalent to a Euclidean diffusion with a position-dependent diffusion tensor, the Riemannian formulation provides a more convenient framework for free energy calculation methods and path-finding algorithms aimed at characterizing the effective conformational dynamics of biomolecules. A simple three-dimensional toy model and a pentapeptide (met-enkephalin) simulated in an explicit solvent environment are used to illustrate the workings of the formalism and its implementation.

Different QM/MM Approaches To Elucidate Enzymatic Reactions: Case Study on ppGalNAcT2.

Hybrid QM/MM computational studies can provide invaluable insight into the mechanisms of enzymatic reactions that can be exploited for rational drug design. Various approaches are available for such studies. However, their strengths and weaknesses may not be immediately apparent. Using the retaining glycosyltransferase ppGalNAcT2 as a case study, we compare different methodologies used to obtain reaction paths and transition state information. Ab Initio MD using CPMD coupled with the String Method is used to derive the minimum free energy reaction path. The geometrical features of the free energy path, especially around the transition state, agree with the minimum potential energy path obtained by the much less computationally expensive Nudged Elastic Band method. The barrier energy, however, differs by 8 kcal/mol. The free energy surface generated by metadynamics provides a rough overview of the reaction and can confirm the physical relevance of optimized paths or provide an initial guess for path optimization methods. Calculations of enzymatic reactions are usually performed at best at the DFT level of theory. A comparison of widely used functionals with high-level DLPNO-CCSD(T)/CBS data on the potential energy profile serves as a validation of the usability of DFT for this type of enzymatic reaction. The M06-2X meta-hybrid functional in particular matches the DLPNO-CCSD(T)/CBS reference extremely well with errors within 1 kcal/mol. However, even pure-GGA functional OPBE provides sufficient accuracy for this type of study.

Heterogeneous nucleation from a supercooled ionic liquid on a carbon surface.

Classical molecular dynamics simulations were used to study the nucleation of the crystal phase of the ionic liquid [dmim+][Cl-] from its supercooled liquid phase, both in the bulk and in contact with a graphitic surface of D = 3 nm. By combining the string method in collective variables [Maragliano et al., J. Chem. Phys. 125, 024106 (2006)], with Markovian milestoning with Voronoi tessellations [Maragliano et al., J. Chem. Theory Comput. 5, 2589-2594 (2009)] and order parameters for molecular crystals [Santiso and Trout, J. Chem. Phys. 134, 064109 (2011)], we computed minimum free energy paths, the approximate size of the critical nucleus, the free energy barrier, and the rates involved in these nucleation processes. For homogeneous nucleation, the subcooled liquid phase has to overcome a free energy barrier of ∼85 kcal/mol to form a critical nucleus of size ∼3.6 nm, which then grows into the monoclinic crystal phase. This free energy barrier becomes about 42% smaller (∼49 kcal/mol) when the subcooled liquid phase is in contact with a graphitic disk, and the critical nucleus formed is about 17% smaller (∼3.0 nm) than the one observed for homogeneous nucleation. The crystal formed in the heterogeneous nucleation scenario has a structure that is similar to that of the bulk crystal, with the exception of the layers of ions next to the graphene surface, which have larger local density and the cations lie with their imidazolium rings parallel to the graphitic surface. The critical nucleus forms near the graphene surface separated only by these layers of ions. The heterogeneous nucleation rate (∼4.8 × 1011 cm-3 s-1) is about one order of magnitude faster than the homogeneous rate (∼6.6 × 1010 cm-3 s-1). The computed free energy barriers and nucleation rates are in reasonable agreement with experimental and simulation values obtained for the homogeneous and heterogeneous nucleation of other systems (ice, urea, Lennard-Jones spheres, and oxide glasses).

Transition-Tempered Metadynamics Is a Promising Tool for Studying the Permeation of Drug-like Molecules through Membranes.

Metadynamics is an important enhanced sampling technique in molecular dynamics simulation to efficiently explore potential energy surfaces. The recently developed transition-tempered metadynamics (TTMetaD) has been proven to converge asymptotically without sacrificing exploration of the collective variable space in the early stages of simulations, unlike other convergent metadynamics (MetaD) methods. We have applied TTMetaD to study the permeation of drug-like molecules through a lipid bilayer to further investigate the usefulness of this method as applied to problems of relevance to medicinal chemistry. First, ethanol permeation through a lipid bilayer was studied to compare TTMetaD with nontempered metadynamics and well-tempered metadynamics. The bias energies computed from various metadynamics simulations were compared to the potential of mean force calculated from umbrella sampling. Though all of the MetaD simulations agree with one another asymptotically, TTMetaD is able to predict the most accurate and reliable estimate of the potential of mean force for permeation in the early stages of the simulations and is robust to the choice of required additional parameters. We also show that using multiple randomly initialized replicas allows convergence analysis and also provides an efficient means to converge the simulations in shorter wall times and, more unexpectedly, in shorter CPU times; splitting the CPU time between multiple replicas appears to lead to less overall error. After validating the method, we studied the permeation of a more complicated drug-like molecule, trimethoprim. Three sets of TTMetaD simulations with different choices of collective variables were carried out, and all converged within feasible simulation time. The minimum free energy paths showed that TTMetaD was able to predict almost identical permeation mechanisms in each case despite significantly different definitions of collective variables.

Generalized Ensemble Sampling of Enzyme Reaction Free Energy Pathways.

Free energy path sampling plays an essential role in computational understanding of chemical reactions, particularly those occurring in enzymatic environments. Among a variety of molecular dynamics simulation approaches, the generalized ensemble sampling strategy is uniquely attractive for the fact that it not only can enhance the sampling of rare chemical events but also can naturally ensure consistent exploration of environmental degrees of freedom. In this review, we plan to provide a tutorial-like tour on an emerging topic: generalized ensemble sampling of enzyme reaction free energy path. The discussion is largely focused on our own studies, particularly ones based on the metadynamics free energy sampling method and the on-the-path random walk path sampling method. We hope that this minipresentation will provide interested practitioners some meaningful guidance for future algorithm formulation and application study.

A general method for molecular modeling of nucleation from the melt.

Crystallization is one of the fundamental phase transition processes, and it is also important practically, for example, in the chemical, food, and pharmaceutical industries. Despite its importance, however, our basic understanding of crystallization, and especially crystal nucleation, at the molecular level is still incomplete. In this work, we present a general molecular simulation approach that can be used to investigate the nucleation of crystals from a subcooled liquid. Our method combines a previously proposed general method to construct structure-based order parameters [E. E. Santiso and B. L. Trout, J. Chem. Phys. 134, 064109 (2011)] with the string method in collective variables [L. Maragliano et al., J. Chem. Phys. 125, 024106 (2006)] to obtain a minimum free energy path connecting the liquid and solid basins. We then use Markovian milestoning with Voronoi tessellations [E. Vanden-Eijnden and M. Venturoli, J. Chem. Phys. 130, 194101 (2009); L. Maragliano et al., J. Chem. Theory Comput. 5, 2589-2594 (2009)] to obtain the free energy profile along the path and the nucleation kinetics. We illustrate the application of this method to the nucleation of Benzene-I crystals from the melt, and compare the results to those previously found using transition path sampling [M. Shah et al., J. Phys. Chem. B 115, 10400-10412 (2011)].

Molecular simulation of homogeneous nucleation of crystals of an ionic liquid from the melt.

The homogeneous nucleation of crystals of the ionic liquid [dmim(+)][Cl(-)] from its supercooled liquid phase in the bulk (P = 1 bar, T = 340 K, representing a supercooling of 58 K) was studied using molecular simulations. The string method in collective variables [Maragliano et al., J. Chem. Phys. 125, 024106 (2006)] was used in combination with Markovian milestoning with Voronoi tessellations [Maragliano et al., J. Chem. Theory Comput. 5, 2589-2594 (2009)] and order parameters for molecular crystals [E. E. Santiso and B. L. Trout, J. Chem. Phys. 134, 064109 (2011)] to sketch a minimum free energy path connecting the supercooled liquid and the monoclinic crystal phases, and to determine the free energy and the rates involved in the homogeneous nucleation process. The physical significance of the configurations found along this minimum free energy path is discussed with the help of calculations based on classical nucleation theory and with additional simulation results obtained for a larger system. Our results indicate that, at a supercooling of 58 K, the liquid has to overcome a free energy barrier of the order of 60 kcal/mol and to form a critical nucleus with an average size of about 3.6 nm, before it reaches the thermodynamically stable crystal phase. A simulated homogeneous nucleation rate of 5.0 × 10(10) cm(-3) s(-1) was obtained for our system, which is in reasonable agreement with experimental and simulation rates for homogeneous nucleation of ice at similar degrees of supercooling. This study represents our first step in a series of studies aimed at understanding the nucleation and growth of crystals of organic salts near surfaces and inside nanopores.

Atomic‐resolution dissection of the energetics and mechanism of isomerization of hydrated ATP‐Mg2+ through the SOMA string method

The atomic mechanisms of isomerization of ATP-Mg(2+) in solution are characterized using the recently developed String Method with Optimal Molecular Alignment (SOMA) and molecular-dynamics simulations. Bias-Exchange Metadynamics simulations are first performed to identify the primary conformers of the ATP-Mg(2+) complex and their connectivity. SOMA is then used to elucidate the minimum free-energy path (MFEP) for each transition, in a 48-dimensional space. Analysis of the per-atom contributions to the global free-energy profiles reveals that the mechanism of these transitions is controlled by the Mg(2+) ion and its coordinating oxygen atoms in the triphosphate moiety, as well as by the ion-hydration shell. Metadynamics simulations in path collective variables based on the MFEP demonstrate these isomerizations proceed across a narrow channel of configurational space, thus validating the premise underlying SOMA. This study provides a roadmap for the examination of conformational changes in biomolecules, based on complementary enhanced-sampling techniques with different strengths. © 2015 Wiley Periodicals, Inc.

Nucleation and interfacial adsorption in ternary systems.

Nucleation is studied in incompressible ternary fluids by examining the topology of the overall landscape of the energy surface. Minimum free energy paths for nucleation (MFEPs) of a single nucleus in an infinite matrix are computed with the string method in the framework of the continuum theory of nucleation for the regular solution. Properties of the critical nucleus are compared with the predictions of the classical nucleation theory. MFEPs are found to exhibit complex nucleation pathways with non-monotonic variations of compositions in the interfacial region, specifically adsorption of a component. In the symmetric regular solution, the minority component is found to segregate at the interface during nucleation with a concomitant depletion of the nucleus core, resulting in unpredicted partition of the non-selective component. Despite increasing the gradient energy, such inhomogeneity in composition is shown to lower the nucleation barrier.

Regulation of the Molecular Brake Region in FGFR2 Kinase

Unregulated activity of tyrosine kinases (TKS) is responsible for numerous developmental musculoskeletal diseases and cancers. The regulatory mechanisms of tks and how they are affected by point mutations are still in need of elucidation. Fibroblast growth factor receptor (FGFR) kinase is a TK whose regulatory element, the “molecular brake,” is thought to prevent constitutive activation of the receptor, despite not being near the kinase's active site. germline mutations in this “molecular brake” region in FGFR cause musculoskeletal disorders such as craniosynostosis and dwarfism, while somatic mutations can lead to uninhibited cell growth, causing an array of cancers. furthermore, it is thought that the majority of pathogenic activating mutations in protein kinases are in “molecular brake”-like regions. However, the mechanism through which the “molecular brake” of FGFR kinase is disengaged, allowing for signal transduction, is still unknown. using path-based molecular dynamics (MD) simulations, we have shown that the inward motion of the kinase's activation loop upon tyrosine phosphorylation is correlated with disengagement of the “molecular brake.” results of the string method in collective variables, which finds the minimum free energy path connecting the loop-out and loop-in states, indicate found that inward motion of loop occurs within 1 frame of “molecular brake” disengagement. Furthermore, umbrella sampling shows that in loop-out conformations, the disengaged state of the “molecular brake” is less stable, while in loop-in conformations, the disengaged state becomes more stable. Our results suggest that the “molecular brake” is allosterically controlled by the activation loop motion. An understanding of the “molecular brake” disengagement mechanism will provide avenues for targeted therapies to counteract pathologic over-activation of fgfr kinase leading to uncontrolled signal transduction and uninhibited cell growth.

A fast tomographic method for searching the minimum free energy path.

Minimum Free Energy Path (MFEP) provides a lot of important information about the chemical reactions, like the free energy barrier, the location of the transition state, and the relative stability between reactant and product. With MFEP, one can study the mechanisms of the reaction in an efficient way. Due to a large number of degrees of freedom, searching the MFEP is a very time-consuming process. Here, we present a fast tomographic method to perform the search. Our approach first calculates the free energy surfaces in a sequence of hyperplanes perpendicular to a transition path. Based on an objective function and the free energy gradient, the transition path is optimized in the collective variable space iteratively. Applications of the present method to model systems show that our method is practical. It can be an alternative approach for finding the state-to-state MFEP.