Comparative Protein Structure Modeling Using MODELLER

Functional characterization of a protein sequence is one of the most frequent problems in biology. This task is usually facilitated by accurate three‐dimensional (3‐D) structure of the studied protein. In the absence of an experimentally determined structure, comparative or homology modeling can sometimes provide a useful 3‐D model for a protein that is related to at least one known protein structure. Comparative modeling predicts the 3‐D structure of a given protein sequence (target) based primarily on its alignment to one or more proteins of known structure (templates). The prediction process consists of fold assignment, target‐template alignment, model building, and model evaluation. This unit describes how to calculate comparative models using the program MODELLER and discusses all four steps of comparative modeling, frequently observed errors, and some applications. Modeling lactate dehydrogenase from Trichomonas vaginalis (TvLDH) is described as an example. The download and installation of the MODELLER software is also described. Curr. Protoc. Bioinform. 47:5.6.1‐5.6.32. © 2014 by John Wiley & Sons, Inc.


INTRODUCTION
Functional characterization of a protein sequence is one of the most frequent problems in biology. This task is usually facilitated by an accurate three-dimensional (3-D) structure of the studied protein. In the absence of an experimentally determined structure, comparative or homology modeling often provides a useful 3-D model for a protein that is related to at least one known protein structure (Marti-Renom et al., 2000;Fiser, 2004;Misura and Baker, 2005;Petrey and Honig, 2005;Misura et al., 2006). Comparative modeling predicts the 3-D structure of a given protein sequence (target) based primarily on its alignment to one or more proteins of known structure (templates).
Comparative modeling consists of four main steps (Marti-Renom et al., 2000;Figure 5.6.1): (i) fold assignment, which identifies similarity between the target and at least one known template structure; (ii) alignment of the target sequence and the template(s); (iii) building a model based on the alignment with the chosen template(s); and (iv) predicting model errors.
There are several computer programs and Web servers that automate the comparative modeling process (Table 5.6.1). The accuracy of the models calculated by many of these servers is evaluated by CAMEO (Haas et al., 2013) and the biannual CASP (Critical Assessment of Techniques for Protein Structure Prediction; Moult, 2005;Moult et al., 2009) experiment.
While automation makes comparative modeling accessible to both experts and nonspecialists, manual intervention is generally still needed to maximize the accuracy of the models in the difficult cases. A number of resources useful in comparative modeling are listed in Table 5  This unit describes how to calculate comparative models using the program MODELLER (Basic Protocol). The Basic Protocol goes on to discuss all four steps of comparative modeling (Figure 5.6.1), frequently observed errors, and some applications. The Support Protocol describes how to download and install MODELLER.

MODELING LACTATE DEHYDROGENASE FROM TRICHOMONAS VAGINALIS (TvLDH) BASED ON A SINGLE TEMPLATE USING MODELLER
MODELLER is a computer program for comparative protein structure modeling (Sali and Blundell, 1993;Fiser et al., 2000). In the simplest case, the input is an alignment of a sequence to be modeled with the template structures, the atomic coordinates of the templates, and a simple script file. MODELLER then automatically calculates a model containing all non-hydrogen atoms, within minutes on a modern PC and with no user intervention. Apart from model building, MODELLER can perform additional auxiliary tasks, including fold assignment, alignment of two protein sequences or their profiles , multiple alignment of protein sequences and/or structures (Madhusudhan et al., 2006;Madhusudhan et al., 2009), calculation of phylogenetic trees, and de novo modeling of loops in protein structures (Fiser et al., 2000).

NOTE:
Further help for all the described commands and parameters may be obtained from the MODELLER Web site (see Internet Resources).

Background to TvLDH
A novel gene for lactate dehydrogenase (LDH) was identified from the genomic sequence of Trichomonas vaginalis (TvLDH). The corresponding protein had higher sequence similarity to the malate dehydrogenase of the same species (TvMDH) than to any other LDH. The authors hypothesized that TvLDH arose from TvMDH by convergent evolution relatively recently (Wu et al., 1999). Comparative models were constructed for TvLDH and TvMDH to study the sequences in a structural context and to suggest site-directed mutagenesis experiments to elucidate changes in enzymatic specificity in this apparent case of convergent evolution. The native and mutated enzymes were subsequently expressed and their activities compared (Wu et al., 1999).

Searching structures related to TvLDH Conversion of sequence to PIR file format
It is first necessary to convert the target TvLDH sequence into a format that is readable by MODELLER (file TvLDH.ali; Fig. 5.6.2). MODELLER uses the PIR format to read and write sequences and alignments. The first line of the PIR-formatted sequence consists of P1; followed by the identifier of the sequence. In this example, the sequence is identified by the code TvLDH. The second line, consisting of ten fields separated by colons, usually contains details about the structure, if any. In the case of sequences with no structural information, only two of these fields are used: the first field should be sequence (indicating that the file contains a sequence without a known structure) and the second should contain the model file name (TvLDH in this case). The rest of the file contains the sequence of TvLDH, with an asterisk (*) marking its end. The standard uppercase single-letter amino acid codes are used to represent the sequence.

Searching for suitable template structures
A search for potentially related sequences of known structure can be performed using the profile.build() command of MODELLER (file build_profile.py). The command uses the local dynamic programming algorithm to identify related sequences (Smith and Waterman, 1981). In the simplest case, the command takes as input the target sequence and a database of sequences of known structure (file pdb_95.pir) and returns a set of statistically significant alignments. The input script file for the command is shown in Figure 5.6.3.
The script, build_profile.py, does the following: 1. Initializes the "environment" for this modeling run by creating a new environ object (called env here). Almost all MODELLER scripts require this step, as the new object is needed to build most other useful objects. 2. Creates a new sequence_db object, calling it sdb, which is used to contain large databases of protein sequences.
3. Reads a file, in text format, containing nonredundant PDB sequences, into the sdb database. The sequences can be found in the file pdb_95.pir. This file is also in the PIR format. Each sequence in this file is representative of a group of PDB sequences that share 95% or more sequence identity to each other and have less than 30 residues or 30% sequence length difference.
4. Writes a binary machine-independent file containing all sequences read in the previous step.
5. Reads the binary format file back in for faster execution.
6. Creates a new "alignment" object (aln), reads the target sequence TvLDH from the file TvLDH.ali, and converts it to a profile object (prf). Profiles contain similar information to alignments, but are more compact and better for sequence database searching. iteration is run, by setting the parameter n_prof_iterations equal to 1. Thus, there is no need to check the profile for deviation (check_profile set to False). Finally, the parameter max_aln_evalue is set to 0.01, indicating that only sequences with E-values smaller than or equal to 0.01 will be included in the output.
Execute the script using the command: python build_profile.py > build_profile.log (or, if Python is not installed on your machine, with mod9.13 build_profile.py). At the end of the execution, a log file is created (build_profile.log). MODELLER always produces a log file. Errors and warnings in log files can be found by searching for the _E> and _W> strings, respectively.

