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

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

No abstract

We describe the application of our recently proposed method of hierarchical optimization of the protein energy landscape to optimize our off-lattice united-residue (UNRES) force field using single training proteins. First, the IgG-binding domain from streptococcal protein G (PDB code 1IGD) was treated; earlier attempts to use this protein to optimize the force field by optimizing the energy gap and Z score between the nativelike and non-native structures failed. The structure of this protein consists of an N-terminal antiparallel β-hairpin, a middle α-helix, and a C-terminal antiparallel β-hairpin, these elements being referred to as β1, α2, and β3, respectively, with the two hairpins forming a parallel β-sheet packed against the α-helix. In our earlier study, one of these elements was assumed to form at level 1, two at level 2, and three at level 3, and higher levels corresponded to the proper packing of two or more elements. This approach resulted in a structure with the wrong packing of the β-sheet, and attempts at further optimization failed. We therefore tried a hierarchy scheme that corresponds to the sequence of folding events deduced from NMR experiments. In this scheme, level 1 corresponds to structures with either β3 or α2, level 2 to structures with both β3 and α2, level 3 to structures with β3, α2, and the N-terminal strand packed against α2 (with β1 still not fully formed), and level 4 to structures with β1, α2, and β3, with β3 being packed to β1, which also implies the packing of β1 and β3 against α2. This optimization was successful and resulted in a reasonably transferable force field that led to well-foldable proteins. This corroborates the conclusion from our model on-lattice studies (Liwo, A.; Arłukowicz, P.; Ołdziej, S.; Czaplewski, C.; Makowski, M.; Scheraga, H. A. J. Phys. Chem. B 2004, 108, 16918) that a proper design of the structural hierarchy is of crucial importance to the foldability with the resulting potential-energy function. Moreover, in the off-lattice approach, the design of the hierarchy also appears to be important to the success of the optimization procedure itself. The next series of calculations was carried out with the LysM domain from the E. coli 1E0G (α + β) protein, which is smaller than 1IGD. In this case, no experimental information about the folding pathway is available; nevertheless, we were able to deduce the appropriate hierarchy by a trial-and-error method. The resulting force field performed worse in tests on α + β- and β-proteins than that derived on the basis of 1IGD with a correct hierarchy, which suggests that the structure of the 1IGD protein encodes more structure-determining interactions common to all proteins than the 1E0G protein does. For 1E0G, we also attempted to carry out a single energy gap and Z-score optimization; this effort resulted in an unsearchable force field. (The nativelike structures could not be found by a global search, although they were the lowest in energy). Technical details of the method, including the maintenance of proper seconda...

The diabatic energy surfaces of terminally-blocked amino-acid residues (modeling the bulk of backbonelocal interaction patterns in proteins: trans-Møller-Plesset (MP2) ab initio level of theory with the 6-31G(d,p) basis set. The dihedral angles for rotation of the peptide groups about the C R -C R virtual-bond axes (λ (1) and λ (2) ) were used as variables; for the proline-derived peptide, only λ (2) was variable, because of the presence of the pyrrolidinering constraint. The resulting energy maps were compared with those obtained with the ECEPP/3 force field. On the basis of the MP2/6-31G(d,p) energy surfaces of terminally blocked single residues, the torsional potentials of mean force for rotation about the C R -C R virtual-bond axes and the double-torsional potentials of mean force for rotations about two consecutive virtual bond axes of all pairs and triplets of the prototypes of L-trans-amino acid residues were determined by numerical integration, fitted to one-and two-dimensional Fourier series in the virtual-bond-dihedral angles γ of the C R trace of a polypeptide chain, and implemented in the united-residue force field.

The 1,4-cycloaddition reactions of the singlet (1Δg) oxygen with s-cis-1,3-butadiene and benzene, with the formation of 3,6-dihydro[1,2]dioxin and 2,3-dioxabicyclo[2.2.2]octa-5,7-diene, respectively, were studied by means of the CAS MCSCF/CAS MCQDPT2 ab initio method with the 6-31G* basis set. In the case of butadiene the reaction was found to be exoenergetic and the product was found to have C 2 symmetry, with the peroxide moiety in the gauche configuration. In the case of benzene the reaction was found to be endoenergetic and the bicyclic product formed was found to have C 2 v symmetry, with the peroxide moiety in the syn configuration. Three possible reaction routes were studied:  (i) concerted cycloaddition, (ii) stepwise cheletropic cycloaddition with the formation of zwitterionic 2,5-dihydrofuran l-oxide as an intermediate, and (iii) stepwise cycloaddition with the formation of a linear intermediate. In the case of butadiene routes (i) and (ii) were excluded, because only second-order saddle points were found on the corresponding reaction pathways. The linear intermediate (I1) found in route (iii) has a biradical character, and its energy relative to that of the separate reactants is 4.1 kcal/mol. The dominant activation barrier corresponds to the transition structure T1 leading to I1 and amounts to 9.9 kcal/mol. The rearrangement of I1 to the product (P) involves only a minor activation barrier of 7.5 kcal/mol (relative to I1). In the case of benzene the reaction occurs in a concerted manner with a single transition structure having C 2 v symmetry; the activation barrier is 25.3 kcal/mol. This difference in binding mechanism can be explained in terms of the configuration of the peroxide moiety in the adduct.

The development of three physics-based energy functions (force fields), designed to simulate the restricted free energy of proteins of the R, β, and R/β structural classes, is described. Each force field corresponds to a particular weighting of the united-residue (UNRES) interactions defined in earlier work. [1][2][3][4][5][6] To find the optimal weights for the R, β, and R/β force fields, both the Z-score and energy gap of the native versus nonnative structures are minimized simultaneously for four benchmark proteins: 1pou (for the R force field), 1tpm (for the β force field), and 1bdd and betanova (for the R/β force field). The simultaneous minimization was carried out by using a novel Monte Carlo method, Vector Monte Carlo (VMC). For R-helical proteins, another weighting of the UNRES interactions (denoted as the R 0 force field) was developed; this fourth force field is described in a companion publication (Lee, J. et al. J. Phys. Chem. B 2001, 105, 7291). The structural implications of the final weights of the four force fields, i.e., the relative contributions of the various UNRES interactions to stabilizing common structural motifs of proteins, are analyzed. The R 0 , R, β, and R/β force fields were used in the CASP4 exercise for ab initio protein-structure prediction with reasonable success. Finally, using a simple model system it was shown that the VMC protocol does not require exhaustive sampling of medium-and high-energy structures in order to optimize the parameters of the potential energy adequately.

We determined the potentials of mean force corresponding to the bending offor use in our united-residue (UNRES) force field. The potentials were determined by integrating the ab initio energy surfaces of terminally-blocked glycine, alanine, and proline calculated in our earlier work at the MP2/6-31G(d, p) level (Ołdziej et al 2003 J. Phys. Chem. A 107 8035), where alanine represents all types of amino-acid residues except for glycine and proline. This resulted in 27 different free-energy surfaces. The potentials were found to depend both on θ and on the two virtual-bond dihedralangles (γ 1 and γ 2 ) whose axes are the edges of θ , as well as on the types of all three consecutive amino-acid residues whose C α atoms define the angle θ . The type of residue at the second and third position of a triad has a major influence on the potentials, while that in the first position is less important. Each surface was fitted well by a three-dimensional Fourier series in the trigonometric functions of multiplicities of θ/2, γ 1 , and γ 2 ; these analytical expressions can readily be implemented in the UNRES force field, thus replacing our earlier knowledge-based virtual-bond valence potentials. M Supplementary data are available from stacks.