On computational efficiency of the hybrid Monte Carlo method applied to the multicanonical ensemble

Abstract

In the present paper, computational efficiency of the hybrid Monte Carlo (HMC) method applied to the multicanonical ensemble is studied; the HMC is an equation of motion guided Monte Carlo method. As in the standard HMC for the canonical ensemble, the multicanonical HMC calculations with high acceptance ratio show better efficiency; about 60% acceptance yields the best performance for the system examined.     

Interaction of Cardiotoxin A5 with Membrane: Role of Conformational Heterogeneity and Hydrophobic Properties

The hypothesis that local conformational differences of the snake venom cardiotoxins (cytotoxins, CT) may play a significant role in their interaction with membrane was tested by molecular modeling of the behavior of the CT A5 from the venom of Naja atra in water and at the water–membrane interface. Two models of the CT A5 spatial structure are known: the first was obtained by X-ray analysis and the second, by NMR studies in solution. A molecular dynamics (MD) analysis demonstrated that loop II of the toxin has a fixed -like shape in water, which does not depend on its initial structure. An interaction of the experimentally derived (X-ray and NMR) conformations and MD simulated conformations of CT A5 with the lipid bilayer was studied by the Monte Carlo method using the previously developed model of the implicit membrane. It is found that: (1) unlike the previously studied CT2 from the venom of cobra Naja oxiana, CT A5 has only loops I and II bound to the membrane with the involvement of a lesser number of hydrophobic residues. (2) A long hydrophobic area is formed on the surface of CT A5 due to the -like shape of loop II and the arrangement of loop I in proximity to loop II. This hydrophobic area favors the toxin embedding into the lipid bilayer. (3) The toxin retains its conformation upon interaction with the membrane. (4) The CT A5 molecule has close values of the potential energy in the membrane and in aqueous environment, which suggests dynamic character of the binding. The results of the molecular modeling indicate a definite configuration of loops I and II and, consequently, a specific character of distribution of polar and apolar properties on the toxin surface, which turns out to be the most energetically favorable.     

Speeding up Monte Carlo molecular simulation by a non-conservative early rejection scheme

Monte Carlo (MC) molecular simulation describes fluid systems with rich information, and it is capable of predicting many fluid properties of engineering interest. In general, it is more accurate and representative than equations of state. On the other hand, it requires much more computational effort and simulation time. For that purpose, several techniques have been developed in order to speed up MC molecular simulations while preserving their precision. In particular, early rejection schemes are capable of reducing computational cost by reaching the rejection decision for the undesired MC trials at an earlier stage in comparison to the conventional scheme. In a recent work, we have introduced a ‘conservative’ early rejection scheme as a method to accelerate MC simulations while producing exactly the same results as the conventional algorithm. In this paper, we introduce a ‘non-conservative’ early rejection scheme, which is much faster than the conservative scheme, yet it preserves the precision of the method. The proposed scheme is tested for systems of structureless Lennard-Jones particles in both canonical and NVT-Gibbs ensembles. Numerical experiments were conducted at several thermodynamic conditions for different number of particles. Results show that at certain thermodynamic conditions, the non-conservative method is capable of doubling the speed of the MC molecular simulations in both canonical and NVT-Gibbs ensembles.     

Sampling efficiency in explicit and implicit membrane environments studied by peptide folding simulations

Efficient exploration of the conformational space of peptides embedded in biological membranes is vital to extract converged thermodynamic and kinetic data from simulation studies. Recently developed implicit membrane models promise vast increases in sampling efficiency compared to explicit membrane simulations, allowing for ab initio structure prediction and functional studies. In this study, a previously developed implicit membrane model, based on the generalized Born method, is compared to an explicit di‐palmitoyl‐phosphatidyl‐choline lipid bilayer and an octane slab membrane mimic. The complete folding process of a synthetic 16‐residue peptide is compared using these three setups. Since the comparison requires the entire folding pathway to be captured, individual simulations ranged up to 3 μs of MD. A quantitative sampling comparison using a wide range of performance metrics reveals that the implicit membrane model is at least 2 orders of magnitude more efficient than the simplest explicit setups. Proteins 2009. © 2008 Wiley‐Liss, Inc.     

Multiple time step diffusive Langevin dynamics for proteins