Selecting a template
An extract (omitting the aligned sequences) from the file build_profile.prf is shown in Figure 5.6.4. The first six commented lines indicate the input parameters used in MODELLER to create the alignments. Subsequent lines correspond to the detected similarities by profile.build(). The most important columns in the output are the second, tenth, eleventh, and twelfth columns. The second column reports the code of the PDB sequence that was aligned to the target sequence. The eleventh column reports the percentage sequence identities between TvLDH and the PDB sequence normalized by the length of the alignment (indicated in the tenth column). In general, a sequence identity value above ß25% indicates a potential template, unless the alignment is too short (i.e., <100 residues). A better measure of the significance of the alignment is given in the twelfth column by the E-value of the alignment (lower the E-value the better). In this example, six PDB sequences show very significant similarities to the query sequence, with E-values equal to 0. As expected, all the hits correspond to malate dehydrogenases (1bdm:A, 5mdh:A, 1b8p:A, 1civ:A, 7mdh:A, and 1smk:A). To select the appropriate template for the target sequence, the alignment.compare_structures() command will first be used to assess the sequence and structure similarity between the six possible templates (file compare.py; Fig. 5.6.5).
In compare.py, the alignment object aln is created and MODELLER is instructed to read into it the protein sequences and information about their PDB files. The command malign()calculates their multiple sequence alignment, which is subsequently used as a starting point for creating a multiple structure alignment by malign3d(). Based on this structural alignment, the compare_structures() command calculates the RMS and DRMS deviations between atomic positions and distances, differences between the main-chain and side-chain dihedral angles, percentage sequence identities, and several other measures. Finally, the id_table() command writes a file (family.mat) with pairwise sequence distances that can be used as input to the dendrogram() command (or the clustering programs in the PHYLIP package; Felsenstein, 1989). dendrogram() calculates a clustering tree from the input matrix of pairwise distances, which helps visualizing differences among the template candidates. Excerpts from the log file (compare.log) are shown in Figure 5.6.6.
The objective of this step is to select the most appropriate single template structure from all the possible templates. The dendrogram in Figure 5.6.6 shows that 1civ:A and 7mdh:A are almost identical, both in terms of sequence and structure. However, 7mdh:A has a better crystallographic resolution than 1civ:A (2.4Å versus 2.8Å). From the second group of similar structures (5mdh:A, 1bdm:A, and 1b8p:A), 1bdm:A has the best resolution (1.8Å). 1smk:A is most structurally divergent among the possible templates. However, it is also the one with the lowest sequence identity (34%) to the target sequence (build_profile.prf). 1bdm:A is finally picked over 7mdh:A as the final template because of its higher overall sequence identity to the target sequence (45%).

Aligning TvLDH with the template
One way to align the sequence of TvLDH with the structure of 1bdm:A is to use the align2d() command in MODELLER (Madhusudhan et al., 2006). Although align2d() is based on a dynamic programming algorithm (Needleman and Wunsch, 1970), it is different from standard sequence-sequence alignment methods because it takes into account structural information from the template when constructing an alignment. This task is achieved through a variable gap penalty function that tends to place gaps in solvent-exposed and curved regions, outside secondary structure segments, and between two positions that are close in space. In the current example, the target-template similarity Weighted pair-group average clustering based on a distance matrix:    is so high that almost any alignment method with reasonable parameters will result in the same alignment.
The MODELLER script shown in Figure 5.6.7 aligns the TvLDH sequence in file TvLDH.ali with the 1bdm:A structure in the PDB file 1bdm.pdb (file align2d.py). In the first line of the script, an empty alignment object aln, and a new model object mdl, into which the chain A of the 1bdm structure is read, are created. append_model() transfers the PDB sequence of this model to aln and assigns it the name of 1bdmA (align_codes). The TvLDH sequence, from file TvLDH.ali, is then added to aln using append(). The align2d() command aligns the two sequences and the alignment is written out in two formats, PIR (TvLDH-1bdmA.ali) and PAP (TvLDH-1bdmA.pap). The PIR format is used by MODELLER in the subsequent model-building stage, while the PAP alignment format is easier to inspect visually. In from modeller import * from modeller.automodel imort * #from modeller import soap_protein_od env ϭ environ() a ϭ automodel(env, alnfileϭ TvLDH-1bdmA.ali , knownsϭ 1bdmA , sequenceϭ TvLDH , assess_methodsϭ(assess.DOPE, #soap_protein_od.Scorer(), assess.GA341)) a.starting_model ϭ 1 a.ending_model ϭ 5 a.make()

