Empirical free energy calculation: comparison to calorimetric data

An effective free energy potential, developed originally for binding free energy calculation, is compared to calorimetric data on protein unfolding, described by a linear combination of changes in polar and nonpolar surface areas. The potential consists of a molecular mechanics energy term calculated for a reference medium (vapor or nonpolar liquid), and empirical terms representing solvation and entropic effects. It is shown that, under suitable conditions, the free energy function agrees well with the calorimetric expression. An additional result of the comparison is an independent estimate of the side-chain entropy loss, which is shown to agree with a structure-based entropy scale. These findings confirm that simple functions can be used to estimate the free energy change in complex systems, and that a binding free energy evaluation model can describe the thermodynamics of protein unfolding correctly. Furthermore, it is shown that folding and binding leave the sum of solute-solute and solute-solvent van der Waals interactions nearly invariant and, due to this invariance, it may be advantageous to use a nonpolar liquid rather than vacuum as the reference medium.

Engineering subunit association of multisubunit proteins: a dimeric streptavidin

A dimeric streptavidin has been designed by molecular modeling using effective binding free energy calculations that decompose the binding free energy into electrostatic, desolvation, and side chain entropy loss terms. A histidine-127 --> aspartic acid (H127D) mutation was sufficient to introduce electrostatic repulsion between subunits that prevents the formation of the natural tetramer. However, the high hydrophobicity of the dimer-dimer interface, which would be exposed to solvent in a dimeric streptavidin, suggests that the resulting molecule would have very low solubility in aqueous media. In agreement with the calculations, a streptavidin containing the H127D mutation formed insoluble aggregates. Thus, the major design goal was to reduce the hydrophobicity of the dimer-dimer interface while maintaining the fundamental structure. Free energy calculations suggested that the hydrophobicity of the dimer-dimer interface could be reduced significantly by deleting a loop from G113 through W120 that should have no apparent contact with biotin in a dimeric molecule. The resulting protein, containing both the H127D mutation and the loop deletion, formed a soluble dimeric streptavidin in the presence of biotin.

Beta-turn new classification and its some features in proteins

An inspection of the phi-psi angle distribution strongly suggests that protein folding is highly constrained. A number of researchers have even suggested that a relatively small set of discrete phi-psi regions might be sufficient to describe most protein conformation. The total of 541 tight turns from 101 non-identical proteins were extracted form Brookhaven DataBank. The dihedral values of tight turns were scattered into the seven regions on the Ramachandran plot. These seven regions were called A1, A2, B1, B2, B22, T1 and T2. A1 and A2 are the traditional alpha-helix regions, B1, B2 and B22 the beta-strand regions, T1 and T2 the beta-turn regions. The A2 and T2 regions were not defined as "discrete" or single points but rather as one dimensional extended states. Based on the geometry of the two central residues of the tight turns, the new classification of beta-turn was defined. This classification of the majority of beta-turns fell into only six of the possible forty nine region combinations and were identifiable with the traditional nomenclature of Venkatachalam(1), but much simpler. The function of beta-turn in the conformation of proteins was studied. The hydrophobicity for different type turns was discussed. It shows that beta-turns have very strong hydrophilic property, so they are usually situated at the folding protein surface. The features of beta-turn and its amino acid distribution in this 541 beta-turn group and different type beta-turn were given.

Empirical potentials and functions for protein folding and binding

Simplified models and empirical potentials are being increasingly used for the analysis of proteins, frequently augmenting or replacing molecular mechanics approaches. Recent folding simulations have employed potentials that, in addition to terms assuring proper polypeptide geometry, include only two noncovalent effects-hydrogen bonding and hydrophobicity, with extremely simple approximations to the latter. The potentials that have been used in the free-energy ranking of protein-ligand complexes have generally been more involved. These potentials have more detailed solvation models and account for both local (hydrophobic and polar) solute-solvent phenomena and long range electrostatic solvation effects. The models of solvation that have been used most frequently are surface area related atomic parameters, knowledge-based models extracted from protein-structure data, and continum electrostatics with an additional area-related parameter. The knowledge-based approaches to solvation, although convenient and accurate enough, are suspect of double counting certain free-energy terms.

