An enhanced partial order curve comparison algorithm and its application to analyzing protein folding trajectories

Related work

Previously, folding simulations analysis is performed mainly for testing various protein folding models [5–7], such as the folding pathway model and the funnel model; and/or for studying energetic aspects of folding kinetics [8–11]. The geometric shapes of the conformations involved in folding trajectories have not been widely explored [4, 12, 13], despite their important role in folding. A particularly interesting work in this direction is by Ota et al. [4], where they investigated the folding trajectories of a mini-protein Trp-cage using phylogenic tree combined with expert knowledge. However in general, an automatic tool to facilitate the folding simulations analysis at large scales is still missing. This paper provides an important step towards this goal by modeling folding trajectories as curves and using a new multiple curve comparison (MCC) algorithm to detect critical folding events.

The closest relative of our MCC problem in computational biology is the multiple structure alignment (MSTA) problem, which aims at aligning a family of protein structures, each modeled as a three dimensional polygonal curve to represent its backbone.

MSTA is a very hard problem. In fact, even the pairwise comparison problem of aligning two structures A and B is believed to be NP-hard since one has to optimize simultaneously both the correspondence between A and B and the relative transformation of one structure with respect to the other. Numerous heuristic-based algorithms have been developed in practice for this fundamental problem [14–20]. If we have a set of k > 0 structures, then even the problem of aligning them optimally without considering transformations becomes intractable – it takes Ω(n ^k ) time using the standard dynamic programming algorithm, where n is the size of each protein involved.

In practice, progressive methods are widely used to attack the MSTA problem [21]. For example, given a set of structures, many approaches start with a seed structure and then progressively align the remaining structures onto it one by one [22–28]. A consensus or core structure is typically built throughout, to maintain the common substructures among the proteins that are already aligned. At each round, usually only pairwise structure comparison is performed to align the current consensus with a new structure.

The above progressive MSTA framework is a greedy approach. Its performance depends on the underlying pairwise comparison methods used, the order of structures that are progressively aligned, as well as the consensus structure maintained. Various heuristics have been exploited to find a good order for the progressive alignments. Note that this order can also be guided by a tree instead of a linear sequence, which removes the need of choosing a seed structure. The progressive procedure may also be iterated several times to locally refine the multiple structure alignments.

Our results

There are two main differences between the MCC problem we are interested in and the traditional MSTA problem. In the case of protein structures, it is usually explicitly or implicitly assumed that the (majority of the) input proteins belong to one family (How to classify a set of input structures into different families is a related problem, and many such classifications exist [17, 29, 30]), or at least share some relations. As such, one can expect that some consensus of the family should exist. However in our case, the set of curves are from a set of simulations including both successful and unsuccessful runs, and we wish to classify this diverse set of curves, and capture common features within as well as across its sub-families. Secondly and more importantly, the level of similarity existing in these folding trajectories is usually much lower than that in a family of related proteins. Hence we aim at an algorithm with high sensitivity, which is able to detect small-scaled partial similarity and handle multi-dimensional curves (trajectories) as well.

In this paper, we propose and develop a sensitive MCC algorithm, called the EPO (enhanced partial order) algorithm, to compare a set of diverse high dimensional curves. Our algorithm follows a similar framework as the POA algorithm [3, 28] to encode the similarities of aligned curves in a partial order graph, instead of in a linear structure used by many traditional MSTA algorithms. This has the advantage that other than similarities among all curves, similarities among a subset of input curves can also be encoded in this graph. See Figure 1 for an example, where nodes in both graphs represent a group of aligned points from input curves.

For the more important problem of sensitivity, we observe that being a greedy approach, the progressive MSTA framework tends to be inherently insensitive to low level of similarities – if one early local decision is wrong, it may completely miss a small-scaled partial similarity. To improve this aspect of the performance of the progressive framework, we first propose a novel two-level scoring function to measure similarity, which, together with a clustering idea, greatly enhances the quality of the local pairwise alignment produced at each round. We then develop an effective merging step to post-process the obtained alignments. This step helps to reassemble vertices (high dimensional points) from input curves that should be matched together, but were scattered in several sub-clusters in the alignments due to some earlier non-optimal decisions. Both techniques are general and can be used to improve the performance of many existing MSTA algorithms. Experimental results show that our MCC algorithm is highly sensitive and able to classify input curves. We also demonstrate the power of our tool in mining critical events from protein folding trajectories using a detailed case study of a miniprotein Trp-cage.

