A restraint molecular dynamics and simulated annealing approach for protein homology modeling utilizing mean angles

Due to the enormous progress that has been made in genomics a large number of DNA sequences including many whole genomes have been published. The evaluation of these data must include the determination of the three-dimensional structures of the proteins encoded. Although the two experimental techniques capable of determining three-dimensional structures of proteins and other biomolecules at atomic resolution, namely nuclear magnetic resonance (NMR) and X-ray crystallography, have seen significant improvements the process of structure determination remains very time-consuming and difficult. Unless unexpected advances of these techniques will occur in future, it is obvious that for the majority of all the primary sequence data available three-dimensional structures cannot be obtained experimentally. Therefore, only computational approaches are capable of filling the gap between existing protein sequences and structures. Although considerable progress has been achieved in ab initio structural prediction strategies [1–3] they are in general still unreliable when atomic resolution is demanded. However, when structures of homologous proteins are available, the prediction of the three-dimensional structure of entire proteins and protein domains is rather successful. In light of the fact that the protein structures elucidated so far only show a remarkably limited number of folds it would be desirable to accelerate the structure determination process especially for proteins possessing a fold already known. According to the SCOP classification [4, 5] (release 1.65, 1 August 2003) 20619 protein structures stored in the Protein Data Bank share only 800 different folds. Comparison of different proteins with similar amino acid sequences showed that they quite often display very similar tertiary structures [6–9].

In the past several different homology modeling approaches were published which range from strongly interactive methods (model building) to fully automated methods (for reviews see e. g. [10] and [11]). Generally the starting point in these approaches is a search in structure databases such as the Protein Data Bank [12] or CATH [13] for all protein structures that are related to the target sequence and then to select those 3D structures that will be used as templates. For searching the structural databases one can employ pairwise sequence-sequence comparisons using for example programs such as FASTA [14] and BLAST [15]. When increased search sensitivity or a larger number of homologs are demanded methods which are based on multiple sequence alignments prove to be particularly efficient. Such an algorithm is implemented in the program PSI-BLAST [16]. An alternative strategy for homolog identification relies on so-called threading methods, which predict whether the target sequence adopts any of the known 3D folds. Threading methods should be useful in cases when no sequences can be found which are clearly related to the target [17].

When a list of related protein structures has been obtained the appropriate templates have to be chosen from these. In this procedure usually factors such as high overall sequence similarity between target and template sequences, quality of the template structure and conditions under which the template structure was obtained are taken into account. Then the selected templates have to be optimally aligned with the target sequence. Since the search methods mentioned above are usually optimized for detecting remote homologs they are not optimal for target-template alignment. A program often used for the latter type of alignments is CLUSTALX [18], which is also used within PERMOL.

Using the template-target alignment a variety of methods has been published for 3D model building. The group of methods which were developed first and are still frequently used were modeling by rigid body assembly [19–21]. Another group of methods use segment matching [22–25]. In the third, most recent group of methods spatial restraints obtained from the template structures are used in distance geometry calculations or energy optimization procedures to obtain the target model [26–31].

The PERMOL approach described presently also uses spatial restraints but in contrast to most other programs mainly dihedral angle restraints as opposed to restraints derived from inter-atomic distances are employed. These restraints enter molecular dynamics calculations in torsion angle space. In the following we will describe this method in more detail and mark differences to existing programs that have been published before.

In ab initio molecular dynamics (MD) simulations in addition to the applied force field only information about the amino acid sequence of the protein in study enters the calculations. While for small molecules such methods show results that are in very good agreement with the experimental data they mostly fail for more complex molecules. On the other hand restrained molecular dynamics calculations based on simulated annealing protocols are routinely and successfully used for the determination of solution NMR structures – in that case strong experimental information is available. Especially effective with regard to computational effort are calculations in torsion angle space as implemented in the programs DYANA [32] and CNS [33].

In this contribution we propose a method which combines the well-developed torsion angle dynamics calculations of DYANA or CNS with structural information extracted from three-dimensional structures of homologous proteins. This information is translated into conformational restraints. Local structural restraints are obtained by a weighted average of the backbone dihedral angles using an algorithm proposed by Döker et al. [34] These averaged dihedral angle restraints are usually well preserved within the local secondary structure elements and therefore are especially well suited for the modeling of these. The program MODELLER [28] for example also uses dihedral angle restraints in an optimization procedure but expresses them as so-called probability density functions which are derived from structural features in several families of homologous proteins. Global structural restraints are obtained from distance relations between carefully selected atoms of amino acids well separated in the primary structure. In contrast to other programs the distance restraints are mainly used for the global arrangement of the secondary structure elements which are defined by the dihedral angle restraints. The efficient structure calculations performed with DYANA or CNS allow calculating a large number of structural models in a relatively short amount of time. From the resulting ensemble of structures the best in terms of the DYANA target function or total energy (CNS) can be selected for further analysis. As has also been shown in NMR spectroscopy, it is useful to describe the target structure by an ensemble of model structures.

It should be noted that the PERMOL approach described here is related to the method detailed by Zhang et al. [35], which uses a combination of torsion angle dynamics and dihedral angle and distance restraints to predict the fold of helical proteins. In contrast to PERMOL the program from Zhang et al. uses methods for secondary structure and contact prediction to derive spatial restraints.

To benchmark the PERMOL approach we used it to determine a homology structure for the histidine-containing phosphocarrier protein (HPr) from E. coli of which the structure has been solved experimentally both by NMR [36] (PDB entry: 1HDN) and X-ray crystallography [37] (1POH). The homology model was compared to the target structures and to a homology model calculated with the program MODELLER [28]. To also investigate the performance of PERMOL on larger proteins that contain substantial disordered loop regions the human peroxisome proliferator activated receptor γ (Ppar γ) was used as a test case. Its structure has been determined previously by X-ray crystallography [38] (3PRG).

In this contribution we have presented a new program for homology modeling of protein structures. Using restraint molecular dynamics simulations together with spatial restraints derived from template structures we calculated homology structures of HPr from E. coli and of human Ppar γ. An advantage of the proposed method is the use of spatial restraints with individual upper and lower limits depending on the local structural conservation in the template structures. This becomes especially evident for the obtained bundle of Ppar γ model structures where one can easily distinguish between the mostly well-defined secondary structure elements and less ordered regions e.g. some of the larger loop regions.

At first glance it appears to be a disadvantage of the proposed method that not a unique, seemingly perfect structure is the result of the calculations as in the case of threading methods. However, the structure bundle produced by our approach gives an idea of the conformational subspace determined by the available experimental basis and the physical model. This is a safeguard against typical over-interpretations of model structures where data in badly predictable regions are used for the detailed interpretation of functional data or are used during the drug design process.

An additional advantage of the simulated annealing approach is that restraint violations are not treated explicitly but contribute to the overall "energy" which is minimized. In contrast to other methods in the approach used in PERMOL the mean torsion angles and their errors provide the main information. A few distance restraints are used to define the long-range relations which cannot be described sufficiently well by the local data. Accordingly, details of the selection of these restraints are not critical. Thus, the selection of pairwise restraints between all conserved residues seems to be plausible. The same is true for conserved hydrogen bonds. However, the PERMOL software also allows to define a custom selection of restraints and thus an adaptation to specific needs. As an example all hydrophobic contacts between amino acid residues observed in the template structures could be selected to serve as restraints. The automated calculation of individual weighting factors during the calculation of the expectation values and standard errors of the individual restraints would permit to introduce information about the local and global sequence conservation and the precision of the used structures. Currently, we are undertaking efforts to address this question. The high quality of the structure models generated with PERMOL illustrates that the same MD programs used for the determination of NMR structures can also be utilized for homology modeling. The programs and strategies developed for NMR structure determination have evolved to efficient optimizers even when only limited information (i. e. small number of structural restraints) is available. This has been recognized for example by Dominguez et al. [47] who use restrained molecular dynamics together with the ARIA protocol [48] for solving the docking problem. While in the case of NMR structure determination the restraints that enter the molecular dynamics simulation are derived from experimental observables like NOE cross-peaks, J- couplings, and residual dipolar couplings, in the case of homology modeling synthetic restraints are generated from previously determined structures of homologous template proteins. The use of standard MD programs and protocols also has a disadvantage since it is not possible to directly introduce properties in the calculation which are not provided for by the programs. An example would be the use of specific potential forms with multiple minima which describe the homology-derived information in more detail as it is done e. g. by MODELLER [28].

We compared the HPr homology structure we obtained with PERMOL to a structural model of HPr from E. coli calculated using MODELLER (version 6v2). When the same alignment file and template structures were used, homology models of similar quality were obtained with the two programs.

A specific advantage of the approach presented here is that it can be well used in the context of standard structure determination by NMR. The restraint files generated by PERMOL are editable and can be easily combined with other data and be adapted for use with different programs. As the same MD programs are used both for modeling with PERMOL and for NMR structure determination, incomplete experimental data can be conveniently combined with spatial restraints derived from homologous template proteins. The validity of the resulting structure models can be checked by calculating NMR R-factors [46]. Different force fields and annealing protocols which are available for the NMR MD programs can also be utilized for homology modeling. In this way recent advances like the structure refinement in explicit solvent [49, 50] can be readily exploited to derive more accurate homology structures.

In summary, we have presented a new method for homology modeling capable of producing high-quality structure models. Compared to many other homology structure prediction programs it is based on a different philosophy since its aim is not to predict a unique best structure but a bundle of structures representing the locally different degrees of reliability of the structure prediction. Since the homology-derived restraints are mainly used to reduce the conformational space to be searched by the MD calculation, their relative importance for obtaining a correct homology model is expected to decrease in future time as the physical model employed in these calculations is improved. Another advantage of the approach described here is its flexibility, conveniently allowing several template structures to be included as sources of structural restraints. Furthermore, the PERMOL software permits to determine which kinds of structural restraints enter the molecular dynamics calculation in a controlled fashion. We demonstrated that the standard MD programs used in the course of structure determination by NMR can also be well utilized for the purpose of homology modeling. Prediction on the basis of averaged torsion angles is a powerful tool which efficiently makes use of the structural information available in the protein data base and leads to well-defined structures.

Recently, a homology model determined with PERMOL was used in the resonance assignment [51] and structure determination process of a mutant form of HPr from S. carnosus [52] and to obtain an initial estimate for the molecular alignment tensor describing the partial orientation of the HPr molecule in anisotropic solution [53, 54]. PERMOL has also been integrated in the NMR structure determination package AUREMOL [39]. In this molecule-centered top-down approach one starts with a trial structure e.g. a homology model obtained by PERMOL that is iteratively refined until it fits the experimental data sufficiently as verified by the calculation of NMR R-factors.