Based on a metabolic network, described by the stoichiometric matrix N ∈ ℝ m × r , and given measurements of metabolite pool sizes, gathered in a diagonal matrix P ∈ ℝ k × k , for all metabolites with available fractional amount of measured MIDs, denoted by x ∈ ℝ k × 1   whose appropriate combinations are gathered in a matrix S ∈ ℝ k × r , one can write the system of ODEs for the fractional amount of measured MIDs as follows:

where v ∈ ℝ r × 1 is the vector of fluxes we aim to estimate. For an example of how S looks like for the two CBC models, consult the Supplementary Material and implementation available on GitHub, the matrix S for a toy example is shown in Figure 1. This ODEs can be discretized following Forward Euler approximation as follows:

corresponding to a system of linear equations. The approximation error by the Forward Euler method is in the order of Δ t i and hence, sufficiently small time steps should be considered in the analyzed labeling experiment. The flux distribution, v , does not depend on time, since we assume, as in other MFA approaches, that the network is in a metabolic steady state.

Let ε M I D ( t i ) ∈ ℝ k × 1 denote the deviation between measured and predicted change in concentration of isotopomers obtained by multiplication of metabolite pool size P and fractional amount of MIDs x ( t i ) . Further, let Q M I D ( t i ) ∈ ℝ k × k represent a diagonal matrix that contains the variance of the change in isotopomer concentration between two time points t i and t i + Δ t i , i.e. Y ( t i ) = P x ( t i + Δ t i ) − P x ( t i ) . Assuming that P x ( t i + Δ t i ) and P x ( t i ) are independent random variables, both of variance σ 2 , it follows that the variance of Y ( t i ) is 2 ⋅ σ 2 . In this study, we assume that measurement error, σ , does not depend on time (see Materials and Methods); therefore, Q M I D will not differ between time points. In addition to these constraints, we also use measurements of particular exchange fluxes to further constrain the feasible flux distributions. Let ε V ∈ ℝ 2 × 1 denote the deviation between measured and estimated exchange fluxes and Q V ∈ ℝ 2 × 2 be a diagonal matrix that contains the variance for the measured exchange fluxes.

Based on these simplifying assumptions, flux estimation corresponds to the following quadratic programming problem (QP) that can be regarded as a constraint-regression problem in which we minimize the chi-square goodness-of-fit measure for the change in MID concentration over time as well as external fluxes. The minimization is performed under the constraint that the entire system, represented by stoichiometric matrix N , is at metabolic steady state (see Figure 1 for a graphical overview of the workflow for the presented approach):

with v m a x = 10 6   n m o l   g D W − 1 s − 1 , v m i n = − 10 6   n m o l   g D W − 1 s − 1 for reversible reactions, and 0 otherwise.

To compare the findings of our SFCR approach with those from the state-of-the-art tool INCA, using an EMU-based global approach, we employed data on MIDs from three different algae previously used for flux estimation (Treves et al., 2021). Since we have access to the majority of MIDs of metabolites that participate in the CBC cycle, we formulate constraints for the MIDs of seven metabolites measured in Chlamydomonas reinhardtii and Chlorella sorokiniana under low light (LL, 100 µmol photons m^-2 s^-1), and in Chlorella ohadii under LL and extreme illumination levels (EIL, 3000 µmol photons m^-2 s^-1). Like in Treves et al. (Treves et al., 2021), we also assume that: (i) the structure of the CBC in the considered algae is the same, (ii) there are no microcompartmentation effects on the MIDs of ribulose-1,5-bisophosphate (RuBP), (iii) dihydroxyacetonphosphate (DHAP) is the predominant form of triose phosphates (T3P), hence, its labeling follows the DHAP labeling measured, and (iv) the labeling of ribulose-5-phosphate follows the labeling of pentose phosphates (PP) measured.

