On the identifiability of metabolic network models

A major problem for the identification of metabolic network models is parameter identifiability, that is, the possibility to unambiguously infer the parameter values from the data. Identifiability problems may be due to the structure of the model, in particular implicit dependencies between the parameters, or to limitations in the quantity and quality of the available data. We address the detection and resolution of identifiability problems for a class of pseudo-linear models of metabolism, so-called linlog models. Linlog models have the advantage that parameter estimation reduces to linear or orthogonal regression, which facilitates the analysis of identifiability. We develop precise definitions of structural and practical identifiability, and clarify the fundamental relations between these concepts. In addition, we use singular value decomposition to detect identifiability problems and reduce the model to an identifiable approximation by a principal component analysis approach. The criterion is adapted to real data, which are frequently scarce, incomplete, and noisy. The test of the criterion on a model with simulated data shows that it is capable of correctly identifying the principal components of the data vector. The application to a state-of-the-art dataset on central carbon metabolism in Escherichia coli yields the surprising result that only \(4\) out of \(31\) reactions, and \(37\) out of \(100\) parameters, are identifiable. This underlines the practical importance of identifiability analysis and model reduction in the modeling of large-scale metabolic networks. Although our approach has been developed in the context of linlog models, it carries over to other pseudo-linear models, such as generalized mass-action (power-law) models. Moreover, it provides useful hints for the identifiability analysis of more general classes of nonlinear models of metabolism.

1. Strictly speaking, a better version of Definition 2 would require that condition (18) holds for almost all \(p^*\) in \(P\). This would automatically rule out trivial definitions of \(\hat{B}_{C(i),i}\) such as \(\hat{B}_{C(i),i}\triangleq B_{C(i),i}^*\) (which makes the reaction identifiable for any \(\alpha \) and \(\fancyscript{B}_i\) but cannot be built without the knowledge of \(B_{C(i),i}^*\) itself). Unfortunately, this is not a good choice in general, in that estimation uncertainty may severely depend on \(p^*\) itself, as we shall see later on in Example 4. For simplicity, here we stick to Definition 2 with the understanding that any such triviality is avoided.

2. This holds as an equality in the sense of expectation, and can also be motivated by asymptotic arguments as \(q\rightarrow +\infty \).

Authors and Affiliations

1. INRIA Grenoble-Rhône-Alpes, Montbonnot, France

   Sara Berthoumieux, Hidde de Jong & Eugenio Cinquemani

2. Laboratoire de Biométrie et Biologie Evolutive, CNRS UMR 5558, Université Lyon 1, INRA, Villeurbanne, France

   Matteo Brilli & Daniel Kahn

Corresponding author

Correspondence to Eugenio Cinquemani.

This work was supported by the Agence Nationale de la Recherche under project MetaGenoReg (ANR-06-BYOS-0003).

Appendix A: Notation and terminology

\(\mathbb R , \mathbb R _{> 0}, \mathbb Z \) and \(\mathbb N \) denote the sets of real, positive real, integer and positive natural numbers, respectively. For an index \(n\in \mathbb N , \mathbb R ^n\) and \(\mathbb R _{> 0}^n\) denote the \(n\) dimensional versions of \(\mathbb R \) and \(\mathbb R _{> 0}\). \(I\) denotes an identity matrix of dimension fixed by the context.

Let \(M\) be any matrix. For two indices \(i\) and \(j\) and a vector of indices \(C\) compatible with the dimensions of \(M, M_i\) denotes the \(i\)th column of \(M, M_{C}\) denotes the submatrix of \(M\) formed by the columns of \(M\) with indices \(C, M_{j,i}\) denotes the element of \(M\) in row \(j\) and column \(i\), and \(M_{j,C}\) denotes the row vector formed by the elements of \(M\) in row \(j\) and columns indexed by \(C\). \([M]_{C,C}\) denotes the minor formed by the rows and columns indexed by \(C\). When convenient, notation \(M_{i:i^{\prime }}\) with \(i^{\prime }\ge i\) is used instead of \(M_C\) with \(C=\begin{bmatrix}i,i+1,\ldots ,i^{\prime }\end{bmatrix}\). For vectors, a subscript \(i\) refers to the \(i\)th element of the vector.

For a square matrix \(\Sigma , \Sigma > 0\) (resp. \(\Sigma \ge 0\)) means that \(\Sigma \) is positive definite (resp. semidefinite). For a vector \(\mu \) of suitable dimension, \(\varepsilon \sim \fancyscript{N}(\mu ,\Sigma )\) means that \(\varepsilon \) is a Gaussian random vector with mean \(\mu \) and covariance matrix \(\Sigma \). For a scalar quantity observed in Gaussian noise with standard deviation \(\nu \), a noise level of \(N~\%\), with \(N\in \mathbb R _{\ge 0}\), will stand for a value of \(\nu \) such that \(\nu \) is equal to \(N~\%\) of the value of that quantity (such that \(\simeq \)99 % of the noise outcomes fall within \(\pm 3\cdot N~\%\) of the observed quantity).

