BWM*: A Novel,Provable,Ensemble-based Dynamic Programming Algorithm for Sparse Approximations of Computational Protein Design

1.1. Design with sparse energy functions

Because protein design is computationally expensive, many protein design algorithms use sparse energy functions that omit interaction energy between sufficiently distant atoms (Desmet et al., 2002; Fleishman et al., 2011; Zhang and Lange, 2013; Robertson and Varani, 2007; Privett et al., 2012; Leaver-Fay et al., 2011; Lazaridis and Karplus, 1999; Jiang et al., 2000; Jones, 1994; Kaufmann et al., 2010; Kilambi and Gray, 2012; King et al., 2014; Kingsford et al., 2005; Koehl and Delarue, 1994; Kortemme et al., 2003; Krivov et al.; 2009). These sparse energy functions not only reduce the time to compute conformational energy, but also define a different energy landscape: the omitted terms eliminate energy differences between conformations, introducing optimal substructure to the energy landscape. Dynamic programming algorithms use concepts such as tree decomposition and tree width to exploit this optimal substructure to compute the corresponding GMEC more efficiently (Leaver-Fay et al., 2005; Xu and Berger, 2006; Krivov et al., 2009). However, these algorithms compute only the GMEC and do not enumerate a gap-free list of conformations. Naïve extensions to do this are worst-case exponential time in enumerating additional conformations.

We propose a novel, dynamic programming algorithm called Branch-Width Minimization* (BWM*) to efficiently and provably enumerate a gap-free ensemble of conformations, in order of increasing sparse energy. To exploit the optimal substructure, our algorithm uses the concepts of branch-decomposition with recursive heaps. Like tree decompositions, branch-decompositions have also been used in dynamic programming approaches for discrete optimization problems (Fomin, 2003; Hicks et al., 2005; Hlinẽný et al., 2008). Our algorithm treats the protein backbone as rigid and models side-chain flexibility using frequently observed, low-energy discrete conformations called rotamers (Lovell et al., 2000; Donald, 2011). Given a branch-decomposition of branch-width w for an n-residue design with at most q rotamers per residue, our algorithm computes the corresponding GMEC, called the sparse GMEC, in O( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$nw^2 q^ { \frac { 3 } { 2 } w } $$ \end{document} ) time and O( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$nq^ { \frac { 3 } { 2 } w } $$ \end{document} ) space, and enumerates each additional conformation in merely O(n log q) time and O(n) space. Because BWM* enumerates conformations in order of increasing sparse energy, the sparse ensemble of all conformations within an energy window E_w of the sparse GMEC may contain different conformations from those in the full ensemble of all conformations within E_w of the GMEC. For a given sparse energy function, we can compute a bound ɛ_B on the difference in energy between the sparse GMEC and GMEC. Using this bound, we prove that a sparse ensemble of all conformations within an expanded energy window E_w + ɛ_B of the sparse GMEC contains all conformations within E_w of the GMEC. Thus, BWM* is guaranteed to compute a sparse ensemble which contains the full ensemble.

Since BWM* is a dynamic programming algorithm, we should be able to distinguish between problems for which BWM* is more efficient, and those for which A* is more efficient. To test this hypothesis, we defined a new measure, Total Effective Search Space (TESS), which can be computed in polynomial time a priori before BWM* is run. Next, we ran BWM* and A* on 67 different protein design problems and found that TESS perfectly discriminates between BWM*-efficient cases and A*-efficient cases. For the 73% of cases in which BWM* is predicted by TESS to outperform A*, it computes the full ensemble or a close approximation faster than A*, enumerating each additional conformation in milliseconds, and computes the full ensemble in seconds, up to 315 times faster than A*. Of the other 27%, BWM* could enumerate the sparse ensemble for some problems by using a smaller energy window and sparser energy functions, while other problems were, as predicted, more suitable for a GMEC-based method or A*. In practice, we found that branch-width w (used to calculate TESS) can be small irrespective of n, making TESS much smaller than the worst-case bounds of O(qⁿ) for A*.

This article makes the following contributions:

1. A new dynamic programming algorithm called BWM*, which exploits branch-decompositions for protein designs with sparse energy functions, and proof of its asymptotic time and space complexity bounds.

2.1. Sparse residue interaction graphs

Let G = (V, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal E}$$ \end{document} ) be a residue interaction graph (Fig. 1B) corresponding to a protein design problem (Fig. 1A), with a vertex for every mutable residue, and an edge e for every pairwise residue interaction. The energy of a conformation c can be computed as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$E ( c ) = \mathop \sum \limits_{e \in {\cal E}} E_e ( c )$$ \end{document} where E_e(c) denotes the interaction energy of the residue pair in c represented by e. By omitting interaction energies between certain residue pairs, the edge set \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal E}$$ \end{document} ′ can be deleted from G, producing a sparse graph G′ = (V, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal E}$$ \end{document} − \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${\cal E}$$ \end{document} ′), shown in Fig. 1C. The energy function E′(c) corresponding to G′ can then be defined to be \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$E^{ \prime} ( c ) = \mathop \sum \limits_{e \in {\cal E}} E_e ( c ) - \mathop \sum \limits_{e \in {\cal E}^{ \prime}} E_e ( c )$$ \end{document} . In this way, the sparse energy of c is represented as the difference between the full energy E(c) and the energy terms missing from E′(c).

OPEN IN VIEWER

FIG. 1.

(A) A sample protein design problem represented as a residue interaction graph (B), with residues as vertices and pairwise interactions as edges. (C) The sparse residue interaction graph generated by deleting (a, d), (a, e), (b, d), and (b, e), shown as red crosses in (A) and (B).

2.2. Branch-decomposition in protein design

Let the branch-decomposition of G′ be an unrooted binary tree T = (D, P) with tree nodes D and tree edges P, such as the one in Fig. 2A. T satisfies the following property: For every edge \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$(v_i, v_j) \in{\cal E}-{\cal E}^\prime$$ \end{document} , there exists a leaf node in T that corresponds to (v_i, v_j). This can be seen explicitly in the rooted branch-decomposition, shown in Fig. 2B. For every edge p ∈P, we arbitrarily define the two subtrees separated by removing p as S_L and S_R, and the respective sparse graph vertices contained in each subtree as S_L(p) and S_R(p). This separation divides the vertex set V into three sets L(p), M(p), and R(p): M(p) = S_R(p) ∩ S_L(p), L(p) = S_L(p) − M(p), and R(p) = S_R(p) − M(p). By definition, this means that there are no edges between L(p) and R(p) in G′. Finally, the branch-width w of a branch-decomposition is defined to be \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$w = \max \limits_{p \in P} \mid M ( p ) \mid$$ \end{document} . Branch-decompositions and similar concepts have previously been used in dynamic programming algorithms for other discrete optimization problems (Fomin, 2003; Hicks et al., 2005; Hlinẽný et al., 2008; Nilsson, 1998; Dechter and Mateescu, 2007). The following section describes how BWM* uses the branch-decomposition T for protein design with sparse energy functions.

OPEN IN VIEWER

FIG. 2.

An example branch-decomposition. (A) The edges of the sparse graph correspond to a node in the branch-decomposition tree. Along the highlighted edge of the branch-decomposition, the mutable residues are separated into three sets: the L-set, which exists only in leaves to the left of the edge; the R-set, which exists only in leaves to the right of that edge; and the M-set, which can be found on both sides of the edge. (B) This tree is arbitrarily rooted for use by BWM*.

3.1. Total Effective Search Space