Although our EPO algorithm is developed with the goal of comparing folding trajectories, the algorithm is general and can be applied to other domains as well, such as protein structures or pedestrian trajectories extracted from surveillance videos [31]. In particular, in this paper, we demonstrate this generality by showing that the EPO algorithm can be used to improve the results of existing multiple protein structure alignment algorithms, especially when input proteins share low structural similarity.

Experimental Results on trajectory data

In this section, we report a systematic performance study on a biological dataset that contains 200 molecular dynamics simulations. The experiments achieve the following goals: First, we show that the quality of the alignments produced by our EPO algorithm is significantly better than that of the original POA algorithm. Second, we demonstrate the effectiveness of our algorithm by applying it to real protein simulation data and obtaining biologically meaningful results that are consistent with previous discoveries [4].

Background of dataset

Our input dataset includes 200 simulated folding trajectories for a particular protein called Trp-cage. The dataset is provided by the Ota's Lab [4]. The folding simulations were performed at 325 K by using the AMBER99 force field with a small modification and the generalized Born implicit solvent model. Trp-cage (see Figure 2) is a mini-protein consisting of 20 amino acids. It has been widely used for folding study because of its short, simple sequence and its quick folding kinetics. Following the definition from [32], a successful folding event has to satisfy the following two criteria:

• The RMSD for a conformation from the native NMR structure [1] is less than 2 Å.

• A subsequence of such near-native conformations holds for at least 200 ps.

In [4], 58 successful folding trajectories reaching successful folding events are identified, and each trajectory includes 101 successive conformations sampled at 20 ps interval. Furthermore, there are two crucial observations in [4] that we will examine in the our experiment. First, before moving to the native conformation, a "ring" sub-structure (see Figure 2) has to be formed. Second, the distinction of native and pseudonative confirmations heavily relies on side-chain position of "ring" sub-structure. Ota et al. [4] obtained the above results by aligning each pair of trajectories first and then applying a neighbor joining method to group similar trajectories together. However this semi-automatic approach requires dedicated expert knowledge. The following experiments applied on the same dataset show that our EPO algorithm can automatically detect the above folding events with little prior knowledge.

Investigation on entire protein structure

In the first set of experiments, we convert each conformation to a high dimensional point (i.e, a 20 × 20 = 400 dimensional point), based on the distance matrix between all of the alpha-carbon atoms. Figure 3 compares the quality of the alignments of the SuccData by performing the POA algorithm, our EPO algorithm without the merging procedure (EPO-NoMerge), and the EPO algorithm. It shows the number of aligned nodes (y-axis) versus the size of aligned nodes (x-axis). Note that EPO-NoMerge is essentially POA with a clustering preprocessing and the new two-level scoring function.

The similarity level between these trajectories is low (i.e, the number of aligned nodes with large size is small). It is clear from this histogram that our EPO algorithm significantly outperforms the other two by producing more aligned nodes with large sizes. The comparison between EPO and EPO-NoMerge demonstrates the effectiveness of our merging procedure, and that EPO-NoMerge is better than POA shows that the two-level scoring function as well as the clustering preprocessing greatly enhances the performance. We have also performed experiments which show that compared to the POA algorithm, EPO-NoMerge is much less sensitive to the order of curves aligned. Comparisons of the three algorithms over the MixData produces a similar result, and majority of points aligned to heavy nodes (i.e, |o| ≥ 40) are from successful runs (the results are not shown in this report).