We present an algorithm for simulating the long time scale dynamics of proteins and other macromolecules. Our method applies the concept of multiple time step integration to the diffusive Langevin equation, in which short time scale dynamics are replaced by friction and noise. The macromolecular force field is represented at atomic resolution. Slow motions are modeled by constrained Langevin dynamics with very large time steps, while faster degrees of freedom are kept in local thermal equilibrium. In the limit of a sufficiently large molecule, our algorithm is shown to reduce the CPU time required by two orders of magnitude. We test the algorithm on two systems, alanine dipeptide and bovine pancreatic trypsin inhibitor (BPTI), and find that it accurately calculates a variety of equilibrium and dynamical properties. In the case of BPTI, the CPU time required is reduced by nearly a factor of 60 compared to a conventional, unconstrained Langevin simulation using the same force field. Proteins 30:215–227, 1998. © 1998 Wiley-Liss, Inc.     

New version of Monte Carlo expanded ensemble method for precise calculations of free energy difference

A new version of Monte Carlo (MC) expanded ensemble (EE) method is proposed for the calculations of free energy difference (FED) between two different systems with close values of the free energy. In order to check the method the FED between simple model systems (fluid of hard spheres and freely jointed polymer chain of hard spheres) was calculated. The free energy of the mentioned above systems was also calculated by a standard MC EE method in order to compare the results of two simulations. It was shown that the accuracy of a new algorithm is the same as of a standard one. At the same time new version of EE allows us to obtain FED between two systems having quite different structures, but similar free energies, during one simulation run.     

Gibbs Ensemble Calculations with an Equation of State: An Application to Vapor-Liquid Equilibria

Abstract

A new modification of the Gibbs ensemble Monte Carlo computer simulation method for fluid phase equilibria is described. The modification is based on a thermodynamic model for the vapor phase, and uses an equation of state to account for the weak interactions between the vapor phase molecules. Reductions in the computational time by 30–40% as compared to the original Gibbs ensemble method are obtained. The algorithm is applied to Lennard-Jones - (12,6) fluids and their mixtures and the results are in good agreement with results obtained from simulations using the full Gibbs ensemble method.     

Low-mass molecular dynamics simulation: A simple and generic technique to enhance configurational sampling

CLN025 is one of the smallest fast-folding proteins. Until now it has not been reported that CLN025 can autonomously fold to its native conformation in a classical, all-atom, and isothermal–isobaric molecular dynamics (MD) simulation. This article reports the autonomous and repeated folding of CLN025 from a fully extended backbone conformation to its native conformation in explicit solvent in multiple 500-ns MD simulations at 277 K and 1 atm with the first folding event occurring as early as 66.1 ns. These simulations were accomplished by using AMBER forcefield derivatives with atomic masses reduced by 10-fold on Apple Mac Pros. By contrast, no folding event was observed when the simulations were repeated using the original AMBER forcefields of FF12SB and FF14SB. The results demonstrate that low-mass MD simulation is a simple and generic technique to enhance configurational sampling. This technique may propel autonomous folding of a wide range of miniature proteins in classical, all-atom, and isothermal–isobaric MD simulations performed on commodity computers—an important step forward in quantitative biology.     

Application of statistical potentials to protein structure refinement from low resolution ab initio models

Recently ab initio protein structure prediction methods have advanced sufficiently so that they often assemble the correct low resolution structure of the protein. To enhance the speed of conformational search, many ab initio prediction programs adopt a reduced protein representation. However, for drug design purposes, better quality structures are probably needed. To achieve this refinement, it is natural to use a more detailed heavy atom representation. Here, as opposed to costly implicit or explicit solvent molecular dynamics simulations, knowledge-based heavy atom pair potentials were employed. By way of illustration, we tried to improve the quality of the predicted structures obtained from the ab initio prediction program TOUCHSTONE by three methods: local constraint refinement, reduced predicted tertiary contact refinement, and statistical pair potential guided molecular dynamics. Sixty-seven predicted structures from 30 small proteins (less than 150 residues in length) representing different structural classes (alpha, beta, alpha;/beta) were examined. In 33 cases, the root mean square deviation (RMSD) from native structures improved by more than 0.3 A; in 19 cases, the improvement was more than 0.5 A, and sometimes as large as 1 A. In only seven (four) cases did the refinement procedure increase the RMSD by more than 0.3 (0.5) A. For the remaining structures, the refinement procedures changed the structures by less than 0.3 A. While modest, the performance of the current refinement methods is better than the published refinement results obtained using standard molecular dynamics.     

Study of the Villin headpiece folding dynamics by combining coarse-grained Monte Carlo evolution and all-atom molecular dynamics

