Dynamic 13C Flux Analysis Captures the Reorganization of Adipocyte Glucose Metabolism in Response to Insulin.

Cellular metabolism is dynamic, but quantifying non-steady metabolic fluxes by stable isotope tracers presents unique computational challenges. Here, we developed an efficient 13C-tracer dynamic metabolic flux analysis (13C-DMFA) framework for modeling central carbon fluxes that vary over time. We used B-splines to generalize the flux parameterization system and to improve the stability of the optimization algorithm. As proof of concept, we investigated how 3T3-L1 cultured adipocytes acutely metabolize glucose in response to insulin. Insulin rapidly stimulates glucose uptake, but intracellular pathways responded with differing speeds and magnitudes. Fluxes in lower glycolysis increased faster than those in upper glycolysis. Glycolysis fluxes rose disproportionally larger and faster than the tricarboxylic acid cycle, with lactate a primary glucose end product. The uncovered array of flux dynamics suggests that glucose catabolism is additionally regulated beyond uptake to help shunt glucose into appropriate pathways. This work demonstrates the value of using dynamic intracellular fluxes to understand metabolic function and pathway regulation.

Introduction

Cellular metabolism is dynamic, responding to an ever-changing environment (Kotte et al., 2010, Wegner et al., 2015). The cell's metabolic network is a complex ensemble of fast and slow interacting processes, including feedback mechanisms, but these dynamic behaviors are overlooked in a steady-state step-change experiment (Link et al., 2015) (Figure 1A). Thus we need the ability to quantify metabolic responses over time.

Figure 1

13C-DFMA Performed Using B-Splines to Investigate Acute Insulin Metabolic Response

(A) Trajectories of metabolic response and information missed by steady-state measurements.

(B) The trajectory of metabolite P enrichment is a measure of pathway flux contribution and turnover rate. P is produced from unlabeled S1 and labeled S2. Relative to baseline (black), increasing v2 flux (red) increases the speed and amount of P enrichment. Increasing P abundance (blue) reduces the speed of P enrichment.

(C) 3T3-L1 adipocyte tracer experiment showing simultaneous addition of insulin and uniformly labeled glucose and the subsequent sampling time points.

(D) Examples of metabolomics data collected, isotopologues of glucose 6-phosphate (G6P), and malate. Optimized solutions shown as trend lines. Data are represented as mean ±2×SEM of 3 biological replicates.

(E) The optimization strategy used in 13C-DMFA to solve the inverse problem by iterative simulation of isotopologues by estimating system's parameters.

(F) The shape of B-splines is manipulated by controlling the coordinates of the control points (CP). More complex curves have more CPs and are produced by increasing knots and order.

See also Figure S1.

