How does the chemistry of a protein relate to its 1D sequence? In Chapter 8, we discussed techniques for detecting characteristic conserved patterns, called motifs, in families of protein sequences. We can find these sequence patterns in 1D data because although the 3D structure of a protein is complex, it is somehow determined by the invariant sequence of amino acids that makes up the protein. Motifs that are conserved in sequence often are related to important structural or functional features of a protein family, and those features often can be understood by their roles in the protein structure.

When amino acids come together in sequence to form a polymer, they do so by forming a peptide bond between the basic amino group and the acidic carboxyl group of each amino acid (Figure 9-1). This results in a long chain of amino acids that has a repeating backbone structure.

Figure 9-1. Peptide bond, peptide chain (chemical notation)

The variable group of each amino acid protrudes from the repeating backbone and is referred to in the protein structure business as a sidechain (Figure 9-2). Each of the 20 amino acid sidechains is chemically different from the others in some respect.

Figure 9-2. The amino acid sidechains (chemical notation)

The sidechains can be classified in many ways. Some are relatively large, while others are tiny or in one case nonexistent. Some have a positive or negative charge. Some are oily, or hydrophobic (water-fearing), meaning that it's energetically unfavorable for them to be solvated in water. Others are hydrophilic (water-loving), and they solvate easily in water. Some have bulky ringlike structures, while others are straight carbon chains. Some are acids, others are bases. Amino acids are conserved through evolution at specific locations in a protein sequence because they are needed there, whether to stabilize the protein structure, to form a specific binding site, or to catalyze a reaction. You can detect that particular amino acids in a protein are conserved by looking at sequence data, but to develop a hypothesis about why they are conserved, it's helpful to examine the 3D protein structure. Figure 9-3 shows the 20 amino acids classified into chemically similar groups. Note that many of the amino acids fall into more than one category. An amino acid sidechain can be both "nonpolar" and "basic," for instance, like lysine, which has a long aliphatic sidechain that terminates in an amino group. Because the relationship between chemical characteristics and amino acids isn't one-to-one, but rather many-to-many, it's not always simple to predict the effects of an amino acid substitution.

Figure 9-3. The amino acid sidechains (classification in a Venn diagram)

Interatomic forces aren't responsible only for specific interactions that form binding and interaction sites; they also are responsible for the formation of certain standard patterns that are consistently observed in protein structure. The amino acid backbone is sterically constrained—restricted from moving in certain ways because atoms will bump into each other—to follow only certain pathways. You may already be familiar with the alpha helix and beta sheet structures that commonly occur in protein structures; the reason that alpha helices and beta sheets are common is the steric restrictions on the protein backbone.

From the known structures of amino acids, Pauling and Corey first predicted the existence of alpha helices and beta sheets as a component of protein structure. Ramachandran first described exactly what range of conformations are available to amino acids in a peptide chain. Peptide chain conformation is simply described by the values of the dihedral angles in the protein backbone (i.e., the angle described by the four atoms surrounding the N-C_α bond and the angle described by the four atoms surrounding the C_α-C bond). These angles are referred to as Φ and Ψ, respectively. The chain isn't free to rotate around the third kind of bond in the protein backbone, the peptide bond, because it is a partial double bond and hence chemically constrained to be planar, so the values of Φ and Ψ for each amino acid provide a complete description of the protein backbone. A Ramachandran map is simply a plot of Φ versus Ψ for an entire protein structure. One means of evaluating a protein structure model is to compare its individual Ramachandran map with the general Ramachandran map of allowed values of Φ and Ψ.

Figure 9-4 is a general Ramachandran map that shows the allowed combinations of Φ and Ψ values for amino acids in protein structures. The small shaded region in the lower left quadrant of the map is the standard conformation of an amino acid in an alpha helix. The larger shaded region in the upper left quadrant of the map is the standard conformation of an amino acid in a beta sheet, or extended structure.

Figure 9-4. Ramachandran map of allowed conformation for protein backbones

It's apparent from the Ramachandran map that steric interactions are very important determinants of the general features of protein structure. Steric interactions instantly eliminate a large fraction of possible conformations for proteins and leave relatively few options for how a compact structure can form from a linear chain of amino acids.

The sequence of a protein is called its primary structure; the most basic level of organization in a protein is the sequence of amino acids. Alpha helix and beta sheet structures, shown in Figure 9-5, are known collectively as secondary structures and are the next level of organization. Interactions between multiple secondary structure elements give rise to supersecondary structure and tertiary structure—helices and sheets contacting each other to form larger characteristic structures, which can be described by their topology.

Figure 9-5. Alpha helix and beta strand structures

To create a functional protein, the sequence of amino acids in the protein chain must give rise to the proper 3D fold for the protein, and it must also place individual amino acids at appropriate points on that scaffold to carry out the protein's chemistry. Finding ways to extract those chemical instructions from the sequences of known proteins, formulating them as rules, and using those rules to predict the structure of other proteins is one of the biggest open research problems in bioinformatics.

Since the form that a protein structure can take and its chemical characteristics are governed by interatomic interactions, it is important to have at least a basic understanding of the interatomic interactions that play a role in protein structure. Interactions between atoms are physically complicated and to describe them in detail would require a whole other book, which fortunately has already been written by someone else: see the Bibliography. What we hope to give you is a rudimentary knowledge of these forces, to help you understand why computer methods have been developed to measure and calculate particular structural properties of proteins.

Understanding these forces gives us a basis for designing evaluative and predictive methods. Threading methods rely on the ability to discriminate between an amino acid that is in a favorable chemical environment and one that isn't. Homology modeling and structure optimization methods rely on rules for spacing between atoms, bond lengths, bond angles, and other values. These rules can be derived from chemical experiments on small molecules or from the distribution of observed values in known protein structures. However these rules are constructed, though, they reflect energetically favorable interactions between atoms.

Hydrogen bonds

Hydrogen bonds arise when two polar groups interact. The two polar groups must be of specific types. One must be a proton donor, a chemical group in which a proton (hydrogen atom) is covalently bonded to a strongly electronegative atom such as oxygen. The bond between the proton and the electronegative atom is polarized, giving the proton a partial positive charge and the electronegative atom a partial negative charge. The other group must be a proton acceptor, an electronegative atom with a partial negative charge and no attached proton. The positively polarized proton in the first group is attracted to the negatively polarized second group, and the two form a bond that isn't covalent, but is nonetheless, much shorter and stronger than a normal nonbonded interaction. Hydrogen bonds are unusual among nonbonded and electrostatic interactions because they are strongly directional; they weaken if the angle described by the three atoms involved is too large or too small.

Hydrogen bond interactions are one of the most important stabilizing forces in protein structure. The protein backbone contains a proton donor, in its N-H group, and a proton acceptor, in its carbonyl oxygen, spaced at regular intervals along the chain (Figure 9-6). The interaction of these groups stabilizes the two major types of secondary structure, the alpha helix and the beta sheet (Figure 9-7). Therefore, some structure prediction methods attempt to use the presence of potential hydrogen bond pairs to improve the accuracy of predictions.

Figure 9-6. Proton donor and acceptor in the protein backbone

Figure 9-7. Hydrogen bonding in alpha helices and beta sheets

Hydrophobic and hydrophilic interactions

A much-discussed (and frequently wrongly used) concept in protein structure analysis is that of the hydrophobic force. We've already mentioned in passing that amino acids can be classified as hydrophobic or hydrophilic. What exactly does this mean?

Proteins, except for those bound within cell membranes, always exist in aqueous solution. They constantly interact with water molecules. Water is a solution that has some interesting properties, and these properties contribute to the stability of the compact globular structures that characterize cellular proteins.

Water is a polar molecule. Individual water molecules in liquid water can each form four hydrogen bonds with neighboring water molecules. Liquid water is an essentially uninterrupted lattice of hydrogen bonded molecules, as seen in Figure 9-8. This unusual property contributes to the high melting and boiling points of water, as well as to such properties as low compressibility and high surface tension. It also results in interesting interactions of water with soluble proteins.

Figure 9-8. Hydrogen bonding in water

A nonpolar molecule dissolved in water interrupts the regular hydrogen bond lattice of liquid water. Individual water molecules can reorient around a small nonpolar molecule to preserve their network of hydrogen bonds, but this reorientation has a cost in terms of free energy (which is how cost is measured in chemistry). The presence of a nonpolar solute forces water molecules into a more ordered conformation than they would ordinarily assume. Instead of being able to face any which way and rotate freely, water molecules near the surface of a nonpolar solute have to work around it and form a cage. This is entropically unfavorable.

The larger a nonpolar solute gets, the more water molecules need to reorient to accommodate it, and the higher the energy cost of solvating the molecule becomes. Of course, if the nonpolar solute has some polar groups on its surface, water molecules can use those groups as hydrogen bonding partners instead of other water molecules, and the water lattice is less disturbed. Globular proteins, which exist in aqueous solution even though they are composed substantially of nonpolar groups, must present a good hydrogen-bonding surface to the world. Hydrophilic amino acids are those whose sidechains offer hydrogen bonding partners to the surrounding medium, while hydrophobic amino acids' sidechains don't. The surface of a globular protein is usually anywhere from 50%-75% polar atoms, and deviations in this pattern can suggest binding or complexation sites.

Solvent accessibility and hydrophobicity play an important role in evaluating model structures. Threading methods for protein fold recognition use amino acid environments in evaluating models. When many hydrophobic amino acids are found in solvent-exposed structural environments or hydrophilic amino acids buried in the protein interior, it is considered unlikely that the protein model is folded correctly.

Relative strength of interatomic forces

The interaction between atoms can be described by a pair potential, such as the Lennard-Jones potential (Figure 9-9), which includes both an attractive and a repulsive term. The form of the potential shows that atoms tend to repel each other at very short range (positive potential energy indicating an unfavorable interaction) but to attract each other at slightly longer range. The strength of the attraction decays with distance, depending on the forces modeled.

Figure 9-9. Plot of Lennard-Jones potential

When making inferences about structural stability or function based on intermolecular interactions, it is important to understand the relative strengths of these interactions, and how they scale with distance (Table 9-1).

Table 9-1. How Interatomic Forces Scale with Distance

Type of Bond	Range of Interaction
Covalent	Complicated short range
Hydrogen bond	Roughly 1/r 2
Charge-charge	Scales with 1/r
Charge-fixed dipole	Scales with 1/r 2
Charge-rotating dipole	Scales with 1/r 4
Fixed dipole-fixed dipole	Scales with 1/r 3
Rotating dipole-rotating dipole	Scales with 1/r 6
Charge-nonpolar	Scales with 1/r 4
Dipole-nonpolar	Scales with 1/r 6
Nonpolar-nonpolar	Scales with 1/r 6

