Molecular dynamics simulation in the grand canonical ensemble

An extended system Hamiltonian is proposed to perform molecular dynamics (MD) simulation in the grand canonical ensemble. The Hamiltonian is similar to the one proposed by Lynch and Pettitt (Lynch and Pettitt, J Chem Phys 1997, 107, 8594), which consists of the kinetic and potential energies for real and fractional particles as well as the kinetic and potential energy terms for material and heat reservoirs interacting with the system. We perform a nonlinear scaling of the potential energy parameters of the fractional particle, as well as its mass to vary the number of particles dynamically. On the basis of the equations of motion derived from this Hamiltonian, an algorithm has been proposed for MD simulation at constant chemical potential. The algorithm has been tested for the ideal gas, for the Lennard-Jones fluid over a wide range of temperatures and densities, and for water. The results for the low-density Lennard-Jones fluid are compared with the predictions from a truncated virial equation of state. In the case of the dense Lennard-Jones fluid and water our predicted results are compared with the results reported using other available methods for the calculation of the chemical potential. The method is also applied to the case of vapor-liquid coexistence point predictions.     

Dissipative particle dynamics simulations in the grand canonical ensemble: applications to polymer brushes.

We have used the dissipative particle dynamics (DPD) method in the grand canonical ensemble to study the compression of grafted polymer brushes in good solvent conditions. The force-distance profiles calculated from DPD simulations in the grand canonical ensemble are in very good agreement with the self-consistent field (SCF) theoretical models and with experimental results for two polystyrene brush layers grafted onto mica surfaces in toluene.     

On‐the‐fly reconstruction of free‐energy profiles using logarithmic mean‐force dynamics

Mean‐force dynamics (MFD), which is a fictitious dynamics for a set of collective variables on a potential of mean‐force, is a powerful algorithm to efficiently explore free‐energy landscapes. Recently, we have introduced logarithmic MFD (LogMFD) (Morishita et al., Phys. Rev. E 2012, 85, 066702) which overcomes difficulties encounterd in free‐energy calculations using standard approaches such as thermodynamic integration. Here, we present a guide to implementing LogMFD calculations paying attention to the practical issues in choosing the parameters in LogMFD. A primary focus is given to the effect of the parameters on the accuracy of the reconstructed free‐energy profiles. A recipe for reducing the errors due to energy dissipation is presented. We also demonstrate that multidimensional free‐energy landscapes can be reconstructed on‐the‐fly using LogMFD, which cannot be accomplished using any other free‐energy calculation techniques. © 2013 Wiley Periodicals, Inc.     

Detailed considerations for a balanced and broadly applicable force field: a study of substituted benzenes modeled with OPLS-AA

Modern classical force fields have been traditionally parameterized by attempting to maximize agreement to any number of experimental and/or quantum mechanical target properties. As these force fields are pushed towards obtaining quantitative estimates of often subtle energetic differences, stringent and consistent parameterization criteria, particularly in regard to charge distributions, are required to ensure that systematic errors cancel, that parameters are transferable between molecules, and that performance does not significantly deteriorate when using more approximate methods, such as with continuum solvent models. Relative free energies of hydration are presented here for 40 mono- and disubstituted benzenes modeled with the OPLS-AA force field; heats of vaporization and pure liquid densities at standard conditions are presented when experimental data is available. Overall agreement between OPLS-AA and experiment is remarkable (average error = 0.5 kcal/mol for DeltaDeltaG(hydration), 1.0 kcal/mol for DeltaH(vap) (0), 0.02 g/mL for densities), yet several functional groups are identified as having consistent and correctable errors (alkyl-, nitro-, and thiobenzenes). Relative free energies of hydration obtained with rigorous free energy perturbations using explicit solvent are also compared with energies from minimizations using a generalized Born model (GB). There is high correlation between these estimates (R = 0.99), and as demonstrated here, reparameterization of the aforementioned groups can be guided with rapid GB calculations.     

Theory for trivial trajectory parallelization of multicanonical molecular dynamics and application to a polypeptide in water

Trivial trajectory parallelization of multicanonical molecular dynamics (TTP-McMD) explores the conformational space of a biological system with multiple short runs of McMD starting from various initial structures. This method simply connects (i.e., trivially parallelizes) the short trajectories and generates a long trajectory. First, we theoretically prove that the simple trajectory connection satisfies a detailed balance automatically. Thus, the resultant long trajectory is regarded as a single multicanonical trajectory. Second, we applied TTP-McMD to an alanine decapeptide with an all-atom model in explicit water to compute a free-energy landscape. The theory imposes two requirements on the multiple trajectories. We have demonstrated that TTP-McMD naturally satisfies the requirements. The TTP-McMD produces the free-energy landscape considerably faster than a single-run McMD does. We quantitatively showed that the accuracy of the computed landscape increases with increasing the number of multiple runs. Generally, the free-energy landscape of a large biological system is unknown a priori. The current method is suitable for conformational sampling of such a large system to reduce the waiting time to obtain a canonical ensemble statistically reliable.     

