Bridging coarse-grained models by jump-in-sample simulations

We present an efficient method to construct coarse-grained (CG) models from models of finer resolution. The method estimates the free energies in a generated sample of the CG conformational space and then fits the entire effective potential surface in the high-dimensional CG conformational space. A jump-in-sample algorithm that uses a random jumping walk in the CG sample is used to iteratively estimate the free energies. We test the method in a tetrahedral molecular fluid where we construct the intermolecular effective potential and evaluate the CG molecular model. Our algorithm for calculating the free energy involves an improved Wang-Landau (WL) algorithm, which not only works more efficiently than the standard WL algorithm, but also can work in high-dimensional spaces.     

Essential energy space random walk via energy space metadynamics method to accelerate molecular dynamics simulations

To overcome the possible pseudoergodicity problem, molecular dynamic simulation can be accelerated via the realization of an energy space random walk. To achieve this, a biased free energy function (BFEF) needs to be priori obtained. Although the quality of BFEF is essential for sampling efficiency, its generation is usually tedious and nontrivial. In this work, we present an energy space metadynamics algorithm to efficiently and robustly obtain BFEFs. Moreover, in order to deal with the associated diffusion sampling problem caused by the random walk in the total energy space, the idea in the original umbrella sampling method is generalized to be the random walk in the essential energy space, which only includes the energy terms determining the conformation of a region of interest. This essential energy space generalization allows the realization of efficient localized enhanced sampling and also offers the possibility of further sampling efficiency improvement when high frequency energy terms irrelevant to the target events are free of activation. The energy space metadynamics method and its generalization in the essential energy space for the molecular dynamics acceleration are demonstrated in the simulation of a pentanelike system, the blocked alanine dipeptide model, and the leucine model.     

MC(JBW): Simple but smart Monte Carlo algorithm for free energy simulations of multiconformational molecules

Many of the most common molecular simulation methods, including Monte Carlo (MC) and molecular or stochastic dynamics (MD or SD), have significant difficulties in sampling the space of molecular potential energy surfaces characterized by multiple conformational minima and significant energy barriers. In such cases improved sampling can be obtained by special techniques that lower such barriers or somehow direct search steps toward different low energy regions of space. We recently described a hybrid MC/SD algorithm [MC(JBW)/SD] incorporating such a technique that directed MC moves of selected torsion and bond angles toward known low energy regions of conformational space. Exploration of other degrees of freedom was left to the SD part of the hybrid algorithm. In the work described here, we develop a related but simpler simulation algorithm that uses only MC to sample all degrees of freedom (e.g., stretch, bend, and torsion). We term this algorithm MC(JBW). Using simulations on various model potential energy surfaces and on simple molecular systems (n-pentane, n-butane, and cyclohexane), MC(JBW) is shown to generate ensembles of states that are indistinguishable from the canonical ensembles generated by classical Metropolis MC in the limit of very long simulations. We further demonstrate the utility of MC(JBW) by evaluating the room temperature free energy differences between conformers of various substituted cyclohexanes and the larger ring hydrocarbons cycloheptane, cyclooctane, cyclononane, and cyclodecane. The results compare favorably with available experimental data and results from previously reported MC(JBW)/SD conformational free energy calculations. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1736–1745, 1998     

A Discrete,Space Variation Model for Studying the Kinetics of Shape Deformation of Vesicles Coupled with Phase Separation

Summary: The evolution dynamics of phase separation, coupled with shape deformation of vesicles is described by using dissipative dynamic equations, specifically the time‐dependent Ginzburg‐Landau (TDGL) equations. In order to improve the numerical stability and thus to efficiently deal with a large deformation of vesicles, a new algorithm, namely the discrete space variation model (DSVM) has been developed for the first time. The algorithm is based on the variation of the discretized free‐energy functional, which is constructed in discrete membrane space, in contrast to the commonly used continuous free‐energy functional. For the sake of numerical tractability, only the cylindrical vesicles (2D), with two components, are taken into consideration to illustrate the efficiency and validity of new algorithm. The simulation results, based on the DSVM algorithm have been compared with those from both linear analysis and strong segregation theory using the continuous space free‐energy functional. It is found that the DSVM algorithm can correctly describe the coupling between the lateral phase‐separation on the vesicle membrane and the vesicle shape deformation, both for early and late stages.

A flower‐like vesicle obtained by DSVM simulation.     

Multiple free energies from a single simulation: extending enveloping distribution sampling to nonoverlapping phase-space distributions

