Adam Liwo

A united-residue force field for off-lattice protein-structure simulations. I. Functional forms and parameters of long-range side-chain interaction potentials from protein crystal data

A two‐stage procedure for the determination of a united‐residue potential designed for protein simulations is outlined. In the first stage, the long‐range and local‐interaction energy terms of the total energy of a polypeptide chain are determined by analyzing protein‐crystal data and averaging the all‐atom energy surfaces. In the second stage (described in the accompanying article), the relative weights of the energy terms are optimized so as to locate the native structures of selected test proteins as the lowest energy structures. The goal of the work in the present study is to parameterize physically reasonable functional forms of the potentials of mean force for side‐chain interactions. The potentials are of both radial and anisotropic type. Radial potentials include the Lennard‐Jones and the shifted Lennard‐Jones potential (with the shift parameter independent of orientation). To treat the angular dependence of side‐chain interactions, three functional forms of the potential that were designed previously to describe anisotropic systems are evaluated: Berne‐Pechukas (dilated Lennard‐Jones); Gay‐Berne (shifted Lennard‐Jones with orientation‐dependent shift parameters); and Gay‐Berne‐Vorobjev (the same as the preceding one, but with one more set of variable parameters). These functional forms were used to parameterize, within a short‐distance range, the potentials of mean force for side‐chain pair interactions that are related by the Boltzmann principle to the pair correlation functions determined from protein‐crystal data. Parameter determination was formulated as a generalized nonlinear least‐squares problem with the target function being the weighted sum of squares of the differences between calculated and “experimental” (i.e., estimated from protein‐crystal data) angular, radial‐angular, and radial pair correlation functions, as well as contact free energies. A set of 195 high‐resolution nonhomologous structures from the Protein Data Bank was used to calculate the “experimental” values. The contact free energies were scaled by the slope of the correlation line between side‐chain hydrophobicities, calculated from the contact free energies, and those determined by Fauchere and Pliška from the partition coefficients of amino acids between water and n‐octanol. The methylene group served to define the reference contact free energy corresponding to that between the glycine methylene groups of backbone residues. Statistical analysis of the goodness of fit revealed that the Gay‐Berne‐Vorobjev anisotropic potential fits best to the experimental radial and angular correlation functions and contact free energies and therefore represents the free‐energy surface of side‐chain‐side‐chain interactions most accurately. Thus, its choice for simulations of protein structure is probably the most appropriate. However, the use of simpler functional forms is recommended, if the speed of computations is an issue.

Cumulant-based expressions for the multibody terms for the correlation between local and electrostatic interactions in the united-residue force field

A general method to derive site-site or united-residue potentials is presented. The basic principle of the method is the separation of the degrees of freedom of a system into the primary and secondary ones. The primary degrees of freedom describe the basic features of the system, while the secondary ones are averaged over when calculating the potential of mean force, which is hereafter referred to as the restricted free energy (RFE) function. The RFE can be factored into one-, two-, and multibody terms, using the cluster-cumulant expansion of Kubo. These factors can be assigned the functional forms of the corresponding lowest-order nonzero generalized cumulants, which can, in most cases, be evaluated analytically, after making some simplifying assumptions. This procedure to derive coarse-grain force fields is very valuable when applied to multibody terms, whose functional forms are hard to deduce in another way (e.g., from structural databases). After the functional forms have been derived, they can be parametrized based on the RFE surfaces of model systems obtained from all-atom models or on the statistics derived from structural databases. The approach has been applied to our united-residue force field for proteins. Analytical expressions were derived for the multibody terms pertaining to the correlation between local and electrostatic interactions within the polypeptide backbone; these expressions correspond to up to sixth-order terms in the cumulant expansion of the RFE. These expressions were subsequently parametrized by fitting to the RFEs of selected peptide fragments, calculated with the empirical conformational energy program for peptides force field. The new multibody terms enable not only the heretofore predictable α-helical segments, but also regular β-sheets, to form as the lowest-energy structures, as assessed by test calculations on a model helical protein A, as well as a model 20-residue polypeptide (betanova); the latter was not possible without introducing these new terms.

Recent improvements in prediction of protein structure by global optimization of a potential energy function.

