Development and use of an atomistic CHARMM‐based forcefield for peptoid simulation

Peptoids are positional isomers of peptides: peptoid sidechains are attached to backbone nitrogens rather than α‐carbons. Peptoids constitute a class of sequence‐specific polymers resistant to biological degradation and potentially as diverse, structurally and functionally, as proteins. While molecular simulation of proteins is commonplace, relatively few tools are available for peptoid simulation. Here, we present a first‐generation atomistic forcefield for peptoids. Our forcefield is based on the peptide forcefield CHARMM22, with key parameters tuned to match both experimental data and quantum mechanical calculations for two model peptoids (dimethylacetamide and a sarcosine dipeptoid). We used this forcefield to demonstrate that solvation of a dipeptoid substantially modifies the conformations it can access. We also simulated a crystal structure of a peptoid homotrimer, H‐(N‐2‐phenylethyl glycine)₃‐OH, and we show that experimentally observed structural and dynamical features of the crystal are accurately described by our forcefield. The forcefield presented here provides a starting point for future development of peptoid‐specific simulation methods within CHARMM. © 2013 Wiley Periodicals, Inc.     

Parameters of monovalent ions in the AMBER-99 forcefield: assessment of inaccuracies and proposed improvements

The monovalent ion parameters used by the AMBER-99 forcefield are shown to exhibit physically inaccurate behavior in molecular dynamics simulations of strong 1:1 electrolytes. These errors arise from an ad hoc adaptation of Aqvist's cation parameters. The result is the rapid formation of large, unphysical clusters at concentrations that are well below solubility limits. The observed unphysical behavior poses a serious challenge for simulating ions around highly charged polymers such as nucleic acids. In this communication, we explain the source of this unphysical behavior. To facilitate the continued use of the popular AMBER parameters, we prescribe a simple fix whereby Aqvist's cations and anions are used in conjunction with the AMBER forcefield for nucleic acids. A preliminary test of this strategy suggests that the proposed fix is reasonable and is likely to be generalizable for simulating diffuse and specific ion binding to nucleic acids.     

The Importance of Relative Reaction Rates in the Optimization of Detailed Kinetic Models

Numerous mathematical tools intended to adjust rate constants employed in complex detailed kinetic models to make them consistent with multiple sets of experimental data have been reported in the literature. Application of such model optimization methods typically begins with the assignment of uncertainties in the absolute rate constants in a starting model, followed by variation of the rate constants within these uncertainty bounds to tune rate parameters to match model outputs to experimental observations. The present work examines the impact of including information on relative reaction rates in the optimization strategy, which is not typically done in current implementations. It is shown that where such rate constant data are available, the available parameter space changes dramatically due to the correlations inherent in such measurements. Relative rate constants are typically measured with greater relative accuracy than corresponding absolute rate constant measurements. This greater accuracy further reduces the available parameter space, which significantly affects the uncertainty in the model outcomes as a result of kinetic parameter uncertainties. We demonstrate this effect by considering a simple example case emulating an ignition event and show that use of relative rate measurements leads to a significantly smaller uncertainty in the output ignition delay time in comparison with results based on absolute measurements. This is true even though the same range of absolute rate constants is sampled in each case. Implications of the results with respect to the maintenance of physically realistic kinetics in optimized models are discussed, and suggestions are made for the path forward in the refinement of detailed kinetic models.     

Development of glycyl radical parameters for the OPLS-AA/L force field

On the basis of quantum chemical calculations C(alpha)-glycyl radical parameters have been developed for the OPLS-AA/L force field. The molecular mechanics hypersurface was fitted to the calculated quantum chemical surface by minimizing their molecular mechanics parameter dependent sum-of-squares deviations. To do this, a computer program in which the molecular mechanics energy derivatives with respect to the parameters were calculated analytically was developed, implementing the general method of Lifson and Warshel (J Chem Phys 1968, 49, 5116) for force field parameter optimization. This program, in principle, can determine the optimal parameter set in one calculation if enough representative value points on the quantum chemical potential energy surface are available and there is no linear dependency between the parameters. Some of the parameters in quantum calculations, including several new torsion types around a bond as well as angle parameters at a new central atom type, are not completely separable. Consequently, some restrictions and/or presumptions were necessary during parameter optimization. The relative OPLS-AA energies reproduced those calculated quantum chemically almost perfectly.     