Predicted and trifluoroethanol-induced alpha-helicity of polypeptides

The alpha-helix stabilizing solvent 2,2,2-trifluoroethanol (TFE) is frequently used as a medium for determining the average alpha-helicity of polypeptides by CD spectroscopy. CD spectra measured in solutions containing 10, 15, 20, 50, and 90% (vol/vol) TFE are presented for 5 peptides that were selected to demonstrate possible variations in the effect of TFE concentration on the secondary structure. The analysis is extended to 6 further peptides whose CD spectra as measured in TFE are documented in the literature. The observed alpha-helicity at a high TFE concentration is compared with the alpha-helicity determined by a structure prediction method that combines conformational filtering [S. Vajda, (1993) Journal of Molecular Biology, Vol. 229, pp. 125-145], and a Monte Carlo simulation [J. Figge et al. (1993) Protein Science, Vol. 2, pp. 155-164]. For the set of 11 peptides we find a correlation of 0.84 between the predicted [theta]222 values and the corresponding values observed by CD spectroscopy in a high concentration of TFE (p < 0.01). Although we generally find a good correlation at high TFE concentration between observed and predicted alpha-helicity, there are several peptides that do not follow the predicted behavior. An analysis of the TFE titration curves in one such case revealed that TFE can induce a sharp transition from a partial beta-sheet conformation to an alpha-helical conformation as the TFE concentration is increased above a critical value.

Streptavidins with intersubunit crosslinks have enhanced stability

Natural tetrameric streptavidin has two subunit interfaces; one is a strong interface between subunits in a tightly associated dimer, and the other is a weak interface between a pair of such dimers (dimer-dimer interface). To test whether strengthening the weak dimer-dimer interface could provide streptavidin with additional structural stability, covalent crosslinks were introduced between adjacent subunits through the dimer-dimer interface. Specific crosslinking sites were designed by site-directed mutations of His-127 residues that are in close proximity in natural streptavidin. The first and second streptavidin constructs have a disulfide bond and an irreversible covalent bond, respectively, between two Cys-127 residues across the dimer-dimer interface. The third variant is a hybrid tetramer consisting of two different streptavidin species, one having lysine and the other aspartic acid at position 127, which are covalently crosslinked. All streptavidin constructs with intersubunit crosslinks showed higher biotin-binding ability than natural core streptavidin after heat treatment. All of these crosslinked streptavidins retained bound biotin more stably than natural core streptavidin in guanidine hydrochloride at very acidic pH. These results suggest that the introduction of covalent bonds across the dimer-dimer interface enhances the overall stability of streptavidin.

Free energy mapping of class I MHC molecules and structural determination of bound peptides

Free energy maps of the binding site are constructed for class I major histocompatibility complex (MHC) proteins, by rotating and translating amino acid probes along the cleft, and performing a side-chain conformational search at each position. The free energy maps are used to determine favorable residue positions that are then combined to form docked peptide conformations. Because the generic backbone structural motif of peptides bound to class I MHC is known, the mapping is restricted to appropriate regions of the site, but allows for the sometimes substantial variations in backbone and side-chain conformations. In a test demonstrating the quality of predictions for a known MHC site using only a rotational and conformational search, we started from the crystal structure of the HIV-1 gp120/HLA-A2 complex, and predicted the HLA-A2 bound structures of peptides from the influenza matrix protein, the HIV-1 reverse transcriptase, and the human T cell leukemia virus. The calculated peptides are at 1.6, 1.3, and 1.4 A all-atom RMSDs from their respective crystal structures (Madden DR, Garboczi DN, Wiley DC, 1993). A further test, which also included a local translational search, predicted structures across MHCs. In particular, we obtained the Kb/SEV-9 complex (Fremont DH et al., 1992, Science 257:919-927) starting with the complex between HLA-B27 and a generic peptide (Madden DR, Gorga JC, Strominger JL, Wiley DC, 1991, Nature (Lond) 353:321-325), with an all-atom RMSD of 1.2 A, indicating that the docking procedure is essentially as effective for predictions across MHCs as it is for determinations within the same MHC, although at substantially greater computational cost. The requirements for further improvement in accuracy are identified and discussed briefly.

