Computers in molecular biology: current applications and emerging trends.

The rate of generation of molecular sequence data is forcing the use of computers as a central tool in molecular biology. Current use of computers is limited largely to data management and sequence comparisons, but rapid growth in the volume of data is generating pressure for the development of high-speed analytical methods for deciphering the codes connecting nucleotide sequence with protein structure and function.

Analysis of breast cancer progression using principal component analysis and clustering.

We develop a new technique to analyse microarray data which uses a combination of principal components analysis and consensus ensemble k-clustering to find robust clusters and gene markers in the data. We apply our method to a public microarray breast cancer dataset which has expression levels of genes in normal samples as well as in three pathological stages of disease; namely, atypical ductal hyperplasia or ADH, ductal carcinoma in situ or DCIS and invasive ductal carcinoma or IDC. Our method averages over clustering techniques and data perturbation to find stable, robust clusters and gene markers. We identify the clusters and their pathways with distinct subtypes of breast cancer (Luminal,Basal and Her2+). We confirm that the cancer phenotype develops early (in early hyperplasia or ADH stage) and find from our analysis that each subtype progresses from ADH to DCIS to IDC along its own specific pathway, as if each was a distinct disease.

Method for quantitating tumor cell removal and tumor cell-invasive capacity in experimental metastases.

An improved method is presented for analyzing the decay of i.v.-injected labeled tumor cells in the lung. A simple compartment analysis yields quantitative values for (a) the rate at which tumor cells are lost or cleared from the lung, and (b) the "invasion" rate at which tumor cells colonize the interstitial space. Statistical methods are outlined for testing the significance of the contribution of each rate to the overall shape of the decay curve and for testing whether or not a given rate is altered significantly during an experiment. The usefulness of the method is demonstrated in analysis of experiments involving perturbations of the host or of the injected tumor cells.

Prediction of protein function from sequence properties. Discriminant analysis of a data base.

The protein superfamilies in the National Biomedical Research Foundation sequence data base cluster into six groups that can be distinguished on the basis of four variables characterizing amino acid composition and local sequence properties. The variables are average hydrophobicity, net charge, sequence length and periodic variation in hydrophobic residues along the chain. The clusters they distinguish are: globins; chromosomal proteins; contractile system proteins and respiratory proteins other than cytochromes; enzyme inhibitors and toxins; enzymes except hydrolases; and all other proteins. The overall probability of correctly allocating a given protein to one of these functional groups is 0.76, with the allocation reliability being highest for globins (0.97) and for chromosomal proteins (0.93).

The detection and classification of membrane-spanning proteins.

Discriminant analysis can be used to precisely classify membrane proteins as integral or peripheral and to estimate the odds that the classification is correct. Specifically, using 102 membrane proteins from the National Biomedical Research Foundation (NBRF) database we find that discrimination between integral and peripheral membrane proteins can be achieved with 99% reliability. Hydrophobic segments of integral membrane proteins can also be distinguished from interior segments of globular soluble proteins with better than 95% reliability. We also propose a procedure for determining boundaries of membrane-spanning segments and apply it to several integral membrane proteins. For the limited data available (such as on transplantation antigens), the residues at the boundaries of a membrane-spanning segment are predictable to within the error inherent in the concept of boundary. As a specific indication of resolution, seven membrane-spanning segments of bacteriorhodopsin are resolved with no information other than sequence, and the predicted boundary residues agree with the experimental data on proteolytic cleavage sites. Several definitive but yet to be tested predictions are also made, and the relation to other predictive methods is briefly discussed. A computer program in FORTRAN for prediction of membrane-spanning segments is available from the authors.

Ligand binding to multiple equivalent sites with steric hindrance.

A general model is developed for the binding of ligands to multiple receptor sites in which steric blockage of sites is taken into account. Analytical expressions for extent of ligand binding as a function of ligand concentration are derived employing a stochastic matrix approach. Using simple numerical techniques to evaluate these expressions it is possible to obtain the inherent affinity of each site for a ligand, the number of sites available and the number of sites excluded when a ligand binds to any site. The expressions derived here are contrasted with other expressions based on simple equilibrium considerations in which there are no interactions between sites and no interactions between ligands in different sites. It is shown that the expressions derived here predict significant departures from linearity in Scatchard plots and a strong negative cooperatively, especially toward the saturation limit.

