Parameter Estimation and Identifiability in Kinetic Flux Profiling Models of Metabolism

The method of Kinetic Flux Profiling includes experimental data collection, conversion of a network diagram to a model consisting of a system of ordinary differential equations and algebraic constraints, and a computational step for the estimation of parameters from the data. The experiment begins by allowing the cells to reach a metabolic steady state in unlabeled media. At the metabolic steady state, all fluxes and all metabolite concentrations are constant. Next, the unlabeled media is switched out for its stable-isotope-labeled equivalent. In most cases, one specific metabolite in the media will be labeled with the carbon isotope ¹³C in place of all its carbons while the rest of the media remains unchanged. After the switch to isotope-labeled media, the ¹³C carbons will be incorporated into metabolites within the network. Samples from the cells are taken at multiple time points, and the proportions of labeled and unlabeled metabolites are quantified using MS or NMR. Ideally, these steps are performed without disrupting the metabolic steady state. That is, without perturbing the metabolic concentrations or fluxes.

Next, the architecture of the metabolic pathway under study is used to write a system of ordinary differential equations that describe the change in isotope labeling. The metabolic pathway is represented as a directed weighted graph where the nodes represent metabolites and the edges represent reactions. The pathway also includes additional inputs and outputs, edges that do not begin or do not terminate at a node. To ensure non-trivial solutions, the pathway will contain at least one labeled input and at least one output. For simplicity, the diagram should only include metabolites that may become labeled. That is, there must be a path from a labeled input to each metabolite in the diagram. In a metabolic steady state the concentrations of the metabolites are not changing. Therefore the total flux into and out of each metabolite must be equal. This flux balance condition gives an additional set of algebraic constraints on the fluxes.

The resulting system contains parameters for each flux in the network as well as for the total concentration of each metabolite. In a standard KFP protocol, it is presumed that these total concentrations are all known. The final step is to use a parameter fitting method to fit the parameters for the metabolic fluxes to the isotope labeling time courses. Since the ODE system is linear, if we had access to perfect, noise-free, continuous time data the parameters could, in theory, be computed. However, we address the accuracy and reliability of the parameter estimates when they are made with limited and noisy data. In the following sections, we seek to answer three crucial questions about the estimation of fluxes in KFP, motivated by realistic limitations on data availability. What can we learn without estimates of the total concentrations? What can we learn with only steady-state labeled proportions? What can we learn when there is a separation of time scales between metabolites?

We begin by introducing a general method for formulating the ODE system and flux parameter constraints. We next suggest a scaling of the system that does not include explicit dependence on the total concentrations of each metabolite. We finish section 2 with a condition on the graph that allows us to estimate relative fluxes (proportion of each flux to the total through a metabolite) and turnover rates of each metabolite without measuring metabolite concentrations. In section 3, we investigate the steady-state problem. Often metabolic changes occur on a time scale that make it technically challenging to collect sufficient data before isotopic steady state is reached. While it is not possible, in general, to compute turnover rates from steady-state data, we give a simple necessary condition for recovering the relative fluxes from steady-state only. This is illustrated with two examples using Bayesian parameter estimation on simulated data to show both the accuracy and reliability of the estimates. Finally, in section 4 we address what can happen if one of the turnover rates is substantially faster than the others. We show that the rapid turnover rate parameter is not identifiable unless the fast metabolite is the one that receives labeled input. Again, this is illustrated using Bayesian parameter estimation on simulated data from specific examples.

Suppose the diagram of the metabolic pathway is represented by the graph $G$. Let $G=(\mathcal{N},\mathcal{E})$ be a special directed weighted graph where the set of nodes $\mathcal{N}$ represent metabolites and the set of edges $\mathcal{E}$ have weights that represent metabolic flux in the direction of the edge. The set $\mathcal{E}$ contains edges that do not originate or terminate at a node (i.e. they enter and exit the scope of the graph). It is therefore convenient to partition $\mathcal{E}$ into four disjoint subsets

