Matthew R. Pincus

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

MOLECULAR MODELING OF MAMMALIAN CYTOCHROME P450s

Abstract. The cytochrome P450s are a superfamily of hemoprotein enzymes responsible for the metabolism of a wide variety of xenobiotic and endogenous compounds. The individual P450s exhibit unique substrate specificity and stereoselectivity profiles which reflect corresponding differences in primary sequence and tertiary structure. In the absence of an experimental structure, models for mammalian P450s have been generated by their homology with bacterial P450s of known structure. The rather low sequence identity between target and template proteins renders P450 modeling a challenging task. However, the substrate recognition properties of several P450s are consistent with recently developed working models. This review summarizes the major concepts and current approaches of molecular modeling of P450s.

Inhibition of human cytochrome P450 1A2 by flavones: a molecular modeling study.

Cytochrome P450 1A2 metabolizes a number of important drugs, procarcinogens, and endogenous compounds. Several flavones, a class of phytochemicals consumed in the human diet, have been shown to differentially inhibit human P450 1A2-mediated methoxyresorufin demethylase. A molecular model of this P450 was constructed in order to elucidate the molecular basis of the P450-flavone interaction. Flavone and its 3,5,7-trihydroxy and 3,5,7-trimethoxy derivatives were docked into the active site to assess their mode of binding. The site is hydrophobic and includes several residues that hydrogen bond with substituents on the flavone nucleus. The binding interactions of these flavones in the modeled active side are consistent with their relative inhibitory potentials, namely 3,5,7-trihydroxylflavone > flavone >3,5,7-trimethoxylflavone, toward P450 1A2-mediated methoxyresorufin demethylation.

Molecular Modeling of Cytochrome P450 2B1: Mode of Membrane Insertion and Substrate Specificity

A molecular model of a mammalian membrane-bound cytochrome P450, rat P450 2B1, was constructed in order to elucidate its mode of attachment to the endoplasmic reticulum and the structural basis of substrate specificity. The model was primarily derived from the structure of P450BM-3, which as a class II P450 is the most functionally similar P450 of known structure. However, model development was also guided by the conserved core regions of P450cam and P450terp. To optimally align the P450 2B1 and P450BM-3 sequences, multiple alignment was performed using sequences of five P450s in the II family, followed by minor adjustments on the basis of secondary structure predictions. The resulting P450 2B1 homology model structure was refined by molecular dynamics heating, equilibration, simulation, and energy minimization. The model suggests that the F–G loop serves as both a hydrophobic membrane anchor and entrance channel for hydrophobic substrates from the membrane to the P450 active site. To assess the mode of substrate binding, benzphetamine, testosterone, and benzo[a]pyrene were docked into the active site. The hydrophobic substrate-binding pocket is consistent with the preferences of this P450 toward hydrophobic substrates, while the presence of an acidic Glu-105 in this pocket is consistent with the preference of this P450 for the cationic substrate benzphetamine. This model is thus consistent with several known experimental properties of this P450, such as membrane attachment and substrate selectivity.

Comparison of the low energy conformations of an oncogenic and a non-oncogenic p21 protein, neither of which binds GTP or GDP