Estimating the number of protein folds.

A number of fundamental questions in structural biology concern the diversity of protein architectures (or folds). Here, we address two of them, the size of the universe of folds, and the distribution of sequence families among them, using an analysis based on a new and rigorous statistical sampling method. In particular we show that the number of known non-transmembrane protein folds is approximately one half of the total that exist, and that certain superfolds should exist, which accommodate dozens of non-homologous sequence families.

Structural principles that govern the peptide-binding motifs of class I MHC molecules.

The peptides that bind class I MHC molecules are restricted in length and often contain key amino acids, anchor residues, at particular positions. The side-chains of peptide anchor residues interact with the polymorphic complementary pockets in MHC peptide-binding grooves and provide the molecular basis for allele-specific recognition of antigenic peptides. We establish correlations between class I MHC specificities for anchor residues and class I MHC sequence markers that occur at the polymorphic positions lining the structural pockets. By analyzing the pocket structures of nine crystallized class I MHC molecules and the modeled structures of another 39 class I MHC molecules, we show that class I pockets can be classified into families that are distinguishable by their common physico-chemical properties and peptide side-chain selectivities. The identification of recurrent structural principles among class I pockets makes it possible to greatly expand the repertoire of known peptide-binding motifs of class I MHC molecules. The evolutionary strategies underlying the emergence of pocket families is briefly discussed.

Two complementary methods for predicting peptides binding major histocompatibility complex molecules.

Peptides that bind to major histocompatibility complex products (MHC) are known to exhibit certain sequence motifs which, though common, are neither necessary nor sufficient for binding: MHCs bind certain peptides that do not have the characteristic motifs and only about 30% of the peptides having the required motif, bind. In order to develop and test more accurate methods we measured the binding affinity of 463 nonamer peptides to HLA-A2.1. We describe two methods for predicting whether a given peptide will bind to an MHC and apply them to these peptides. One method is based on simulating a neural network and another, called the polynomial method, is based on statistical parameter estimation assuming independent binding of the side-chains of residues. We compare these methods with each other and with standard motif-based methods. The two methods are complementary, and both are superior to sequence motifs. The neural net is superior to simple motif searches in eliminating false positives. Its behavior can be coarsely tuned to the strength of binding desired and it is extendable in a straightforward fashion to other alleles. The polynomial method, on the other hand, has high sensitivity and is a superior method for eliminating false negatives. We discuss the validity of the independent binding assumption in such predictions.

TcR recognition of the MHC-peptide dimer: structural properties of a ternary complex.

We have developed a method that utilizes site-specific mutation data, sequence analysis, immunological data and free-energy minimization, to determine structural features of the ternary complex formed by the T-cell receptor (TcR) and the class I major histocompatibility complex (MHC) molecule bound by peptide. The analysis focuses on the mouse Kd MHC system, for which a large set of clones with sequenced T-cell receptors is available for specific peptides. The general philosophy is to reduce the uncertainties and computation time in a free-energy minimization procedure by identifying and imposing experimental constraints. In addition to assessing compatibility with various kinds of immunological data, we are particularly interested in differentiating the structural features peculiar to this particular system from generic features, and in ascertaining the robustness of the structure; i.e. determining, in so far as possible, the variations in the structure that leave its compatibility with experiment unaltered from those that do not. This last is equivalent to recognizing that certain features of the model are presented with a reasonable degree of confidence, while others remain highly tentative. The central conclusion in the former category is a placement of the TcR on the Kd peptide complex, which has its beta 2, beta 3 and alpha 3 loops (i.e. the second and third complementarity-determining region of the TcR beta chain, and the third complementarity-determining region of the alpha chain) covering the peptide; the alpha 1 and alpha 2 loops covering the MHC alpha 1 helix; the alpha 2 loop interacting with residues on the MHC beta sheet; and the beta 1 and (part of) the beta 2 loops covering the alpha 2 MHC helix. More specifically, our findings include the following. (1) A highly conserved histidine residue in the first complementarity-determining region of the TcR beta chain (beta:CDR1) points outward and interacts with highly conserved side-chains on the MHC alpha 2 helix. (2) The amino-terminal portion of the beta 2 loop interacts with the carboxyl portion of the peptide. A particularly important interaction is K4 of the loop interacting with E8 of the peptide. (3) Charged side-chains of the 11-residue TcR alpha 2 loop interact with conserved charged side-chains at positions 44, 58, 61 and 68 on the MHC. (4) The TcR beta 3 loop interacts with the amino-terminal part of the peptide, up through position 4. (5) the TcR alpha 3 loop interacts with the central portion of the peptide and stacks against the beta 2 loop. (6) Because of the interaction between the beta 2 loop and the peptide, and stacking of beta 2 on alpha 3, alpha 3 gene and V beta gene selection can be correlated. (7) Using the topology of the recently solved TcR alpha chain we predict that the alpha 2 loop interacts with the loop on the MHC beta sheet floor, which encompasses residues 42 to 44.

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.