where $\mathcal{E}_{\mathcal{L}}$ is the set of edges that provide labeled metabolite to the pathway, $\mathcal{E}_{\mathcal{U}}$ is the set of edges that provide unlabeled metabolite to the pathway, $\mathcal{E}_{\mathcal{V}}$ is the set of edges that exit the scope of the pathway, and $\mathcal{E}_{\mathcal{W}}$ are the edges that connect nodes within the pathway. We will consider only connected graphs in which there is a path from the labeled input(s) to every node. This is because any node that can not receive label will not provide any additional information. In some cases where a system has edges exiting or entering, another node is added to be the target or origin of those edges. We have chosen not to do that in this case both because it would require separate handling of inputs with labeled and unlabeled metabolites and because we feel the notion of entering and exiting the pathway is useful for the interpretation of the experiments.

We will use the following notation for the number of metabolites and reactions in the pathway, $|\mathcal{N}|=N$ and $|\mathcal{E}|=R$. Because we only consider the case where there is at least one labeled influx $1\leq|\mathcal{E}_{\mathcal{L}}|\leq N$. Each metabolite may have an unlabeled influx and there must be at least one outflux so $0\leq|\mathcal{E}_{\mathcal{U}}|\leq N$ and $1\leq|\mathcal{E}_{\mathcal{V}}|\leq N$. Finally, since the graph is connected, $N-1\leq|\mathcal{E}_{\mathcal{W}}|\leq N(N-1)$. In Figure 1 we have $N=3$ metabolites and $R=10$ reactions. $\mathcal{E}_{\mathcal{L}}$ corresponds to the labeled inflow with rate $f_{1}$. $\mathcal{E}_{\mathcal{U}}$ corresponds to the unlabeled inflows with rates $f_{2},f_{5}$ and $f_{8}$. $\mathcal{E}_{\mathcal{V}}$ corresponds to the outflows with rates $f_{3},f_{6}$ and $f_{9}$. $\mathcal{E}_{\mathcal{W}}$ corresponds to the flows between nodes with rates $f_{4},f_{7}$ and $f_{10}$.

By considering only the directed weighted graph with edges in $\mathcal{E}_{\mathcal{W}}$ we can define the weighted adjacency matrix. The weighted adjacency matrix, $\boldsymbol{W}\in\mathbb{R}^{N\times N}_{\geq 0}$, is the matrix composed of the weights of the directed edges between nodes, $\mathcal{E}_{\mathcal{W}}$. The matrix entry $w_{i,j}$ is the weight of the edge from node $i$ to node $j$. If no connection between node $i$ and node $j$ exists $w_{ij}=0$ [13]. We further define the total concentration matrix $\boldsymbol{X}^{\mathcal{T}}$ as the diagonal matrix containing the total concentration of each metabolite ${[\boldsymbol{X}^{\mathcal{T}}]}_{i,i}=x_{i}^{\mathcal{T}}$. Then the terms in eq. 4 corresponding to the influx into each node from the other nodes is $\boldsymbol{W}^{T}{\boldsymbol{X}^{\mathcal{T}}}^{-1}\boldsymbol{X}^{\mathcal{U}}$. In our example,

Another commonly used matrix in graph theory is the degree matrix. The degree matrix, $\boldsymbol{D}$, is a diagonal matrix of the degree of each node [12]. Since we have a directed weighted graph, we will define two variations of the degree matrix: the weighted out-degree matrix and the weighted in-degree matrix.

The terms in eq. 4 that correspond to outfluxes to other metabolites can be written. $\boldsymbol{D}_{\text{out}}{\boldsymbol{X}^{\mathcal{T}}}^{-1}\boldsymbol{X}^{\mathcal{U}}$. In our example,

Our graph $G$ also includes edges in sets $\mathcal{E}_{\mathcal{L}}$, $\mathcal{E}_{\mathcal{U}}$, and $\mathcal{E}_{\mathcal{V}}$ that enter and exit the scope of the graph so we define matrices corresponding to these sets of fluxes.

