Parameter estimation in tree graph metabolic networks

Parameter estimation in general networks

We first consider the parameter estimation problem in general linear time-invariant ordinary differential equation (LTI-ODE) systems. For convenience, we briefly sketch the approach when the catalytic rates are constants over time as in our previous work (Astola et al., 2011).

We recall that any network can be represented as a graph, where nodes are connected by edges when there is some interaction between these nodes. In a metabolic network a node represents a substrate or a product, and a directed edge from node i to node j means that i can be converted to j by enzymatic activity. To an edge from node i to j, we assign a weight, i.e., the catalytic rate k_ij ≥ 0, which represents the rate of product formation. In parameter inference one estimates the k_ij from data.

Denoting the concentration of substrate i at time t as X_i(t), a general time-invariant linear ODE model with a constant nonhomogeneous term, satisfying the mass conservation law, can be written as (1)${\displaystyle {\dot{X}}_i\left(t\right)=-\sum _{j\ne i}{k}_{ij}{X}_i\left(t\right)+\sum _{j\ne i}{k}_{ji}{X}_j\left(t\right)+b_i,}$ for i = 1, …, n, with (2)${\displaystyle b_i=\left\{\begin{array}{cc}\text{constant}\hfill & \text{if}i=1\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}\right.}$ The first summation in (1) stands for the edges leaving X_i, the second for the incoming edges and b_i for the possible in or outflow to the system. To simplify the notation, we introduce a matrix A with components given by (3)${\displaystyle \left\{\begin{array}{c}{A}_{ij}=k_{ji},i\ne j\hfill \\ A_{ii}=-\sum _{j\ne i}{k}_{ij}.\hfill \end{array}\right.}$ Then, (1) becomes (4)${\displaystyle {\dot{X}}_i\left(t\right)=\sum _{j=1}^n{A}_{ij}{X}_j\left(t\right)+b_i,i=1,\dots ,n.}$ Equation (4) can be rewritten in a compact matrix form as (5)${\displaystyle \dot{X}\left(t\right)=\left(\begin{array}{c}{\dot{X}}_1\left(t\right)\hfill \\ {\dot{X}}_2\left(t\right)\hfill \\ ⋮\hfill \\ {\dot{X}}_n\left(t\right)\hfill \end{array}\right)=\left(\begin{array}{cccc}-\sum _{j\ne 1}{k}_{1j}\hfill & \hfill & k_{n1}\hfill & b_1\hfill \\ k_{12}\hfill & \hfill & k_{n2}\hfill & 0\hfill \\ ⋮\hfill & \hfill & \hfill & ⋮\hfill \\ k_{1n}\hfill & \hfill & -\sum _{j\ne n}{k}_{nj}\hfill & 0\hfill \end{array}\right)\left(\begin{array}{c}{X}_1\left(t\right)\hfill \\ X_2\left(t\right)\hfill \\ ⋮\hfill \\ X_n\left(t\right)\hfill \\ 1\hfill \end{array}\right)=\tilde{A}\cdot \tilde{X}\left(t\right),}$ where $\tilde{X}\left(t\right)$ is obtained from X(t) by appending an extra 1 and matrix $\tilde{A}$ is obtained from A by extending it with an extra column containing the constant b₁.

To reconstruct a metabolic network from time-series measurements, we have to estimate the reaction rates k_ij, i.e., the weights of the edges in the network and the flow terms b_i. In view of (5), it is sufficient to estimate the (n + 1) × (n + 1) matrix $\tilde{A}$. We denote the data, i.e., measured concentrations of substrate i at time points t_j, j = 1, …, m, as an (n × m) matrix 𝕏. Estimates of the derivatives of the data curves we will store in a matrix $\dot{\mathbb{X}}$. To compute these estimates we may proceed in two ways. First, construct two n × m data matrices 𝕏₀, 𝕏₁ as follows (6)${\displaystyle {\mathbb{X}}_0=\left(\begin{array}{cccc}{\mathbb{X}}_{1,m-1}\hfill & {\mathbb{X}}_{1,m-2}\hfill & \dots \hfill & {\mathbb{X}}_{1,0}\hfill \\ {\mathbb{X}}_{2,m-1}\hfill & {\mathbb{X}}_{2,m-2}\hfill & \dots \hfill & {\mathbb{X}}_{2,0}\hfill \\ ⋮\hfill & \hfill & \dots \hfill & ⋮\hfill \\ {\mathbb{X}}_{n,m-1}\hfill & {\mathbb{X}}_{n,m-2}\hfill & \dots \hfill & {\mathbb{X}}_{n,0}\hfill \end{array}\right),{\mathbb{X}}_1=\left(\begin{array}{cccc}{\mathbb{X}}_{1,m}\hfill & {\mathbb{X}}_{1,m-1}\hfill & \dots \hfill & {\mathbb{X}}_{1,1}\hfill \\ {\mathbb{X}}_{2,m}\hfill & {\mathbb{X}}_{2,m-1}\hfill & \dots \hfill & {\mathbb{X}}_{2,1}\hfill \\ ⋮\hfill & \hfill & \dots \hfill & ⋮\hfill \\ {\mathbb{X}}_{n,m}\hfill & {\mathbb{X}}_{n,m-1}\hfill & \dots \hfill & {\mathbb{X}}_{n,1}\hfill \end{array}\right),}$ where m is the number of measurements. The matrix (7)${\displaystyle \dot{\mathbb{X}}\equiv \frac{1}{\mathrm{\Delta }t}\left({\mathbb{X}}_1-{\mathbb{X}}_0\right),}$ could then be used as an approximation for $\dot{X}$. For simplicity we assume the time grid to be equidistant with time step Δt. If this is not the case, the necessary modifications are easily implemented.

Secondly, we may use an alternative and often better approach to obtain approximations for ${\dot{\mathbb{X}}}_i$ by fitting splines to the time series data 𝕏_i (Zhan & Yeung, 2011). To obtain curves that interpolate the data faithfully, we require that the distances between the curves and the measurements are minimal and that at the same time the curves are smooth. To achieve this we fit P-splines, which are B-splines with a penalization for non-smoothness (Eilers & Marx, 1996). From these splines, we evaluate the derivative estimates at time points t_j. These estimates are then used as entries in the matrix $\dot{\mathbb{X}}$. Having at hand an estimate for matrix $\dot{\mathbb{X}}$, the problem of network inference comes down to finding the matrix $\tilde{A}$ from the equation (8)${\displaystyle \dot{\mathbb{X}}=\tilde{A}\tilde{\mathbb{X}},}$ in which $\dot{\mathbb{X}}$ is known and $\tilde{\mathbb{X}}$ is obtained from the data matrix 𝕏 by extending this with an extra row of ones. However, solving $\tilde{A}$ directly from (8) often results in over-fitting, since all possible edges are included in the modeled network. Another serious shortcoming of such a matrix (pseudo-) inversion approach is the fact that we cannot control the positivity of the reaction rates. Although in Schmidt, Cho & Jacobsen (2005) negative coefficients were interpreted as inhibition of the compounds, in many biological pathways, negative coefficients are not permitted. Thus we take a more general approach in which one can exclude all edges that are biologically not acceptable, and in which one can constrain the reaction rates to be positive, without substantially compromising computation time.

To this end, we reformulate the equation as a minimization problem: (9)${\displaystyle \underset{\tilde{A}}{argmin}\left(\parallel \dot{\mathbb{X}}-\tilde{A}\tilde{\mathbb{X}}\parallel \right).}$ The matrix norm used here is the Frobenius norm: (10)${\displaystyle \parallel \tilde{A}\parallel =\sqrt{\sum _{i=1}^n\sum _{j=1}^m{\tilde{A}}_{ij}^2}.}$ This alternative formulation allows inclusion of expert knowledge in a simple way. We put ${\tilde{A}}_{ij}=0$, when an edge from node i to node j cannot exist. Nearly all mathematical software packages (Mathematica, Matlab, Maple, etc.) can numerically find the minimizer $\tilde{A}$ (and thus the reaction rates k_ij and the flow term b₁) with the constraint that k_ij ≥ 0.

Parameter estimation in tree networks

As described in the introduction, tree networks are networks, whose graphs resemble trees in that they branch away from the root and the directions of the edges always point from the root towards the leaves. In Fig. 1 we presented, using an example, the difference between a tree and a non-tree graph. In a kinetic reaction system with a tree network, the parameters can be uniquely estimated even when they are time dependent. We could write this down in general. However, the proof is based on one central idea. We feel that the reader gains more insight if we simply show this idea through an example. To that end we use as example the network in the left hand side of Fig. 1. The extension to the general is straightforward.

For the network on Fig. 1A, we have the following kinetic mass balance model: (11)${\displaystyle \left(\begin{array}{c}{\dot{X}}_1\hfill \\ {\dot{X}}_2\hfill \\ {\dot{X}}_3\hfill \\ {\dot{X}}_4\hfill \\ {\dot{X}}_5\hfill \\ {\dot{X}}_6\hfill \\ {\dot{X}}_7\hfill \end{array}\right)=\left(\begin{array}{ccccccccc}-\left(k_{1,2}+k_{1,3}\right)\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & b_1\hfill \\ k_{1,2}\hfill & -\left(k_{2,4}+k_{2,5}\right)\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ k_{1,3}\hfill & 0\hfill & -\left(k_{3,6}+k_{3,7}\right)\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & k_{2,4}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & k_{2,5}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & k_{3,6}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & k_{3,7}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right)\left(\begin{array}{c}{X}_1\hfill \\ X_2\hfill \\ X_3\hfill \\ X_4\hfill \\ X_5\hfill \\ X_6\hfill \\ X_7\hfill \\ 1\hfill \end{array}\right),}$ where the constant b₁ represents the influx into the system and the k_i,j are the catalytic rates. Note that there are as many unknown parameters (k_i,j, b₁) as there are measured variables X_i(t_j). Therefore, as can be directly verified, we can rewrite the previous matrix equation by exchanging the X_i and k_i,j as follows: (12)${\displaystyle \left(\begin{array}{c}{\dot{X}}_1\hfill \\ {\dot{X}}_2\hfill \\ {\dot{X}}_3\hfill \\ {\dot{X}}_4\hfill \\ {\dot{X}}_5\hfill \\ {\dot{X}}_6\hfill \\ {\dot{X}}_7\hfill \end{array}\right)=\underset{\text{matrix}B}{\underbrace{\left(\begin{array}{ccccccc}1\hfill & -X_1\hfill & -X_1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & X_1\hfill & 0\hfill & -X_2\hfill & -X_2\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & X_1\hfill & 0\hfill & 0\hfill & -X_3\hfill & -X_3\hfill \\ 0\hfill & 0\hfill & 0\hfill & X_2\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & X_2\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & X_3\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & X_3\hfill \end{array}\right)}}\left(\begin{array}{c}{b}_1\hfill \\ k_{1,2}\hfill \\ k_{1,3}\hfill \\ k_{2,4}\hfill \\ k_{2,5}\hfill \\ k_{3,6}\hfill \\ k_{3,7}\hfill \end{array}\right).}$

We immediately see that B is an upper triangular matrix since the entries below the diagonal are zero. This implies that the determinant of the matrix B in (12) is the product of the entries on the diagonal: $X_1^2\cdot X_2^2\cdot X_3^2$, and thus unequal to 0 since X_i ≠ 0, ∀i = 1, …, n. So, B is invertible and the system of equations has the unique solution. (13)${\displaystyle \left(\begin{array}{c}{b}_1\hfill \\ k_{1,2}\hfill \\ k_{1,3}\hfill \\ k_{2,4}\hfill \\ k_{2,5}\hfill \\ k_{3,6}\hfill \\ k_{3,7}\hfill \end{array}\right)=\underset{\text{matrix}{B}^{-1}}{\underbrace{\left(\begin{array}{ccccccc}{X}_0^{-1}\hfill & X_0^{-1}\hfill & X_0^{-1}\hfill & X_0^{-1}\hfill & X_0^{-1}\hfill & X_0^{-1}\hfill & X_0^{-1}\hfill \\ 0\hfill & X_1^{-1}\hfill & 0\hfill & X_1^{-1}\hfill & X_1^{-1}\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & X_1^{-1}\hfill & 0\hfill & 0\hfill & X_1^{-1}\hfill & X_1^{-1}\hfill \\ 0\hfill & 0\hfill & 0\hfill & X_2^{-1}\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & X_2^{-1}\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & X_3^{-1}\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & X_3^{-1}\hfill \end{array}\right)}}\left(\begin{array}{c}{\dot{X}}_1\hfill \\ {\dot{X}}_2\hfill \\ {\dot{X}}_3\hfill \\ {\dot{X}}_4\hfill \\ {\dot{X}}_5\hfill \\ {\dot{X}}_6\hfill \\ {\dot{X}}_7\hfill \end{array}\right).}$

Time varying kinetic rates

In earlier work we developed a fast method to reconstruct metabolic networks (Astola et al., 2011). The idea in this approach was to substitute the measurements directly into the model equations and not only in the objective function. This approach had as a limitation that all parameters were assumed to be constant in time. Here we extend our previous approach by allowing the catalytic rates to be time dependent, to better reflect the real situation, since in practice the enzyme concentrations are fluctuating in time. This has also immediately resulted in reconstructions that better fit the observed data as can be seen in Fig. 2. While the standard practice in enzyme kinetics is to either use constant catalytic rates in mass balance equation or to model product formation through a Hill function (Goutelle et al., 2008) such as in the Michaelis–Menten equation (Savageau, 1995), none of these take into account the fact that the enzyme concentration is also changing in time. Since we also want to study the relation of gene expression and enzyme concentration in time, we need to capture their dynamics.

As the catalytic rate is now modeled as a function in time, and not as a constant, it is no longer possible to infer this with the standard procedure of solving for those parameters that fit the ordinary differential equations to data in the sense of maximum likelihood. We cannot clearly separate the substrate/product, enzyme concentrations and noise, since we have no measurements of the enzyme concentrations. To solve them, we would have to impose a model on them, which we don’t have a priori. A reasonable approach in this situation is to first estimate a model for the metabolite concentrations for which we have several measurements. By fixing the concentrations first using spline approximations, we may then estimate the trends in the enzyme concentrations. This method assumes that the solutions are rather smooth. If this is not the case and the sampling frequency is low, the derivatives obtained by fitting splines can introduce errors that distort the reconstruction. The inference method proposed here is by no means restricted to tree networks, but in case the network has a tree structure, the parameters can be estimated in an unambiguous way. We summarize the general work flow for the proposed parameter inference in the schematic diagram in Fig. 3.

Time dependent parameter estimation

In this section we present three different schemes to estimate the k_ij(t) in model (4). In (9) we used the data at all time points simultaneously to estimate the time independent parameters. However, a remarkable feature of tree structured networks is that the data at one time point is already enough to calculate unique estimates for the parameter values at that particular time point. This is immediately clear from (13): as soon as we have estimates for the time derivatives $\dot{X}\left(t_k\right)$ available, we may calculate estimates for the k_ij(t_k).

• To estimate the derivatives at some time point one still needs the data of neighboring time points. So, the first step in this scheme is to fit, e.g., P-splines to the data time series (O’Sullivan, 1986; Eilers & Marx, 1996). From these splines we calculate estimates for the time derivatives ${\dot{X}}_i\left(t_k\right)$. Then by substituting these estimates as well as the measurements into Eq. (4), we are left with a set of linear equations to solve k_ij(t_k) and b₁ at all times t_k. Finally, for smooth and continuous catalytic rates, one may fit, e.g., a second order polynomial through these estimates.

• An alternative approach in which the number of parameters is smaller than in scheme 1, is to assume that the functions k_ij(t) can be adequately represented as polynomials in time of some order. In practice order 2 is often sufficient. With this choice we have then: (14)${\displaystyle k_{ij}\left(t\right)={\alpha }_{ij}{t}^2+{\beta }_{ij}t+{\gamma }_{ij}.}$

This implies that per k_ij we have 3 parameters to be estimated using the whole time series data. By substituting (14) into matrix $\tilde{A}$ in (9) we then obtain estimates for α_ij, β_ij and γ_ij, and thus for k_ij(t).

• As in the previous scheme, we assume (14). We construct an objective function like the following: (15)${\displaystyle \sum _k\parallel X\left(t_k\right)-\mathbb{X}\left(t_k\right)\parallel ,}$ which is the sum of the distances between X(t_k) and the measurements. We look for a matrix $\tilde{A}$, such that the solutions X_i(t) to (4) minimize this objective function. Using suitable optimization algorithm we simultaneously estimate X_i, k_ij, and b₁.

To compare the fit, accuracy and speed of these three schemes we applied them using as test networks random tree networks that have equal numbers of nodes and edges as the network on Fig. 1A.

In these networks, we simulated time series data with time varying catalytic rates. To generate artificial data, we assigned random values to α_ij, β_ij and γ_ij in a range, such that the resulting solutions have approximately the same range as the metabolite concentration data for quercetin glycosides measured in tomato seedlings (cf. Fig. 2). To assess the reconstruction power of the three schemes, we also tested them on networks that are not trees. The corresponding data generation process is the same but the network models contain cycles. In the third set of simulations we added ±10% uniformly distributed noise to tree structured network data.

Parameter inference as a mean to select active genes

In addition, as a potentially powerful application, we show how we may infer the gene candidates likely to be involved in the enzymatic reactions. This can be done by comparing estimated time dependent catalytic rates with simultaneously measured gene expression data. If, according to the model, the formation of a metabolite necessitates higher/lower enzyme concentration, this should be also observable in the expression level of the gene that codes for this enzyme. Using this heuristic, we were able to select from a large set of potential genes the most likely candidate genes for further experimental validation of their functioning in particular reactions. In view of this application, small inaccuracies in parameters are not detrimental, since here we are mainly interested the dynamic trends of the catalytic rates instead of their precise numeric values.

As an example, we take the quercetin glycosylation pathway in cotyledons, occurring during the development of tomato seedlings (Koes, Quattrocchio & Mol, 1994). Quercetin glycosides are a subset of flavonoids, which are plant secondary metabolites naturally produced by plants. Flavonoids are being intensively studied for their proposed beneficial effects on prevention of chronic diseases (Bovy, Schijlen & Hall, 2007; Rein et al., 2006; Moon, Wang & Morris, 2006).

We have measured the concentrations of several quercetin derivative compounds accumulating in cotyledon- and hypocotyl tissues. We have daily measurements from day 5 after sowing up to day 9. The same sample used for the metabolite analysis with liquid chromatography mass spectrometry were used for gene expression analysis. The expression levels of genes, putatively involved in the glycosylation of quercetin, were quantified using microarray analysis. Glycosyltransferases (GTs) are members of the multigene superfamily in plants that can transfer single or multiple sugars to various plant molecules, resulting in the glycosylation of these compounds (Wang, 2009). To date, it is not known exactly which GTs catalyze each glycosylation reaction. With more than 200 GT candidates an experimental validation of every single GT is costly. Therefore we wanted to make a pre-selection of the potentially strongest gene candidates, using mathematical modeling and simulations. We use the heuristics that if the kinetic ODE model describes the system of enzymatic reactions reasonably well, the estimated catalytic rates should reflect the real enzymatic activity. This in turn should correlate with the expression trends of the GTs observed using the time series microarray analysis.

Our procedure for the GT inference is as follows:

• 1.

  Given the time series metabolite concentration data, estimate the time dependent parameters using all biologically relevant networks. Select the network that gives the best fit to measurements with respect to residual or goodness of fit etc. Save the estimated catalytic rates corresponding to the best network as reference.

• 2.

  Compute correlations between the time series of expression levels of each GT and the previously saved series of catalytic rates.

• 3.

  Select those GTs whose dynamics correlate best with catalytic dynamics for further experimental validation.

Comparison of parameter inference schemes

As can be seen from Figs. 4A, 4D and 4G, concerning the fitting errors, all schemes give similar results and their box-plots have some overlap. In principle they are solving the same optimization problem, only scheme 1 first solves the point wise rate values and then fits a polynomial, whereas scheme 2 searches for a polynomial-valued rates that fit to the whole series of data and scheme 3 tries simultaneously estimate the parameters as well as the derivatives. We measured the accuracy of the parameter estimation by computing the Frobenius norm (10) of the difference between the original timevariant kinetic rates used in simulation and the reconstructed rates. Besides the actual estimation accuracy, also computation times are relevant. In terms of computation time, scheme 2 is the fastest and scheme 3 is slowest, although the differences are not large. Notice that the comparisons in Fig. 4 were done in a setting where equal parameter constraints (k_ij > 0) were given to the solvers and the parameters were estimated using constrained non-linear global optimization (NMinimize in Mathematica) choosing for the fast Nelder–Mead algorithm (with option “PostProcess” → False).

This result is more or less to be expected, since when the data is reasonably accurate, it does not always make sense to re-estimate the data by using it as an unknown variable in the equations of the system. Rather, it may pay off to substitute the data directly into the equations reducing the number of unknown elements. Also it is logical that schemes 1 and 2 perform less well on non-tree graph networks, since the assumption on unique point-wise estimability is not valid anymore. Since our method is based on initial fitting of splines, the major sensitivity is indeed with respect to data. This was also confirmed by the sensitivity analysis we conducted.

Our network models, although relatively small, belong to the general group of the so-called sloppy biochemical models (Gutenkunst, 2008), despite of which the parameters still may be identifiable. For a separate discussion and more background on this subject, please see ‘Discussion’. The range of eigenvalues of the Hessian of the residual (between predicted and measured values) varies from 10⁻⁴ to 10⁵. For the sensitivity analysis numerical derivatives need to be computed. Since we are considering time varying parameters, we have taken time-averages of point-wise derivatives. Eigenvectors corresponding to very small eigenvalues, implying sloppiness in sensitivity, all point towards those parameters that are associated with network nodes where the measured metabolite concentrations are very low. This is logical since the parameters associated with concentration values close to zero have little effect on the residual, because our objective function does not contain the standard deviation term in the denominator. By this choice we explicitly wanted to avoid that those measurements that are close to noise level shall have equal weight with the more abundant ones.

Enzyme inference from microarray data

In Fig. 5 we illustrate the results of the analysis as described in ‘Parameter inference as a mean to select active genes’. These are the expression levels of best matching six GTs together with the estimated catalytic rates for the reactions that corresponds to the conversions from node X_i to X_j exactly as in Fig. 1A. We have standardized, i.e., subtracted the mean and divided by standard deviation both predicted and measured expressions for visual comparison. As can be seen from Fig. 5, the deviation of the expression levels between samples can vary from gene to gene. One could also weight the correlation according to this variation so that more precise observations are favored. For accurate reconstruction of both the kinetic rates as well as the selection of appropriate genes, a time series with more data points is desired. What exactly the minimal sample number and sampling method should be depends on the data and the system model, but a rule of thumb from experienced modelers would be a minimum of 15 data points. To test experimentally whether the inferred genes are actually related to the enzymes that glycosylate the flavonols, a set of selected genes are currently being cloned.

As a computational validation of the selection procedure, we tested whether substituting the (scaled) expression levels of the selected genes into the model will result in a decreased residual (better likelihood of observing the measurements). The reason we want to do this post-analysis is two-fold. First of all, our GT candidates are ranked according to their correlation with the predicted enzymatic trends, but it may happen that several candidates have almost equal correlation coefficients. This makes it difficult to distinguish between the candidates, especially because the initial GT-population is already a result of an ontology-based selection. Another point is that the selection of the most likely GT’s is based on individual matchings with single dynamic parameters whose magnitudes are unknown. It is not absolutely clear, say, whether the combination of the very best candidates will always give better results than when, for example, one candidate is actually the second best one (in terms of correlation). In each network combination, at most seven GTs are considered, but the number of all possible combinations is still very high. Also the expression levels need to be scaled to match the metabolic model.

To ensure a rich set of gene combinations, we ran a Markov Chain Monte Carlo-algorithm (MCMC) (Calvetti & Somersalo, 2007). To address the question of whether the differences in correlations are significant enough, we first ordered the genes into a sequence according to their correlation with the predicted enzyme concentration levels and took two sets of genes according to their order number in the sequence: 1, 2, …, 10 and 11, 12, …, 20. We tested whether the residuals, obtained after 200 iterations of 1,000 samples with MCMC algorithm using the data of these two sets, have equal means and variances. For the mean test we obtained a P-value less than 0.00001 and for the variance test a P-value of less than 0.006. We may conclude that in the context of a dynamic kinetic reaction model, those genes with expression levels highly correlating to the predicted enzyme dynamics, are significantly more likely to be responsible for the observations.

Structural identifiability

A general problem in parameter estimation is that it is difficult and sometimes even impossible to be sure that the estimated parameters are unique. If the model is structurally unidentifiable, there is an infinite number of parameter sets that give equal results. This is a substantial challenge, especially when the network structure is not known, since an overly complex network can result in over-fitting. This problem is not present in any of the (biologically) potential networks as sketched in Fig. 6, since as tree graphs these all turn out to be locally structurally identifiable as they can be embedded in an upper triangular matrix as discussed in the preceding section.

Practical identifiability

Structural identifiability does not imply practical identifiability and therefore we have studied the practical identifiability of the parameters in our system by means of profile likelihood (Raue et al., 2009). We learned that all the kinetic parameters connecting substrates and products with concentrations above detection limit show also practical identifiability (see Supplementary Data). Another observation is that if we allow a product to decay without constraints, the practical identifiability as well as the tree structure of the graph is lost.