A recently proposed method to obtain free energy differences for multiple end states from a single simulation of a reference state which was called enveloping distribution sampling (EDS) [J. Chem. Phys. 126, 184110 (2007)] is expanded to situations where the end state configuration space densities do not show overlap. It uses a reference state Hamiltonian suggested by Han in 1992 [Phys. Lett. A 165, 28 (1992)] in a molecular dynamics implementation. The method allows us to calculate multiple free energy differences "on the fly" from a single molecular dynamics simulation. The influence of the parameters on the accuracy and precision of the obtained free energy differences is investigated. A connection is established between the presented method and the Bennett acceptance ratio method. The method is applied to four two-state test systems (dipole inversion, van der Waals perturbation, charge inversion, and water to methanol conversion) and two multiple-state test systems [dipole inversion with five charging states and five (dis-)appearing water molecules]. Accurate results could be obtained for all test applications if the parameters of the reference state Hamiltonian were optimized according to a given algorithm. The deviations from the exact result or from an independent calculation were at most 0.6 kJ/mol. An accurate estimation of the free energy difference is always possible, independent of how different the end states are. However, the convergence times of the free energy differences are longer in cases where the end state configuration space densities do not show overlap [charge inversion, water to methanol conversion, (dis-)appearing water molecules] than in cases where the configuration space densities do show some overlap [(multiple) dipole inversion and van der Waals perturbation].     

A hybrid recursion method to robustly ensure convergence efficiencies in the simulated scaling based free energy simulations

Recently, we developed an efficient free energy simulation technique, the simulated scaling (SS) method [H. Li et al., J. Chem. Phys. 126, 024106 (2007)], in the framework of generalized ensemble simulations. In the SS simulations, random walks in the scaling parameter space are realized so that both phase space overlap sampling and conformational space sampling can be simultaneously enhanced. To flatten the distribution in the scaling parameter space, in the original SS implementation, the Wang-Landau recursion was employed due to its well-known recursion capability. In the Wang-Landau recursion based SS free energy simulation scheme, at the early stage, recursion efficiencies are high and free energy regions are quickly located, although at this stage, the errors of estimated free energy values are large; at the later stage, the errors of estimated free energy values become smaller, however, recursions become increasingly slow and free energy refinements require very long simulation time. In order to robustly resolve this efficiency problem during free energy refinements, a hybrid recursion strategy is presented in this paper. Specifically, we let the Wang-Landau update method take care of the early stage recursion: the location of target free energy regions, and let the adaptive reweighting method take care of the late stage recursion: the refinements of free energy values. As comparably studied in the model systems, among three possible recursion procedures, the adaptive reweighting recursion approach is the least favorable one because of its low recursion efficiency during free energy region locations; and compared to the original Wang-Landau recursion approach, the proposed hybrid recursion technique can be more robust to guarantee free energy simulation efficiencies.     

Constructing a multidimensional free energy surface like a spider weaving a web

Complete free energy surface in the collective variable space provides important information of the reaction mechanisms of the molecules. But, sufficient sampling in the collective variable space is not easy. The space expands quickly with the number of the collective variables. To solve the problem, many methods utilize artificial biasing potentials to flatten out the original free energy surface of the molecule in the simulation. Their performances are sensitive to the definitions of the biasing potentials. Fast‐growing biasing potential accelerates the sampling speed but decreases the accuracy of the free energy result. Slow‐growing biasing potential gives an optimized result but needs more simulation time. In this article, we propose an alternative method. It adds the biasing potential to a representative point of the molecule in the collective variable space to improve the conformational sampling. And the free energy surface is calculated from the free energy gradient in the constrained simulation, not given by the negative of the biasing potential as previous methods. So the presented method does not require the biasing potential to remove all the barriers and basins on the free energy surface exactly. Practical applications show that the method in this work is able to produce the accurate free energy surfaces for different molecules in a short time period. The free energy errors are small in the cases of various biasing potentials. © 2017 Wiley Periodicals, Inc.     

Calibration of force-field dependency in free energy landscapes of peptide conformations by quantum chemical calculations

