Quantitative flux coupling analysis

Abstract

Flux coupling analysis (FCA) aims to describe the functional dependencies among reactions in a metabolic network. Currently studied coupling relations are qualitative in the sense that they identify pairs of reactions for which the activity of one reaction necessitates the activity of the other one, but without giving any numerical bounds relating the possible activity rates. The potential applications of FCA are heavily investigated, however apart from some trivial cases there is no clue of what bottleneck in the metabolic network causes each dependency. In this article, we introduce a quantitative approach to the same flux coupling problem named quantitative flux coupling analysis (QFCA). It generalizes the current concepts as we show that all the qualitative information provided by FCA is readily available in the quantitative flux coupling equations of QFCA, without the need for any additional analysis. Moreover, we design a simple algorithm to efficiently identify these flux coupling equations which scales up to the genome-scale metabolic networks with thousands of reactions and metabolites in an effective way. Furthermore, this framework enables us to quantify the “strength” of the flux coupling relations. We also provide different biologically meaningful interpretations, including one which gives an intuitive certificate of precisely which metabolites in the network enforce each flux coupling relation. Eventually, we conclude by suggesting the probable application of QFCA to the metabolic gap-filling problem, which we only begin to address here and is left for future research to further investigate.

Appendix

1.1 Derivation of the dual problem for (10)

For the standard definitions from the theory of Lagrange duality used in this appendix, we refer the reader to the fifth chapter of (Boyd and Vandenberghe 2004).

By definition, the Lagrangian for the LP (10) is equal to

$$\begin{aligned} \mathcal {L}(u, v, \lambda _1, \lambda _2, \lambda _3, \nu ) = -{\mathbf {1}}^T u + \nu ^TS v+\lambda _1^T(u-v_I)+\lambda _2^T(u-{\mathbf {1}})+\lambda _3^T(-u), \end{aligned}$$

where $\lambda _1, \lambda _2, \lambda _3\in \mathbb {R}^k$, and $\nu \in \mathbb {R}^m$ are the Lagrange dual variables. The Lagrange dual function for this Lagrangian is defined as

$$\begin{aligned} g(\lambda _1, \lambda _2, \lambda _3, \nu ) = \inf _{u \in \mathbb {R}^k,v \in \mathbb {R}^n} \mathcal {L}(u, v, \lambda _1, \lambda _2, \lambda _3, \nu ). \end{aligned}$$

Let $p^\star$ denote the optimal objective value and assume $\lambda _i \ge 0$ for $i = 1,2,3$. If $u^\star$ and $v^\star$ are optimal, then

$$\begin{aligned} \begin{aligned} g(\lambda _1, \lambda _2, \lambda _3, \nu )&= \inf _{u \in \mathbb {R}^k,v \in \mathbb {R}^n} \mathcal {L}(u, v, \lambda _1, \lambda _2, \lambda _3, \nu )\\&\le \mathcal {L}\big (u^\star , v^\star , \lambda _1, \lambda _2, \lambda _3, \nu \big )\\&= -{\mathbf {1}}^T u^\star + \nu ^T S v^\star +\lambda _1^T\big (u^\star -v_I^\star \big )+\lambda _2^T(u^\star -{\mathbf {1}}) +\lambda _3^T(-u^\star )\\&= -{\mathbf {1}}^T u^\star +\lambda _1^T\big (u^\star -v_I^\star \big )+\lambda _2^T(u^\star -{\mathbf {1}}) +\lambda _3^T(-u^\star )\\&\le -{\mathbf {1}}^Tu^\star \\&= -p^\star . \end{aligned} \end{aligned}$$

Therefore, for any nonnegative $\lambda _1, \lambda _2, \lambda _3$, and any $\nu$, the negative Lagrange dual function yields an upper bound on $p^\star$. In order to find the tightest bound, we should solve the following Lagrange dual problem

$$\begin{aligned} \begin{array}{ll} \text{ maximize } &{} {\mathbf {1}}^T\lambda _2 \\ \text{ subject } \text{ to } &{} S^T \nu = \left[ \begin{array}{c} \lambda _1\\ 0 \end{array} \right] \\ &{} \lambda _1+\lambda _2 \ge {\mathbf {1}}\\ &{} \lambda _1 \ge 0 \\ &{} \lambda _2 \ge 0, \end{array} \end{aligned}$$

over the Lagrange dual variables. The proof follows from rewriting the Lagrangian as

$$\begin{aligned} \mathcal {L}(u, v, \lambda _1, \lambda _2, \lambda _3, \nu ) = \left( S^T\nu -\left[ \begin{array}{c} \lambda _1\\ 0 \end{array} \right] \right) ^T v-\lambda _2^T {\mathbf {1}}+(-{\mathbf {1}}+\lambda _1+\lambda _2-\lambda _3)^T u. \end{aligned}$$

Again, the dual LP is always both feasible (e.g., $\lambda _1=\lambda _2=0$, $\nu =0$ is feasible) and bounded (${\mathbf {1}}^T \lambda _2 \le k$).

For this primal-dual pair of LPs strong duality holds which means that the gap between the dual optimal objective and the primal optimal objective $p^\star$ is zero. In other words, the bound given by the optimal dual variables is sharp.