Dynamic lattice searching methods for optimization of clusters

Global optimization of clusters is a subject of intense interest in computational chemistry. Especially for large clusters, locating the global minima is a challenging problem. Two strategies are generally used for the problem, i.e., the stochastic optimization and the static modeling strategy. The former is known as unbiased global optimization method, while the latter is more efficient but biased. This review describes the development of a dynamic lattice searching (DLS) approach. In DLS, the lattices are constructed dynamically and optimization is achieved by searching these lattices. Therefore, DLS possesses the characteristics of both the stochastic and static methods. With the aim of improving the efficiency of DLS for optimization of large clusters, several variants of the method have been developed. The results show that DLS methods may be promising tools for fast modeling of large clusters. With this review, greater interests are expected for global optimization of atomic or molecular clusters.     

Global optimization for parameter estimation of differential-algebraic systems

The estimation of parameters in semi-empirical models is essential in numerous areas of engineering and applied science. In many cases, these models are described by a set of ordinary-differential equations or by a set of differential-algebraic equations. Due to the presence of non-convexities of functions participating in these equations, current gradient-based optimization methods can guarantee only locally optimal solutions. This deficiency can have a marked impact on the operation of chemical processes from the economical, environmental and safety points of view and it thus motivates the development of global optimization algorithms. This paper presents a global optimization method which guarantees ɛ-convergence to the global solution. The approach consists in the transformation of the dynamic optimization problem into a nonlinear programming problem (NLP) using the method of orthogonal collocation on finite elements. Rigorous convex underestimators of the nonconvex NLP problem are employed within the spatial branch-and-bound method and solved to global optimality. The proposed method was applied to two example problems dealing with parameter estimation from time series data.     

Solvation parameters. Part 5: physicochemical interpretation of experimental solvent values for stationary phases of gas-liquid chromatography

It has been demonstrated for a long time that in the particular case of gas-liquid chromatography (GLC), a linear free energy relationship (LFER) of five terms can be established, each term including a parameter of solute and a parameter of solvent. The nature of some of these parameters has been quite clearly identified, even if not always well predicted from the molecular structure. First of all, the five solute parameters: two involved in the hydrogen bonding and three in the Van der Waals forces; secondly, the two solvent parameters involved in hydrogen bonding. It was remaining an uncertainty concerning the nature of the solvent parameters named D, W and E, respectively associated with the solute parameters of dispersion, orientation and induction/polarizability. This uncertainty has been solved using experimental chromatographic data of McReynolds (56 phases) and of the Kováts group (11 phases). The parameter W appears as of polar nature strictly speaking. The parameters D and E can be expressed by two opposite bilinear functions of 1/V (inverse of molecular volume) and PSA/V (ratio of the polar surface area over the molecular volume). These results are in agreement with previous studies limited to alkanes by the Kováts group.     

SODOCK: swarm optimization for highly flexible protein-ligand docking

Protein-ligand docking can be formulated as a parameter optimization problem associated with an accurate scoring function, which aims to identify the translation, orientation, and conformation of a docked ligand with the lowest energy. The parameter optimization problem for highly flexible ligands with many rotatable bonds is more difficult than that for less flexible ligands using genetic algorithm (GA)-based approaches, due to the large numbers of parameters and high correlations among these parameters. This investigation presents a novel optimization algorithm SODOCK based on particle swarm optimization (PSO) for solving flexible protein-ligand docking problems. To improve efficiency and robustness of PSO, an efficient local search strategy is incorporated into SODOCK. The implementation of SODOCK adopts the environment and energy function of AutoDock 3.05. Computer simulation results reveal that SODOCK is superior to the Lamarckian genetic algorithm (LGA) of AutoDock, in terms of convergence performance, robustness, and obtained energy, especially for highly flexible ligands. The results also reveal that PSO is more suitable than the conventional GA in dealing with flexible docking problems with high correlations among parameters. This investigation also compared SODOCK with four state-of-the-art docking methods, namely GOLD 1.2, DOCK 4.0, FlexX 1.8, and LGA of AutoDock 3.05. SODOCK obtained the smallest RMSD in 19 of 37 cases. The average 2.29 A of the 37 RMSD values of SODOCK was better than those of other docking programs, which were all above 3.0 A.     