The folding mechanism of the Villin headpiece (HP36) is studied by means of a novel approach which entails an initial coarse-grained Monte Carlo (MC) scheme followed by all-atom molecular dynamics (MD) simulations in explicit solvent. The MC evolution occurs in a simplified free-energy landscape and allows an efficient selection of marginally-compact structures which are taken as viable initial conformations for the MD. The coarse-grained MC structural representation is connected to the one with atomic resolution through a "fine-graining" reconstruction algorithm. This two-stage strategy is used to select and follow the dynamics of seven different unrelated conformations of HP36. In a notable case the MD trajectory rapidly evolves towards the folded state, yielding a typical root-mean-square deviation (RMSD) of the core region of only 2.4 A from the closest NMR model (the typical RMSD over the whole structure being 4.0 A). The analysis of the various MC-MD trajectories provides valuable insight into the details of the folding and mis-folding mechanisms and particularly about the delicate influence of local and nonlocal interactions in steering the folding process.     

Gibbs ensemble Monte Carlo simulations of binary vapour–liquid–liquid equilibrium: application to n-hexane–water and ethane–ethanol systems

Gibbs ensemble Monte Carlo (GEMC) simulations in the isochoric–isothermal (NVT) ensemble were used to simulate vapour–liquid–liquid equilibrium (VLLE) for binary n-hexane–water and ethane–ethanol mixtures. The GEMC simulation of binary VLLE data proved to be extremely difficult and that is probably the reason why the open literature is so sparse with simulations for these types of systems. The results presented in this paper are to our knowledge the first successful binary three-phase GEMC simulations of non-idealised fluid systems. This paper also shows that the isobaric–isothermal (NPT) ensemble is unsuitable for the simulation of phase equilibria of binary three-phase systems.     

RASPA: molecular simulation software for adsorption and diffusion in flexible nanoporous materials

A new software package, RASPA, for simulating adsorption and diffusion of molecules in flexible nanoporous materials is presented. The code implements the latest state-of-the-art algorithms for molecular dynamics and Monte Carlo (MC) in various ensembles including symplectic/measure-preserving integrators, Ewald summation, configurational-bias MC, continuous fractional component MC, reactive MC and Baker's minimisation. We show example applications of RASPA in computing coexistence properties, adsorption isotherms for single and multiple components, self- and collective diffusivities, reaction systems and visualisation. The software is released under the GNU General Public License.     

A Monte Carlo method to account for sampling error in multi-species indicators

The usefulness of biodiversity indicators strongly increases if accompanied by measures of uncertainty. In the case of indicators that combine population indices of species, however, the inclusion of the uncertainty of the species indices has shown to be hard to realize, usually due to imperfections in monitoring programmes. Missing values and time series of different lengths preclude the use of analytical approaches, whereas bootstrapping across sites requires the raw abundance data on the site level, which may not always be available. Sometimes bootstrapping across species rather than sites is opted for, but this approach ignores the uncertainty attached to species indices. We developed a method to account for sampling error of species indices in the calculation of multi-species indicators based on Monte Carlo simulation of annual species indices. The construction of confidence intervals enables various trend assessments, like testing for linear or smooth trends, testing for changes between two time points, testing the significance of a suspected change-point and testing for differences between two multi-species indicators. Here, we compare our method with conventional methods and illustrate the benefits of our approach using Dutch breeding bird indicators.     

Conformational space exploration of Met- and Leu-enkephalin using the MOLS method, molecular dynamics, and Monte Carlo simulation--a comparative study

We report here a comparative study of the molecular conformational energy landscape generated using the mutually orthogonal Latin squares (MOLS) method, molecular dynamics (MD), and Monte Carlo (MC) simulation. The MOLS method, as described earlier from our laboratory, uses an experimental design technique to rapidly and exhaustively sample the low energy conformations of a molecule. MD and MC simulations have been used to perform similar tasks. In the comparison reported here, the three methods were applied to a pair of neuropeptides, namely Met- and Leu-enkephalin. A set of 1500 conformations of these enkephalins were generated using these methods with CHARMM22 force field, and the resulting samples were analyzed to determine the extent and nature of coverage of the conformational space. The results indicate that the MOLS method samples a larger number of possible conformations and identifies conformations closer to the experimental structures than the MD and MC simulations.     

An efficient conformational sampling method for homology modeling