The terms in eq. 4 which correspond to outfluxes exiting the graph can be written $\boldsymbol{D}_{\mathcal{V}}{\boldsymbol{X}^{\mathcal{T}}}^{-1}\boldsymbol{X}^{\mathcal{U}}$. In our example,

The last set of terms in eq. 4 are the constant terms for the unlabeled inputs and are the entries on the diagonal of $\boldsymbol{D}_{\mathcal{U}}$. Therefore, we can define the vector $\hat{\boldsymbol{b}}$ in eq. 4 as having entries $\hat{b}_{i}={[\boldsymbol{D}_{\mathcal{U}}]}_{i,i}$.

Finally, we can define $\boldsymbol{\hat{A}}$ in eq. 4 and $\boldsymbol{F}$ in eq. 5. The matrix $\boldsymbol{\hat{A}}$ is defined by the equation

where the diagonal matrix $\boldsymbol{F}_{\text{out}}$ is defined by the equation $\boldsymbol{F}_{\text{out}}=\boldsymbol{D}_{\text{out}}+\boldsymbol{D}_{\mathcal{V}}$ and gives the total rate of flux out of each metabolite. Similarly, we define the matrix $\boldsymbol{F}_{\text{in}}$ to be the diagonal matrix of all fluxes entering each metabolite, $\boldsymbol{F}_{\text{in}}=\boldsymbol{D}_{\text{in}}+\boldsymbol{D}_{\mathcal{L}}+\boldsymbol{D}_{\mathcal{U}}$. As stated above, the algebraic constraints for constant total metabolite concentrations is that $\boldsymbol{F}_{\text{in}}=\boldsymbol{F}_{\text{out}}$. When these constraints are satisfied, we call this diagonal matrix $\boldsymbol{F}$ for simplicity. The diagonal entries of $\boldsymbol{F}$ will be denoted $F_{i}$ for the flux through the $i^{th}$ metabolite.

Measurements of labeled and unlabeled metabolite concentrations are typically reported as a “percent enrichment”, the percentage of the total pool of a metabolite that is labeled with at least one ¹³C. This is because the mass spectrometer measures a signal intensity at specific masses. The intensity of the signal at the mass of the labeled metabolite is divided by the total signal intensity for the metabolite to compute a labeled relative signal intensity. Conversion of signal intensity to concentration requires extensive additional experiments to develop calibration curves for each metabolite in the experiment. However, the relative signal intensity for metabolite $i$ is a reasonable approximation of the relative concentration for metabolite $i$ whether labeled ($x^{\mathcal{L}}_{i}/x^{\mathcal{T}}_{i}$) or unlabeled ($x^{\mathcal{U}}_{i}/x^{\mathcal{T}}_{i}$). Therefore, we scale the model so that the new variables are the proportions of each metabolite that are unlabeled, $\bar{x}_{i}=x^{\mathcal{U}}_{i}/x^{\mathcal{T}}_{i}$ or in matrix form

In this scaled format, it is useful to also write these equations in terms of new parameters. Let us define the vector $\vec{\boldsymbol{\alpha}}=\boldsymbol{F}_{\text{in}}^{-1}\hat{\boldsymbol{b}}$ and the matrices $\boldsymbol{B}=\boldsymbol{F}_{\text{in}}^{-1}\boldsymbol{W}^{T}$, $\boldsymbol{K}_{\text{in}}={\boldsymbol{X}^{\mathcal{T}}}^{-1}\boldsymbol{F}_{\text{in}}$, and $\boldsymbol{K}_{\text{out}}={\boldsymbol{X}^{\mathcal{T}}}^{-1}\boldsymbol{F}_{\text{out}}$. The entries of $\vec{\boldsymbol{\alpha}}$ are $\alpha_{i}=b_{i}/{[\boldsymbol{F}_{\text{in}}]}_{i,i}$, the proportion of the total flux into metabolite $i$ that comes from an unlabeled input originating outside of the diagram. Similarly, for the fluxes originating from other metabolites, the entries of $\boldsymbol{B}$ are $\beta_{i,j}=w_{j,i}/{[\boldsymbol{F}_{\text{in}}]}_{i,i}$, the proportion of the total flux into metabolite $i$ that comes from metabolite $j$. This new notation for the fluxes has been incorporated into a diagram of the network for our example and is shown in fig. 2. The entries of $\vec{\boldsymbol{\alpha}}$ and $\boldsymbol{B}$ are all the dimensionless proportions of fluxes into each metabolite. Finally, $\boldsymbol{K}_{\text{in}}$ and $\boldsymbol{K}_{\text{out}}$ are diagonal matrices with entries that can be viewed as turnover rates of the metabolites in units $1/$time. The scaled system can be written