Model building
Once a target-template alignment is constructed, MODELLER calculates a 3-D model of the target completely automatically, using its automodel class. The script in Figure 5.6.9 will generate five different models of TvLDH based on the 1bdm:A template structure and the alignment in file TvLDH-1bdmA.ali (file model-single.py).
The first line ( Fig. 5.6.9) loads the automodel class and prepares it for use. An automodel object is then created and called "a" and parameters are set to guide the model-building procedure. alnfile names the file that contains the target-template alignment in the PIR format. knowns defines the known template structure(s) in alnfile (TvLDH-1bdmA.ali) and sequence defines the code of the target sequence. starting_model and ending_model define the number of models that are calculated (their indices will run from 1 to 5). The last line in the file calls the make method that actually calculates the models. The most important output files are model-single.log, which reports warnings, errors and other useful information including the input restraints used for modeling that remain violated in the final model, and TvLDH.B9999000 [1][2][3][4][5].pdb, which contain the coordinates of the five produced models, in the PDB format. The models can be viewed by any program that reads the PDB format, such as Chimera (http://www.cgl.ucsf.edu/chimera/) or RasMol (http://www.rasmol.org).

Evaluating a model
If several models are calculated for the same target, the best model can be selected by picking the model with the lowest value of the MODELLER objective function or the DOPE (Shen and Sali, 2006) or SOAP (Dong et al., 2013) assessment scores, which are reported at the end of the log file. (To calculate the SOAP score, download the SOAP-Protein library file from http://salilab.org/SOAP/ and uncomment the two SOAPrelated lines in model-single.py by removing the '#' characters.) In this example, the second model (TvLDH.B99990002.pdb) has the lowest objective function and is selected. None of these scores are absolute measures, meaning that they can only be used to rank models calculated from the same alignment.
Once a final model is selected, there are many ways to further assess it. In this example, the DOPE potential in MODELLER is used to evaluate the fold of the selected model. Links to other programs for model assessment can be found in Table 5.6.1. However, before any external evaluation of the model, one should check the log file from the modeling run for runtime errors (model-single.log) and restraint violations (see the MODELLER manual for details).
The script, evaluate_model.py ( Fig. 5.6.10) evaluates the model with the DOPE potential. In this script, the atomic coordinates of the PDB file are read in (using com-plete_pdb()) to a model object, mdl. This is necessary for MODELLER to correctly calculate the energy, and additionally allows for the possibility of the PDB file having atoms in a nonstandard order, or having different subsets of atoms (e.g., all atoms including hydrogens, while MODELLER uses only heavy atoms, or vice versa). The DOPE energy is then calculated using assess_dope(). An energy profile is additionally requested, smoothed over a 15-residue window, and normalized by the number of restraints acting on each residue. This profile is written to a file TvLDH.profile, which can be used as input to a graphing program such as GNUPLOT.
Similarly, the profile can be calculated for the template structure (see the scripts eval-uate_template.py and plot_profiles.py in the zipfile). A comparison of the two profiles is shown in Figure 5.6.11. It can be seen that the DOPE score profile shows clear differences between the two profiles for the long active-site loop between residues 90 and 100 and the long helices at the C-terminal end of the target sequence. This long loop interacts with region 220 to 250, which forms the other half of the active site. This latter region is well resolved in both the template and the target structure. However, probably due to the unfavorable nonbonded interactions with the 90 to 100 region, it is reported to be of high energy by DOPE. It is to be noted that a region of high energy indicated by DOPE may not always necessarily indicate actual error, especially when it highlights an active site or a protein-protein interface. However, in this case, the same active-site loops have a better profile in the template structure, which strengthens the argument that the model is probably incorrect in the active-site region.
Resolution of such problems is beyond the scope of this unit, but is described in a more advanced modeling tutorial available at http://salilab.org/modeller/tutorial/advanced .html.

OBTAINING AND INSTALLING MODELLER
MODELLER is written in Fortran 90 and uses Python for its control language. All input scripts to MODELLER are, hence, Python scripts. While knowledge of Python is not necessary to run MODELLER, it can be useful in performing more advanced tasks. Precompiled binaries for MODELLER can be downloaded from http://salilab.org/ modeller. gunzip modeller-9.13.tar.gz tar -xvf modeller-9.13.tar 6. The files needed for the installation can be found in a newly created directory called modeller-9.13. Move into that directory and start the installation with the following commands:

Necessary Resources
cd modeller-9.13 ./Install 7. The installation script will prompt the user with several questions and suggest default answers. To accept the default answers, press the Enter key. The various prompts are briefly discussed below: a. For the prompt below, choose the appropriate combination of the machine architecture and operating system. For this example, choose the default answer by pressing the Enter key.
3) x86_64 (Opteron/EM64T) box (Linux). 4) Alternative Linux x86 PC binary (e.g., for FreeBSD). Select the type of your computer from the list above [1]: b. For the prompt below, tell the installer where to install the MODELLER executables. The default choice will place it in the directory indicated, but any directory to which the user has write permissions may be specified. Full directory name for the installed MODELLER9.13 [<YOUR-HOME-DIRECTORY>/bin/modeller9.13]: c. For the prompt below, enter the MODELLER license key obtained in step 3.
KEY_MODELLER9v13, obtained from our academic license server at http://salilab.org/modeller/registration .html: 8. The installer will now confirm the answers to the above prompts. Press Enter to begin the installation. The mod9.13 script installed in the chosen directory can now be used to invoke MODELLER. The installer will also provide information on how to set up MODELLER to work with your operating system's copy of Python.

Background Information
As stated earlier, comparative modeling consists of four main steps: fold assignment, target-template alignment, model building and model evaluation (Marti-Renom et al., 2000;Fig. 5.6.1).

Fold assignment and target-template alignment
Although fold assignment and sequencestructure alignment are logically two distinct steps in the process of comparative modeling, in practice, almost all fold-assignment methods also provide sequence-structure alignments. In the past, fold-assignment methods were optimized for better sensitivity in detecting remotely related homologs, often at the cost of alignment accuracy. However, recent methods simultaneously optimize both the sensitivity and alignment accuracy. Therefore, in the following discussion, fold assignment and sequence-structure alignment will be treated as a single procedure, explaining the differences as needed.

Fold assignment
The primary requirement for comparative modeling is the identification of one or more known template structures with detectable similarity to the target sequence. The identification of suitable templates is achieved by scanning structure databases, such as PDB (Berman et al., 2000), SCOP (Andreeva et al., 2004), DALI, UNIT 5.5 (Dietmann et al., 2001), and CATH (Pearl et al., 2005), with the target sequence as the query. The detected similarity is usually quantified in terms of sequence identity or statistical measures such as E-value or z-score, depending on the method used.

Three regimes of the sequence-structure relationship
The sequence-structure relationship can be subdivided into three different regimes in the sequence similarity spectrum: (i) the easily detected relationships, characterized by >30% sequence identity; (ii) the "twilight zone" (Rost, 1999), corresponding to relationships with statistically significant sequence similarity, with identities in the 10% to 30% range; and (iii) the "midnight zone" (Rost, 1999), corresponding to statistically insignificant sequence similarity.

Pairwise sequence alignment methods
For closely related protein sequences with identities higher than 30% to 40%, the alignments produced by all methods are almost always largely correct. The quickest way to search for suitable templates in this regime is to use simple pairwise sequence alignment methods such as SSEARCH (Pearson, 1994), BLAST (Altschul et al., 1997), and FASTA (Pearson, 1994). Brenner et al. (1998) showed that these methods detect only ß18% of the Modeling Structure from Sequence homologous pairs at less than 40% sequence identity, while they identify more than 90% of the relationships when sequence identity is between 30% and 40% (Brenner et al., 1998). Another benchmark, based on 200 reference structural alignments with 0% to 40% sequence identity, indicated that BLAST is able to correctly align only 26% of the residue positions (Sauder et al., 2000).