We compared flux estimates resulting from integration of data about MIDs in two model variants of different complexity. The model denoted as ‘simple CBC model’, includes reactions of the CBC with only a branch towards starch synthesis; further, metabolites in the simple CBC model are not compartmentalized between the chloroplast and cytosol. The ‘simple CBC model’ included ODEs that describe the nonstationary ¹³C labeling of RuBP, F6P, FBP, G1P, G6P, and ADPG by using the MIDs of PP, RuBP, DHAP, FBP, F6P, G1P, G6P, and ADPG (see Supplementary Information 1.1 for a full list of model reactions and ODEs as well as the full names of the metabolites, and Figure 2 for a graphical model representation). The second model denoted as ‘compartmented CBC model’, includes reactions of the CBC, photorespiration, starch and sucrose syntheses, and considers metabolite compartmentation between chloroplast and cytosol. Based on the structure of the ‘compartmented CBC model’, in addition to the ODEs for MIDs of RuBP, F6P, FBP, G1P, G6P, ADPG, we also formulated ODEs for the MIDs of UDPG (see Supplementary Information 1.2 for a full list of model reactions and ODEs and Figure 2 for a graphical model representation).

Since the proposed estimation of fluxes using SFCR, like other MFA approaches, assumes metabolic steady state, we next determined if this assumption held for the metabolites whose MIDs are described in the model of the investigated algae species. To this end, we conducted ANOVA with the metabolite pool sizes over time. Our findings indicated that RuBP, FBP and ADPG were not at steady state over the investigated time course of C. reinhardtii. For C. ohadii under EIL, ADPG was the only metabolite not at steady state, while for C. sorokiniana all metabolites were found to be at steady state (Supplementary Table S1). Since steady state was ensured for the majority of metabolites, we enforced this constraint in flux estimation. In contrast, for C. ohadii under LL, we found RuBP and PP to be the only metabolite pools at steady state over the entire time course of the measurement. However, if we take out the last time point (40s) FBP, ADPG and UDPG were the only metabolites not at steady state (Supplementary Table S1). Therefore, we only used the first four time points (0s, 5s, 10s and 20s) for flux estimation in C. ohadii under LL.

With these assumptions, we cast the problem of flux estimation as a quadratic optimization problem in which we minimized the chi-squared goodness-of-fit statistic ( χ 2 ) for the considered MIDs along four time points (three time points in case of C. ohadii under LL) as well as flux values for starch and sucrose synthesis reactions (see Materials and Methods). This approach was possible since all MIDs included in the ODEs were measured (Supplementary Table S2). To determine robustness of flux estimates we calculate confidence intervals (CI) for all reaction fluxes (see Methods Section 4.4, Supplementary Table S3, S4, S8) and report χ 2 for the best fit and CI.

Simple CBC model At a significance level α = 0.05 , fitting the MIDs for six metabolites and starch synthesis rate using the ‘simple CBC model’ resulted in a system with 149 degrees of freedom for which the expected range of χ 2 , the range where χ 2 is considered statistically acceptable (Young, 2014), is [117.1, 184.7]. For C. reinhardtii and C. sorokiniana the χ 2 was 87.4 and 85.7 with CI [84.7, 108.6] and [78, 96.5], respectively. Since the obtained χ 2 values were smaller than the lower bound of the interval for the expected χ 2 , there may be an issue of over-fitting. Nevertheless, since we wanted to compare flux estimates across algae under matching set of assumptions, we used the corresponding flux estimates in the comparative analysis. The χ 2 for C. ohadii under LL and EIL were 159.7 (with 108 degrees of freedom, since the last time point was not considered) and 520.6 (149 degrees of freedom), respectively, and were not considered statistically acceptable. Since the relative error in measurements under EIL is not significantly lower than in the other algae, there is no experimental basis for omitting the last time point, as there was for C. ohadii under LL. Our previous modeling efforts with INCA did not find acceptable fit for the data of C. ohadii under EIL (Treves et al., 2021). Possible differences in the model structure for C. ohadii under EIL provides one explanation for the discrepancy.

Compartmented model of CBC and related pathways At a significance level α = 0.05 , fitting the MIDs for seven metabolites using the ‘compartmented CBC model’ variant as well as starch and sucrose synthesis rates resulted in a system with 166 degrees of freedom for which in the expected range of the χ 2 is [132.2, 203.6]. Parameters modeling compartmentation were fixed to random numbers sampled from the interval [0.1, 0.9], with 50 repetitions. We reported the parameterization of c that resulted in the best fit; this may not correspond to the global optima, over all possible parameterizations for c , which is anyhow difficult to ensure for such nonlinear optimization problems. Following this procedure, we found parameterization of c for which the χ 2 CI for C. reinhardtii and C. sorokiniana were not above the expected range with [87, 111] and [110, 132.3], respectively. In case of C. ohadii under LL, we fit a system with 118 degrees of freedom; therefore, at the significance level α = 0.05 , the expected range of χ 2 that would allow to accept model fit is [89.8, 150]. However, the CI for the χ 2 of fit to labeling patterns of C. ohadii was [121.1, 163.7], i.e. the upper bound of the CI was above the expected range. However, the majority (81.6%) of bootstrap samples were within the acceptable range (Supplementary Figure S1); and therefore, we used the flux estimates for comparison across algae. The χ 2 for the fit of ¹³C labeling data from C. ohadii under EIL was 598.6 and hence is not statistically acceptable. Hence, like above, it is likely that modification of the model structure may be required to model labeling dynamics in C. ohadii under EIL.