Determination of atomic desolvation energies from the structures of crystallized proteins.

We estimated effective atomic contact energies (ACE), the desolvation free energies required to transfer atoms from water to a protein's interior, using an adaptation of a method introduced by S. Miyazawa and R. L. Jernigan. The energies were obtained for 18 different atom types, which were resolved on the basis of the way their properties cluster in the 20 common amino acids. In addition to providing information on atoms at the highest resolution compatible with the amount and quality of data currently available, the method itself has several new features, including its reference state, the random crystal structure, which removes compositional bias, and a scaling factor that makes contact energies quantitatively comparable with experimentally measured energies. The high level of resolution, the explicit accounting of the local properties of protein interiors during determination of the energies, and the very high computational efficiency with which they can be assigned during any computation, should make the results presented here widely applicable. First we used ACE to calculate the free energies of transferring side-chains from protein interior into water. A comparison of the results thus obtained with the measured free energies of transferring side-chains from n-octanol to water, indicates that the magnitude of protein to water transfer free energies for hydrophobic side-chains is larger than that of n-octanol to water transfer free energies. The difference is consistent with observations made by D. Shortle and co-workers, who measured differential free energies of protein unfolding for site-specific mutants in which Ala or Gly was substituted for various hydrophobic side-chains. A direct comparison (calculated versus observed free energy differences) with those experiments finds slopes of 1.15 and 1.13 for Gly and Ala substitutions, respectively. Finally we compared calculated and observed binding free energies of nine protease-inhibitor complexes. This requires a full free energy function, which is created by adding direct electrostatic interactions and an appropriate entropic component to the solvation free energy term. The calculated free energies are typically within 10% of the observed values. Taken collectively, these results suggest that ACE should provide a reasonably accurate and rapidly evaluatable solvation component of free energy, and should thus make accessible a range of docking, design and protein folding calculations that would otherwise be difficult to perform.

Hydrophobicity scales and computational techniques for detecting amphipathic structures in proteins.

