De novo protein design. II. Plasticity in sequence space.

It is generally accepted that many different protein sequences have similar folded structures, and that there is a relatively high probability that a new sequence possesses a previously observed fold. An indirect consequence of this is that protein design should define the sequence space accessible to a given structure, rather than providing a single optimized sequence. We have recently developed a new approach for protein sequence design, which optimizes the complete sequence of a protein based on the knowledge of its backbone structure, its amino acid composition and a physical energy function including van der Waals interactions, electrostatics, and environment free energy. The specificity of the designed sequence for its template backbone is imposed by keeping the amino acid composition fixed. Here, we show that our procedure converges in sequence space, albeit not to the native sequence of the protein. We observe that while polar residues are well conserved in our designed sequences, non-polar amino acids at the surface of a protein are often replaced by polar residues. The designed sequences provide a multiple alignment of sequences that all adopt the same three-dimensional fold. This alignment is used to derive a profile matrix for chicken triose phosphate isomerase, TIM. The matrix is found to recognize significantly the native sequence for TIM, as well as closely related sequences. Possible application of this approach to protein fold recognition is discussed.     

Torsional relaxation for biopolymers.

We describe a method for making natural, physical movements in a chained polymer by sequentially adjusting a few neighboring torsion angles in the polymer backbone. In addition to being very fast and easy to implement, the method is also very general. It applies equally well to proteins and nucleic acids. This method is then used to design a local refinement procedure. We test the refinement procedure on the minimization of a simple energy function for proteins. The energy function has a simplified potential for hydrophobic interaction, a hydrogen-bond term, and a term for van der Waals interaction. There is considerable current interest in such simple energy functions for protein folding. When applied to refine structures found by a global search method, the refinement is able to produce large reduction in the hydrogen-bond term and the van der Waal term of the energy. We conclude that the method is particularly effective in finding good "packing" of residues in an initially compact conformation.     

Construction of side-chains in homology modelling. Application to the C-terminal lobe of rhizopuspepsin

A detailed and rule-based side-chain modelling procedure for globular proteins is presented. It uses the conformational information contained in a homologous (template) structure as a starting point and includes recipes for atom placement and for checking and improving the atomic positions. The scheme does not rely on intuitive judgements or visual examination of the model during construction or refinement. It comprises four stages; the first three are relatively simple and the fourth is more complex. In the first stage, initial conformations for as many atoms as possible are transferred from the template structure based on the application of trends reported previously. Second, these trends are used to correct poor van der Waals overlaps. Third, the remaining side-chains atoms (those for which no information is contained in the template) are placed by evaluating their rigid rotation, van der Waals surfaces. The fourth stage consists of a hierarchial series of conformational checks. They involve the evaluation of individual residue energies in the absence and presence of the rest of the protein relative to statistical trends observed in the template structure, the comparison of hydrogen-bonding patterns and side-chain accessibilities in the model and template and brief energy minimization followed by an evaluation of the rigid rotation potential energy surfaces of each side-chain. The checks pinpoint "incorrectly" modelled side-chains, suggest conformational changes and provide a means for determining the portions of the model that are likely to be correct and those likely to be in error. The procedure developed in the paper is tested by modelling the side-chains of the C-terminal lobe of the aspartyl proteinase rhizopuspepsin, using the rhizopuspepsin backbone and the homologous protein, penicillopepsin, as a template for the side-chains. The resultant model was compared to the high-resolution X-ray structure of rhizopuspepsin. Using penicillopepsin data only (stage I), 58% of the chi 1 dihedrals and 44% of the chi 2 dihedrals were modelled correctly. Once poor van der Waals overlaps had been corrected and all of the atoms had been placed (stages II and III), 86% of the chi 1 dihedrals and 75% of the chi 2 dihedrals were correct. After the refinement had been completed (stage IV), 92% of the chi 1 dihedrals and 81% of the chi 2 dihedrals were correctly positioned.(ABSTRACT TRUNCATED AT 400 WORDS)     

Energy functions for protein design: adjustment with protein-protein complex affinities, models for the unfolded state, and negative design of solubility and specificity