Empirical free energy as a target function in docking and design: application to HIV-1 protease inhibitors

Structure-based drug design requires the development of efficient computer programs for exploring the structural compatibility of various flexible ligands with a given receptor. While various algorithms are available for finding docked conformations, selecting a target function that can reliably score the conformations remains a serious problem. We show that the use of an empirical free energy evaluation method, originally developed to characterize protein-protein interactions, can substantially improve the efficacy of search algorithms. In addition to the molecular mechanics interaction energy, the function takes account of solvation and side chain conformational entropy, while remaining simple enough to replace the incomplete target functions used in many drug docking and design procedures. The free energy function is used here in conjunction with a simple site mapping-fragment assembly algorithm, for docking the MVT-101 non-peptide inhibitor to HIV-1 protease. In particular, we predict the bound structure with an all atom RMSD of 1.21 A, compared to 1.69 A using an energy target function, and also accurately predict the free energy shifts obtained with a series of five trimeric hydroxyethylene isostere analogs.

Prediction of protein complexes using empirical free energy functions

A long sought goal in the physical chemistry of macromolecular structure, and one directly relevant to understanding the molecular basis of biological recognition, is predicting the geometry of bimolecular complexes from the geometries of their free monomers. Even when the monomers remain relatively unchanged by complex formation, prediction has been difficult because the free energies of alternative conformations of the complex have been difficult to evaluate quickly and accurately. This has forced the use of incomplete target functions, which typically do no better than to provide tens of possible complexes with no way of choosing between them. Here we present a general framework for empirical free energy evaluation and report calculations, based on a relatively complete and easily executable free energy function, that indicate that the structures of complexes can be predicted accurately from the structures of monomers, including close sequence homologues. The calculations also suggest that the binding free energies themselves may be predicted with reasonable accuracy. The method is compared to an alternative formulation that has also been applied recently to the same data set. Both approaches promise to open new opportunities in macromolecular design and specificity modification.

Extracting hydrophobicity parameters from solute partition and protein mutation/unfolding experiments

Hydrophobicity values for amino acids obtained from protein unfolding experiments are about twice as large as those obtained from data on the partitioning of amino acids between water and octanol. Quantitative analyses of several data sets, presented here, indicate that the difference is best explained by the most direct hypothesis, i.e. that the environment of hydrophobic groups in the interior of a protein is poorly modeled by octanol. Instead, we propose--and provide supporting evidence--that hydrocarbons are a more suitable model. First, we reanalyze data from both solute partitioning and protein unfolding experiments, taking account of the effects that were omitted previously, by introducing a volume dependence in the former and a full free energy analysis in the latter. Both changes in evaluation methodology decrease the discrepancy, but the differences remain substantial. The hydrophobicity parameter obtained from side-chain transfers between octanol and water increases from 16.7 to 22 cal/mol/Angstrom2, while that obtained from protein unfolding decreases from 34.9 to 31.2 cal/mol/Angstrom2. On the other hand, our analysis of the solubilities of pure hydrocarbons in water provides a hydrophobicity parameter of 30.8 cal/mol/Angstrom2. This apparent hydrocarbon-like environment of a protein's interior is also suggested more directly by an analysis of the contact environment of hydrophobic side chains in mutation/unfolding experiments, which have polar contact areas that are <2% of the total.

Flexible docking and design

Docking and design are the major computational steps toward understanding and affecting receptor-ligand interactions. The flexibility of many ligands makes these calculations difficult and requires the development and use of special methods. The need for such tools is illustrated by two examples: the design of protease inhibitors and the analysis and design of peptide antigens binding to specific MHC receptors. We review the computational concepts that have been extended from rigid-body to flexible docking, as well as the following important strategies for flexible docking and design: (a) Monte Carlo/molecular dynamics docking, (b) in-site combinatorial search, (c) ligand build-up, and (d) site mapping and fragment assembly. The use of empirical free energy as a target function is discussed. Due to the rapid development of the methodology, most new methods have been tested on only a limited number of applications and are likely to improve results obtained by more traditional computational or graphic tools.