Both model variants allowed us to fit labeling kinetics for C. reinhardtii and C. sorokiniana. Next, we compared the flux estimates to investigate the effects of model structure used in the fitting. We found that CI for the flux estimates of the reactions in the CBC overlap for the majority of reactions (Figure 3, Supplementary Table S3). In addition, mean flux values were highly correlated with r P = 0.995 (Spearman r S = 0.998 ) for C. reinhardtii and r P = 1 (Spearman r S = 1 ) for C. sorokiniana. As a result, we concluded that the complexity of the model had only a slight effect on the estimated flux ranges with the labeling patterns of C. reinhardtii and C. sorokiniana.

Given that the complexity of the model resulted in small effects on the flux estimates and the goodness-of-fit values, here we provide a more detailed comparative analysis using the more complex, compartmented CBC model across the three algal species under LL. C. ohadhii is the fastest growing green algae known to date, with high photoprotection capacity that facilitates growth under extreme light conditions (Ananyev et al., 2017). Hence, identification of differential fluxes in comparison to algae with slow growth rates provides the potential to understand the underlying molecular mechanisms that underpin the capacity for fast growth (Treves et al., 2021). We found that the estimated mean flux of the carboxylation reaction of RuBisCO under LL condition was the largest in C. ohadii ( 6.7   n m o l   g D W − 1 s − 1 ), followed by C. sorokiniana (4.5   n m o l   g D W − 1 s − 1 ) and C. reinhardtii ( 1.9   n m o l   g D W − 1 s − 1 ) (Figures 4A, B). The mean values for RuBisCO carboxylation were in the expected order across algae, as were the upper bounds of the 95% CI. In line with this observation, previous studies found high photosynthetic rates for C. ohadii (Treves et al., 2016).

However, it must be also noted that the CI for C. sorokiniana and C. ohadii overlap (Figures 4A, B) and the difference between the fluxes was not statistically significant (t-test p-value=0.89). The mean estimated flux through RuBisCO oxygenation reaction was zero for all algae. In line with the observation for RuBisCO carboxylation, mean flux through SBPase is highest in C. ohadii followed by C. sorokiniana and C. reinhradtii (Figures 4A, D). However, the increase in flux for C. ohadii in comparison to C. sorokiniana was not significant (t-test p-value=0.08). For reactions PFP/FBA (Figures 4A, C), we found the smallest flux for C. reinhardtii. While the estimated mean flux of PFP/FBA in C. sorokiniana and C. ohadii were similar, flux sampled by the bootstrap procedure indicated that PFP flux could have significantly larger values in C. sorokiniana than C. ohadii, as observed in the CI. Measured starch synthesis rates for the three algae under LL (Treves et al., 2021) showed that starch synthesis rate was increasing with increasing photosynthesis rate and hence, growth. Moreover, the measured sucrose rates indicated a lack in sucrose production for C. reinhardtii, while both C. sorokiniana and C. ohadii showed similar levels of sucrose synthesis. Using these measurements to fit external fluxes in the QP, we would expect reactions PGM and PGI, the first steps in starch synthesis from CBC intermediates, to show the highest flux in C. ohadii, followed by C. sorokiniana and C. reinhadtii (Figures 4A, E). However, we observed the highest flux towards starch synthesis in C. sorokiniana and zero mean flux in C. ohadii. Flux towards sucrose synthesis was found not to be zero only for C. ohadii (Figures 4A, F). While the absence of sucrose synthesis flux is expected form experimental data integrated in the QP for C. reinhardtii, this is not the case for C. sorokiniana. In addition, it must be noted that the estimated flux values for starch and sucrose synthesis rate are two orders of magnitude smaller in comparison to the measured. The work of Treves et al. (2021) indicated that flux through PEP is differential in C. ohadii compared with the other algae and might also explain differences in growth. However, our analysis was not able to predict these differences. The reason for this lies in the missing coverage of MIDs in the TCA cycle as well as MIDs of PEP needed to correctly predict this flux.