The free energy landscapes of peptide conformations were calibrated by ab initio quantum chemical calculations, after the enhanced conformational diversity search using the multicanonical molecular dynamics simulations. Three different potentials of mean force for an isolated dipeptide were individually obtained by the multicanonical molecular dynamics simulations using the conventional force fields, AMBER parm94, AMBER parm96, and CHARMm22. Each potential of mean force was then calibrated based upon the umbrella sampling algorithm from the adiabatic energy map that was calculated separately by the ab initio molecular orbital method, and all of the calibrated potentials of mean force coincided well. The calibration method was also applied to the simulations of a peptide dimer in explicit water models, and it was shown that the calibrated free energy landscapes did not depend on the force field used in the classical simulations, as far as the conformational space was sampled well. The current calibration method fuses the classical free energy calculation with the quantum chemical calculation, and it should generally make simulations for biomolecular systems much more reliable when combining with enhanced conformational sampling.     

Combination of molecular dynamics method and 3D-RISM theory for conformational sampling of large flexible molecules in solution

We have developed an algorithm for sampling the conformational space of large flexible molecules in solution, which combines the molecular dynamics (MD) method and the three-dimensional reference interaction site model (3D-RISM) theory. The solvent-induced force acting on solute atoms was evaluated as the gradient of the solvation free energy with respect to the solute-atom coordinates. To enhance the computation speed, we have applied a multiple timestep algorithm based on the RESPA (Reversible System Propagator Algorithm) to the combined MD/3D-RISM method. By virtue of the algorithm, one can choose a longer timestep for renewing the solvent-induced force compared with that of the conformational update. To illustrate the present MD/3D-RISM simulation, we applied the method to a model of acetylacetone in aqueous solution. The multiple timestep algorithm succeeded in enhancing the computation speed by 3.4 times for this model case. Acetylacetone possesses an intramolecular hydrogen-bonding capability between the hydroxyl group and the carbonyl oxygen atom, and the molecule is significantly stabilized due to this hydrogen bond, especially in gas phase. The intramolecular hydrogen bond was kept intact during almost entire course of the MD simulation in gas phase, while in the aqueous solutions the bond is disrupted in a significant number of conformations. This result qualitatively agrees with the behavior on a free energy barrier lying upon the process for rotating a torsional degree of freedom of the hydroxyl group, where it is significantly reduced in aqueous solution by a cancellation between the electrostatic interaction and the solvation free energy.     

Monte Carlo sampling algorithm for searching a scale-transformed energy space of polypeptides

A Monte Carlo sampling algorithm for searching a scale-transformed conformational energy space of polypeptides is presented. This algorithm is based on the assumption that energy barriers can be overcome by a uniform sampling of the logarithmically transformed energy space. This algorithm is tested with Met-enkephalin. For comparison, the entropy sampling Monte Carlo (ESMC) simulation is performed. First, the global minimum is easily found by the optimization of a scale-transformed energy space. With a new Monte Carlo sampling, energy barriers of 3000 kcal/mol are frequently overcome, and low-energy conformations are sampled more efficiently than with ESMC simulations. Several thermodynamic quantities are calculated with good accuracy.     

Bennett's acceptance ratio and histogram analysis methods enhanced by umbrella sampling along a reaction coordinate in configurational space

Free energy perturbation, a method for computing the free energy difference between two states, is often combined with non-Boltzmann biased sampling techniques in order to accelerate the convergence of free energy calculations. Here we present a new extension of the Bennett acceptance ratio (BAR) method by combining it with umbrella sampling (US) along a reaction coordinate in configurational space. In this approach, which we call Bennett acceptance ratio with umbrella sampling (BAR-US), the conditional histogram of energy difference (a mapping of the 3N-dimensional configurational space via a reaction coordinate onto 1D energy difference space) is weighted for marginalization with the associated population density along a reaction coordinate computed by US. This procedure produces marginal histograms of energy difference, from forward and backward simulations, with higher overlap in energy difference space, rendering free energy difference estimations using BAR statistically more reliable. In addition to BAR-US, two histogram analysis methods, termed Bennett overlapping histograms with US (BOH-US) and Bennett-Hummer (linear) least square with US (BHLS-US), are employed as consistency and convergence checks for free energy difference estimation by BAR-US. The proposed methods (BAR-US, BOH-US, and BHLS-US) are applied to a 1-dimensional asymmetric model potential, as has been used previously to test free energy calculations from non-equilibrium processes. We then consider the more stringent test of a 1-dimensional strongly (but linearly) shifted harmonic oscillator, which exhibits no overlap between two states when sampled using unbiased Brownian dynamics. We find that the efficiency of the proposed methods is enhanced over the original Bennett's methods (BAR, BOH, and BHLS) through fast uniform sampling of energy difference space via US in configurational space. We apply the proposed methods to the calculation of the electrostatic contribution to the absolute solvation free energy (excess chemical potential) of water. We then address the controversial issue of ion selectivity in the K(+) ion channel, KcsA. We have calculated the relative binding affinity of K(+) over Na(+) within a binding site of the KcsA channel for which different, though adjacent, K(+) and Na(+) configurations exist, ideally suited to these US-enhanced methods. Our studies demonstrate that the significant improvements in free energy calculations obtained using the proposed methods can have serious consequences for elucidating biological mechanisms and for the interpretation of experimental data.     

