Improved Pruning algorithms and Divide-and-Conquer strategies for Dead-End Elimination, with application to protein design

MOTIVATION: Structure-based protein redesign can help engineer proteins with desired novel function. Improving computational efficiency while still maintaining the accuracy of the design predictions has been a major goal for protein design algorithms. The combinatorial nature of protein design results both from allowing residue mutations and from the incorporation of protein side-chain flexibility. Under the assumption that a single conformation can model protein folding and binding, the goal of many algorithms is the identification of the Global Minimum Energy Conformation (GMEC). A dominant theorem for the identification of the GMEC is Dead-End Elimination (DEE). DEE-based algorithms have proven capable of eliminating the majority of candidate conformations, while guaranteeing that only rotamers not belonging to the GMEC are pruned. However, when the protein design process incorporates rotameric energy minimization, DEE is no longer provably-accurate. Hence, with energy minimization, the minimized-DEE (MinDEE) criterion must be used instead. RESULTS: In this paper, we present provably-accurate improvements to both the DEE and MinDEE criteria. We show that our novel enhancements result in a speedup of up to a factor of more than 1000 when applied in redesign for three different proteins: Gramicidin Synthetase A, plastocyanin, and protein G. AVAILABILITY: Contact authors for source code.     

Relating side-chain mobility in proteins to rotameric transitions: Insights from molecular dynamics simulations and NMR

The dynamic aspect of proteins is fundamental to understanding protein stability and function. One of the goals of NMR studies of side-chain dynamics in proteins is to relate spin relaxation rates to discrete conformational states and the timescales of interconversion between those states. Reported here is a physical analysis of side-chain dynamics that occur on a timescale commensurate with monitoring by ²H spin relaxation within methyl groups. Motivated by observations made from tens-of-nanoseconds long MD simulations on the small protein eglin c in explicit solvent, we propose a simple molecular mechanics-based model for the motions of side-chain methyl groups. By using a Boltzmann distribution within rotamers, and by considering the transitions between different rotamer states, the model semi-quantitatively correlates the population of rotamer states with ‘model-free’ order parameters typically fitted from NMR relaxation experiments. Two easy-to-use, analytical expressions are given for converting S²_axis’ values (order parameter for C–CH₃ bond) into side-chain rotamer populations. These predict that S²_axis’ values below 0.8 result from population of more than one rotameric state. The relations are shown to predict rotameric sampling with reasonable accuracy on the ps–ns timescale for eglin c and are validated for longer timescales on ubiquitin, for which side-chain residual dipolar coupling (RDC) data have been collected.     

Determinants of protein side-chain packing.

The problem of protein side-chain packing for a given backbone trace is investigated using 3 different prediction models. The first requires an exhaustive search of all possible combinations of side-chain conformers, using the dead-end elimination theorem. The second considers only side-chain-backbone interactions, whereas the third neglects side-chain-backbone interactions and instead keeps side-chain-side-chain interactions. Predictions of side-chain conformations for 11 proteins using all 3 models show that removal of side-chain-side-chain interactions does not cause a large decrease in the prediction accuracy, whereas the model having only side-chain-side-chain interactions still retains a significant level of accuracy. These results suggest that the 2 classes of interactions, side-chain-backbone and side-chain-side-chain, are consistent with each other and work concurrently to stabilize the native conformations. This is confirmed by analyses of energy spectra of the side-chain conformations derived from the fourth prediction model, the Independent model, which gives almost the same quality of the prediction as the dead-end elimination. The analyses indicate that the 2 classes of interactions simultaneously increase the energy difference between the native and nonnative conformations.     

Side-chain rotamers in protein α-α hairpins and a mechanism of their selection