The development of the EGAD program and energy function for protein design is described. In contrast to most protein design methods, which require several empirical parameters or heuristics such as patterning of residues or rotamers, EGAD has a minimalist philosophy; it uses very few empirical factors to account for inaccuracies resulting from the use of fixed backbones and discrete rotamers in protein design calculations, and describes the unfolded state, aggregates, and alternative conformers explicitly with physical models instead of fitted parameters. This approach unveils important issues in protein design that are often camouflaged by heuristic-emphasizing methods. Inter-atom energies are modeled with the OPLS-AA all-atom forcefield, electrostatics with the generalized Born continuum model, and the hydrophobic effect with a solvent-accessible surface area-dependent term. Experimental characterization of proteins designed with an unmodified version of the energy function revealed problems with under-packing, stability, aggregation, and structural specificity. Under-packing was addressed by modifying the van der Waals function. By optimizing only three parameters, the effects of >400 mutations on protein-protein complex formation were predicted to within 1.0 kcal mol(-1). As an independent test, this modified energy function was used to predict the stabilities of >1500 mutants to within 1.0 kcal mol(-1); this required a physical model of the unfolded state that includes more interactions than traditional tripeptide-based models. Solubility and structural specificity were addressed with simple physical approximations of aggregation and conformational equilibria. The complete energy function can design protein sequences that have high levels of identity with their natural counterparts, and have predicted structural properties more consistent with soluble and uniquely folded proteins than the initial designs.     

Protein docking using continuum electrostatics and geometric fit

The computer program DOT quickly finds low-energy docked structures for two proteins by performing a systematic search over six degrees of freedom. A novel feature of DOT is its energy function, which is the sum of both a Poisson-Boltzmann electrostatic energy and a van der Waals energy, each represented as a grid-based correlation function. DOT evaluates the energy of interaction for many orientations of the moving molecule and maintains separate lists scored by either the electrostatic energy, the van der Waals energy or the composite sum of both. The free energy is obtained by summing the Boltzmann factor over all rotations at each grid point. Three important findings are presented. First, for a wide variety of protein-protein interactions, the composite-energy function is shown to produce larger clusters of correct answers than found by scoring with either van der Waals energy (geometric fit) or electrostatic energy alone. Second, free-energy clusters are demonstrated to be indicators of binding sites. Third, the contributions of electrostatic and attractive van der Waals energies to the total energy term appropriately reflect the nature of the various types of protein-protein interactions studied.     

Computing van der Waals energies in the context of the rotamer approximation

The rotamer approximation states that protein side-chain conformations can be described well using a finite set of rotational isomers. This approximation is often applied in the context of computational protein design and structure prediction to reduce the complexity of structural sampling. It is an effective way of reducing the structure space to the most relevant conformations. However, the appropriateness of rotamers for sampling structure space does not imply that a rotamer-based energy landscape preserves any of the properties of the true continuous energy landscape. Specifically, because the energy of a van der Waals interaction can be very sensitive to small changes in atomic separation, meaningful van der Waals energies are particularly difficult to calculate from rotamer-based structures. This presents a problem for computational protein design, where the total energy of a given structure is often represented as a sum of precalculated rigid rotamer self and pair contributions. A common way of addressing this issue is to modify the van der Waals function to reduce its sensitivity to atomic position, but excessive modification may result in a strongly nonphysical potential. Although many different van der Waals modifications have been used in protein design, little is known about which performs best, and why. In this paper, we study 10 ways of computing van der Waals energies under the rotamer approximation, representing four general classes, and compare their performance using a variety of metrics relevant to protein design and native-sequence repacking calculations. Scaling van der Waals radii by anywhere from 85 to 95% gives the best performance. Linearizing and capping the repulsive portion of the potential can give additional improvement, which comes primarily from getting rid of unrealistically large clash energies. On the other hand, continuously minimizing individual rotamer pairs prior to evaluating their interaction works acceptably in native-sequence repacking, but fails in protein design. Additionally, we show that the problem of predicting relevant van der Waals energies from rotamer-based structures is strongly nonpairwise decomposable and hence further modifications of the potential are unlikely to give significant improvement.     

The impact of <Emphasis Type="Italic">Trichoderma reesei</Emphasis> Cel7A carbohydrate binding domain mutations on its binding to a cellulose surface: a molecular dynamics free energy study