We can now define a new measure Total Effective Search Space (TESS), which predicts BWM* performance. Let the set of mutable residues be R, and the number of unpruned rotamers for a residue r ∈R be q_r. For each nonempty λ-set λ and its corresponding M-set M_λ, M_λ∪λ ⊆ R and the number of all unpruned conformations for a particular subproblem can be computed as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\prod \limits_{r \in M_{ \lambda} \cup \lambda} q_r$$ \end{document} . This value corresponds to the total number of conformations enumerated by BWM* for each M_λ ∪λ, and summing over all nonempty λ-sets is a deterministic measure of the time and space complexity of BWM*. We define the sum of products \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\mathop \sum \limits_{ \lambda} \prod \limits_{r \in M_{ \lambda} \cup \lambda} q_r$$ \end{document} as the Total Effective Search Space. Because computing the smallest-width branch-decomposition is NP-complete, polynomial-time approximation algorithms that typically do not return the optimal branch-decomposition are used in practice. We now show that when using an algorithm that provably computes a branch-decomposition in polynomial time β_t = n^O(1), TESS can be computed in polynomial time. The proof of Theorem 1 is provided in section A.1 of the supplementary information (SI) in Jou et al. (2015).

Theorem 1. TESS can be computed in O(n² + β_t) time, where β_t is the time taken to compute the branch-decomposition of a sparse graph.

3.2. Algorithm: preprocessing and enumeration

The preprocessing phase of BWM* computes the energy of the sparse GMEC and constructs a data structure for efficient enumeration of conformations in order of increasing sparse energy. This data structure is called a recursive heap. Let a recursive heap H be a canonical min heap satisfying the heap property, with the following two additional properties:

1. For every heap node in H, there are zero, one, or two child heaps. These child heaps are also recurisve heaps, and their heap nodes have the same properties.

During preprocessing, BWM* performs a post-order traversal of T. For each edge with a nonempty λ-set, the following operations are performed:

1. Enumerate all possible rotamer assignments to the residues in M ∪ λ, and for each assignment look up optimal assignments for the residues in L − λ, which were previously computed in its child edges.

For every assignment in the M-set, we enumerate all assignments of the λ-set and look up all remaining assignments to the (L − λ)-set from the child edges. Therefore, at the end of this procedure we have the energy of the optimal assignment to the L-set for each assignment to the M-set in each recursive heap. Since the residues of the L-set interact only with each other and the residues of the M-set, the optimal assignment to the L-set is determined at the end of this procedure. As the L-set of the root edge contains all n mutable residues, once the traversal has returned to the root edge the optimal solution for every residue has been calculated, and the energy of the sparse GMEC is contained in the root node of the recursive heap constructed at the root edge. The sparse GMEC can then be calculated by finding the lowest-energy partial conformation in the λ-set at the root, and recursively looking up the optimal assignments in its children to reassemble the corresponding full conformation. After the sparse GMEC has been returned, the heap must be updated to return the next best conformation. This procedure has two steps:

1. Call this procedure recursively on its two child heaps, and update the energy of the root node with the new energies at the roots of its two child heaps.

After these steps are completed, the resulting root of the heap now contains the energy of the next best conformation. This procedure is called repeatedly to enumerate additional conformations.

We now give an upper bound on the size of |M ∪λ|, and analysis of the time and space complexity of the preprocessing and enumeration steps of BWM*. The proofs of the theorems, along with details of the construction of recursive heaps and enumeration phase are provided in section A.2 and section A.3 of the SI in Jou et al. (2015).

Theorem 2. The maximum subproblem M ∪ λ is bounded by the relation \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\mid M \cup \lambda \mid \leq \frac { 3 } { 2 } w$$ \end{document} .

Theorem 3. BWM* takes O( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$nw^2 q^ { \frac { 3 } { 2 } w } $$ \end{document} ) time and O( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$nq^ { \frac { 3 } { 2 } w } $$ \end{document} ) space to preprocess and return the sparse GMEC.

Theorem 4. Enumerating the next best conformation takes O(n log q) time and O(n) space to remove the minimum-energy conformation from H.

3.3. Sparse error bounds