Several regularities were observed for the distribution of side-chain rotamers in α-α hairpins of globular proteins. In left-turned α-α hairpins, most side chains adopt t rotamers in d-positions and g^? rotamers in g-positions. In right-turned α-α hairpins, most side-chains adopt g^? rotamers in a-positions and t rotamers in e-positions. Analysis of these regularities suggested that selection of the side-chain conformation in α-α hairpins depends on two main factors, the mode of α-helix packing and the positions of side chains in α-helices. The regularities were explained by the squeezing mechanism: interhelical interactions bring the α-helices close together so that the side chains are squeezed out of the helix-helix interface and adopt unique conformations.     

Analysis of side-chain rotamers in transmembrane proteins

We measured the frequency of side-chain rotamers in 14 alpha-helical and 16 beta-barrel membrane protein structures and found that the membrane environment considerably perturbs the rotamer frequencies compared to soluble proteins. Although there are limited experimental data, we found statistically significant changes in rotamer preferences depending on the residue environment. Rotamer distributions were influenced by whether the residues were lipid or protein facing, and whether the residues were found near the N- or C-terminus. Hydrogen-bonding interactions with the helical backbone perturbs the rotamer populations of Ser and His. Trp and Tyr favor side-chain conformations that allow their side chains to extend their polar atoms out of the membrane core, thereby aligning the side-chain polarity gradient with the polarity gradient of the membrane. Our results demonstrate how the membrane environment influences protein structures, providing information that will be useful in the structure prediction and design of transmembrane proteins.     

Rotamer strain as a determinant of protein structural specificity

We present direct evidence for a change in protein structural specificity due to hydrophobic core packing. High resolution structural analysis of a designed core variant of ubiquitin reveals that the protein is in slow exchange between two conformations. Examination of side-chain rotamers indicates that this dynamic response and the lower stability of the protein are coupled to greater strain and mobility in the core. The results suggest that manipulating the level of side-chain strain may be one way of fine tuning the stability and specificity of proteins.     

Side-chain flexibility in protein-ligand binding: the minimal rotation hypothesis

The goal of this work is to learn from nature about the magnitudes of side-chain motions that occur when proteins bind small organic molecules, and model these motions to improve the prediction of protein-ligand complexes. Following analysis of protein side-chain motions upon ligand binding in 63 complexes, we tested the ability of the docking tool SLIDE to model these motions without being restricted to rotameric transitions or deciding which side chains should be considered as flexible. The model tested is that side-chain conformational changes involving more atoms or larger rotations are likely to be more costly and less prevalent than small motions due to energy barriers between rotamers and the potential of large motions to cause new steric clashes. Accordingly, SLIDE adjusts the protein and ligand side groups as little as necessary to achieve steric complementarity. We tested the hypothesis that small motions are sufficient to achieve good dockings using 63 ligands and the apo structures of 20 different proteins and compared SLIDE side-chain rotations to those experimentally observed. None of these proteins undergoes major main-chain conformational change upon ligand binding, ensuring that side-chain flexibility modeling is not required to compensate for main-chain motions. Although more frugal in the number of side-chain rotations performed, this model substantially mimics the experimentally observed motions. Most side chains do not shift to a new rotamer, and small motions are both necessary and sufficient to predict the correct binding orientation and most protein-ligand interactions for the 20 proteins analyzed.     

Role of the structural context in selection of hydrophobic side-chain rotamers in a- and d-positions of alpha-helices

The study shows that in coiled-coil proteins the distribution of hydrophobic side-chain rotamers in a- and d-positions of alpha-helices is strongly dependent on the mutual arrangement of the a-helices. In coiled-coil dimers, where a-helices are packed "face-to-face", most side chains occupying a-positions adopt t-rotamers, and those in d-positions adopt g- -rotamers. In tetramers, where alpha-helices are packed "side-by-side", most side chains in a-positions adopt g- -rotamers and those in d-positions adopt t-rotamers. These features can be used for prediction of side-chain rotamers in protein modeling and design.     

Excluded volume in protein side-chain packing

The excluded volume occupied by protein side-chains and the requirement of high packing density in the protein interior should severely limit the number of side-chain conformations compatible with a given native backbone. To examine the relationship between side-chain geometry and side-chain packing, we use an all-atom Monte Carlo simulation to sample the large space of side-chain conformations. We study three models of excluded volume and use umbrella sampling to effectively explore the entire space. We find that while excluded volume constraints reduce the size of conformational space by many orders of magnitude, the number of allowed conformations is still large. An average repacked conformation has 20 % of its chi angles in a non-native state, a marked reduction from the expected 67 % in the absence of excluded volume. Interestingly, well-packed conformations with up to 50 % non-native chi angles exist. The repacked conformations have native packing density as measured by a standard Voronoi procedure. Entropy is distributed non-uniformly over positions, and we partially explain the observed distribution using rotamer probabilities derived from the Protein Data Bank database. In several cases, native rotamers that occur infrequently in the database are seen with high probability in our simulation, indicating that sequence-specific excluded volume interactions can stabilize rotamers that are rare for a given backbone. In spite of our finding that 65 % of the native rotamers and 85 % of chi(1) angles can be predicted correctly on the basis of excluded volume only, 95 % of positions can accommodate more than one rotamer in simulation. We estimate that, in order to quench the side-chain entropy observed in the presence of excluded volume interactions, other interactions (hydrophobic, polar, electrostatic) must provide an additional stabilization of at least 0.6 kT per residue in order to single out the native state.     

Cooperative charge fluctuations by migrating protons in globular proteins

A review of the hydrogen bonded network on the protein surface shows the presence of a charged complex system with parallel and competitive interactions, including ionizable side-chains, migrating protons, bound water and nearby backbone peptides. This system displays cooperative effects of dynamical nature, reviewed for lysozyme as a case. By increasing the water coverage of the protein powder, the bound water cluster exhibits a percolative transition, detectable by the onset of large water-assisted displacements of migrating protons, with a parallel emergence of protein mobility and biological function. By lowering the temperature, migrating protons exhibit a glassy dielectric relaxation in the low frequency range, pointing to a frustration by competing interactions similar to that observed in spin glasses and fragile glass forming liquids. The observation of these dissipative processes implies the occurrence of spontaneous charge fluctuations. A simplified model of the protein surface, where conformational and ionizable side-chain fluctuations are averaged out, is used to discuss the statistical physics of these cooperative effects. Some biological implications of this dynamical cooperativity for enzymatic activity are briefly suggested at the end.     

Role of the structural context in selection of hydrophobic side-chain rotamers in a- and d-positions of α-helices

It was demonstrated for the first time that the distribution of side-chain rotamers in the a-and d-positions of α-helices of coiled-coil (cc) proteins follows a certain trend, rather then being random. For instance, most side chains adopt t rotamers in the a-positions and g^? rotamers in the d-positions of helical dimers. Vice versa, most side chains adopt g^? rotamers in the a-positions and t rotamers in the d-positions of tetramers. It was concluded that selection of the side-chain rotamers depends on the packing of α-helices and, consequently, depends on the structural context.     

A graph-theory algorithm for rapid protein side-chain prediction

Fast and accurate side-chain conformation prediction is important for homology modeling, ab initio protein structure prediction, and protein design applications. Many methods have been presented, although only a few computer programs are publicly available. The SCWRL program is one such method and is widely used because of its speed, accuracy, and ease of use. A new algorithm for SCWRL is presented that uses results from graph theory to solve the combinatorial problem encountered in the side-chain prediction problem. In this method, side chains are represented as vertices in an undirected graph. Any two residues that have rotamers with nonzero interaction energies are considered to have an edge in the graph. The resulting graph can be partitioned into connected subgraphs with no edges between them. These subgraphs can in turn be broken into biconnected components, which are graphs that cannot be disconnected by removal of a single vertex. The combinatorial problem is reduced to finding the minimum energy of these small biconnected components and combining the results to identify the global minimum energy conformation. This algorithm is able to complete predictions on a set of 180 proteins with 34342 side chains in <7 min of computer time. The total chi(1) and chi(1 + 2) dihedral angle accuracies are 82.6% and 73.7% using a simple energy function based on the backbone-dependent rotamer library and a linear repulsive steric energy. The new algorithm will allow for use of SCWRL in more demanding applications such as sequence design and ab initio structure prediction, as well addition of a more complex energy function and conformational flexibility, leading to increased accuracy.     