Recent improvements of a hierarchical ab initio or de novo approach for predicting both α and β structures of proteins are described. The united-residue energy function used in this procedure includes multibody interactions from a cumulant expansion of the free energy of polypeptide chains, with their relative weights determined by Z-score optimization. The critical initial stage of the hierarchical procedure involves a search of conformational space by the conformational space annealing (CSA) method, followed by optimization of an all-atom model. The procedure was assessed in a recent blind test of protein structure prediction (CASP4). The resulting lowest-energy structures of the target proteins (ranging in size from 70 to 244 residues) agreed with the experimental structures in many respects. The entire experimental structure of a cyclic α-helical protein of 70 residues was predicted to within 4.3 Å α-carbon (Cα) rms deviation (rmsd) whereas, for other α-helical proteins, fragments of roughly 60 residues were predicted to within 6.0 Å Cα rmsd. Whereas β structures can now be predicted with the new procedure, the success rate for α/β- and β-proteins is lower than that for α-proteins at present. For the β portions of α/β structures, the Cα rmsds are less than 6.0 Å for contiguous fragments of 30–40 residues; for one target, three fragments (of length 10, 23, and 28 residues, respectively) formed a compact part of the tertiary structure with a Cα rmsd less than 6.0 Å. Overall, these results constitute an important step toward the ab initio prediction of protein structure solely from the amino acid sequence.

Principal Component Analysis for Protein Folding Dynamics

Protein folding is considered here by studying the dynamics of the folding of the triple beta-strand WW domain from the Formin-binding protein 28. Starting from the unfolded state and ending either in the native or nonnative conformational states, trajectories are generated with the coarse-grained united residue (UNRES) force field. The effectiveness of principal components analysis (PCA), an already established mathematical technique for finding global, correlated motions in atomic simulations of proteins, is evaluated here for coarse-grained trajectories. The problems related to PCA and their solutions are discussed. The folding and nonfolding of proteins are examined with free-energy landscapes. Detailed analyses of many folding and nonfolding trajectories at different temperatures show that PCA is very efficient for characterizing the general folding and nonfolding features of proteins. It is shown that the first principal component captures and describes in detail the dynamics of a system. Anomalous diffusion in the folding/nonfolding dynamics is examined by the mean-square displacement (MSD) and the fractional diffusion and fractional kinetic equations. The collisionless (or ballistic) behavior of a polypeptide undergoing Brownian motion along the first few principal components is accounted for.

A united-residue force field for off-lattice protein-structure simulations. II. Parameterization of short-range interactions and determination of weights of energy terms by Z-score optimization

Continuing our work on the determination of an off‐lattice united‐residue force field for protein‐structure simulations, we determined and parameterized appropriate functional forms for the local‐interaction terms, corresponding to the rotation about the virtual bonds (Utor), the bending of virtual‐bond angles (Ub), and the energy of different rotameric states of side chains (Urot). These terms were determined by applying the Boltzmann principle to the distributions of virtual‐bond torsional and virtual‐bond angles and side‐chain rotameric states, respectively, calculated from a data base of 195 high‐resolution nonhomologous proteins. The complete energy function was constructed by combining the individual energy terms with appropriate weights. The weights were determined by optimizing the so‐called Z‐score value (which is the normalized difference between the energy of the native structure and the mean energy of non‐native structures) of the histidine‐containing phosphocarrier protein from Streptococcus faecalis (1PTF; an 88‐residue α + β protein). To accomplish this, a database of Cα patterns was created using high‐resolution nonhomologous protein structures from the Protein Data Bank, and the distributions of energy components of 1PTF were obtained by threading its sequence through ∼500 randomly chosen Cα‐patterns from the X‐ray structures in the PDB, followed by energy minimization, with the energy function incorporating initially guessed weights. The resulting minimized energies were used to optimize the Z‐score value of 1PTF as a function of the weights of the various energy terms, and the new weights were used to generate new energy‐component distributions. The process was iterated, until the weights used to generate the distributions and the optimized weights were self‐consistent. The potential function with the weights of the various energy terms obtained by optimizing the Z‐score value for 1PTF was found to locate the native structures of other test proteins (within an average RMS deviation of 3 Å): calcium‐binding protein (4ICB), ubiquitin (1UBQ), α‐spectrin (1SHG), major cold‐shock protein (1MJC), and cytochrome b5 (3B5C) (which included α and β structures) as distinctively lowest in energy in similar threading experiments.

United-residue force field for off-lattice protein-structure simulations: III. Origin of backbone hydrogen-bonding cooperativity in united-residue potentials