Profile-sequence alignment methods
The sensitivity of the search and accuracy of the alignment become progressively difficult as the relationships move into the twilight zone (Saqi et al., 1998;Rost, 1999). A significant improvement in this area was the introduction of profile methods by Gribskov et al. (1987). The profile of a sequence is derived from a multiple sequence alignment and specifies residue-type occurrences for each alignment position. The information in a multiple sequence alignment is most often encoded as either a position-specific scoring matrix (PSSM; Henikoff and Henikoff, 1994;Altschul et al., 1997) or as a Hidden Markov Model (HMM; Krogh et al., 1994;Eddy, 1998). In order to identify suitable templates for comparative modeling, the profile of the target sequence is used to search against a database of template sequences. The profilesequence methods are more sensitive in detecting related structures in the twilight zone than the pairwise sequence-based methods; they detect approximately twice the number of homologs under 40% sequence identity (Park et al., 1998;Lindahl and Elofsson, 2000;Sauder et al., 2000). The resulting profilesequence alignments correctly align approximately 43% to 48% of residues in the 0% to 40% sequence identity range (Sauder et al., 2000;Marti-Renom et al., 2004); this number is almost twice as large as that of the pairwise sequence methods. Frequently used programs for profile-sequence alignment are PSI-BLAST (Altschul et al., 1997), SAM , HMMER (Eddy, 1998), HHsearch (Soding, 2005), HHBlits (Remmert et al., 2012), and BUILD_PROFILE (part of MODELLER; Sali and Blundell, 1993).

Profile-profile alignment methods
As a natural extension, the profile-sequence alignment methods have led to profile-profile alignment methods that search for suitable template structures by scanning the profile of the target sequence against a database of template profiles as opposed to a database of template sequences. These methods have proven to include the most sensitive and accurate fold assignment and alignment protocols to date (Edgar and Sjolander, 2004;Marti-Renom et al., 2004;Ohlson et al., 2004;Wang and Dunbrack, 2004). Profile-profile methods detect ß28% more relationships at the superfamily level and improve the alignment accuracy by 15% to 20%, compared to profile-sequence methods Zhou and Zhou, 2005). There are a number of variants of profile-profile alignment methods that differ in the scoring functions they use (Pietrokovski, 1996;Rychlewski et al., 1998;Yona and Levitt, 2002;Panchenko, 2003;Sadreyev and Grishin, 2003;von Ohsen et al., 2003;Edgar and Sjolander, 2004;Marti-Renom et al., 2004;Zhou and Zhou, 2005). However, several analyses have shown that the overall performances of these methods are comparable (Edgar and Sjolander, 2004;Marti-Renom et al., 2004;Ohlson et al., 2004;Wang and Dunbrack, 2004). Some of the programs that can be used to detect suitable templates are FFAS (Jaroszewski et al., 2005), SP3 (Zhou and Zhou, 2005), SALIGN , HHBlits (Remmert et al., 2012), HHsearch (Soding, 2005), and PPSCAN (part of MODELLER; Sali and Blundell, 1993).

Sequence-structure threading methods
As the sequence identity drops below the threshold of the twilight zone, there is usually insufficient signal in the sequences or their profiles for the sequence-based methods discussed above to detect true relationships (Lindahl and Elofsson, 2000). Sequencestructure threading methods are most useful in this regime, as they can sometimes recognize common folds even in the absence of any statistically significant sequence similarity (Godzik, 2003). These methods achieve higher sensitivity by using structural information derived from the templates. The accuracy of a sequence-structure match is assessed by the score of a corresponding coarse model and not by sequence similarity, as in sequencecomparison methods (Godzik, 2003). The scoring scheme used to evaluate the accuracy is either based on residue substitution tables dependent on structural features such as solvent exposure, secondary structure type, and hydrogen-bonding properties (Shi et al., 2001;Karchin et al., 2003;McGuffin and Jones, 2003;Zhou and Zhou, 2005), or on statistical potentials for residue interactions implied by the alignment (Sippl, 1990(Sippl, , 1995Bowie et al., 1991;Skolnick and Kihara, 2001 al., 2003). The use of structural data does not have to be restricted to the structure side of the aligned sequence-structure pair. For example, SAM-T08 makes use of the predicted local structure for the target sequence to enhance homolog detection and alignment accuracy . Commonly used threading programs are GenTHREADER (Jones, 1999;McGuffin and Jones, 2003), 3D-PSSM (Kelley et al., 2000), FUGUE (Shi et al., 2001), SP3 (Zhou and Zhou, 2005), SAM-T08 multi-track HMM Karplus et al., 2003), and MUSTER (Wu and Zhang, 2008).

