Continuum and atomistic modeling of ion partitioning into a peptide nanotube.

Continuum and atomistic descriptions of the partitioning of ions into a self-assembled (D,L)-octapeptide nanotube, cyclo[-(L-Ala-D-Ala)(4)-], are presented. Perturbation free energy calculations, including Ewald electrostatics, are used to estimate the electrostatic component of the excess free energy of charging Li(+), Na(+), Rb(+), and Cl(minus sign) ions inside the nanotube. The radial density and orientational distribution of water around the ion is calculated for the ion at two different positions inside the tube; it is seen that the calculated distributions are sensitive to the location of the ions. Two different continuum electrostatic models are formulated to describe the ion solvation inside the nanotube. When enhanced orientational structuring of water dipoles is evidenced, explicitly including the first solvation shell as part of the low dielectric nanotube environment provides good agreement with molecular dynamics simulations. When water orientational structuring is as in the reference bulk solvent, we find that treating the first shell water explicitly or as a high dielectric continuum leads to similar results. These results are discussed, and their importance for continuum electrostatic modeling of ion channels are highlighted.

pK(a) Calculations suggest storage of an excess proton in a hydrogen-bonded water network in bacteriorhodopsin.

Calculations of protonation states and pK(a) values for the ionizable groups in the resting state of bacteriorhodopsin have been carried out using the recently available 1.55 A resolution X-ray crystallographic structure. The calculations are in reasonable agreement with the available experimental data for groups on or near the ion transport chain (the retinal Schiff base; Asp85, 96, 115, 212, and Arg82). In contrast to earlier studies using lower-resolution structural data, this agreement is achieved without manipulations of the crystallographically determined heavy-atom positions or ad hoc adjustments of the intrinsic pK(a) of the Schiff base. Thus, the theoretical methods used provide increased reliability as the input structural data are improved. Only minor effects on the agreement with experiment are found with respect to methodological variations, such as single versus multi-conformational treatment of hydrogen atom placements, or retaining the crystallographically determined internal water molecules versus treating them as high-dielectric cavities. The long-standing question of the identity of the group that releases a proton to the extracellular side of the membrane during the L-to-M transition of the photocycle is addressed by including as pH-titratable sites not only Glu204 and Glu194, residues near the extracellular side that have been proposed as the release group, but also an H(5)O(2)(+) molecule in a nearby cavity. The latter represents the recently proposed storage of the release proton in an hydrogen-bonded water network. In all calculations where this possibility is included, the proton is stored in the H(5)O(2)(+) rather than on either of the glutamic acids, thus establishing the plausibility on theoretical grounds of the storage of the release proton in bacteriorhodopsin in a hydrogen-bonded water network. The methods used here may also be applicable to other proteins that may store a proton in this way, such as the photosynthetic reaction center and cytochrome c oxidase.

Solvation energy density occlusion approximation for evaluation of desolvation penalties in biomolecular interactions.

In calculations involving many displacements of an interacting pair of biomolecules, such as brownian dynamics, the docking of a substrate/ligand to an enzyme/receptor, or the screening of a large number of ligands as prospective inhibitors for a particular receptor site, there is a need for rapid evaluation of the desolvation penalties of the interacting pair. Although continuum electrostatic treatments with distinct dielectric constants for solute and solvent provide an account of the electrostatics of solvation and desolvation, it is necessary to re-solve the Poisson equation, at considerable computational cost, for each displacement of the interacting pair. We present a new method that uses a formulation of continuum electrostatic solvation in terms of the solvation energy density and approximates desolvation in terms of the occlusion of this density. We call it the SEDO approximation. It avoids the need to re-solve the Poisson equation, as desolvation is now estimated by an integral over the occluded volume. Test calculations are presented for some simple model systems and for some real systems that have previously been studied using the Poisson equation approach: MHC class I protein-peptide complexes and a congeneric series of human immunodeficiency virus type 1 (HIV-1) protease--ligand complexes. For most of the systems considered, the trends and magnitudes of the desolvation component of interaction energies obtained using the SEDO approximation are in reasonable correlation with those obtained by re-solving the Poisson equation. In most cases, the error introduced by the SEDO approximation is much less than that of the often-used test-charge approximation for the charge-charge components of intermolecular interactions. Proteins 2001;43:12-27.

Generalized born models of macromolecular solvation effects.

