The minimized dead-end elimination criterion and its application to protein redesign in a hybrid scoring and search algorithm for computing partition functions over molecular ensembles

One of the main challenges for protein redesign is the efficient evaluation of a combinatorial number of candidate structures. The modeling of protein flexibility, typically by using a rotamer library of commonly-observed low-energy side-chain conformations, further increases the complexity of the redesign problem. A dominant algorithm for protein redesign is dead-end elimination (DEE), which prunes the majority of candidate conformations by eliminating rigid rotamers that provably are not part of the global minimum energy conformation (GMEC). The identified GMEC consists of rigid rotamers (i.e., rotamers that have not been energy-minimized) and is thus referred to as the rigid-GMEC. As a postprocessing step, the conformations that survive DEE may be energy-minimized. When energy minimization is performed after pruning with DEE, the combined protein design process becomes heuristic, and is no longer provably accurate: a conformation that is pruned using rigid-rotamer energies may subsequently minimize to a lower energy than the rigid-GMEC. That is, the rigid-GMEC and the conformation with the lowest energy among all energy-minimized conformations (the minimized-GMEC) are likely to be different. While the traditional DEE algorithm succeeds in not pruning rotamers that are part of the rigid-GMEC, it makes no guarantees regarding the identification of the minimized-GMEC. In this paper we derive a novel, provable, and efficient DEE-like algorithm, called minimized-DEE (MinDEE), that guarantees that rotamers belonging to the minimized-GMEC will not be pruned, while still pruning a combinatorial number of conformations. We show that MinDEE is useful not only in identifying the minimized-GMEC, but also as a filter in an ensemble-based scoring and search algorithm for protein redesign that exploits energy-minimized conformations. We compare our results both to our previous computational predictions of protein designs and to biological activity assays of predicted protein mutants. Our provable and efficient minimized-DEE algorithm is applicable in protein redesign, protein-ligand binding prediction, and computer-aided drug design.     

Preprocessing of rotamers for protein design calculations

We have developed a process that significantly reduces the number of rotamers in computational protein design calculations. This process, which we call Vegas, results in dramatic computational performance increases when used with algorithms based on the dead-end elimination (DEE) theorem. Vegas estimates the energy of each rotamer at each position by fixing each rotamer in turn and utilizing various search algorithms to optimize the remaining positions. Algorithms used for this context specific optimization can include Monte Carlo, self-consistent mean field, and the evaluation of an expression that generates a lower bound energy for the fixed rotamer. Rotamers with energies above a user-defined cutoff value are eliminated. We found that using Vegas to preprocess rotamers significantly reduced the calculation time of subsequent DEE-based algorithms while retaining the global minimum energy conformation. For a full boundary design of a 51 amino acid fragment of engrailed homeodomain, the total calculation time was reduced by 12-fold.     

Generate: a program for 3-D structure generation and conformational analysis of peptides and peptidomimetics

The program Generate, aimed at generating 3-D structures for peptides and peptidomimetics, is presented. The algorithm is based on a build-up procedure, using a library of conformations of amino acid residues. This library is built from conformational analysis of amino acids placed in a di- or tripeptide environment to mimic the surroundings of the amino acid in a true peptide, considering different positions of the residue in the peptide chain (peptidyl fragment, NH(+)(3)-terminus or COO(-)-terminus). Cis-trans isomerism in the amide bonds is taken into account by construction of rotamer libraries for different isomers. Water solvation is included through the GB/SA model. New amino acid residues can easily be added to the libraries, making it possible to generate conformations of peptidomimetics.     

Molecular modelling prediction of ligand binding site flexibility