Flexible docking of peptides to class I major-histocompatibility-complex receptors

We present a new method for docking flexible peptides to class I Major-Histocompatibility-Complex (MHC) receptors. Docking is performed in two steps: (a) The charged terminal peptide residues are located by randomly distributing multiple copies of each in volumes of approximately 150 A at either end of the binding groove, and then minimizing the system energy using a modified multiple-copy search algorithm. This is followed by (b) construction of the intervening chain using the multiple-copy bond-scaling-relaxation loop closure algorithm. In both steps, the copies tend to cluster and the size of the resulting clusters is proportional to the basin of attraction of the corresponding energy well. We show that native MHC-bound peptides have broad minima and, consequently, that misfolded, low-energy peptide conformations can be eliminated by restricting consideration to groups of peptides which cluster into broad minima. The accuracy of the method is assessed by comparing the predictions with crystallographic data for three different MHC peptide systems, at various degrees of stringency: (a) the extent to which we can determine side chain function (anchor vs. T-cell epitopes); (b) the extent to which we can determine the peptide-receptor orientation; and (c) the accuracy with which we can predict atomic coordinates. We find the method correct on (a) for 19 of the 22 non-Gly positions; the failures appearing to be a consequence of omitting solvation. Predictions related to (b) are also very encouraging, with the overall orientation of the predicted peptides being very similar to the crystal conformation, when measured by the hydrogen bonding pattern between the two. The degree of success in predicting atomic coordinates varied considerably, however, from 1.4 A for the HLA-A2 peptide to 2.7 A for the Kb peptide. The inaccuracy of the latter appears to reflect an incomplete target function, most likely the ommission of solvation. The calculations thus define the current limits of accuracy in docking flexible peptides to Class I receptors and identify the methodological improvements that must be made for the next advance in accuracy.

Effect of conformational flexibility and solvation on receptor-ligand binding free energies

A coherent framework is presented for determining the free energy change accompanying ligand binding to protein receptors. The most important new feature of the method is the contribution of the flexibility of the free ligand, and hence its conformational change on binding, to the free energy. Flexibility introduces two additional terms in the free energy difference: the internal energy difference between the ligand in the bound and free states and the backbone entropy loss. The former requires taking explicit account of the difference in solvation of the various forms of the free ligand. The solvation free energy change is estimated using an atomic solvation parameter model [Eisenberg & Mclachlan (1986) Nature 319, 199-203], with an improved parameter set. In order to evaluate the method, we applied it to three data sets for which increasingly general methods are required. The set to which the most restrictive theory can be applied consists of eight crystallized endopeptidase--protein inhibitor complexes which do not change conformation on binding and for which the major contribution to the solvation free energy is entropic. The results are in good agreement with the measured values and somewhat better than those previously reported in the literature. The second data set compares the relative binding free energies of biotin and its analogs for streptavidin. In this case the structures are also rigid, but solvation free energy must include both enthalpic and entropic components. We find that differential free energy predictions are approximately the same as those obtained by free energy perturbation techniques. The final application is an analysis of the measured stabilities of 13 different MHC receptor-peptide complexes. In this case we show that flexibility contributes 30-50% of the free energy change and find a correlation of 0.88 between our predicted free energies and peptide dissociation times.

Identifiability and indistinguishability of nonlinear pharmacokinetic models

Three nonlinear model structures of interest in pharmacokinetics are analyzed to determine whether the unknown, independent, model parameters can be deduced if perfect input-output data were available. This is the problem of identifiability. The method used is based on the local state isomorphism theorem. In certain circumstances, the modeler may be undecided between several model structures and it is then of interest to determine whether different model structures can be distinguished from perfect input-output data. This is the problem of model indistinguishability. The technique used, again based on the local state isomorphism theorem, parallels the similarity transformation approach for linear systems described previously in this journal. The analysis is performed on three two-compartment examples having one linear and one nonlinear (Michaelis-Menten) elimination pathway. In each model there is, on physiological and other grounds, some uncertainty over the precise location (central compartment or peripheral compartment) of one of the elimination pathways.

Indistinguishability for a class of nonlinear compartmental models