层柱状微孔材料吸附存储天然气的Monte Carlo模拟

采用巨正则系综MonteCarlo方法模拟了天然气中主要成分甲烷在层柱状微孔材料中T=300K下的吸附存储,在模拟中层柱状微孔采用Yi等人建立的柱子均匀分布在两炭孔墙之间的模型来表征。甲烷分子采用Lennard-Jones球型分子模型,炭孔墙采用Steele的10-4-3模型,对孔宽为1.36nm的层柱微孔,模拟了四个不同孔率的层柱材料吸附甲烷的情形。得到了孔中流体的局部密度分布以及吸附等温线,对比不同孔率下甲烷的吸附量,得到了此情形吸附甲烷的较佳孔率为0.94。     

Stockmayer流体在活性炭孔中的吸附的分子模拟

应用巨正则系综monteCarlo方法模拟Stockmayer流体[以一氯二氟甲烷(R22)为代表]在活性炭孔中的吸附。模拟中R22分子采用等效Stockmayer势能模型,狭缝碳孔墙采用10-4-3模型。通过模拟得到了最佳孔径,并在最佳孔径下,针对不同的主体压力及活性基团密度,得到了吸附等温线、孔中流体的局部密度分布图和较为直观的孔内流体分子的瞬时构象,分析了吸附等温线的特征及孔内流体的吸附结构,认为在0.0,1.0sites/nm^2的活性基团密度下的碳孔内分别发生物理及化学吸附,并确定了最佳操作压力,为工业设计合适的催化剂提供依据。     

MCM-22型分子筛中苯分子吸附行为的蒙特卡罗模 拟研究

用巨正则统计系综蒙特卡罗模拟方法研究了纯硅MCM-22型分子筛(ITQ-1)中苯分子的吸附行为。结果表明苯分子在ITQ-1型分子筛中主要存在4个吸附位点。从苯分子粒子分布云图上可以看到苯分子的扩散和吸附主要在12元环超笼内发生。在苯分子的扩散过程中,S2位置附近的苯分子分布较为集中,而S3和S4附近的苯分子分布则较为离散。苯分子通过10元环窗口的运动路径势能面的计算结果表明,苯分子在12元环超笼内可以较为自由迁移,而通过10元环窗口从一个超笼扩散到附近的超笼时则需要较高的激发能量,这个能量大约为100kJ/mol。     

Optimal designs for pairwise calculation: An application to free energy perturbation in minimizing prediction variability

Pairwise-based methods such as the free energy perturbation (FEP) method have been widely deployed to compute the binding free energy differences between two similar host–guest complexes. The calculated pairwise free energy difference is either directly adopted or transformed to absolute binding free energy for molecule rank ordering. We investigated, through both analytic derivations and simulations, how the selection of pairs in the experiment could impact the overall prediction precision. Our studies showed that (1) the estimated absolute binding free energy () derived from calculated pairwise differences (ΔΔG) through weighted least squares fitting is more precise in prediction than the pairwise difference values when the number of pairs is more than the number of ligands and (2) prediction precision is influenced by both the total number of pairs and the specifically selected pairs, the latter being critically important when the number of calculated pairs is limited. Furthermore, we applied optimal experimental design in pair selection and found that the optimally selected pairs can outperform randomly selected pairs in prediction precision. In an illustrative example, we showed that, upon weighing ligand structure similarity into design optimization, the weighted optimal designs are more efficient than the literature reported designs. This work provides a new approach to assess retrospective pairwise-based prediction results, and a method to design new prospective pairwise-based experiments for molecular lead optimization. © 2019 Wiley Periodicals, Inc.     

Potential of mean force of water–proton bath and molecular dynamic simulation of proteins at constant pH

An advanced implicit solvent model of water–proton bath for protein simulations at constant pH is presented. The implicit water–proton bath model approximates the potential of mean force of a protein in water solvent in a presence of hydrogen ions. Accurate and fast computational implementation of the implicit water–proton bath model is developed using the continuum electrostatic Poisson equation model for calculation of ionization equilibrium and the corrected MSR6 generalized Born model for calculation of the electrostatic atom–atom interactions and forces. Molecular dynamics (MD) method for protein simulation in the potential of mean force of water–proton bath is developed and tested on three proteins. The model allows to run MD simulations of proteins at constant pH, to calculate pH‐dependent properties and free energies of protein conformations. The obtained results indicate that the developed implicit model of water–proton bath provides an efficient way to study thermodynamics of biomolecular systems as a function of pH, pH‐dependent ionization‐conformation coupling, and proton transfer events. © 2012 Wiley Periodicals, Inc.     