In line with the findings on individual reaction level (e.g. RuBisCO carboxylation, PFP/FBA) the sum of flux through the entire CBC was the highest for C. sorokiniana and the smallest for C. reinhardtii (Supplementary Table S4). However, photosynthesis rates measured for the three algae suggest highest flux in C. ohadii (Treves et al., 2021). For a detailed examination of possible reasons, we point the reader to the Discussion.

Considering the ordering of fluxes based on their magnitudes, we found a very high correlation (Spearman, r S = 1 ) between mean flux of C. reinhardtii and C. sorokiniana. Interestingly, the Spearman correlation of the flux distribution for C. ohadii to that of either C. reinhardtii or C. sorokiniana was smaller ( r S = 0.75 , in both cases). The decrease in correlation for the flux distribution of C. ohadii with those of the other algae can be explained by the complementary pattern of carbon partitioning between starch and sucrose syntheses; specifically, in C. ohadii there was only positive flux towards sucrose synthesis, but no flux through starch synthesis, in contrast to the findings for C. reinhardtii and C. sorokiniana (Figures 4E, F). For a detailed examination of possible reasons, we point the reader to the Discussion.

We also determined the turnover times for three modeled metabolites (Figures 4G–I; Supplementary Table S5). Using the measured pool sizes, we found that the turnover of RuBP was the fastest in C. ohadii (40s), followed by C. reinhardtii (66s) and C. sorokiniana (6min). In line with the observed turnover of RuBP, we found that turnover of FBP was the fastest in C. ohadii (15s), followed by C. reinhardtii (22s) and C. sorokiniana (31s). However, while SBP exhibited the fastest turnover in C. ohadii (18s), C. sorokiniana (24s) showed similar turnover as C. ohadii, followed by C. reinhadtii (49s) with the slowest turnover for SBP. Together, these findings indicated that C. ohadii has more efficient photosynthesis in comparison to C. reinhardtii, indicated by higher carboxylation rates in comparison to C. reinhardtii, while metabolites pool sizes are in the same order of magnitude for both algae. In addition, it can be also concluded that C. ohadii has more efficient photosynthesis in comparison to C. sorokiniana, indicated by similar carboxylation rates for C. ohadii in comparison with C. sorokiniana, despite pool sizes in C. ohadii which are one order of magnitude smaller than those observed for C. sorokiniana (Supplementary Table S7). In addition, C. ohadii has a faster turnover of metabolite pool sizes in comparison to the other algae.

To compare the estimates from the SFCR approach to those from global MFA flux estimation, we used the state-of-the-art tool INCA (Young, 2014) with the ‘compartmented CBC model’. In addition to MID entering the SFCR formulation, the underlying model structure in INCA allowed us to use the MIDs for 3PGA, SBP, S7P, R5P, PEP and 2PG in flux estimation. Therefore, we estimated fluxes using INCA in two scenarios using: (1) MIDs fitted in the SFCR formulation and (2) all MIDs that can be used in INCA. We opted to include MIDs for more metabolites with the aim of increasing the precision of the estimates (expected due to the usage of more data). Altogether, the flux estimation included fitting of pool size, starch and sucrose synthesis rates as well as data on the aforementioned metabolite MIDs for both approaches. For a comparison of fitted labelling pattern obtained by SFCR and INCA see Supplementary Figures S3–S5.