It would often be useful in computer simulations to use a simple description of solvation effects, instead of explicitly representing the individual solvent molecules. Continuum dielectric models often work well in describing the thermodynamic aspects of aqueous solvation, and approximations to such models that avoid the need to solve the Poisson equation are attractive because of their computational efficiency. Here we give an overview of one such approximation, the generalized Born model, which is simple and fast enough to be used for molecular dynamics simulations of proteins and nucleic acids. We discuss its strengths and weaknesses, both for its fidelity to the underlying continuum model and for its ability to replace explicit consideration of solvent molecules in macromolecular simulations. We focus particularly on versions of the generalized Born model that have a pair-wise analytical form, and therefore fit most naturally into conventional molecular mechanics calculations.

Protonation states and pH titration in the photocycle of photoactive yellow protein.

Photoactive yellow protein (PYP) undergoes a light-driven cycle of color and protonation states that is part of a mechanism of bacterial phototaxis. This article concerns functionally important protonation states of PYP and the interactions that stabilize them, and changes in the protonation state during the photocycle. In particular, the chromophore pK(a) is known to be shifted down so that the chromophore is negatively charged in the ground state (dark state) even though it is buried in the protein, while nearby Glu46 has an unusually high pK(a). The photocycle involves changes of one or both of these protonation states. Calculations of pK(a) values and protonation states using a semi-macroscopic electrostatic model are presented for the wild-type and three mutants, in both the ground state and the bleached (I(2)) intermediate state. Calculations allowing multiple H-bonding arrangements around the chromophore also have been carried out. In addition, ground-state pK(a) values of the chromophore have been measured by UV-visible spectroscopy for the wild-type and the same three mutants. Because of the unusual protonation states and strong electrostatic interactions, PYP represents a severe test of the ability of theoretical models to yield correct calculations of electrostatic interactions in proteins. Good agreement between experiment and theory can be obtained for the ground state provided the protein interior is assumed to have a relatively low dielectric constant, but only partial agreement between theory and experiment is obtained for the bleached state. We also present a reinterpretation of previously published data on the pH-dependence of the recovery of the ground state from the bleached state. The new analysis implies a pK(a) value of 6.37 for Glu46 in the bleached state, which is consistent with other available experimental data, including data that only became available after this analysis. The new analysis suggests that signal transduction is modulated by the titration properties of the bleached state, which are in turn determined by electrostatic interactions. Overall, the results of this study provide a quantitative picture of the interactions responsible for the unusual protonation states of the chromophore and Glu46, and of protonation changes upon bleaching.

Electrostatic coupling to pH-titrating sites as a source of cooperativity in protein-ligand binding.

This paper describes an alternative mechanism for the cooperative binding of charged ligands to proteins. The ligand-binding sites are electrostatically coupled to protein side chains that can undergo protonation and deprotonation. The binding of one ligand alters the protein's protonation equilibrium in a manner that makes the the binding of the second ligand more favorable. This mechanism requires no conformational change to produce a cooperative effect, although it is not exclusive of conformational change. We present a theoretical description of the mechanism, and calculations on three kinds of systems: A model system containing one protonation site and two ligand-binding sites; a model system containing two protonation sites and two ligand-binding sites; and calbindin D9k, which contains two Ca2+-binding sites and 30 protonation sites. For the one-protonation-site model, it is shown that the influence of the protonation site can only be cooperative. The competition of this effect with the anticooperative effect of ligand-ligand repulsion is studied in detail. For the two-protonation site model, the effect can be either cooperative or, in special cases, anticooperative. For calbindin D9k, the calculations predict that six protonation sites in or near the ligand-binding sites make a cooperative contribution that approximately cancels the anticooperative effect of Ca2+-Ca2+ repulsion, accounting for more than half of the total cooperative effect that is needed to overcome repulsion and produce the net cooperativity observed experimentally. We argue that cooperative mechanisms of the kind described here are likely when there is more than one ligand-binding site in a protein domain.

Calculations of electrostatic interactions and pKas in the active site of Escherichia coli thioredoxin.