Replica state exchange metadynamics for improving the convergence of free energy estimates

Metadynamics (MTD) is a powerful enhanced sampling method for systems with rugged energy landscapes. It constructs a bias potential in a predefined collective variable (CV) space to overcome barriers between metastable states. In bias‐exchange MTD (BE‐MTD), multiple replicas approximate the CV space by exchanging bias potentials (replica conditions) with the Metropolis–Hastings (MH) algorithm. We demonstrate that the replica‐exchange rates and the convergence of free energy estimates of BE‐MTD are improved by introducing the infinite swapping (IS) or the Suwa‐Todo (ST) algorithms. Conceptually, IS and ST perform transitions in a replica state space rather than exchanges in a replica condition space. To emphasize this, the proposed scheme is called the replica state exchange MTD (RSE‐MTD). Benchmarks were performed with alanine polypeptides in vacuum and water. For the systems tested in this work, there is no significant performance difference between IS and ST. © 2015 Wiley Periodicals, Inc.     

Enhanced conformational sampling method for proteins based on the TaBoo SeArch algorithm: Application to the folding of a mini‐protein,chignolin

The conformational samplings are indispensible for obtaining reliable canonical ensembles, which provide statistical averages of physical quantities such as free energies. However, the samplings of vast conformational space of biomacromolecules by conventional molecular dynamics (MD) simulations might be insufficient, due to their inadequate accessible time‐scales for investigating biological functions. Therefore, the development of methodologies for enhancing the conformational sampling of biomacromolecules still remains as a challenging issue in computational biology. To tackle this problem, we newly propose an efficient conformational search method, which is referred as TaBoo SeArch (TBSA) algorithm. In TBSA, an inverse energy histogram is used to select seeds for the conformational resampling so that states with high frequencies are inhibited, while states with low frequencies are efficiently sampled to explore the unvisited conformational space. As a demonstration, TBSA was applied to the folding of a mini‐protein, chignolin, and automatically sampled the native structure (C_α root mean square deviation < 1.0 Å) with nanosecond order computational costs started from a completely extended structure, although a long‐time 1‐µs normal MD simulation failed to sample the native structure. Furthermore, a multiscale free energy landscape method based on the conformational sampling of TBSA were quantitatively evaluated through free energy calculations with both implicit and explicit solvent models, which enable us to find several metastable states on the folding landscape. © 2015 Wiley Periodicals, Inc.     

Free energy of solvation from molecular dynamics simulation applying Voronoi-Delaunay triangulation to the cavity creation

The free energy of solvation for a large number of representative solutes in various solvents has been calculated from the polarizable continuum model coupled to molecular dynamics computer simulation. A new algorithm based on the Voronoi-Delaunay triangulation of atom-atom contact points between the solute and the solvent molecules is presented for the estimation of the solvent-accessible surface surrounding the solute. The volume of the inscribed cavity is used to rescale the cavitational contribution to the solvation free energy for each atom of the solute atom within scaled particle theory. The computation of the electrostatic free energy of solvation is performed using the Voronoi-Delaunay surface around the solute as the boundary for the polarizable continuum model. Additional short-range contributions to the solvation free energy are included directly from the solute-solvent force field for the van der Waals-type interactions. Calculated solvation free energies for neutral molecules dissolved in benzene, water, CCl4, and octanol are compared with experimental data. We found an excellent correlation between the experimental and computed free energies of solvation for all the solvents. In addition, the employed algorithm for the cavity creation by Voronoi-Delaunay triangulation is compared with the GEPOL algorithm and is shown to predict more accurate free energies of solvation, especially in solvents composed by molecules with nonspherical molecular shapes.     

Mapping potential energy surfaces