Protein segments that form amphipathic alpha-helices in their native state have periodic variation in the hydrophobicity values of the residues along the segment, with a 3.6 residue per cycle period characteristic of the alpha-helix. The assignment of hydrophobicity values to amino acids (hydrophobicity scale) affects the display of periodicity. Thirty-eight published hydrophobicity scales are compared for their ability to identify the characteristic period of alpha-helices, and an optimum scale for this purpose is computed using a new eigenvector method. Two of the published scales are also characterized by eigenvectors. We compare the usual method for detecting periodicity based on the discrete Fourier transform with a method based on a least-squares fit of a harmonic sequence to a sequence of hydrophobicity values. The two become equivalent for very long sequences, but, for shorter sequences with lengths commonly found in alpha-helices, the least-squares procedure gives a more reliable estimate of the period. The analog to the usual Fourier transform power spectrum is the "least-squares power spectrum", the sum of squares accounted for in fitting a sinusoid of given frequency to a sequence of hydrophobicity values. The sum of the spectra of the alpha-helices in our data base peaks at 97.5 degrees, and approximately 50% of the helices can account for this peak. Thus, approximately 50% of the alpha-helices appear to be amphipathic, and, of those that are, the dominant frequency at 97.5 degrees rather than 100 degrees indicates that the helix is slightly more open than previously thought, with the number of residues per turn closer to 3.7 than 3.6. The extra openness is examined in crystallographic data, and is shown to be associated with the C terminus of the helix. The alpha amphipathic index, the key quantity in our analysis, measures the fraction of the total spectral area that is under the 97.5 degrees peak, and is a characteristic of hydrophobicity scales that is consistent for different sets of helices. Our optimized scale maximizes the amphipathic index and has a correlation of 0.85 or higher with nine previously published scales. The most surprising feature of the optimized scale is that arginine tends to behave as if it were hydrophobic; i.e. in the crystallographic data base it has a tendency to be on the hydrophobic face of teh amphipathic helix. Although the scale is optimal only for predicting alpha-amphipathicity, it also ranks high in identifying beta-amphipathicity and in distinguishing interior from exterior residues in a protein.(ABSTRACT TRUNCATED AT 400 WORDS)

Computational determination of the structure of rat Fc bound to the neonatal Fc receptor.

The available crystal structure for the complex between the Fc fragment of immunoglobulin G (IgG) and the neonatal Fc receptor (FcRn) was determined at low resolution and has no electron density for a large portion of the CH2 domain of the Fc. Here, we use a well validated computational docking algorithm in conjunction with known crystallographic data to predict the orientation of CH2 when bound to FcRn, and validate the predicted structure with data from site-specific mutagenesis experiments. The predicted Fc structure indicates that the CH2 domain moves upon binding FcRn , such that the end-to-end distance of the bound Fc fragment is greater than it is in the crystal structure of isolated Fc. The calculated orientation of the bound CH2 domain is displaced by an average of 6 A from the CH2 orientation in the structure of Fc alone, and shows improved charge complementarity with FcRn. The predicted effects of 11 specific mutations in Fc and FcRn are calculated and the results are compared with experimental measurements. The predicted structure is consistent with all reported mutagenesis data, some of which are explicable only on the basis of our model. The current study predicts that FcRn-bound Fc is asymmetric due to reorientation of the CH2 domain upon FcRn binding, a rearrangement that would be likely to interfere with optimal binding of FcRn at the second binding site of the Fc homodimer.

A method for determining whether the descending limb of a biphasic histamine release curve reflects insufficient cross-linking.

We introduce a kinetic method for determining whether the descending limb of a biphasic histamine release dose-response curve is the result of insufficient cross-linking, and delineate conditions under which it is applicable. The method involves examining kinetic curves showing for various fixed antigen concentrations the cumulative amount of histamine release as a function of time. From the slope of the kinetic curves measured at some fixed time one determines how the rate of release depends on concentration. We show under very general conditions that if the dose-response curve for histamine release reaches a peak, and then decreases over a concentration interval in which the rate of release does not decline, then the decline in the dose-response curve cannot be due to insufficient cross-linking. Consequently, a characteristic feature of antigen-excess inhibition of histamine release due to mechanisms other than insufficient cross-linking is the crossing of a kinetic curve generated at a suboptimal antigen concentration by a kinetic curve generated at a supraoptimal antigen concentration. We show that the technique is easily executed experimentally and provide kinetic evidence-suggesting that the rabbit basophils the antigen-excess inhibition of histamine release by bis-benzylpenicilloyl-1,6-diaminohexane (BPO2) is due to insufficient cross-linking, whereas the antigen-excess inhibition observed with ovalbumin probably is due to more complete desensitization mechanisms.

Identifying weak signals in the presence of noise: a new method of locating potential ligand contact residues in immunoglobulin-related molecules.