As before, if the metabolite concentrations are unchanging then we have the algebraic constraints $\boldsymbol{F}_{\text{in}}=\boldsymbol{F}_{\text{out}}=\boldsymbol{F}$ which gives $\boldsymbol{K}_{\text{in}}=\boldsymbol{K}_{\text{out}}=\boldsymbol{K}$. This simplifies the system further to

where $\vec{\boldsymbol{1}}$ is the vector with all entries equal to one.

This formulation of the problem is powerful because it allows us to write the model variables entirely in terms of proportions unlabeled. This proportion can be measured directly as one minus the isotopic enrichment. Additionally, the parameters in $\vec{\boldsymbol{\alpha}}$ and $\boldsymbol{B}$ are also ratiometric. They are proportions of each flux to total flux through the metabolites so we will refer to them as flux ratio parameters. The other parameters in $\boldsymbol{K}$ are the turnover rates for each metabolite. There is no explicit dependence on the total concentrations. This means we can construct the model and fit the parameters in this format without measuring metabolite concentrations. To convert the resulting values for the parameters back to fluxes in units of concentration per unit of time would still require the costly measurement of total concentrations. However, in many instances functional changes in metabolic networks could be gleaned from changes in the turnover rates and relative flux values which are the parameters in equation eq. 9.

Since the computational time for parameter estimation can quickly become impractical, it is essential to reduce the number of parameters as much as possible. In this section, we use the algebraic constraints to find the number of independent parameters which must be fit. We begin with a general statement of the number of parameters needed. Next, we investigate the special case where each metabolite has a non-zero flux exiting the graph. Finally, we discuss the more difficult case, where some exiting fluxes are missing from the reaction diagram.

Before we turn to the special case, let us begin by counting the parameters in our scaled system eq. 9. This model has $N$ parameters corresponding to the non-zero entries of $\boldsymbol{K}$, $k_{i}$ for the turnover rate of each metabolite, $|\mathcal{E}_{\mathcal{U}}|$ parameters for the non-zero entries of $\vec{\boldsymbol{\alpha}}$ for the unlabeled inputs, and $|\mathcal{E}_{\mathcal{W}}|$ parameters corresponding to the non-zero entries of $\boldsymbol{B}$ for the fluxes between metabolites. Since $\alpha_{i}$ and $\beta_{i,j}$ are all of the dimensionless proportions of fluxes into metabolite $i$, for any metabolite which has no labeled inputs, we know that

We can therefore, without loss of generality, replace the first non-zero entry in these rows

There are therefore $|\mathcal{E}_{\mathcal{W}}|-(N-|\mathcal{E}_{\mathcal{L}}|)$ remaining independent $\beta$ values. In our example, this means we can replace $\beta_{2,1}$ with $(1-\alpha_{2})$ and $\beta_{3,2}$ with $(1-\alpha_{3})$ so we have written both forms in fig. 2. Now we can write our system in terms of the seven remaining parameters.