The active-site protonation state is crucial to the reductive mechanism of Escherichia coli thioredoxin, which involves a nucleophilic attack by the thiolate form of Cys32. We have calculated the titration properties of the active-site residues using a continuum electrostatic model, the X-ray structure of the oxidized protein, and ensembles of NMR structures of the oxidized and reduced protein. Protein dipoles, especially the SH dipole of Cys35, can provide sufficient stabilization of the Cys32 thiolate to account for its low experimental pKa (approximately 7.4), but this effect is very sensitive to local conformational variations. The experimental finding that Cys32 titrates at a lower pH than Cys35 is explained by the latter's deeper burial from solvent exposure, and stronger interaction with the carboxylate of Asp26, and not by helix dipoles or positively charged side chains. The calculated very strong interaction between Cys32 and Cys35 in their thiolate forms implies that their titration must occur in two widely pH-separated steps and that the thiolate groups must move apart in the second step. The calculations are very consistent with the experimental Asp26 pKa value of 7.5 for the oxidized X-ray structure. Both the oxidized and reduced NMR structures fall into two categories: "tight" structures in which the Asp26 and Lys57 side chains are in direct contact, and for which the calculations predict unreasonably low pKas; and "loose" structures, which resemble the oxidized X-ray structure in that these side chains are farther apart, and for which the calculations are in very good agreement with experiment. We propose that the calculations over the NMR ensemble can be used as a test of the alternative structural models provided by NMR.

All-atom empirical potential for molecular modeling and dynamics studies of proteins.

New protein parameters are reported for the all-atom empirical energy function in the CHARMM program. The parameter evaluation was based on a self-consistent approach designed to achieve a balance between the internal (bonding) and interaction (nonbonding) terms of the force field and among the solvent-solvent, solvent-solute, and solute-solute interactions. Optimization of the internal parameters used experimental gas-phase geometries, vibrational spectra, and torsional energy surfaces supplemented with ab initio results. The peptide backbone bonding parameters were optimized with respect to data for N-methylacetamide and the alanine dipeptide. The interaction parameters, particularly the atomic charges, were determined by fitting ab initio interaction energies and geometries of complexes between water and model compounds that represented the backbone and the various side chains. In addition, dipole moments, experimental heats and free energies of vaporization, solvation and sublimation, molecular volumes, and crystal pressures and structures were used in the optimization. The resulting protein parameters were tested by applying them to noncyclic tripeptide crystals, cyclic peptide crystals, and the proteins crambin, bovine pancreatic trypsin inhibitor, and carbonmonoxy myoglobin in vacuo and in crystals. A detailed analysis of the relationship between the alanine dipeptide potential energy surface and calculated protein φ, χ angles was made and used in optimizing the peptide group torsional parameters. The results demonstrate that use of ab initio structural and energetic data by themselves are not sufficient to obtain an adequate backbone representation for peptides and proteins in solution and in crystals. Extensive comparisons between molecular dynamics simulations and experimental data for polypeptides and proteins were performed for both structural and dynamic properties. Energy minimization and dynamics simulations for crystals demonstrate that the latter are needed to obtain meaningful comparisons with experimental crystal structures. The presented parameters, in combination with the previously published CHARMM all-atom parameters for nucleic acids and lipids, provide a consistent set for condensed-phase simulations of a wide variety of molecules of biological interest.

Use of 1H NMR spectroscopy and computer simulations To analyze histidine pKa changes in a protein tyrosine phosphatase: experimental and theoretical determination of electrostatic properties in a small protein.

In bovine low Mr protein tyrosine phosphatase, the pKa values of His-66 and His-72 are 8.3 and 9.2, respectively. These unusually high values were hypothesized to be caused by electrostatic interactions with several nearby negatively charged groups. To test this, mutant enzymes were made in which one or more carboxylate side chains were removed or introduced near the histidines. Michaelis kinetic parameters, measured using p-nitrophenyl phosphate as a substrate, indicated that all mutant enzymes retained approximately 50% or more of the activity of wild-type enzyme. The effect that each mutation had on the pKa of the nearby histidine was monitored by 1H NMR spectroscopy using the MLEV-17 pulse sequence to filter out the broad interfering amide resonances in the spectra. Independently, computer simulations of the pKas were obtained using the finite difference method to solve the linear Poisson-Boltzmann equation. The proximity of a charged residue to the titrating histidine imidazole largely determined the extent of the pKa perturbation. The change in pKa for His-72 in the mutant enzymes was -1.69 units for D42A, -2.36 units for E23A, -2.99 units for E23A/D42A, and unchanged for E139A and Q143E. Thus, the pKa of His-72 in the double mutant E23A/D42A decreased to nearly that of a free histidine imidazole group. The His-66 pKa change was -1.25 units for E139A and was not significant for the other mutants. His-66, Glu-139, and Gln-143 are at the protein surface and much more exposed to the higher solvent dielectric compared to His-72, Glu-23, and Asp-42. These structural characteristics explain the smaller decrease in the observed pKa of His-66 for the E139A mutant compared to the decrease in the pKa of His-72 when a single nearby carboxylate was removed. These observations were adequately predicted by theoretical electrostatic calculations using the Poisson-Boltzmann equation as a model for a solute molecule of low dielectric in solution of high dielectric.