Derivation and testing of pair potentials for protein folding. When is the quasichemical approximation correct?

Many existing derivations of knowledge-based statistical pair potentials invoke the quasichemical approximation to estimate the expected side-chain contact frequency if there were no amino acid pair-specific interactions. At first glance, the quasichemical approximation that treats the residues in a protein as being disconnected and expresses the side-chain contact probability as being proportional to the product of the mole fractions of the pair of residues would appear to be rather severe. To investigate the validity of this approximation, we introduce two new reference states in which no specific pair interactions between amino acids are allowed, but in which the connectivity of the protein chain is retained. The first estimates the expected number of side-chain contracts by treating the protein as a Gaussian random coil polymer. The second, more realistic reference state includes the effects of chain connectivity, secondary structure, and chain compactness by estimating the expected side-chain contrast probability by placing the sequence of interest in each member of a library of structures of comparable compactness to the native conformation. The side-chain contact maps are not allowed to readjust to the sequence of interest, i.e., the side chains cannot repack. This situation would hold rigorously if all amino acids were the same size. Both reference states effectively permit the factorization of the side-chain contact probability into sequence-dependent and structure-dependent terms. Then, because the sequence distribution of amino acids in proteins is random, the quasichemical approximation to each of these reference states is shown to be excellent. Thus, the range of validity of the quasichemical approximation is determined by the magnitude of the side-chain repacking term, which is, at present, unknown. Finally, the performance of these two sets of pair interaction potentials as well as side-chain contact fraction-based interaction scales is assessed by inverse folding tests both without and with allowing for gaps.     

Relationship between ion pair geometries and electrostatic strengths in proteins

The electrostatic free energy contribution of an ion pair in a protein depends on two factors, geometrical orientation of the side-chain charged groups with respect to each other and the structural context of the ion pair in the protein. Conformers in NMR ensembles enable studies of the relationship between geometry and electrostatic strengths of ion pairs, because the protein structural contexts are highly similar across different conformers. We have studied this relationship using a dataset of 22 unique ion pairs in 14 NMR conformer ensembles for 11 nonhomologous proteins. In different NMR conformers, the ion pairs are classified as salt bridges, nitrogen-oxygen (N-O) bridges and longer-range ion pairs on the basis of geometrical criteria. In salt bridges, centroids of the side-chain charged groups and at least a pair of side-chain nitrogen and oxygen atoms of the ion-pairing residues are within a 4 A distance. In N-O bridges, at least a pair of the side-chain nitrogen and oxygen atoms of the ion-pairing residues are within 4 A distance, but the distance between the side-chain charged group centroids is greater than 4 A. In the longer-range ion pairs, the side-chain charged group centroids as well as the side-chain nitrogen and oxygen atoms are more than 4 A apart. Continuum electrostatic calculations indicate that most of the ion pairs have stabilizing electrostatic contributions when their side-chain charged group centroids are within 5 A distance. Hence, most (approximately 92%) of the salt bridges and a majority (68%) of the N-O bridges are stabilizing. Most (approximately 89%) of the destabilizing ion pairs are the longer-range ion pairs. In the NMR conformer ensembles, the electrostatic interaction between side-chain charged groups of the ion-pairing residues is the strongest for salt bridges, considerably weaker for N-O bridges, and the weakest for longer-range ion pairs. These results suggest empirical rules for stabilizing electrostatic interactions in proteins.     

Protein design using continuous rotamers