To show the bounds on error introduced by sparse graphs, we first define two terms to bound the largest possible difference between sparse and full energy for any conformation c from the set of all allowed conformations C. Let \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$E^{ \prime}_{max} = \mathop \sum \limits_{e \in {\cal E}^{ \prime}} \max \limits_{c \in C}E_e ( c )$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$E^{ \prime}_{min} = \mathop \sum \limits_{e \in {\cal E}^{ \prime}} \min \limits_{c \in C} E_e ( c )$$ \end{document} bound the positive and negative energy change for c, respectively. While these bounds are loose, they have the benefit of being computable in polynomial time. With these two terms, we can bound the total energy difference between the GMEC c* and the sparse GMEC c′, and show the relationship between the sparse ensemble and full ensemble. The following lemma bounds the sparse energy difference between the sparse GMEC and GMEC. The proof for lemma 1 is provided in section A.4 of the SI in Jou et al. (2015).

Lemma 1. The difference in sparse energy between the sparse GMEC c′ and the GMEC c* is bounded by the relationship |E′(c*) − E′(c′)| < E′ _max − E′_min.

With these results we can now show that a sparse ensemble S′ of conformations within E_w + E′ _max − E′ _min of the sparse GMEC is guaranteed to contain the full ensemble S* of all conformations within E_w of the GMEC. The proof for Theorem 5 is provided in section A.4 of the SI in Jou et al. (2015).

Theorem 5. For any conformation c, if E(c) ≤ E(c*) + E_w then E′(c) ≤ E′(c′) + E_w + E′ _max  − E′ _min

We can now show the total complexity to enumerate a gap free list of k conformations, starting from the sparse GMEC. The proof for Theorem 6 is provided in section A of the SI in Jou et al. (2015).

Theorem 6. BWM* provably computes the sparse GMEC and an ensemble of the top k conformations in O( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$nw^2 q^ { \frac { 3 } { 2 } w } + kn \log q$$ \end{document} ) time and O( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$nq^ { \frac { 3 } { 2 } w } + kn$$ \end{document} ) space, which is guaranteed to contain all conformations within E_w of the GMEC when all conformations within E_w + E′ _max − E′ _min have been enumerated.

4.1. Experimental methods

We implemented BWM* in our laboratory's open source osprey (Gainza et al., 2013) protein design package. To test the hypothesis that TESS could predict whether BWM* or A* would be faster, and also determine empirical enumeration times, we ran BWM* and A* on 67 different design problems taken from Gainza et al., (2012). Each protein design problem consists of a rigid backbone design with 4–16 mutable residues and 5–10 allowed amino acid mutations per residue. Side-chain flexibility was modeled using a rigid rotamer library with 153 modal rotamers (Lovell et al., 2000). For each protein design problem, DEE was run with energy windows E_w of both 0.5 kcal/mol and 1 kcal/mol, followed by either A* or BWM* with different distance and energy cutoffs: distance cutoffs (deleting edges whose minimum distance between all unpruned rotamer pairs exceeds the cutoff distance) of 7 Å, 8 Å, and energy cutoffs (deleting edges whose maximum absolute interaction energy between all unpruned rotamer pairs is less than cutoff energy) of 0.1 kcal/mol and 0.2 kcal/mol. Then the performance was compared with the TESS predictions. Information on the 67 protein design problems, along with details of the computational experiments performed, are provided in section C of the SI in Jou et al. (2015).

4.2. Total Effective Search Space predictions

To determine if TESS could be used to predict BWM* performance, we plotted TESS for all 67 protein design problems against the number of unpruned conformations after DEE, as shown in Fig. 3. Shown in green dots is the TESS for 39 out of 67 design problems for which BWM* enumerated the sparse ensemble for an energy cutoff of 0.1 kcal/mol and an energy window of 1 kcal/mol. Ten more design problems (shown as light green dots) were solved efficiently by increasing the energy cutoff to 0.2 kcal/mol to make the residue interaction graph sparser, and decreasing the energy window to 0.5 kcal/mol to prune additional conformations from the design search space. The TESS for these 49 problems in which BWM* outperforms A* is less than 8.97 × 10⁸ conformations (shown with a line in Fig. 3). TESS remains larger for the remaining 18 problems for which BWM* exceeded memory limits during preprocessing. Using a larger energy cutoff (0.3–0.5 kcal/mol) reduced the TESS for 6 of these 18 problems (yellow triangles) below 8.97 × 10⁸ (blue triangles), and BWM* was able to compute the sparse ensemble, as seen in Fig. 3B. For 11 additional problems (orange diamonds), reducing the energy window to 0 kcal/mol reduced their TESS below 1.57 × 10¹¹ conformations (blue diamonds), which was sufficient to compute the sparse GMEC but not the sparse ensemble. In one case neither the sparse ensemble nor the sparse GMEC could be computed, which is shown as a red dot and red cross, respectively. In this case the periplasmic copper/silver-binding protein (PDB id: 2qcp) contains a set of residues with large interaction energy greater than 0.5 kcal/mol and many unpruned rotamers, producing a dense subgraph with a large TESS.

OPEN IN VIEWER

FIG. 3.

Total Effective Search Space predicts BWM* performance. x-axis shows TESS for each design problem using their respective cutoff and energy window. y-axis shows unpruned conformations for each design problem after DEE pruning with their respective energy window. (A) Dark green points: full ensemble with energy cutoff 0.1 kcal/mol, energy window 1 kcal/mol. Light green points: full ensemble with energy cutoff 0.2 kcal/mol, energy window 0.5 kcal/mol. The yellow triangles, orange diamonds, and red dot have x values of TESS using unpruned conformations with energy cutoff 0.1 kcal/mol, energy window 1 kcal/mol for cases where larger energy cutoffs or smaller energy windows were necessary to finish preprocessing. (B) Blue triangles: x values are TESS for calculating only the sparse ensemble (energy cutoff 0.3 kcal/mol or higher, energy window 0.5 kcal/mol). (C) Blue diamonds: x values are TESS for calculating only the sparse GMEC (energy cutoff between 0.1 and 0.5 kcal/mol, energy window 0 kcal/mol). Red cross: x value is TESS for computing the sparse GMEC of 2qcp (see text).

Based on our analysis, the plot contains three distinct regions: all protein design problems with TESS below 8.97 × 10⁸ conformations are tractable for BWM*, and in this region it will always outperform A*. For protein design problems with TESS values above 1.57 × 10¹¹ conformations, BWM* attempts to enumerate all conformations of a prohibitively large subproblem even when computing only the sparse GMEC, and A* will have the advantage. For the problems that lie between these two thresholds, increasing the energy cutoff and reducing the energy window allows BWM* to return the sparse GMEC and possibly the sparse ensemble as well. This shows that the performance of BWM* can be reliably predicted with TESS. The next section compares the performance of BWM* and A* for the 49 problems on which BWM* is predicted to perform better.

4.3. Ensemble enumeration time

BWM* was able to provably compute the full ensemble in 45 out of 49 protein design problems for at least one combination of energy window, distance, and energy cutoffs. For the 36 design problems where BWM* enumerates the full ensemble with E_w = 1 kcal/mol, energy cutoff 0.1 kcal/mol, Fig. 4A shows the maximum observed problem size |M ∪ λ| and branch-width w compared to the number of mutable residues n. The maximum observed problem size never exceeds 6, and its upper bound never goes beyond 7.5 (w ≤ 5). The fact that the maximum observed subproblem is small irrespective of n makes TESS much smaller than the search space of A* (O(qⁿ)), as shown by the solid line in Fig. 4C. This reduced effective search space allows BWM* to search more efficiently than A*, and as plotted in Fig. 4B, BWM* enumerates not only the sparse ensemble but also the full ensemble faster than A*. In cases where BWM* cannot fully enumerate the full ensemble (because the error bounds calculated according to Lemma 1 are large), the sparse ensemble contains more than 93% of all conformations in the full ensemble. One such case is the protein YciI from H. influenzae (PDB id: 1mwq) with energy window of 1 kcal/mol and an energy cutoff of 0.1 kcal/mol, where the number of conformations enumerated by BWM* exceeded 10,000 before it provably enumerated the full ensemble. Fig. 5 compares the sequence-specific ensembles for this case, and shows that the sparse ensemble contains almost all low-energy conformations in the full ensemble. The missing conformations introduce less than a 3% difference in the Boltzmann-weighted partition function (Roberts et al., 2012; Georgiev et al., 2008b). A detailed description of the high similarity between sparse and full ensemble, as well as additional runtime comparisons for the remaining nine cases in which BWM* computes the full ensemble can be found in section B of the SI in Jou et al. (2015).