Thermodynamics of a reverse turn motif. Solvent effects and side-chain packing.

The linear pentapeptide, Ala-Tyr-cis-Pro-Tyr-Asp-NMA (AYPYD) is known to have a significant population of type VI turn conformers in aqueous solvent. We have carried out theoretical studies of the conformational energetics of this peptide using a potential of mean force (PMF) consisting of the AMBER/OPLS empirical potential energy function, a macroscopic electrostatic model of polar solvation, and a surface area-based model of non-polar solvation. Conformers were taken from molecular dynamics simulations reported elsewhere, or generated by a random search method reported here. The chain entropy of folding was calculated by a systematic search of accessible dihedral angle space. The intra-peptide component was found to strongly favor folding and was nearly cancelled by the polar solvation term which disfavored folding. The non-polar solvation term had little effect. Fluctuations about the average value of the PMF were small and in accord with estimates from a simple harmonic model. When applied to conformers generated by a random search, the PMF selected a conformer close to the NMR-determined structure as the lowest energy conformer. The conformer with the second-lowest energy was extended, but was found to fold rapidly to the turn state in a subsequent molecular dynamics study, and may be an important state on the folding-unfolding pathway. Averages of the PMF were combined with the entropy estimates to provide an estimate of the free energy of folding that is in reasonable agreement with experimental results. In terms of the interplay between backbone electrostatic interactions and the packing of apolar side-chains, this peptide provides a model for the energetics of protein folding, and therefore makes a useful test case for calculations.

Dynamics of a type VI reverse turn in a linear peptide in aqueous solution.

BACKGROUND: Peptide sequences with aromatic groups flanking a cis-proline residue are known to have a high propensity for adopting compact structures in which the aromatic sidechains pack against the proline ring. In particular, the sequence Ser-Tyr-Pro-Phe-Asp-Val (and variants of this) is known by NMR to form a high proportion of type VI turns in aqueous solution. We set out to explore the energetic and dynamic features of such sequences using molecular dynamics simulation techniques. RESULTS: The conformation properties of the linear pentapeptide NH3(+)-Ala-Tyr-cisPro-Tyr-Asp-NMA (cis-AYPYD) have been explored in three solvated molecular dynamics simulations. The first began from an NMR-derived model structure containing a type VIa turn and close-stacking interactions between the tyrosine and proline sidechains. During 20 ns of simulation, the peptide made transitions between type VIa and VIb turns, but did not 'unfold' to more extended conformers, consistent with the unusual stability for folded forms observed by NMR for this sequence. Distances monitored by nuclear Overhauser peaks and sidechain rotamer populations in the trajectory are in good agreement with NMR data. Two additional 5 ns trajectories were begun from more extended conformers. The first folded into a conformer much like the NMR-derived structure within 3 ns and remained folded for the remainder of the trajectory. The second was begun from a structure in which the sidechain orientations were deliberately misfolded relative to that required for turn formation; this structure did not make a transition to a turn-like state. CONCLUSIONS: The kinetic stability of folded forms of AYPYD, along with the observation of spontaneous folding from an extended conformation, indicates that the special stability seen experimentally is reflected in computer simulations. The results provide new information about the stabilization of secondary structure in short peptides, particularly by aromatic-proline interactions, and offer a description of pathways of interconversion of type VIa and VIb turns.

Conformation and hydrogen ion titration of proteins: a continuum electrostatic model with conformational flexibility.