We also observe that most of the heavily aligned nodes are close to the end of the trajectories for the SuccData. In fact, many aligned points have conformation IDs around and greater than 90, which is indeed the time that the folding starts to get stabilized. More specifically, consider the set of aligned nodes of size greater than 40 for the SuccData. Among all points aligned to these nodes, 67.2% has a conformation ID greater than 90, and 24.4% has an ID between 80 and 90. This implies that our algorithm has the potential to detect the stabilization of successful folding events in an automatic manner.

This also implies that using the entire protein structure may be too coarse to detect critical folding events, as they are usually induced by small key motifs. In what follows, we map only a substructure of the input protein into a multi-dimensional point and provide more detailed analysis of this folding data.

Ring-substructure

The second motif corresponds to the neck of a ring structure. In particular, it consists of the sub-chains of No. 2 – 5 and No. 16 – 19 amino acids. The following results demonstrate that EPO can automatically not only find but also track the formation of such fingerprint sub-structures (critical motif).

First, we observe from Table 1 that when applying the EPO algorithm to the MixData (with the sidechain weight factor α = 0.9), signifificant alignments involve mainly trajectories from SuccData. For example, the last row of Table 1 shows that among the 62 points (from 62 trajectories) aligned to a particular node, 58 are from SuccData, with the remaining 4 from FailData. Hence this motif is potentially critical to the success of the folding of Trp-cage. It also suggests that we can automatically classify the MixData into SuccData and FailData with few false positives based on this ring-neck motif, while previously, the classification in the input data was obtained by a few expert defined rules.

		Classification
Aligned Node ID	|o i | ≥ τ	SuccData	FailData	D(o) Å
1	49	27	22	1.852
2	45	28	17	1.798
3	41	29	12	2.189
4	40	31	9	1.447
5	48	31	16	1.761
6	40	32	8	1.322
7	47	34	13	1.133
8	42	35	7	1.923
9	44	36	8	0.873
10	49	42	7	1.428
11	54	48	6	1.020
12	59	50	9	1.294
13	60	51	9	0.932
14	56	52	4	1.255
15	62	56	6	1.782
16	62	58	4	1.503

Column 2 – 4 shows the size of an aligned node (i.e, the number of points aligned to this node) from MixData, SuccData, and FailData, respectively. Column 4 shows the diameter of this node (note that the distance threshold ε = 1 Å means that the diameter of a node can be up to 2 Å).

Second, when the side-chain weighting factor α = 0.9, it turns out that 49.6% of significant aligned nodes are formed before the conformation ID 85 (compared to results from Section). For example, there are two aligned nodes from the successful runs, where 80% of points (i.e. trajectories) aligned to them has a conformation ID between 75 – 85. This implies that the complete formation of this ring-neck usually immediately precedes the stabilization of the folding structure (which is roughly at conformation ID 90 for successful trajectories).

If reducing the side-chain weighting factor α to 0.5, naturally, we found more aligned nodes. In particular, other than the cluster with conformations of IDs around 80, we observe more significant clusters involving conformations with IDs from 50–70. By comparing the conformations of the ring-neck motif in these clusters with those in the aligned nodes around 80, we found that the backbone structures are rather similar, but the side-chains are of different orientations. In other words, the shape of the ring-neck motif is first stabilized by the backbone structure, and then the side-chains gradually move into right position. There are a few trajectories where the side-chains eventually move to the mirror image of their correct positions, and lead to pseudo-native conformations which can be detected when considering the side-chains.

Figure 4 displays several groups of critical events identified by the EPO algorithm (corresponding to the aligned nodes as shown in Table 1). In particular, 4(a) includes two closely occurred events during the early stage of the folding procedure (one of the conformations is selected from the aligned node 1, and the other one is from aligned-node 2). At this time, the sequence has started to fold and one can observe the helical structure, but the ring is not yet formed. 4(b) presents two conformations representing aligned nodes 3 and 4, respectively. We observe that the ring has started to form at this point, but is still not stable. 4(c) shows several conformations (one each from aligned nodes 5 to 15) occurred in order in most successful runs. At this stage, the ring is adjusted and stabilized. The adjustment mainly happens around the turn area and for side chains. 4(d) shows the final successful structure.

The above results are consistent with the results from [4], where such a ring-shaped substructure was discovered semi-automatically by pairwise structure comparisons together with expert knowledge.

Experimental Results on protein structure data

In summary, the above results are similar to the case of comparing multiple folding trajectories. In particular, the two-stage scoring function exploited by our EPO algorithm is effective and alleviates the ordering issue in progressive MSTA. The merging stage helps to produce better correspondences between input structures and makes the algorithm robust. Hence the EPO algorithm can be used as a post-processing tool for current MSTA softwares to refine and optimize their alignment results.

Notations and algorithm overview

Before we formally define the MCC problem, we introduce some necessary notations. Given a set of elements V = {v₁,..., v _l }, a relation $\prec$ over V is transitive if v _i $\prec$v _j and v _j $\prec$v _k imply that v _i $\prec$v _k . In this paper, we also refer to v _i $\prec$v _j as a partial order constraint. A partial order graph (POG) G = (V, E) is a directed acyclic graph with V = {v₁,..., v _l }, where v _i $\prec$v _j if there is an edge (v _i , v _j ). Note that by the transitivity of this relation, two nodes may have a partial order constraint even when there is no edge between them in G. Let R be the set of partial order constraints induced by G. We say that a permutation Π(V) of V is a partial order list w.r.t. G if for any v _i $\prec$v _j ∈ R, we have that vi appears before v _j in the permutation Π(V). In other words, the linear order in Π(V) is a total order satisfying all partial order constraints induced from G. See Figure 5 for an example.

Let $\mathcal{T}$ = {T₁,..., T _N } be a set of N trajectories in ℝ ^d , where each trajectory T _i is an ordered sequence of n points $p_1^i,\mathrm{...},p_n^i$. (For simplicity, we assume without loss of generality that all T _i s have the same length n.) The goal of the MCC algorithm is to find aligned sub-sequences from $\mathcal{T}$.

More formally, an aligned node o is a collection of vertices from T _i s, with at most one point from each T _i . Given a 3-tuple ($\mathcal{T}$, τ, ε), where τ and ε are input thresholds, an alignment of $\mathcal{T}$ is a POG G with the corresponding set of partial order constraints R and a partial order list of aligned nodes $\mathcal{O}$ = {o₁,..., o _L } such that the following three criteria are satisfied:

C1. |o _k | ≥ τ, for any k ∈ [1, L];

C2. for any $p_j^i,p_{j^{\prime }}^{i^{\prime }}\in o_k$, $\left|\right|p_j^i-p_{j^{\prime }}^{i^{\prime }}\left|\right|\le \epsilon$;

C3. if $p_j^i\in o_{k_1}$ and $p_{j^{\prime }}^{i^{\prime }}\in o_{k_2}$ with $o_{k_1}\prec o_{k_2}$, then j <j'.

(C1) indicates that the number of vertices of input curves aligned to each aligned node o _k is greater than a size threshold τ . (C2) means that these aligned points are tightly clustered together (i.e, the diameter of them is bounded by a distance threshold ε). (C3) enforces that points in different aligned nodes still maintain their partial order along their respective trajectory, which means that o _k s are inherited and thus consistent to the points in each trajectory. Our goal is to maximize L, the size of such an alignment $\mathcal{O}$. See Figure 6(b) for an example of an alignment graph.

Initial POG construction

Standard dynamic programming (DP) [35, 36] is an effective method for pairwise comparison between sequences. It produces an optimal alignment between two sequences with respect to a given scoring function. One can perform multiple sequences alignment progressively based on this DP pairwise comparison method. Roughly speaking, in the i th round of the algorithm, the alignment of the first i – 1 sequences is represented in a consensus sequence. The algorithm then updates this consensus by aligning it with the i th sequence S _i using the standard DP algorithm. Information from S _i that is not aligned to the consensus sequence isessentially lost. See Figure 1(a).