In line with observations from the regression approach, using INCA we obtained a statistically acceptable fit for C. ohadii under LL only when considering time points 0s to 20s. For C. ohadii under EIL no acceptable fit was obtained for both sets of integrated MID data. Using the mean flux estimates from INCA as well as mean flux estimated from SFCR (Supplementary Table S4), we calculated Spearman correlation and found a high qualitative agreement between flux estimates in each algae (C. reinhardtii r s = 1 ,   p − v a l u e < 10 − 10 , C. sorokiniana r s = 0.91 ,   p − v a l u e < 10 − 10 , C. ohadii r s = 0.89 ,   p − v a l u e < 10 − 10 , Supplementary Table S9). In contrast to SFCR, INCA was able to match measured starch and sucrose synthesis rates resulting in a two order increase of absolute flux values in comparison to the SFCR approach (Figure 5). In addition, this discrepancy resulted in different conclusions that can be drawn from the two approaches. SFCR indicated that the main difference across algae was the sugar produced, i.e. starch in C. reinhardtii and C. sorokiniana and, in contrast, sucrose in C. ohadii; in comparison, the results from INCA showed that, in line with the fitted experimental evidence, starch is the main sugar produced across all algae, with increased production for C. sorokiniana and highest for C. ohadii (see Supplementary Table S 4 for a list of mean flux as well as CI of all model reactions obtained from SFCR and INCA). While the solution obtained from INCA included the correct values for starch and sucrose synthesis rates, it should be also noted that CI obtained from INCA often span several orders of magnitude. In contrast, SFCR was able to predict considerably smaller CIs (Supplementary Table S4).

Therefore, we next examined possible reasons for the two order of magnitude difference in flux estimates between INCA and SFCR. The main difference between the modeling with INCA and SFCR is the consideration of pool size as variable or constant (corresponding to mean measurements). Hence, we investigated if the difference of flux estimates may be due to variance in metabolite pool sizes used in INCA and measured pool sizes integrated in the SFCR approach. While the estimated total pool size from INCA matched the measured content for ADPG, RuBP and UDPG, estimates of pool sizes for F6P, FBP, G1P and G6P were different from mean measurements by at least one order of magnitude (Supplementary Table S7), underlining that the used metabolic pool size might be the origin of discrepancy in estimated fluxes. Therefore, as a last check, we repeated the estimation of fluxes with SFCR by fixing the pool sizes to those estimated from INCA. This modeling scenario indicated that pool size from INCA cannot explain the differences of two orders of magnitude to estimates from SFCR. In addition, for C. sorokiniana we even obtained no acceptable fit to experimental data based on the χ 2 = 53812 . Since no flux towards starch and sucrose synthesis is predicted for C. sorokiniana (Figure 5B), we also found a decrease in qualitative agreement ( r s = 0.52 ) between the estimates from SFCR with measured and estimated pool size from INCA (Supplementary Table S9). Similarly, also for C. reinhardtii no flux towards starch or sucrose synthesis is predicted (Figure 5A), although, in contrast to C. sorokiniana the χ 2 is acceptable ( χ 2 = 128.8 ). The only flux distribution that remains unchanged considering the ordering of fluxes is the one of C. ohadii ( r s = 1 ) (Figure 5C). Therefore, pool sizes cannot explain the order difference in estimated flux between INCA and the proposed SFCR approach.

As shown in Section 2.1. the SFCR approach makes use of the Forward Euler approximation to discretize the ODEs. Hence, flux predictions might be affected by large approximation errors as a result of large time steps considered in the experimental setup. To investigate the extent the experimental time steps affect flux predictions from SFCR we used the interpolation of the measured MID data to compare previous results with flux estimation from denser time series (see Methods). We found that the flux estimates obtained from the original time series and the time series with · t being 1 sec are highly correlated ( r = 0.97, Supplementary Table S8) and had no effect on the order of magnitude difference in the estimated fluxes between SFCR and INCA.

We used two variants of a metabolic model for the CBC differing in their complexity. The first model version, denoted as ‘simple CBC model’, considers only CBC with a branch towards starch synthesis and no compartmentation of metabolite pools. The second model version, denoted as ‘compartmented CBC model’, includes reactions for the CBC and related processes, i.e. photorespiration, starch and sucrose synthesis. The model considers metabolite compartmentation into chloroplast and cytosol. Under the assumptions indicated in the main text, ODEs for the MIDs of RuBP, F6P, FBP, G6P, G1P, and ADPG could be written for the ‘simple CBC model’ and ODEs for the MIDs of RuBP, F6P, FBP, G1P, G6P, UDPG, and ADPG could be written using the ‘compartmented CBC model’. A full list of model reactions and ordinary differential equations can be found in Supplementary Information 1.1 and 1.2 for the ‘simple CBC model’ and ‘compartmented CBC model’, respectively. Combined with MID measurements for metabolites RuBP, DHAP (used as MID for T3P), F6P, FBP, G1P, G6P, UDPG and ADPG we used the ODEs to find reaction flux distributions that allow to fit measured labeling patterns.

