A Mathematical Model for Enzyme Clustering in Glucose Metabolism

We have recently demonstrated that the rate-limiting enzymes in human glucose metabolism organize into cytoplasmic clusters to form a multienzyme complex, the glucosome, in at least three different sizes. Quantitative high-content imaging data support a hypothesis that the glucosome clusters regulate the direction of glucose flux between energy metabolism and building block biosynthesis in a cluster size-dependent manner. However, direct measurement of their functional contributions to cellular metabolism at subcellular levels has remained challenging. In this work, we develop a mathematical model using a system of ordinary differential equations, in which the association of the rate-limiting enzymes into multienzyme complexes is included as an essential element. We then demonstrate that our mathematical model provides a quantitative principle to simulate glucose flux at both subcellular and population levels in human cancer cells. Lastly, we use the model to simulate 2-deoxyglucose-mediated alteration of glucose flux in a population level based on subcellular high-content imaging data. Collectively, we introduce a new mathematical model for human glucose metabolism, which promotes our understanding of functional roles of differently sized multienzyme complexes in both single-cell and population levels.

describe the importance of bistable behavior in glycolysis with strong emphasis on the known allosteric regulation of the rate-limiting enzymes 12 . Particularly, recent mathematical models for human liver or cancer cells cover a complicated network beyond glycolysis by including the pentose phosphate pathway, serine biosynthesis, glutamine metabolism, and/or mitochondrial metabolism 13,14 . However, mathematical models for glucose metabolism have barely explained yet how a multienzyme complex influences glucose flux in cancer cells in a cluster size-dependent manner.
In this work, we have constructed a mathematical model to understand how the formation of a multienzyme complex in glucose metabolism affects glucose flux in cancer cells. Since the glucosome is composed of at least four cytoplasmic rate-limiting enzymes of glucose metabolism 4 , we have formulated a model to accommodate all the members of the glucosome. Our model indeed reveals the impact of their spatial complexation in silico on temporal changes of glycolytic metabolites at a subcellular level in single cells. We further show that our model is capable of predicting metabolic outcomes of a population of cells based on the relative distribution of glucosome clusters in the population. In conclusion, our mathematical model includes not only enzyme kinetics and their allosteric regulation but also spatial compartmentalization of the rate-limiting enzymes of the pathway, thus allowing us to understand the metabolic contributions of glucose flux at both single-cell and population levels.