A graph which is an arborescence will have $|\mathcal{E}_{\mathcal{L}}|=1$ and $|\mathcal{E}_{\mathcal{W}}|=N-1$. Therefore, such a graph will have no remaining independent $\beta$ parameters. On the other hand, a graph that is not an arborescence will have either $|\mathcal{E}_{\mathcal{L}}|>1$ or $|\mathcal{E}_{\mathcal{W}}|>N-1$ or both. Therefore, any graph which is not an arborescence will have non-zero $\beta$ parameters in this scaling.

So far we have counted the turnover rate parameters, $k_{i}$, and the proportional flux parameters, $\alpha_{i}$ and $\beta_{i,j}$. This brings us to a total number of parameters in eq. 9 of $N+|\mathcal{E}_{\mathcal{U}}|+|\mathcal{E}_{\mathcal{W}}|-(N-|\mathcal{E}_{\mathcal{L}}|)=R-|\mathcal{E}_{\mathcal{V}}|$. In developing this parameter set, we have labeled the edges in $\mathcal{E}_{\mathcal{U}}$ with $\alpha$s and those in $\mathcal{E}_{\mathcal{W}}$ with $\beta$s. Let us now turn our attention to $\mathcal{E}_{\mathcal{L}}$ and $\mathcal{E}_{\mathcal{V}}$. Because of our choice to scale all inputs to the total input for each metabolite, a labeled input edge in $\mathcal{E}_{\mathcal{L}}$ can also scale to the influx into its node $i$ as $(1-\alpha_{i}-\Sigma_{j=1}^{N}\beta_{i,j})F_{i}$. In our example, this means that the labeled input will be $f_{1}=(1-\alpha_{1}-\beta_{1,3})F_{1}$. In practice, it is typical to have only one labeled input so each of our examples has one edge in $\mathcal{E}_{\mathcal{L}}$ which is an input into the first metabolite.

The labeling of the edges in $\mathcal{E}_{\mathcal{V}}$ may pose a new complication because of our algebraic constraints. Notice first that in the special case $|\mathcal{E}_{\mathcal{V}}|=N$, the number of parameters already defined after the substitution in eq. 11 will be exactly the minimum number, $R-N$. This is because in this case, every node has an output vector that can be used to satisfy the algebraic constraint that the total output from each node must equal the flux through that node. Since the total of all fluxes out of node $i$ must be $F_{i}$, the flux from metabolite $i$ exiting the diagram will be $F_{i}-\sum_{j\neq i}\beta_{j,i}F_{j}$. In our example, this means that $f_{3}=F_{1}-(1-\alpha_{2})F_{2}$, $f_{6}=F_{2}-(1-\alpha_{3})F_{3}$, and $f_{10}=F_{3}-\beta_{1,3}F_{1}$. Since every node in the example has an output, we have completed our parameterization in a way that gives the minimum number of parameters and does not depend on any concentrations which would require additional measurements.

We will finish this section with a discussion about the case where $\mathcal{E}_{\mathcal{V}}<N$. In this case, there is at least one metabolite that does not have an outflux leaving the pathway. This means that we still have too many parameters in our formulation. However, our algebraic constraints now imply that the equation for the missing outflux must be zero. That is, for any node $i$ which does not have a flux out of the pathway, $F_{i}=\Sigma_{j\neq i}\beta_{j,i}F_{j}$. This constraint is challenging because it is a relationship among the total fluxes through multiple metabolites. It can be solved for one of the $\beta$ parameters or one of the $k$ parameters, however, the resulting equation will depend on the total concentrations of node $i$ and any immediately downstream nodes. In our example, suppose there is no outflux from metabolite $3$. Then $f_{9}=0$ implies $F_{3}=\beta_{1,3}F_{1}$ To write this in our preferred parameters we have $\beta_{1,3}=\frac{k_{3}}{k_{1}}\frac{x_{3}^{\mathcal{T}}}{x_{1}^{\mathcal{T}}}$ or $k_{3}=\beta_{1,3}k_{1}\frac{x_{1}^{\mathcal{T}}}{x_{3}^{\mathcal{T}}}$. In general, any metabolite that does not have an outflow will create a constraint on the parameters which requires reintroduction of the total concentrations of the metabolites. In designing KFP experiments, this leads to a decision that must be made. In some cases, the experimenter may already plan to conduct additional experiments to measure the concentrations of the relevant metabolites, either as part of their larger study to understand mechanistic changes or because they wish to convert the relative fluxes and turnover rates in equation eq. 9 back to the original fluxes in eq. 4. In this case, the constraint induced by the missing outflux can be used to further reduce the number of parameters. If the total concentrations will not be measured, we will need to include a flux out of the pathway for every metabolite. The additional parameter not only increases computation time, but it may also increase the number of time points needed to accurately and confidently estimate relative fluxes and turnover rates.