Indistinguishability, as applied to nonlinear compartmental models, is analyzed by means of the local state isomorphism theorem. The method of analysis involves the determination of all local, diffeomorphic transformations connecting the state variables of two models. This is then applied to two two-compartment models, in the first instance with linear eliminations, and then with the addition of eliminations with Michaelis-Menten kinetics. In the nonlinear example, the state transformation turns out to be linear or possibly affine. It is found that the nonlinear analysis could be eased by splitting the state isomorphism equations into those of the initial linear models together with extra equations due to the nonlinearities.

Computing the structure of bound peptides. Application to antigen recognition by class I major histocompatibility complex receptors

The ability to accurately compute the atomic positions of substrate-bound ligands is central to understanding biological recognition. Although substantial progress has been made in docking small, relatively rigid ligands, the problem of docking flexible peptides remains open. In this communication we present a new method that allows configurational flexibility of peptides, and apply it to predict the conformation of peptides bound to two class-I major histocompatibility complex receptors: human HLA-A2, and murine H-2Kb. Using only the approximate locations of the amino and carboxyl-terminal residues of the bound peptide, our calculations yield structures with backbone conformations that are similar to structures reported crystallographically.

Toward computational determination of peptide-receptor structure

We introduce a method for docking small flexible ligands of the size of dipeptides and phosphocholine and test it against crystallographic complexes. We then show how the method can be used as the basis for a strategy for solving the much more difficult problem of docking fully flexible peptides in the 8-10-residue size range. After developing the method we apply it to peptide-MHC class I systems and find that the predictions are in accord with biological and crystallographic data.

Determining protein loop conformation using scaling-relaxation techniques

We recently developed a rapid loop closure algorithm in which bond lengths are scaled to constrain the ends of a segment to match a known distance and then gradually relaxed to their standard values, with boundary constraints maintained. Although the algorithm predicted the Zif286 zinc-finger loop to within approximately 2 A, it had a serious limitation that made its more general use tentative: it omitted the atomic environment of the loop. Here we report an extension of the algorithm to take into account the protein environment surrounding a given loop from the outset of the conformational search and show that it predicts structure with an efficiency and accuracy that could not be achieved without continuous environmental inclusion. The algorithm should be widely applicable to structure determination when complete experimental information is unavailable.

The binding domain structure of retinoblastoma-binding proteins

The retinoblastoma gene product (Rb), a cellular growth suppressor, complexes with viral and cellular proteins that contain a specific binding domain incorporating three invariant residues: Leu-X-Cys-X-Glu, where X denotes a nonconserved residue. Hydrophobic and electrostatic properties are strongly conserved in this segment even though the nonconserved amino acids vary considerably from one Rb-binding protein to another. In this report, we present a diagnostic computer pattern for a high-affinity Rb-binding domain featuring the three conserved residues as well as the conserved physico-chemical properties. Although the pattern encompasses only 10 residues (with only 4 of these explicitly defined), it exhibits 100% sensitivity and 99.95% specificity in database searches. This implies that a certain pattern of structural and physico-chemical properties encoded by this short sequence is sufficient to govern specific Rb binding. We also present evidence that the secondary structural conformation through this region is important for effective Rb binding.

Conformational filtering in polypeptides and proteins

We present a method for assigning an ensemble of conformational states to each amino acid residue of a sequence. The states are defined as regions in the (phi, psi) map. The procedure is based on the use of conformational filters. In each filter we use a different set of approximations to estimate the probability of conformational states, and retain only the ones whose probability exceeds an acceptance probability. The resulting state assignment is not necessarily unique, but provides information that can be further exploited in searches for the tertiary structure. This conformational filtering approach to the de novo analysis of a sequence has a number of advantages over traditional structure prediction. First, it is possible to select acceptance probabilities such that the true conformational state is retained for up to 87% of residues, while substantially reducing the number of potential conformations. Second, in solution most linear peptides are present as ensembles of rapidly interconverting conformers, and such ensembles can be well predicted by filtering. Third, we can use Markov chains instead of a statistical mechanical (Ising) treatment, and avoid the need for estimating statistical weight matrices valid for the molecule as a whole. Markov models can use local transition matrices that are assumed to be independent of the rest of the chain, and are directly calculated from pairwise data. We show here that the locally identifiable transition matrices are transferable from the crystal structures of proteins to the solution structures of short peptides, and the ensembles of filtered conformations are in good agreement with nuclear magnetic resonance data. When applied to proteins, the filters retain several conformational states for most residues, and provide a measure of conformational variability. Small variability means that the segment is well defined by local interactions alone, and hence is likely to preserve its structure when isolated from the rest of the chain. Conversely, the structure of a segment with above-average conformational variability is likely to be significantly affected by its protein environment.