A new method for including local conformational flexibility in calculations of the hydrogen ion titration of proteins using macroscopic electrostatic models is presented. Intrinsic pKa values and electrostatic interactions between titrating sites are calculated from an ensemble of conformers in which the positions of titrating side chains are systematically varied. The method is applied to the Asp, Glu, and Tyr residues of hen lysozyme. The effects of different minimization and/or sampling protocols for both single-conformer and multi-conformer calculations are studied. For single-conformer calculations it is found that the results are sensitive to the choice of all-hydrogen versus polar-hydrogen-only atomic models and to the minimization protocol chosen. The best overall agreement of single-conformer calculations with experiment is obtained with an all-hydrogen model and either a two-step minimization process or minimization using a high dielectric constant. Multi-conformational calculations give significantly improved agreement with experiment, slightly smaller shifts between model compound pKa values and calculated intrinsic pKa values, and reduced sensitivity of the intrinsic pKa calculations to the initial details of the structure compared to single-conformer calculations. The extent of these improvements depends on the type of minimization used during the generation of conformers, with more extensive minimization giving greater improvements. The ordering of the titrations of the active-site residues, Glu-35 and Asp-52, is particularly sensitive to the minimization and sampling protocols used. The balance of strong site-site interactions in the active site suggests a need for including site-site conformational correlations.

Computational studies of the early intermediates of the bacteriorhodopsin photocycle.

Starting from a refined model of bacteriorhodopsin's ground state, alternative models of the K and L intermediates with retinal in either 13-cis or 13-14-dicis configuration have been generated by molecular dynamics simulations. All models have been submitted to electrostatic calculations in order to determine the pK1/2 values of particular residues of interest in the active site. Our pK1/2 calculations for the refined ground state can reestablish our former results, this time without adjusting the intrinsic pK of the Schiff base. For the K intermediate the electrostatic calculations show no significant change in the pK1/2 values compared to the ground state for most of the titrating groups in the active site. For the L intermediate where retinal possesses a 13-cis configuration, we found that electrostatic factors decrease the pK1/2 value of the Schiff base by 4-5 pK-units compared to the ground state. The calculations suggest that changes of the electrostatic environment via a pure 13-cis model are sufficient to produce a pK reduction of the Schiff base that will promote subsequent proton transfer steps.

Electrostatic calculations of side-chain pK(a) values in myoglobin and comparison with NMR data for histidines.

Site-specific titration curves for 12 histidine residues in carbon monoxy sperm whale myoglobin (MbCO) have been determined from two-dimensional (2D) double quantum NMR experiments. Eight of these histidine residues are observed to titrate over the accessible pH range, and pK(a) values have been determined; bounds on the titration midpoints of the remaining four histidines are also reported. Results for residues 48, 81, and 119 differ significantly from those estimated from earlier, one-dimensional studies, but they are in good agreement with values recently determined for metaquomyoglobin. These experimental values (plus those determined earlier for tyrosine titrations) are compared to predictions from crystal structures of myoglobin using a numerical Poisson-Boltzmann model and a Monte Carlo treatment of the multiple-site titration. An extension of existing models is described that accounts for alternate tautomers for histidines. Calculations are reported using several choices for radii and charges, and for five crystal structures, in order to assess the sensitivity of the results to details of the calculations. In general, the agreement between calculated and observed titration behavior suggests that this theoretical model captures much of the electrostatic behavior in this system, even though it ignores conformational fluctuations and the differences in mean structures that may exist between crystal and solution. Interactions among titrating groups are often important; in general, these interactions lead to more gradual individual site titrations (the mean Hill coefficient is about 0.8), and in several cases the interactions are so strong that two side chains need to be considered as a unit and single residues may participate in two-step titrations. It is suggested that histidines involved in such two-step titrations and carboxylic acid residues with abnormally low pK(a) values in the native conformation may be involved in the acid-induced partial unfolding of MbCO.

Electrostatic calculations of the pKa values of ionizable groups in bacteriorhodopsin.

The effects of solvation and charge-charge interactions on the pKa of ionizable groups in bacteriorhodopsin have been studied using a macroscopic dielectric model with atom-level detail. The calculations are based on the atomic model for bacteriorhodopsin recently proposed by Henderson et al. Even if the structural data are not resolved at the atomic level, such calculations can indicate the quality of the model, outline some general aspects of electrostatic interactions in membrane proteins, and predict some features. The effects of structural uncertainties on the calculations have been investigated by conformational sampling. The results are in reasonable agreement with experimental measurements of several unusually large pKa shifts (e.g. the experimental findings that Asp96 and Asp115 are protonated in the ground state over a wide pH range). In general, we find that the large unfavorable desolvation energies of forming charges in the protein interior must be compensated by strong favorable charge-charge interactions, with the result that the titrations of many ionizable groups are strongly coupled to each other. We find several instances of complex titration behavior due to strong electrostatic interactions between titrating sites, and suggest that such behavior may be common in proton transfer systems. We also propose that they can help to resolve structural ambiguities in the currently available density map. In particular, we find better agreement between theory and experiment when a structural ambiguity in the position of the Arg82 side-chain is resolved in favor of a position near the Schiff base.