For two vectors \(v\) and \(e\) of equal size, both \(v/e\) and \(\frac{v}{e}\) indicate the vector obtained by element-wise division. Given a vector sequence \(v^1,\ldots , v^q\), \(\bar{v}\) is the mean \((1/q)\sum _{k=1}^q v^k\). For vectors and sets, \(|\cdot |\) denotes vector dimension and set cardinality, respectively. For a function \(f:A\rightarrow B\), \(f|_D\) indicates its restriction on \(D\subseteq A\).

Appendix B: Mathematical proofs

Proposition 1

A reaction \(i\) of \(\fancyscript{M}_p\) is structurally identifiable at \(p^*\) if and only if there exists \(D=\{ (e^1,u^1), \ldots , (e^q,u^q)\} \subseteq E\times U\) such that the solution of the equation \(W_i^*=Y^* B_i^*\), with

is unique in the parameters \(B_i^*=\big ([{B^x}^* \; {B^u}^*]^T\big )_{i}\).

Proof

(If) Assume that, for a given \(D \subseteq E\times U\), the solution of \(W_i^*=Y^* B_i^*\) is unique. We need to prove that \(\big ((J_p)_i,x_p\big )|_D = \big ((J_{p^*})_i,x_{p^*}\big )|_D\) implies \(p_i=p_i^*\). For simplicity, here we drop index \(i\) from subscripts.

Given any two parameters \(p^*=\begin{bmatrix}a^*&{B^*}^T\end{bmatrix}^T\) and \(p=\begin{bmatrix}a&B^T\end{bmatrix}^T\), for which \(\fancyscript{M}_p:(e,u)\mapsto (J_p,x_p)\) and \(\fancyscript{M}_{p^*}:(e,u)\mapsto (J_{p^*},x_{p^*})\), it holds by construction that \(W=YB\) and \(W^*=YB^*\). If \((J_p,x_p)|_D = (J_{p^*},x_{p^*})|_D\), then it also holds that \(Y=Y^*\) and \(W=W^*\), therefore we can write \(W^*=Y^*B\). Because the solution in \(B\) of the latter is unique and one solution is \(B^*\), it follows that \(B=B^*\). To conclude that \(p=p^*\), we are left with showing that \(a=a^*\). This follows from

(Only if) Here the hypothesis is that, for a given \(D \subseteq E\times U\), \(\big ((J_p)_i,x_p\big )|_D = \big ((J_{p^*})_i,x_{p^*}\big )|_D\) implies \(p_i=p^*_i\), and we need to show that the solution in \(B_i\) of \(W_i^* = Y^* B_i\) is unique. For simplicity, we will again drop \(i\) from the subscripts.

For the sake of contradiction, assume that \(W^* = Y^* B\) admits distinct solutions. Since \(B^*\) is a solution, all solutions are of the form \(B=B^*+z\), with \(z\) in the nontrivial kernel of \(Y^*\). For any such \(z\) we can write \(Y^*B^*=Y^*(B^*+z)\), i.e., \(\forall (e,u)\in D\),

Let \(p^*=\begin{bmatrix}a^*&{B^*}^T\end{bmatrix}^T\). For any \((e,u)\in D, J_*=J_{p^*}(e,u)\) and \(x_*=x_{p^*}(e,u)\) are given by the solution of

(which is unique by virtue of Assumption 1). Using (28), term \(\begin{bmatrix} \ln x_*^T&\ln u_*^T\end{bmatrix}B^*\) can be rewritten as \(\begin{bmatrix} \ln x_*^T&\ln u_*^T \end{bmatrix}(B^*+z) - \begin{bmatrix} \overline{\ln x_*^T}&\overline{\ln u_*^T}\end{bmatrix}z\). Replacing this into (29) yields

From this we see that \(p=\begin{bmatrix}a&{B}^T\end{bmatrix}^T\), with \(a\) defined as above, is different from \(p^*\) but is such that \(\big (J_p(e,u),x_p(e,u)\big )=\big (J_{p^*}(e,u),x_{p^*}(e,u)\big )\) for all \((e,u)\in D\), which contradicts the hypothesis. \(\square \)

Corollary 1

A reaction \(i\) of \(\fancyscript{M}_p\) is structurally identifiable at \(p^*\) if and only if there exists \(D=\{ (e^1,u^1) \,, \, \cdots \,, (e^q,u^q)\}\subseteq E\times U\) such that \(Y^*_{C(i)}\) is full column-rank.

Proof

