Computing H/D-Exchange rates of single residues from data of proteolytic fragments

Abstract

Background

Protein conformation and protein/protein interaction can be elucidated by solution-phase Hydrogen/Deuterium exchange (sHDX) coupled to high-resolution mass analysis of the digested protein or protein complex. In sHDX experiments mutant proteins are compared to wild-type proteins or a ligand is added to the protein and compared to the wild-type protein (or mutant). The number of deuteriums incorporated into the polypeptides generated from the protease digest of the protein is related to the solvent accessibility of amide protons within the original protein construct.

Results

In this work, sHDX data was collected on a 14.5 T FT-ICR MS. An algorithm was developed based on combinatorial optimization that predicts deuterium exchange with high spatial resolution based on the sHDX data of overlapping proteolytic fragments. Often the algorithm assigns deuterium exchange with single residue resolution.

Conclusions

With our new method it is possible to automatically determine deuterium exchange with higher spatial resolution than the level of digested fragments.

Background

In the solution-phase Hyrdogen/Deuterium ex-change (sHDX) experiment, protein surface accessibility is probed by exchange of labile hydrogen for deuterium. Simply speaking, hydrogens located at solvent exposed sites exchange at a higher rate with deuteriums from the solution than others. From these exchange rates one can therefore deduce information about protein solvent accessibility and thus protein conformation.

There is controversy surrounding the effect of D₂O solvent on the conformation of proteins. Sheu et al. [1] used molecular dynamic modeling of a small peptide to illustrate compaction of the peptide conformation in D₂O versus H₂O. This small compaction of the conformation occurs when the pep-tide is fully deuterated (which is never observed in the sHDX experiments). Since sHDX monitors the incorporation of deuterium over time the resulting slight compaction of the structure is minimized. Other methods used for the study of protein/protein interaction or protein conformation such as cross-linking [2, 3] or hydroxyl radical addition [4–6] result in large conformational change of the protein structure; leaving sHDX as the method of choice for probing protein conformational changes in solution.

NMR spectroscopy has been the gold standard for determination of protein structure, but it has limitations on protein solubility and molecular weight (<50 kD). Solution-phase HDX with mass spectrometry analysis has higher sensitivity and is not limited by molecular weight, but sHDX is hampered with a major difficulty. One only obtains exchange data for peptic fragments and assigning exchange rates to single residues has to be done by manual interpretation.

We provide an automated method to resolve this problem. More precisely, we present an algorithm that enumerates all possible exchange rates for single residues that explain the observed data of the peptic fragments. As the number of possibilities is often very large, we combine sets of assignments to equivalence classes which are easily interpreted such that the number of equivalence classes is typically very small.

The assignment of exchange rates to single residues from the data of the peptic fragments is a combinatorial problem. Hence, we apply methods from combinatorial optimization to it, i.e. we show how to formalize the problem as an integer linear program and propose methods to solve the problem.

Biochemical Background

Concerning the determination of protein-protein interaction, X-ray crystal diffraction and NMR [7] pro-vide the highest resolution of the sites of interaction. On the downside, both methods require large (milligram) quantities of protein. Other techniques rely on chemical or photo-induced reactions with MS analysis [8, 9] to reveal functional groups that are ex-posed to the solvent. These methods also suffer from physical limitations.

Another method utilizes hydroxyl radical reactions with alkyl CH bonds. The OH tends to re-act mainly with surface-exposed residues providing a good footprint of the solvent exposed surface of the protein(s) [4, 6]. The modification is covalent and thus irreversible, but each modification can potentially change the conformation of the protein, thus skewing results.

Exchange of labile hydrogens for deuteriums (sHDX) as a probe of protein surface accessibility does not change the conformation of the protein. Advantages of MS over NMR and X-ray crystallography structural determination are the ability to work at low concentration and high molecular weight.

The experiment is initiated by dilution of the protein solution into a biological buffer made with D₂O. Solvent accessible hydrogens are exchanged with deuterium. The exchange is quenched (greatly slowed) by dropping the pH to between pH 2.3 and pH 2.5 and lowering the temperature to approximately 0°C. The protein complex is digested with a protease that is active under quench conditions (such as pepsin) and on-line liquid chromatography is performed directly to the FT-ICR MS. Deuterium in corporation is monitored by the increase in mass of each peptic fragment as the deuterons are added.

These data sets are large, often with many over-lapping proteolytic fragments. From these data, the exchange rate is easily determined for the same peptic fragments from the protein and the protein/protein complex [10] (all other fragments are disregarded). When peptic fragments are not directly comparable, but are overlapping (Figure 1) manual interpretation must be performed to assign exchange rate to single residues. Data analysis is the greatest bottleneck in sHDX experiments; thus automated data analysis is necessary. Furthermore, we are interested in all such assignments, as averaging over all solutions gives better results in practice.

Figure 1