Configuration-specific kinetic theory applied to the elastic collisions of hard spherical molecules

Classical trajectory simulations can be used to glean a wealth of information on the geometric details of gas-phase molecular collision events for which the standard theoretical treatment lacks the ability to predict. For instance, the standard treatment gives no information on configuration-specific collision parameters. A configuration-specific parameter is defined here as the average value for a collision parameter that is exclusive to either an ensemble of front-end or an ensemble of rear-end molecular collisions. This paper presents statistical results of simulation "measurements" on several configuration-specific parameters, including the configuration-specific collision frequencies. The simulations use single-component systems of hard spherical molecules confined within a spherical boundary. To complement the simulation effort, a systematic mathematical analysis for the configuration-specific parameters is presented. This analysis uses the Maxwell-Boltzmann distribution of molecular speeds as usual, but exploits the distinction between front-end and rear-end collision space, and uses the line-of-centers speed rather than the relative speed. The configuration-specific expressions derived from this analysis are in very good agreement with the simulation measurements for every molecular collision parameter studied in this work.     

An improved nucleic acid parameter set for the GROMOS force field

Over the past decades, the GROMOS force field for biomolecular simulation has primarily been developed for performing molecular dynamics (MD) simulations of polypeptides and, to a lesser extent, sugars. When applied to DNA, the 43A1 and 45A3 parameter sets of the years 1996 and 2001 produced rather flexible double-helical structures, in which the Watson-Crick hydrogen-bonding content was more limited than expected. To improve on the currently available parameter sets, the nucleotide backbone torsional-angle parameters and the charge distribution of the nucleotide bases are reconsidered based on quantum-chemical data. The new 45A4 parameter set resulting from this refinement appears to perform well in terms of reproducing solution NMR data and canonical hydrogen bonding. The deviation between simulated and experimental observables is now of the same order of magnitude as the uncertainty in the experimental values themselves.     

Systematic comparison of empirical forcefields for molecular dynamic simulation of insulin

The use of atomistic simulation methodologies based on empirical forcefields has enhanced our understanding of many physical processes governing protein structure and dynamics. However, the forcefields used in classical modeling studies are often designed for a particular class of proteins and rely on continuous improvement and validation by comparison of simulations with experimental data. We present a comprehensive comparison of five popular forcefields for simulating insulin. The effect of each forcefield on the conformational evolution and structural properties of the peptide is analyzed in detail and compared with available experimental results. In this study we observed that different forcefields favor different structural trends. However, the all-atom forcefield CHARMM27 and the united-atom forcefield GROMOS 43A1 delivered the best representation of the experimentally observed dynamic behavior of chain B of insulin.     

Regression of binary interaction parameters for thermodynamic models using an inside-variance estimation method (IVEM)

An inside-variance estimation method (IVEM) for binary interaction parameter regression in thermodynamic models is proposed. This maximum likelihood method involves the re-computation of the variance for each iteration of the optimization procedure, automatically re-weighting the objective function. Most of the maximum likelihood approaches currently used to regress the parameters of thermodynamic models fix the variances, converting the problem into a traditional weighted least squares minimization. However, such approaches lead to residual variances (between measured and calculated values) that are inconsistent with the fixed variances and, thus, do not necessarily produce optimum parameters for prediction purposes. The new method (IVEM) substantially improves fluid phase equilibria predictions (as shown by the examples presented) by maintaining consistency between the residual variances and the variance used in the objective function. This results in better parameter estimation and to a direct measure of the uncertainty in the model prediction.     

Optimization of parameters for molecular dynamics simulation using smooth particle-mesh Ewald in GROMACS 4.5

Based on our critique of requirements for performing an efficient molecular dynamics simulation with the particle-mesh Ewald (PME) implementation in GROMACS 4.5, we present a computational tool to enable the discovery of parameters that produce a given accuracy in the PME approximation of the full electrostatics. Calculations on two parallel computers with different processor and communication structures showed that a given accuracy can be attained over a range of parameter space, and that the attributes of the hardware and simulation system control which parameter sets are optimal. This information can be used to find the fastest available PME parameter sets that achieve a given accuracy. We hope that this tool will stimulate future work to assess the impact of the quality of the PME approximation on simulation outcomes, particularly with regard to the trade-off between cost and scientific reliability in biomolecular applications.     