A critical role of the Family 7 cellobiohydrolase (Cel7A) carbohydrate binding domain (CBD) is to bind to a cellulose surface and increase the enzyme concentration on the surface. Several residues of Trichoderma reesei Cel7A CBD, including Y5, N29, Y31, Y32 and Q34, contribute to cellulose binding, as revealed by early experimental studies. To investigate the interactions between these important residues and cellulose, we applied a thermodynamic integration method to calculate the cellulose–Cel7A CBD binding free energy changes caused by Y5A, N29A, Y31A, Y32A and Q34A mutations. The experimental binding trend was successfully predicted, proving the effectiveness of the complex model. For the two polar residue mutants N29A and Q34A, the changes in the electrostatics are comparable to those of van der Waals, while for three Y to A mutants, the free energy differences mainly come from van der Waals interactions. However, in both cases, the electrostatics dominates the interactions between individual residues and cellulose. The side chains of these residues are rigidified after the complex is formed. The binding free energy changes for the two mutants Y5W and Y31W were also determined, and for these the van der Waals interaction was strengthened but the electrostatics was weakened.     

Solvation thermodynamics of biopolymers. I. Separation of the volume and surface interactions with estimates for proteins

The present paper is a systematic first approach to the problem of solvation thermodynamics of biomolecules. Most previous approaches have been only crude estimates of solvent contributions, and have simply assessed solvation free energy as proportional to surface areas. Here we estimate the various contributions and divide them into (a) hard-core interactions dependent upon the entire volume of solute and (b) the remainder of interactions manifested through surfaces, such as van der Waals, charge-charge, or hydrogen bonds. We have estimated the work to create a cavity with scaled-particle theory (SPT), the van der Waals interactions on the surface, and hydrogen bonds between the surface and the solvent. The conclusion here is that this latter term is the largest component of the solvation free energy of proteins. From estimates on nine diverse proteins, it is clear that the larger the protein, the more dominant is the hydrogen-bond term. In the next paper, we indicate that correlations between hydrogen-bonding groups on the surfaces could increase the magnitude of the hydrogen-bond contribution.     

Developing optimal non-linear scoring function for protein design

Motivation. Protein design aims to identify sequences compatible with a given protein fold but incompatible to any alternative folds. To select the correct sequences and to guide the search process, a design scoring function is critically important. Such a scoring function should be able to characterize the global fitness landscape of many proteins simultaneously. RESULTS: To find optimal design scoring functions, we introduce two geometric views and propose a formulation using a mixture of non-linear Gaussian kernel functions. We aim to solve a simplified protein sequence design problem. Our goal is to distinguish each native sequence for a major portion of representative protein structures from a large number of alternative decoy sequences, each a fragment from proteins of different folds. Our scoring function discriminates perfectly a set of 440 native proteins from 14 million sequence decoys. We show that no linear scoring function can succeed in this task. In a blind test of unrelated proteins, our scoring function misclassfies only 13 native proteins out of 194. This compares favorably with about three-four times more misclassifications when optimal linear functions reported in the literature are used. We also discuss how to develop protein folding scoring function.     

A method for determining overall protein fold from NMR distance restraints

Summary We describe a simple method for determining the overall fold of a polypeptide chain from NOE-derived distance restraints. The method uses a reduced representation consisting of two particles per residue, and a force field containing pseudo-bond and pseudo-angle terms, an electrostatic term, but no van der Waals or hard shell repulsive terms. The method is fast and robust, requiring relatively few distance restraints to approximate the correct fold, and the correct mirror image is readily determined. The method is easily implemented using commercially available molecular modeling software.     

Prediction of fluorine chemical shifts in proteins

Molecular dynamics calculations have been used in an effort to estimate the change in fluorine nmr shielding when a fluorine nucleus enters the tertiary structure of a protein. Considerations of the possible interactions that can define the shift parameter change suggest that van der Waals interactions are the leading determinant of fluorine shifts in proteins, although aromatic ring currents, other magnetic anisotropies, and electrostatic field effects could result in shift distinctions of 1 ppm or smaller. Results of our studies of a fluorine-containing analogue of the ribonuclease A S-protein/S-peptide complex indicate that static structures such as those implied by crystallographic data lead to overestimates of the magnitude of the van der Waals shielding term; molecular dynamics simulations provide indications of the effects of conformational averaging in defining this term. The treatment used predicts the correct direction of the shift change when the fluorine enters this protein environment from aqueous solution and, with an experimentally supported choice of adjustable parameters, gives agreement with the magnitude of the shift.     