The partial order alignment (POA) algorithm [3] alleviates this problem by encoding the consensus in a POG instead of a linear sequence (see Figure 1(b)). That is, the alignment of S₁,..., S_i-1is encoded in a partial order graph G _i , which is then updated to G_i+1by aligning it with S _i . Due to the partial order in a POG, the alignment between G _i and S _i can still be achieved by a DP algorithm. The POA algorithm reduces the influence of the order of the sequences aligned, and is able to capture alignments between a subset of sequences. More details of the POA algorithm and its variants can be found in [2, 3].

In our case, each trajectory is mapped to an ordered sequence of points (i.e, a polygonal curve), and a similar algorithm can be applied to our trajectory data: Instead of the usual 1D sequences, we now have dD sequences, where d is the dimension of each point. Note that since each point corresponds to the distance map of a conformation, no transformation is needed when comparing such curves. The first stage of our EPO algorithm constructs a POG G with respect to the input set of trajectories $\mathcal{T}$ using a modified POA algorithm. Below we explain the main differences between them.

Scoring function

The choice of the scoring function when aligning G _i = (V _i , E _i ) with T _i , is in general a crucial aspect of an alignment algorithm. A good scoring function will align as many points as possible globally. Given a point p ∈ T _i and a node o ∈ G _i , let δ(o, p) be the similarity between p and o, the definition of which will be described shortly. The score of aligning p with o is usually defined as:

$\begin{array}{l}\text{Score}(\text{o},\text{p})=\\ \begin{array}{ll}\max \hfill & \{\underset{(o^{\prime },o)\in E_i}{\max }(\text{Score}(\text{o}\prime ,\text{q})+\delta (\text{o},\text{p})),\hfill \\ \begin{array}{cc}\underset{(o^{\prime },o)\in E_i}{\max }\text{Score}(\text{o}\prime ,\text{p}),& \text{Score}(\text{o},\text{q})\},\end{array}\hfill \end{array}\end{array}$

where q is the parent of the point p along T _i , and o' ranges over all immediate predecessors of o in the POG G _i . It is easy to verify that such scores can be computed by a dynamic programming procedure due to the inherent order existing in both the trajectory and the POG.

A common way to define δ(o, p), the similarity between o and p, is as follows. Assume that each node o is associated with a node center ω (o) to represent all the points aligned to this node. Then

