Simple Energy Function Research Articles

New correlation energy functionals with explicit dependence on the number of electrons

Based upon the idea of effective number of electrons, we develop simple but accurate correlation energy functionals to be used for density functional theory calculations. We derive both a spin-independent and a spin-dependent functional. The spin-dependent one, used in conjunction with Becke’s exchange functional [A. D. Becke, Phys. Rev. A 38, 3098 (1988)], yields excellent results for thermochemistry calculations, giving an average absolute error of 2.9 kcal/mol for a test set comprised of the enthalpies of formation of the 148 molecules in the extended G2 set [L. A. Curtiss, K. Raghavachari, P. C. Redfern, and J. A. Pople, J. Chem. Phys. 106, 1063 (1997); L. A. Curtiss, P. C. Redfern, K. Raghavachari, and J. A. Pople, ibid. 109, 42 (1998)] plus the total energies of the atoms H through Ar. We also discuss the problem of fractional occupation number, and we show that the corresponding principle of integer preference can be fulfilled by the procedure that we propose to build correlation energy functionals.

The Mechanism and Pathway of pH Induced Swelling in Cowpea Chlorotic Mottle Virus

Normal mode analysis based on a simplified energy function was used to study the swelling process of the icosahedral virus, cowpea chlorotic mottle virus (CCMV). Native state virus particles (coat proteins) of this T=3 icosahedral virus have been shown to undergo a large conformational change to a swollen state when metal ions are removed or the pH is raised. A normal mode analysis based on the native state capsid showed one preferential direction, a breathing mode, that explains the majority of the structural rearrangement necessary to bring the native structure close to the swollen state. From the native form of CCMV, the structure can be displaced along the direction of a single breathing mode by different amounts to create several candidate swollen structures and a putative pathway for virus expansion. The R-factor between these predicted swollen capsid structures and experimental electron density from cryoelectron microscopy (cryo-EM) measurements is then calculated to indicate how well each structure satisfies the experimental measurements on the swollen capsid state. A decrease of the crystallographic R-factor value from ∼72% to ∼49% was observed for these simple incremental displacements along the breathing mode. The simultaneous displacement of the native structure along other relevant (symmetric, non-degenerate) modes produce a structure with an R-factor of 45%, which is further reduced to 43.9% after minimization: a value in good accord with models based on the EM data at 28 Å resolution. Based on the incrementally expanded structures, a pathway for the swelling process has been proposed. Analysis of the intermediate structures along this pathway indicates a significant loss of interactions at the quasi-3-fold interfaces occurs in the initial stages of the swelling process and this serves as a trigger for the compact to swollen transition. Furthermore, the pH dependent swelling appears to be triggered by the titration of a single residue with an anomalous p K a value in the unswollen particle.

Realistic modelling of earthquake-isolation bearings

The physical properties of rubber components are greatly influenced by the techniques used in their manufacture. Normally, high-order strain energy (density) functions are required in finite-element analyses (FEAs) to represent elastomeric behaviour. Natural rubber is bonded to steel in alternate layers to produce earthquake-isolation bearings and these bearings are predominantly subjected to compression and shear. In compression, rubber is thought to behave in compliance with a simple neo-Hookean strain energy function, whilst tests in shear provide only information about low-order functions. This suggests that it may be possible to model the bearings with a comparatively simple equation. Natural rubber test-pieces were manufactured by three different processes over a range of hardness. Initially, a Flory–Rhener test was applied to each compound to establish the variation in cross-link density. Uniaxial, compression and tension tests were carried out on specimens produced by compression, transfer and injection moulding. Neo-Hookean and two-term Ogden models were determined for each compound for input to FEA. Simultaneous shear and compression loads were applied to a rectangular block comprising rubber and steel layers. Results from this test were compared with a plane-stress finite-element model. In conclusion, recommendations are made for the manufacture, computer modelling and testing of the bearings.