Although metabolomic data are proxy for metabolic responses (Liu and Locasale, 2017), they are not definitive measures of pathway activities (Lee et al., 2017). For instance, when a metabolite decreases, this could be due to reduced synthesis and/or increased degradation (Yugi et al., 2019). Thus, we need to perform flux analysis to quantify the flow of metabolites through metabolic pathways and subsequently reconstruct how the activity of these pathways is altered. This is now commonly accomplished with stable isotope 13C tracer experiments because they generate metabolite enrichment data that have quantitative flow information embedded (Cordes and Metallo, 2016, O'Brien et al., 2015). Fluxes that best explain the observed enrichment patterns can be deciphered by modeling, recovering turnover rates and routes of synthesis from isotopically non-stationary metabolite enrichments (Figure 1B).

However, 13C flux analysis is challenging because we need to track substrate atoms, often through an intertwining network of reactions in central carbon metabolism (Buescher et al., 2015). Network-level deconvolution by mainstream 13C flux analysis packages is largely limited to steady-state metabolism, where fluxes and metabolite levels are constant. More recent software can cater for data sampled from isotopic non-steady state conditions (i.e., metabolism remains constant but with label enrichment increasing over time), but these systems are still studied in a metabolic steady-state (Kogadeeva and Zamboni, 2016, Weitzel et al., 2013, Young, 2014). Techniques like Kinetic Flux Profiling can be used to estimate non-steady fluxes (Alves et al., 2015, Horl et al., 2013, Yuan et al., 2008), but their scope is limited to metabolic pathways with simple label mixing. The current state of the art to infer network-level non-steady fluxes from 13C enrichment data is to use linear piecewise affine approximations of fluxes (PWA) (Abate et al., 2012, Schumacher and Wahl, 2015).

To model curves, PWA's workflow first breaks time-dependent fluxes into linear pieces; then the fluxes are estimated in segments. The first task of optimizing the number and the placement of these breakpoints uses only metabolite abundances, but not enrichments (Schumacher and Wahl, 2015), as currently there is no elegant way to use both. This is not ideal because, in balanced flux scenarios, the metabolite levels will appear constant despite the underlying fluxes being dynamic (Figures S1A–S1D). When optimized breakpoints are too sparse to capture the underlying curvature, the latter flux estimation will struggle to produce a good fit.

Fluxes are estimated by least-square optimization (Wiechert et al., 2001), but this task is computationally intensive when dealing with a regular-sized central carbon metabolism model. The tracking of tracer atom propagations involves solving isotope balance equations. These ordinary differential equations (ODE) form a stiff problem owing to the coexistence of very slow and very fast reactions. Implementing non-steady fluxes worsens the situation because metabolite pool sizes are changing at different paces as well. Any task that involves running optimization repeatedly becomes overwhelmingly slow, such as model and data troubleshooting and sensitivity analyses by resampling. Thus, we aimed to make non-steady-state metabolic flux analysis of stable isotope tracer data more generalized and implementable.

Here, we present a generalized computational framework for performing 13C dynamic metabolic flux analysis (13C-DMFA). The crux of our solution was using B-splines to model time-dependent fluxes (Martínez et al., 2015), analogous to drawing curves in vector-based graphics. B-splines made it possible to (1) use higher-order curves to reconstruct a variety of flux trajectories and (2) apply a more stable algorithm to speed up flux estimation.

As a proof of principle, 13C-DMFA was applied to characterize how cultured adipocytes process glucose in response to insulin, which stimulates protein phosphorylation within seconds and glucose uptake within minutes (Humphrey et al., 2013, Krycer et al., 2017). We found that flux rates varied substantially between pathways, responding to insulin with differing speeds and to differing magnitudes. The modular organization-portrayed pathway flux dynamics suggests that glucose catabolism was not simply an autonomous cascade driven by glucose influx. Overall, this work provides a MATLAB-based package and a workflow to perform 13C-tracer flux analysis on non-steady-state systems.

Results

Network-Level Modeling Is Required to Extract Dynamic Flux Information

We previously examined how the adipocyte metabolome changed in response to acute insulin treatment (Krycer et al., 2017). We used 3T3-L1 cells (Green and Kehinde, 1975), a mouse embryonic fibroblast cell line that was differentiated into adipocytes before treatment. Briefly, 3T3-L1 adipocytes were exchanged to media containing uniformly labeled 13C6-glucose and immediately treated with or without insulin (INS and CTRL, respectively) (Figure 1C). Lactate in the media and 27 intracellular metabolites were then measured at six time points (1, 5, 10, 20, 40, and 60 min) after treatment by mass spectrometry (capillary electrophoresis-time of flight and electrospray ionization-tandem mass spectrometry) and then normalized to protein content (units: pmol/mg protein). It was obvious that insulin altered pathway fluxes, if not, metabolite abundances (Figure 1D). Insulin increased glucose uptake rapidly, and glucose propagated into various metabolic pathways at different speeds (Krycer et al., 2017). This was assessed by examining the levels and 13C-enrichment of individual metabolites, treating them as separate entities.

Here, we sought to apply network-level flux modeling to understand how adipocytes process glucose in response to insulin. Enrichment states of the metabolites were increasing in a cascading manner, with distal metabolites like tricarboxylic acid (TCA) cycle intermediates enriching slower than metabolites proximal to glucose uptake. The lag observed in the metabolite enrichments is intricately linked to how metabolic fluxes are changing with insulin stimulation and how the metabolites are connected. Thus we applied mass conservation principles, namely, via metabolite balance equations and isotopologue balance equations (Antoniewicz et al., 2007, Quek et al., 2010), to infer these fluxes by modeling the transfer of carbon atoms between metabolites. As the steady-state assumption is no longer valid, this raised the challenge of estimating time-dependent fluxes by optimization, which we address in the following sections.

Non-steady Fluxes and Near-Zero Isotopologues Make the 13C-DMFA Problem Hard to Optimize

Generalizing for a network of reactions, the isotopologues of metabolites over time X(t) can be expressed as a function of the metabolic fluxes v(t) and the system's boundary conditions, namely, the initial metabolite concentrations c(t=0) and the enrichment of input substrates Xsubstrate. X(t) is often expressed in the intensive form x(t)=X(t)/c(t) as isotopologue fractions; note that and . This forward problem is a system of non-linear ODEs (Dai and Locasale, 2017), which means the inverse problem of estimating v(t) from x(t) is solved by least-square optimization (Figure 1E). Mainstream techniques are limited to constant v(t) and x(t), i.e., steady-state metabolism (Buescher et al., 2015, Wiechert and de Graaf, 1997), due to the computational difficulties faced when optimizing non-steady fluxes.

The crux of dynamic flux analysis is to parameterize each v(t) independently, but this creates computational issues during the optimization. When a metabolite's total influx and efflux are not equal (i.e., imbalanced), the metabolite's concentration and isotopologue abundance will increase (influx > efflux) or decrease (influx < efflux). The allowance for imbalanced fluxes significantly increases the encounter of infeasible points where c(t)<0 and x(t)<0 during optimization, especially when effluxes can exceed influxes. Core to this problem is that numerically ensuring c(t)>0 does not guarantee x(t)>0, as the latter is estimated numerically with smaller steps.

The risk of encountering negative is ever present despite having established c(t)>0 numerically, because (1) most x(t) are trending toward or away from zero in a pulse tracer experiment, (2) the initial guesses for optimization are almost always infeasible, and (3) underdetermined parameters often settle just within the feasible envelope. Negative c(t) and x(t) are not caused by numerical or truncation errors; in an iterative method they are still mathematically valid outputs of v(t) because efflux exceeding influx can lead to a deficit. However, as x(t) is an intensive variable, adaptive ODE solvers cannot progress past the point of reaching negative x(t) due to subsequent numerical instability, causing the solver to take even smaller steps.

Currently there is no way to avoid breaching the feasible envelope when ODE solvers are used during optimization. In particular, infeasible points are difficult to anticipate numerically because negative values may only emerge when smaller steps are taken (Figure 2A). The abrupt stalling of ODE solvers causes optimizations to terminate prematurely, which occurs more frequently with larger models. Thus, we considered B-splines as an alternative, to efficiently (1) contain the optimization in the feasible space and (2) guarantee x(t) is positive when evaluated in the feasible space.

Figure 2

The SBR Algorithm Overcomes Stability Issue Encountered by Conventional ODE Solvers

(A) As SBR utilizes predetermined time steps, negative isotopologues are avoided when x(t) is evaluated at the given time steps, even though negative values are part of the true solution. ODE solvers are unable to forecast the infeasibility because time steps are determined dynamically.

(B) Isotopologue propagation in a metabolic network is abstracted as a series of tanks and modeled by mass conservation equations.

(C) Schematic of the Sequencing Batch Reactors (SBR) algorithm. Each tank follows the cyclic sequence of fill, mix, and drain. The drain step is not computed by the SBR algorithm.

(D) The schedule of the sequence shows that the “drain” step can be omitted because x(t+1) is an intensive value calculated during the “mix” step.

(E) The toy TCA cycle model used to test performance of the SBR algorithm.

(F) SBR is generally faster and converges more consistently at the theoretical solution than MATLAB's ODE45 and ODE15s solvers. Manhattan distances are used to measure deviations from the true solution, with bars showing median values. Error bars show interquartile range of 100 replicates.

(G) Benchmarking SBR's speed-accuracy trade-off at the scale of a real (i.e., adipocyte) model. EMUs generated by SBR are as accurate as ODE solvers', but will require as much computational time as ODE15s without the Jacobian, which is slower than with the Jacobian. SSRs measure the weighted deviations of EMUs from the “true” measurements generated by ODE15s (with Jacobian) using prototyping CTRL and INS solutions. Bars show median values. SBRx2.5 and SBRx20 refer to increasing SBR's number of steps by 2.5- and 20-fold relative to the optimized SBR time steps showed in Figure 3D (SBRx1). Error bars show interquartile range of 100 replicates.

(H) The expected error in the CTRL and INS datasets when corrupted randomly with one standard deviation noise.

(I) When optimized again using ODE15, five randomly chosen CTRL and INS solutions generated by SBR did not change significantly relative to the rest of the SBR solutions.

Figure 3

A 13C-DMFA Workflow that Incorporates a Prototyping Phase to Improve Model Fit

(A) Overview of the 13C-DMFA workflow, showing iterative input revision and prototyping, followed by optimization and results interrogation.

(B) An abridged diagram of adipocytes central metabolism showing placement of key inputs and outputs.

(C) AIC was used to help choose the B-spline orders and number of knots used for parameter estimation chosen based on the lowest AIC for INS. Placement of the first (blue), second (red), and third (orange) knots shown, with correspondingly colored cubic curves marking knots' lower envelope.

(D) SBR time step configuration of [500, 400, 200, 75, 75, 50] steps between sampling points manually chosen at the elbow point as a good compromise between numerical accuracy and computation speed. As reference for SSR calculations, a synthetic “error-free” dataset was generated by SBR using uniform 30,000-step and the best prototyping INS solution.

(E) Optimized SBR time step (red cross) still showed a good compromise between accuracy and speed when retested against actual INS solution and synthetic data generated by ODE15s (with Jacobian). Numerical error was less than theoretical error estimated from measurement noise.

See also Figures S2 and S3.

B-Splines Provide Means to Keep Optimization in the Feasible Space

Briefly, B-splines are piecewise polynomial functions that can be expressed as basis functions (Figure 1F). They can be used to draw any curve, directly addressing our need to parameterize time-dependent flux trajectories in a generalized manner. B-spline computation is fast because the system can be rendered linear, enabling the use of matrix algebra.

Here, time-dependent trajectories of fluxes v(t) and metabolites c(t) are expressed as a linear combination of B-spline basis functions N and ∫N, respectively (Equations 1 and 2). N and ∫N are coefficient matrices generated by the Cox-de Boor recursion formulation, using a prescribed set of time steps, knot positions, and order of splines (e.g., 1 = constant; 2 = linear; 3 = quadratic) (Figures S1E–S1H) (Martínez et al., 2015, Vercammen et al., 2014). Using the stoichiometric matrix S, a matrix representation of the metabolite balance equations (Quek et al., 2010), c(t) is calculated analytically by integrating v(t), given the initial metabolite concentrations c0. The B-spline construction and integration are further elaborated in Transparent Methods.for R reactions, M metabolites, T time steps, and j = (order + number of internal knots).

Equations 1 and 2 show that we can control v(t) using just one very simple system of parameterization: the number of control points (CP) and their XY-coordinates and order of the B-spline (Figure S1G). This system resembles how we construct B-splines in graphic software. Note the number of knots is the number of CP minus the B-spline order; c(t) is determined from v(t) and c0 (Equation 2). With reference to PWA, knots are equivalent to breakpoints, the order is 2 (i.e., linear), and the formulation for v(t) and c(t) are exactly the same. The main difference is that to model flux trajectories with greater curvature, the order of B-splines can be increased as well instead of just the number of knots/breakpoints.

Most importantly, once the knot sequence and order of splines are prescribed, v(t) and c(t) are linear with respect to the estimated parameters CP. This circumvents solving c(t) numerically, thus allowing the incorporation of c(t)>0 as linear constraints to define the feasible bounds for optimization.

Euler Method Improvised to Eliminate Negative Isotopologues

Next, we sought to eliminate negative isotopologues x(t)<0. To achieve this, first we will describe how X(t) is calculated. Cell metabolism can be abstracted as a network of continuous stirred tank reactors, with the contents of a tank and the flow between tanks representing the metabolite pool and their inter-conversion, respectively (Figure 2B) (Cordes and Metallo, 2016). This abstraction can be dissected further at the sub-molecular level, treating atoms or group of atoms of metabolites as separate tanks, for instance, to deal with cleavage reactions.

This is accomplished by the Elementary Metabolite Unit (EMU) framework (Antoniewicz et al., 2007). Within this framework, X(t) denotes isotopologue abundances at time t. An EMU of a metabolite is a distinct group of the metabolite's atoms. For example, three-carbon pyruvate can be decomposed into seven EMUs sized from 1 to 3 (PYR1, PYR2, PYR3, PYR12, PYR23, PYR13 and PYR123, where PYR denotes carbon X as a member of that particular EMU). The propagation of 13C atoms through a network of isotopologues can be tracked by mass conservation. Equation 3 represents a generalized form of the EMU balance equation, with X and x now representing EMUs and their fractions. The EMU framework is elegant because it invokes the least number of balancing equations and is compatible with mass spectrometry data.where and

Equation 3 is often solved numerically using ODE solvers. The forward Euler Method (Equation 4) is the simplest numerical approach but is sufficient to demonstrate that X can become negative when the EMU efflux term is too large. X<0 is difficult to avoid numerically because adaptive solvers reduce the step size dynamically to suppress numerical error, but the sudden occurrence of negative X(t) would destabilize the equation system and amplify numerical error as a result. This causes the ODE solver to take even smaller steps, until the solver stalls completely because negative X(t) is indeed the solution. In our hands, rounding negative values up to zero did not improve convergence rate of optimization. In contrast, B-splines uses predetermined time steps and are not adjusted dynamically.

To eliminate the negative term altogether, we simply treated metabolite conversions as “sequencing batch reactors” (SBR) (Figure 2C). The operation of SBR is composed of three stages, namely fill, mix, and drain, cycled repeatedly in discrete steps to model the continuous conversion of metabolites over time (Figure 2D). In each step, the SBR algorithm cycles through the three stages, first accounting for the bolus addition of new metabolites to the existing pool (fill), followed by mixing of the new and existing metabolites (mix), and finally reducing the post-mix pool to produce the existing pool for the next cycle (drain). The “mixing” equation sums the normalized EMU contributions of source substrates and the existing pool at a given time step t (Equation 5). Note that x is an intensive variable carried over to the next cycle and is calculated before the “drain” step.

The discretization in SBR takes advantage of how B-splines are calculated using predetermined time steps. With the efflux term absent in Equation 5, x is always positive when Equation 5 is evaluated at time steps satisfying the constraint c>0. Thus we circumvented negative values during optimization as long as c are positive at the time steps taken (Equation 2), even if solutions for c(t) or x(t) include negative values. Despite eliminating the efflux term, mass conservation is still enforced because the size of the existing pool c is calculated analytically with respect to v(t) by B-spline integration.

The SBR Algorithm Is Fast, Accurate, and Numerically Stable

The SBR approach was benchmarked against ODE optimization using an eight-reaction toy model representing the citric acid cycle (TCA) (Figure 2E, see Transparent Methods). Optimization solutions were generated by solving Equation 3 using MATLAB's ODE45 and ODE15s, and by the SBR approach using Equation 5. All methods converged at the theoretical solution. The Manhattan distance between the optimized values and the theoretical solution showed that SBR excelled by consistently converging at the optimum, whereas ODE solvers suffered failure rates of at least 60% (Figure 2F). The inability of the ODE solvers to converge consistently at this toy scale, indicated by a large spread in the calculated Manhattan distances when compared with the SBR approach, demonstrated potential optimization challenges for real models.

SBR time steps with decreasing sizes were tested to evaluate the trade-off between the numerical accuracy of calculating x(t) and runtime. For an SBR approach to be as accurate as ODE45, 35 steps were required, but computation was still 10-fold faster based on the median runtime of optimizations that converged (Figure 2F).

The above-mentioned results indicated that we have a good margin for calibrating SBR time steps such that high/comparable computational accuracy is obtained without being slower than ODE solvers. We next made a similar comparison, using INS and CTRL datasets. We found that to achieve the same level of accuracy as ODE solvers, the SBR algorithm took as long as the ODE15s without the Jacobian (Figure 2G). Matching the time taken by ODE15s with Jacobian, which was as accurate but 7- to 10-fold faster than ODE15s without the Jacobian, increased SBR errors (Figure 2G), but they were still within the tolerance of noisy data (error <500, Figure 2H). As expected, the numerical errors incurred for CTRL were greater than for INS because SBR time steps were calibrated using INS dataset. We struggled to complete the optimizations when we used ODE solvers, but we showed that re-optimizing SBR solutions using ODE15s with Jacobian neither shifted the SBR solutions significantly nor improved the fit (Figure 2I). Thus the performance of the SBR algorithm is at least on par with ODE solvers, but is significantly more stable.

Using the Adipocyte Data to Design a Practical 13C-DMFA Workflow

Another advantage of the SBR approach is that it makes our workflow highly amenable to prototyping, the iterative process of testing model-to-data fit followed by ad hoc modifications (Figure 3A). This takes advantage of the ease to manipulate optimizations for speed over accuracy to achieve quicker learning cycles to improve the model. For prototyping, the SBR time steps taken were coarse and arbitrary, 20 uniform steps between sampling time points.

The optimization objection function (Equation 6) was defined based on isotopologue concentrations for 27 intracellular metabolites and extracellular lactate over time and the incorporation of radiolabeled 14C-glucose into glycogen and fatty acids at 60 min (Data S1). The construction and components of Equation 6 are described in greater detail in Figure S2A and Transparent Methods. The optimization procedure implemented largely adhered to existing 13C-MFA methodologies (Quek et al., 2009) (see Transparent Methods, Data S4). Briefly, CP and c0 (from Equation 2) were estimated by minimizing the sum of squared residuals (SSR) calculated by the objective function (Equation 6), using MATLAB's constrained non-linear solver fmincon. Residual errors between the measured (data) and simulated (data) were weighted by the respective measurement errors (error weight). The objective function was scripted dynamically to accommodate frequent changes during prototyping, as well as to deal with missing values and gaps in the metabolite data (Figure S2B).subjected to c(t)>0, v(t)≥ 0

We used the adipocyte metabolite data to dictate the reconstruction of the central carbon model (Figures 3B and 4A, see Transparent Methods). Reaction pathways were assembled such that metabolites showing significant enrichment were traceable to glucose. Redundant or alternative routes were included at first, and subsequently marginalized by examining the fit of the candidate model against the adipocyte data (results not shown). For example, lactate had m+3 and m+2 label, but there was a lack of substantial m+3 and m+2 labeling for glucose 6-phosphate (G6P) and phosphoenolpyruvate (PEP) respectively, suggesting that gluconeogenic (or glyceroneogenic) flux was not significant relative to glycolysis; thus we eliminated PEP carboxykinase. In addition, the high m+3 label in malate indicated that pyruvate anapleurosis should be included. A single compartment model was adopted because, for adipocytes, metabolites are largely synthesized or consumed in a single compartment. There is a general lack of parallel pathways that can give rise to compartment-specific enrichment profiles. For example, as glucose and branched-chain amino acids are the major substrates for de novo lipogenesis (Green et al., 2016), and reductive carboxylation is negligible in adipocytes (Koh et al., 2004, Liu et al., 2016), the cytosolic synthesis of citrate from oxoglutarate was omitted. Furthermore, due to the lack of compartment-specific metabolite data, it was more conservative to use a single-compartment model to avoid overfitting. The availability of compartment-specific metabolite data in the future would enable the modeling of additional pathways such as the malate-aspartate shuttle, which likely operates in our cell system but could not be resolved by our tracer dataset.

Figure 4

Temporal Fluxes Showed Differences in Dynamics and Magnitude

(A) Network map showing median temporal net fluxes and metabolite abundances.

(B) Isotopologue data and fitted solutions generated by Monte-Carlo (left column). Experimental data are represented as mean ±2×SEM of 3 biological replicates. Temporal fluxes generated to recapitulate the isotopologue data showed differences in flux dynamics (middle column). Line shows median flux and shaded area shows interquartile range from 50 Monte-Carlo replicates. The distributions of average fluxes show summarized differences in flux magnitude (right column). Average fluxes were calculated by integrating the area under the curve from 10 to 60 min.

See also Figure S4.

We also considered metabolic inputs and outputs to the reaction network, which are important to consolidate the metabolite data (Figure 3B). The former causes dilution of label at specific points of the network, and the latter draws carbon out of the system to create flow through a pathway. Major inputs (glycolysis, glutamine anapleurosis) and outputs (lactate production, aspartate catapleurosis) were specified upfront, then minor inputs and outputs were deduced from the data, either by inspection of the data or by modeling trial and error. For example, pyruvate and AcCoA influxes were found to be essential to create entry of unlabeled substrates upstream of the TCA cycle in the CTRL dataset (Data S1). Unlabeled pyruvate can be produced by alanine transaminase, and AcCoA, by fatty acid oxidation in the mitochondria.

The final adipocyte metabolic model used for flux estimation contained 66 reactions and 34 internal metabolites; the same model was used for both INS and CTRL datasets to facilitate comparison (Data S1). Initially, the model could not recapitulate the unlabeled (m0) fractions, which were unexpectedly high and stable (Figure S3A). Thus, nuisance parameters were introduced ad hoc to represent “stagnant” partitions (Antoniewicz, 2018), which improved data fitting significantly as shown by a reduction in Akaike information criterion (AIC) (Figure S3B).

B-splines can model any flux trajectory, but careful parameterization of order and knot is needed to avoid overfitting. Thus, in parallel to model reconstruction, we tested configurations of B-splines up to orders of three, in combination with up to three internal knots (Figure 3C, see Transparent Methods). Using AIC to balance overfitting and underfitting, we found that the final adipocyte model was best modeled using quadratic fluxes (order = 3) and one internal knot. This meant four CP values were used to model each v(t). The flux model had a total of 329 estimated parameters: 66 fluxes × 4, 37 initial concentrations, and 28 stagnant pools. In comparison, INS and CTRL dataset each has a total of 750 and 671 data points, respectively; the difference is due to missing values. Underfitting was not an immediate issue as the calculated SSRs were larger than the chi-square limit.

The final prototyping task was to calibrate the SBR time steps (see Transparent Methods). The aim of this calibration task is to reduce numerical error incurred by our SBR algorithm; commercial ODE solvers implement adaptive step size to achieve this. Smaller time steps are needed to more accurately model faster rates of enrichment, thus the calibration was performed on INS rather than CTRL. Briefly, based on the best prototyping solution of INS, we substituted data with c(t) and x(t) to calculate data. We started with a uniform 30,000 step as upper limit and tested various SBR time steps generated by a combination of grid search and random sampling, with a uniform 20-step set as the lower limit. Employing a strategy analogous to the Elbow method for clustering optima, SSRs here were used to show that [500,400,200,75,75,50] steps between sampling time points achieved a good balance between speed and accuracy (Figure 3D). Indeed, this configuration reflected that metabolite enrichments in the INS dataset were rapidly increasing within the first four time points (1–20 min). When final flux results were generated, we verified that the error incurred by this optimized set of SBR time steps was small relative to the error associated to measurement noise (Figures 2H and 3E). Compared with the uniform 20-step used for prototyping, this time step configuration was 10-fold slower; a single optimization generally took 30 h to complete (CPU 3.4 GHz).

To generate flux results in parameter estimation (Figure 3A), a multi-start procedure—repeating the optimization using random initial values—was implemented to address the local optimum problem of non-linear optimization (see Transparent Methods). In parallel, the multi-start procedure also sampled for different knot positions, because this parameter was set before the optimization. Thus, knot optimization incorporated both metabolite abundances and enrichments, unlike PWA. SSRs were sensitive to knot position, with CTRL knots placed early and INS knots more spread out for the lowest 20% SSRs (Figure S3C).

Despite our best efforts, we could not bring the minimum SSR for both INS and CTRL below the chi-square limits. Visually, the estimated fluxes have reproduced the respective dataset (Data S3), but the normal probability plots of the residual errors showed a slight departure from normality (Figure S3E). However, the data points with high residuals did not belong to specific pathway or metabolites (see Transparent Methods). We may have underestimated measurement uncertainties because increasing the minimum error threshold from 1% to 5% can bring SSRs down to within the chi-square limits without significantly changing the optimized solution, as the optimum solutions with 5% minimum error were still within proximity of the original Monte-Carlo solutions (Figure S3I). The application of a blanket 1% minimum error in our modeling strategy mainly reflects the expected instrument accuracy of determining relative enrichment of isotopologues (<1.5%) (Buescher et al., 2015), but this does not capture biological variance. Without any apparent systematic bias as leads to troubleshoot, we proceeded without making any further modification. Overall, this issue raised the caveat that although estimated fluxes were generally correct, they still contained systematic biases. Nonetheless, we have validated that our optimization approach was able to converge at the global optimum when challenged with a synthetic INS and CTRL dataset (Figure S3D).

Last, Monte-Carlo resampling was used to evaluate the spread of the estimated parameters; this included knot positions (Figures S3F and S3H). We generated 50 corrupted sets of metabolite data per condition, CTRL and INS, using one standard deviation of Gaussian noise, and multi-start optimization was again performed on each set. By resampling, we verified that 50 iterations were sufficient because the medians and standard deviations of the estimated parameters stabilized by the first 30 iterations of sampling (Figure S3G). By 50 samples, the standard deviations of 245 of 329 parameters have stabilized within 1% (ΔSD/SD). Subsequent interrogations of flux results and their distributions were based on these 2×50 solutions.

Insulin Stimulates a Heterogeneous Metabolic Response beyond Glucose Uptake

We previously established that insulin treatment rapidly increased glucose influx and directed glucose flow in a pathway-specific manner in cultured adipocytes (Krycer et al., 2017). This was corroborated here by an elevation in metabolic fluxes across central carbon metabolism by insulin, with this effect sustained over the time course (Figure 4A). On the other hand, CTRL-treated cells exhibited initially high fluxes that decayed rapidly without insulin (Figures 4A and 4B), likely a result of media exchange (discussed in the next section).

Examining insulin-responsive changes in flux rates (INS-CTRL), we found that glycolytic flux only took 10 min to reach a new, higher steady state at about 40 nmol/mg/h (Figure 4B, Data S2). The acceleration rate of glycolysis was consistent with published extracellular acidification rates of insulin-stimulated 3T3-L1 cells (Schreiber et al., 2017). If glucose influx was the sole determinant of central carbon fluxes, we would expect that a majority of glucose-dependent (13C-labeled) pathway fluxes rise concomitantly with glycolysis in response to insulin. However, branching pathways reached differing magnitudes (Figures 4B and S4A). For instance, pyruvate dehydrogenase (PDH) flux only increased from 0.21 to 0.72 nmol/mg/h (Figure S4A). These pathways peripheral to glycolysis also showed a more gradual ramping in response to insulin. The flux of pyruvate carboxylase (PC), transketolase (TK), and malic enzyme (ME) took about 30–40 min to reach maximum velocity (Figures 4B and S4A). Overall, these flux differences suggest that glucose metabolism is driven by more than just glucose influx.

Furthermore, these flux differences were not always evident from the metabolite data itself. For instance, metabolite abundances in both oxidative and non-oxidative arms of pentose phosphate pathway rose concurrently in response to insulin (Krycer et al., 2017), but it was TK flux that substantially increased from near zero to 4 nmol/mg/h, with G6P dehydrogenase (G6PDH) flux remaining constant at 0.2 nmol/mg/h (Figures 4B and S4A). Parenthetically, in contrast, ME flux increased 3.7-fold to 41 nmol/mg/h by insulin (Figure S4A), reinforcing previous findings that ME is the primary site of NADPH generation in adipocytes (Liu et al., 2016). Overall, this demonstrates the value of calculating fluxes to quantify pathway activity.

To explore these flux changes on a broader scale, we mapped insulin-dependent changes onto the metabolic network to understand how adipocytes handle imported glucose. This revealed that adjacent steps within pathways responded similarly to insulin but differed markedly from other pathways or even other segments within the same pathway. We reasoned that flux-controlling steps regulated this modularity. Thus, we identified modules by unsupervised K-means clustering, generating eight optimum clusters (Figures 5A, S5A, and S5B, see Transparent Methods). When overlayed onto the pathway map, this analysis revealed that segments of the glycolytic pathway responded differently to insulin (Figure 5B). For instance, the flux rates in lower glycolysis (#7) increased more rapidly than upper glycolysis (#5), the latter showing a potential overshoot. At the border between upper and lower glycolysis, the reversible aldolase and triose phosphate isomerase steps (#4) were very active. The punctuation of flux dynamics at GAPDH reinforced this enzyme as a flux-controlling step of glycolysis (Shestov et al., 2014). Thus, the enzymes downstream of glucose transport were not simply driven by glucose influx. The modular system suggested a feedforward mechanism, whereby lower glycolysis pre-emptively increased in flux before glucose transport and upper glycolysis.

Figure 5

Adipocyte Fluxes Responded Rapidly to Insulin in a Modular Fashion

(A) Eight clusters of insulin responses grouped based on the control points of temporal flux differences (INS minus CTRL). The control points summarize the trajectory of each flux difference into five points in time, with the vertical axis normalized to [-1,1] using the maximum flux and the horizontal axis representing time in minutes. Black lines represent centroid location.

(B) Network reactions were colored using cluster information based on majority. Dashed arrows represent net flux. Metabolites and enzymes in the glycolysis pathway outlined.

See also Figure S5.

Similarly, the TCA cycle was segmented, with pyruvate situated at the interface (Faubert et al., 2017). The series of enzymes converting 2-oxoglutarate to fumarate were clustered together (#2), demonstrating that the TCA cycle still operated in the canonical oxidative fashion (Figure 5B). Superimposed were the large PC and ME fluxes that were increasing together (#5), albeit more gradually, in greater magnitude than the oxidative half of the TCA cycle. The extracellular lactate m+2 fraction was unusually large (Figure S4B) and directly constrained the model to produce high cyclic pyruvate carboxylation-decarboxylation flux, a process that is coupled to pyruvate carbon being scrambled/inverted at symmetrical fumarate to produce the m+2 lactate signal (10% for CTRL, 19% for INS). This futile cyclic flux around PC and ME was insulin responsive and a potential avenue for cell to generate NADPH for lipogenesis from surplus NADH by consuming ATP.

Overall, the formation of clusters based on flux dynamics and magnitude demonstrated that adipocytes' central carbon metabolism is more complex than just a glucose overflow driven by uptake. This suggests that metabolic pathways are organized into modules to better coordinate diversion and utilization of glucose.

The Effect of Media Exchange Depends on Insulin Pre-conditioning

We expected the CTRL cells to remain in metabolic steady state, but we observed non-steady fluxes in these cells (Figure 4B). “Constant flux” B-splines (order = 1, knots = 0) gave high SSRs (≥5,000) and AIC scores (Figure 3C). Indeed, “hockey stick” fluxes, characterized by a high but rapidly declining flux in the first 10 min, were required to achieve fast and early enrichment of CTRL metabolites (Data S3).

We speculated that this was due to media exchange between the pre-incubation and labeling/treatment periods (Figure 6A) (Krycer et al., 2017), involving PBS and DMEM washes to maximize glucose labeling at ∼100%. Based on CTRL fluxes, the effects of these perturbations dissipated within 10 min (Figures 4B and S4A). We explored this further, comparing media exchange to the effects of insulin stimulation. The abundances of central carbon metabolites under CTRL and INS conditions were analyzed using principal-component analysis (PCA). Principal component 1 (PC1) appeared to have captured the insulin-dependent evolution of metabolite levels over time (Figure 6B). Unstimulated cells CTRL-10 and CTRL-60 formed a tight cluster relative to the plotted trajectory of INS-1 to INS-60 and were proximal to INS-1 and INS-5. This indicated that metabolite abundances in the CTRL dataset appeared to have normalized quickly after media exchange, much like the transient initial fluxes (Figure 4B). Thus the effects of media exchange were neither significant nor prolonged.

Figure 6

The Effects of Media Exchange Exacerbated by Insulin Preconditioning

(A) Media exchange experimental design, showing the timings of media exchange and harvest, contrasted against the main 13C-DMFA experiment. Conditioned adipocytes denoted by the prime symbol (′).

(B) PCA plot of total abundances of central carbon metabolites (in triplicates). The ellipses demarcate the two groups of adipocytes, naive and conditioned. INS-10 did not deviate from the insulin response trajectory captured by PC1, whereas metabolite abundances of INS-10′ and INS-60′ were markedly different.

(C) Heatmap compares fold changes of metabolite levels between naive and conditioned adipocytes, showing build-up of metabolites above GAPDH and mononucleotides in conditioned adipocytes. 2,3-DPG: 2,3-bisphosphoglycerate.

(D) Energy charge ratios summarize changes of nucleotide levels, showing a significant drop in energy charge when media exchange was performed on conditioned adipocyte (INS-10′). Data are represented as mean ± SD of 3 biological replicates.

(E) Relative enrichment of glycolytic metabolites at 10 min. INS-10′ enrichments suggest that flux constriction at GAPDH was partial and glucose was still catabolized to pyruvate.

We tested whether media exchange perturbed insulin-conditioned adipocytes more than untreated (naive) cells. Here, adipocytes were first exposed to insulin for 50 min, and then the culture media was exchanged from 12C- to 13C-glucose-containing media 10 min before harvesting (INS-10′). Effectively, except for INS-10′, all media exchanges performed were on adipocytes in the naive state (Figure 6A). The distance between INS-60′ and INS-10′ indicated that the metabolic perturbations by media exchange were significant in adipocytes with an established insulin response (Figure 6B). The direction of shift along PC2 indicated that the perturbations were very different from insulin responses observed in naive adipocytes. These differences, and the proximity of INS-1, INS-5, and INS-10 to the cluster of CTRL data points, provide evidence that the flux differences between INS and CTRL data were principally insulin responses. This validates our pathway clustering approach (Figure 5A).

Examining the glycolytic metabolites, media exchange prominently increased fructose 1,6-bisphosphate, glyceraldehyde 3-phosphate (G3P), and dihydroxyacetone phosphate levels in insulin-conditioned, but not naive, adipocytes (Figure 6C). These metabolites are upstream of GAPDH and are enriched by exogenous glucose (Figure 6E). This suggested that media exchange caused a constriction of glucose catabolism at GAPDH, and the effects were more apparent when glycolytic flux was high. The reciprocal depletion and accumulation of nucleotide tri- and monophosphates, respectively, (Figure 6C) decreased the energy charge (Figure 6D), supporting an imbalance between upper and lower glycolysis. This resembles the “phosphate stress” phenomenon seen in yeast (van Heerden et al., 2014). In contrast, naive adipocytes did not show these responses to media exchange (Figures 6C and 6D). Overall, this further demonstrated that GAPDH is a flux-controlling step of glycolysis.

Insulin Switched Adipocytes' Substrates to Glucose

Last, we complemented our flux analysis with carbon book-keeping, enabling us to determine the primary fates of glucose. While insulin increased the incorporation of glucose carbon into most metabolites (Figure 7A), the majority ended up as lactate at 58%, followed by glycogen at 17%. Despite being fat cells, the conversion of glucose into lipogenic acetyl-CoA, determined by radiolabeling methods, was trivial at less than 2%, with or without insulin. Lipogenesis did not appear to be a significant sink for NADPH to balance out the large ME flux, although our estimates appeared consistent with previous studies (Crown et al., 2015, Green et al., 2016). Furthermore, glucose oxidation rates (incorporation of 14C-labeled glucose into CO2) were significantly smaller than the lactate production rate (Figure 7B).

Figure 7

Insulin Altered the Profile of Substrates' Fates

(A) The fates of glucose carbon, showing the increased conversion of glucose into lactate and glycogen. Pie charts report median percentages of glucose carbon contained in the various end products.

(B) The partitioning of radiolabeled glucose into end products CO2, glycerol (TAG-Gly), and fatty acid. These pools constituted a small fraction of the glucose consumed compared with lactate efflux determined by enzymatic assay. Data are represented as mean ± SEM of 3 biological replicates.

(C) The simulated redistribution profile of adipocytes' major substrates at 60 min. The partial penetration of glucose into TCA cycle metabolites compared with glycolytic metabolites shows that the increase in glucose oxidation by the TCA cycle was relatively small compared with aerobic glycolysis.

(D) Unstimulated adipocytes showed a mixed substrate profile, with insulin subsequently shifting carbon sources away from glutamine and pyruvate to glucose. The endogenous pool represents all the metabolites present before the addition of glucose label and is a significant carbon source. Pie charts report medians of carbon uptakes.

In contrast, in the unstimulated state, 52% glucose carbon consumed lingered as intracellular glucose (22%), G3P (12%), and other metabolites (18%), with a further 33% converted into lactate (Figure 7A). We then used the flux models to extrapolate for the fates of non-glucose substrates (see Transparent Methods). The results showed that TCA metabolites and lactate came from a variety of sources in unstimulated adipocytes, but these sources were displaced by glucose upon insulin stimulation (Figure 7C). Insulin constricted the uptake of non-glucose substrates, namely, glutamine and pyruvate (Figure 7D). For instance, the shift from glutaminolysis to pyruvate anaplerosis by insulin was indicated by the increase of net pyruvate carboxylation flux (PC minus ME flux) from −1.8 to 1.5 nmol/mg/h (Figure 4B). The large interconversion between pyruvate, malate, and oxaloacetate caused a significant assimilation of unlabeled CO2. Overall, the analysis of substrate fates suggested that insulin treatment reprogrammed adipocytes to preferentially catabolize glucose.

Discussion

In this study, we have developed a 13C-DMFA workflow to estimate non-steady-state fluxes from temporal tracer metabolite data. The incorporation of B-splines and SBR was instrumental to overcome the computational challenges faced when modeling fluxes in central carbon metabolism. The built-in algorithm can model a variety of flux trajectories (e.g., ramp, decline, sigmoidal, and even no change). The provided MATLAB package simplifies the user's task to input manipulation and prototyping (Figures 1E and 3A). Here, we used 13C-DMFA to mechanistically interpret adipocyte metabolite data, showcasing that the adipocyte dynamically coordinates numerous metabolic pathways to efficiently process glucose in response to insulin.

Our analysis revealed a heterogeneous response to insulin, with segments of metabolic pathways being differentially regulated by insulin (Figure 5B). The different flux dynamics between upper and lower glycolysis implicated the role of GAPDH in regulating glycolytic fluxes in an insulin-dependent manner. This corroborates our previous observations that insulin exerts a coordinated response at multiple sites beyond GLUT4 to regulate adipocyte metabolism (Humphrey et al., 2013, Krycer et al., 2017, Ma et al., 2015). It extends the concept of anabolic priming, complementing the demand-driven diversion of glucose into appropriate sinks through fast protein phosphorylation events independently of glucose supply. Altogether, this work demonstrates the value of acute tracer experiments in generating testable hypotheses pertinent to metabolic regulation.

In particular, 13C-DMFA enabled us to compare the speed and magnitude of pathway fluxes and facilitated detailed carbon book-keeping (Figures 4A and 7A). This revealed that glycolysis operated faster than the TCA cycle and that a major fate of glucose was lactate, together suggesting that aerobic glycolysis is a major metabolic trait of adipocyte metabolism. Indeed, a preference for lactate production has been well documented in primary adipocytes (DiGirolamo et al., 1992). Furthermore, aerobic glycolysis is a cancer hallmark (Pavlova and Thompson, 2016, Vander Heiden and DeBerardinis, 2017), and may also be utilized here in terminally differentiated adipocytes to facilitate insulin-stimulated anabolism.

Another rationale is that lactate production enables the rate of glycolysis to be rapidly elevated independently of other pathways (Hui et al., 2017). In our data, glycolysis achieved substantially higher fluxes than the oxidative TCA cycle and peripheral pathways (Figure 4A), with an acceleration matching the dynamics of GLUT4 trafficking (Burchfield et al., 2013) and glucose uptake (Krycer et al., 2017) in response to insulin. Thus, aerobic glycolysis serves to cater for the large glucose influx upon insulin stimulation.

Exploring this further, we challenged insulin-stimulated adipocytes with media exchange, finding a significant drop in energy charge. This demonstrates the importance of prioritizing ATP harvesting before ramping up glycolysis because glycolysis is an autocatalytic cycle that has inherent potential to be unstable (Barenholz et al., 2017). GAPDH regulation is a safeguard against glycolysis stalling (van Heerden et al., 2014), but GAPDH itself was prone to disruption when glycolytic flux was high (Figures 6C and 6D). Parenthetically, although fluxes were perturbed by media exchange, the effect has no bearing on previous observations of metabolic priming based on metabolite data (Krycer et al., 2017), because naive adipocytes 10 min after media exchange stayed on course (Figure 6B).

Collectively, this suggests that adipocyte metabolism is rewired in response to insulin to cater for the rapid influx of glucose, diverting glucose into relevant pathways, preventing ATP stress, and perhaps utilization of lactate dehydrogenase to maintain cytosolic redox status. Future investigations should explore how these metabolic responses by the adipocyte change with varying glucose availability, including insulin resistance where glucose uptake is impaired.

Overall, 13C-DMFA caters for tracer experiments examining non-steady-state metabolism, its utility demonstrated through the temporal resolution of insulin responses at the intracellular flux level. We envisage that 13C-DMFA can be integrated or cross-validated with other datasets, such as gene expression and protein modification data, to dissect how enzyme levels, allostery, and post-translational modifications interact to acutely regulate metabolism.

Limitations of the Study

Our approach provides a substantial improvement on existing approaches to quantify dynamic flux rates. Nonetheless, there are future avenues for improvement. For instance, we did not achieve statistically acceptable fit as the minimum SSRs were greater than the chi-square critical value. This indicated the presence of systematic biases in the model and/or data that were not apparent, or we have underestimated the uncertainties of the experimental data. A potential source of bias was the adoption of a single-compartment model to avoid overfitting, which can impact interpretation of labeling dynamics when metabolite production has significant compartment dependency (Buescher et al., 2015). However, the methodology itself is sound because it can find true solutions (Figure S3D), and the use of AIC to optimize B-spline order and number of knots provided countermeasure against overfitting (Figure 3C). Furthermore, our model can be adapted to improvements in metabolite detection, such as the incorporation of compartmentalization in metabolomics data.

Methods

All methods can be found in the accompanying Transparent Methods supplemental file.