Proof of convexity

Disciplined Convex Programming (DCP)

DCP is a methodology that imposes a set of conventions that one must follow when constructing convex programs. The conventions are simple and teachable, taken from basic principles of convex analysis, and inspired by the practices of those who regularly study and apply convex optimization today. The conventions do not limit generality; but they do allow much of the manipulation and transformation required to analyze and solve convex programs to be automated.

To read more about DCP, we recommend this paper by M. Grant, S. Boyd, and Y. Ye.

For our implementation of Convex Model Balancing, we used CVXPY - a Python-embedded modeling language for convex optimization problems.

The CVXPY function that are required for Model Balancing are:

  • exp\((x)\) - elementwise \(e^x\)

  • log1p\((x)\) - elementwise \(\ln(1 + x)\)

  • logistic\((x)\) - elementwise \(\ln(1 + e^x)\)

  • log_sum_exp\((x)\) - alias for \(\ln(\Sigma_i e^{x_i})\)

  • quad_form\((x,P)\) - alias for \(x^\top P x\)

  • pos\((x,P)\) - alias for elementwise \(max(x, 0)\)

The enzyme concentration term is a convex function

In this section, we will show that the most complex term in Model Balancing (i.e. the enzyme abundance score) is a convex function using DCP. This proof only works when dropping the negative part of the z-score (i.e. taking the square only for the cases where the enzyme concentrations are higher than the geometric mean). According to the notation of Model Balancing, this is the case when one defines α to be 0.

The following formulae show how we calculate the vector of log enzyme concentrations \(\mathbf{z}\) for a single condition. If the model contains more than one condition, the total score will simply be the sum over the conditions.

In log-scale, \(\mathbf{z}\) is a combination of four convex expressions:

  • \(\mathbf{z}^{min}\) - representing the minimal enzyme log-concentrations (when working in full capacity)

  • \(\ln(\mathbf{η}^{thr})\) - a thermodynamic efficiency term for being close to equilibrium (i.e. low driving forces)

  • \(\ln(\mathbf{η}^{kin})\) - a kinetic efficiency term for saturation effects of substrates and products

  • \(\ln(\mathbf{η}^{reg})\) - a regulation efficiency term for allosteric regulations (inhibitors and activators)

Each one of these expressions is convex (according to DCP) and therefore their sum is convex as well. We show this by writing down the functions that calculate each term, and ensure that CVXPY is able to verify that they are convex.

\[\begin{split}\mathbf{z}^{min} &= \ln(\mathbf{v}) - \ln(\mathbf{k^{cat+}}) \\ \ln(\mathbf{η}^{thr}) &= -\text{log1p} \left( - \exp \left( \mathbf{S}^\top \mathbf{x} - \ln(\mathbf{k^{eq}}) \right) \right) \\ \forall l\;\;\ln(η^{kin}_l) &= -\text{log_sum_exp} \left[\mathbf{B}_l (\mathbf{x} - \ln(\mathbf{k^M}_l)) \right]\\ \forall l\;\;\ln(η^{reg}_l) &= -\mathbf{h^I}_l^\top~\text{logistic}(\mathbf{x} - \ln(\mathbf{k}^I_l)) - \mathbf{h^A}_l^\top~\text{logistic}(\ln(\mathbf{k}^A_l) - \mathbf{x})\end{split}\]

where \(\mathbf{v}\) is the flux vector (in M/sec), \(\mathbf{x}\) is the vector of metabolite log-concentrations (in log(M)), \(\mathbf{S}\) is the stoichiometric matrix, \(\mathbf{k^{cat+}}\) is the vector of turnover numbers in the direction of positive flux (in sec), \(\mathbf{k^M}_l\) is the vector of the reactant affinity coefficients for reaction \(l\) (in M), \(\mathbf{k^I}_l\) is the same for inhibitors and \(\mathbf{k^A}_l\) for activators. The vectors \(\mathbf{h^I}_l\) and \(\mathbf{h^A}_l\) hold the Hill coefficients for inhibitors and activators (respectively) or zero for all other metabolites. Note that the formulae for \(\ln(\mathbf{η}^{kin})\) and \(\ln(\mathbf{η}^{reg})\) are computed by combining the expressions for all the single reactions into one vector.

