MORPH-PRO: a novel algorithm and web server for protein morphing

Background Proteins are known to be dynamic in nature, changing from one conformation to another while performing vital cellular tasks. It is important to understand these movements in order to better understand protein function. At the same time, experimental techniques provide us with only single snapshots of the whole ensemble of available conformations. Computational protein morphing provides a visualization of a protein structure transitioning from one conformation to another by producing a series of intermediate conformations. Results We present a novel, efficient morphing algorithm, Morph-Pro based on linear interpolation. We also show that apart from visualization, morphing can be used to provide plausible intermediate structures. We test this by using the intermediate structures of a c-Jun N-terminal kinase (JNK1) conformational change in a virtual docking experiment. The structures are shown to dock with higher score to known JNK1-binding ligands than structures solved using X-Ray crystallography. This experiment demonstrates the potential applications of the intermediate structures in modeling or virtual screening efforts. Conclusions Visualization of protein conformational changes is important for characterization of protein function. Furthermore, the intermediate structures produced by our algorithm are good approximations to true structures. We believe there is great potential for these computationally predicted structures in protein-ligand docking experiments and virtual screening. The Morph-Pro web server can be accessed at http://morph-pro.bioinf.spbau.ru.


Background
The number of solved protein structures in PDB [1] has grown enormously in recent years. However, the function of many proteins is highly correlated with their movement. X-Ray crystallography, which contributes most of the structures in PDB, gives us only a static view of protein structure. Recent developments in computational protein morphing [2][3][4] provide visualization of a molecule transitioning from one conformation to another by producing a series of intermediate conformations. In this paper we present a novel, computationally efficient algorithm for generating intermediate structures between two solved conformations of the same protein. In addition, we *Correspondence: kira@math.spbu.ru 4 Algorithmic Biology Laboratory, Saint Petersburg Academic University, Saint Petersburg, Russia Full list of author information is available at the end of the article explore the possibility that intermediate structures generated in the morphing procedure may also represent realistic approximations of the actual protein conformational change, including the structures of the intermediate conformations.
Various attempts to predict the trajectory of proteins through conformational space have been made. Some success has been achieved through the use of elastic network models [5,6]. However, the accuracy of these methods depends on the chosen starting conformation (either apo-or holo-) and collectivity of the atoms in the motion [7]. Other attempts require numerous iterations of energy-minimization [8], which can be computationally expensive. Molecular dynamics simulations [9] may also be useful in determining the nature of conformational changes, but currently require significant computing power. Furthermore, motion planning techniques can http://www.almob.org/content/8/1/19 be adapted to model molecular motions [10][11][12], providing an attractive alternative to the mentioned approaches due to their efficiency.
The most widely-used application to produce protein morphs is the Morph Server developed by Krebs and Gerstein [8]. The goal of the Morph Server is to provide visualization and classification of protein movements. Our emphasis is on the fast generation of intermediate structures that represent realistic conformations.
Given two aligned proteins as input, our MORPH-PRO algorithm produces a series of intermediate conformations. We use linear interpolation, so that at each step every residue will move along the straight line between its current position and its ending position. Unfortunately, this can lead to biologically infeasible intermediate structures with atoms occupying the same space, incorrect bond lengths, and incorrect bond angles. Therefore, we use the atom positions generated by linear interpolation as a first approximation to the correct solution, and use a dynamic programming algorithm to ensure that certain biological constraints are satisfied. This produces structures which better resemble real proteins. Because these techniques are very efficient, our algorithm can produce many intermediate structures very quickly.
The intermediate structures produced by morphing algorithms show great promise in molecular docking [13]. Molecular docking, which uses computer simulations to model and score protein-ligand binding, is a critical tool for drug discovery. Protein flexibility is believed to play a significant role in ligand binding [14]. One method for including flexibility in the docking experiment is to perform ensemble docking [15], which uses multiple conformations of the protein for evaluation. Performing docking against several conformations of a protein has been shown to provide better screening results, than against a single static structure [16]. The intermediate structures produced by morphing algorithms may improve our ability to detect these ligands, and therefore aide in the development of drug-like molecules [17].

Methods
In this section we analyze the simplest form of the morphing problem and present our MORPH-PRO algorithm. We designate P start and P end as the sequences of 3-D coordinates of the Cα atoms for the starting and ending conformations. For simplicity, we assume that proteins P start and P end have an equal number of residues, and are aligned in 3-D. Later we will discuss the situation where P start and P end do not meet these conditions and will address various extensions to the simplest model of the protein morphing problem.