Effects of system net charge and electrostatic truncation on all‐atom constant pH molecular dynamics

Constant pH molecular dynamics offers a means to rigorously study the effects of solution pH on dynamical processes. Here, we address two critical questions arising from the most recent developments of the all‐atom continuous constant pH molecular dynamics (CpHMD) method: (1) What is the effect of spatial electrostatic truncation on the sampling of protonation states? (2) Is the enforcement of electrical neutrality necessary for constant pH simulations? We first examined how the generalized reaction field and force‐shifting schemes modify the electrostatic forces on the titration coordinates. Free energy simulations of model compounds were then carried out to delineate the errors in the deprotonation free energy and salt‐bridge stability due to electrostatic truncation and system net charge. Finally, CpHMD titration of a mini‐protein HP36 was used to understand the manifestation of the two types of errors in the calculated pK_a values. The major finding is that enforcing charge neutrality under all pH conditions and at all time via cotitrating ions significantly improves the accuracy of protonation‐state sampling. We suggest that such finding is also relevant for simulations with particle mesh Ewald, considering the known artifacts due to charge‐compensating background plasma. © 2014 Wiley Periodicals, Inc.     

Adaptive lambda square dynamics simulation: An efficient conformational sampling method for biomolecules

A novel, efficient sampling method for biomolecules is proposed. The partial multicanonical molecular dynamics (McMD) was recently developed as a method that improved generalized ensemble (GE) methods to focus sampling only on a part of a system (GEPS); however, it was not tested well. We found that partial McMD did not work well for polylysine decapeptide and gave significantly worse sampling efficiency than a conventional GE. Herein, we elucidate the fundamental reason for this and propose a novel GEPS, adaptive lambda square dynamics (ALSD), which can resolve the problem faced when using partial McMD. We demonstrate that ALSD greatly increases the sampling efficiency over a conventional GE. We believe that ALSD is an effective method and is applicable to the conformational sampling of larger and more complicated biomolecule systems. © 2013 Wiley Periodicals, Inc.     

自洽平均场正则系综下计算嵌段共聚物/均聚物共混相图

在自洽平均场中计算聚合物宏观相分离体系时,需要将正则系综与巨正则系综结合使用。 通过将正则系综与巨正则系综之间的变量进行转化,只在正则系综下计算即可得到巨正则系综下的相应变量的值,在很大程度上减少了计算量。 本文利用这种简化方法计算了A-b-B两嵌段共聚物与均聚物A在不同均聚物聚合度下随着均聚物含量变化的相图,其结果与巨正则系综下的计算结果相同。 该结果表明,在嵌段共聚物与均聚物的共混体系中,增加嵌段共聚物组成fA或者减小均聚物的聚合度,将有效阻止体系发生宏观相分离。     

Surveying implicit solvent models for estimating small molecule absolute hydration free energies

Implicit solvent models are powerful tools in accounting for the aqueous environment at a fraction of the computational expense of explicit solvent representations. Here, we compare the ability of common implicit solvent models (TC, OBC, OBC2, GBMV, GBMV2, GBSW, GBSW/MS, GBSW/MS2 and FACTS) to reproduce experimental absolute hydration free energies for a series of 499 small neutral molecules that are modeled using AMBER/GAFF parameters and AM1-BCC charges. Given optimized surface tension coefficients for scaling the surface area term in the nonpolar contribution, most implicit solvent models demonstrate reasonable agreement with extensive explicit solvent simulations (average difference 1.0-1.7 kcal/mol and R(2)=0.81-0.91) and with experimental hydration free energies (average unsigned errors=1.1-1.4 kcal/mol and R(2)=0.66-0.81). Chemical classes of compounds are identified that need further optimization of their ligand force field parameters and others that require improvement in the physical parameters of the implicit solvent models themselves. More sophisticated nonpolar models are also likely necessary to more effectively represent the underlying physics of solvation and take the quality of hydration free energies estimated from implicit solvent models to the next level.     

Free energy perturbation approach to the critical assessment of selective cyclooxygenase-2 inhibitors