A comparative study of particle swarm optimization and its variants for phase stability and equilibrium calculations in multicomponent reactive and non-reactive systems

Particle swarm optimization is a novel evolutionary stochastic global optimization method that has gained popularity in the chemical engineering community. This optimization strategy has been successfully used for several applications including thermodynamic calculations. To the best of our knowledge, the performance of PSO in phase stability and equilibrium calculations for both multicomponent reactive and non-reactive mixtures has not yet been reported. This study introduces the application of particle swarm optimization and several of its variants for solving phase stability and equilibrium problems in multicomponent systems with or without chemical equilibrium. The reliability and efficiency of a number of particle swarm optimization algorithms are tested and compared using multicomponent systems with vapor–liquid and liquid–liquid equilibrium. Our results indicate that the classical particle swarm optimization with constant cognitive and social parameters is a reliable method and offers the best performance for global minimization of the tangent plane distance function and the Gibbs energy function in both reactive and non-reactive systems.     

Parameterization of a reactive force field using a Monte Carlo algorithm

Parameterization of a molecular dynamics force field is essential in realistically modeling the physicochemical processes involved in a molecular system. This step is often challenging when the equations involved in describing the force field are complicated as well as when the parameters are mostly empirical. ReaxFF is one such reactive force field which uses hundreds of parameters to describe the interactions between atoms. The optimization of the parameters in ReaxFF is done such that the properties predicted by ReaxFF matches with a set of quantum chemical or experimental data. Usually, the optimization of the parameters is done by an inefficient single‐parameter parabolic‐search algorithm. In this study, we use a robust metropolis Monte‐Carlo algorithm with simulated annealing to search for the optimum parameters for the ReaxFF force field in a high‐dimensional parameter space. The optimization is done against a set of quantum chemical data for MgSO₄ hydrates. The optimized force field reproduced the chemical structures, the equations of state, and the water binding curves of MgSO₄ hydrates. The transferability test of the ReaxFF force field shows the extend of transferability for a particular molecular system. This study points out that the ReaxFF force field is not indefinitely transferable. © 2013 Wiley Periodicals, Inc.     

Peptide structure prediction using distributed volunteer computing networks

Recent investigations to develop novel antimicrobial, antibiotical drugs have focused on the development of artificial protein peptides. As short peptides are naturally involved in many important biological processes in the cell and therefore target many kinds of cells. To functionalize peptides it is vital to design peptides, which can differentially target bacterial and eucariotic cells. Although the length of the peptides investigated in this study was limited to 16 amino acids, the number of possible peptide sequences is still too large to synthesize them in a trial- and error manner, therefore requiring a method for directed, but also high-througput peptide design. By predicting the structure of peptide proteins, this design process can be supported through structure-function analysis and peptide-membrane interaction simulation. In this investigation we could predict peptide structures de-novo, i.e. with the sequence information alone, using a massively parallel simulation scheme. We sample a sizable fraction of the peptide’s conformational space using Monte-Carlo simulations in the free-energy forcefield PFF02 on the volunteer computing network POEM@HOME. This forcefield models the protein’s native conformation as the global minimum of the free-energy. We could identify peptides of different topologies in a completely automated manner, which allows for the high-throughput screening of large peptide databases for their structural features, which would allow the rapid protopying of peptides needed for novel peptide design.     

Optimal pairwise and non-pairwise alchemical pathways for free energy calculations of molecular transformation in solution phase

We estimate the global minimum variance path for computing the free energy insertion into or deletion of small molecules from a dense fluid. We perform this optimization over all pair potentials, irrespective of functional form, using functional optimization with a two-body approximation for the radial distribution function. Surprisingly, the optimal pairwise path obtained via this method is almost identical to the path obtained using a optimized generalized "soft core" potential reported by Pham and Shirts [J. Chem. Phys. 135, 034114 (2011)]. We also derive the lowest variance non-pairwise potential path for molecular insertion or deletion and compare its efficiency to the pairwise path. Under certain conditions, non-pairwise pathways can reduce the total variance by up to 60% compared to optimal pairwise pathways. However, optimal non-pairwise pathways do not appear generally feasible for practical free energy calculations because an accurate estimate of the free energy, the parameter that is itself is desired, is required for constructing this non-pairwise path. Additionally, simulations at most intermediate states of these non-pairwise paths have significantly longer correlation times, often exceeding standard simulation lengths for solvation of bulky molecules. The findings suggest that the previously obtained soft core pathway is the lowest variance pathway for molecular insertion or deletion in practice. The findings also demonstrate the utility of functional optimization for determining the efficiency of thermodynamic processes performed with molecular simulation.     