In Table 9-1, r represents the distance between two atoms in angstroms. Interactions that decrease in strength with 1/r are effective at a much longer range than those that decrease in strength with higher powers of r. Covalent interactions and hydrogen bonds are strong, and very energetically significant at short distances. Charge-charge interactions have some of the longest-range effects; electrostatic effects on protein activity have been experimentally shown at over 15-angstrom distance, a substantial range in molecular terms. A concentration of charges on a protein surface can create a powerful electrostatic steering effect that can attract ligand molecules or other proteins at even longer range. Hydrogen bonds and charge-dipole interactions are also relatively strong. The effects of these interactions are modeled by computing electrostatic potentials and using the computed potentials as the basis for calculating other molecular properties such as binding constants (via Brownian dynamics) or pKa values.

On the other hand, interactions between noncharged and nonpolar atoms are very weak and effective only at short range. However, the effects of these interactions can be cumulative, stabilizing structure and making intermolecular associations more favorable. The effects of these interactions are addressed when you compute the size of intermolecular contact surfaces or enumerate interactions between neighboring interactions in a protein. In the remainder of this chapter, we discuss various methods for measuring and evaluating atomic structures of proteins, all of which can be used together to add to your understanding of protein chemistry.

One type of molecular structure viewers are lightweight applications that can be set up to work with your web browser. When properly configured, they will display molecular data as you access it on the Web. RasMol and CnD3 are two of the most popular viewers.

# DEPTHDEF = -DTHIRTYTWOBIT
DEPTHDEF = -DSIXTEENBIT
# DEPTHDEF = -DEIGHTBIT

Description:		Brookhaven PDB 
MIMEType:		chemical/x-pdb 
Suffixes:	.pdb 
Application:		/usr/local/bin/rasmol

Description:		NCBI ASN.1 
MIMEType:		chemical/ncbi-asn1-binary 
Suffixes:	.prt 
Application:		/usr/local/cn3d/Cn3D

Heavy-duty molecular structure viewers tend to have many more features than web applications such as RasMol and Cn3D. The most popular examples are MolMol, MidasPlus, and VMD. These programs run on your desktop machine, and to use them you need copies of the PDB files you're interested in using already stored on your computer.

MolScript has a completely different purpose from the other visualization packages we have discussed. It is designed to produce high-quality graphics for print publication, as you can see in Figure 9-10. It can be configured to run from the command line and to produce PostScript, Raster3D, and VRML output only; it can also be configured to run interactively in its own window, using OpenGL, and to produce output in many additional image file formats.^[*]

Figure 9-10. A sample image generated by molscript

Setting up interactive MolScript with OpenGL on a Linux workstation isn't straightforwRasMolard; it requires the installation of Mesa (open source OpenGL) libraries and customization of the Makefile that comes with the distribution. However, the basic MolScript installation is quick and simple and can produce visually appealing line drawings of molecular structure cartoons in color or black and white, in a style that is uniquely elegant and appropriate for print media. To install the basic version of MolScript, simply follow the directions in the install file. Copy the resulting executables (molscript and molauto) to your /usr/local/bin directory or to another directory in your default path. Here's what molscript and molauto do:

MolScript takes two input files: a MolScript command file and a PDB coordinate file. Here's the MolScript input file that produced the images in Figure 9-10:

! MolScript v2.1 input file 
! generated by MolAuto v1.1.1
title "MYOGLOBIN  (FERRIC IRON - METMYOGLOBIN)"
plot
  read mol "1MBN.pdb";
  transform atom * by centre position atom *;
  set segments 2;
  set planecolour hsb 0.6667 1 1
coil from 1 to 3
set planecolour hsb 0.619 1 1
helix from 3 to 18
set planecolour hsb 0.5714 1 1
coil from 18 to 20
set planecolour hsb 0.5238 1 1
helix from 20 to 35
...
coil from 94 to 100
set planecolour hsb 0.1429 1 1
helix from 100 to 118
set planecolour hsb 0.09524 1 1
coil from 118 to 125
set planecolour hsb 0.04762 1 1
helix from 125 to 148
set planecolour hsb 0 1 1
coil from 148 to 153;

set colourparts on
bonds in require residue 1 and type HEM;

end_plot

The MolScript scripting language is unique and not really based on any standard computer language. The only way to learn it is to decide what you want to do, study the manual and examples, and learn the language. The example just shown is a simple MolScript command file; it reads in a single molecule, centers it on the molecule's center of mass, defines the locations of the various secondary structure elements and shades them through the spectrum from red to blue. MolScript can produce much more complex figures than this, however. MolScript plots can be scaled and multiple plots shown on a single page. Subsets of atoms in the molecule can be turned on, displayed in different formats, and custom colored. Labels can be added to figures.

Fortunately, the molauto program automatically produces simple input files for the molscript program, which can help you get started using the MolScript command language. molauto does the most tedious part of input file setup for you—assigning helix, sheet, or coil drawing styles, and colors, to each segment of secondary structure. molauto has a variety of command-line options, which you can access by entering molauto -h. molauto reads input in the standard PDB file format, and writes to standard output unless a redirector is used.

The following are some of the most useful command line options for molauto:

The output of the molauto program is an input for the main molscript program. Command-line options for molscript include:

