On the Uniqueness of the Reverse Monte Carlo Simulation for Molecular Liquids

Abstract

A careful analysis of the three dimensional structures of liquid Chlorine produced by the Reverse Monte Carlo (RMC) and Molecular Dynamics (MD) techniques is presented. The analysis allows us to measure the degree of uniqueness between the potential and the atom-atom distribution functions, g _aa(r), in the case of pairwise potentials formed by isotropic and anisotropic site-site interactions. The g _aa(r) obtained from MD simulations are used as ‘experimental’ input data in the RMC procedure and the constraint of rigid molecules is imposed. The particle configurations produced by RMC are then studied by using a recently proposed general method for analysing the local order in liquids. The same analysis applied to the particle configurations produced by the conventional MD simulation yields a set of partial distribution functions which relates the main features of the g _aa(r) to microscopic pair geometries. The comparison between the partial centre-centre g _cc(r) shows that the three dimensional structures, produced by MD and RMC simulations, agree very well when only isotropic site-site interactions act. In this case RMC produces the same radial distribution function g(r, ω₁, ω₂) as that obtained from the original MD configurations; it is therefore a valid tool for deriving a complete information on the physical properties of a fluid. For anisotropic site-site interactions the partial g _cc(r) of MD and RMC differ significantly and show that the three dimensional structures, produced by MD and RMC simulations, differ too. The discrepancies are particularly evident for the T shaped configurations and affect the values of the potential energy. Therefore, even if the potential is purely pairwise additive, the use of the atomic radial distribution function as input data and the imposition of atomic constraints which model the molecules as hard dumbbells are not sufficient to bring the RMC procedure towards the ‘true’ microscopic structure of the liquid; the presence of non central forces between sites disrupts the bijective correspondence between the potential and the g _aa(r).     

Introducing a novel method based on the imperialistic competitive algorithm to determine fluorine intermolecular potential from ab initio calculations and calculation of some properties via MD simulations

In this work, the possibility of obtaining an accurate site-site potential model suitable for use in molecular dynamics (MD) simulations of fluorine from ab initio calculations has been explored. The exploration was made on ab initio calculations. To reduce the ab initio pair potentials into a site-site potential, a higher significance was assigned to the configuration which is more stable. For this purpose, the imperialistic competitive algorithm (ICA) was implemented as a powerful optimisation tool. The calculated second virial coefficients were compared to the experimental values to test the quality of the presented intermolecular potential. The relative error for the calculated second virial coefficient ranged from 0.1 to 5.6%. MD simulations were used to evaluate the ability of the proposed intermolecular potential function. The relative error for the MD simulations ranged from 0.5 to 5.2%. The results are in good agreement with experimental data.     

Molecular Dynamics on an Excel Spreadsheet

Abstract

We discuss some of the problems that have frustrated the development of reliable model intermolecular potentials for polyatomic molecules. In particular, the usual assumption of an isotropic atom-atom model potential is analysed, and evidence for its inadequacies is presented. A new approach to designing model potentials, an anisotropic site—site model, is introduced by describing several applications to both small and organic molecules, including molecular dynamics and Monte Carlo simulations. The anisotropy required in an atom—atom potential can be directly linked to the non-spherical features in the valence electron distribution, such as lone pairs and π electrons. An accurate electrostatic model for these effects can be constructed from a distributed multipole analysis of the ab initio wavefunction. The empirically required forms of anisotropy in the repulsion potential can also be qualitatively linked to the molecular electron density difference map. Thus, consideration of the molecular bonding can be a useful indication of how to construct adequate model intermolecular pair potentials.     

Coarse-grained simulations for organic molecular liquids based on Gay-Berne and electric multipole potentials

Coarse-grained studies of CH₃SH, CH₃CHO and CHCl₃ liquids, based on anisotropic Gay-Berne (GB) and electric multipole potentials (EMP), demonstrate that the coarse-grained model is able to qualitatively reproduce the results obtained from the atomistic model (AMOEBA polarizable force field) and allows for significant saving in computation time. It should be pointed out that the accuracy of the coarse-grained model is very sensitive to how well the anisotropic GB particle is defined and how satisfactorily the EMP sites are chosen.