We have investigated the efficacy of generating multiple sidechain conformations using a rotamer library in order to find the experimentally observed ligand binding site conformation of a protein in the presence of a bound ligand. We made use of a recently published algorithm that performs an exhaustive conformational search using a rotamer library to enumerate all possible sidechain conformations in a binding site. This approach was applied to a dataset of proteins whose structures were determined by X-ray and NMR methods. All chosen proteins had two or more structures, generally involving different bound ligands. By taking one of these structures as a reference, we were able in most cases to successfully reproduce the experimentally determined conformations of the other structures, as well as to suggest alternative low-energy conformations of the binding site. In those few cases where this procedure failed, we observed that the bound ligand had induced a high-energy conformation of the binding site. These results suggest that for most proteins that exhibit limited backbone motion, ligands tend to bind to low energy conformations of their binding sites. Our results also reveal that it is possible in most cases to use a rotamer search-based approach to predict alternative low-energy protein binding site conformations that can be used by different ligands. This opens the possibility of incorporating alternative binding site conformations to improve the efficacy of docking and structure-based drug design algorithms.     

Importance of chirality and reduced flexibility of protein side chains: a study with square and tetrahedral lattice models

Side chains of amino acid residues are the determining factor that distinguishes proteins from other unstable chain polymers. In simple models they are often represented implicitly (e.g., by spin states) or simplified as one atom. Here we study side chain effects using two-dimensional square lattice and three-dimensional tetrahedral lattice models, with explicitly constructed side chains formed by two atoms of different chirality and flexibility. We distinguish effects due to chirality and effects due to side chain flexibilities, since residues in proteins are L residues, and their side chains adopt different rotameric states. For short chains, we enumerate exhaustively all possible conformations. For long chains, we sample effectively rare events such as compact conformations and obtain complete pictures of ensemble properties of conformations of these models at all compactness region. This is made possible by using sequential Monte Carlo techniques based on chain growth method. Our results show that both chirality and reduced side chain flexibility lower the folding entropy significantly for globally compact conformations, suggesting that they are important properties of residues to ensure fast folding and stable native structure. This corresponds well with our finding that natural amino acid residues have reduced effective flexibility, as evidenced by statistical analysis of rotamer libraries and side chain rotatable bonds. We further develop a method calculating the exact side chain entropy for a given backbone structure. We show that simple rotamer counting underestimates side chain entropy significantly for both extended and near maximally compact conformations. We find that side chain entropy does not always correlate well with main chain packing. With explicit side chains, extended backbones do not have the largest side chain entropy. Among compact backbones with maximum side chain entropy, helical structures emerge as the dominating configurations. Our results suggest that side chain entropy may be an important factor contributing to the formation of alpha helices for compact conformations.     

Rotamer optimization for protein design through MAP estimation and problem‐size reduction

The search for the global minimum energy conformation (GMEC) of protein side chains is an important computational challenge in protein structure prediction and design. Using rotamer models, the problem is formulated as a NP‐hard optimization problem. Dead‐end elimination (DEE) methods combined with systematic A* search (DEE/A*) has proven useful, but may not be strong enough as we attempt to solve protein design problems where a large number of similar rotamers is eligible and the network of interactions between residues is dense. In this work, we present an exact solution method, named BroMAP (branch‐and‐bound rotamer optimization using MAP estimation), for such protein design problems. The design goal of BroMAP is to be able to expand smaller search trees than conventional branch‐and‐bound methods while performing only a moderate amount of computation in each node, thereby reducing the total running time. To achieve that, BroMAP attempts reduction of the problem size within each node through DEE and elimination by lower bounds from approximate maximum‐a‐posteriori (MAP) estimation. The lower bounds are also exploited in branching and subproblem selection for fast discovery of strong upper bounds. Our computational results show that BroMAP tends to be faster than DEE/A* for large protein design cases. BroMAP also solved cases that were not solved by DEE/A* within the maximum allowed time, and did not incur significant disadvantage for cases where DEE/A* performed well. Therefore, BroMAP is particularly applicable to large protein design problems where DEE/A* struggles and can also substitute for DEE/A* in general GMEC search. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

Approaches to solving the rigid receptor problem by identifying a minimal set of flexible residues during ligand docking