Optimizing amino acid conformation and identity is a central problem in computational protein design. Protein design algorithms must allow realistic protein flexibility to occur during this optimization, or they may fail to find the best sequence with the lowest energy. Most design algorithms implement side-chain flexibility by allowing the side chains to move between a small set of discrete, low-energy states, which we call rigid rotamers. In this work we show that allowing continuous side-chain flexibility (which we call continuous rotamers) greatly improves protein flexibility modeling. We present a large-scale study that compares the sequences and best energy conformations in 69 protein-core redesigns using a rigid-rotamer model versus a continuous-rotamer model. We show that in nearly all of our redesigns the sequence found by the continuous-rotamer model is different and has a lower energy than the one found by the rigid-rotamer model. Moreover, the sequences found by the continuous-rotamer model are more similar to the native sequences. We then show that the seemingly easy solution of sampling more rigid rotamers within the continuous region is not a practical alternative to a continuous-rotamer model: at computationally feasible resolutions, using more rigid rotamers was never better than a continuous-rotamer model and almost always resulted in higher energies. Finally, we present a new protein design algorithm based on the dead-end elimination (DEE) algorithm, which we call iMinDEE, that makes the use of continuous rotamers feasible in larger systems. iMinDEE guarantees finding the optimal answer while pruning the search space with close to the same efficiency of DEE. Availability: Software is available under the Lesser GNU Public License v3. Contact the authors for source code.     

Breaking symmetry in protein dimers: designs and functions

Symmetry, and in particular point group symmetry, is generally the rule for the global arrangement between subunits in homodimeric and other oligomeric proteins. The structures of fragments of tropomyosin and bovine fibrinogen are recently published examples, however, of asymmetric interactions between chemically identical chains. Their departures from strict twofold symmetry are based on simple and generalizable chemical designs, but were not anticipated prior to their structure determinations. The current review aims to improve our understanding of the structural principles and functional consequences of asymmetric interactions in proteins. Here, a survey of >100 diverse homodimers has focused on the structures immediately adjacent to the twofold axis. Five regular frameworks in alpha-helical coiled coils and antiparallel beta-sheets accommodate many of the twofold symmetric axes. On the basis of these frameworks, certain sequence motifs can break symmetry in geometrically defined manners. In antiparallel beta-sheets, these asymmetries include register slips between strands of repeating residues and the adoption of different side-chain rotamers to avoid steric clashes of bulky residues. In parallel coiled coils, an axial stagger between the alpha-helices is produced by clusters of core alanines. Such simple designs lead to a basic understanding of the functions of diverse proteins. These functions include regulation of muscle contraction by tropomyosin, blood clot formation by fibrin, half-of-site reactivity of caspase-9, and adaptive protein recognition in the matrix metalloproteinase MMP9. Moreover, asymmetry between chemically identical subunits, by producing multiple equally stable conformations, leads to unique dynamic and self-assembly properties.     

Role of the protein side-chain fluctuations on the strength of pair-wise electrostatic interactions: comparing experimental with computed pK(a)s

The effect of the protein side-chain fluctuations on the strength of electrostatic interactions was studied. The effect was modeled on 7 different crystal structures on the same enzyme as well as on 20 molecular dynamics snapshot structures. It was shown that the side-chain flexibility affects predominantly the magnitude of the strong pair-wise interactions, that is, the pair-wise interaction among ion pairs, and practically does not affect the interactions with the rest of the protein. This was used to suggest a correction function that should be applied to the original pair-wise electrostatic interaction to mimic the effects of the fluctuations. The procedure is applied on three ion pairs identified in lysozyme. It was shown that sampling different side-chain rotamers and modifying the strength of the pair-wise interaction energies makes calculated pK(a)s less sensitive to the fluctuations of the structure and improves the prediction accuracy.     

SIRIUS. An automated method for the analysis of the preferred packing arrangements between protein groups