The structural refinement of protein models is a challenging problem in protein structure prediction (Moult et al., Proteins 2003;53(Suppl 6):334-339). Most attempts to refine comparative models lead to degradation rather than improvement in model quality, so most current comparative modeling procedures omit the refinement step. However, it has been shown that even in the absence of alignment errors and using optimal templates, methods based on a single template have intrinsic limitations, and that refinement is needed to improve model accuracy. It is thought that failure of current methods originates on one hand from the inaccuracy of the effective free energy functions adopted, which do not represent properly the energetic balance in the native state, and on the other hand from the difficulty to sample the high dimensional and rugged free energy landscape of protein folding, in the search for the global minimum. Here, we address this second issue. We define the evolutionary and vibrational armonics subspace (EVA), a reduced sampling subspace that consists of a combination of evolutionarily favored directions, defined by the principal components of the structural variation within a homologous family, plus topologically favored directions, derived from the low frequency normal modes of the vibrational dynamics, up to 50 dimensions. This subspace is accurate enough so that the cores of most proteins can be represented within 1 A accuracy, and reduced enough so that Replica Exchange Monte Carlo (Hukushima and Nemoto, J Phys Soc Jpn 1996;65:1604-1608; Hukushima et al., Int J Mod Phys C: Phys Comput 1996;7:337-344; Mitsutake et al., J Chem Phys 2003;118:6664-6675; Mitsutake et al., J Chem Phys 2003;118:6676-6688) (REMC) can be applied. REMC is one of the best sampling methods currently available, but its applicability is restricted to spaces of small dimensionality. We show that the combination of the EVA subspace and REMC can essentially solve the optimization problem for backbone atoms in the reduced sampling subspace, even for rather rugged free energy landscapes. Applications and limitations of this methodology are finally discussed.     

All‐atom structural models of insulin binding to the insulin receptor in the presence of a tandem hormone‐binding element

Insulin regulates blood glucose levels in higher organisms by binding to and activating insulin receptor (IR), a constitutively homodimeric glycoprotein of the receptor tyrosine kinase (RTK) superfamily. Therapeutic efforts in treating diabetes have been significantly impeded by the absence of structural information on the activated form of the insulin/IR complex. Mutagenesis and photo‐crosslinking experiments and structural information on insulin and apo‐IR strongly suggest that the dual‐chain insulin molecule, unlike the related single‐chain insulin‐like growth factors, binds to IR in a very different conformation than what is displayed in storage forms of the hormone. In particular, hydrophobic residues buried in the core of the folded insulin molecule engage the receptor. There is also the possibility of plasticity in the receptor structure based on these data, which may in part be due to rearrangement of the so‐called CT‐peptide, a tandem hormone‐binding element of IR. These possibilities provide opportunity for large‐scale molecular modeling to contribute to our understanding of this system. Using various atomistic simulation approaches, we have constructed all‐atom structural models of hormone/receptor complexes in the presence of CT in its crystallographic position and a thermodynamically favorable displaced position. In the “displaced‐CT” complex, many more insulin–receptor contacts suggested by experiments are satisfied, and our simulations also suggest that R‐insulin potentially represents the receptor‐bound form of hormone. The results presented in this work have further implications for the design of receptor‐specific agonists/antagonists. Proteins 2013; © 2012 Wiley Periodicals, Inc.     

The Lennard-Jones Fluid Revisited: Computer Simulation Results

Abstract

Pseudoexperimental data of high accuracy on the pressure and the internal energy of the Lennard-Jones fluid have been generated both by the Monte Carlo and molecular dynamics methods for five subcritical and three supercritical isotherms. Values of the chemical potential of the Lennard-Jones fluid computed by a new version of the gradual insertion particle method for two isotherms up to very high densities are also reported and discussed, and compared with existing data.     

Analytical expressions on harmonic oscillators for variational path integrals

In this study, we derive analytical expressions for one-dimensional harmonic oscillators for variational path integrals (VPIs). A Gaussian-type trial wavefunction is adopted. Total and potential energies of the system are analytically expressed both for the continuous time and an approximate discretised VPIs. Obtained expressions are numerically verified using molecular dynamics calculations. Convergence properties regarding the projection time and the Trotter number are discussed.     

Gas phase infrared spectra via the phase integration quasi-classical method

We apply the recently developed phase integration method (PIM) (Monteferrante et al. Mol Phys. 2011;109:3015–3027) to the calculation of infrared spectra of gas phase molecules. The PIM combines a generalised Monte Carlo sampling of the exact thermal density of the system with classical molecular dynamics to obtain approximate time quantum correlation functions. To describe the molecules, we adopt very simple analytical potentials that have, however, proved interesting, and surprisingly challenging, benchmarks for approximate quantum dynamical schemes. We show that, in contrast with two other commonly applied methods, our spectra do not exhibit spurious features or unphysical shifts depending on the temperature. Identifying the positions of the peaks requires only a few tens of trajectories, while an accurate evaluation of the relative intensities of the peaks is computationally more demanding.