An all atom energy based computational protocol for predicting binding affinities of protein-ligand complexes

We report here a computationally fast protocol for predicting binding affinities of non-metallo protein-ligand complexes. The protocol builds in an all atom energy based empirical scoring function comprising electrostatics, van der Waals, hydrophobicity and loss of conformational entropy of protein side chains upon ligand binding. The method is designed to ensure transferability across diverse systems and has been validated on a heterogenous dataset of 161 complexes consisting of 55 unique protein targets. The scoring function trained on a dataset of 61 complexes yielded a correlation of r=0.92 for the predicted binding free energies against the experimental binding affinities. Model validation and parameter analysis studies ensure the predictive ability of the scoring function. When tested on the remaining 100 protein-ligand complexes a correlation of r=0.92 was recovered. The high correlation obtained underscores the potential applicability of the methodology in drug design endeavors. The scoring function has been web enabled at as binding affinity prediction of protein-ligand (BAPPL) server.     

Van der Waals locks: loop-n-lock structure of globular proteins

In a globular protein the polypeptide chain returns to itself many times, making numerous chain-to-chain contacts. The stability of these contacts is maintained primarily by van der Waals interactions. In this work we isolated and analysed van der Waals contacts that stabilise spatial structures of nine major folds. We suggest a specific way to identify the tightest contacts of prime importance for the stability of a given crystallized protein and introduce the notion of the van der Waals lock. The loops closed by the van der Waals interactions provide a basically novel view of protein globule organization: the loop-n-lock structure. This opens a new perspective in understanding protein folding as well: the consecutive looping of the polypeptide chain and the locking of the loop ends by tight van der Waals interactions.     

Van der Waals forces. Special characteristics in lipid-water systems and a general method of calculation based on the Lifshitz theory

A practical method for examining and calculating van der Waals forces is derived from Lifshitz'' theory. Rather than treat the total van der Waals energy as a sum of pairwise interactions between atoms, the Lifshitz theory treats component materials as continua in which there are electromagnetic fluctuations at all frequencies over the entire body. It is necessary in principle to use total macroscopic dielectric data from component substances to analyze the permitted fluctuations; in practice it is possible to use only partial information to perform satisfactory calculations. The biologically interesting case of lipid-water systems is considered in detail for illustration. The method gives good agreement with measured van der Waals energy of interaction across a lipid film. It appears that fluctuations at infrared frequencies and microwave frequencies are very important although these are usually ignored in preference to UV contributions. “Retardation effects” are such as to damp out high frequency fluctuation contributions; if interaction specificity is due to UV spectra, this will be revealed only at interactions across <200 angstrom (A). Dependence of van der Waals forces on material electric properties is discussed in terms of illustrative numerical calculations.     

Monte Carlo refinement of rigid-body protein docking structures with backbone displacement and side-chain optimization

Structures of hitherto unknown protein complexes can be predicted by docking the solved protein monomers. Here, we present a method to refine initial docking estimates of protein complex structures by a Monte Carlo approach including rigid-body moves and side-chain optimization. The energy function used is comprised of van der Waals, Coulomb, and atomic contact energy terms. During the simulation, we gradually shift from a novel smoothed van der Waals potential, which prevents trapping in local energy minima, to the standard Lennard-Jones potential. Following the simulation, the conformations are clustered to obtain the final predictions. Using only the first 100 decoys generated by a fast Fourier transform (FFT)-based rigid-body docking method, our refinement procedure is able to generate near-native structures (interface RMSD <2.5 A) as first model in 14 of 59 cases in a benchmark set. In most cases, clear binding funnels around the native structure can be observed. The results show the potential of Monte Carlo refinement methods and emphasize their applicability for protein-protein docking.     

Calculation of binding affinities of HIV-1 RT and β-secretase inhibitors using the linear interaction energy method with explicit and continuum solvation approaches