The discovery of selective cyclooxygenase-2 (COX-2) inhibitors represents a major achievement of the efforts over the past few decades to develop therapeutic treatments for inflammation. To gain insights into designing new COX-2-selective inhibitors, we address the energetic and structural basis for the selective inhibition of COX isozymes by means of a combined computational protocol involving docking experiment, force field design for the heme prothetic group, and free energy perturbation (FEP) simulation. We consider both COX-2- and COX-1-selective inhibitors taking the V523I mutant of COX-2 to be a relevant structural model for COX-1 as confirmed by a variety of experimental and theoretical evidences. For all COX-2-selective inhibitors under consideration, we find that free energies of binding become less favorable as the receptor changes from COX-2 to COX-1, due to the weakening and/or loss of hydrogen bond and hydrophobic interactions that stabilize the inhibitors in the COX-2 active site. On the other hand, COX-1-selective oxicam inhibitors gain extra stabilization energy with the change of residue 523 from valine to isoleucine because of the formations of new hydrogen bonds in the enzyme-inhibitor complexes. The utility of the combined computational approach, as a valuable tool for in silico screening of COX-2-selective inhibitors, is further exemplified by identifying the physicochemical origins of the enantiospecific selective inhibition of COX-2 by -substituted indomethacin ethanolamide inhibitors.     

In silico analysis of carbohydrate-binding pockets in the lectin genes from various species of Canavalia

Legumes are endowed with an opulent class of proteins called lectins that can detect tenuous variations in carbohydrate structures and bind them reversibly with high affinity and specificity. The genus Canavalia, in the family of Leguminosae, is considered to be an affluent source of lectin. An effort has been made to analyse the sequences encoded by the lectin gene and its carbohydrate binding pockets from three species of Canavalia, including C. virosa, C. rosea, and C. pubescens. Crude seed extract showed highest haemagglutination titer against buffalo RBCs and has high affinity to mannose and trehalose. Amplification of the lectin gene by gene-specific primers showed the presence of an 870 bp amplicon. Physicochemical characterization using various bioinformatic tools showed that the isoelectric point was below 7, suggesting that lectin molecules were acidic. A high aliphatic index and high instability index were observed, which indicated that lectin molecules were stable towards a wide range of temperatures. The occurrence of N-glycosylation sites at two sites was also identified in all three species. Prediction of secondary structure showed that approximately 59.05 %, 56.76 % and 54.88 % of the elements were random coils in the case of C. virosa, C. pubescens and C. rosea, respectively. Comparative modelling of the proteins and docking of hypothetical models with sugar moieties that inhibited the agglutination activity suggested that asparagine, serine, alanine, valine, tyrosine and threonine were the major residues involved in hydrogen bonding and other stacking interactions. This can further provide insights on its prospective antibiosis property.     

Analysis of low-abundance proteins using the proteomic reactor with pH fractionation

We developed a new method consisting of the proteomic reactor coupled with step pH fractionation for the analysis of low-abundance proteins from minute amount of sample. These new reactors were implemented using both SAX and SCX materials. The pH fractions from the SAX reactor provided higher peptide and protein identification than SCX reactor and conventional solution digestion. Interestingly, the physical characteristics (pI, molecular weight, missed cleavage site and grand average hydrophobicity (GRAVY) index, and number of acid and basic amino acid) of the peptides obtained from the SAX and SCX proteomic reactors are drastically different. Furthermore, nearly half of the peptides observed from the pH fractionations from the SAX reactor are of low abundance while only 22% low-abundance proteins are observed with conventional in-solution digestion following 2D LC-MS/MS analysis.     

Monte carlo simulations of proteins at constant pH with generalized born solvent,flexible sidechains,and an effective dielectric boundary

Titratable residues determine the acid/base behavior of proteins, strongly influencing their function; in addition, proton binding is a valuable reporter on electrostatic interactions. We describe a method for pK_a calculations, using constant‐pH Monte Carlo (MC) simulations to explore the space of sidechain conformations and protonation states, with an efficient and accurate generalized Born model (GB) for the solvent effects. To overcome the many‐body dependency of the GB model, we use a “Native Environment” approximation, whose accuracy is shown to be good. It allows the precalculation and storage of interactions between all sidechain pairs, a strategy borrowed from computational protein design, which makes the MC simulations themselves very fast. The method is tested for 12 proteins and 167 titratable sidechains. It gives an rms error of 1.1 pH units, similar to the trivial “Null” model. The only adjustable parameter is the protein dielectric constant. The best accuracy is achieved for values between 4 and 8, a range that is physically plausible for a protein interior. For sidechains with large pK_a shifts, ≥2, the rms error is 1.6, compared to 2.5 with the Null model and 1.5 with the empirical PROPKA method. © 2013 Wiley Periodicals, Inc.