Sample data set. Overlap of peptic fragments obtained from our sHDX on myoglobin. The table on the right shows the number of amide hydrogens predicted to be either slow, medium or fast (based on MEM). The vertical lines show the decomposition of the sequence into parts induced by the fragments. We want to automatically draw conclusions on the exchange rates of single amino acids, like the one that the second D residue has to have medium exchange rate, concluded from the restrictions imposed by fragments 3 and 5. Notice that we already deleted the N-terminal amino acid as we cannot see them exchange.

Mathematical Abstraction

In this section, we present our mathematical model for the assignment of exchange rates to residues. A brief overview of the introduced terms and symbols can be found in Table 1. In an idealized setting, we consider the following problem. We sequentially number the n residues of a protein from1 to n, beginning at the N-terminal residue and ending at the C-terminal residue. The set of peptic fragments resulting from the digestion of the protein is captured by a set ℱ of integer intervals (i, j) := {k ∈ ℕ | i ≤ k ≤ j}, for two positive integers i, j with i ≤ j, representing the endpoints of the corresponding fragment. In other words, the peptic fragment represented by (i, j) spans residues i, i+1, ..., j. Furthermore, K denotes the number of different classes of exchange rates, arising from the discretization of the experimentally measured deuterium uptake rates [11]. The K distinct classes of exchange rates, to which we simply refer as colors in the following, are represented by set \mathcal{S}. To simplify notation we number the colors from 1 to K and identify in the following the colors by their respective number. The experimentally found bulk information of how many residues within each fragment (i, j) ∈ ℱ fall into each of the exchange rate categories is given by "requirement" integers b_{\left(i,j\right)}^k, for each fragment (i, j) ∈ ℱ and each color k ∈ \mathcal{S}. We call the vector b^kof "requirements" with respect to color k, indexed by fragments from, ℱ the right hand side for color k. In our experimental data, exactly three different colors are distinguished (interpreted as slow, medium, and fast exchange rates), i.e. K = 3. However, our method is not restricted to this case.

Table 1 Overview of terms and symbols

The mathematical notion introduced above is illustrated in Figure 1. There the residues, numbered from 1 to 28, are spanned by 9 peptic fragments, i.e. |ℱ = 9|. The third peptic fragment "VWGKVEAD" will then be represented by the integer interval (12; 19). From the experimental data we know that 5 out of the 8 residues contained in this fragment exchanged slowly (s), two at medium rate (m), and the last remaining residue exchanged fast (f) and thus b_{(12,19)}^s=5,b_{(12,19)}^m=2, and b_{(12,19)}^f=1.

Determining the exchange rate of single residues from the experimentally found data for the peptic fragments then translates into finding a "consistent" assignment of colors from \mathcal{S} to the integer points from {1, ..., n}, representing the residues of the protein, that complies with the constraints imposed by the "requirements" b_{\left(i,j\right)}^k. More precisely, we have to determine an assignment π : {1, ..., n} ↦ \mathcal{S} such that |{i ≤ l ≤ j : π(l) = k}| = b_{\left(i,j\right)}^k for all given fragments (i, j) ∈ ℱ and all possible colors k ∈ \mathcal{S}. We call such an assignment feasible.

We say that two fragments (i, j) and (i', j') over-lap, if they share at least one common residue, i.e. (i, j) ∩ (i', j') ≠ ∅. The partition of the set of fragments ℱ into a maximum number of subsets, such that no two fragments from different subsets overlap, defines independent subproblems; an assignment of exchange rates to the residues spanned by the fragments of one subset does not affect the solution of a subproblem corresponding to any other subset of fragments.

Furthermore, we denote by \mathcal{P} the partition of the set of residues {1, ..., n} into maximal subsets such that residues from the same subset are spanned by exactly the same set of fragments. More precisely, for all residues i and j in the same part of \mathcal{P} and for all fragments f ∈ ℱ it holds i ∈ f ⇔ j ∈ f. Hence, for each part p ∈ \mathcal{P} and each fragment f ∈ ℱ either p ⊆ f or p ∩ f = ∅. In Figure 1 for example, residues number 7 (Q), 8 (Q), and 9 (V), are all contained in fragments number 1, 2, 6, and 8 and thus form an element \stackrel{ª}{p} of partition \mathcal{P}. Note that for the two neighboring residues the set of containing fragments differs from {1, 2, 6, 8} and therefore part\stackrel{ª}{p} = {7, 8, 9} is maximal.

However, data collected in real experiments usually contain some noise, such that no feasible assignment of exchange rates as defined above exists. Therefore, the goal is to compute all assignments that minimize the total sum of errors. Here, the error of an assignment π in fragment (i, j) ∈ ℱ with respect to color k is defined as the surplus, respectively the shortage, of residues in (i, j) that are assigned color k, compared to the number of such residues suggested by the experimental data. That is,

e_{(i,j)}^k=|b_{(i,j)}^k-|\{i\le l\le j:\pi (l)=k\}\left|\right|

(1)

and thus the objective is to minimize the sum of this deviations over all fragments and colors, i.e.

