Insights into the three-dimensional (3D) structure of a protein are of great assistance when planning experiments aimed at the understanding of protein function and during the drug design process. The experimental elucidation of the 3D-structure of proteins is however often hampered by difficulties in obtaining sufficient protein, diffracting crystals and many other technical aspects. Therefore the number of solved 3D-structures increases only slowly compared to the rate of sequencing of novel cDNAs, and no structural information is available for the vast majority of the protein sequences registered in the SWISS-PROT database (nearly 60'000 entries in release 34). In this context it is not surprising that predictive methods have gained much interest.

Proteins from different sources and sometimes diverse biological functions can have similar sequences, and it is generally accepted that high sequence similarity is reflected by distinct structure similarity. Indeed, the relative mean square deviation (rmsd) of the alpha-carbon co-ordinates for protein cores sharing 50% residue identity is expected to be around 1Å. This fact served as the premise for the development of comparative protein modelling (also often called modelling by homology or knowledge-based modelling), which is presently the most reliable method. Comparative model building consist of the extrapolation of the structure for a new (target) sequence from the known 3D-structure of related family members (templates).

While the high precision structures required for detailed studies of protein-ligand interaction can only be obtained experimentally, theoretical protein modelling provides the molecular biologists with "low-resolution" models which hold enough essential information about the spatial arrangement of important residues to guide the design of experiments. The rational design of many site-directed mutagenesis experiments could therefore be improved if more of these "low-resolution" theoretical model structures were available.

Comparative protein modelling requires at least one sequence of known 3D-structure with significant similarity to the target sequence. In order to determine if a modelling request can be carried out, one compares the target sequence with a database of sequences derived from the Brookhaven Protein Data Bank (PDB), using programs such as FastA and BLAST. Sequences with a FastA score 10.0 standard deviations above the mean of the random scores or a Poisson unlikelyhood probability P(N) lower than 10^-5 (BLAST) can be considered for the model building procedure. The choice of template structures can be further restricted to those which share at least 30% residue identity as determined by SIM.

The above procedure might allow the selection of several suitable templates for a given target sequence, and up to ten templates are used in the modelling process. The best template structure - the one with the highest sequence similarity to the target - will serve as the reference. All the other selected templates will be superimposed onto it in 3D. The 3D match is carried out by superimposing corresponding Ca atom pairs selected automatically from the highest scoring local sequence alignment determined by SIM. This superposition can then be optimised by maximising the number of Ca pairs in the common core while minimising their relative mean square deviation. Each residue of the reference structure is then aligned with a residue from every other available template structure if their Ca atoms are located within 3.0 Å. This generates a structurally corrected multiple sequence alignment.

The target sequence now needs to be aligned with the template sequence or, if several templates were selected, with the structurally corrected multiple sequence alignment. This can be achieved by using the best-scoring diagonals obtained by SIM. Residues which should not be used for model building, for example those located in non-conserved loops, will be ignored during the modelling process. Thus, the common core of the target protein and the loops completely defined by at least one supplied template structure will be built.

The next step is the construction of a framework, which is computed by averaging the position of each atom in the target sequence, based on the location of the corresponding atoms in the template. When more than one template is available, the relative contribution, or weight, of each structure is determined by its local degree of sequence identity with the target sequence.

Following framework generation, loops for which no structural information was available in the template structures are not defined and therefore must be constructed. Although most of the known 3D-structures available share no overall similarity with the template, there may be similarities in the loop regions, and these can be inserted as loop structure in the new protein model. Using a "spare part" algorithm, one searches for fragments which could be accommodated onto the framework among the Brookhaven Protein Data Bank (PDB) entries determined with a resolution better than 2.5 Å. Each loop is defined by its length and its "stems", namely the alpha carbon (Ca) atom co-ordinates of the four residues preceding and following the loop. The fragments which correspond to the loop definition are extracted from the PDB entries and rejected if the relative mean square deviation (rmsd) computed for their "stems" is greater than a specified cut-off value. Furthermore, only fragments which do not overlap with neighbouring segments should be retained. The accepted "spare parts" are sorted according to their rmsd, and a Ca framework based on the five best fragments can be added to the model. In order to ensure that the best possible fragments are used for loop rebuilding, the rmsd cut-off can be incremented from 0.2 onwards until all loops are rebuilt.

Since the loop building only adds Ca atoms, the backbone carbonyl and nitrogens must be completed in these regions. This step can be performed by using a library of pentapeptide backbone fragments derived from the PDB entries determined with a resolution better than 2.0 Å. These fragments are then fitted to overlapping runs of five Ca atoms of the target model. The co-ordinates of each central tripeptide are then averaged for each target backbone atom (N, C, O) and added to the model. This process yields modelled backbones that differ from experimental co-ordinates by approx. 0.2 Å rms.

For many of the protein side chains there is no structural information available in the templates. These cannot therefore be built during the framework generation and must be added later. The number of side chains that need to be built is dictated by the degree of sequence identity between target and template sequences. To this end one uses a table of the most probable rotamers for each amino acid side chain depending on their backbone conformation. All the allowed rotamers of the residues missing from the structure are analysed to see if they are acceptable by a van der Waals exclusion test. The most favoured rotamer is added to the model. The atoms defining the c1 and c2 angles of incomplete side chains can be used to restrict the choice of rotamers to those fitting these angles. If some side chains cannot be rebuilt in a first attempt, they will be assigned initially in a second pass. This allows some side chains to be rebuilt even if the most probable allowed rotamer of a neighbouring residue already occupies some of this portion of space. The latter may then switch to a less probable but allowed rotamer. In case that not all of the side chains can be added, an additional tolerance of 0.15 Å can be introduced in the van der Waals exclusion test and the procedure repeated.

Idealisation of bond geometry and removal of unfavourable non-bonded contacts can be performed by energy minimisation with force fields such as CHARMM, AMBER or GROMOS. The refinement of a primary model should be performed by no more than 100 steps of steepest descent, followed by 200-300 steps of conjugate gradient energy minimisation. Experience has shown models optimised that energy minimisation (or molecular dynamics) usually move away from a control structure. It is thus necessary to keep the number of minimisation steps to a minimum. Constraining the positions of selected atoms (such as Ca, or using a B-factor based function) in each residue generally helps avoiding excessive structural drift during force field computations.