Necessary conditions for avoiding incorrect polypeptide folds in conformational search by energy minimization

Low energy conformations have been generated for melittin, pancreatic polypeptide, and ribonuclease S-peptide, both in the vicinity of x-ray structures by energy refinement and by an unconstrained search over the entire conformational space. Since the structural polymorphism of these medium-sized peptides in crystal and solution is moderate, comparing the calculated conformation to x-ray and nmr data provides information on local and global behavior of potential functions. Local analysis includes standardization calculations, which show that models with standard geometry can approximate good resolution x-ray data with less than 0.5 A rms deviation (RMSD). However, the atomic coordinates are shifted up to 2 A RMSD by local energy minimization, and thus 2 A is generally the smallest RMSD value one can target in a conformational search using the same energy evaluation models. The unconstrained search was performed by a buildup-type method based on dynamic programming. To accelerate the generation of structures in the conformational search, we used the ECEPP potential, defined in terms of standard polypeptide geometry. A number of low energy conformations were further refined by relaxing the assumption of standard bond lengths and bond angles through the use of the CHARMM potential, and the hydrophobic folding energies of Eisenberg and McLachlan were calculated. Each conformation is described in terms of the RMSD from the native, hydrogen-bonding structure, solvent-accessible surface area, and the ratio of surfaces corresponding to nonpolar and polar residues. The unconstrained search finds conformations that are different from the native, sometimes substantially, and in addition, have lower conformational energies than the native. The origin of deviations is different for each of the three peptides, but in all examples the refined x-ray structures have lower energies than the calculated incorrect folds when (1) the assumption of standard bond lengths and bond angles is relaxed; (2) a small and constant effective dielectric permittivity (epsilon < 10) is used; and (3) the hydrophobic folding energy is incorporated into the potential.

A membrane model for cytosolic calcium oscillations. A study using Xenopus oocytes

Cytosolic calcium oscillations occur in a wide variety of cells and are involved in different cellular functions. We describe these calcium oscillations by a mathematical model based on the putative electrophysiological properties of the endoplasmic reticulum (ER) membrane. The salient features of our membrane model are calcium-dependent calcium channels and calcium pumps in the ER membrane, constant entry of calcium into the cytosol, calcium dependent removal from the cytosol, and buffering by cytoplasmic calcium binding proteins. Numerical integration of the model allows us to study the fluctuations in the cytosolic calcium concentration, the ER membrane potential, and the concentration of free calcium binding sites on a calcium binding protein. The model demonstrates the physiological features necessary for calcium oscillations and suggests that the level of calcium flux into the cytosol controls the frequency and amplitude of oscillations. The model also suggests that the level of buffering affects the frequency and amplitude of the oscillations. The model is supported by experiments indirectly measuring cytosolic calcium by calcium-induced chloride currents in Xenopus oocytes as well as cytosolic calcium oscillations observed in other preparations.

Nanosecond fluorescence of tryptophans in cytochrome P-450scc (CYP11A1): effect of substrate binding

Fluorescence of eight tryptophan residues in cytochrome P-450scc with bound endogenous cholesterol could be fitted with a two component model: a single exponential and a "top-hat" distribution of lifetimes as the second component. The short-lived component (tau 1 about 700 ps) does not change significantly upon binding of substrate (22R-hydroxycholesterol). The parameters of the long-lived component (central lifetime tau m about 3.4 ns) change upon binding of carbon monoxide and substrate. 22R-hydroxycholesterol binding broadens the distribution of the long-lived component; that is the heterogeneity of the Trp environment is increased when this substrate displaces the endogenous cholesterol.