\text{minimize}{\displaystyle \sum _{k\in \mathcal{S}}{\displaystyle \sum _{\left(i,j\right)\in ℱ}{e}_{\left(i,j\right)}^k}.}

(2)

In Figure 1 the colors green, yellow and red encode an optimal assignment π* of the exchange rates slow, medium and fast, with respect to objective (2). Under this assignment, fragment 3 contributes an error of 1 both w.r.t. color yellow (medium exchange rate)and color red (fast rate) to the total error of 17, while it satisfies the requirement for color green (slow rate, numbered 1) b_{\left(12,19\right)}^1=5 exactly.

Results and Discussion

In the following, we present different approaches to tackle the assingment problem that we have derrived from the mathematical abstraction mentioned before.

Integer linear programming formulation

First, we formulate the idealized version of the problem assuming error-free experimental data as an integer linear program (ILP). That is, we give an ILP whose feasible solutions correspond one-to-one to the feasible assignments of colors to residues.

Let π : {1, ..., n} ↦ \mathcal{S} be an assignment of colors to residues. A binary variable x_i^k for every color k ∈\mathcal{S} and every residue i ∈ {1, ..., n} indicates whether residue i is assigned color k or not, i.e.

x_i^k=\{\begin{array}{l}1\text{if}\pi (i)=k\hfill \\ 0\text{otherwise}.\hfill \end{array}

We denote by x^k=\left(x_1^k,x_2^k,\mathrm{...},x_n^k\right) the vector of binary variables modeling the assignment of color k and let x={(x^k)}_{k\in \mathcal{S}},x\in {\{0,1\}}^{Kn}.

Since every residue is assigned exactly one color, it must hold {\displaystyle {\sum }_{k\in \mathcal{S}}{x}_i^k}=1 for all i ∈ {1, ..., n}. Conversely, every 0-1 assignment to variables x satisfying {\displaystyle {\sum }_{k\in \mathcal{S}}{x}_i^k}=1 for all i ∈ {1, ..., n} corresponds to an assignment of colors to residues. A 0-1 assignment to x corresponds to a feasible color assignment π, if and only if furthermore {\displaystyle {\sum }_{l=i}^j{x}_l^k}=b_{\left(i,j\right)}^k holds for all (i, j) ∈ ℱ and k ∈ \mathcal{S}.

Now consider the problem of computing an assignment with minimum total error. Translating the definition of the error that we make when assigning color k (or not) to residues in fragment (i, j) (see equation (1)) to the context of 0-1 assignments to variables x^k, the problem of minimizing (2) becomes

\text{minimize}{\displaystyle \sum _{k\in \mathcal{S}}{\displaystyle \sum _{\left(i,j\right)\in ℱ}|b_{\left(i,j\right)}^k-{\displaystyle \sum _{l=i}^j{x}_l^k}|.}}

Concerning the formulation of a minimum sum of absolute values in terms of a linear objective function and linear constraints, observe that |b_{\left(i,j\right)}^k-{\displaystyle {\sum }_{l=i}^j{x}_l^k}| is the smallest number e_{\left(i,j\right)}^k that satisfies

\begin{array}{c}{e}_{\left(i,j\right)}^k\ge {\displaystyle \sum _{l=i}^j{x}_l^k}-b_{\left(i,j\right)}^k\text{and}\\ e_{\left(i,j\right)}^k\ge -{\displaystyle \sum _{l=i}^j{x}_l^k}+b_{\left(i,j\right)}^k\end{array}

Hence, after introducing a variable e_{\left(i,j\right)}^k for every color k ∈ \mathcal{S} and every fragment (i, j) ∈ ℱ, the integer linear program we are looking at is

\begin{array}{ll}\min \hfill & {\displaystyle \sum _{k\in \mathcal{S}}{\displaystyle \sum _{(i,j)\in ℱ}{e}_{(i,j)}^k}}\hfill \\ \text{s}.\text{t}.\hfill & e_{(i,j)}^k\ge {\displaystyle \sum _{l=i}^j{x}_l^k-b_{(i,j)}^k}\hfill & \text{for}\text{all}k\in \mathcal{S},(i,J)\in ℱ\hfill \\ e_{(i,j)}^k\ge -{\displaystyle \sum _{l=i}^j{x}_l^k+b_{(i,j)}^k}\hfill & \text{for}\text{all}k\in \mathcal{S},(i,J)\in ℱ\hfill \\ {\displaystyle \sum _{k\in \mathcal{S}}{x}_l^k}=1\hfill & \text{for}\text{all}1\le l\le n\hfill \\ x\in {\{0,1\}}^{Kn}\hfill \end{array}

We refer to this integer linear program as basic-ILP.

In our experiments, it turns out that finding a single solution is very fast, whereas enumerating all solutions takes quite some time due to their large number. This large number can be explained as follows: Recall that \mathcal{P} is the partition of {1, ..., n} into a minimal number of parts, such that for each element p ∈ \mathcal{P} and each fragment f ∈ F either p ⊆ f or p ∩ f = ∅. In other words, no fragment starts or ends within such a part. Therefore, from an assignment π we can derive further assignments π' exhibiting the same total error, by simply permuting the colors within these parts, i.e. if i, j ∈ p for p ∈ \mathcal{P} and the total error of an assignment π is e₁, than π' with π' (i) = π (j), π' (j) = π (i) and π' (l) = π (l) for l ≠ i, j has total error e₂ with e₂ = e₁. We call two assignments equivalent, if one can be obtained from the other by iteratively applying this rule.

In order to enumerate equivalent solutions only once, we modify our integer linear program as follows: For k ∈ \mathcal{S} and p ∈ \mathcal{P}, we replace the binary variables {(x_l^k)}_{l\in p} by a single integer variable y_p^k with y_p^k:={\displaystyle {\sum }_{l\in p}{x}_l^k}. Moreover, let A be the |ℱ| × |\mathcal{P}| inclusion matrix, i.e. for every f ∈ ℱ and p ∈ \mathcal{P}, the corresponding entry is given by

a_{f,p}=\{\begin{array}{cc}1& \text{if }p\subseteq f\\ 0& \text{otherwise}.\end{array}

We denote by e^k={(e_{(i,j)}^k)}_{(i,j)\in ℱ} the vector of errors with respect to color k and by y^k={(y_p^k)}_{p\in } the number of residues colored k. In matrix notation the constraints are then of the form

\begin{array}{c}-A{y}^k+e^k\ge -b^k\\ A{y}^k+e^k\ge b^k\end{array}

for all k ∈ \mathcal{S}. Hence our integer linear program be-comes

\begin{array}{ll}\min \hfill & {\displaystyle \sum _{k\in \mathcal{S}}{\displaystyle \sum _{f\in ℱ}{e}_f^k}}\hfill \\ \text{s}.\text{t}\text{.}\hfill & -A{y}^k+e^k\ge -b^k\hfill & \text{for}\text{all}k\in \mathcal{S}\hfill \\ A{y}^k+e^k\ge b^k\hfill & \text{for}\text{all}k\in \mathcal{S}\hfill \\ {\displaystyle \sum _{k\in \mathcal{S}}{y}^k=P}\hfill \\ y\ge 0,\text{integer}\hfill \end{array}

(3)

where P is the vector that contains |p| for each component p ∈ \mathcal{P}and y={(y^k)}_{k\in \mathcal{S}} We refer to this integer linear program as improved-ILP. We compute all solutions within a certain error bound by following basically the same approach as described above. However, the number of solutions now is just a fraction of the number of solutions of the original basic-ILP yielding a significant speed-up

Although there is commercial software for integer programming which quickly solves instances of reasonable size, there is no algorithm that is guaranteed to find an optimum solution in polynomial time, since integer programming is NP-complete in general. However, the problem of assigning exchange rates to residues in a way that is conform with the experimentally found bulk data exhibits a certain combinatorial structure. In the next section, we exploit this fact to derive an exact polynomial-time algorithm for the case of two colors and use it as a building block for approximation algorithms for more than two colors subsequently.

A Combinatorial Approach

First, let us consider the special case of two colors, i.e. K = 2 and thus \mathcal{S} = {1,2}. That is, we have constraints of the form y_p^1+y_p^2=\left|p\right| for all p ∈ \mathcal{P}. This allows us to simplify the linear program considerably. We replace y_p^2=\left|p\right|-y_p^1 and omit the superscript of the y-variables in the following. This yields

\begin{array}{cc}-Ay+e^1\ge -b^1& Ay+e^2\ge F-b^2\\ Ay+e^1\ge b^1& -Ay+e^2\ge -F+b^2\end{array}

where F is the vector of fragment sizes. We may get rid of half of the constraints by the following observation. Let b := max {b¹, F - b²} and \stackrel{ª}{b}:=\min \{b^1,F-b^2\} where the maximum is taken component-wise. Let y be an arbitrary feasible solution with minimum total error {\sum }_{f\in ℱ}{e}_f^1+e_f^2. We may consider the contribution of each fragment independently for that particular y. We may rename the error variables e¹ and e² component-wise according to b and \stackrel{ª}{b}, i.e.

\begin{array}{cc}{e}_f:=\{\begin{array}{ll}{e}_f^1\hfill & \text{if}{b}_f=b_f^1\hfill \\ e_f^2\hfill & \text{otherwise}\hfill \end{array}& {\overline{e}}_f:=\{\begin{array}{ll}{e}_f^1\hfill & \text{if}{\overline{b}}_f=b_f^1\hfill \\ e_f^2\hfill & \text{otherwise}\hfill \end{array}\end{array}

(4)

For each f ∈ ℱ with {\stackrel{ª}{b}}_f\le a_f^{T}y\le b_f, we have e_f^1+e_f^2=b_f-{\stackrel{ª}{b}}_f. If a_f^{T}y>b_f, we get e_f^1+e_f^2=2e_f+b_f-{\stackrel{ª}{b}}_f. Analogously, we get e_f^1+e_f^2=2{\stackrel{ª}{e}}_f+b_f-{\stackrel{ª}{b}}_f if . Hence, it is sufficient to optimize the following linear program

\begin{array}{ll}\min \hfill & {\displaystyle \sum _{f\in ℱ}{e}_f+{\overline{e}}_f}\hfill \\ \text{s}.\text{t}.\hfill & -Ay+e\ge -b\hfill \\ Ay+\overline{e}\ge \overline{b}\hfill \\ -y\ge -P\hfill \\ y,e,\overline{e}\ge 0\hfill \end{array}

(5)

which is integral if b and \stackrel{ª}{b} are integral since the constraint matrix is totally unimodular. The corresponding dual LP is given by

\begin{array}{ll}\max \hfill & -b^T{f}^1-{\overline{b}}^T{f}^2+P^T{f}^3\hfill \\ \text{s}.\text{t}.\hfill & -A^T{f}^1+A^T{f}^2-f^3\le 0\hfill \\ 0\le f^{1,2}\le 1\hfill \\ 0\le f^3\hfill \end{array}

(6)

which is equivalent to (multiplying the objective function by -1 and introducing slack variables)

\begin{array}{ll}-\min \hfill & -b^T{f}^1-{\overline{b}}^T{f}^2+P^T{f}^3\hfill \\ \text{s}.\text{t}.\hfill & -A^T{f}^1+A^T{f}^2-f^3+f^4=0\hfill \\ 0\le f^{1,2}\le 1\hfill \\ 0\le f^{3,4}\hfill \end{array}

(7)

We will show next that this LP is a Minimum Cost Circulation Problem. To this end, let M be the matrix of the equality constraints, i.e.

M:=(-A^T{A}^T-II).

Note that this matrix has the column-wise consecutive-ones property. By row operations like in Gaussian elimination, we can easily transform M such that each column contains exactly one +1 and one -1, as follows. We add the dummy constraint 0 = 0 at the end and subtract from each row its predecessor. The resulting matrix, say \stackrel{ª}{M}, can be considered as the node-arc-incidence matrix of a directed graph. Since the right hand side remains unchanged, we get a Minimum Cost Circulation problem on a graph with |\mathcal{P}| + 1 nodes and O(||+|ℱ|) arcs [12]. As a matter of fact, we have for each variable y _p two arcs corresponding to the constraint 0 ≤ y _p ≤ |p| and for each fragment (i, j) the arcs (i, j + 1) and (j + 1, i) as depicted in Figure 2.

Figure 2

Fragment graph example. Example of a fragment graph with |\mathcal{P}| = 7 The corresponding fragments are (1, 3), (2, 5), (3, 6), and (5, 7).

For three or more colors the complexity is open. The totally unimodularity of the constraint matrix is destroyed, i.e. there are instances with fractional vertices, e.g. the one from Figure 2 with the appropriate right hand sides. Moreover, there is an instance which has a positive error, but the value of the LP is 0. Hence the integrality gap is infinite. If the number of colors is not fixed but part of the input, the problem is NP-complete [13].

A Simple and Efficient Heuristic for the General Case

We present an algorithm that uses our combinatorial approach for the 2-color case (K = 2) from previous section as a subroutine to provide solutions that approximate (without performance guarantee) a coloring, i.e. an assignment of colors to residues, with minimum total error for instances with arbitrary but fixed number of colors. The general idea is to reduce the problem to the 2-color case by merging all but one color, say color i, to a single color and solve the resulting problem by an algorithm for the minimum cost circulation problem, as described in the section about the Combinatorial Approach. We remove residues colored i by the obtained solution and solve the coloring problem on the remaining residues using K - 1 colors recursively.

Our approach works as follows. Consider an arbitrary color k ∈ \mathcal{S}. We compute a subset of the residues that are assigned color k such that the total error with respect to color k and the sum of all remaining colors is minimized, i.e. we solve the two color problem with requirements (= right hand sides)

(b^k,{\displaystyle \sum _{k^{\prime }\ne k}{b}^{k^{\prime }}}).

Residues assigned color k in an optimal solution to this problem will be colored k in the final solution too, the assignment of the remaining colors \mathcal{S}\{k} to the remaining residues is computed recursively.

Note that the order in which colors are selected to be the next fixed color k in the recursive computation can be arbitrary. Nevertheless, they might lead to solutions of different total error. As we have only three different colors in our experimental data, we evaluate all six orderings and return the best solution found.

In the next section we present a Lagrangian relaxation method to compute, based on our combinatorial approach for the 2-color case, a bound on the minimum total error, which is exploited in a branch-&-bound manner to determine all optimal colorings.

A Lagrangian Relaxation Approach

In this section we propose a Lagrangian relaxation approach for the problem, which is particularly suit-able for finding all optimal solutions. It is based on the improved-ILP formulation:

\min {\displaystyle \sum _{k\in \mathcal{S}}{\displaystyle \sum _{f\in ℱ}{e}_f^k}}

(8)

\begin{array}{ccc}\text{s}.\text{t}.& -A{y}^k+e^k\ge -b^k& \text{for}\text{all}k\in \mathcal{S}\end{array}

(9)

\begin{array}{ll}A{y}^k+e^k\ge b^k\hfill & \text{for}\text{all}k\in \mathcal{S}\hfill \end{array}

(10)

\begin{array}{l}{\displaystyle \sum _{k\in \mathcal{S}}{y}^k=P}\hfill \\ y\ge 0,\text{integer}\hfill \end{array}

(11)

where P is the vector that contains the length of parts in \mathcal{P}. The problem can be considered to contain independent structures for each color k ∈ \mathcal{S}, namely the set of positive integer vectors y^ksatisfying (9) and (10) under the objective (8), that are linked by constraints (11). Therefore, dualizing the linking constraints (11), with Lagrangian multipliers λ, splits the problem into an independent problem for each color k ∈ \mathcal{S}:

\begin{array}{lll}\min \hfill & {\displaystyle \sum _{k\in \mathcal{S}}{\displaystyle \sum _{f\in ℱ}{e}_f^k+{\lambda }^T\left({\displaystyle \sum _{k\in \mathcal{S}}{y}^k-P}\right)}}\hfill & (IP(\lambda ))\hfill \\ \text{s}.\text{t}.\hfill & -A{y}^k+e^k\ge -b^k\hfill & \text{for}\text{all}k\in \mathcal{S}\hfill \\ A{y}^k+e^k\ge b^k\hfill & \text{for}\text{all}k\in \mathcal{S}\hfill \\ 0\le y^k\le P\hfill & \text{for}\text{all}k\in \mathcal{S}\hfill \\ y\text{interger}\hfill \end{array}

Neglecting the constant term -λ^TP in the objective function and replacing error variable e by e + ē we have to determine, for every color k ∈ \mathcal{S}, an optimal integral solution to the following linear program:

\begin{array}{l}\min \sum _{f\in ℱ}\left(e_f^k+{\stackrel{ª}{e}}_f^k\right)+{\lambda }^T{y}^k\\ \text{s}.\text{t}.-A{y}^k+e^k\ge -b^k\end{array}

(12)

A{y}^k+{\stackrel{ª}{e}}^k\ge b^k

(13)

\begin{array}{l}{e}^k,{\stackrel{ª}{e}}^k\ge 0\\ 0\le y^k\le P\end{array}

(14)

Note that we added constraint (14) to enforce e_f^k or {\stackrel{ª}{e}}_f^k to be zero if {\stackrel{ª}{e}}_f^k, respectively e_f^k, corresponds to the absolute value of the error, i.e. if the constraint (13), respectively the constraint (12), for fragment f is tight. Note that we have to enforce e_f^k and {\stackrel{ª}{e}}_f^k to be nonnegative. In every optimum solution either e or ē (or both) will be zero for each fragment f. Similar as for linear program (5), its dual is given by (omitting the color superscript k):

\begin{array}{ll}-\min \hfill & b^T{f}^1-b^T{f}^2+P^T{f}^3\hfill \\ \text{s}\text{.t}\text{.}\hfill & -A^T{f}^1+A^T{f}^2-f^3+f^4=\lambda \hfill \\ 0\le f^{1,2}\le 1\hfill \\ 0\le f^{3,4}\hfill \end{array}

(15)

This linear program differs from LP (7) only in the right-hand sides of the equality constraints.

Conclusions

We applied our methods to process data from typical biochemical experiments. We report our results for four proteins: Calcium-binding protein (Cabin), Cytochrome P450 (CytoC), FK506 binding protein(FKBP), with two different digests (pepsin and XIII), and myoglobin. As a preprocessing step, the single fragments were analyzed with our integer linear programming based technique [14], except for FKBP V2 (MEM) which was analyzed with the MEM-method [15] and is based on the same data as FKBP V1 (ILP). We analyzed FKBP with only the xiii digestion (V3) and combined the datasets from the two digestions (V1, V2 and V4). The number discretized exchange rates per fragment obtained in this preprocessing step serves as input to the algorithm.

The instances have between 74 and 152 residues and between 18 and 49 fragments. The solutions with a minimal number of errors could be computed in less than 0.11 seconds for all instances. Computing all (non-equivalent) solutions with a minimal number of errors, from 96 up to almost 20 million in number, took less than 7 minutes, where the running time greatly depends on the number of solutions (see Table 2). Computing all solutions using the basic-ILP takes much longer as with the improved-ILP. For all instances, the heuristic computes a solution with the minimal error.

Table 2 Results & Runtime

Where available, we compared our assignments of exchange rates to the results obtained by NMR-analysis (FKBP and CytoC [16]). The error measure is based on a comparison per part. Within each part, the rates assigned by the algorithm are compared to the ones from NMR. Table 3 summarizes the results. The table also shows the importance of taking all solutions into account, as averaging typically yields better results than a single solution. The assignments coincide to 60 - 75% to the ones obtained by NMR, when choosing the optimal ordering with in the parts of equivalent residues. Figure 3 provides the results for FK506 and Cytochrome P450 at single residue resolution for manual inspection.

Table 3 Comparison with NMR

Figure 3

Single residue results for FKBP and cytochrome C. Results for FKBP (top) and cytochrome C (bottom) at single residue resolution. Rates from NMR are colored in red, our results from sHDX are colored in green. Vertical Black lines denote the boundary between consecutive parts. Legend for Symbols: '0' Rates agree on position, 'X' Rates disagree,'?' NMR data is missing, '-' sHDX data is missing. The horizontal lines indicate the range of the discretized values. slow ≤ 0.1 h^-1, and fast ≥ 8 h^-1. Secondary structure is indicated by the horizontal bars on the top: alpha helical in blue, beta sheets in green and loops in yellow.

A structural view on the results for FK506 and myoglobin is given in Figure 4. For myoglobin we do not have NMR data at hand. Nevertheless does the figure nicely agree with the expected out come, as buried parts of the protein show on the average lower exchange rates than exposed parts. The two figures have been produced by use of PyMOL [17].

Figure 4

Structural view of results. The exchange rates (fast, medium slow) are color-coded (red, yellow, green) in the structural view. The color blue means that we have no fragments covering this part and hence we do not predict anything. We present FKBP on the left and Myoglobin on the right.

In our solutions, the resolution is significantly increased compared to the input data, i.e., the length of fragments obtained from the sHDX experiments. The parts are typically small (see Table 2), between2 and 4 residues. 75% of the parts are smaller than 8 residues. For 46% of the amino-acids, we get single-residue resolution on the data.

The results for the real instances are very promising, as the small number of easily interpretable classes of equivalent solutions can be used in protein structure prediction tools and for manual inspection.

Methods

In this Section, we describe the computational methods, which we use to solve the different formulation, as well as the biochemical methods to obtain the experimental data.

Solution of the integer linear program

We implemented our approach using the C++-Library SCIL [18] to solve integer linear programs. SCIL uses the libraries LEDA [19] and SCIP [20].SCIP uses CPLEX [21] or SoPlex [22] as solver for linear programs. The underlying solution method is branch-&-bound, that is described in detail in [23].

In order to find all solutions within a given error bound e, the constraint

{\displaystyle \sum _{k\in S}{\displaystyle \sum _{(i,j)\in ℱ}{e}_{(i,j)}^k\le e}}

is added to the integer linear program and hence we are faced with the problem of computing all feasible solutions of an integer linear program. We do this with a branching-approach similar to the classical branch-&-bound method for finding an optimal solution: First, the linear relaxation is solved. If the linear relaxation is infeasible, the search on this branch terminates. If the solution is integral, it is stored (provided the solution was not found yet). If there is a binary variable which was not fixed so far (i.e. not set to 0 or 1), one such variable x_l^k is picked and the two subproblems, in which the variable is fixed o 0 or 1 recursively, are solved. Notice that it is possible that we branch on a variable which already has an integral value. In this case, the solution of the linear relaxation of the subproblem will be the same as in problem itself. Nevertheless, we will terminate, as there are only a finite number of variables to branch on.

Solving the Combinatorial Problem

We may use any algorithm that solves the Minimum Cost Circulation problem, e.g. Cycle Canceling or Successive Shortest Path (see [12] for further reference). Both approaches have their advantages. The former always maintains a feasible circulation, i.e. we start with the zero flow and augment flow along negative cycles in the residual network until no negative cycle remains. Since the residual network with respect to an optimal circulation does not contain a directed negative circuit, we can find node potentials, i.e. a corresponding dual solution, using the Bellman-Ford algorithm in O\left(\left|\mathcal{P}\right|\cdot \left|ℱ\right|\right) time. The difference between the potential of two neighboring nodes then yields the value of the corresponding y-variable. The errors are determined straight forward. If there is a solution without error this approach yields a solution within the running time of Bellman-Ford. On the other hand, the Successive Shortest Path algorithm maintains similar node potentials such that the arc-weights remain non-negative. Since the total excess is bounded by |\mathcal{P}| in our case, the running time of that algorithm is O\left(\left|\mathcal{P}\right|·\left|ℱ\right|+\left|\mathcal{P}{|}^2\log \right|\mathcal{P}|\right).

Solving the Lagrangian Dual

Instead of a minimum cost circulation problem(right-hand side is 0), we have to solve the more general minimum cost flow problem [12] where the supplies and demands \stackrel{ª}{\lambda } of the nodes are determined by the difference of Lagrangian multipliers, i.e. \stackrel{ª}{\lambda } is of dimension |\mathcal{P}| + 1 and {\stackrel{ª}{\lambda }}_i={\lambda }_i-{\lambda }_{i-1} for 2\le i\le \left|\mathcal{P}\right|,{\stackrel{ª}{\lambda }}_1={\lambda }_1 and {\stackrel{ª}{\lambda }}_{\left|\mathcal{P}\right|+1}=-{\lambda }_{\left|\mathcal{P}\right|}. A feasible flow of minimum cost can be computed efficiently by, e.g., the cycle-canceling algorithm and the successive shortest path algorithm, as well as variants of them, like the capacity scaling algorithm[12]. In our implementation (C++) we used the LEDA library [19] to solve the Lagrangian subproblem by analgorithm based on capacity scaling and successive shortest path computation [12].

We improve the resulting bounds by the subgradient optimization method described in the following and incorporate the overall approach into a branch-&-bound algorithm as the lower bounding scheme.

Let v(IP (λ)) denote the optimal value of IP (λ). Then for any vector λ of Lagrangian multipliers, the (non-differentiable) Lagrangian function

z(\lambda )=v(IP(\lambda ))

provides a lower bound on the minimum total error. To benefit from the sharpest possible bound in the branch-&-bound framework we are interested in solving the Lagrangian dual problem

z^{*}=\underset{\lambda }{\max }z(\lambda ).

We apply the subgradient method to obtain near-optimal Lagrangian multipliers. Following the approach by Held and Karp [24] we iteratively determine values λ^ℓ+1for ℓ = 0,1, ..., of the Lagrangian multipliers by moving in the direction of a subgradient with "step length" μ _ℓ :

{\lambda }^{\ell +1}={\lambda }^{\ell }+{\mu }_{\ell }\left({\displaystyle \sum _{k\in S}{y}^k}(\ell )-P\right),

where {(y^k(\ell ))}_{k\in S} is any optimal solution to IP(λ^ℓ ). The step length is computed according to formula

{\mu }_{\ell }=\frac{{\theta }_{\ell }\left(UB-z\left({\lambda }^k\right)\right)}{{\Vert {\Sigma }_{k\in S}{y}^k\left(\ell \right)-P\Vert }^2},

where UB is a previously computed upper bound on z* and θ is a step size parameter assuming values in {x ∈ ℝ | 0 < x ≤ 2}. In the experiments it turns out, that initializing the vector of Lagrangian multipliers λ⁰ to the length P of the corresponding intervals in \mathcal{P} increases the convergence rate dramatically. We also experienced a fast convergence to near-optimal Lagrangian multipliers when following the classical Held-Karp method to choose the step size scalar θ: We start with θ₀ = 2 and half θ_ℓ whenever the best Lagrangian bound v(IP(λ)) found so far has not increased in a certain number of iterations. As soon as the step size scalar falls below a specified threshold or the number of iterations exceeds a certain limit (which is adaptive with respect to the depth of the branch-&-bound node), we branch on a variable y_p^k,k\in S,p\in \mathcal{P}, such that{\stackrel{ª}{y}}_p^k-\lfloor {\stackrel{ª}{y}}_p^k\rfloor is close to 0.5, where {\stackrel{ª}{y}}_p^k is the average value of variable y_p^k in the last h = 10 Lagrangian solutions. Since we aim to find all optimal colorings, we also branch on variables that are integral. Incorporating the Lagrangian approach as a lower bounding scheme into a branch-&-bound frame work gives an alternative algorithm that does not depend on commercial software packages.

Experimental Setting

The entire sHDX experiment was automated with a LEAP robot (HTS PAL, Leap Technologies, Carrboro, NC).

Automation of the experiment reduces human error and reduces deuterium for hydrogen back-exchange. All time points where interlaced and performed in triplicate to ensure experimental reproducibility. After digestion, the protein digest was injected from a 10 μ L loop to either a 1 mm × 50 mm C5 column (Phenomenex) or a Pro-Zap Pro-sphere HP C18 HR 1.5u 10 mm × 2.1 mm (All-tech). A rapid gradient 2% B to 95% B in 1.5 min(A: acetonitrile/H₂O/formic acid 5/94.5/0.5, B: ace-tonitrile/H₂O/formic acid 95/4.5/0.5) was used to elute peptides. The eluent was post-column split and infused by microelectrospray ionization into a custom built 14.5 T LTQ FT-ICR mass spectrometer. The extraction of the peptic fragments and their deuterium uptakes from these data was done by an in-house analysis package [25]. Then we compute the cumulative exchange rates from the deuterium uptakes with either the MEM-method [15] or a new approach based on integer linear programming [14].

A current limitation for implementation of this software is back exchange of deuterium-to-hydrogen during the separation of the samples. It has been reported that different peptides have a different percentage of back exchange due to the sequence of amino acids [26, 27]. Furthermore, the peptide sequence overlap will limit the ability to map single amino acid rate kinetics. Thus, reduction of backexchange has been investigated [28, 29], along with multiple acid proteases to increase sequence coverage [30]. The sHDX experiment is continually being improved, but in its current state the sHDX experiment does not take away from the integrity of the algorithm to discern single amino acid exchange kinetics.