Investigation of the effect of correlated uncertain rate parameters via the calculation of global and local sensitivity indices

Applications of global uncertainty methods for models with correlated parameters are essential to investigate chemical kinetics models. A global sensitivity analysis method is presented that is able to handle correlated parameter sets. It is based on the coupling of the Rosenblatt transformation with an optimized Random Sampling High Dimensional Model Representation method. The accuracy of the computational method was tested on a series of examples where the analytical solution was available. The capabilities of the method were also investigated by exploring the effect of the uncertainty of rate parameters of a syngas–air combustion mechanism on the calculated ignition delay times. Most of the parameters have large correlated sensitivity indices and the correlation between the parameters has a high influence on the results. It was demonstrated that the values of the calculated total correlated and final marginal sensitivity indices are independent of the order of the decorrelation steps. The final marginal sensitivity indices are meaningful for the investigation of the chemical significance of the reaction steps. The parameters belonging to five elementary reactions only, have significant final marginal sensitivity indices. Local sensitivity indices for correlated parameters were defined which are the linear equivalents of the global ones. The results of the global sensitivity analysis were compared with the corresponding results of local sensitivity analysis for the case of the syngas–air combustion system. The same set of reactions was indicated to be important by both approaches.     

Uncertainty analysis of correlated stability constants

Every attempt of using a computer to model reality has two main uncertainties: the conceptual uncertainty and the data uncertainty. The conceptual uncertainty deals with the choice of model selected for the simulation and the data uncertainty is about the precision and accuracy of the input data. They are often determined experimentally and may thus be encumbered by a number of uncertainties. Normally when treating uncertainties in input data these data are treated as independent variables. However, since many of these parameters are determined together they are actually correlated. This paper focuses on chemical stability constants, a most important parameter for chemical calculations based on speciation. Commonly in the literature they are at best given with an uncertainty interval. We propose to also give the covariance matrix thus giving the opportunity to really assess correlations. In addition we discuss the effect of these correlations on speciations.     

Uncertainty analysis of correlated parameters in automated reaction mechanism generation

Uncertainty analysis is a useful tool for inspecting and improving detailed kinetic mechanisms because it can identify the greatest sources of model output error. Owing to the very nonlinear relationship between kinetic and thermodynamic parameters and computed concentrations, model predictions can be extremely sensitive to uncertainties in some parameters while uncertainties in other parameters can be irrelevant. Error propagation becomes even more convoluted in automatically generated kinetic models, where input uncertainties are correlated through kinetic rate rules and thermodynamic group values. Local and global uncertainty analyses were implemented and used to analyze error propagation in Reaction Mechanism Generator (RMG), an open-source software for generating kinetic models. A framework for automatically assigning parameter uncertainties to estimated thermodynamics and kinetics was created, enabling tracking of correlated uncertainties. Local first-order uncertainty propagation was implemented using sensitivities computed natively within RMG. Global uncertainty analysis was implemented using adaptive Smolyak pseudospectral approximations as implemented in the MIT Uncertainty Quantification Library to efficiently compute and construct polynomial chaos expansions to approximate the dependence of outputs on a subset of uncertain inputs. Cantera was used as a backend for simulating the reactor system in the global analysis. Analyses were performed for a phenyldodecane pyrolysis model. Local and global methods demonstrated similar trends; however, many uncertainties were significantly overestimated by the local analysis. Both local and global analyses show that correlated uncertainties based on kinetic rate rules and thermochemical groups drastically reduce a model's degrees of freedom and have a large impact on the determination of the most influential input parameters. These results highlight the necessity of incorporating uncertainty analysis in the mechanism generation workflow.