Folding funnels: the key to robust protein structure prediction.

Natural proteins fold because their free energy landscapes are funneled to their native states. The degree to which a model energy function for protein structure prediction can avoid the multiple minima problem and reliably yield at least low-resolution predictions is also dependent on the topography of the energy landscape. We show that the degree of funneling can be quantitatively expressed in terms of a few averaged properties of the landscape. This allows us to optimize simplified energy functions for protein structure prediction even in the absence of homology information. Here we outline the optimization procedure in the context of associative memory energy functions originally introduced for tertiary structure recognition and demonstrate that even partially funneled landscapes lead to qualitatively correct, low-resolution predictions.

Hierarchy of simulation models in predicting molecular recognition mechanisms from the binding energy landscapes: Structural analysis of the peptide complexes with SH2 domains

Computer simulations using the simplified energy function and simulated tempering dynamics have accurately determined the native structure of the pYVPML, SVLpYTAVQPNE, and SPGEpYVNIEF peptides in the complexes with SH2 domains. Structural and equilibrium aspects of the peptide binding with SH2 domains have been studied by generating temperature-dependent binding free energy landscapes. Once some native peptide-SH2 domain contacts are constrained, the underlying binding free energy profile has the funnel-like shape that leads to a rapid and consistent acquisition of the native structure. The dominant native topology of the peptide-SH2 domain complexes represents an extended peptide conformation with strong specific interactions in the phosphotyrosine pocket and hydrophobic interactions of the peptide residues C-terminal to the pTyr group. The topological features of the peptide-protein interface are primarily determined by the thermodynamically stable phosphotyrosyl group. A diversity of structurally different binding orientations has been observed for the amino-terminal residues to the phosphotyrosine. The dominant native topology for the peptide residues carboxy-terminal to the phosphotyrosine is tolerant to flexibility in this region of the peptide-SH2 domain interface observed in equilibrium simulations. The energy landscape analysis has revealed a broad, entropically favorable topology of the native binding mode for the bound peptides, which is robust to structural perturbations. This could provide an additional positive mechanism underlying tolerance of the SH2 domains to hydrophobic conservative substitutions in the peptide specificity region.

Molecular docking of α-cyclodextrin inclusion complexes by genetic algorithm and empirical binding free energy function

A molecular docking method that predicts the lowest energy geometries of inclusion complexes between host and guests was developed and tested, in combination with a new simple empirical function that estimates the free energy of binding. The total interaction energies of the host–guest inclusion complexes were optimized using a genetic algorithm (GA). The docking method was applied to 43 complexes of α-cyclodextrin (α-CD) and mono- or 1,4-disubstituted benzenes with known binding constants. The new simple empirical free energy function was calibrated by the 43 docked complex structures and gave a good relationship between the predicted binding constants and the observed values.

Conformational propagation with prion-like characteristics in a simple model of protein folding.

Protein refolding/misfolding to an alternative form plays an aetiologic role in many diseases in humans, including Alzheimer's disease, the systemic amyloidoses, and the prion diseases. Here we have discovered that such refolding can occur readily for a simple lattice model of proteins in a propagatable manner without designing for any particular alternative native state. The model uses a simple contact energy function for interactions between residues and does not consider the peculiarities of polypeptide geometry. In this model, under conditions where the normal (N) native state is marginally stable or unstable, two chains refold from the N native state to an alternative multimeric energetic minimum comprising a single refolded conformation that can then propagate itself to other protein chains. The only requirement for efficient propagation is that a two-faced mode of packing must be in the ground state as a dimer (a higher-energy state for this packing leads to less efficient propagation). For random sequences, these ground-state dimeric configurations tend to have more beta-sheet-like extended structure than almost any other sort of dimeric ground-state assembly. This implies that propagating states (such as for prions) are beta-sheet rich because the only likely propagating forms are beta-sheet rich. We examine the details of our simulations to see to what extent the observed properties of prion propagation can be predicted by a simple protein folding model. The formation of the alternative state in the present model shows several distinct features of amyloidogenesis and of prion propagation. For example, an analog of the phenomenon of conformationally distinct strains in prions is observed. We find a parallel between 'glassy' behavior in liquids and the formation of a propagatable state in proteins. This is the first report of simulation of conformational propagation using any heteropolymer model. The results imply that some (but not most) small protein sequences must maintain a sequence signal that resists refolding to propagatable alternative native states and that the ability to form such states is not limited to polypeptides (or reliant on regular hydrogen bonding per se) but can occur for other protein-like heteropolymers.