$\delta (o,p)=\{\begin{array}{ll}\epsilon -\left|\right|p-\omega (o)\left|\right|\hfill & \text{if}\left|\right|p-\omega (o)\left|\right|<\epsilon \hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$

(1)

An alternative way to view this is that each node o has an influence region of radius ε around its center. A point p can be aligned to a node o only if it lies within the influence region of o.

Natural choices for the node center ω(o) of o include using an earlier computed canonical cluster center, or the center of the minimum enclosing ball of points already aligned to this node (or some weighted variants of it), which is a dynamic point relying on alignment order. The advantage of the former is that canonical cluster centers tend to spread apart, which helps to increase coverage of aligned nodes. Furthermore, the canonical cluster centers serve as good candidates for node centers as we already know that there are many points around them. The disadvantage is that it does not consider the distribution of points already aligned to this node. See Figure 7, where without considering the distribution of points aligned to o _a and o _b , the new point p will be aligned to o _b even though o _a is a better choice. Using the center of the minimum enclosing ball alleviates this problem. However, the influence regions of nodes produced this way tend to overlap much more than using the canonical cluster centers and the position of these centers also depend heavily on the order of curves aligned. We combine the advantages of both approaches into the following two-level scoring function for measuring the similarity δ(o, p).

Specifically, for a node o, let q be the first point aligned to this node. This means that at the time we were examining q, q cannot be aligned to any existing node in the POG. Let c _k ∈ C be the nearest canonical cluster center of q – recall that the node o was created because ||q – c _k || ≤ ε. We add c _k as a point aligned to this node, and at any time, the center of the minimum enclosing ball of currently aligned points, including c _k , will be used as the node center ω(o). Now let

$D(o)=\underset{q,q^{\prime }\in o}{\max }\left|\right|q-q^{\prime }\left|\right|$

be the diameter of points currently aligned to o. We define that:

$\delta (o,p)=\{\begin{array}{ll}2\epsilon \hfill & \text{if}\left|\right|p-\omega (o)\left|\right|<D(o)\hfill \\ \epsilon \hfill & \text{else if}\left|\right|p-\omega (o)\left|\right|<\epsilon \hfill \\ 0\hfill & \text{else}\hfill \end{array}$

(2)

In other words, the new scoring function prefers centering points to be around previously computed cluster centers, thus tending to reduce overlaps between the influence regions of different nodes. Furthermore, it gives higher similarity score for points that are more tightly grouped together with those already aligned at current node, addressing the problem shown in Figure 7. Our experimental tests have shown that this two level scoring function significantly outperforms the ones using either only the canonical centers or only the centers of minimal enclosing balls. We remark that it is possible to use variants of the above two-level scoring function, such as making it continuous (instead of being a step function). We choose the current form for its simplicity. Furthermore, experiments show that there is only marginal difference if we use the continuous version.

Merging stage

In the first stage, we have applied a progressive method to align each trajectory onto an alignment graph one by one. In the i th iteration, a point from T _i is either aligned to the best matched node in the current POG G _i , or a new node is created containing this point and the corresponding canonical cluster center, or discarded. After processing all of the N trajectories in order, we return a POG G = G _N . In the second stage of our EPO algorithm, we further improve the quality of the alignment in G by using a novel merging process.

Given the greedy nature of the POA algorithm, the alignment obtained in G is not optimal and depends on the alignment order. Furthermore, since the influence regions of different aligned-nodes may overlap, no matter how we improve the scoring function, sometimes it is simply ambiguous to decide locally where to align a new in-coming point, and a wrong decision may have grave consequence later.

For example, see Figure 8, where the set of points P (enclosed in the dotted circle) should have been aligned to one node. However, suppose the nodes o _a and o _b already exist before any point in P is inserted.

Then as points from P come in, it is rather likely that they are distributed evenly into both o _a and o _b . This problem becomes much more severe in higher dimensions, where P can be distributed to several nodes whose centers are well-separated around P, but whose influence regions still cover some points from P (the number of such regions grows exponentially w.r.t. the dimension d). Hence instead of being captured in one heavily aligned node, P is broken into many nodes with small size. Our experimental tests confirm that this is happening rather commonly in both standard and our modified POA algorithms.

To address this problem, we propose a novel postprocessing on G. The goal is to merge qualified points from neighboring less-aligned nodes to augment more heavily loaded nodes. In particular, the following two invariants are maintained during the merging process:

(I1) At any time, the diameter of the target node is still bounded by the distance threshold ε;

(I2) The partial order constraints induced by the POG graph are always consistent with the order of points along each trajectory.

The second criterion means that at any time in the POG graph G', if p ∈ o₁, q ∈ o₂, p, q ∈ T _i and p precedes q along the trajectory T _i , then either o₁ $\prec$o₂, or there is no partial order relation between them. In other words, the resulting POG still corresponds to a valid alignment of $\mathcal{T}$ with respect to the same thresholds.

As an example, see Figure 6, where the point $p_2^1$ (i.e, the second point of the trajectory T₁) in the node b in (b) is moved to the node e in (c). Note that the graph is also updated to reflect the change (the dashed edge in (c)), in order to maintain the invariants (I1) and (I2). When all points aligned to a specific node o are merged (thus moved) to other nodes (i.e, o becomes empty), we delete o, and its successors in the POG will then become the successors of its parent.

A high level pseudocode of the merging process is shown in Figure 9. It augments better aligned nodes from the current POG G by processing first the nodes with larger size. We perform this procedure a few times till there is no significant increase in the quality of the resulting alignment. In practice, to speed up the algorithm, we merge neighbors to a node o only if its size is greater than some threshold (fixed at half of the size threshold, i.e, τ/2, in our experiments), as otherwise, there is low probability that o will become a heavy node later.