Protein structural variation in computational models and crystallographic data

Normal mode analysis offers an efficient way of modeling the conformational flexibility of protein structures. We use anisotropic displacement parameters from crystallography to test the quality of prediction of both the magnitude and directionality of conformational flexibility. Normal modes from four simple elastic network model potentials and from the CHARMM force field are calculated for a data set of 83 diverse, ultrahigh-resolution crystal structures. While all five potentials provide good predictions of the magnitude of flexibility, all-atom potentials have a clear edge at prediction of directionality, and the CHARMM potential has the highest prediction quality. The low-frequency modes from different potentials are similar, but those computed from the CHARMM potential show the greatest difference from the elastic network models. The comprehensive evaluation demonstrates the costs and benefits of using normal mode potentials of varying complexity.     

Reduced C(beta) statistical potentials can outperform all-atom potentials in decoy identification

We developed a series of statistical potentials to recognize the native protein from decoys, particularly when using only a reduced representation in which each side chain is treated as a single C(beta) atom. Beginning with a highly successful all-atom statistical potential, the Discrete Optimized Protein Energy function (DOPE), we considered the implications of including additional information in the all-atom statistical potential and subsequently reducing to the C(beta) representation. One of the potentials includes interaction energies conditional on backbone geometries. A second potential separates sequence local from sequence nonlocal interactions and introduces a novel reference state for the sequence local interactions. The resultant potentials perform better than the original DOPE statistical potential in decoy identification. Moreover, even upon passing to a reduced C(beta) representation, these statistical potentials outscore the original (all-atom) DOPE potential in identifying native states for sets of decoys. Interestingly, the backbone-dependent statistical potential is shown to retain nearly all of the information content of the all-atom representation in the C(beta) representation. In addition, these new statistical potentials are combined with existing potentials to model hydrogen bonding, torsion energies, and solvation energies to produce even better performing potentials. The ability of the C(beta) statistical potentials to accurately represent protein interactions bodes well for computational efficiency in protein folding calculations using reduced backbone representations, while the extensions to DOPE illustrate general principles for improving knowledge-based potentials.     

Distance-scaled,finite ideal-gas reference state improves structure-derived potentials of mean force for structure selection and stability prediction

The distance-dependent structure-derived potentials developed so far all employed a reference state that can be characterized as a residue (atom)-averaged state. Here, we establish a new reference state called the distance-scaled, finite ideal-gas reference (DFIRE) state. The reference state is used to construct a residue-specific all-atom potential of mean force from a database of 1011 nonhomologous (less than 30% homology) protein structures with resolution less than 2 A. The new all-atom potential recognizes more native proteins from 32 multiple decoy sets, and raises an average Z-score by 1.4 units more than two previously developed, residue-specific, all-atom knowledge-based potentials. When only backbone and C(beta) atoms are used in scoring, the performance of the DFIRE-based potential, although is worse than that of the all-atom version, is comparable to those of the previously developed potentials on the all-atom level. In addition, the DFIRE-based all-atom potential provides the most accurate prediction of the stabilities of 895 mutants among three knowledge-based all-atom potentials. Comparison with several physical-based potentials is made.     

MD Simulations of Liquids with Td,Oh Molecular Symmetry

Abstract

A new method is proposed for the calculation of intermolecular interactions in Molecular Dynamics simulations of liquids with T_d, O_h molecular symmetry. The new algorithm is based on the separation of the pair potential into a short-range and a long-range contribution described by a site-site and a spherical centre-centre potential model respectively using an additional cutoff distance. Test calculations for the Lennard-Jones fluids CCl₄ and SF₆ show significant savings in CPU time. We compare thermodynamic properties, pair correlation functions and a few dynamic autocorrelation functions obtained with the novel strategy with results of the commonly used algorithm for systems containing 864 molecules. Since no significant differences appear the new algorithm may be suggested as a useful contribution to the area of Molecular Dynamics simulation of liquids with these rather high molecular symmetries.     