BACKGROUND: Using fixed receptor sites derived from high-resolution crystal structures in structure-based drug design does not properly account for ligand-induced enzyme conformational change and imparts a bias into the discovery and design of novel ligands. We sought to facilitate the design of improved drug leads by defining residues most likely to change conformation, and then defining a minimal manifold of possible conformations of a target site for drug design based on a small number of identified flexible residues. RESULTS: The crystal structure of thymidylate synthase from an important pathogenic target Pneumocystis carinii (PcTS) bound to its substrate and the inhibitor, BW1843U89, is reported here and reveals a new conformation with respect to the structure of PcTS bound to substrate and the more conventional antifolate inhibitor, CB3717. We developed an algorithm for determining which residues provide 'soft spots' in the protein, regions where conformational adaptation suggests possible modifications for a drug lead that may yield higher affinity. Remodeling the active site of thymidylate synthase with new conformations for only three residues that were identified with this algorithm yields scores for ligands that are compatible with experimental kinetic data. CONCLUSIONS: Based on the examination of many protein/ligand complexes, we develop an algorithm (SOFTSPOTS) for identifying regions of a protein target that are more likely to accommodate plastically to regions of a drug molecule. Using these indicators we develop a second algorithm (PLASTIC) that provides a minimal manifold of possible conformations of a protein target for drug design, reducing the bias in structure-based drug design imparted by structures of enzymes co-crystallized with inhibitors.     

W2466.48 Opens a Gate for a Continuous Intrinsic Water Pathway during Activation of the Adenosine A2A Receptor

The question how G‐protein‐coupled receptors transduce an extracellular signal by a sequence of transmembrane conformational transitions into an intracellular response remains to be solved at molecular detail. Herein, we use molecular dynamics simulations to reveal distinct conformational transitions of the adenosine A_2A receptor, and we found that the conserved W246^6.48 residue in transmembrane helix TM6 performs a key rotamer toggle switch. Agonist binding induces the sidechain of W246^6.48 to fluctuate between two distinct conformations enabling the diffusion of water molecules from the bulk into the center of the receptor. After passing the W246^6.48 gate, the internal water molecules induce another conserved residue, Y288^7.53, to switch to a distinct rotamer conformation establishing a continuous transmembrane water pathway. Further, structural changes of TM6 and TM7 induce local structural changes of the adjacent lipid bilayer.     

Toward deciphering the code to aminergic G protein-coupled receptor drug design

The trace amine-associated receptor 1 (TAAR(1)) is a biogenic amine G protein-coupled receptor (GPCR) that is potently activated by 3-iodothyronamine (1, T(1)AM) in vitro. Compound 1 is an endogenous derivative of the thyroid hormone thyroxine which rapidly induces hypothermia, anergia, and bradycardia when administered to mice. To explore the role of TAAR(1) in mediating the effects of 1, we rationally designed and synthesized rat TAAR(1) superagonists and lead antagonists using the rotamer toggle switch model of aminergic GPCR activation. The functional activity of a ligand is proposed to be correlated to its probable interactions with the rotamer switch residues; agonists allow the rotamer switch residues to toggle to their active conformation, whereas antagonists interfere with this conformational transition. These agonist and antagonist design principles provide a conceptual model for understanding the relationship between the molecular structure of a drug and its pharmacological properties.     

Statistically based reduced representation of amino acid side chains

Preferred conformations of amino acid side chains have been well established through statistically obtained rotamer libraries. Typically, these provide bond torsion angles allowing a side chain to be traced atom by atom. In cases where it is desirable to reduce the complexity of a protein representation or prediction, fixing all side-chain atoms may prove unwieldy. Therefore, we introduce a general parametrization to allow positions of representative atoms (in the present study, these are terminal atoms) to be predicted directly given backbone atom coordinates. Using a large, culled data set of amino acid residues from high-resolution protein crystal structures, anywhere from 1 to 7 preferred conformations were observed for each terminal atom of the non-glycine residues. Side-chain length from the backbone C(alpha) is one of the parameters determined for each conformation, which should itself be useful. Prediction of terminal atoms was then carried out for a second, nonredundant set of protein structures to validate the data set. Using four simple probabilistic approaches, the Monte Carlo style prediction of terminal atom locations given only backbone coordinates produced an average root mean-square deviation (RMSD) of approximately 3 A from the experimentally determined terminal atom positions. With prediction using conditional probabilities based on the side-chain chi(1) rotamer, this average RMSD was improved to 1.74 A. The observed terminal atom conformations therefore provide reasonable and potentially highly accurate representations of side-chain conformation, offering a viable alternative to existing all-atom rotamers for any case where reduction in protein model complexity, or in the amount of data to be handled, is desired. One application of this representation with strong potential is the prediction of charge density in proteins. This would likely be especially valuable on protein surfaces, where side chains are much less likely to be fixed in single rotamers. Prediction of ensembles of structures provides a method to determine the probability density of charge and atom location; such a prediction is demonstrated graphically.     