Simple Landau model of the smectic- A-isotropic phase transition

A simple Landau-type free energy function is presented to describe the smectic-A-isotropic phase transition. Varying the coupling between orientational and positional order parameters, a smectic-A-isotropic or a nematic-isotropic phase transition occurs. Within this model the smectic-A-isotropic phase transition is found to be always more strongly first order than the nematic-isotropic phase transition. The theoretical results are found to be in good agreement with all published experimental results.

A Statistical Mechanical Model for Predicting B-DNA Curvature and Flexibility

A statistical mechanical model taking into account the symmetric twisting, tilting, sliding fluctuations and asymmetric rolling fluctuations has been proposed to predict the macroscopic curvature and flexibility of B-DNA. Based on the statistical data of structural parameters of double helix in nucleic acid database and the related theoretical analysis, the equilibrium angular parameters (Ω, ρ and τ) describing the orientation of successive base-pair planes, the translation parameters (Dy) along the long axis of neighboring base-pair step and the corresponding force constants are arranged for ten dimers appropriately. Under the assumption of independent angular parameters, independent base-pair steps and a simple energy function, we can calculate the macroscopic curvature and the flexibility of DNA sequences through the transformation matrix 〈R〉 and the Boltzmann ensemble average. The predictions on curvature and flexibility of DNA have been compared with the corresponding experimental data. The agreement is remarkably good. It is demonstrated that the lowering of the temperature does increase the DNA curvature.

Optimal annealing schedules for two-, three-, and four-level systems using a genetic algorithm approach

An annealing schedule, T(t), is the temperature as function of time whose goal is to bring a system from some initial low-order state to a final high-order state. We use the probability in the lowest energy level as the order parameter, so that an ideally annealed system would have all its population in its ground-state. We consider a model system comprised of discrete energy levels separated by activation barriers. We have carried out annealing calculations on this system for a range of system parameters. In particular, we considered the schedule as a function of the energy level spacing, of the height of the activation barriers, and, in some cases, as a function of degeneracies of the levels. For a given set of physical parameters, and maximum available time, tm, we were able to obtain the optimal schedule by using a genetic algorithm (GA) approach. For the two-level system, analytic solutions are available, and were compared with the GA-optimized results. The agreement was essentially exact. We were able to identify systematic behaviors of the schedules and trends in final probabilities as a function of parameters. We have also carried out Metropolis Monte Carlo (MMC) calculations on simple potential energy functions using the optimal schedules available from the model calculations. Agreement between the model and MMC calculations was excellent.

Ab initio prediction of the solution structures and populations of a cyclic pentapeptide in DMSO based on an implicit solvation model.