Multiscale Coarse-Graining of the Protein Energy Landscape

A variety of coarse-grained (CG) models exists for simulation of proteins. An outstanding problem is the construction of a CG model with physically accurate conformational energetics rivaling all-atom force fields. In the present work, atomistic simulations of peptide folding and aggregation equilibria are force-matched using multiscale coarse-graining to develop and test a CG interaction potential of general utility for the simulation of proteins of arbitrary sequence. The reduced representation relies on multiple interaction sites to maintain the anisotropic packing and polarity of individual sidechains. CG energy landscapes computed from replica exchange simulations of the folding of Trpzip, Trp-cage and adenylate kinase resemble those of other reduced representations; non-native structures are observed with energies similar to those of the native state. The artifactual stabilization of misfolded states implies that non-native interactions play a deciding role in deviations from ideal funnel-like cooperative folding. The role of surface tension, backbone hydrogen bonding and the smooth pairwise CG landscape is discussed. Ab initio folding aside, the improved treatment of sidechain rotamers results in stability of the native state in constant temperature simulations of Trpzip, Trp-cage, and the open to closed conformational transition of adenylate kinase, illustrating the potential value of the CG force field for simulating protein complexes and transitions between well-defined structural states.     

Conformational gating in ammonia lyases

BackgroundAmmonia lyases are enzymes of industrial and biomedical interest. Knowledge of structure-dynamics-function relationship in ammonia lyases is instrumental for exploiting the potential of these enzymes in industrial or biomedical applications.MethodsWe investigated the conformational changes in the proximity of the catalytic pocket of a 3-methylaspartate ammonia lyase (MAL) as a model system. At this scope, we used microsecond all-atom molecular dynamics simulations, analyzed with dimensionality reduction techniques, as well as in terms of contact networks and correlated motions.ResultsWe identify two regulatory elements in the MAL structure, i.e., the β5-α2 loop and the helix-hairpin-loop subdomain. These regulatory elements undergo conformational changes switching from ‘occluded’ to ‘open’ states. The rearrangements are coupled to changes in the accessibility of the active site. The β5-α2 loop and the helix-hairpin-loop subdomain modulate the formation of tunnels from the protein surface to the catalytic site, making the active site more accessible to the substrate when they are in an open state.ConclusionsOur work pinpoints a sequential mechanism, in which the helix-hairpin-loop subdomain of MAL needs to break a subset of intramolecular interactions first to favor the displacement of the β5-α2 loop. The coupled conformational changes of these two elements contribute to modulate the accessibility of the catalytic site.General significanceSimilar molecular mechanisms can have broad relevance in other ammonia lyases with similar regulatory loops. Our results also imply that it is important to account for protein dynamics in the design of variants of ammonia lyases for industrial and biomedical applications.     

Computer Simulation of the Dielectric Properties of Liquid Water

Abstract

Monte Carlo simulations of water in the NVT ensemble using three models (SPC, TIP4P and TIPS2) are reported. The internal energy, dielectric constant, and the site-site radial distribution functions of liquid water (temperature 300 K and mass density 1 gm cc^?1) were calculated and compared with experiment. It was found that of the three intermolecular potential models, SPC gives the best dielectric constant. Since SPC also yields acceptable results for the energy and structure, it is judged to be the best among the three models studied.     

Distinguishing native conformations of proteins from decoys with an effective free energy estimator based on the OPLS all-atom force field and the Surface Generalized Born solvent model