Automated methods have been developed to determine the preferred packing arrangement between interacting protein groups. A suite of FORTRAN programs, SIRIUS, is described for calculating and analysing the geometries of interacting protein groups using crystallographically derived atomic co-ordinates. The programs involved in calculating the geometries search for interacting pairs of protein groups using a distance criterion, and then calculate the spatial disposition and orientation of the pair. The second set of programs is devoted to analysis. This involves calculating the observed and expected distributions of the angles and assessing the statistical significance of the difference between the two. A database of the geometries of the 400 combinations of side-chain to side-chain interaction has been created. The approach used in analysing the geometrical information is illustrated here with specific examples of interactions between side-chains, peptide groups and particular types of atom. At the side-chain level, an analysis of aromatic-amino interactions, and the interactions of peptide carbonyl groups with arginine residues is presented. At the atomic level the analyses include the spatial disposition of oxygen atoms around tyrosine residues, and the frequency and type of contact between carbon, nitrogen and oxygen atoms. This information is currently being applied to the modelling of protein interactions.     

Trading accuracy for speed: A quantitative comparison of search algorithms in protein sequence design

Finding the minimum energy amino acid side-chain conformation is a fundamental problem in both homology modeling and protein design. To address this issue, numerous computational algorithms have been proposed. However, there have been few quantitative comparisons between methods and there is very little general understanding of the types of problems that are appropriate for each algorithm. Here, we study four common search techniques: Monte Carlo (MC) and Monte Carlo plus quench (MCQ); genetic algorithms (GA); self-consistent mean field (SCMF); and dead-end elimination (DEE). Both SCMF and DEE are deterministic, and if DEE converges, it is guaranteed that its solution is the global minimum energy conformation (GMEC). This provides a means to compare the accuracy of SCMF and the stochastic methods. For the side-chain placement calculations, we find that DEE rapidly converges to the GMEC in all the test cases. The other algorithms converge on significantly incorrect solutions; the average fraction of incorrect rotamers for SCMF is 0.12, GA 0.09, and MCQ 0.05. For the protein design calculations, design positions are progressively added to the side-chain placement calculation until the time required for DEE diverges sharply. As the complexity of the problem increases, the accuracy of each method is determined so that the results can be extrapolated into the region where DEE is no longer tractable. We find that both SCMF and MCQ perform reasonably well on core calculations (fraction amino acids incorrect is SCMF 0.07, MCQ 0.04), but fail considerably on the boundary (SCMF 0.28, MCQ 0.32) and surface calculations (SCMF 0.37, MCQ 0.44).     

Side-chain modeling with an optimized scoring function

Modeling side-chain conformations on a fixed protein backbone has a wide application in structure prediction and molecular design. Each effort in this field requires decisions about a rotamer set, scoring function, and search strategy. We have developed a new and simple scoring function, which operates on side-chain rotamers and consists of the following energy terms: contact surface, volume overlap, backbone dependency, electrostatic interactions, and desolvation energy. The weights of these energy terms were optimized to achieve the minimal average root mean square (rms) deviation between the lowest energy rotamer and real side-chain conformation on a training set of high-resolution protein structures. In the course of optimization, for every residue, its side chain was replaced by varying rotamers, whereas conformations for all other residues were kept as they appeared in the crystal structure. We obtained prediction accuracy of 90.4% for chi(1), 78.3% for chi(1 + 2), and 1.18 A overall rms deviation. Furthermore, the derived scoring function combined with a Monte Carlo search algorithm was used to place all side chains onto a protein backbone simultaneously. The average prediction accuracy was 87.9% for chi(1), 73.2% for chi(1 + 2), and 1.34 A rms deviation for 30 protein structures. Our approach was compared with available side-chain construction methods and showed improvement over the best among them: 4.4% for chi(1), 4.7% for chi(1 + 2), and 0.21 A for rms deviation. We hypothesize that the scoring function instead of the search strategy is the main obstacle in side-chain modeling. Additionally, we show that a more detailed rotamer library is expected to increase chi(1 + 2) prediction accuracy but may have little effect on chi(1) prediction accuracy.