Iterative sequence-structure alignment and model building
Yet another strategy is to optimize the alignment by iterating over the process of calculating alignments, building models, and evaluating models. Such a protocol can sample alignments that are not statistically significant and identify the alignment that yields the best model. Although this procedure can be time consuming, it can significantly improve the accuracy of the resulting comparative models in difficult cases (John and .

Importance of an accurate alignment
Regardless of the method used, searching in the twilight and midnight zones of the sequence-structure relationship often results in false negatives, false positives, or alignments that contain an increasingly large number of gaps and alignment errors. Improving the performance and accuracy of methods in this regime remains one of the main tasks of comparative modeling today (Moult, 2005). It is imperative to calculate an accurate alignment between the target-template pair, as comparative modeling can almost never recover from an alignment error (Sanchez and Sali, 1997a).

Template selection
After a list of all related protein structures and their alignments with the target sequence have been obtained, template structures are prioritized depending on the purpose of the comparative model. Template structures may be chosen based purely on the target-template sequence identity, or on a combination of several other criteria, such as experimental accuracy of the structures (resolution of X-ray structures, number of restraints per residue for NMR structures), conservation of activesite residues, holo-structures that have bound ligands of interest, and prior biological information that pertains to the solvent, pH, and quaternary contacts. It is not necessary to select only one template. In fact, the use of several templates approximately equidistant from the target sequence generally increases the model accuracy (Srinivasan and Blundell, 1993;Sanchez and Sali, 1997b).

Modeling by assembly of rigid bodies
The first and still widely used approach in comparative modeling is to assemble a model from a small number of rigid bodies obtained from the aligned protein structures (Browne et al., 1969;Greer, 1981;Blundell et al., 1987). The approach is based on the natural dissection of the protein structures into conserved core regions, variable loops that connect them, and side chains that decorate the backbone. For example, the following semiautomated procedure is implemented in the computer program COMPOSER (Sutcliffe et al., 1987a). First, the template structures are selected and superposed. Second, the "framework" is calculated by averaging the coordinates of the C α atoms of structurally conserved regions in the template structures. Third, the main-chain atoms of each core region in the target model are obtained by superposing the core segment, from the template whose sequence is closest to the target, on the framework. Fourth, the loops are generated by scanning a database of all known protein structures to identify the structurally variable regions that fit the anchor core regions and have a compatible sequence (Topham et al., 1993). Fifth, the side chains are modeled based on their intrinsic conformational preferences and on the conformation of the equivalent side chains in the template structures (Sutcliffe et al., 1987b). Finally, the stereochemistry of the model is improved either by a restrained energy minimization or a molecular dynamics refinement. The accuracy of a model can be somewhat increased when more than one template structure is used to construct the framework and when the templates are averaged into the framework using weights corresponding to their sequence similarities to the target sequence (Srinivasan and Blundell, 1993). Possible future improvements of modeling by rigid-body assembly include incorporation of rigid body shifts, such as the relative shifts in the packing of a helices and β-sheets (Nagarajaram et al., 1999). Three other programs that implement this method are 3D-JIGSAW (Bates et al., 2001), RosettaCM (Song et al., 2013), and SWISS-MODEL (Schwede et al., 2003). Modeling by segment matching or coordinate reconstruction The basis of modeling by coordinate reconstruction is the finding that most hexapeptide segments of protein structure can be clustered into only 100 structurally different classes (Jones and Thirup, 1986;Claessens et al., 1989;Unger et al., 1989;Levitt, 1992;Bystroff and Baker, 1998). Thus, comparative models can be constructed by using a subset of atomic positions from template structures as guiding positions to identify and assemble short, allatom segments that fit these guiding positions. The guiding positions usually correspond to the C α atoms of the segments that are conserved in the alignment between the template structure and the target sequence. The all-atom segments that fit the guiding positions can be obtained either by scanning all known protein structures, including those that are not related to the sequence being modeled (Claessens et al., 1989;Holm and Sander, 1991), or by a conformational search restrained by an energy function (Bruccoleri and Karplus, 1987;van Gelder et al., 1994). This method can construct both main-chain and side-chain atoms, and can also model unaligned regions (gaps). It is implemented in the program SegMod (Levitt, 1992). Even some side-chain modeling methods (Chinea et al., 1995) and the class of loop-construction methods based on finding suitable fragments in the database of known structures (Jones and Thirup, 1986) can be seen as segment-matching or coordinatereconstruction methods.

Modeling by satisfaction of spatial restraints
The methods in this class begin by generating many constraints or restraints on the structure of the target sequence, using its alignment to related protein structures as a guide. The procedure is conceptually similar to that used in determination of protein structures from NMR-derived restraints. The restraints are generally obtained by assuming that the corresponding distances between aligned residues in the template and the target structures are similar. These homology-derived restraints are usually supplemented by stereochemical restraints on bond lengths, bond angles, dihedral angles, and nonbonded atom-atom contacts that are obtained from a molecular mechanics force field. The model is then derived by minimizing the violations of all the restraints. This optimization can be achieved either by distance geometry or real-space optimization. For example, an elegant distance geometry approach constructs all-atom models from lower and upper bounds on distances and dihedral angles (Havel and Snow, 1991).
Comparative protein structure modeling by MODELLER. MODELLER, the authors' own program for comparative modeling, belongs to this group of methods (Sali and Blundell, 1993;Sali and Overington, 1994;Fiser et al., 2000;. MODELLER implements comparative protein structure modeling by satisfaction of spatial restraints. The program was designed to use as many different types of information about the target sequence as possible. Homology-derived restraints. In the first step of model building, distance and dihedral angle restraints on the target sequence are derived from its alignment with template 3-D structures. The form of these restraints was obtained from a statistical analysis of the relationships between similar protein structures. The analysis relied on a database of 105 family alignments that included 416 proteins of known 3-D structure (Sali and Overington, 1994). By scanning the database of alignments, tables quantifying various correlations were obtained, such as the correlations between two equivalent C α -C α distances, or between equivalent main-chain dihedral angles from two related proteins (Sali and Blundell, 1993). These relationships are expressed as conditional probability density functions (pdf's), and can be used directly as spatial restraints. For example, probabilities for different values of the main-chain dihedral angles are calculated from the type of residue considered, from main-chain conformation of an equivalent residue, and from sequence similarity between the two proteins. Another example is the pdf for a certain C α -C α distance given equivalent distances in two related protein structures. An important feature of the method is that the form of spatial restraints was obtained empirically, from a database of protein structure alignments.
Stereochemical restraints. In the second step, the spatial restraints and the CHARMM22 force-field terms enforcing proper stereochemistry (MacKerell et al., 1998) are combined into an objective function. The general form of the objective function is similar to that in molecular dynamics programs, such as CHARMM22 (MacKerell et al., 1998). The objective function depends on the Cartesian coordinates of ß10,000 atoms (3-D points) that form the modeled molecules. For a 10,000-atom system, there can be on the order of 200,000 restraints. The functional form of each term is simple; it includes a quadratic function, harmonic lower and upper bounds, cosine, a weighted sum of a few Gaussian functions, Coulomb's law, Lennard-Jones potential, and cubic splines. The geometric features presently include a distance, an angle, a dihedral angle, a pair of dihedral angles between two, three, four, and eight atoms, respectively, the shortest distance in the set of distances, solvent accessibility, and atom density that is expressed as the number of atoms around the central atom. Some restraints can be used to restrain pseudo-atoms, e.g., the gravity center of several atoms.
Optimization of the objective function. Finally, the model is obtained by optimizing the objective function in Cartesian space. The optimization is carried out by the use of the variable target function method (Braun and Go, 1985), employing methods of conjugate gradients and molecular dynamics with simulated annealing (Clore et al., 1986). Several slightly different models can be calculated by varying the initial structure, and the variability among these models can be used to estimate the lower bound on the errors in the corresponding regions of the fold.
Restraints derived from experimental data. Because the modeling by satisfaction of spatial restraints can use many different types of information about the target sequence, it is perhaps the most promising of all comparative modeling techniques. One of the strengths of modeling by satisfaction of spatial restraints is that restraints derived from a number of different sources can easily be added to the homology-derived restraints. For example, restraints could be provided by rules for secondary-structure packing (Cohen et al., 1989), analyses of hydrophobicity (Aszodi and Taylor, 1994) and correlated mutations , empirical potentials of mean force (Sippl, 1990), nuclear magnetic resonance (NMR) experiments (Sutcliffe et al., 1992), cross-linking experiments, fluorescence spectroscopy, image reconstruction in electron microscopy, site-directed mutagenesis (Boissel et al., 1993), and intuition, among other sources. Especially in difficult cases, a comparative model could be improved by making it consistent with available experimental data and/or with more general knowledge about protein structure.
Relative accuracy, flexibility, and automation. Accuracies of the various model-building methods are relatively similar when used optimally (Marti-Renom et al., 2002). Other factors such as template selection and align-ment accuracy usually have a larger impact on the model accuracy, especially for models based on low sequence identity to the templates. However, it is important that a modeling method allow a degree of flexibility and automation to obtain better models more easily and rapidly. For example, a method should allow for an easy recalculation of a model when a change is made in the alignment. It should also be straightforward enough to calculate models based on several templates, and should provide tools for incorporation of prior knowledge about the target (e.g., cross-linking restraints, predicted secondary structure) and allow ab initio modeling of insertions (e.g., loops), which can be crucial for annotation of function.