Using a recently developed statistical mechanics methodology, the solution structures and populations of the cyclic pentapeptide cyclo(D-Pro(1)-Ala(2)-Ala(3)-Ala(4)-Ala(5)) in DMSO are obtained ab initio, i.e., without using experimental restraints. An important ingredient of this methodology is a novel optimization of implicit solvation parameters, which in our previous publication [Baysal, C.; Meirovitch, H. J Am Chem Soc 1998, 120, 800-812] has been applied to a cyclic hexapeptide in DMSO. The molecule has been described by the simplified energy function E(tot) = E(GRO) + summation operator(k) sigma(k)A(k), where E(GRO) is the GROMOS force-field energy, sigma(k) and A(k) are the atomic solvation parameter (ASP) and the solvent accessible surface area of atom k. This methodology, which relies on an extensive conformational search, Monte Carlo simulations, and free energy calculations, is applied here with E(tot) based on the ASPs derived in our previous work, and for comparison also with E(GRO) alone. For both models, entropy effects are found to be significant. For E(tot), the theoretical values of proton-proton distances and (3)J coupling constants agree very well with the NMR results [Mierke, D. F.; Kurz, M.; Kessler, H. J Am Chem Soc 1994, 116, 1042-1049], while the results for E(GRO) are significantly worse. This suggests that our ASPs might be transferrable to other cyclic peptides in DMSO as well, making our methodology a reliable tool for an ab initio structure prediction; obviously, if necessary, parts of this methodology can also be incorporated in a best-fit analysis where experimental restraints are used.

Calculation of Protein Surface Loops Using Monte-Carlo Simulated Annealing Simulation

A modified program for protein loop modelling is presented with significant improvements on our previous study (Zhang et al. 1997, Biopolymers, Vol. 41, pp. 61-72), which is capable of sampling the entire conformational space and identifying the low-energy candidates by Monte-Carlo simulated annealing simulation and a cluster analysis method. Twenty flexible surface loops connecting different secondary structures were selected to test the efficiency of this program. The averaged deviations of backbone heavy atoms for four to eight-residue-loops are 0.19, 0.27, 0.46, 0.41 and 0.87A respectively. High speed of calculation is achieved with a simplified energy function and a grid-mapping method. As a comparison of single simulation, it takes only four seconds for the simplest four-residue loop and forty-three seconds for the most complex eight-residue loop on a PII-350-Linux platform.

On the transferability of atomic solvation parameters: Ab initio structural prediction of cyclic heptapeptides in DMSO.

A statistical mechanics methodology for predicting the solution structures and populations of peptides developed recently is based on a novel method for optimizing implicit solvation models, which was applied initially to a cyclic hexapeptide in DMSO (C. Baysal and H. Meirovitch, Journal of American Chemical Society, 1998, vol. 120, pp. 800-812). Thus, the molecule has been described by the simplified energy function E(tot) = E(GRO) + summation operator(k) sigma(k)A(k), where E(GRO) is the GROMOS force-field energy, sigma(k) and A(k) are the atomic solvation parameter (ASP) and the solvent accessible surface area of atom k, respectively. In a more recent study, these ASPs have been found to be transferable to the cyclic pentapeptide cyclo(D-Pro(1)-Ala(2)-Ala(3)-Ala(4)-Ala(5)) in DMSO (C. Baysal and H. Meirovitch, Biopolymers, 2000, vol. 53, pp. 423-433). In the present paper, our methodology is applied to the cyclic heptapeptides axinastatin 2 [cyclo(Asn(1)-Pro(2)-Phe(3)-Val(4)-Leu(5)-Pro(6)-Val(7))] and axinastatin 3 [cyclo(Asn(1)-Pro(2)-Phe(3)-Ile(4)-Leu(5)-Pro(6)-Val(7))], in DMSO, which were studied by nmr by Mechnich et al. (Helvetica Chimica Acta, 1997, vol. 80, pp. 1338-1354). The calculations for axinastatin 2 show that special ASPs should be optimized for the partially charged side-chain atoms of Asn while the rest of the atoms take their values derived in our previous work; this suggests that similar optimization might be needed for other side chains as well. The solution structures of these peptides are obtained ab initio (i.e., without using experimental restraints) by an extensive conformational search based on E(GRO) alone and E(*)(tot), which consists of the new set of ASPs. For E(*)(tot), the theoretical values of proton-proton distances, (3)J coupling constants, and other properties are found to agree very well with the nmr results, and they are always better than those based on E(GRO).

Deciphering common failures in molecular docking of ligand-protein complexes.