Oncogenic p21 protein, encoded by theras-oncogene, that causes malignant transformation of normal cells and many human tumors, is almost identical in sequence to its normal protooncogene-encoded counterpart protein, except for the substitution of arbitrary amino acids for the normally occurring amino acids at critical positions such as Gly 12 and Gin 61. Since p21 is normally activated by the binding of GTP in place of GDP, it has been postulated that oncogenic forms must retain bound GTP for prolonged time periods. However, two multiply substituted p21 proteins have been cloned, neither of which binds GDP or GTP. One of these mutant proteins with Val for Gly 10, Arg for Gly 12, and Thr for Ala 59 causes cell transformation, while the other, similar protein with Gly 10, Arg 12, Val for Gly 13 and Thr 59 does not transform cells. To define the critical conformational changes that occur in the p21 protein that cause it to become oncogenic, we have calculated the low energy conformations of the two multiply substituted mutant p21 proteins using a new adaptation of the electrostatically driven Monte Carlo (EDMC) technique, based on the program ECEPP. We have used this method to explore the conformational space available to both proteins and to compute the average structures for both using statistical mechanical averaging. Comparison of the average structures allows us to detect the major differences in conformation between the two proteins. Starting structures for each protein were calculated using the recently deposited x-ray crystal coordinates for the p21 protein, that was energy-refined using ECEPP, and then perturbed using the EDMC method to compute its average structure. The specific amino acid substitutions for both proteins were then generated into the lowest energy structure generated by this procedure, subjected to energy minimization and then to full EDMC perturbations. We find that both mutant proteins exhibit major differences in conformation in specific regions, viz., residues 35–47, 55–78, 81–93, 96–110, 115–126, and 123–134, compared with the EDMC-refined x-ray structure of the wild-type protein. These regions have been found to be the most flexible in the p21 protein bound to GDP from prior molecular dynamics calculations (Dykeset al., 1993). Comparison of the EDMC-average structure of the transforming mutant with that of the nontransforming mutant reveals major structural differences at residues 10–16, 32–40, and 60–68. These structural differences appear to be the ones that are critical in activation of the p21 protein. Analysis of the correlated motions of the different regions of the two mutant proteins reveals that changes in the conformation of regions in the carboxyl half of the protein are caused by changes in conformation around residues 10–16 and are transmitted by means of residues around Gln 61. The latter region therefore constitutes a “molecular switch” unit, in agreement with conclusions from prior work.

Conformation of the transmembrane domain of the c-erbB-2 oncogene-encoded protein in its monomeric and dimeric states

The human c-erbB-2 oncogene is homologous to the ratneu oncogene, both encoding transmembrane growth factor receptors. Overexpression and point mutations in the transmembrane domain of the encoded proteins in both cases have been implicated in cell transformation and carcinogenesis. In the case of theneu protein, it has been proposed that these effects are mediated by conformational preferences for anα-helix in the transmembrane domain, which facilitates receptor dimerization, an important step in the signal transduction process. To examine whether this is the case for c-erbB-2 as well, we have used conformational energy analysis to determine the preferred three-dimensional structures for the transmembrane domain of the c-erbB-2 protein from residues 650 to 668 with Val (nontransforming) and Glu (transforming) at position 659. The global minimum energy conformation for the Val-659 peptide from the normal, nontransforming protein was found to contain several bends, whereas the global minimum energy conformation for Glu-659 peptide from the mutant, transforming protein was found to beα-helical. Thus, the difference in conformational preferences for these transmembrane domains may explain the difference in transforming ability of these proteins. The presence of higher-energyα-helical conformations for the transmembrane domain from the normal Val-659 protein may provide an explanation for the presence of a transforming effect from overexpression of c-erbB-2. In addition, docking of the oncogenic sequences in theirα-helical and bend conformations shows that the all-α-helical dimer is clearly favored energetically over the bend dimer.

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.

Common Conformational Effects of p53 Mutations

The tumor suppressor gene p53 has been identified as the most frequent target of genetic alterations in human cancers. Most of these mutations occur in highly conserved regions in the DNA-binding core domain of the p53 protein, suggesting that the amino acid residues in these regions are critical for maintaining normal p53 structure and function. We previously used molecular dynamics calculations to demonstrate that several amino acid substitutions in these regions that are induced by environmental carcinogens and found in human tumors produce certain common conformational changes in the mutant proteins that differ substantially from the wild-type structure. In order to determine whether these conformational changes are consistent for other p53 mutants, we have now used molecular dynamics to determine the structure of the DNA-binding core domain of seven other environmentally induced, cancer-related p53 mutants, namely His 175, Asp 245, Asn 245, Trp 248, Met 249, Ser 278, and Lys 286. The results indicate that all of these mutants differ substantially from the wild-type structure in certain discrete regions and that some of these conformational changes are similar for these mutants as well as those determined previously. The changes are also consistent with experimental evidence for alterations in structure in p53 mutants determined by epitope detectability using monoclonal antibodies directed against these regions of predicted conformational change.