A recently proposed dynamical method [A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U.S.A. 99, 12562 (2002)] allows us to globally sample the free energy surface. This approach uses a coarse-grained non-Markovian dynamics to bias microscopic atomic trajectories. After a sufficiently long simulation time, the global free energy surface can be reconstructed from the non-Markovian dynamics. Here we apply this scheme to study the T=0 free energy surface, i.e., the potential energy surface in coarse-grained space. We show that the accuracy of the reconstructed potential energy surface can be dramatically improved by a simple postprocessing procedure with only minor computational overhead. We illustrate this approach by conducting conformational analysis on a small organic molecule, demonstrating its superiority over traditional unbiased approaches in sampling potential energy surfaces in coarse-grained space.     

Generating conformer ensembles using a multiobjective genetic algorithm

The task of generating a nonredundant set of low-energy conformations for small molecules is of fundamental importance for many molecular modeling and drug-design methodologies. Several approaches to conformer generation have been published. Exhaustive searches suffer from the exponential growth of the search space with increasing degrees of conformational freedom (number of rotatable bonds). Stochastic algorithms do not suffer as much from the exponential increase of search space and provide a good coverage of the energy minima. Here, the use of a multiobjective genetic algorithm in the generation of conformer ensembles is investigated. Distance geometry is used to generate an initial conformer, which is then subject to geometric modifications encoded by the individuals of the genetic algorithm. The geometric modifications apply to torsion angles about rotatable bonds, stereochemistry of double bonds and tetrahedral chiral centers, and ring conformations. The geometric diversity of the evolving conformer ensemble is preserved by a fitness-sharing mechanism based on the root-mean-square distance of the atomic coordinates. Molecular symmetry is taken into account in the distance calculation. The geometric modifications introduce strain into the structures. The strain is relaxed using an MMFF94-like force field in a postprocessing step that also removes conformational duplicates and structures whose strain energy remains above a predefined window from the minimum energy value found in the set. The implementation, called Balloon, is available free of charge on the Internet ( http://www.abo.fi/~mivainio/balloon/).     

Shape minimisation problems in liquid crystals

ABSTRACT

We consider a class of liquid crystal free-boundary problems for which both the equilibrium shape and internal configuration of a system must simultaneously be determined, for example in films with air–liquid or fluid–liquid crystal interfaces and elastomers. We develop a finite element algorithm to solve such problems with dynamic mesh control, achieved by supplementing the free energy with an auxiliary functional that promotes mesh quality and is minimised in the null space of the energy. We apply this algorithm to a flexible capacitor, as well as to determine the shape of liquid crystal tactoids as a function of the surface tension and elastic constants. These are compared with theoretical predictions and experimental observations of tactoids from the literature.     

One- and multidimensional conformational free energy simulations

A new thermodynamic integration approach to conformational free energy simulations is presented. The method is applicable both to one-dimensional cases (reaction coordinates) and multidimensional situations (free energy surfaces). Analysis of the properties of the thermodynamic integration algorithm is used to formulate methods of calculating multidimensional free energy gradients. The method is applied to calculate the free energy profile for rotation around the central C—C bond of n-butane in the gas and liquid phase and to generate maps of the 18-dimensional free energy gradient with respect to all nine ϕ and nine ψ dihedrals of the decaalanine and deca-α-methylalanine peptides in vacuum. For n-butane essentially no change in the gauche–trans equilibrium between the gas and liquid is predicted within the CHARMM explicit hydrogen model, with the thermodynamic integration, thermodynamic perturbation, and direct simulation methods yielding free energy profiles that are identical within errors. For the decapeptides the right-handed helical region of conformational space is investigated. For decaalanine a minimum on the free energy surface is found in the vicinity of (ϕ, ψ) = (-64.5°, -42.5°) in the α-helix region; no minimum exists for 3₁₀-helix-type conformers. For deca-α-methylalanine free energy minima corresponding to both the α-helix at ( - 55.5°, - 51.5°) and the 3₁₀-helix at ( - 54°, - 29°) are found; the α-helix state is favored by about 4 kcal/mol and the barrier for the concerted 3₁₀-helix → α-helix transition is about 3 kcal/mol. The α-methylation also considerably increases the rigidity of the α-helix with respect to deformations. The computational efficiency, ease of generalization to calculations of multidimensional gradients, and analytical capability due to component analysis of free energy differences make the method a novel, powerful tool to improve the basic understanding of conformational equilibria of flexible molecules in condensed phases. A related scheme for energy minimization in the presence of holonomic constraints is also presented, allowing generation of adiabatic energy surfaces in constrained systems. © 1996 by John Wiley & Sons, Inc.