A normal mode-based geometric simulation approach for exploring biologically relevant conformational transitions in proteins

A three-step approach for multiscale modeling of protein conformational changes is presented that incorporates information about preferred directions of protein motions into a geometric simulation algorithm. The first two steps are based on a rigid cluster normal-mode analysis (RCNMA). Low-frequency normal modes are used in the third step (NMSim) to extend the recently introduced idea of constrained geometric simulations of diffusive motions in proteins by biasing backbone motions of the protein, whereas side-chain motions are biased toward favorable rotamer states. The generated structures are iteratively corrected regarding steric clashes and stereochemical constraint violations. The approach allows performing three simulation types: unbiased exploration of conformational space; pathway generation by a targeted simulation; and radius of gyration-guided simulation. When applied to a data set of proteins with experimentally observed conformational changes, conformational variabilities are reproduced very well for 4 out of 5 proteins that show domain motions, with correlation coefficients r > 0.70 and as high as r = 0.92 in the case of adenylate kinase. In 7 out of 8 cases, NMSim simulations starting from unbound structures are able to sample conformations that are similar (root-mean-square deviation = 1.0-3.1 ?) to ligand bound conformations. An NMSim generated pathway of conformational change of adenylate kinase correctly describes the sequence of domain closing. The NMSim approach is a computationally efficient alternative to molecular dynamics simulations for conformational sampling of proteins. The generated conformations and pathways of conformational transitions can serve as input to docking approaches or as starting points for more sophisticated sampling techniques.     

Joint neighbors approximation of macromolecular solvent accessible surface area

A new method for approximate analytical calculations of solvent accessible surface area (SASA) for arbitrary molecules and their gradients with respect to their atomic coordinates was developed. This method is based on the recursive procedure of pairwise joining of neighboring atoms. Unlike other available methods of approximate SASA calculations, the method has no empirical parameters, and therefore can be used with comparable accuracy in calculations of SASA in folded and unfolded conformations of macromolecules of any chemical nature. As shown by tests with globular proteins in folded conformations, average errors in absolute atomic surface area is around 1 A2, while for unfolded protein conformations it varies from 1.65 to 1.87 A2. Computational times of the method are comparable with those by GETAREA, one of the fastest exact analytical methods available today.     

Maintaining solvent accessible surface area under rotamer substitution for protein design

Although quantities derived from solvent accessible surface areas (SASA) are useful in many applications in protein design and structural biology, the computational cost of accurate SASA calculation makes SASA-based scores difficult to integrate into commonly used protein design methodologies. We demonstrate a method for maintaining accurate SASA during a Monte Carlo search of sequence and rotamer space for a fixed protein backbone. We extend the fast Le Grand and Merz algorithm (Le Grand and Merz, J Comput Chem, 14, 349), which discretizes the solvent accessible surface for each atom by placing dots on a sphere and combines Boolean masks to determine which dots are exposed. By replacing semigroup operations with group operations (from Boolean logic to counting dot coverage) we support SASA updates. Our algorithm takes time proportional to the number of atoms affected by rotamer substitution, rather than the number of atoms in the protein. For design simulations with a one hundred residue protein our approach is approximately 145 times faster than performing a Le Grand and Merz SASA calculation from scratch following each rotamer substitution. To demonstrate practical effectiveness, we optimize a SASA-based measure of protein packing in the complete redesign of a large set of proteins and protein-protein interfaces.     