From Proposition 1, we know that identifiability is equivalent to the uniqueness of the solution in \(B_i\) of \(W_i^* = Y^* B_i\), i.e. of the solution in \(B_{C(i),i}\) of \(W_i^* = Y^*_{C(i)} B_{C(i),i}\) (the elements of \(B_i\) not included in \(B_{C(i),i}\) are set to zero by definition). Uniqueness holds if and only if \(\ker (Y^*_{C(i)})=\{0\}\), i.e. \(Y^*_{C(i)}\) is full column-rank or equivalently \(\mathrm{{rank}}(Y^*_{C(i)})=n_i\). \(\square \)

Proposition 2

If a reaction \(i\) of \(\fancyscript{M}_p\) is structurally identifiable at \(p^*\) in the sense of Definition 1 then, for every \(\alpha \in (0,1)\), it is practically identifiable in the sense of Definition 2 with confidence level at least \(1-\alpha \) for any uncertainty set \(\fancyscript{B}_i\supseteq \fancyscript{E}_{\widehat{\Sigma }}(\alpha )\), where \(\fancyscript{E}_{\widehat{\Sigma }}(\alpha )\) denotes the \((1-\alpha )\)-confidence ellipsoid of a zero-mean Gaussian distribution with variance \(\widehat{\Sigma }=(Y_{C(i)}^T\Sigma _{\varepsilon _{i}}^{-1} Y_{C(i)})^{-1}\).

Proof

The definition of \(\fancyscript{M}_p\) ensures that \(W_i^* = Y_{C(i)}^* B_{C(i),i}^*\). Using the fact that \(Y^*=Y\) and \(W^*=W\), given a noisy dataset \(\tilde{W}_i = W_i + \varepsilon _i\) and the errorless dataset \(Y\), the regression problem becomes

From Sect. 3.1, identifiability of \(\fancyscript{M}_p\) at \(p^*\) for the given input set \(D\) is equivalent to \(Y_{C(i)}^*=Y_{C(i)}\) being full column-rank. Thus, if \(D\) ensures structural identifiability of \(\fancyscript{M}_p\) at \(p^*\), \(Y_{C(i)}\) is full column-rank and the weighted pseudoinverse of \(Y_{C(i)}\), defined as \(Y^\dag \triangleq \left(Y_{C(i)}^T\Sigma _{\varepsilon _{i}}^{-1}Y_{C(i)}\right)^{-1} Y^T\Sigma _{\varepsilon _{i}}^{-1}\), is well-defined. This enables us to define the minimum variance estimator of \(B^*_{C(i),i}\), \(\hat{B}_{C(i),i}= Y^\dag W_i\). From the linearity of the estimator in the Gaussian noise \(\varepsilon _{i}\), after simple calculations of first and second-order moments, one gets \(\hat{B}_{C(i),i}\sim \fancyscript{N}\left(B^*_{C(i),i},\widehat{\Sigma }\right)\) (also compare Ljung 1999, Appendix II). Thus, from the definition of \(\fancyscript{B}_i\), \(\mathbb P _{p^*} [\hat{B}_{C(i),i}-B^*_{C(i),i}\in \fancyscript{B}_i]\ge \mathbb P _{p^*} [\hat{B}_{C(i),i}-B^*_{C(i),i}\in \fancyscript{E}_{\widehat{\Sigma }}(\alpha )] = 1-\alpha \). \(\square \)

Proposition 3

Let index \(j\) be an element of \(C(i)\), i.e., for some \(\ell \!=\!1,\ldots , n_i\), \(j\!=\!C_\ell (i)\). Suppose that the entries of \(V_{\ell ,r+1:n_i}\) are all zero (that is, the \(\ell \)th entry of all vectors in \(K_Y\) is identically zero). If \(\breve{B}_i\) is the (unique) solution to (22), then \(B_{j,i}\!=\!V_{\ell ,1:r}\breve{B}_{i}\) is uniquely determined.

Proof

Given the (unique) solution \(\breve{B}_i\) to Eq. (22), recall that all possible solutions in the original parameter space have the form \(B_{C(i),i}=V_{1:r}\cdot \breve{B}_i+k_Y\), for some \(k_Y\in K_Y=\mathrm{{range}}\{V_{r+1:n_i}\}\). That is, \(B_{C(i),i}=V_{1:r}\cdot \breve{B}_i+V_{r+1:n_i}h\) for some (column) vector \(h\in \mathbb R ^{n_i-r}\). Suppose that a row index \(\ell \) exists such that \(V_{\ell ,r+1:n_i}\) is a null row. Clearly, for any \(h\), \(V_{\ell ,r+1:n_i}h=0\). It follows that the \(\ell \)th row of \(B_{C(i),i}\) is given by \(V_{\ell ,1:r}\cdot \breve{B}_i+V_{\ell ,r+1:n_i}h=V_{\ell ,1:r} \cdot \breve{B}_i\), i.e. it is uniquely determined (independent of \(h\)). Since the \(\ell \)th row of \(B_{C(i),i}\) corresponds to \(B_{j,i}\) by the definition of \(j\), the assertion is proven. \(\square \)

Reprints and permissions