The only symbols left to define are the \(\mathbf{B}_l\), which are a set of constant auxiliary matrices, each one depending only on the stoichiometry of the reaction \(l\).

Defining the B matrix

For the sake of an example, let us assume a substrate list of length 4. We define matrices with rows corresponding to all single items, pairs, triples, and quadruples (and columns corresponding to the elements of \(\mathbf{x}\)):

\[\begin{split}\mathbf{K}_{1}^{(4)} &=& \left(\begin{array}{llll}1 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1\\\end{array}\right)\\ \mathbf{K}_{2}^{(4)} &=& \left(\begin{array}{llll}1 & 1 & 0 & 0\\1 & 0 & 1 & 0 \\1 & 0 & 0 & 1\\ 0 &1&1&0 \\ 0& 1&0&1\\0&0&1&1\end{array}\right)\\ \mathbf{K}_{3}^{(4)} &=& \left(\begin{array}{llll}1 & 1& 1& 0 \\1 & 1 & 0 &1 \\1&0&1 &1 \\0&1&1&1\end{array}\right)\\ \mathbf{K}_{4}^{(4)} &=& \left(\begin{array}{llll}1 & 1&1&1 \\\end{array}\right)\end{split}\]

and combine them in a matrix:

\[\begin{split}\mathbf{K}^{(4)} = \left( \begin{array}{l} \mathbf{K}_{1}^{(4)}\\ \mathbf{K}_{2}^{(4)}\\ \mathbf{K}_{3}^{(4)}\\ \mathbf{K}_{4}^{(4)}\\ \end{array} \right)\,.\end{split}\]

We can generalise this from the case \(n_s = 4\) to substrate lists of any length, and obtain the respective matrix \(\mathbf{K}^{n_s}\). Altogether, for a reaction with substrate list length \(n_s\) and product list length \(n_p\), we define another auxiliary matrix \(\mathbf{A}\):

\[\begin{split}\mathbf{A} = \left( \begin{array}{l} \mathbf{0}\\ \mathbf{K}^{n_s} \,\mathbf{S}^s\\ \mathbf{K}^{n_p} \,\mathbf{S}^p \end{array} \right)\end{split}\]

where \(\mathbf{S}^s\) is a matrix (with \(n_s\) rows and with columns corresponding to all metabolites in the model) where each row is an indicator vector for the index of one of the substrates in reaction \(l\). Likewise, \(\mathbf{S}^p\) is a matrix of indicator vectors for the products of this reaction.

Finally, we can define the \(\mathbf{B}\) matrix:

\[\mathbf{B} = \mathbf{A} - (\mathbf{1}\, \mathbf{1}^\top) \,\mathbf{S}^s.\]

Here, for simplicity, we did not use the subscript \(l\), but the size and values of \(\mathbf{A},\mathbf{B},\mathbf{S}^s,\mathbf{S}^p\) is dependent on the reaction. This procedure must to be repeated for every reaction in the model in order to generate the set of auxiliary matrices.

Calculating the z-score for enzyme concentration

So, we now have an expression for the enzyme log-concentrations \(\mathbf{z}\):

\[\mathbf{z} = \mathbf{z}^{min} - \ln(\mathbf{η}^{thr}) - \ln(\mathbf{η}^{kin}) - \ln(\mathbf{η}^{reg})\]

which is convex according to CVXPY. To calculate the log likelihood, we need to take the square difference from the mean (normalized by the covariance), but the square of a convex function will only be convex if the inner function is strictly positive. Since there is no guarantee that the calculated values are always higher than the mean, we need to apply the pos function which keeps only positive values and discard all negative values (replacing them with 0):

\[L = \text{quad_form}(\text{pos}(\mathbf{z} - \mathbf{\mu}_z), \mathbf{\Sigma}_z^{-1})\,.\]

According to CVXPY, this function is convex.