It is easily seen that

$$\begin{aligned} \lambda _2^\star =\max ({\mathbf {1}}-\lambda _1^\star ,0), \end{aligned}$$

has either zero or one entries just like the optimal $u^\star$. However, from zero duality gap

$$\begin{aligned} {\mathbf {1}}^T u^\star ={\mathbf {1}}^T\lambda _2^\star , \end{aligned}$$

hence $\lambda _2^\star$ and $u^\star$ have the same number of ones. Also by complementary slackness, $\lambda _2^\star$ is zero wherever $u_i^\star \ne 1$. Altogether, $\lambda _2^\star$ and $u^\star$ are 0-1 vectors with the same sparsity pattern, thus they are equal.

Ultimately, we can also rewrite the dual problem as (13) in analogy to the primal problem (9), by substituting

$$\begin{aligned} \lambda =\left[ \begin{array}{c} \lambda _1\\ 0 \end{array} \right] . \end{aligned}$$

1.2 Proofs of Sect. 2.4

Proof of Theorem 2.1

From the Lemma 2.2, we can assume the seemingly stronger but equivalent right hand side of (4), which by definition means that $\lambda \in \mathrm {ker}(S)^\perp$. From rank-nullity theorem we know that $\mathrm {ker}(S)^\perp = \mathrm {range}(S^T)$. Thus, $\lambda \in \mathrm {range}(S^T)$ which in turn implies the desired result. The converse also holds because in the reverse direction we have $\lambda = S^T \nu$ and hence, for any $v\in \mathcal {C},$

$$\begin{aligned} \lambda ^Tv = \nu ^T Sv = \nu ^T 0 = 0. \end{aligned}$$

$\square$

Proof of Lemma 2.2

($\Leftarrow$) is immediate from $\mathcal {C}\subseteq \mathrm {ker}(S)$. For the other direction, suppose that the left-hand side is true.

The proof goes by contradiction. Assume to the contrary that there exists $u\in \mathrm {ker}(S)$ such that $\lambda ^Tu\ne 0$. Let

$$\begin{aligned} v = u+\sum _{i:R_i\in \mathcal {I}}|u_i |v^i, \end{aligned}$$

where for any $R_i\in \mathcal {I}$, $v^i$ is an arbitrary feasible flux distribution with its ith flux coefficient equal to one, namely $v^i_i=1$. One can easily show that the feasibility constraints (1) and (2) hold for v, and hence $v\in \mathcal {C}$.

On the other hand,

$$\begin{aligned} \lambda ^Tv - \sum _{i:R_i\in \mathcal {I}}|u_i |\lambda ^Tv^i = \lambda ^{T} u \ne 0, \end{aligned}$$

by the way we constructed v. Therefore, at least one of the terms on the left hand side is nonzero. The proof is complete since this is in contradiction with the assumption that $\lambda ^Tv=0$ for all $v\in \mathcal {C}$. $\square$

1.3 Reversibility type correction

Earlier in Sect. 2.5, we have derived a naive method for reversibility type correction by solving $2n_r$ LPs. Recall that for $R_j \notin \mathcal {I}$, our task is to figure out if both its forward and reverse directions are unblocked, and as always we assume that all the blocked reactions are already removed.

Following the same fashion, this time we try to search for positive certificates proving that either the forward or reverse direction of $R_j$ becomes blocked when $R_j$ is added to $\mathcal {I}$ (for the reverse direction we should also replace $S^j$ by $-S^j$). Applying (13) to the resulting modified metabolic network, either $\lambda ^\star = 0$ or $\lambda ^\star _j\ne 0$, otherwise the same $\lambda ^\star$ can be considered as a positive certificate for the original metabolic network where $R_j \notin \mathcal {I}$ proving that some reactions other than $R_j$ are blocked, which is in contradiction to the assumption that we have already removed all the blocked reactions.

If $\lambda ^\star _j < 0$, then $\frac{\lambda ^\star }{-\lambda ^\star _j}$ is in fact a DCE which clearly shows $v_j > 0$ for all $v\in \mathcal {C}$. If $\lambda ^\star _j > 0$, then $\frac{\lambda ^\star }{\lambda ^\star _j}$ shows $v_j < 0$ for all $v\in \mathcal {C}$ and becomes a DCE if we replace $S^j$ by $-S^j$ permanently. If $\lambda ^\star _j = 0$ and hence $\lambda ^\star = 0$ for both directions, then $R_j$ is truly reversible. As a consequence, we have just shown that in the modified metabolic network either $R_j$ is not blocked in the selected direction or some other irreversible reactions should also become blocked which means that $R_j$ is effectively irreversible because some irreversible reactions are directionally coupled to it.

In this way, DCEs can also be used for reversibility type correction with two major advantages over the naive method. The first one is that in order to derive DCEs with either $\lambda ^\star _j < 0$ or $\lambda ^\star _j > 0$, it is enough to solve (17) as there is no constraint on the sign of $\lambda _j$. This is a twofold decrease in the number of required LPs and brings the total number of them down to $n_r$. The second advantage is that we also get the $\mathcal {D}_j$ for free by the same LP.