Insights into the mobility of methyl-bearing side chains in proteins from (3)J(CC) and (3)J(CN) couplings

Side-chain dynamics in proteins can be characterized by the NMR measurement of (13)C and (2)H relaxation rates. Evaluation of the corresponding spectral densities limits the slowest motions that can be studied quantitatively to the time scale on which the overall molecular tumbling takes place. A different measure for the degree of side-chain order about the C(alpha)-C(beta) bond (chi(1) angle) can be derived from (3)J(C)(')(-)(C)(gamma) and (3)J(N)(-)(C)(gamma) couplings. These couplings can be measured at high accuracy, in particular for Thr, Ile, and Val residues. In conjunction with the known backbone structures of ubiquitin and the third IgG-binding domain of protein G, and an extensive set of (13)C-(1)H side-chain dipolar coupling measurements in oriented media, these (3)J couplings were used to parametrize empirical Karplus relationships for (3)J(C)(')(-)(C)(gamma) and (3)J(N)(-)(C)(gamma). These Karplus curves agree well with results from DFT calculations, including an unusual phase shift, which causes the maximum (3)J(CC) and (3)J(CN) couplings to occur for dihedral angles slightly smaller than 180 degrees, particularly noticeable in Thr residues. The new Karplus curves permit determination of rotamer populations for the chi(1) torsion angles. Similar rotamer populations can be derived from side-chain dipolar couplings. Conversion of these rotamer populations into generalized order parameters, S(J)(2) and S(D)(2), provides a view of side-chain dynamics that is complementary to that obtained from (13)C and (2)H relaxation. On average, results agree well with literature values for (2)H-relaxation-derived S(rel)(2) values in ubiquitin and HIV protease, but also identify a fraction of residues for which S(J,D)(2) < S(rel)(2). This indicates that some of the rotameric averaging occurs on a time scale too slow to be observable in traditional relaxation measurements.     

A search algorithm for fixed-composition protein design

We present a computational protein design algorithm for finding low-energy sequences of fixed amino acid composition. The search algorithms used in protein design typically do not restrict amino acid composition. However, the random energy model of Shakhnovich suggests that the use of fixed-composition sequences may circumvent defects in the modeling of the denatured state. Our algorithm, FC_FASTER, links fixed-composition versions of Monte Carlo and the FASTER algorithm. As proof of principle, FC_FASTER was tested on an experimentally validated, full-sequence design of the beta1 domain of protein G. For the wild-type composition, FC_FASTER found a lower energy sequence than the experimentally validated sequence. Also, for a different composition, FC_FASTER found the hypothetical lowest-energy sequence in 14 out of 32 trials.     

Peptide backbone reconstruction using dead-end elimination and a knowledge-based forcefield

A novel, yet simple and automated, protocol for reconstruction of complete peptide backbones from C(alpha) coordinates only is described, validated, and benchmarked. The described method collates a set of possible backbone conformations for each set of residue triads from a structural library derived from the PDB. The optimal permutation of these three residue segments of backbone conformations is determined using the dead-end elimination (DEE) algorithm. Putative conformations are evaluated using a pairwise-additive knowledge-based forcefield term and a fragment overlap term. The protocol described in this report is able to restore the full backbone coordinates to within 0.2-0.6 A of the actual crystal structure from C(alpha) coordinates only. In addition, it is insensitive to errors in the input C(alpha) coordinates with RMSDs of 3.0 A, and this is illustrated through application to deliberately distorted C(alpha) traces. The entire process, as described, is rapid, requiring of the order of a few minutes for a typical protein on a typical desktop PC. Approximations enable this to be reduced to a few seconds, although this is at the expense of prediction accuracy. This compares very favorably to previously published methods, being sufficiently fast for general use and being one of the most accurate methods. Because the method is not restricted to the reconstruction from only C(alpha) coordinates, reconstruction based on C(beta) coordinates is also demonstrated.     