In this section, we have used the structure of the directed weighted graph representation of a metabolic pathway to formulate differential equations for the unlabeled metabolite concentrations. Scaling this model in terms of isotope enrichment allows us to write the model entirely in terms of relative fluxes and turnover rates. Using the algebraic constraints that the flux in and out of each metabolite must be equal, we can further reduce the number of parameters. In the case where every metabolite has a flux out of the pathway ($|\mathcal{E}_{\mathcal{V}}|=N$) the model can be written in terms of the minimal number of parameters using only relative fluxes and turnover rates. In the case where not all metabolites have a flux out of the pathway, the resulting minimal set of parameters will include dependence on the total concentrations.

Note that in the special case where each metabolite receives at least one input from $\mathcal{E}_{\mathcal{L}}$ or $\mathcal{E}_{\mathcal{U}}$ we can show that $(\boldsymbol{I}-\boldsymbol{B})$ has strictly positive row sums. This is equivalent to the condition that $\alpha_{i}>0$ for all metabolites with no labeled input. Consider a metabolite $i$ with no labeled input. Recall from eq. 10 that the entries of row $i$ of matrix $\boldsymbol{B}$ must add to $1-\alpha_{i}$. This implies that the row sum for the $i^{th}$ row of $(\boldsymbol{I}-\boldsymbol{B})$ is $\alpha_{i}$. In the case of a row corresponding to a metabolite which does have labeled input, the row sum of $\boldsymbol{B}$ must be strictly less than one since we require a non-zero labeled influx. Therefore, $(\boldsymbol{I}-\boldsymbol{B})$ is a Z-Matrix with positive row sums and is invertible [4].

We have chosen the scaled parameters so that the equilibrium equation eq. 12 will have no dependence on the turnover rates in $\boldsymbol{K}$, therefore these turnover rates can not be estimated without some data collected during the transient time. While we have defined the function from $\alpha$ and $\beta$ parameters to the steady state values, this function is not, in general, invertible. Recall from section 2.3 that the number of nonzero parameters in $\vec{\boldsymbol{\alpha}}$ is $|\mathcal{E}_{\mathcal{U}}|$ and that the number of non-zero parameters in $\boldsymbol{B}$ is $|\mathcal{E}_{\mathcal{W}}|$. After substitution using eq. 11, there will be $|\mathcal{E}_{\mathcal{L}}|+|\mathcal{E}_{\mathcal{U}}|+|\mathcal{E}_{\mathcal{W}}|-N$ proportional flux parameters in eq. 12. The simple condition that there must be at least as many steady state values as parameters leads to the following necessary condition for recovery of the $\alpha$ and $\beta$ parameters

Note that the number of arrows leaving the pathway, $\mathcal{E}_{\mathcal{V}}$ is not included in this condition. As we saw in section 2.3, if all nodes have an output vector in $\mathcal{E}_{\mathcal{V}}$ then these output vectors are determined by the algebraic constraints. If one of the nodes does not have such an output vector, the algebraic constraint gives a new condition on parameters. However, this can not reduce the number of parameters in $\vec{\boldsymbol{\alpha}}$ or $\boldsymbol{B}$ so can be ignored for this discussion.