Morphing algorithm
We represent a sequence of n points in 3-D (n-tuple) as a 3 · n matrix (p ij ), where p ij is the i-th coordinate of the j-th point. Let n be the number of residues in P start and P end . Given a parameter α, we define the α-intermediate of proteins P and P as (1 − α) · P + α · P . The simplest way to morph P start into P end is to generate intermediate reconstructions (1−α)·P start +α·P end for 0 < α < 1. However, some α-intermediates may not look like real proteins, for example they may consist of consecutive Cα atoms at biologically impossible distances. Below we show how to solve the protein morphing problem thereby transforming every intermediate reconstruction (being a sequence of n points) into a protein-like sequence of points. At each iteration, every point first moves by an appropriate distance towards its ending position, and then the obtained sequence of points is adjusted to become protein-like.
The pseudo code of the algorithm for generating K protein-like sequences P 1 . . . , P K of points is as follows: procedure Morph(P start ,P end , K ) Below we describe the algorithm for transforming a sequence of points P into a protein-like structure Proteinize(P).

Optimal equidistant sequence problem
Given a sequence P of n points, we define d j (P) as the distance between the (j)-th and the (j + 1)-th points in P: d j (P) = Protein structures exhibit a strict distance constraint between consecutive Cα atoms that are 3.8 Å apart within an error margin of 0.1 Å. A sequence of points is protein-like if it is (3.8,0.1)-equidistant. We note that the consecutive Cα atoms in cis-proline do not adhere to this distance rule, and these cases are not handled by our algorithm. We define the distance d(P, P ) between two sequences P and P , of n points each, as n An (a, )-equidistant sequence P is called an optimal (a, )equidistant approximation of P if d(P, P ) is minimum among all possible (a, )-equidistant sequences P . Below we describe an approximate solution to the following problem: Optimal Equidistant Sequence Problem (OESP): Given a sequence of points, find its optimal equidistant approximation. http://www.almob.org/content/8/1/19

Solving OESP
Here we describe an approximate OESP algorithm that assumes the space of possible solutions is discretized. For each point from the sequence P, we construct a lattice of 3-D points centered around it, as shown at Figure 1. Thus, each lattice is local to its corresponding point from P, which distinguishes our approach from naive and outdated attempts to understand protein folding which utilize a global lattice [18][19][20]. The selection of the number of points in the lattice and the edge length is discussed later.
Let v i,j be the i th vertex in the lattice constructed around the j th point. Let v 0,j be the vertex corresponding to the j-th point in P. Let Q be the number of vertices in each lattice.
We construct a directed edge from a vertex v i,j to a vertex v g,j+1 for 1 ≤ i, g ≤ Q and 1 ≤ j ≤ n − 1. The score of an edge is defined as: We also assign a score to each vertex, v i,j , Finding a protein-like sequence P of points which minimizes d(P, P ) translates into finding the path with the minimum score through the graph starting in the first lattice and ending in the n th lattice. The score of a path is defined as the sum of the scores of its edges and vertices. Let PATH(v i,j ) be the value of the minimum scoring path among those that start in the first lattice and end at vertex v i,j . Variable PATH(v i,j ) can be computed using the following recurrence: The score of the protein-like sequence of points which is closest to our original approximation is then The solution of OESP can be determined by backtracking. The time complexity of generating a protein-like conformation of Cα atoms from a collection of n points, if one exists, is O(nQ 2 ).

Angle and proximity constraints
The above approach solves OESP and produces a (3.8,0.1)equidistant sequence. There is more, however, to consider when defining a protein-like structure than consecutive residue distance. We now redefine the notion of a proteinlike sequence of points to take into account consecutive residue angles and proximity constraints.
Given 3-D points q 1 , q 2 , and q 3 , a function ang(q 1 , q 2 , q 3 ) is defined as the minor angle in degrees created by the lines through q 1 and q 2 and through q 2 and q 3 , respectively. Given a sequence P of n points, we let ang j (P) = ang( We observed that most Cα angles in real proteins fall in the range of 70°to 120°.
Furthermore, a sequence P of points is z-distance consistent if the distance between any two non-consecutive points in P is at least z Å. We determined that a distance of 2.0 Å was typical in real proteins.
We introduce a new score to evaluate the angle defined by three vertices, v 1 , v 2 , and v 3 .
In order to incorporate angles into our algorithm, we must use a more complex recurrence which relies on both the current vertex, v i,j , and a preceding vertex, v h,j−1 . We define PATH(v i,j , v h,j−1 ) as the path with minimum score among all paths that start in the first lattice, end in v i,j , and pass through v h,j−1 . We replace (2) with the following for To determine the score of the protein-like sequence of points which is closest to our original approximation, we find: This construction does not force the sequence of points to be 2.0-distance consistent. For this, we apply a heuristic, which increases the VScore of vertices which are close to other lattices. We replace (1) with We chose the multiplier 100 because it worked well to prevent Cα clashes in our morphs. The addition of the angle and distance constraints requires O(n 2 Q 3 ).
However, the advanced strategy described above may be impractical if the proteins being examined are large or the conformational change is dramatic. Therefore, we also considered a simplified strategy which can significantly improve the running time. In the simplified strategy, (2) is replaced with is the vertex preceding v h,j−1 in the best path ending at v h,j−1 , the score of which is determined by the value of PATH(v h,j−1 ). Similar to the the basic method, the score of the optimal protein-like sequence of points is and thus, the time complexity of the simplified strategy is also O(nQ 2 ). The simplified strategy may provide a sub-optimal intermediate structure. However, if a structure is produced, it obeys both the angle and proximity constraints. It should be noted that the simplified strategy may fail to find a solution to OESP instances, even when a solution can be found by the advanced algorithm. The advanced algorithm looks for an optimal path among all feasible ones stretching from the first to the last lattice, while the former takes into consideration only a subset of paths. In addition, the simplified strategy may require an increase of the lattice size (see Parameter Selection), thus reducing the difference in the running time in practice of the algorithms.
Our experiments described in detail below were carried out using the simplified strategy.

Preprocessing
Our algorithm only interpolates intermediate positions for residues which are aligned. Therefore, if the input proteins have different lengths we use the Needleman-Wunsch global sequence alignment algorithm [21] to align them, and reduce our starting and ending conformations to include only positions that are aligned. We chose to use a sequence-based alignment method because P start and P end are likely related proteins and will have similar sequences. The output of this phase of the algorithm is a set of coordinates of aligned Cα's for P start and P end . In this situation, the i th residue in the alignment may not correspond to the i th residue in P start . If the i th and (i + 1) st residues produced from the alignment are not consecutive in P start then EScore for the edge connecting them is 0. Similarly, if either the (i−1) th and i th or the i th and (i+1) st residues are not consecutive in P start then AScore for the angle at the i th residue is 0.
In order for the morphing algorithm to work, the proteins should be aligned in 3-D using a structure alignment program. In the implementation we used for the experiments described in this paper, this task is accomplished by Kabsch's algorithm [22] (also see [23]). Our server uses the Quaternion Characteristic Polynomial (QCP) method recently proposed by [24].

Parameter selection
For our experiments we set the number of intermediate structures, K, to be the rounded displacement of the largest Cα movement. For example, if the greatest movement of any Cα from the starting conformation and the ending conformation is 15.2 Å, then K = 15. This results in only small differences between consecutive structures.
We selected the edge length and point density for the lattices based on experimental evidence. Increasing the density of vertices in the lattice allows for a finer grained set of possible coordinates, but we found that a density higher than 6 points per Å (216 points per Å 3 ) does not produce significantly better intermediate structures. Consequently, we fixed the density at 6 points per Å. The length of the lattice edge is set initially to 1 Å. However, if OESP solution cannot be found at this lattice size, we increase the lattice edge length (to 1.5 Å and then to 2.0 Å). If an OESP solution cannot be found with lattice edge length of 2.0 Å then our algorithm will not produce a morph.

Server implementation
We implemented the MORPH-PRO server using an open source web framework Ruby on Rails and SQLite3 database engine, and a new 3D graphics standard WebGL. The algorithm for protein morphing was implemented in ANSI C. We used BioRuby [25] -an open source bioinformatics library for Ruby -for parsing PDB files, and the QCProt 1.3 realization of the QCP algorithm for aligning proteins in 3D, distributed under a BSD open source license.
The interface allows a user to upload two PDB files containing the starting and the ending conformations, and either to explicitly indicate the number of intermediate http://www.almob.org/content/8/1/19 conformations or to let it be determined automatically (based on the maximum Cα displacement, as described in Parameter Selection). After the intermediate conformations are computed, the morphing process can be visualized either as a movie or step-by-step. A transformation between two consecutive conformations is accomplished via linear interpolation. A 3D chain representing a conformation can be rotated, and zoomed in and out. In addition, a user can choose an appropriate level of detail for rendering and elect to use the full algorithm or the simplified version. A publicly available archive of submitted morph requests is stored on the server in an SQLite3 database, making it easier to re-run the algorithm on the same input.

Results and discussion
We evaluate our morphs by looking at both the biological feasibility of each individual structure, as well as the series of structures as a whole. We evaluate our morphs by comparing to proteins which have 3 or more solved structures in PDB, as proposed by [26]. In many instances, multiple conformations of the same protein are not available. Instead, we used proteins from the same family with nearly identical sequences as endpoints in our morph.

Pyrophosphokinases
We created a morph between two members of the pyrophosphokinase family (PDB codes: 1DY3, 1RAO). The alignment produced 158 residues with a maximum Cα displacement of 22 Å. The RMSD between the starting structure and the ending structure is 4.07 Å.
We examined each intermediate structure produced from this morph, and looked for clashing Cα atoms. None of the intermediate structures had atoms within 2 Å of another atom. We also looked at torsion angles created by Cα atoms. The Ramachandran plot of phi versus psi angles of the intermediate structure, which occurs halfway through the morph, is shown in Figure 2. The majority of the points in the plot fall within a region that is observed in real proteins. This indicates that our structure exhibits characteristics of real proteins.
It is also beneficial to look at the intermediate structures in the context of the entire morph. We have shown that our intermediates are protein-like, and we now demonstrate that the series of intermediate structures closely mimics the series of conformations a protein would visit. If multiple conformations of the same protein are known, then we can compare our predicted trajectory to the solved trajectory by calculating the RMSD between our  [28], and energy minimization was performed using Swiss-PDB viewer [29]. http://www.almob.org/content/8/1 /19 intermediates and the experimentally solved intermediates. However, alternate conformations were not available for these proteins, so instead we used solved structures for proteins in the pyrophosphokinase family.
We chose two additional pyrophosphokinases to act as 'experimental' intermediates (PDB codes: 1RB0, 1HKA). We chose these proteins because they can be ordered by their RMSD between 1DY3 and 1RAO, and therefore are likely to be similar to the trajectory the morph should take. We plot the RMSD of our intermediate structures against each of these four proteins in Figure 3.
Intermediates which are produced early in the morph are closest to the starting protein, 1DY3, while those that are produced late in the morph are closest to the ending protein, 1RAO, as expected. Our intermediates from the middle of the morph become close to both 'experimental' intermediates, 1RB0 and 1RAO, suggesting that our movement closely follows the evolutionary changes which occurred between the two proteins. In addition, the intermediate structures generated by our algorithm come roughly as close, if not closer, to the known homologs as those produced by Morph Server, as demonstrated in Table 1. A direct speed test with the Morph Server was not possible because a fully functional standalone tool was not available.

F1-ATPase
The technique of looking at RMSD of the intermediate structures to known structures is most useful when X-Ray structures of actual intermediate conformations are available. There are three conformations solved for the F1-ATPase molecular motor (PDB code: 1E79) which exhibit a subtle change. The RMSD between the starting and ending conformations is 1.78 Å. The protein has 492 residues and the largest movement of a Cα is 11 Å. We produce a morph of 11 total structures from 1E79A to 1E79C.
The intermediate structures are very similar to all of the known structures, with RMSD consistently less than 2 Å. We do, however, see our intermediate structures become closer to the known intermediate 1E79B.
One intermediate structure comes as close as 1.61 Å, while the starting structure (1E79A) is 1.85 Å and the ending structure (1E79C) is 1.73 Å. Figure 4 demonstrates how the predicted intermediates are similar to the starting structure early in the morph, become more similar to the known intermediate structure in the middle of the moprh, and then finally become similar to the ending structure. In Figure 4, we generated 30 intermediate structures to better illustrate this point.

GroEL
Our algorithm also performs well on large proteins. GroEL proteins chaperon the folding of other proteins. Two GroEL proteins (PDB codes: 1GRL and 1AON) exhibit a simple morph on 515 aligned residues, changing from a closed conformation to an open conformation. The RMSD between these two structures is 12.36 Å while the largest movement of a single Cα is 34.8 Å. Despite the large number of atoms and the significant movement, the morph took only a couple minutes to run. Figure 5 shows the initial conformation, the final conformation and 2 out of 34 intermediate structures produced in this morph.

Virtual screening
Virtual screening [30] is a technique which simulates the binding of a protein and a ligand, in order to determine the best ligand candidates from a large database. Most often, virtual screening is used as part of a drug development pipeline, guiding the selection of likely drug candidates. The predicted binding affinity of a ligand for a protein is determined by a docking algorithm, which finds the orientation and location of the ligand with respect to the protein. Modeling protein flexibility is very difficult due to the large degrees of freedom of a protein structure [13,31]. One promising approach to implicitly incorporating protein flexibility is to dock against an ensemble of static protein structures [32]. If multiple conformations of the target protein are solved using NMR or X-Ray studies, these are good candidates for ensemble docking. However, in the more common case of unknown intermediate conformations a computational method can provide accurate models more quickly. Use of computationally-produced intermediates in virtual screening has shown promising results [33].
To test the potential for our intermediate structures in virtual screening we examined docking scores of our structures versus those solved experimentally against a small database of ligands. First, we produced a morph of the c-Jun N-terminal kinase 1 (JNK1). The starting conformation of this protein (1UKH) was solved complexed with a peptide (pepJIP1) derived from the binding portion of the scaffolding protein JIP1. The ending conformation (1UKI) was solved complexed with pepJIP1 and the ATP mimic SP600125. The binding of pepJIP1 to the JIP1 binding site on JNK1 causes a small conformational change at the ATP site. Though the movement is small, it produces a morph of 3 intermediates (P 2 , P 3 , P 4 ) in addition to the starting and ending conformations. The absent backbone atoms and side chains of each intermediate structure were reconstructed using Maxsprout [28], and energy minimization was performed using Swiss-PDB viewer [29]. As a basis for comparison, the X-Ray structures of 1UKH and 1UKI were also reduced to their Cα's and then reconstructed in the same manner to produce P 1 and P 5 , respectively.
Next, we performed docking with GOLD [34], a commonly used docking program and scoring scheme, on four ligands (extracted from PDB) known to bind to JNK1, as well as SP600125. Table 2 shows the rankings of the binding affinities from highest to lowest based on the GoldScore. The headings are the PDB codes for the solved structures of JNK1 complexed with each ligand.
The first column behaves as expected. The structure which has the highest binding affinity for SP600125 is 1UKI which is the structure of JNK1 complexed with SP600125. The X-Ray structures docked with SP600125 rank significantly higher than the reconstructed P 1 and P 5 . This suggests that better side chain reconstruction could greatly improve the docking results.
For three of the other ligands, the second intermediate structure, P 2 scores higher than any other intermediate structure as well as any X-Ray structure. This demonstrates that our intermediate structures would be more likely to identify ligands which bind to JNK1 than either of the two X-Ray structures.

Conclusions
It is clear that there is much to learn about the nature of protein structure dynamics that is not addressed in the static information contained in PDB. The intermediate structures representing a protein as it moves from one conformation to another may yield much information about how a protein functions. Experimental techniques are inadequate for this task due to practical and technological limitations. For this reason, structural biology The morphed structures also show promise in the area of virtual screening. Most techniques limit protein flexibility to the side chain atoms, and may allow limited flexibility of the substrate. Our morph produces intermediate structures which are hypotheses for possible backbone movements. For this reason, some ligands bound more favorably to our intermediate structures than the solved structures. These are strong implications for the potential of morphs in guiding drug development.
Like all other approaches, our algorithm also has limitations. Linear interpolation, with only small corrections, prevents our method from correctly producing a morph for proteins with very large or complex movements. Many of these morphs could be solved by allowing a larger movement from the first approximation (a larger lattice), or allowing higher granularity of possible Cα positions (more points in each lattice) but the time cost would be significant. Clearly, in protein morphing there is a trade-off between speed and accuracy. http://www.almob.org/content/8/1/19 The manuscript was prepared by NEC. PAP and KV contributed to finalizing the manuscript. All authors read and approved the final manuscript.