Common failures in predicting crystal structures of ligand-protein complexes are investigated for three ligand-protein systems by a combined thermodynamic and kinetic analysis of the binding energy landscapes. Misdocked predictions in ligand-protein docking are classified as 'soft' and 'hard' failures. While a soft failure arises when the search algorithm is unable to find the global energy minimum corresponding to the crystal structure, a hard failure results from a flaw of the energy function to qualify the crystal structure as the predicted lowest energy conformation in docking simulations. We find that neither the determination of a single structure with the lowest energy nor finding the most common binding mode is sufficient to predict crystal structures of the complexes, which belong to the category of hard failures. In a proposed hierarchical approach, structural similarity clustering of the conformations, generated from equilibrium simulations with the simplified energy function, is followed by energy refinement with the AMBER force field. This protocol, that involves a hierarchy of energy functions, resolves some common failures in ligand-protein docking and detects crystallographic binding modes that were not found during docking simulations.

Folding protein models with a simple hydrophobic energy function: the fundamental importance of monomer inside/outside segregation.

The present study explores a "hydrophobic" energy function for folding simulations of the protein lattice model. The contribution of each monomer to conformational energy is the product of its "hydrophobicity" and the number of contacts it makes, i.e., E(h,c) = -Sigma N/i=1 c(i)h(i) = -(h.c) is the negative scalar product between two vectors in N-dimensional cartesian space: h = (h1,., hN), which represents monomer hydrophobicities and is sequence-dependent; and c = (c(1),., c(N)), which represents the number of contacts made by each monomer and is conformation-dependent. A simple theoretical analysis shows that restrictions are imposed concomitantly on both sequences and native structures if the stability criterion for protein-like behavior is to be satisfied. Given a conformation with vector c, the best sequence is a vector h on the direction upon which the projection of c - c is maximal, where c is the diagonal vector with components equal to c, the average number of contacts per monomer in the unfolded state. Best native conformations are suggested to be not maximally compact, as assumed in many studies, but the ones with largest variance of contacts among its monomers, i.e., with monomers tending to occupy completely buried or completely exposed positions. This inside/outside segregation is reflected on an apolar/polar distribution on the corresponding sequence. Monte Carlo simulations in two dimensions corroborate this general scheme. Sequences targeted to conformations with large contact variances folded cooperatively with thermodynamics of a two-state transition. Sequences targeted to maximally compact conformations, which have lower contact variance, were either found to have degenerate ground state or to fold with much lower cooperativity.

Averaging interaction energies over homologs improves protein fold recognition in gapless threading

Protein structure prediction is limited by the inaccuracy of the simplified energy functions necessary for efficient sorting over many conformations. It was recently suggested (Finkelstein, Phys Rev Lett 1998;80:4823-4825) that these errors can be reduced by energy averaging over a set of homologous sequences. This conclusion is confirmed in this study by testing protein structure recognition in gapless threading. The accuracy of recognition was estimated by the Z-score values obtained in gapless threading tests. For threading, we used 20 target proteins, each having from 20 to 70 homologs taken from the HSSP sequence base. The energy of the native structures was compared with the energy from 34 to 75 thousand of alternative structures generated by threading. The energy calculations were done with our recently developed Calpha atom-based phenomenological potentials. We show that averaging of protein energies over homologs reduces the Z-score from approximately -6.1 (average Z-score for individual chains) to approximately -8.1. This means that a correct fold can be found among 3 x 10(9) random folds in the first case and among 3 x 10(15) in the second. Such increase in selectivity is important for recognition of protein folds.

Computer simulations of ligand-protein binding with ensembles of protein conformations: A Monte Carlo study of HIV-1 protease binding energy landscapes

