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

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.

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}$.

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.

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

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.

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.

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

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

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].

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.

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.