Loop modeling
Loop modeling is an especially important aspect of comparative modeling in the range of 30% to 50% sequence identity. In this range of overall similarity, loops among the homologs vary while the core regions are still relatively conserved and aligned accurately. Loops often play an important role in defining the functional specificity of a given protein, forming the active and binding sites. Loop modeling can be seen as a mini protein folding problem, because the correct conformation of a given segment of a polypeptide chain has to be calculated mainly from the sequence of the segment itself. However, loops are generally too short to provide sufficient information about their local fold. Even identical decapeptides in different proteins do not always have the same conformation (Kabsch and Sander, 1984;Mezei, 1998). Some additional restraints are provided by the core anchor regions that span the loop and by the structure of the rest of the protein that cradles the loop. Although many loopmodeling methods have been described, it is still challenging to correctly and confidently model loops longer than ß10 to 12 residues (Fiser et al., 2000;Jacobson et al., 2004;Zhu et al., 2006).
There are two main classes of loopmodeling methods: (i) database search approaches that scan a database of all known protein structures to find segments fitting the anchor core regions (Jones and Thirup, 1986;Chothia and Lesk, 1987); (ii) conformational search approaches that rely on optimizing a scoring function (Moult and James, 1986;Bruccoleri and Karplus, 1987;Shenkin et al., 1987). There are also methods that combine these two approaches (van Vlijmen and Karplus, 1997;Deane and Blundell, 2001 Loop modeling by database search. The database search approach to loop modeling is accurate and efficient when a database of specific loops is created to address the modeling of the same class of loops, such as βhairpins (Sibanda et al., 1989), or loops on a specific fold, such as the hypervariable regions in the immunoglobulin fold (Chothia and Lesk, 1987;Chothia et al., 1989). There are attempts to classify loop conformations into more general categories, thus extending the applicability of the database search approach (Ring et al., 1992;Oliva et al., 1997;Rufino et al., 1997;Fernandez-Fuentes and Fiser, 2006). However, the database methods are limited because the number of possible conformations increases exponentially with the length of a loop, and until the late 1990s only loops up to 7 residues long could be modeled using the database of known protein structures (Fidelis et al., 1994;Lessel and Schomburg, 1994). However, the growth of the PDB in recent years has largely eliminated this problem (Fernandez-Fuentes and Fiser, 2006).
Loop modeling by conformational search. There are many such methods, exploiting different protein representations, objective functions, and optimization or enumeration algorithms. The search algorithms include the minimum perturbation method (Fine et al., 1986), dihedral angle search through a rotamer library (Zhu et al., 2006;Sellers et al., 2008), molecular dynamics simulations (Bruccoleri and Karplus, 1990;van Vlijmen and Karplus, 1997), genetic algorithms (Ring et al., 1993), Monte Carlo and simulated annealing (Higo et al., 1992;Collura et al., 1993;Abagyan and Totrov, 1994), multiple copy simultaneous search (Zheng et al., 1993), self-consistent field optimization (Koehl and Delarue, 1995), robotics-inspired kinematic closure (Mandell et al., 2009), and enumeration based on graph theory (Samudrala and Moult, 1998). The accuracy of loop predictions can be further improved by clustering the sampled loop conformations and partially accounting for the entropic contribution to the free energy (Xiang et al., 2002). Another way to improve the accuracy of loop predictions is to consider the solvent effects. Improvements in implicit solvation models, such as the Generalized Born solvation model, motivated their use in loop modeling. The solvent contribution to the free energy can be added to the scoring function for optimization, or it can be used to rank the sampled loop conformations after they are generated with a scoring function that does not include the solvent terms (Fiser et al., 2000;Felts et al., 2002;DePristo et al., 2003).
Loop modeling in MODELLER. The loopmodeling module in MODELLER implements the optimization-based approach (Fiser et al., 2000;Fiser and Sali, 2003b). The main reasons for choosing this implementation are the generality and conceptual simplicity of scoring function minimization. Loop prediction by optimization is applicable to simultaneous modeling of several loops and loops interacting with ligands, which is not straightforward with the database-search approaches. Loop optimization in MODELLER relies on conjugate gradients and molecular dynamics with simulated annealing. The pseudo energy function is a sum of many terms, including some terms from the CHARMM22 molecular mechanics force field (MacKerell et al., 1998) and spatial restraints based on distributions of distances (Sippl, 1990;Melo et al., 2002) and dihedral angles in known protein structures. The method was tested on a large number of loops of known structure, both in the native and nearnative environments (Fiser et al., 2000).

Comparative model building by iterative alignment, model building, and model assessment
Comparative or homology protein structure modeling is severely limited by errors in the alignment of a modeled sequence with related proteins of known three-dimensional structure. To ameliorate this problem, one can use an iterative method that optimizes both the alignment and the model implied by it (Sanchez and Sali, 1997a;Miwa et al., 1999). This task can be achieved by a genetic algorithm protocol that starts with a set of initial alignments and then iterates through realignment, model building, and model assessment to optimize a model assessment score (John and Sali, 2003). During this iterative process: (1) new alignments are constructed by the application of a number of genetic algorithm operators, such as alignment mutations and crossovers; (2) comparative models corresponding to these alignments are built by satisfaction of spatial restraints, as implemented in the program MODELLER; and (3) the models are assessed by a composite score, partly depending on an atomic statistical potential (Melo et al., 2002). When testing the procedure on a very difficult set of 19 modeling targets sharing only 4% to 27% sequence identity with their template structures, the average final alignment accuracy increased from 37% to 45% relative to the initial alignment (the alignment accuracy was measured as the percentage of positions in the tested alignment that were identical to the reference structurebased alignment). Correspondingly, the average model accuracy increased from 43% to 54% (the model accuracy was measured as the percentage of the C α atoms of the model that were within 5Å of the corresponding C α atoms in the superimposed native structure).

Errors in comparative models
As the similarity between the target and the templates decreases, the errors in the model increase. Errors in comparative models can be divided into five categories (Sanchez and Sali, 1997a,b;Fig. 5.6.12), as follows: Errors in side-chain packing (Fig. 5.6.12A). As the sequences diverge, the packing of side chains in the protein core changes. Sometimes even the conformation of identical side chains is not conserved, a pitfall for many comparative modeling methods. Side-chain errors are critical if they occur in regions that are involved in protein function, such as active sites and ligand-binding sites.
Distortions and shifts in correctly aligned regions (Fig. 5.6.12B). As a consequence of sequence divergence, the main-chain conformation changes, even if the overall fold remains the same. Therefore, it is possible that in some correctly aligned segments of a model the template is locally different (>3Å) from the target, resulting in errors in that region. The structural differences are sometimes not due to differences in sequence, but are a consequence of artifacts in structure determination or structure determination in different environments (e.g., packing of subunits in a crystal). The simultaneous use of several templates can minimize this kind of error (Srinivasan and Blundell, 1993;Sanchez and Sali, 1997a,b).
Errors in regions without a template (Fig. 5.6.12C). Segments of the target sequence that have no equivalent region in the template structure (i.e., insertions or loops) are the most difficult regions to model. If the insertion is relatively short, <9 residues long, some methods can correctly predict the conformation of the backbone (van Vlijmen and Karplus, 1997;Fiser et al., 2000;Jacobson et al., 2004). Conditions for successful prediction are the correct alignment and an accurately modeled environment surrounding the insertion.
Errors due to misalignments (Fig. 5.6.12D). The largest single source of errors in comparative modeling is misalignments, especially when the target-template sequence identity decreases below 30%. However, alignment errors can be minimized in two ways. First, it is usually possible to use a large number of sequences to construct a multiple alignment, even if most of these sequences do not have known structures. Multiple alignments are generally more reliable than pairwise alignments (Barton and Sternberg, 1987;Taylor et al., 1994). The second way of improving the alignment is to iteratively modify those regions in the alignment that correspond to predicted errors in the model (Sanchez and Sali, 1997a,b;John and Sali, 2003).
Incorrect templates (Fig. 5.6.12E). This is a potential problem when distantly related proteins are used as templates (i.e., <25% sequence identity). Distinguishing between a model based on an incorrect template and a model based on an incorrect alignment with a correct template is difficult. In both cases, the evaluation methods will predict an unreliable model. The conservation of the key functional or structural residues in the target sequence increases the confidence in a given fold assignment.

Predicting the model accuracy
The accuracy of the predicted model determines the information that can be extracted from it. Thus, estimating the accuracy of a model in the absence of the known structure is essential for interpreting it.
Initial assessment of the fold. As discussed earlier, a model calculated using a template structure that shares more than 30% sequence identity is indicative of an overall accurate structure. However, when the sequence identity is lower, the first aspect of model evaluation is to confirm whether or not a correct template was used for modeling. It is often the case, when operating in this regime, that the fold-assignment step produces only false positives. A further complication is that at such low similarities the alignment generally contains many errors, making it difficult to distinguish between an incorrect template on one hand and an incorrect alignment with a correct template on the other hand. There are several methods that use 3-D profiles and statistical potentials (Sippl, 1990;Luthy et al., 1992;Melo et al., 2002) to assess the compatibility between the sequence and modeled structure by evaluating the environment of each residue in a model with respect to the expected environment as found in native high-resolution experimental structures. These methods can be used to assess whether or not the correct is compared with its model (green) and with the template fatty acid binding protein (blue). (C) Errors in regions without a template. The C α trace of the 112-117 loop is shown for the X-ray structure of human eosinophil neurotoxin (red), its model (green), and the template ribonuclease A structure (residues 111-117; blue). (D) Errors due to misalignments. The N-terminal region in the crystal structure of human eosinophil neurotoxin (red) is compared with its model (green). The corresponding region of the alignment with the template ribonuclease A is shown. The red lines show correct equivalences, that is, residues whose C α atoms are within 5Å of each other in the optimal leastsquares superposition of the two X-ray structures. The "a" characters in the bottom line indicate helical residues and "b" characters, the residues in sheets. (E) Errors due to an incorrect template. The X-ray structure of α-trichosanthin (red) is compared with its model (green) that was calculated using indole-3-glycerophosphate synthase as the template.
template was used for the modeling. They include VERIFY3D (Luthy et al., 1992), Prosa2003 (Sippl, 1993;Wiederstein and Sippl, 2007), HARMONY (Topham et al., 1994), ANOLEA (Melo and Feytmans, 1998), DFIRE (Zhou and Zhou, 2002), DOPE (Shen and Sali, 2006), SOAP (Dong et al., 2013), QMEAN local (Benkert et al., 2011), and TSV-Mod (Eramian et al., 2008). Even when the model is based on alignments that have >30% sequence identity, other factors, including the environment, can strongly influence the accuracy of a model. For instance, some calcium-binding proteins undergo large conformational changes when bound to calcium. If a calcium-free template is used to model the calcium-bound state of the target, it is likely that the model will be incorrect irrespective of the target-template similarity or accuracy of the template structure (Pawlowski et al., 1996).
Evaluations of self-consistency. The model should also be subjected to evaluations of self-consistency to ensure that it satisfies the restraints used to calculate it. Additionally, the stereochemistry of the model (e.g., bondlengths, bond-angles, backbone torsion angles, and nonbonded contacts) may be evaluated using programs such as PROCHECK (Laskowski et al., 1993) and WHATCHECK (Hooft et al., 1996). Although errors in stereochemistry are rare and less informative than errors detected by statistical potentials, a cluster of stereochemical errors may indicate that there are larger errors (e.g., alignment errors) in that region.

Applications
Comparative modeling is often an efficient way to obtain useful information about the protein of interest. For example, comparative models can be helpful in designing mutants to test hypotheses about the protein's function (Wu et al., 1999;Vernal et al., 2002); in identifying active and binding sites (Sheng et al., 1996); in searching for, designing, and improving ligand binding strength for a given binding site (Ring et al., 1993;Li et al., 1996;Selzer et al., 1997;Enyedy et al., 2001;Que et al., 2002); in modeling substrate specificity (Xu et al., 1996); in predicting antigenic epitopes (Sali and Blundell, 1993); in simulating protein-protein docking (Vakser, 1995); in inferring function from calculated electrostatic potential around the protein (Matsumoto et al., 1995); in facilitating molecular replacement in X-ray structure determination (Howell et al., 1992); in refining models based on NMR constraints (Modi et al., 1996); in testing and improving a sequence-structure alignment (Wolf et al., 1998); in annotating single nucleotide polymorphisms (Mirkovic et al., 2004;Karchin et al., 2005); in structural characterization of large complexes by docking to low-resolution cryo-electron density maps (Spahn et al., 2001;Gao et al., 2003); and in rationalizing known experimental observations. Fortunately, a 3-D model does not have to be absolutely perfect to be helpful in biology, as demonstrated by the applications listed above. The type of a question that can be addressed with a particular model does depend on its accuracy (Fig. 5.6.13).
At the low end of the accuracy spectrum, there are models that are based on less than 25% sequence identity, and that sometimes have less than 50% of their C α atoms within 3.5Å of their correct positions. However, such models still have the correct fold, and even knowing only the fold of a protein may sometimes be sufficient to predict its approximate biochemical function. Models in this low range of accuracy, combined with model evaluation, can be used for confirming or rejecting a match between remotely related proteins (Sanchez and Sali, 1997a;. In the middle of the accuracy spectrum are the models based on approximately 35% sequence identity, corresponding to 85% of the C α atoms modeled within 3.5Å of their correct positions. Fortunately, the active and binding sites are frequently more conserved than the rest of the fold, and are thus modeled more accurately (Sanchez and Sali, 1998). In general, medium-resolution models frequently allow a refinement of the functional prediction based on sequence alone, because ligand binding is most directly determined by the structure of the binding site rather than its sequence. It is frequently possible to correctly predict important features of the target protein that do not occur in the template structure. For example, the location of a binding site can be predicted from clusters of charged residues (Matsumoto et al., 1995), and the size of a ligand may be predicted from the volume of the binding-site cleft (Xu et al., 1996). Medium-resolution models can also be used to construct site-directed mutants with altered or destroyed binding capacity, which in turn could test hypotheses about the sequence-structure-function relationships. Other problems that can be addressed with medium-resolution comparative models include designing proteins that have compact structures, without long tails, loops, and exposed hydrophobic residues, for better crystallization, or designing proteins with added disulfide bonds for extra stability.
The high end of the accuracy spectrum corresponds to models based on 50% sequence identity or more. The average accuracy of these models approaches that of low-resolution X-ray structures (3Å resolution) or medium-resolution NMR structures (10 distance restraints per residue; Sanchez and Sali, 1997b). The alignments on which these models are based generally contain almost no errors. Models with such high accuracy have been shown to be useful even for refining crystallographic structures by the method of molecular replacement (Howell et al., 1992;Baker and Sali, 2001;Jones, 2001;Claude et al., 2004;Schwarzenbacher et al., 2004).

Conclusion
Over the past few years, there has been a gradual increase in both the accuracy of comparative models and the fraction of protein sequences that can be modeled with useful accuracy (Marti-Renom et al., 2000;Baker and Sali, 2001;Pieper et al., 2011). The magnitude of errors in fold assignment, alignment, and modeling of side-chains and loops have decreased considerably. These improvements are a consequence both of better techniques and a larger number of known protein sequences and structures. Nevertheless, all the errors remain  Figure 5.6.13 Accuracy and application of protein structure models. The vertical axis indicates the different ranges of applicability of comparative protein structure modeling, the corresponding accuracy of protein structure models, and their sample applications. (A) The docosahexaenoic fatty acid ligand (violet) was docked into a high-accuracy comparative model of brain lipid-binding protein (right), modeled based on its 62% sequence identity to the crystallographic structure of adipocyte lipid-binding protein (PDB code 1adl). A number of fatty acids were ranked for their affinity to brain lipid-binding protein consistently with site-directed mutagenesis and affinity chromatography experiments (Xu et al., 1996), even though the ligand specificity profile of this protein is different from that of the template structure. Typical overall accuracy of a comparative model in this range of sequence similarity is indicated by a comparison of a model for adipocyte fatty acid binding protein with its actual structure (left). (B) A putative proteoglycan binding patch was identified on a medium-accuracy comparative model of mouse mast cell protease 7 (right), modeled based on its 39% sequence identity to the crystallographic structure of bovine pancreatic trypsin (2ptn) that does not bind proteoglycans. The prediction was confirmed by site-directed mutagenesis and heparin-affinity chromatography experiments (Matsumoto et al., 1995). Typical accuracy of a comparative model in this range of sequence similarity is indicated by a comparison of a trypsin model with the actual structure. (C) A molecular model of the whole yeast ribosome (right) was calculated by fitting atomic rRNA and protein models into the electron density of the 80S ribosomal particle, obtained by electron microscopy at 15Å resolution (Spahn et al., 2001). Most of the models for 40 out of the 75 ribosomal proteins were based on template structures that were approximately 30% sequence-identical. Typical accuracy of a comparative model in this range of sequence similarity is indicated by a comparison of a model for a domain in L2 protein from B. stearothermophilus with the actual structure (1rl2). significant and demand future methodological improvements. In addition, there is a great need for more accurate modeling of distortions and rigid-body shifts, as well as detection of errors in a given protein structure model. Error detection is useful both for refinement and interpretation of the models.