pKa's of ionizable groups in proteins: atomic detail from a continuum electrostatic model.

A macroscopic electrostatic model is used to calculate the pKa values of the titratable groups in lysozyme. The model makes use of detailed structural information and treats solvation self-energies and interactions arising from permanent partial charges and titratable charges. Both the tetragonal and triclinic crystal structures are analyzed. Half of the experimentally observed pKa shifts (11 out of 21) are well reproduced by calculations for both structures; this includes the unusually high pKa of Glu 35 in the active site. For more than half the pKa's (13 out of 21), there is a large difference (1-3.3 pK units) between the results from the two structures. Many of these correspond to the titrating groups for which the calculations are in error. Since for an ionic strength of 0.1 M the Debye screening between titratable groups leads to a very high effective dielectric constant (the average value for all pairs of titrating groups is approximately 900), near-neighbor interactions dominate the pKa perturbations. Thus, the pKa values are very sensitive to the details of the local protein conformation, and it is likely that side-chain mobility has an important role in determining the observed pKa shifts.

Electrostatic effects of charge perturbations introduced by metal oxidation in proteins. A theoretical analysis.

A macroscopic dielectric model for the interactions between charges in proteins is used to calculate the changes in His residue pKa values induced in azurin by oxidation of the copper. The calculated results agree with nuclear magnetic resonance experiments to within the uncertainty associated with the measurements. It is found that a large apparent dielectric constant can describe the interaction between two protein groups, even if the shortest path between them is through the protein, which is assumed to have a low dielectric constant.

Diffusion-collision model for the folding kinetics of myoglobin.

The diffusion-collision model has been used to analyze the folding kinetics of myoglobin. The microdomains, which are the basic units that coalesce during the folding, are identified with the helices and the stabilizing contacts between helices are determined from the native structure. Both association and dissociation reactions are included and a range of stabilization parameters is investigated to determine the variation in overall rate and the relative contributions made by different intermediates during the folding process. In a comparison of folding to the native state and to the midpoint of the folding transition (i.e., 50% native protein at the completion of the reaction) significant differences in the contributing intermediates are found.

Determinants of a protein fold. Unique features of the globin amino acid sequences.

The three-dimensional structures of globins are known, from crystallographic analyses, to be very similar. Their amino acid sequences, however, differ greatly. Only two residues are absolutely conserved in all sequences, and the residue identities of some pairs of sequences are only 16%. We have determined the nature and exact extent of the sequence variations and the extent to which the conserved features of the globin sequences are unique to this family. The 226 globin sequences now known were aligned and analysed. Because distantly related protein sequences cannot be aligned correctly without the use of structural data, we developed a method that incorporated structural information into the alignment procedure. Analysis of the aligned sequences show that: (1) Although individual chains vary in size between 132 and 157 residues, deletions and insertions result in there being only 102 residue sites common to all globins. These sites form six separate regions. Insertions and deletions between these regions means that their separations can vary in different sequences. (2) Within the conserved regions there are 32 sites that almost always contain hydrophobic residues. In the known structures, these sites are in the protein interior. We measured the variations in the size of the residues that occur in the 226 sequences at these sites. At six sites the residues differ in size by less than 40 A3, at 11 sites they differ by 40 to 100 A3, and at 15 sites they differ by more than 100 A3. There are two other conserved buried sites: one contains the His linked to the haem iron and the other usually contains a His involved with the haem ligand. (3) Within the conserved regions there are another 32 sites that are almost always occupied by charged, polar or small non-polar (Gly or Ala) residues. In the known structures, these sites are on the protein surface. To determine the extent to which the conserved features found for the globin sequences are unique to that protein family, the following procedure was used. The six conserved regions, and the residue restrictions that occur at the 66 sites within these regions, were encoded into two "templates". One was based only on the sequences so far determined; the other was extended to include as yet unobserved substitutions that seemed plausible on the basis of size, hydrophobicity and polarity. Each of the 3286 non-globin sequences in the data bank was then examined by a computer program to see how closely it could be matched to these templates.(ABSTRACT TRUNCATED AT 400 WORDS)