The role of plastic beta-hairpin and weak hydrophobic core in the stability and unfolding of a full sequence design protein

In this study, the thermal stability of a designed alpha/beta protein FSD (full sequence design) was studied by explicit solvent simulations at three moderate temperatures, 273 K, 300 K, and 330 K. The average properties of the ten trajectories at each temperature were analyzed. The thermal unfolding, as judged by backbone root-mean-square deviation and percentage of native contacts, was displayed with increased sampling outside of the native basin as the temperature was raised. The positional fluctuation of the hairpin residues was significantly higher than that of the helix residues at all three temperatures. The hairpin segment displayed certain plasticity even at 273 K. Apart from the terminal residues, the highest fluctuation was shown in the turn residues 7-9. Secondary structure analysis manifested the structural heterogeneity of the hairpin segment. It was also revealed by the simulation that the hydrophobic core was vulnerable to thermal denaturation. Consistent with the experiment, the I7Y mutation in the double mutant FSD-EY (FSD with mutations Q1E and I7Y) dramatically increased the protein stability in the simulation, suggesting that the plasticity of the hairpin can be partially compensated by a stronger hydrophobic core. As for the unfolding pathway, the breathing of the hydrophobic core and the separation of the two secondary structure elements (alpha helix and beta hairpin) was the initiation step of the unfolding. The loss of global contacts from the separation further destabilized the hairpin structure and also led to the unwinding of the helix.     

Comparing three stochastic search algorithms for computational protein design: Monte Carlo,replica exchange Monte Carlo,and a multistart,steepest‐descent heuristic

Computational protein design depends on an energy function and an algorithm to search the sequence/conformation space. We compare three stochastic search algorithms: a heuristic, Monte Carlo (MC), and a Replica Exchange Monte Carlo method (REMC). The heuristic performs a steepest‐descent minimization starting from thousands of random starting points. The methods are applied to nine test proteins from three structural families, with a fixed backbone structure, a molecular mechanics energy function, and with 1, 5, 10, 20, 30, or all amino acids allowed to mutate. Results are compared to an exact, “Cost Function Network” method that identifies the global minimum energy conformation (GMEC) in favorable cases. The designed sequences accurately reproduce experimental sequences in the hydrophobic core. The heuristic and REMC agree closely and reproduce the GMEC when it is known, with a few exceptions. Plain MC performs well for most cases, occasionally departing from the GMEC by 3–4 kcal/mol. With REMC, the diversity of the sequences sampled agrees with exact enumeration where the latter is possible: up to 2 kcal/mol above the GMEC. Beyond, room temperature replicas sample sequences up to 10 kcal/mol above the GMEC, providing thermal averages and a solution to the inverse protein folding problem. © 2016 Wiley Periodicals, Inc.     

Chi1 torsion angle dynamics in proteins from dipolar couplings.

Experiments are presented for the measurement of one-bond carbon-proton dipolar coupling values at CH and CH2 ositions in 13C-labeled, approximately 50% fractionally deuterated proteins. 13Cbeta-1Hbeta dipolar couplings have been measured for 38 of 49 possible residues in the 63-amino-acid B1 domain of peptostreptococcal protein L in two aligning media and interpreted in the context of side-chain chi1 torsion angle dynamics. The beta protons for 18 of the 25 beta-methylene-containing amino acids for which dipolar data are available can be unambiguously stereoassigned, and for those residues which are best fit to a single rotamer model the chi(1) angles obtained deviate from crystal structure values by only 5.2 degrees (rmsd). The results for 11 other residues are significantly better fit by a model that assumes jumps between the three canonical (chi1 approximately -60 degrees, 60 degrees, 180 degrees ) rotamers. Relative populations of the rotamers are determined to within +/-6% uncertainty on average and correlate with dihedral angles observed for the three molecules in the crystal asymmetric unit. Entropic penalties for quenching chi1 jumps are considered for six mobile residues thought to be involved in binding to human immunoglobulins. This study demonstrates that dipolar couplings may be used to characterize both the conformation of static residues and side-chain motion with high precision.