Protein decoy data sets provide a benchmark for testing scoring functions designed for fold recognition and protein homology modeling problems. It is commonly believed that statistical potentials based on reduced atomic models are better able to discriminate native-like from misfolded decoys than scoring functions based on more detailed molecular mechanics models. Recent benchmark tests on small data sets, however, suggest otherwise. In this work, we report the results of extensive decoy detection tests using an effective free energy function based on the OPLS all-atom (OPLS-AA) force field and the Surface Generalized Born (SGB) model for the solvent electrostatic effects. The OPLS-AA/SGB effective free energy is used as a scoring function to detect native protein folds among a total of 48,832 decoys for 32 different proteins from Park and Levitt's 4-state-reduced, Levitt's local-minima, Baker's ROSETTA all-atom, and Skolnick's decoy sets. Solvent electrostatic effects are included through the Surface Generalized Born (SGB) model. All structures are locally minimized without restraints. From an analysis of the individual energy components of the OPLS-AA/SGB energy function for the native and the best-ranked decoy, it is determined that a balance of the terms of the potential is responsible for the minimized energies that most successfully distinguish the native from the misfolded conformations. Different combinations of individual energy terms provide less discrimination than the total energy. The results are consistent with observations that all-atom molecular potentials coupled with intermediate level solvent dielectric models are competitive with knowledge-based potentials for decoy detection and protein modeling problems such as fold recognition and homology modeling.     

Predicting Liquid–Vapour Equilibria for Water Using an ab-initio Potential from Histogram Reweighting Monte Carlo Simulations

Abstract

The coexisting densities for an ab-initio model for water have been calculated using grand canonical Monte Carlo simulations with the histogram reweighting technique. Although good agreement with experimental data is found for the radial distribution function at room temperature, the predicted critical density and temperature are well below both the experimental value as well as predictions from semi-empirical potentials. Improvement in the repulsive part of the ab-initio potential is suggested as a way to obtain better agreement with experiment.     

Ion transport and displacement currents with membrane-bound carriers

Summary The standard carrier model for ion transport by a one-to-one mechanism is developed to predict the time-dependent currents for systems that are symmetrical at zero applied potential. The complete solution for ions and carriers bearing any charge is derived by assuming that the concentration of ions in the membrane is low and either that the applied potential is small or that the applied potential affects equally all of the association and dissociation reactions between the ions and the carriers. The response to an abruptly applied potential is then given by the sum of a constant and two declining exponential terms. The time constants of these relaxations are described by the equations derived for neutral carriers by Stark, Ketterer, Benz and Läuger in 1971 (Biophys. J. 11:981). The sum of the amplitudes of the exponentials for small applied potentials obeys a relation like that first derived by Markin and Liberman in 1973 (Biofizika 18:453). For small applied potentials expressions are also provided for the voltage transients in charge-pulse experiments and for the membrane admittance.     

Correlation between knowledge-based and detailed atomic potentials: application to the unfolding of the GCN4 leucine zipper.

The relationship between the unfolding pseudo free energies of reduced and detailed atomic models of the GCN4 leucine zipper is examined. Starting from the native crystal structure, a large number of conformations ranging from folded to unfolded were generated by all-atom molecular dynamics unfolding simulations in an aqueous environment at elevated temperatures. For the detailed atomic model, the pseudo free energies are obtained by combining the CHARMM all-atom potential with a solvation component from the generalized Born, surface accessibility, GB/SA, model. Reduced model energies were evaluated using a knowledge-based potential. Both energies are highly correlated. In addition, both show a good correlation with the root mean square deviation, RMSD, of the backbone from native. These results suggest that knowledge-based potentials are capable of describing at least some of the properties of the folded as well as the unfolded states of proteins, even though they are derived from a database of native protein structures. Since only conformations generated from an unfolding simulation are used, we cannot assess whether these potentials can discriminate the native conformation from the manifold of alternative, low-energy misfolded states. Nevertheless, these results also have significant implications for the development of a methodology for multiscale modeling of proteins that combines reduced and detailed atomic models.     

Chemical Potential,Partial Enthalpy and Partial Volume of Mixtures by NPT Molecular Dynamics

Abstract