The linear interaction energy (LIE) approach has been applied to estimate the binding free energies of representative sets of HIV-1 RT and β-Secretase inhibitors, using both molecular dynamics (MD) and tethered energy minimization sampling protocols with the OPLS-AA potential, using a range of solvation methodologies. Generalized Born (GB), ‘shell’ and periodic boundary condition (PBC) solvation were used, the latter with reaction field (RF) electrostatics. Poisson-Boltzmann (PB) and GB continuum electrostatics schemes were applied to the simulation trajectories for each solvation type to estimate the electrostatic ligand-water interaction energy in both the free and bound states. Reasonable agreement of the LIE predictions was obtained with respect to experimental binding free energy estimates for both systems: for instance, ‘PB’ fits on MD trajectories carried out with PBC solvation and RF electrostatics led to models with standard errors of 1.11 and 1.03 kcal mol⁻¹ and coefficients of determination, r ² of 0.76 and 0.75 for the HIV-1 RT and β-Secretase sets. However, it was also found that results from MD sampling using PBC solvation provided only slightly better fits than from simulations using shell or Born solvation or tethered energy minimization sampling. Figure Evolution of the running averages for compound H11 (binding to HIV-1RT) of the bound state ligand-water and ligand-protein interaction energies. The ligand-water electrostatic terms are twice the corresponding GB and PB electrostatic solvation free energies. The ligand-receptor van der Waals and Coulombic interaction energies are also shown, in addition to the ligand-water van der Waals interaction term. The terms were calculated (without application of a cut-off) from a trajectory sampled under PBC solvation with reaction field electrostatics Electronic supplementary material The online version of this article (doi:) contains supplementary material, which is available to authorized users.     

Free energy landscapes of two model peptides: α-helical and β-hairpin peptides explored with Brownian dynamics simulation

We applied an atomistic Brownian dynamics (BD) simulation with multiple time step method for the folding simulation of a 13-mer α-helical peptide and a 12-mer β-hairpin peptide, giving successful folding simulations. In this model, the driving energy contribution towards folding came from both electrostatic and van der Waals interactions for the α-helical peptide and from van der Waals interactions for the β-hairpin peptide. Although, many non-native structures having the same or lower energy than that of native structure were observed, the folded states formed the most populated cluster when the structures obtained by the BD simulations were subjected to the cluster analysis based on distance-based root mean square deviation of side-chains between different structures. This result indicates that we can predict the native structures from conformations sampled by BD simulation.     

Understanding beta-hairpin formation by molecular dynamics simulations of unfolding

We have studied the mechanism of formation of a 16-residue beta-hairpin from the protein GB1 using molecular dynamics simulations in an aqueous environment. The analysis of unfolding trajectories at high temperatures suggests a refolding pathway consisting of several transient intermediates. The changes in the interaction energies of residues are related with the structural changes during the unfolding of the hairpin. The electrostatic energies of the residues in the turn region are found to be responsible for the transition between the folded state and the hydrophobic core state. The van der Waals interaction energies of the residues in the hydrophobic core reflect the behavior of the radius of gyration of the core region. We have examined the opposing influences of the protein-protein (PP) energy, which favors the native state, and the protein-solvent (PS) energy, which favors unfolding, in the formation of the beta-hairpin structure. It is found that the behavior of the electrostatic components of PP and PS energies reflects the structural changes associated with the loss of backbone hydrogen bonding. Relative changes in the PP and PS van der Waals interactions are related with the disruption of the hydrophobic core of a protein. The results of the simulations support the hydrophobic collapse mechanism of beta-hairpin folding.     

Reintroducing electrostatics into macromolecular crystallographic refinement: application to neutron crystallography and DNA hydration

Most current crystallographic structure refinements augment the diffraction data with a priori information consisting of bond, angle, dihedral, planarity restraints, and atomic repulsion based on the Pauli exclusion principle. Yet, electrostatics and van der Waals attraction are physical forces that provide additional a priori information. Here, we assess the inclusion of electrostatics for the force field used for all-atom (including hydrogen) joint neutron/X-ray refinement. Two DNA and a protein crystal structure were refined against joint neutron/X-ray diffraction data sets using force fields without electrostatics or?with electrostatics. Hydrogen-bond orientation/geometry favors the inclusion of electrostatics. Refinement of Z-DNA with electrostatics leads to?a hypothesis for the entropic stabilization of Z-DNA that may partly explain the thermodynamics of converting the B form of DNA to its Z form. Thus, inclusion of electrostatics assists joint neutron/X-ray refinements, especially for placing and orienting hydrogen atoms.