We develop and apply a new method for estimating the locations of hypervariable residues in immunoglobulin-related molecules. The method differs from the standard introduced by Wu and Kabat in two essential ways: (1) we take explicit account of the type of substitution at a given position, rather than just the total number of substitutions, and (2) we use an explicit statistical decision criterion for classifying a site into either the complementarity determining or framework category. Simulations indicate that the method is reliable with relatively little data, approximately 5% of the sites being misclassified when 10 sequences are aligned. The method is applied to immunoglobulin light chains and to class 1 and class 2 products of the major histocompatibility complex.

Computational determination of side chain specificity for pockets in class I MHC molecules.

We show that a rapidly executable computational procedure provides the basis for a predictive understanding of antigenic peptide side chain specificity, for binding to class I major histocompatibility complex (MHC) molecules. The procedure consists of a combined search to identify the joint conformations of peptide side chains and side chains comprising the MHC pocket, followed by conformational selection, using a target function, based on solvation energies and modified electrostatic energies. The method was applied to the B pocket region of five MHC molecules, which were chosen to encompass the full range of specificities displayed by anchors at peptide position 2. These were a medium hydrophobic residue (Leu or Met) for HLA-A*0201, a basic residue (Arg or Lys) for HLA-B*2705; a small hydrophobic residue (Val) for HLA-A*6801, an acidic residue (Glu) for HLA-B*4001 and a bulky residue (Tyr) for H-2K(d). The observed anchors are correctly predicted in each case. The agreement for HLA-B40 and H-2K(d) is especially promising, since their structures have not yet been determined experimentally. Because the experimental determination of motifs by elution is difficult and these calculations take only hours on a high speed workstation, the results open the possibility of routine determination of motifs computationally.

Characterization of a helper T cell epitope recognized by mice of a low responder major histocompatibility type.

Most known helper T cell (Th) epitopes studied have naturally been immunodominant epitopes recognized by T cells from animals of high responder major histocompatibility complex (MHC) haplotype. We have previously found that most such immunodominant Th epitopes tend to be amphipathic alpha helices, that is, helices with hydrophobic residues on one side and hydrophilic residues on the other, and the corresponding peptide can usually elicit a response to the native protein. However, very few epitopes seen by MHC low responder T cells have been identified. Within the CNBr fragment of residues 1-55 of sperm whale myoglobin (SwMb), a Th epitope is known to exist that stimulates T cells from low responder H-2k mice, but it has not yet been localized to a length of 8-12 residues, the usual length of a Th epitope. To determine whether this low responder epitope would have similar properties, we located it using 10 evenly overlapping 15-residue peptides that span the region. Analysis of this region by the computer program predicted the site covered by two peptides (residues 26-40 and 31-45 which overlap by 10 residues) to be the most likely site for a Th epitope. Of the 10 peptides tested experimentally, only one peptide (residues 26-40) was able to stimulate two low responder Th clones that are specific for the 1-55 region. The peptide was able to prime T cells of low responder B10.BR mice in vivo for in vitro response to the native SwMb as well as to the peptide fragment of residues 1-55. Immunization of low responder mice with SwMb showed that, of the 10 overlapping peptides, the major site of response within the 1-55 region is to the identified peptide. Finally, an extended peptide of residues 24-42 was made to increase the amphipathic score. This extended peptide induced greater proliferation of the clones. Thus, this low responder epitope has properties similar to those of immunodominant epitopes recognized by high responders.

Minimal requirements for peptide mediated activation of CD8+ CTL.

A physical chemical model of T cell stimulation by class I-peptide complexes was developed and used to analyse in vitro studies of gamma-interferon release as a function of the number of peptide and MHC molecules. The analysis provided reasonable estimates of well identified parameters, including equilibrium constants and the minimum number of T cell receptor-class I-peptide ternary complexes on a presenting cell required to activate T cells. The latter number was estimated as 3-5 per T cell. This is in distinct contrast to estimates in the literature of the number of peptide-MHC complexes required for activity, which is necessarily larger. The analysis also predicted that activity is potentiated by interaction between class I molecules, even if one member of the pair is not bound by antigen. The analytical approach used in this paper may be applicable to other activation systems.

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.