Equilibrium NPT molecular dynamics computer simulations have been used to determine the chemical potential, partial enthalpy and partial volume of model Ar-Kr mixtures using newly devised non-intrusive particle insertion and particle swap techniques [P. Sindzingre et al. Chemical Physics, 129 (1989) 213]. In this report we examine, for the first time, in some detail the relative convergence statistics of the particle swap and particle insertion methods for these properties for binary Lennard-Jones (LJ) mixtures. Both species are represented by single-site Lennard-Jones pair potentials with Lorentz-Berthelot rules for the cross-species interactions. We show that, over the whole phase diagram and especially in the vicinity of the fluid-solid coexistence line, the particle swap method gives significantly better statistics than the particle insertion method for the difference in chemical potential of the two species, partial enthalpy and partial volume of each species. Also, we find that, using the particle swap method, the difference in the chemical potential converges more rapidly than the differences in the partial enthalpy and volume.     

Extracting the excitonic Hamiltonian of the Fenna-Matthews-Olson complex using three-dimensional third-order electronic spectroscopy

An accurate scoring function is a key component for successful protein structure prediction. To address this important unsolved problem, we develop a generalized orientation and distance-dependent all-atom statistical potential. The new statistical potential, generalized orientation-dependent all-atom potential (GOAP), depends on the relative orientation of the planes associated with each heavy atom in interacting pairs. GOAP is a generalization of previous orientation-dependent potentials that consider only representative atoms or blocks of side-chain or polar atoms. GOAP is decomposed into distance- and angle-dependent contributions. The DFIRE distance-scaled finite ideal gas reference state is employed for the distance-dependent component of GOAP. GOAP was tested on 11 commonly used decoy sets containing 278 targets, and recognized 226 native structures as best from the decoys, whereas DFIRE recognized 127 targets. The major improvement comes from decoy sets that have homology-modeled structures that are close to native (all within ∼4.0 Å) or from the ROSETTA ab initio decoy set. For these two kinds of decoys, orientation-independent DFIRE or only side-chain orientation-dependent RWplus performed poorly. Although the OPUS-PSP block-based orientation-dependent, side-chain atom contact potential performs much better (recognizing 196 targets) than DFIRE, RWplus, and dDFIRE, it is still ∼15% worse than GOAP. Thus, GOAP is a promising advance in knowledge-based, all-atom statistical potentials. GOAP is available for download at http://cssb.biology.gatech.edu/GOAP.     

Estimation of Solvent Diffusion Coefficients Using Molecular Dynamics Simulations

Abstract

Molecular dynamics simulations have been used to investigate diffusion in two commonly used industrial solvents, toluene and tetrahydrofuran. Several different models for the solvents are compared (flexible vs. rigid, all-atom vs. united atom), and it is found that united atom and all-atom models of the solvents produce very different diffusion coefficients at the experimental density. This disagreement can be explained by the pressure dependence of the diffusion coefficient, which is found to vary in accord with the Chapman-Enskog result for hard spheres. It is recommended that force fields be parametrized carefully to produce reasonable pressures at the experimental densities, or that simulations be carried out at constant pressure, if they are to be used for the purposes of calculating transport coefficients.     

The inconsistencies of present-day mathematical-physical methods pertaining to current generators in volume conductors used in electrocardiography

The present-day practices of electrocardiography and vectorardiography are based upon the theory that the surface potential differences can be assumed to be due to a single dipole inside the body. It is shown in this paper that a dipole cannot account for all the surface potentials due to realistic current generators, and hence the determination of the current generator from surface potential measurements based upon such a theory will lead to inconsistent representations of the heart for one and the same subject. To demonstrate this point two eccentric dipoles of different strengths and locations representing two muscle fibers are taken to be the current generator in a homogeneous spherical conductor. The exact surface potentials are then expressed by means of the “interior sphere theorem” of the authors. With these expressions the magnitude, direction, and location of the resultant dipole are determined by the method of D. Gabor and C. V. Nelson (J. App. Physics,25, 413–16, 1954). The surface potentials due to this resultant dipole are again exactly expressed by means of the “interior sphere theorem” and compared with those due to the eccentric dipoles assumed. It can be seen that the differences can be considerable. It is suggested that the multipole model of the authors (Bull. Math. Biophysics,20, 203–16, 1958) be used as a more accurate and the only unique representation of the heart. This investigation was supported by the National Heart Institute under a research grant H-2263(c).