A few simple examples illustrate why this condition eq. 13 is not sufficient. If $\vec{\boldsymbol{\alpha}}=\vec{\boldsymbol{0}}$ then the only steady state is all zeros. Intuitively, this case has no unlabeled input so at steady state, all metabolites are fully labeled. Clearly no information about the values of the $\beta$ parameters can be deduced from this steady state. Another example is if metabolite $i$ only receives input from metabolite $j$, that is $\alpha_{i}=0$, $\beta_{i,k}=0$ for $k\neq j$. In this case, the metabolites $i$ and $j$ will approach the same steady state and the steady states will provide fewer than $N$ equations for specifying underlying flux ratios.

In the case where the graph is an arborescence, the condition is met and the $\alpha$ proportional flux parameters can always be recovered from steady state values. In this case the steady state of the root node will be its proportion of unlabeled input. Every other node $i$, can be written in terms of only its unique upstream metabolite $j$. $\bar{x}_{i}^{ss}=1-(1-\alpha_{i})(1-\bar{x}^{ss}_{j}$). Therefore $\alpha_{i}=\frac{\bar{x}_{i}^{ss}-\bar{x}_{j}^{ss}}{1-\bar{x}_{j}^{ss}}$.

Next we will look in detail at two simple examples which illustrate this necessary condition. In the first, there are two metabolites and the reaction between them is not reversible. In this case, the condition eq. 13 is met and the parameters can be written as functions of the steady state values. Using Bayesian parameter estimation of simulated noisy data, we see that the $\alpha$ parameters are recovered accurately and confidently. Next, we consider the similar case with a reversible reaction. In this case the eq. 13 condition is not met and the estimation of parameters is problematic.

Given data collected only at the steady state of the system, proportional flux parameters can only be uniquely recovered from steady-state estimates if the map from proportional flux parameter values to steady-state values is one-to-one. In this section we explore two examples, one irreversible and one reversible which demonstrate this condition.

The scaled equation for fig. 3 is

which has a steady state

Suppose we can get estimates for the steady-state values of the metabolites, $(\bar{x}_{1}^{ss},\bar{x}_{2}^{ss})$. The proportional parameters $\alpha_{1}$ and $\alpha_{2}$ can be recovered with the map

To illustrate the ability to recover $\alpha_{1}$ and $\alpha_{2}$ in this example, we have used Bayesian parameter estimation with noisy simulated data as described in Appendix appendix A. In all of our Bayesian parameter estimation examples we have used the following general procedure. We generate solutions to the ODE system using known parameter values. Then at each time point, we produce three noisy data points to simulate technical replicates of the experiment. To explore the role of the noise, we have added normally distributed noise to the true solution with a standard deviation that is 2.5, 5, or 10% of the true value. To explore the dependence on the number of time points we have included 3, 5, or 10 equally spaced time points in each parameter estimation. All parameters, including proportional flux parameters and turnover rates are given naive priors. In this section we show only the posterior distributions for the proportional flux parameters, since we anticipate poor estimates of turnover rates. The simulated data, the estimates for the turnover rates, and the credible intervals for each parameter estimation can all be found in the supplement.

In the violin plots in fig. 4, we see that the method consistently provides accurate recovery of both $\alpha_{1}$ and $\alpha_{2}$ as the mode of the posterior distributions. The spread of the violin plots or box plots, is a visual indicator of the uncertainty in our estimate. The comparison across number of time points shows, not surprisingly, that more time points generally provides a more accurate and less variable estimate. The results are relatively insensitive to the level of noise in the data. Increasing the noise level only slightly increases the error and the uncertainty in the estimate.

Next, we turn to the second example shown in fig. 5. In this reversible reaction, the number of proportional flux parameters is 3, and therefore it does not satisfy the condition eq. 13. The scaled system for this example is

which has a steady state

Suppose we can get estimates for the steady-state values of the metabolites, $(\bar{x}_{1}^{ss},\bar{x}_{2}^{ss})$. Since these values depend on all three parameters, there is no way to recover unique values for the parameter with only two pieces of information.