Model and Methods
A Mathematical Model. We develop a model to investigate how the size of glucosome clusters affects the direction of glucose flux in cancer environment as an indication of the metabolic activity of glucosome clusters in various sizes. Our model involves 7 metabolic intermediates (S i ), 9 enzyme-associated species (E i ) with various forms by which enzyme activities of the rate limiting steps are regulated, 3 metabolic products (P i ), and 28 reactions ( Fig. 1 and Tables 1 and 2). Among the fourteen enzymes participating in glycolysis and/or gluconeogenesis, we only focus on the enzyme activities of the glucosome members (i.e., PFKL, FBPase, PKM2, and PEPCK1), which are hypothesized to have different activity levels in a cluster size-dependent manner. In Fig. 1, we present our simplified metabolic pathway with the enzymes involved in the glucosome. All chemical species and reaction rate constants in our pathway are introduced in Tables 1 and 2. An ordinary differential equations (ODE) model was used to describe temporal dynamics of the metabolic network involving enzymatic reactions associated with glucose metabolism. Temporal concentration changes in each metabolic intermediate are determined by propensity functions of the reactions where the intermediate is produced or consumed. Table 3 summarizes the propensity function R i for the ith reaction. In Table 3, most of the reactions are modeled by the law of mass action, in which the reaction rate is proportional to the substrate concentrations. The allosteric regulations of metabolic enzymes by metabolites are modeled as Michaelis-Menten kinetics, and their propensity functions include rational functions of the substrate concentrations to describe allosteric effects of the corresponding enzymes. Table 4 summarizes a set of all ODEs used in our model.
We also provide reaction rate constants and initial values for the intermediate concentrations used in our model (Tables 1 and 2). We determined the parameters in Tables 1 and 2 using the numerical simulation of the model to be consistent with the qualitative changes in glucose flux in the literature 15 . As shown in Table 1, we assumed that initial concentrations of metabolic intermediates (S i ), products (P i ), PKM2 tetramers (E 3 *) and glycosylated PFKL (E 1 gly ) are at the basal levels and set them as 0.01. We also assumed that no multienzyme complex (E S/M/L ) is formed at the beginning of simulation. An enzyme concentration of each glucosome member is assumed to be conserved, setting as 100. Assuming that all reaction rate constants for the enzymes are in the same order of magnitude, the reaction rate constants of metabolic enzymes are set to 10 as shown in Table 2. However, we slightly modified the values for k 2 , k −2 , k 4 , and k −4 so that glycolysis becomes dominant relative to the metabolic shunt to the pentose phosphate pathway and/or serine biosynthesis in the absence of enzyme clusters (i.e., P 3 ≫ P 1 , P 2 ). Lastly, we assumed that all the metabolic products are produced (k p/s/f ) and degraded (δ p/s/f ) in the same rates, thus setting them to 5 and 0.5, respectively. Note that we provided both forward and backward reaction rate constants in Table 2, which are separated by comma.
Next, different levels of the enzyme activity depending on the size of glucosome clusters are expressed in terms of the parameters c i (for medium-sized enzyme clusters) and e i (for large-sized enzyme clusters) in Table 2. Based on our hypothesis that the medium-sized enzyme clusters promote glucose shunt to the pentose phosphate pathway 4 , we decelerated glycolysis (c 2 , c −d ≪ 1) and accelerated gluconeogenesis (c −2 , c −6 ≫ 1). Similarly, based on our hypothesis that the large-sized enzyme clusters will shunt glucose flux to serine biosynthesis 4 , we changed parameters to direct glucose flux toward serine production (e 2 , e −6 ≫ 1 , and e −2 , e −d ≪ 1). An enzyme activity in the small-sized clusters is assumed similar to the activity of a freestanding enzyme because the small-sized clusters formed by PFKL do not represent a multienzyme complex 4 . Using the initial values and the reaction rate constants shown in Tables 1 and 2, our mathematical model was simulated using MATLAB to compute temporal changes of glucose flux in different scenarios at the subcellular level. Note that our MATLAB simulation codes are provided in the supplementary document.
Importantly, we have also performed sensitivity analysis of the parameters to investigate how the concentrations of three metabolic products (P i ) are affected by small changes in each parameter that we have incorporated in our model. In the sensitivity analysis, we have used the method of partial rank correlation coefficient (PRCC) 16 to evaluate the level of the linear association between the input values of the parameters and the output concentrations of the metabolic products. To perform the sensitivity analysis, we modified the publically available MATLAB codes which are developed by Kirschner and coworker 16 . The default parameter values were set as the values used in our model, and we defined an interval for each parameter from half of the default value to twice of the default value. Then, we used Latin hypercube sampling 16 to choose each parameter value in the corresponding interval by assuming that each parameter is uniformly distributed in the interval. Afterward, the sampled parameter values were used to calculate the concentrations of the metabolic products at the end time of the simulation. We repeated the sampling process described above 20,000 times and calculated the Spearman correlation coefficients, which correspond to PRCCs.
Experimental Methods. Materials. The plasmid expressing human PFKL with a monomeric form of enhanced green fluorescent protein (mEGFP) was prepared previously (PFKL-mEGFP) 4 . We have employed PFKL-mEGFP as a marker of the glucosome. 2-Deoxyglucose was purchased from Sigma.
Transfection. To prepare cells for transfection and subsequent imaging, Hs578T cells were gently removed from the culture flask by replacing the culture medium with Trypsin-EDTA solution (Corning, Cat# 25-053-Cl). Fresh, antibiotic-free growth media was subsequently used to harvest and resuspend the cells, followed by plating on glass-bottomed 35 mm Petri dishes (MatTek). When the confluency became ~70-90% on the following day, the cells were transfected with Lipofectamine 2000 (Invitrogen) using the Opti-MEM-I reduced serum media Figure 1. Simplified glucose metabolism with multienzyme complexes. Seven metabolic intermediates are involved in the pathway: S 1 represents glucose, S 2 is fructose-6-phosphate, S 3 is fructose-1,6-bisphosphate, S 4 is 3-phosphoglycerate, S 5 is phosphoenolpyruvate, S 6 is pyruvate, and S 7 is oxaloacetate. Four rate-limiting enzymes are E 1 (phosphofructokinase 1, PFK), E 2 (fructose-1,6-bisphosphatase, FBPase), E 3 (pyruvate kinase M2 dimer, PKM2), and E 4 (phosphoenolpyruvate carboxykinase 1, PEPCK1). Pyruvate kinase M2 catalyzes conversion from S 5 to S 6 when it becomes a tetramer ( ⁎ E 3 ). On the other hand, phosphofructokinase 1 is inactivated after post-translational glycosylation (E 1 gly ). PFK forms three differently sized clusters: E S , E M , and E L represent small-, medium-and large-sized clusters, where E M and E L are multienzyme complexes. To measure the direction of glucose flux, we denote three metabolic products as P 1 , P 2 , and P 3 , which represent metabolic outcomes of the pentose phosphate pathway, serine biosynthesis and the downstream of glycolysis. All the used parameters are summarized in Tables 1 and 2 (Opti-MEM-I; Gibco, Cat# 11058). The Opti-MEM-I medium was then exchanged with the fresh antibiotic-free growth medium after a 5 h incubation (37 °C, 5% CO 2 , and 95% humidity), followed by ~18-24 h incubation in the CO 2 incubator at 37 °C.
Fluorescence Live-cell Imaging. On the day of imaging (~18-24 h post-transfection), cells were washed with imaging solution (20 mM HEPES (pH 7.4), 135 mM NaCl, 5 mM KCl, 1 mM MgCl 2 , 1.8 mM CaCl 2 , and 5.6 mM glucose) for three 10 min incubations, followed by a ~1-2 h incubation at ambient temperature. All samples were then imaged at ambient temperature (~25 °C) with a 60× 1.45 NA objective (Nikon CFI Plan Apo TIRF) using a Photometrics CoolSnap EZ monochrome CCD camera on a Nikon Eclipse Ti inverted C2 confocal microscope. Epifluorescence imaging was carried out using the following filter set from Chroma Technology; mEGFP detection by a set of Z488/10-HC cleanup, HC TIRF Dichroic and 525/50-HC emission filter.
For cell-based high-content imaging assays, 2-deoxyglucose (25 mM) was added to Hs578T cells after washing three times with the imaging solution. Images showing PFKL-mEGFP clusters were acquired before and after cells were incubated with 2-deoxyglucose at various time points (up to 6 hours). Control experiments were also carried out with 167 μL of vehicle (i.e., water). To ensure reproducibility, our experiments were repeated at least five times over the course of a few months. Statistical analysis was performed using two-sample two-tail t-tests.
Cluster Size Analysis. Cluster size analysis was accomplished using the ImageJ processing software (National Institutes of Health) as we have done before 4 . Briefly, fluorescent wide-field images were processed through ImageJ using a custom script and macro that automates the counting of fluorescent clusters using its built-in module, so-called robust automatic threshold selection (RATS). In this analysis, the images were scaled according to the pixel size of the microscope (i.e., 0.12 µm/pixel) before the default parameters for RATS (i.e., noise threshold = 25, λ factor = 3) were used in this analysis. Once fluorescent clusters were selected from an image, the particle analysis module was applied to attain both the number and area of fluorescent clusters within an image. This process was repeated for all subsequent cell images. The operator then evaluated the original cell images against the particle mask to eliminate data in which more than one cluster was counted as a single particle.

Results
Modeling Glucose Metabolism with Multienzyme Complexes. Glucose metabolism consists of glycolysis and gluconeogenesis where three irreversible reactions are involved in glycolysis, four irreversible reactions in gluconeogenesis, and seven reversible reactions in both pathways. In our model, we simplify the 14-step pathway by condensing a few reversible steps into one conversion (Fig. 1). Accordingly, we omit some of the metabolic intermediates such as glucose-6-phosphate, dihydroxyacetone phosphate, glyceraldehyde-3-phosphate, 1,3-bisphosphoglycerate, and 2-phosphoglycerate because they are not directly shuttled to other metabolic pathways. In addition, three metabolic products (P i ) represent three metabolic fates of glycolytic metabolites. Importantly, we focus on four cytoplasmic enzymes which spatially organize into multienzyme clusters (i.e., glucosome clusters) 4 : PFKL, FBPase, PKM2, and PEPCK1. Note that the other ten enzymes are not cytoplasmic or rate-determining in glucose metabolism.
In addition, we incorporated cancer-relevant mechanisms of the four enzymes forming the glucosome. First, pyruvate kinase catalyzes conversion from phosphoenolpyruvate to pyruvate. The M2 isoform of pyruvate kinase

Variables
Chemical species Values (non-dimensional) (PKM2) is predominantly expressed in cancer cells while suppressing the expression of pyruvate kinase M1 isoform (PKM1) 17 . PKM2 plays an important role in aerobic glycolysis as a dimer by impairing pyruvate production, thus redirecting glucose flux into serine biosynthesis 18,19 . However, a tetrameric form of PKM2 possesses its canonical glycolytic activity in production of pyruvate in cancer cells. We take into account the dimer-tetramer conversion as part of the regulation of glucose metabolism (Fig. 1). Second, PFK catalyzes the conversion of fructose-6-phosphate to fructose-1,6-bisphosphate. Particularly, liver-type PFK (PFKL) is glycosylated in several types of cancer cells under hypoxia conditions, diverting the direction of glucose flux into the pentose phosphate pathway 20 . Our model includes the impact of glycosylation on PFKL as part of the regulation of its metabolic activity in glycolysis. Third, cytoplasmic PEPCK1 plays an essential role in the first step of gluconeogenesis 21 . Recently, a new role of PEPCK1 in tumor cell proliferation was discovered in certain types of cancer cells by controlling carbon metabolic flux through the TCA cycle 22 . However, a molecular level mechanism of PEPCK1 in cancer cells (e.g., post-translational modification-dependent activity or oligomerization-dependent activity) is not firmly established yet to be included in our model. Fourth, FBPase catalyzes the conversion of fructose-1,6-bisphosphate to fructose-6-phosphate. Since FBPase catalyzes the backward reaction of PFK, the activities of PFK and FBPase are reciprocally regulated by various allosteric metabolites. Similarly, an antagonistic role of FBPase for glycolytic flux was validated to play a role in kidney cancer progression 23 . However, other than the well-established allosteric regulation, cancer-specific alteration of FBPase is not established yet to be included in our model. Unlike many previous mathematical models for glycolysis, we take into account the recent discovery of the glucosome being formed in various sizes in cancer cells. Briefly, various sizes of PFKL clusters in cancer cells are categorized into three subclasses; small, medium and large-sized cluster 4 . Small-sized clusters formed by PFKL are single-enzyme assemblies whereas medium-and large-sized clusters represent spatial organizations of multienzyme complexes. When cancer cells display small-sized clusters, the given cells barely include the other  Glucose Flux Analysis at Subcellular Levels. In Fig. 2, metabolic product concentrations are shown in different experimental conditions, which are computed by numerically solving the mathematical model given in Tables 3 and 4. The model results in time-dependent changes of three metabolic products: P 1 represent a metabolic product of the pentose phosphate shunt (green line), P 2 is a metabolic product of serine biosynthesis (red line), and P 3 indicates a metabolic product of the downstream of glycolysis (blue line). In addition, the activity of non-clustering PFKL or clustered PFKL into small sizes are anticipated to be similar in our model. For the simulation of the cases with no cluster or only small-sized PFKL clusters, we set all the parameters related to the formation of medium or large-sized enzyme complexes as zero (i.e., k am , k −am , k al , k −al , c i , and e i ). Please, note that the parameters for small-sized PFKL clusters (i.e., k as and k −as ) were also set to zero when we simulated the case showing no cluster. Consequently, our simulation for no cluster or only small-sized PFKL clusters showed the high level of P 3 relative to those of P 1 and P 2 ( Fig. 2A and B), indicating that most of glucose flux flows to glycolysis to produce pyruvate and beyond. However, when the rate-limiting enzymes in glucose metabolism are spatially organized into medium-sized glucosome clusters, the enzymes in medium-sized clusters may have different levels of catalytic activity. Because the increased flux of the pentose phosphate shunt is correlated with the high level of medium-sized clusters in cancer cells 4 , our assumption of c 2 , c −d ≪ 1 and c −2 , c −6 ≫ 1 resulted in a significant increase of P 1 but decreased the concentrations of P 2 and P 3 (Fig. 2C). In this case, we set the values of k al , k −al , and e i as zero. Lastly, we considered all the rate-limiting enzymes being assembled into large-sized Table 3. Propensities of 28 reactions in the glucose metabolic pathway. , P = (P 1 , P 2 ,    Table 4. A system of ODEs for metabolic intermediates, enzymes and their clusters, and metabolic products.

Equations
glucosome clusters in the presence of smaller clusters. Because the promotion of large-sized clusters in cancer cells is correlated with the increased flux of serine biosynthesis 4 , our assumption of e 2 , e −6 ≫ 1 and e −2 , e −d ≪ 1 along with the same assumption of c i in the previous case resulted in the substantially increased level of P 2 relative to the levels of P 1 and P 3 (Fig. 2D). Collectively, our mathematical model is developed to understand the cluster size-dependent changes of glucose flux in cancer cells.
Sensitivity Analysis. Next, we performed a sensitivity analysis by which we evaluated without bias how the input values of simulation parameters influence the production of specific metabolic flux (P i ). According to the publically available MATLAB codes developed by Kirschner and coworker 16 , we modified to calculate PRCCs in our model ( Figure S1). If the PRCC values are close to −1 or +1, the corresponding parameters have strong correlations with the formation of the given metabolic product. On the other hand, if the correlation coefficient is close to 0, there is no correlation between the input parameters and the given metabolic product. When the absolute values of the PRCCs are greater than ±0.2 and their p-values are less than 0.01, we defined the corresponding parameters as 'sensitive' parameters ( Fig. 3) to determine which simulation parameters are critical for our model and how sensitive each parameter is to understand its contribution to the outcomes of glucose flux. Noticeably, the formation of each metabolic product was sensitive not only to the enzyme activities of the glucosome members but also the efficiencies of enzyme clustering into the certain sizes of glucosomes. In Fig. 4, we further analyzed a few critical parameters that influence the enzyme activities of the glucosome members (i.e., k 2 , k −2 , and k −d ). Each selected parameter was perturbed in one direction while all other parameter values were fixed as default values as shown in Table 2. The changes in the metabolic product concentrations were then presented under four cases; (i) no spatial organization, (ii) small-sized PFKL clusters, (iii) medium-sized multienzyme complexes, and (iv) large-sized multienzyme complexes. Note that for the cases dealing with the medium-sized or large-sized complexes, the smaller-sized clusters are also formed in the model. When k 2 was decreased from 40 to 10 (Fig. 4A), the conversion from S 2 to S 3 was reduced. We thus observed significant increase in P 1 while decrease in P 2 and P 3 . On the other hand, when k −2 was increased from 7 to 10 (Fig. 4B), a similar pattern was observed as in Fig. 4A although the changes were small due to the narrow ranges of the parameter. Finally, when k −d was increased from 1 to 10 (Fig. 4C), glucose flux moved toward the production of P 3 , which resulted in significant increase in P 3 and decrease in P 1 or P 2 . Note that we did not observe much increase of P 3 in the cases where no spatial organization or only small-sized PFKL clusters are dominant (Fig. 4C) because the flux was already committed toward the production of P 3 in these cases. Collectively, the simulated flux changes reflect the changes of enzymatic activities involved in glucosome clusters.
In Fig. 5, we also analyzed more critical parameters that influence the efficiencies of enzyme clustering into the glucosomes (i.e., c −2 , e 2 , and e −d ). Similarly, each selected parameter was perturbed in one direction while fixing the other parameter values. The changes in the metabolic product concentrations were shown under the four cases as well. When c −2 was decreased from 10 to 1 (Fig. 5A), the level of P 1 decreased and that of P 2 increased significantly due to the decreased contribution of the medium-sized cluster to the conversion from S 3 to S 2 . When e 2 was decreased from 2.5 to 1 (Fig. 5B), the resulting effect was quite opposite from the one shown Fig. 5A because the decreased e 2 reduces the contribution of the large-sized clusters to the conversion of S 2 to S 3 . When e −d was increased from 0.05 to 1 (Fig. 5C), the level of P 2 decreased and that of P 3 increased significantly with Figure 3. The partial rank correlation coefficients (PRCCs) between the 'essential' input parameters and the concentrations of metabolic products. The PRCCs at time = 10 are graphed to provide relative strengths of the correlations between the input parameters and the concentrations of metabolic products (P 1 , P 2 , and P 3 ). The horizontal lines at ±0.2 indicate the thresholds we used to distinguish sensitive essential parameters from non-essential parameters. The PRCCs of all parameters are shown in supplementary figure S1. P 1 , P 2 , and P 3 represent metabolic outcomes of the pentose phosphate pathway, serine biosynthesis and the downstream of glycolysis, respectively.
Scientific REPORtS | (2018) 8:2696 | DOI:10.1038/s41598-018-20348-7 slight decrease in P 1 . This is due to the increased efficiency of the large-sized cluster on its contribution to the conversion of E 3 to * E 3 .

Metabolic Flux Analysis at Ensemble Levels.
We have further evaluated our model by quantifying the changes of glucose flux in the presence of methylene blue, fructose-1,6-bisphosphate or epidermal growth factor (EGF), which are known to regulate the direction of glucose flux at ensemble levels [24][25][26][27][28][29] . As shown in Table 5, we have counted the percentages of cells showing no cluster, small, medium, and large clusters in the absence and presence of the listed glucose flux regulators. In combination with our subcellular flux analysis (Fig. 2), we have computed the time-dependent changes of metabolic products (P 1 , P 2 , and P 3 ) at an ensemble level in the four published cases 4 . To predict the levels of metabolic products at an ensemble level, we multiplied the levels of P 1 , P 2 , and P 3 of the four cases shown in Fig. 2 by the corresponding proportions of the cells showing differently sized clusters. In a control population of Hs578T cells without any exogenous stimulus (Fig. 6A), the level of P 3 appears to be slightly greater than those of P 1 and P 2 . However, relative to Fig. 6A, the cancer cells that were treated with methylene blue (Fig. 6B) or fructose-1,6-bisphosphate (Fig. 6C) increase the level of P 1 while decreasing the level of P 3 at an ensemble level. On the other hand, Fig. 6D shows the case where the cancer cells have been treated with EGF. In this case, the level of P 2 increases but the level of P 1 decreases relative to the other cases. It appears clear that the simulation is capable of showing that glucose flux shunt from glycolysis to either the pentose phosphate pathway or the serine biosynthesis at ensemble levels, whose trend is indeed consistent with the already known functions of the glucose flux regulators in populations [24][25][26][27][28][29] . Collectively, our mathematical model with high-content imaging analysis is adequate to predict metabolic consequences of glucose flux regulators in a population of cancer cells.

Mathematical Prediction of Metabolic Flux with 2-Deoxyglucose at Ensemble Levels.
Our model is now applied to assess metabolic consequences of 2-deoxyglucose in Hs578T cells. 2-Deoxyglucose is a competitive inhibitor of phosphoglucose isomerase catalyzing step 2 in glycolysis 30 . However, due to the lack of 2-hydroxyl group, 2-deoxyglucose cannot further metabolize into fructose-6-phosphate in downstream glycolysis, thereby resulting in the inhibition of glycolysis. Accordingly, 2-deoxyglucose appeared to impair cell growth and thus implemented for anti-tumor therapeutics. However, its pharmacological effects on cancer treatment [31][32][33] remain inconclusive, consequently raising a number of questions particularly including potential metabolic alterations by 2-deoxyglucose in cancer cells.
In this work, we have employed our mathematical model to predict whether 2-deoxyglucose can potentially alter the direction of glucose flux in cancer cells beyond glycolysis inhibition. First, we performed cell-based high-content imaging assays to validate the metabolic effect of 2-deoxyglucose at single cell levels (Fig. 7). Briefly, we expressed PFKL-mEGFP as a glucosome marker in Hs578T cells. To be consistent with Table 5, Hs578T cells were cultured in RPMI1640 with 10% dialyzed FBS. When we treated transfected Hs578T cells with 25 mM of 2-deoxyglucose for 6 hours, the percentage of cells showing small-sized clusters was significantly reduced from 58.3 % to 34.7% (Fig. 7A). Meanwhile, the percentages of cells showing medium-and large-sized clusters were significantly increased from 13.4% to 21.2% and from 26.7% to 44.1%, respectively (Fig. 7A), indicating the strong promotion of glucosome formation and thus metabolic shunts to anabolic pathways.
Subsequently, we computed an ensemble-level outcome of glucose flux in the presence 2-deoxyglucose. Similarly to Fig. 6, we calculated the changes of metabolic outcomes at ensemble levels based on our subcellular   flux analysis (Fig. 2) and high-content imaging analysis (Fig. 7A, Table 5). Relative to a control flux in Hs578T cells where glycolytic flux reaches to ~7 arbitrary units at t = 10 (Fig. 6A), glycolytic flux was indeed inhibited at ~ 5 arbitrary units in the presence of 2-deoxyglucose (Fig. 7B). Importantly, our model provided that metabolic shunts from glycolysis to anabolic biosynthetic pathways were promoted and relatively became dominant in the population of Hs578T cells. Although more experimental validation may be necessary, it is clear that 2-deoxyglucose promoted the formation of glucosome clusters in both single-cell (Fig. 7C,D) and ensemble levels (Fig. 7B). Therefore, we propose that along with the inhibitory role of 2-deoxyglucose in glycolysis, its commitment diverting glucose flux into anabolic pathways may explain in part cancer progression observed during its clinical trials 31 .

Discussion
We constructed a simple mathematical model to understand the cluster size-dependent functional contributions of metabolic enzymes and their multienzyme complexes to cancer cell metabolism. Briefly, glycolysis is interconnected with energy metabolism and anabolic biosynthetic pathways, including the pentose phosphate pathway and serine biosynthesis. Three metabolic products (P i ) thus represent downstream pathways of glycolysis, allowing us to investigate how glucose flux changes its direction at metabolic nodes between energy metabolism and anabolic pathways in cancer cells. Based on our experimental data 4 , we developed a model to predict that The graph shows the average percentages (%) of cells displaying the given sized clusters along with their standard deviations (±) in the absence (black bars) and presence (red bars) of 2-deoxyglucose. At least five independent imaging sessions were performed and total 1200 transfected cells were analyzed. Statistical analyses were performed using two-sample two-tail t-tests. *p < 0.01. (B) Metabolic flux analysis were performed using our mathematical model with the high-content imaging data of 2-deoxyglucose. Relative to a control flux ( Fig. 2A), glycolytic flux (P 3 ) decreased, but the metabolic shunts of glucose to the pentose phosphate pathway (P 1 ) and serine biosynthesis (P 2 ) increased. (C and D) Representative images of Hs578T cells show subcellular localization of PFKL-mEGFP before and after treatment of 2-deoxyglucose (25 mM) for 6 hours. Scale bar, 10 µm.
medium-sized clusters of glucosomes shunt glucose flux into the pentose phosphate pathway whereas large-sized clusters of glucosome divert glucose flux into serine biosynthesis at subcellular levels. Importantly, our model supports that the changes of relative ratios of cancer cells displaying small-, medium-and large-sized clusters in a population appear to be significant enough to influence overall net metabolic outcomes of cancer cells at ensemble levels. Moreover, our mathematical model is further evaluated to predict the effect of 2-deoxyglucose on the fate of glucose, thus providing new quantitative insights of how 2-deoxyglucose alters glucose flux in cancer cells.
Collectively, we conclude that our mathematical model supports the hypothesis that glucosomes divert glucose flux into the pentose phosphate pathway and serine biosynthesis in a cluster size-dependent manner.
It is important to emphasize here that our mathematical model accounting for size-dependent metabolic functions of glucosome clusters is significantly different from other models that mostly rely on in vitro enzyme kinetics. Our model explains altered glycolysis in cancer cells that redirects glucose flux into the pentose phosphate pathway and serine biosynthesis in a cluster size-dependent manner. Importantly, we can now predict the direction of glucose flux in the presence of small-molecule drug candidates as long as quantitative high-content imaging data is obtained in single-cell levels. Additionally, we can quantify relative partition ratios of glucose flux between glycolysis, the pentose phosphate pathway and serine biosynthesis in various conditions. Collectively, our mathematical model fully integrates various cancer-associated mechanisms discovered in recent years, thus advancing our understanding of cancer cell metabolism in single cells.