The default input files produced by molauto can be hand-edited to produce various effects. One important thing you might want to do (and can't do automatically unless you have installed the MolScript package with OpenGL support) is to rotate the molecular structure until you achieve a good view.

To rotate the molecule view using the noninteractive version of molscript, add the following lines to your molscript input file, replacing the line that currently reads:

transform atom * by centre position atom *; 

with:

transform atom * by centre position in amino-acids
                 by rotation x    0.0
                 by rotation y    0.0
                 by rotation z    0.0
                                    ;  !Be sure to include this semicolon.

After you generate your first version of the image, open it in a fast PostScript viewer such as gv. To change the view of the molecule, experiment with changing the values of x, y, and z rotation in your input file. Since molscript takes only seconds to run on any protein input file, you can make changes to the input file, save the file, and redisplay the new output several times until you like the view.

Once generated, the molscript image file can be viewed, converted to other file formats, and edited using standard Unix image-manipulation tools. One program you can load when you install most major Linux distributions is GIMP, the freeware package similar to Adobe Photoshop.

Another useful tool for producing graphics for publication is the program LIGPLOT (http://www.biochem.ucl.ac.uk/bsm/ligplot/ligplot.html), which is available from the Structure and Modelling group at University College London (UCL). Given a molecular structure and a specific residue or heteroatom group within the structure as input, LIGPLOT automatically generates a 2D schematic drawing showing hydrogen bonds, interatomic contacts, and solvent accessibility. A sample of LIGPLOT is shown in Figure 9-11

Figure 9-11. A schematic diagram of ligands to the heme cofactor in cytochrome B5, generated with LIGPLOT

To install LIGPLOT on a Linux workstation, simply follow the directions in the README file.

In order for LIGPLOT to find its parameter files and helper programs correctly, you need to add some path information to your .cshrc file:

setenv ligdir /usr/local/ligplot
alias ligplot $ligdir'/ligplot.scr'
alias ligonly $ligdir'/ligonly.scr'
alias dimplot $ligdir'/dimplot.scr'
alias dimonly $ligdir'/dimonly.scr'
setenv hbdir /usr/local/hbplus
alias hbplus $hbdir'/hbplus' 

The values on the command line specify a residue range in a particular protein chain. The program doesn't have to display only interactions with ligands and prosthetic groups; it can also display the network of close interactions with any residue in a protein. This works best when the residue range selected is small.

The dimplot program, a variant of LIGPLOT, displays interactions across an interface between two protein chains or domains. The domain variant works only if your PDB file labels proteins at the domain level of organization.

The painful part of installing the LIGPLOT, hbplus, and naccess programs on some Linux systems is, ironically, not the installation itself, but having the capability to decrypt the encrypted archives you get from UCL. The files are encrypted using the standard Unix crypt command. This sounds straightforward enough, but many Linux vendors don't include crypt in their distributions. In order to use crypt on your system, you may in fact need to reinstall the latest version of glibc-2.0. If you don't want to deal with this, request a decrypted copy of the LIGPLOT tar archive from the authors when you send in your license agreement.

Protein coordinate data sets don't automatically come labeled with alpha-helix and beta-sheet classifiers. Secondary structure features in the protein can be distinguished with reasonable certainty by their hydrogen bonding patterns and their backbone torsion angles.

The standard program for extracting secondary structure from sequence is the DSSP program. DSSP analyzes the geometry and backbone hydrogen bonding partners of each residue in a known protein structure, producing a tabular output that includes residue numbering, sequence, hydrogen bonding, and geometry details. The DSSP database, and DSSP executables derived from the 1995 release of the program, are available from the European Bioinformatics Institute (EBI); these executables may still cause Y2K-related errors on some older Linux systems. Updated DSSP source code is available from the Gerrit Vriend at the Center for Molecular and Biomolecular Informatics at the University of Nijmegen, Netherlands.

ACC  ALA -  143  142 ->  TYR -  146  145  3.3  107.8  125.8   58.5   76.9  1MBN 

ACC  ALA -  143  142 ->  LYS -  147  146  3.2  154.3  113.4    0.1   43.4  1MBN 
DNR  ALA -  144  143 ->  LYS -  140  139  3.0  153.6  109.9   16.4   27.2  1MBN 

ACC  ALA -  144  143 ->  GLU -  148  147  3.0  160.3  109.4   11.6    6.4  1MBN 
DNR  LYS -  145  144 ->  ASP -  141  140  3.2  145.3  119.5    3.7   73.8  1MBN 

ACC  LYS -  145  144 ->  LEU -  149  148  3.0  149.4  128.8    4.7   63.7  1MBN 
DNR  TYR -  146  145 ->  ILE -  142  141  3.2  158.7  121.8   20.1   52.6  1MBN 
DNR  TYR -  146  145 ->  ALA -  143  142  3.3  107.8  125.8   58.5   76.9  1MBN 

ACC  TYR -  146  145 ->  GLY -  150  149  3.0  156.9   96.3   37.1   37.7  1MBN 

ACC  TYR -  146  145 ->  TYR -  151  150  3.1  111.2  118.0    4.2   89.9  1MBN 
DNR  LYS -  147  146 ->  ALA -  143  142  3.2  154.3  113.4     0.1   43.4  1MBN 

Topology cartoons are a 2D notation for depicting the topological arrangement of secondary structural elements in proteins. The cartoons can clarify the spatial relationships and connectivity between secondary structure elements in a protein. These relationships may not be easily seen in a 3D structure, even if only the structural backbone is displayed or a ribbon diagram is drawn. Software for generating your own cartoons may be found on the Protein topology page, http://www.sander.embl-ebi.ac.uk/tops/.

Topology cartoons, as illustrated in Figure 9-12, represent each secondary structural unit as a shape. Circles are helices, and triangles are beta strands. The beginning of the chain is marked with an N, the end with a C. Each element has a directionality, which can be deduced from the way the connecting segment is drawn. If the N-terminal connection is to the edge of the secondary structural element, that element is directed out of the plane of the drawing; if the N-terminal connection is to the center of the secondary structural element, it is directed back into the plane of the drawing.

Figure 9-12. A protein topology cartoon

setenv PATH "/usr/sbin:/sbin:/usr/jdk118/bin:${PATH:."
setenv CLASSPATH "/usr/local/Tops/classes/TOPS.jar:${CLASSPATH"
setenv TOPS_HOME "/usr/local/Tops"

Classification databases are taxonomies of protein structure, and they bear a strong resemblance to the morphology-based taxonomies developed by early biologists. Proteins that "look" grossly the same, in terms of shape and topology, are classified as more closely related than proteins that look substantially different. Protein structure types have whimsical names (like Greek key beta barrel ) based on visual observation and comparison with familiar objects. The classification databases can be envisioned as trees with many branchings at each branch point—very similar to phylogenetic trees, in concept.

The most common parameter that expresses the difference between two protein structures is RMSD, or root mean squared deviation, in atomic positions between the two structures. RMSD can be computed as a function of all the atoms in a protein or as a function of some subset of the atoms, such as the protein backbone or the alpha-carbon positions only. Using a subset of the protein atoms is common, because it is likely that, when two protein structures are compared, they will not be identical to each other in sequence, and therefore the only atoms between which one-to-one comparisons in position can be made will be the backbone atoms.

This is the first context we've discussed in which the orientation of a molecular structure becomes important. Because protein structures are generally described in Cartesian coordinates, they essentially exist within a virtual space, and they come with a built-in orientation with respect to that space. RMSD is a function of the distance between atoms in one structure and the same atoms in another structure. Thus, if one molecule starts out in a different position with respect to the reference coordinate system, the other molecule—the RMSD between the two proteins—will be large whether they are similar or not.

In order to compute meaningful RMSDs, the two structures under consideration must first be superimposed, insofar as that is possible. Superimposition of protein structures usually starts with a sequence comparison. The sequence comparison establishes the one-to-one relationships between pairs of atoms from which the RMSD is computed. Atom-to-atom relationships, for the purpose of structure comparison, may actually occur between residues that aren't in the same relative position in the amino acid sequence. Sequence insertions and deletions can push two sequences out of register with each other, while the core architecture of the two structures remains similar.

Once atom-to-atom relationships between two structures are established, the task of a superposition program is to achieve an optimal superposition between the two programs—that is, the superposition with the smallest possible RMSD. Because protein scaffolds, or cores, can be similar in topology without being identical, it isn't usually possible to achieve perfect overlap in all pairs of atoms in two structures that are being compared. Overlaying one pair of atoms perfectly may push another pair of atoms further apart. Superposition algorithms optimize the orientation and spatial position of the two molecules with respect to each other.

Figure 9-13 shows an optimal alignment between atomic structures of triosephosphate isomerase and beta-mannase, shown in Compare3D. The two structures are similar enough to be classified as structural neighbors, and their chain traces are relatively similar. However, their sequence identity is only 8.5%.

Figure 9-13. An optimal superposition of myoglobin and the 4 chain of hemoglobin, which are structural neighbors

Once optimal superpositions of all pairs of structures have been made, the RMSD values that are computed as a result can be compared with each other, because the structures have been moved to the same frame of reference before making the RMSD calculations.

The DALI Domain Dictionary (DDD) at the EBI is based on an automatic classification of protein domains by sequence identity. Rather than using a human-designed classification scheme, DDD is constructed by clustering protein neighbors within an abstract fold space. Instead of working with whole proteins, DDD classifies structures based on compact, recurring structures (called domains ) that may repeat themselves within, and among, different protein structures. The content of DDD may also be familiar to you as FSSP, the "Fold classification based on Structure-Structure alignment of Proteins" database.

DDD can be searched based on text keywords; it can also be viewed as a tree or a clickable graphical representation of fold space. Views of sequence data for conserved domains are available through the DDD interface, as well as connections to structural neighbors.

The superposition program (SUPPOS) that produces the structural alignments in DALI/FSSP is available within the WHAT IF software package of protein structure analysis tools, which is discussed in Section 9.7.1.2.

The Combinatorial Extension of the Optimal Path (CE) is a sophisticated automatic structure alignment algorithm that uses characteristics of local geometry to "seed" structural alignments and then joins these regions of local similarity into an optimal path for the full alignment. Dynamic programming can then optimize the alignment.

CE is available either as a web server or as source code from the San Diego Supercomputer Center. The web server allows you to upload files for pairwise comparison to each other or to proteins in the PDB, to compare a structure to all structures in the PDB, to compare a structure to a list of representative chains, and review alignments for specific protein families. CE also is fully integrated with the PDB's web site, and CE searches can be initiated directly from the web page generated for any protein you identify in a sequence search. Along with the source code, you can download a current, precomputed pairwise comparison database containing all structures in the PDB. If you're doing only a few comparisons, however, you probably won't even want to do this.

When using the CE server to compute similarities, there are several parameters that you can set, including cutoffs for percent sequence identity, percent of the alignment spanned by gaps, and percent length difference between two chains. You can also set an RMSD cutoff and a Z-score cutoff. The Z-score is a measure of the significance of an alignment relative to a random alignment, analogous to a BLAST E-value. A Z-score of 3.5 or above from CE usually indicates that two proteins have a similar fold.

Along with CE, the SDSC offers the Compound Likeness (CL) server, a suite of tools for probabilistic comparison between protein structures. In CL, you select either an entire protein structure or a structure fragment to use as a probe for searching the PDB. Search features include bond length and angle parameters, surface polarity and accessibility, dihedral angles, secondary structure, shape, and predicted alpha helix and beta sheet coefficients. CL allows you to ask the question "what else is chemically similar to this protein (or fragment) that is of interest to me" and to define chemical similarity very broadly. A full tutorial on CL is available at the CL web site (http://cl.sdsc.edu/cl1.html/ ).

VAST is a pairwise structural alignment tool offered by NCBI. VAST reports slightly different parameters about structural comparison than CE does, and the underlying algorithm differs in significant respects. However, the results tend to be quite similar. VAST searches automatically allow you to view your superimposed protein structures in the Cn3D browser plug-in, with aligned sequences displayed in Cn3D as well. For practical purposes, either CE or VAST is sufficient to give you an idea of how two structures match up; if you are concerned about the algorithmic differences, both groups provide access to detailed explanations at their sites. Unlike CE, the VAST software doesn't appear to be available to download, so if you want to perform a large number of comparisons on your own server, CE may be preferable.

Geometric analysis can show where a model developed from x-ray crystallography data or NMR data violates the laws of chemistry. As mentioned earlier, there are physical laws governing intermolecular interactions: nonbonded atoms can get only so close to each other because as two atoms are forced together beyond the boundary set by their van der Waals radius, the energetics of the contact become very unfavorable. These interactions limit not only the contacts between pairs of atoms in different parts of a protein chain, but also how freely atoms can rotate around the bonds that connect them. The structure of atomic orbitals and the nature of bonds between atoms place natural limits on the position of bonded atoms with respect to each other, so bond angles and dihedral angles are, in practice, restricted to a limited set of values. Tools for geometric analysis generally have been developed by crystallographers to show where their structural models violate these laws of nature; they also can be used by homology modelers or ab-initio structure modelers to evaluate the quality of a structural model.

There are a variety of tools for analyzing structure quality. Some run as standalones; others are incorporated into more comprehensive structure analysis and simulation packages. An exhaustive listing of the best of these tools can be found on the PDB web site.

setenv prodir /usr/local/procheck
alias procheck	'$prodir /procheck.scr'
alias procheck_comp	'$prodir /procheck_comp.scr'
alias procheck_nmr	'$prodir /procheck_nmr.scr'
alias proplot	'$prodir /proplot.scr'
alias proplot_comp	'$prodir /proplot_comp.scr'
alias proplot_nmr	'$prodir /proplot_nmr.scr'
alias aquapro	'$prodir /aquapro.scr'
alias gfac2pdb	'$prodir /gfac2pdb.scr' 
alias viol2pdb	'$prodir /viol2pdb.scr' 

Colour all plots?
-----------------  
Y <- Produce all plots in colour (Y/N)?

Which plots to produce
----------------------
Y	<-  1. Ramachandran plot (Y/N)?
N	<-  2. Gly & Pro Ramachandran plots (Y/N)?
N	<-  3. Chi1-Chi2 plots (Y/N)?
N	<-  4. Main-chain parameters (Y/N)?
N	<-  5. Side-chain parameters (Y/N)?
N	<-  6. Residue properties (Y/N)?
N	<-  7. Main-chain bond length distributions (Y/N)?
N	<-  8. Main-chain bond angle distributions (Y/N)?
N	<-  9. RMS distances from planarity (Y/N)?
N	<- 10. Distorted geometry plots (Y/N)?

Geometric analysis can also be useful in understanding a protein's fold and function. In this case, the geometry of interest isn't the chemical bonding interactions between atoms adjacent to each other in the protein chain, but the nonbonded interactions between atoms that are widely separated in the protein chain. The density of intramolecular contacts in the structural core of a domain may be quite different from the density of contacts in a region between two structural domains. Measuring this density over the whole protein may give clues as to the process by which a protein folds. The patterns of hydrogen bonds that hold a protein together may serve as an identifying signature for a protein fold. And contacts between certain chemically important residues in a protein may suggest hypotheses about the protein's catalytic mechanism or function. Protein engineers may want to examine the intramolecular contacts in a protein to determine where changes are least likely to disrupt the protein's structure.

There are many programs for calculating solvent accessibility by probe-rolling. They are all straightforward and easy to use, requiring a standard PDB file as input and usually giving output in the form of a percentage of accessibility for each amino acid or atom in the protein. One popular program is naccess , which is available from the Biomolecular Structure and Modelling group at UCL. naccess can be used in combination with other programs developed by this group, such as HBPLUS (a program for computing intermolecular interactions and hydrogen bonds) and LIGPLOT, which we covered earlier. It also runs as a standalone. naccess is written in FORTRAN, so you'll need the g77 compiler installed on your machine to compile it.

naccess outputs two files, an .asa file with accessible surface areas for each atom in the molecule and an .rsa file with accessible surface areas and relative accessibilities for each amino acid. It handles both protein and nucleic acid molecules and produces reports of accessibilities for individual molecular chains as well as complete structures. The -h, -w, and -y flags cause the program to ignore hydrogen atoms, water molecules, or heteroatoms, respectively. Run with the -c option, naccess produces intermolecular contact areas rather than accessible areas.

SURFNET is a program developed by Roman Laskowski at UCL to manipulate solvent-accessibility information and display useful representations of surface features, clefts, cavities, and binding sites. SURFNET generates surface output in formats that can be displayed in molecular visualization programs, including RasMol.

The Alpha Shapes software is a mathematically exact alternative to the standard Connolly surface method of computing solvent accessibility. Developed by the research group of mathematician Herbert Edelsbrunner at NCSA (http://www.alphashapes.org/alpha/), the Alpha Shapes software is a general-purpose program for modeling the surfaces of objects. A set of extensions to the Alpha Shapes software, specifically for analyzing protein molecules, is also available.

The Alpha Shapes method constructs what is called a simplicial complex or alpha complex of a structure, based on a rigorous geometrical decomposition of the space surrounding the collection of points that describes the structure. Once the alpha complex is constructed for a protein structure, algorithms for inclusion and exclusion can describe exactly the surface area or volume of the structure as well as cavities, clefts, and regions of contact. The main benefit of using the Alpha Shapes algorithm to compute protein shapes is that the software comes with a sophisticated visualization program called alvis, which can display such geometrical features as the interior shape of an ion channel or the cavities in the interior of a protein.

Several programs make up the Alpha Shapes distribution. These programs must be run in the proper sequence to correctly analyze molecular data:

You can find usage details of each of these programs in a README file that accompanies the Alpha Shapes distribution, or by attempting to run the program with no arguments on the command line.

Using alvis to visualize your Alpha Shapes data can be quite interesting. To do this, you need output from delcx and mkalf, but not from VOLBL. To run alvis on a data set generated from molecule.PDB, where output files molecule.dt and molecule.alf are also present, enter alvis molecule. The visualizer opens with the convex hull of the molecule displayed. The standard atomic structure of the molecule can't be seen from within the current version of alvis, but you can compare your alvis view with another view of the same molecule (perhaps using RasMol or a similar molecular visualization program) side by side.

In the bottom left of the alvis control panel, you'll see a box containing a graph with three colored curves. This graph is called the alpha rank graph, and it can be used to select a desirable view of the molecule. Positioning your cursor at peaks, valleys, or intersections on these graphs gives the most interesting views of the molecular shape.

Using the Pocket Panel, available from the Gizmos pulldown in alvis, you can make selections that shows voids, pockets, and difference areas in a protein. The online alvis tutorial at http://www.alphashapes.org describes in full the settings needed to view these features.

Along with the main Alpha Shapes programs, a number of utility scripts are provided that can postprocess VOLBL output to give specific information. These include:

Analogous scripts for computing volume differences are also included.

Analytical surface potentials based on Alpha Shapes can also be accessed with the CAST-P web server at the University of Illinois at Chicago. At the time of this writing, not all protein structures in the PDB are represented on CAST-P; the site is currently under development. However, it promises to be a useful analytical tool in the future. CAST-P features an integrated Java-based visualization of cavities in protein structures and the amino acids that are in contact with cavities.

A protein molecule can be thought of as a collection of charges in a dielectric medium. In the model that computes electrostatic potentials for protein molecules, each atom is represented as a point with a partial atomic charge. The solvent accessible surface of the protein forms the boundary between the interior medium of the protein and the exterior medium surrounding the protein.

Computing the electrostatic potential of a protein structure allows you to predict quantities such as individual amino acid pKa values, solvation energies, and approximations to intermolecular binding energies. If you are interested in protein modeling, macromolecular electrostatics is a topic that you may wish to explore further. Our review of the subject in the March 2000 issue of Methods provides an entry point to the molecular electrostatics literature.

The University of Houston Brownian Dynamics (UHBD) package is a state-of-the-art software package for computing macromolecular electrostatics. UHBD computes electrostatic potentials and can also use those potentials as parameters in subsequent Brownian Dynamics and Molecular Dynamics simulations. The most recent release of UHBD can be compiled on Linux systems and includes several control scripts that implement UHBD to calculate pKa values for individual titrating amino acids in the protein, as well as theoretical titration curves for the protein as a whole. UHBD is accessed by a scripting language; it requires a protein structure file and a command script as input. It also requires a file containing atomic partial charges for any amino acids and other atoms in the input structure. Detailed scripting examples are provided in the UHBD distribution, along with charge datafiles that allow the program to assign correct partial atomic charges to all but unusual atom groups.

UHBD, and other similar programs such as DelPhi—which overlaps only the electrostatics functionality of UHBD—use numerical approximations to solve the Poisson-Boltzmann equation for the large number of interacting charges that make up a protein. In the finite-difference approach used by UHBD, the irregularly spaced charges in a protein molecule are mapped onto a regular 3D grid, and the Poisson-Boltzmann equation is solved iteratively for each point on the grid until the solution converges to an electrostatic potential for each point.

Other than alvis, which doesn't truly display a molecular surface but rather a mathematically derived pseudosurface, there are several options for displaying the shapes of molecules. Most molecular modeling packages incorporate a molecular surface display feature and allow the surface to be colored according to chemical properties. However, the display schemes in programs not specifically designed for that purpose are too unsophisticated to handle data from macromolecular electrostatics calculations and other representations of physicochemical properties. An exception seems to be the SWISS-PDBViewer (discussed earlier), which can interpret data from external electrostatics calculations and analytical molecular surface calculations.

What are the "ideal" parameters or constraints used in optimization? In some cases, they are based entirely on chemical principles: bond lengths and angles determined by steric restrictions and nonbonded interactions described as Lennard-Jones potentials. In other cases, structural constraints are based on information derived from the database of known protein structures. If a particular amino acid in a particular sequence context always has the same conformation, a higher probability can be assigned to it assuming that conformation again, rather than a different conformation. Secondary structure prediction methods use an information-based approach to predicting likely conformations for the protein backbone. Optimization methods use information to refine atomic structures at the level of individual sidechain atoms once the backbone trace has been worked out.

Rotamer libraries are parameter sets specifically for the optimization of sidechain positions in molecular model building. They are called rotamer libraries because they contain information about allowed rotations of the remote amino acid sidechain atoms around the Cα - Cβ bond, expressed as the allowed values of sidechain dihedral angles.

Because of steric constraints on bond rotation, amino acid sidechains in proteins can assume only a few conformations without unfavorable energetic consequences. Rotamer libraries can be derived using chemical bond and angle constraints, but, they are more likely to be developed by analysis of the conformations assumed by amino acid sidechains in known protein structures. Rotamer libraries can be either backbone-dependent or -independent. Backbone-independent rotamer libraries classify all instances of a particular amino acid as part of the same set, even if one occurrence is within a beta sheet and the other is within an alpha helix. Backbone-dependent rotamer libraries, on the other hand, further classify amino acids according to their occurrence in specific secondary structures.^[†]

SCRWL, available from the Fred Cohen research group at UCSF, is a program that allows you to model sidechain conformations using a backbone-dependent rotamer library.

The derivation of probability density functions (PDFs) is similar in concept to the development of rotamer libraries, although more mathematically rigorous. The essence of a PDF is that a mathematical function is developed to represent a distribution of discrete values. The discrete values that make up the distribution are harvested from occurrences of a situation in a representative database of samples. That mathematical function can then be used to evaluate and optimize (and in some cases even predict) the properties of future occurrences of the same situation.

In protein modeling, PDFs have been used to describe intra- and inter-residue interatomic distances, as well as bond angles, dihedral angles, and other more spatially extensive regions of protein structure. Modeller, which is discussed in Chapter 10, uses a combination of bond angle and dihedral angle PDFs to optimize the protein structure models it builds. Modeller's internal OPTIMIZE routine can be used for PDF-based structure optimization.

The data from which PDFs are generated can be broken down into specific occurrences; for example, all contacts between Cβ in residue i and Cβ in residue i+4 when both residues are leucine but again, trade-offs between classification detail and class population occur. Distance PDFs for proteins have been used by several groups to evaluate and optimize protein structures. Most such work is still in its early stages, and software isn't yet available for public use.

Figure 9-15 shows a plot of a distance probability density function for tertiary interactions between sulfur atoms in cysteine residues generated from known protein structures. The function's peak near 2 angstroms corresponds to the high propensity with which the sulfur atoms form disulfide bridges between cysteine residues. These data are taken from the Biology Workbench at the San Diego Supercomputer Center (http://workbench.sdsc.edu/) and plotted using xmgr. Note that the Workbench PDFs make a distinction between cysteine residues participating in disulfide bridges (pictured here and referred to as CSS residues at the Workbench site) and those cysteines that don't participate in disulfide bonds (which the Workbench site calls CYS).

Figure 9-15. Interatomic distance probability density function

Structure evaluation based on PDFs is implemented in the Structure Tools section of the SDSC Biology Workbench (Figure 9-16). You can upload a PDB structure or a theoretical model and score the structure either on a residue-by-residue or an atom-by-atom basis. Scores can be displayed on a plot, where the Y-axis represents the relative probability of the region of structure that's being evaluated. This can be thought of in terms of the probability that a particular residue or atom is in the "correct" position, given what is known about other occurrences of that residue or atom in similar sequence environments. Regions with low probability are likely to be misfolded or poorly modeled. PDF probability scores can also be written out in a special PDB file, in place of the temperature factor values found in the original PDB file. These special PDB files can then be displayed using a visualization program such as RasMol or Chime. Coloring the molecule by temperature factor maps the PDF probability scores onto the molecular structure, highlighting regions of the structure that score poorly.

Figure 9-16. Comparing PDF scores for an obsolete PDB structure (1B5C) and (1CY0) that superseded it

The GeneCensus project is a broad, sequence-based comparison of the protein content of several genomes. At the time of this writing, GeneCensus contains information from 20 genomes. GeneCensus currently can be searched with a PDB ID or an ORF ID to locate occurrences of specific protein folds in the GeneCensus genomes.

PRESAGE is a database of information about experimental progress with the structures of various proteins. You can search PRESAGE with a TIGR ORF code, NCBI or SWISS-PROT accession number, and a number of other codes to find out if someone is attempting to isolate, crystallize, and solve the structure of that protein, and if so, how far along they are. PRESAGE is relatively new at press time; it currently contains only about 6,000 records and isn't guaranteed to be comprehensive. It can't be searched directly by BLAST or FASTA search with a sequence, so before checking PRESAGE for instances of an unknown sequence, you have to search for matching accession codes. However, in principle, the PRESAGE database promises to be useful, both for crystallographers and their collaborators, and for curious citizens of the molecular biology community who are wondering if the structure of their protein has been solved yet.

The Biomolecular Interaction Network Database (BIND) is another relatively new data repository offered by the Samuel Lunenfeld Research Institute. BIND was designed to be a central deposition site for known information about macromolecular interactions. BIND is a complex database, containing information about interactions between objects in the database, molecular complex information, and metabolic pathway information. The BIND format is designed to contain information about experimental conditions that observe the interactions stored in the database, as well as information about binding site location, biochemical activity, kinetics, and thermodynamics. BIND is still in its beta testing phase, containing only a few hundred interactions. However, BIND has been funded to hire indexers and it is expected to grow rapidly in the near future. One interesting aspect of the BIND data entry process is that methods for automated reading of existing journal literature are being implemented to extract known interactions from their inconvenient location in dusty library stacks and put them more effectively in the public domain.