Ryszard J. Wawak

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.

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

Efficient parallel algorithms in global optimization of potential energy functions for peptides, proteins, and crystals

Global optimization is playing an increasing role in physics, chemistry, and biophysical chemistry. One of the most important applications of global optimization is to find the global minima of the potential energy of molecules or molecular assemblies, such as crystals. The solution of this problem typically requires huge computational effort. Even the fastest processor available is not fast enough to carry out this kind of computation in real time for the problems of real interest, e.g., protein and crystal structure prediction. One way to circumvent this problem is to take advantage of massively parallel computing. In this paper, we provide several examples of parallel implementations of global optimization algorithms developed in our laboratory. All of these examples follow the master/worker approach. Most of the methods are parallelized on the algorithmic (coarse-grain) level and one example of fine-grain parallelism is given, in which the function evaluation itself is computationally expensive. All parallel algorithms were initially implemented on an IBM/SP2 (distributed-memory) machine. In all cases, however, message passing is handled through the standard Message Passing Interface (MPI); consequently the algorithms can also be implemented on any distributed- or shared-memory system that runs MPI. The efficiency of these implementations is discussed.

Comment on “Anti-cooperativity in hydrophobic interactions: A simulation study of spatial dependence of three-body effects and beyond” [J. Chem. Phys. 115, 1414 (2001)]

We address the criticism of our methodology for determination of the three-body cooperative terms in the potential of mean force (PMF) of the hydrophobic interaction of methane molecules in water [Czaplewski et al., Prot. Sci. 9, 1235 (2000)] expressed in the title paper of Shimizu and Chan, as well as their conclusion that hydrophobic association is predominantly anti-cooperative. We demonstrate that their reference two-methane PMF curve is subject to a systematic error, which invalidates their conclusions about the sign of the cooperative PMF.

Prediction of conformation of rat galanin in the presence and absence of water with the use of monte Carlo methods and the ECEPP/3 force field

The conformation of the 29-residue rat galanin neuropeptide was studied using the Monte Carlo with energy minimization (MCM) and electrostatically driven Monte Carlo (EDMC) methods. According to a previously elaborated procedure, the polypeptide chain was first treated in a united-residue approximation, in order to enable extensive exploration of the conformational space to be carried out (with the use of MCM), Then the low-energy united-residue conformations were converted to the all-atom representations, and EDMC simulations were carried out for the all-atom polypeptide chains, using the ECEPP/3 force field with hydration included. In order to estimate the effect of environment on galanin conformation, the low-energy conformations obtained as a result of these simulations were taken as starting structures for further EDMC runs that did not include hydration. The lowest-energy conformation obtained in aqueous solution calculations had a nonhelical N-terminal part packed against the nonpolar face of a residual helix that extended from Pro13 toward the C-terminus. One next lowest-energy structure was a nearly-all-helical conformation, but with a markedly higher energy. In contrast, all of the low-energy conformations in the absence of water were all-helical differing only by the extent to which the helix was kinked around Pro13. These results are in qualitative agreement with the available NMR and CD data of galanin in aqueous and nonaqueous solvents.

Gradient discontinuities in calculations involving molecular surface area

The free energy of solvation of a polypeptide or a protein can be expressed in terms of the accessible surface area of the molecule. Algorithms for energy minimization or for molecular dynamics, which involve the first derivatives of the energy, including the free energy of solvation, are commonly used in the conformational analysis of proteins. Discontinuities of the first derivatives, which occur in the accessible surface area and, hence, in the solvation energy, can cause serious numerical problems. In this paper, we describe all the situations in which the gradient of the molecular surface area becomes discontinuous.