We present the results of molecular docking simulations with HIV-1 protease for the sb203386 and skf107457 inhibitors by Monte Carlo simulated annealing. A simplified piecewise linear energy function, the standard AMBER force field, and the AMBER force field with solvation and a soft-core smoothing component are employed in simulations with a single-protein conformation to determine the relationship between docking simulations with a simple energy function and more realistic force fields. The temperature-dependent binding free energy profiles of the inhibitors interacting with a single protein conformation provide a detailed picture of relative thermodynamic stability and a distribution of ligand binding modes in agreement with experimental crystallographic data. Using the simplified piecewise linear energy function, we also performed Monte Carlo docking simulations with an ensemble of protein conformations employing preferential biased sampling of low-energy protein conformations, and the results are analyzed in connection with the free energy profiles. ©1999 John Wiley & Sons, Inc. Int J Quant Chem 72: 73–84, 1999

Predicting the structures of 18 peptides using Geocore.

We describe an extensive test of Geocore, an ab initio peptide folding algorithm. We studied 18 short molecules for which there are structures in the Protein Data Bank; chains are up to 31 monomers long. Except for the very shortest peptides, an extremely simple energy function is sufficient to discriminate the true native state from more than 10(8) lowest energy conformations that are searched explicitly for each peptide. A high incidence of native-like structures is found within the best few hundred conformations generated by Geocore for each amino acid sequence. Predictions improve when the number of discrete phi/psi choices is increased.

Knowledge-based structure prediction of MHC class I bound peptides: a study of 23 complexes

The binding of T-cell antigenic peptides to MHC molecules is a prerequisite for their immunogenicity. The ability to identify binding peptides based on the protein sequence is of great importance to the rational design of peptide vaccines. As the requirements for peptide binding cannot be fully explained by the peptide sequence per se, structural considerations should be taken into account and are expected to improve predictive algorithms. The first step in such an algorithm requires accurate and fast modeling of the peptide structure in the MHC-binding groove. We have used 23 solved peptide-MHC class I complexes as a source of structural information in the development of a modeling algorithm. The peptide backbones and MHC structures were used as the templates for prediction. Sidechain conformations were built based on a rotamer library, using the 'dead end elimination' approach. A simple energy function selects the favorable combination of rotamers for a given sequence. It further selects the correct backbone structure from a limited library. The influence of different parameters on the prediction quality was assessed. With a specific rotamer library that incorporates information from the peptide sidechains in the solved complexes, the algorithm correctly identifies 85% (92%) of all (buried) sidechains and selects the correct backbones. Under cross-validation, 70% (78%) of all (buried) residues are correctly predicted and most of all backbones. The interaction between peptide sidechains has a negligible effect on the prediction quality. The structure of the peptide sidechains follows from the interactions with the MHC and the peptide backbone, as the prediction is hardly influenced by sidechain interactions. The proposed methodology was able to select the correct backbone from a limited set. The impairment in performance under cross-validation suggests that, currently, the specific rotamer library is not satisfactorily representative. The predictions might improve with an increase in the data.

Structure and property prediction of sub- and supercritical water

In this report two approaches are presented to predict the structure and PVT behavior of associating fluids with emphasis on water. One approach is the development of equations of state based on the analytic chain association theory (ACAT). An associating fluid is assumed to be a mixture of monomers, dimers, trimers, etc., for which the composition distribution is obtained. The resulting equations are simple enough to be used for PVT calculations. The second approach is the development of an effective Kihara pair potential for water which incorporates the hydrogen bonding using the ACAT. This potential function has been used in an analytical expression to predict the first shell of the radial distribution function (RDF) for water. The expression for RDF which was initially applied to simple potential energy functions, such as the Lennard–Jones and Kihara functions satisfies the general functionality of the RDF with respect to intermolecular potential, temperature and density as well as all the limiting values of RDF at high temperature and dilute gas, and infinite separation. The effective potential parameters are determined to predict the first shell of the RDF data for water at various subcritical temperatures and densities. The predicted results for water at near-critical and supercritical conditions are shown to be in agreement with the data obtained by neutron diffraction experiments and with the simulation data.