Providing Bayesian parameter estimation with simulated data of this system as described in Appendix A shows the inability to recover all parameter values accurately. In the violin plots in fig. 6, we see that we only get moderately accurate estimates for $\alpha_{2}$. As before, as we increase the number of time points, the uncertainty for $\alpha_{2}$ estimates decreases and the results are relatively insensitive to the magnitude of the noise. The estimates for $\alpha_{1}$ and $\beta_{1,2}$ on the other hand are completely unsatisfactory. The estimates of $\alpha_{1}$ are much lower than the true value, while the estimates for $\beta_{1,2}$ are much higher than the true value. The sample distributions for both $\alpha_{1}$ and $\beta_{1,2}$ show that there is a high uncertainty in the estimates. Since the two steady-state values depend on the values of $\alpha_{1}$, $\beta_{1,2}$, and $\alpha_{2}$, several combinations of parameter values give the same steady-state value. While the method may be able to correctly determine combinations of these parameters that give the correct steady-state values for both $\bar{x}_{1}$ and $\bar{x}_{2}$, it cannot discern the correct parameter values.

We further hypothesize that $\alpha_{1}$ and $\beta_{1,2}$ are more difficult to estimate because of the structure of our steady state solutions. The true values for $\alpha_{1}$, $\beta_{1,2}$, and $\alpha_{2}$ are all fairly small, so if we ignore the quadratic terms in eq. 17 we obtain approximate steady states $\left(\frac{\alpha_{1}}{1-\beta_{1,2}},\frac{\alpha_{1}+\alpha_{2}}{1-\beta_{1,2}}\right)$. In this expression, $\alpha_{1}$ and $\beta_{1,2}$ always occur in the same combination and therefore cannot be distinguished.

In this section we have shown that the estimation of parameters may be severely limited if only steady state data is obtained. In this case, the values for the turnover rates can not be estimated as they do not appear in the steady state equation. It may, however, be possible to estimate the proportional flux parameters if the number of measured steady state exceeds the number of parameters. The specific example with irreversible and reversible reactions illustrates this necessary condition. The examples also makes clear the utility of using a Bayesian parameter estimation. A classic least squares parameter estimate will be able to replicate the limited data but will not necessarily give accurate estimates of the underlying fluxes. The Bayesian approach, however, allows the researcher to simultaneously quantify the best estimate and the certainty in that estimate. This allows for iterative improvement of the experimental protocol for a particular pathway and lowers the chances of misinterpretation.

Here we show that if the metabolite with the fast turnover rate is not the metabolite receiving the labeled nutrient input, the solution is entirely on the slow manifold and the parameter for the fast turnover rate cannot be estimated.

As a result of theorem 3, fast dynamics will not be present in the system and the experimenter will not be able to estimate the faster turnover rate $k_{n}$. In this situation, the experimenter may still be able to estimate all the slower turnover rates and the proportional flux paramters.

In the special case where the metabolite with the quickest turnover rate $X_{n}$ is a target of an edge in $\mathcal{E}_{\mathcal{L}}$, the fast subsystem is defined by

Fast dynamics are present in this system but estimating $k_{n}$ relies on the ability of the experimenter to collect time points on the scale of $k_{n}$.

In the following example, we show that data collected from a chain reaction network with a disparity of time scales does not contain fast dynamics if the slow turnover rate is upstream of a more rapid turnover rate. We will begin by formulating the general reversible chain reaction network model and its fast and slow subsystems.

Consider the reversible chain reaction network in fig. 7. The scaled ODE model can be written in form eq. 9 with the following variables and parameters. For $n=1,\ldots,N$,

and for $n=1,\ldots,N-1$,

Using the reversible two-metabolite example, we will illustrate both cases in section 4.2.1. Consider the reversible two-metabolite reaction network in fig. 5 and its scaled ODE model eq. 16.