Based on the dipole model of peptide groups developed in our earlier work [Liwo et al., Prot. Sci., 2, 1697 (1993)], a cumulant expansion of the average free energy of the system of freely rotating peptide‐group dipoles tethered to a fixed α‐carbon trace is derived. A graphical approach is presented to find all nonvanishing terms in the cumulants. In particular, analytical expressions for three‐ and four‐body (correlation) terms in the averaged interaction potential of united peptide groups are derived. These expressions are similar to the cooperative forces in hydrogen bonding introduced by Koliński and Skolnick [J. Chem. Phys., 97, 9412 (1992)]. The cooperativity arises here naturally from the higher order terms in the power‐series expansion (in the inverse of the temperature) for the average energy. Test calculations have shown that addition of the derived four‐body term to the statistical united‐residue potential of our earlier work [Liwo et al., J. Comput. Chem., 18, 849, 874 (1997)] greatly improves its performance in folding poly‐l‐alanine into an α‐helix. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 259–276, 1998

Computational techniques for efficient conformational sampling of proteins

In this review, we summarize the computational methods for sampling the conformational space of biomacromolecules. We discuss the methods applicable to find only lowest energy conformations (global minimization of the potential-energy function) and to generate canonical ensembles (canonical Monte Carlo method and canonical molecular dynamics method and their extensions). Special attention is devoted to the use of coarse-grained models that enable simulations to be enhanced by several orders of magnitude.

A method for optimizing potential-energy functions by a hierarchical design of the potential-energy landscape: Application to the UNRES force field

A method for optimizing potential-energy functions of proteins is proposed. The method assumes a hierarchical structure of the energy landscape, which means that the energy decreases as the number of native-like elements in a structure increases, being lowest for structures from the native family and highest for structures with no native-like element. A level of the hierarchy is defined as a family of structures with the same number of native-like elements (or degree of native likeness). Optimization of a potential-energy function is aimed at achieving such a hierarchical structure of the energy landscape by forcing appropriate free-energy gaps between hierarchy levels to place their energies in ascending order. This procedure is different from methods developed thus far, in which the energy gap and/or the Z score between the native structure and all non-native structures are maximized, regardless of the degree of native likeness of the non-native structures. The advantage of this approach lies in reducing the number of structures with decreasing energy, which should ensure the searchability of the potential. The method was tested on two proteins, PDB ID codes 1FSD and 1IGD, with an off-lattice united-residue force field. For 1FSD, the search of the conformational space with the use of the conformational space annealing method and the newly optimized potential-energy function found the native structure very quickly, as opposed to the potential-energy functions obtained by former optimization methods. After even incomplete optimization, the force field obtained by using 1IGD located the native-like structures of two peptides, 1FSD and betanova (a designed three-stranded β-sheet peptide), as the lowest-energy conformations, whereas for the 46-residue N-terminal fragment of staphylococcal protein A, the native-like conformation was the second-lowest-energy conformation and had an energy 2 kcal/mol above that of the lowest-energy structure.

A general method for the determination of the stoichiometry of unknown species in multicomponent systems from physicochemical measurements

Abstract A general algorithm is presented for the determination of the stoichiometry of unknown species present in chemical equilibrium systems. It is based on processing data from physicochemical measurements. The algorithm consists of two parts. In the first stage the stoichiometric coefficients are treated as real numbers, which are evaluated together with other parameters by a nonlinear least-squares method. Then integral models lying in the neighbourhood of the resultant “nonintegral” stoichiometry are examined, and that model which best fits the experimental data is chosen. Usually the stoichiometric coefficients obtained in the first step are close to integral numbers, which values are then established unambiguously in the search step. In the analysis of the experimental data errors in the composition of the solutions (e.g. in the titrant volumes) are taken into account as well as the errors in the measured quantities [e.g. in electromotive force (e.m.f.) or absorbance]. Numerical examples of calculations are presented for both simulated and real potentiometric spectrophotometric, and osmometric data. Problems of determinability of an equilibrium model, of deciding on the most probable model, and of estimating the number of species present are discussed.

Calculation of protein conformation by global optimization of a potential energy function.

A novel hierarchical approach to protein folding has been applied to compute the unknown structures of seven target proteins provided by CASP3. The approach is based exclusively on the global optimization of a potential energy function for a united‐residue model by conformational space annealing, followed by energy refinement using an all‐atom potential. Comparison of the submitted models for five globular proteins with the experimental structures shows that the conformations of large fragments (∼60 aa) were predicted with rmsds of 4.2–6.8 Å for the Cα atoms. Our lowest‐energy models for targets T0056 and T0061 were particularly successful, producing the correct fold of approximately 52% and 80% of the structures, respectively. These results support the thermodynamic hypothesis that protein structure can be computed solely by global optimization of a potential energy function for a given amino acid sequence. Proteins Suppl 1999;3:204–208.