OPEN IN VIEWER

FIG. 4.

Both actual and worst-case bounds for BWM* are better than A*. (A) Comparison of the number of mutable residues considered. The maximum observed subproblem size is |M ∪ λ| and worst-case problem size is \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \frac { 3 } { 2 } w$$ \end{document} , which bounds |M ∪ λ| from above. (B) Times to enumerate full ensemble for all design problems with energy window of 1 kcal/mol and energy cutoff of 0.1 kcal/mol. (C) Actual search space size vs. worst-case bounds. x-Axis is the 36 design problems in (A). Problems are ranked by A* worst-case problem size (qⁿ); q and n measure the size of the design problem and are shown below the line; n increases from 4 to 16, shown in orange, and q is between 3 and 41, shown in red. Solid green line: Actual number ( \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ \prod \limits_r q_r$$ \end{document} ) of unpruned conformations after DEE, where q_r is the number of unpruned rotamers for mutable residue r and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$q = \max \limits_r q_r$$ \end{document} . Solid blue: TESS. Dotted blue: BWM worst-case problem size bounds q^|M∪λ|. Dotted green: A* worst-case problem size bounds qⁿ.

OPEN IN VIEWER

FIG. 5.

The sparse ensemble closely approximates the full ensemble. Figures shown for YciI protein from H. influenzae (PDB id: 1mwq) with energy cutoff of 0.1 kcal/mol and energy window of 1 kcal/mol. For each sequence, the x-axis plots the full energy of each conformation of the ensemble, and the y-axis plots the percent contribution of each conformation to the partition function Z (called q in Georgiev et al., 2008b). Conformations in both the full and sparse ensembles are colored green. Conformations in the full ensemble that are not in the sparse ensemble are in red. Conformations in the sparse ensemble that are not in the full ensemble are in blue. GMEC is circled in green.

Fig. 6 shows the individual and cumulative time taken by BWM* and A* to return each additional conformation in the gap-free list for two representative design problems. For the aortic preferentially expressed protein-1 (PDB id:1u2h), BWM* computes the full ensemble before A* even returns the first conformation. For cases like the third KH domain of heterogeneous nuclear ribonucleoprotein K (PDB id: 1zzk), the preprocessing time of BWM* can be significant, but once BWM* begins enumerating conformations it rapidly enumerates more conformations than A* in less time. This is because BWM* requires only O(n log q) time (milliseconds, empirically) for each additional conformation generated, as shown in Fig. 6B and D. Simply modifying A* to use the sparse energy function does not change the worst-case complexity of A*, and BWM* is still faster. A comparison of runtimes showing superior performance of BWM* compared to A* with a sparse energy function are provided as Fig. 3 in section B of the SI in Jou et al. (2015).

OPEN IN VIEWER

FIG. 6.

BWM* enumeration is orders of magnitude faster than A*. Figures shown for aortic preferentially expressed protein-1 (PDB id:1u2h) with energy cutoff 0.1 kcal/mol, and third KH domain of heterogeneous nuclear ribonucleoprotein K (PDB id: 1zzk) with distance cutoff of 7 Å. Energy window was 1 kcal/mol for both. Blue: BWM* time. Green: A* time. (A) and (C) show cumulative time against total conformations enumerated. (B) and (D) show enumeration time per conformation.