To assess the quality of the fit we used: (i) the χ 2 statistic with degrees of freedom corresponding to the difference of the number of measurements and the number of modelled parameters (corresponding to the rank of stoichiometric matrix N ); (ii) the distribution of residuals over the individual measurements, and assessed whether it follows a normal distribution with a mean of zero and standard deviation of one (Supplementary Figure S2). Using the χ 2 statistic, the model fit strongly depends on the underlying measurement error. For MID data the standard deviation is often zero for high labeled isotopologues at early time points preventing the model to fit the data. To overcome this problem, we assumed that measurement error, σ , does not depend on time and considered the maximum standard deviation for each isotopologue observed over time as measurement error. In case an isotopologue M m + j had standard deviation of zero across all measured time points, we considered the standard deviation of M m + j − 1 as measurement error for M m + j in data fitting. The code to obtain the measurement error used in model fit, that was calculated based on the standard deviation of measured MID data, can be found on GitHub (https://github.com/ankueken/MFA_Reg).

Having measurements of absolute pool size concentrations and MIDs at five time points 0, 5, 10, 20 and 40s, t i and Δ t i being {0, 5, 10, 20} and {5, 5, 10, 20}, respectively. Note, that for flux estimation in C. ohadii under LL we did not consider the last time point (see main text). Data on metabolite pool size, MID level and starch/sucrose synthesis rates that entered in the program above can be found in Supplementary Table S2. To avoid describing differences that only arise due to reactions that carry flux below solver tolerance, we consider reactions with flux below 10^-6 n m o l   g D W − 1 s − 1 to carry no flux. We used B-splines as implemented in Python package Scipy to interpolate the MIDs using each individual replicate.

To obtain CI and χ 2 for the model fit, we used non-parametric bootstrap by running the following procedure: we resampled the residuals ε , obtained from the QP above, with replacement to generate a new set of residuals ε * . To generate a new bootstrap data set, ε * is added to the measured MIDs. Then the bootstrap dataset was used as an independent replicate experiment, and the QP was solved to calculate new estimates of model parameters. The procedure was repeated 1000 times and bootstrap model parameters as well as bootstrap χ 2 were stored. To obtain 95% CI we computed the 97.5th and the 2.5th percentile of the obtained bootstrap distribution for model parameters and χ 2 .

The fit between model prediction and data can be accepted if the obtained χ 2 for model fit falls within the interval [ χ α 2 2 ( d f ) ,   χ 1 − α 2 2 ( d f ) ] , where α is the significance level, here we use 0.05, and d f is the degree of freedom, corresponding to the difference of the number of measurements and the number of modelled parameters (i.e. the rank of the stoichiometric matrix N ) (Young, 2014).

For a metabolite present in both the chloroplast [ p ] and cytosol [ c ] , the dynamics of the ¹³C labeling of the total pool are the summation of the labeling dynamics observed in the plastid and those observed in the cytosol.

Hence, in the QP above matrix S = S [ p ] + S [ c ] , with S [ p ] combining MID data for plastidial labeling dynamics and S [ c ] combining MID data for cytosolic labeling dynamics. However, compartment specific MID data were not available. Therefore, we substitute MID data in S that model cytosolic labeling dynamics by M m + j [ c ] = c ⋅ M m + j [ p ] , with c denoting the fraction of M m + j [ c ] M m + j [ p ] . Assuming that the labeling dynamics of the total pool are predominantly driven by the reactions in the chloroplast, we sampled values for c from the interval [0.1, 0.9]. Since we assume parameters c to show no drastic change for bootstrapped data, we limit the deviation of parameters c chosen during bootstrap sampling to be at most ±2.5% from the parameterization of c used in the best fit.

The SFCR approach also allows to treat the pool size parameter P as an unknown variable. This assumption will result in the following change in variance σ 2 entering Q M I D in the QP. Assuming that P and x ( t i ) are independent variables the variance of Y ( t i ) is σ P 2 ⋅ σ x ( t i ) 2 + σ P 2 ( μ x ( t i ) ) 2 + σ x ( t i ) 2 ( μ P ) 2 . Moreover, the QP solved included ε P ∈ ℝ n × 1 denote the deviation between measured and estimated pool size for metabolites n whose MIDs enter S and Q P ∈ ℝ n × n be a diagonal matrix that contains the variance for the measured pool size. Based on these simplifying assumptions, flux estimation with variable pool size corresponds to the following QP in which we minimize chi-square goodness-of-fit measure: