Abstract
Background
Metabolic flux profiling based on the analysis of distribution of stable isotope tracer in metabolites is an important method widely used in cancer research to understand the regulation of cell metabolism and elaborate new therapeutic strategies. Recently, we developed software Isodyn, which extends the methodology of kinetic modeling to the analysis of isotopic isomer distribution for the evaluation of cellular metabolic flux profile under relevant conditions. This tool can be applied to reveal the metabolic effect of proapoptotic drug edelfosine in leukemia Jurkat cell line, uncovering the mechanisms of induction of apoptosis in cancer cells.
Results
The study of ^{13}C distribution of Jukat cells exposed to low edelfosine concentration, which induces apoptosis in ≤5% of cells, revealed metabolic changes previous to the development of apoptotic program. Specifically, it was found that low dose of edelfosine stimulates the TCA cycle. These metabolic perturbations were coupled with an increase of nucleic acid synthesis de novo, which indicates acceleration of biosynthetic and reparative processes. The further increase of the TCA cycle fluxes, when higher doses of drug applied, eventually enhance reactive oxygen species (ROS) production and trigger apoptotic program.
Conclusion
The application of Isodyn to the analysis of mechanism of edelfosineinduced apoptosis revealed primary druginduced metabolic changes, which are important for the subsequent initiation of apoptotic program. Initiation of such metabolic changes could be exploited in anticancer therapy.
Background
The characterization of metabolic flux profile in living cells is an important issue in understanding the regulation of normal metabolism and the development of disease processes. Such characterization is then necessary for the development of novel therapeutic strategies. Stable isotope tracing using [1,2^{13}C_{2}]glucose as a source of carbon, has been described as a very powerful tool for metabolic flux profiling [16]. The specific pattern of various ^{13}C isotopic isomers (isotopomers) fractions measured using mass spectrometry or nuclear magnetic resonance techniques characterized the distribution of metabolic fluxes in the cells under the studied conditions. To evaluate the flux distribution from measured isotopomer distribution a special software tool is necessary.
Classical ^{13}C metabolic flux analysis (MFA) evaluated steady state metabolic fluxes based on isotopomer fractions measured under the conditions of isotopic steady state [7]. For nonstationary metabolic flux analysis we developed a tool called "Isodyn" (from "isotopomer dynamics") [810] that simulates ^{13}C redistribution in metabolites by automatically constructing and solving large systems of differential equations for isotopomers. Although intracellular metabolites could reach isotopic steady state in a range of minutes, the existence of intracellular stores essentially delays the time necessary for establishing isotopic steady state. Such stores as glycogen, aminoacids and lipids, which intensively exchange with intermediates of central carbohydrate metabolism, could prolong the presteady state phase for all isotopomers. Of course, there is always a possibility of measuring the labeling of such stores and apply classical ^{13}CMFA for the "fast" intermediates of central metabolism considering that they are in quasisteady state. However, the simulation of such "slow" variables provides additional information for the determination of the characteristics of the system.
Moreover, there is another reason for using nonstationary analysis based on a kinetic model of considered pathways: it allows, when experimental data is enough, a more profound analysis of kinetic characteristics and regulation in the pathway. Such advantages stimulated the development of other bioinformatic tools for nonstationary flux analysis [11]. Here, an application of Isodyn for revealing the characteristics of cancer cell metabolism and their change induced by a proapoptotic agent edelfosine is described.
Apoptosis is a programmed cell death and the evasion of apoptotic programm is one of the most fundamental characteristics of cancer cells [12]. However, transformed cells still possess the components of apoptotic mechanism, and it could be induced by various agents. The strategy of selectively killing tumor cells by inducing apoptosis could be used for cancer therapy [13,14], and the presented analysis provides information for the development of such strategy.
Apoptotic process is a complex sequence of signaling events and metabolic changes. The cascade of signaling events resulting in cell death is well studied. However, the signals to apoptosis could be seen as a result of severe distortions in metabolism. In this way, the metabolic changes could be primary events that activate or inhibit apoptotic process. For example, the stimulation of mitochondrial metabolism related to reactive oxygen species (ROS) production [1517] or the inhibition of glycolisis [18] has been linked with activation of apoptotic cascade. Our objective was to understand whether relevant metabolic changes precede the development of apoptosis, or they just follow the progression of the apoptotic signaling program. To reveal the early metabolic changes, the metabolic effects of very low doses of edelfosine, which induce apoptosis in less than 5% of cellular population, were studied.
Synthetic antitumour ether phospholipid edelfosine (1Ooctadecyl2Omethylracglycero3phosphocholine, ET18OCH_{3}) selectively induces apoptosis in cancer cells [1925]. The cellkilling mechanism of edelfosine is mediated by signalling events such as blocking some protein kinases [26] or activation of specific apoptotic receptors [21]. Also edelfosine induces the increase in mitochondrial reactive oxygen species (ROS) production [20,27], which could be a consequence of certain metabolic distortions. If metabolic changes are primary with respect to the development of apoptotic program, it could be expected that essential changes in cell metabolism could take place at low doses of such drug, which hardly induce apoptosis.
In order to find the metabolic changes caused by the low doses of edelfosine, Isodyn simulated the isotopomer distribution using the available enzyme kinetic information and the experimentally acquired ^{13}C isotopomer distribution data. This approach allowed us to obtain sets of fluxes in the modeled metabolic network, which were consistent with the experimental distribution of mass isotopomers derived from labeled glucose in the presence of edelfosine. The determination of the metabolic conditions that promote apoptosis could be an essential contribution to the therapy of cancer.
Results and Discussion
Glucose consumption and lactate production
The rates of glucose consumption normalized per cell volume (for the convenience of comparison with the other intracellular fluxes) were defined taking into account the change of cell number and concentrations of medium glucose and lactate as it is described in "Methods". These fluxes are summarized in Table 1. These values calculated directly from experimental data were used in the fitting of isotopic isomer distribution described below.
Table 1. The metabolic fluxes of glucose consumption (J_{Glc}) and lactate production (J_{lac}).
Measured isotopomer distributions
Since the objective of this work was to register metabolic changes that precede the development of apoptosis, very low doses of apoptosis inducing drug edelfosine were used (0.5 and 1 μg × mL^{1}) for ^{13}C metabolite distribution experiments. After 48 hours of treatment, the dose of 0.5 μg × mL^{1 }induced less than 1% of apoptosis whereas the dose of 1 μg × mL^{1 }induced between 4 to 5% of apoptosis (measured as subG1 population with respect total number of cells). Incubation with [1,2^{13}C_{2}]glucose as a tracer resulted in a specific isotopomer (^{13}C isotopic isomer) distributions in metabolites, which was measured by GC/MS techniques in medium lactate and ribose5phosphate (r5p) derived from intracellular RNA. The applied technique determined the fractions of different mass isotopomers: m0 (without any ^{13}C labels), m1 (with one ^{13}C label), m2 (with two ^{13}C labels), etc. These fractions for control and edelfosinetreated cells are shown in Table 2.
Table 2. Isotopomer distribution in lactate secreted into the medium and RNA ribose.
Simulation of measured isotopomer distributions
The distribution of isotopomer fractions of a metabolite contains information about the fluxes through the metabolic pathways of its formation. Roughly, when [1,2^{13}C_{2}]glucose is metabolized, anaerobic glycolysis produces mainly m2 lactate, whereas passage of labeled glucose through the TCA cycle (including anaplerotic pyruvate carboxylation, pyr→oaa) or the oxidative and non oxidative branches of the pentose phosphate pathway results in m1 and m3 lactate. Thus, the fractions of m1 and m3 with respect to m2 characterize the contribution of the TCA cycle and pentose phosphate pathway. The simulation and fitting the measured isotopomer distribution allows the determination of a set of metabolic fluxes, which correspond to the measured isotopomer distribution. The details of such determination are described in method section and in [810] and Additional file 1.
Additional file 1. Supplementary materials. PDF containing supplemental figures and other information.
Format: PDF Size: 478KB Download file
This file can be viewed with: Adobe Acrobat Reader
The results of data fitting by Isodyn are also presented in Table 2. The quality of fit is characterized by χ^{2}, the sum of squares of differences between experimental and simulated data normalized by the standard deviations of experimental data. The ribose used for the analysis was extracted from cellular RNA. Thus, the isotopomer distribution in ribose contains information on both the label isotopomer distribution of the de novo synthesized nucleotides and also on the fraction of initial nonlabeled nucleotides that were reused. The program calculates this initial fraction with respect to the one synthesized de novo during the treatment (as described in Additional file 1; it is referred in the tables as "dilution" and characterized RNA synthesis de novo. According to the data of Table 2, in edelfosinetreated cells dilution decreased, which indicates that a greater fraction of RNA was synthesized de novo.
Analysis of metabolic flux profiles
Table 3 show the fluxes corresponding to the best fit shown in Table 2 and indicates the fluxes for which the difference between treated and nontreated cells are statistically significant. According to the table, to fit the measured isotopomer distribution in cell population where edelfosine induced ≤5% of apoptosis, glucose consumption (via hk, hexokinase activity) must increase, the TCA cycle (via pdh, pyruvate dehydrogenase activity) must be activated and pentose phosphate pathways must be inhibited with respect to the control.
Table 3. Simulated metabolic fluxes corresponding to the best fit presented in Table 2.
Table 3 illustrates that a small change in the distribution of mass isotopomers shown in Table 2 could be a consequence of large changes in metabolic fluxes. Specifically, the flux through the TCA cycle (pdh, citrate synthase, flux from citrate to malate) increases almost three folds although the stimulation of apoptotic program can be measured in only 45% of cells. These fluxes normalized per respective glucose uptake are increased also, although not so tremendously. Thus, the low doses of edelfosine activate the whole central metabolism and even more activate the TCA cycle. In fact, it is not so simple to decide what value, normalized or not normalized, characterizes the TCA cycle activation better. Although glycolysis provides substrates for the TCA cycle, it is known that activation of glycolysis is not necessary coupled with the activation of the TCA cycle. For instance, in muscle cells starting active contractions, a hundred fold increase in glycolysis hardly activativates TCA cycle [28]. Glycolysis has much more capacity for activation, while the activation of the TCA cycle coupled directly with mitochondrial bioenergetics requires much more structural changes. If the TCA cycle is activated without the respective increase in the volume occupied by mitochondria, this activation probably could have negative consequences for cell survival.
Despite the changes in isotopomer distribution induced by a low dose of edelfosine is small, χ^{2 }criterion is sufficiently sensitive to them. The program fits the data for control with very small deviations (χ^{2 }= 0.18). Howener, if this "control" set of parameters is used to simulate the data for treated cells, χ^{2 }increases to 60, which indicates that the model well accepted as a simulator of metabolic fluxes in control cells becomes unacceptable for the edelfosinetreated cells. Subsequent fitting procedure decreases χ^{2 }to 5.75 (Table 2) just by increase of metabolic fluxes through the TCA cycle and changes in the pentose phosphate pathway as Table 3 indicates.
The activation of the TCA cycle, revealed for a low dose of edelfosine, could be a reason for the activation of apoptotic process when higher doses of drug are applied. The main function of the TCA cycle in energy metabolism is to produce substrates for mitochondrial respiration. Therefore, the activation of TCA cycle favors the reduction of electron transporters. This is a factor that could switch the mitochondrial respiratory chain to the state of damaging generation of reactive oxygen species (ROS) [29], which is an important component of apoptotic process. The low doses used for the metabolomic studies did not increase ROS production, although their metabolic effect was significant. The increase of dose of edelfosine eventually produce significant increase of ROS, as it is shown in Figure 1. This observation, which is in line with other experimental studies reported elsewhere [20,27], validates the conclusion based on the isotopomer distribution data analysis.
Figure 1. Increase in DCF fluorescence in edefosinetreated Jurkat cells. The DCF fluorescence as a result of ROS production was measured in a fluorescence microplate reader (top panel) and by fluorescence microscopy (left images, bottom two panels). Also shown are the phasecontrast images superimposed with the fluorescence images of the Jurkat cells (right images, bottom two panels).
Another important distinction in the revealed flux profiles of control and edelfosinetreated Jurkat cells is the difference in the fluxes of pentose phosphate pathway. In control cells all pentose phosphate fluxes were higher than those in edelfosinetreated cells (Table 3). Moreover, the model predicted an increase in ribose consumption for RNA synthesis in edelfosinetreated cells, which led to a decrease in ribose 5phosphate concentration, such that its synthesis became practically unidirectional. Thus, newly synthesized ribose 5phosphate was consumed for RNA synthesis to a greater extent in edelfosinetreated cells than control (see flux r5p, Table 3), in accordance with the value of ^{13}C dilution in r5p (see Table 2). In contrast, in control cells, the newly synthesized ribose 5phosphate is not massively used for RNA synthesis, but reutilized in the central energetic metabolism.
Conclusions
The simulation of changes in isotopomer distribution induced by low doses of edelfosine revealed essential activation of the TCA cycle (including anaplerosis). Such direction of changes in metabolism could enhance ROS production when the higher doses of drug are applied. Indeed, when doses of edelfosine causing >5% of apoptosis rate were used (from 1 to 2 μg × mL^{1}) the increase of ROS production becomes measurable experimentally.
Another effect of the drug is the stimulation of RNA synthesis. Both effects of the drug are consistent: the increase of the TCA cycle fluxes are aimed in the increase of ATP energy production necessary for biosynthetic processes induced by such stress factor as edelfosine.
Methods
Structure of program and algorithms for isotopomer distribution analysis
The computer program "Isodyn", which we developed in C++, represents a simulation environment for the dynamics of metabolite labeling by ^{13}C isotopes in metabolic reactions of living cells. For such simulations it uses a classical kinetic model of metabolic pathways linked with a module that computes the distribution of ^{13}C isotopic isomers of metabolites. For the case of metabolic steady state it uses following algorithm for the simulation of dynamics of isotopomer distribution in metabolites.
1. To simulate reaching steady state in the kinetic model for total metabolite concentrations and fluxes for a given set of parameters.
2. To decompose the combined fluxes of kinetic model to the isotopeexchange fluxes, which differently affect isotopomer distribution.
3. To simulate the distribution of isotopomers using the total metabolite concentrations and decomposed fluxes obtained in steps 1 and 2.
Each simulation, performed through the steps 13, gives the distribution of isotopomers. The computed distribution is compared with the measured one using χ^{2 }criterion (see below) and a procedure of optimization is applied, which changes parameters and performs calculations each time passing through steps 13 with the objective to decrease χ^{2}. The steps 13 and the procedure of optimization are described next.
Step 1
The classical kinetic model is represented by a system of ordinary differential equations (ODE) simulating the evolution of concentrations of metabolites:
here c = {c_{1}, c_{2},..., c_{m}} is the vector of m metabolite concentrations described by the model. The metabolites considered in this particular case are presented in Figure 2. This scheme reflects the simplifications accepted for the model, in particular, pentoses represented by r5p and xu5p, which are in fast equilibrium. Oxidative branch of pentose phosphate pathway is represented by one reaction converting g6p into pentose phosphates. One lumped reaction connects citrate with malate, g3p goes directly to pep. We described in detail the most important reactions where carbon skeleton changes and label is redistributed, and combined the reactions which either do not change the carbon skeleton or only release CO_{2}.
N = {n_{ij}}, i = 1,..., m, j = 1,..., k is stoichiometric matrix for m metabolites and k reactions and v(c,p) = {v_{1}(c,p), v_{2}(c,p),...,v_{k}(c,p)}, is the vector of k reaction rates (metabolic fluxes), which are functions of concentrations and parameters such as Vmax. Concentrations vary with time in accordance with (1), and parameters are constant in each simulation, but they are unknown and could be found by fitting the experimental data. The rates for the most of enzymatic reactions are described by Michaelis' equation, although the reactions performing splitting/reformation of carbon skeleton of substrate molecule are considered in more detail taking into account the known elementary steps of reaction mechanism, as exemplified below for aldolase reaction. Matrix N is organized in such a way that production (input) and consumption (output) for the metabolites correspond to the arrows shown in Figure 2. Additional file 1 describes the whole set of differential equations of the model, which included the reactions of glycolysis, the TCA cycle, anaplerotic reactions, pentose phosphate pathway (PPP), exchange with extracellular glucose, lactate, and glutamate, and also the biosynthetic fluxes to nucleic acids and fatty acids. By solving these equations numerically, the program computes evolution of total metabolite concentrations and metabolic fluxes for a given set of parameters (Vmax, Km, rate constants).
Step 2
According to (1) differential equation for the concentration of a metabolite s (c_{s}) is:
The reactions j, which change concentrations c_{s}, change also the concentrations of isotopomers x_{s }= {x_{s1}, x_{s2},...,x_{sI}}. However, the net reaction rate v_{i }accounted for kinetic model (1) could be composed of several isotope exchange processes, which differently affect isotopomer distribution and have to be accounted separately. Below we describe an example of such decomposition for aldolase reaction. Thus, reactions v_{i }should be decomposed to several reactions u_{i}: v_{i }→ u_{i1 }, u_{i2 }.... This decomposition depends on reaction mechanism and is specific for each particular reaction. For instance, if a reaction does not produce any change in carbon skeleton of substrates, the decomposition implies only that the rates of transformation of substrate (u_{f}) into product and reverse transformation of product into substrate (u_{r}) must be calculated separately. The change of isotopic composition of reactants depends not only on net reaction rate (v_{i }= u_{f } u_{r}), but on the forward and reverse rates taken separately. If a reaction performs splitting/reformation of carbon skeleton of substrate molecule, additional isotopeexchange fluxes, different from forward and reverse reaction rates, could take place. Specifically, for aldolase reaction: fbp ↔ g3p + dhap, which proceeds through several elementary steps:
1. e + fbp ↔ efbp
2. efbp ↔ edhap + g3p
3. edhap ↔ e + dhap,
the fluxes through the whole reaction cycle (reactions 13) in the forward and reverse directions (u_{f}, u_{r}) and also the exchange flux of a half of fbp molecule with g3p (fbp ↔ g3p, steps 1 and 2 without entering to step 3, u_{e}) should be evaluated. The latter does not change the total metabolite concentrations but it characterizes the interchange of first three atoms of fbp with g3p pool. u_{f }and u_{r }change also total concentrations, their values indicate the transition of last three atoms of f6p to and from dhap. The algorithms for the calculations of such isotope exchange fluxes for given parameters and metabolite concentrations are described elsewhere (Selivanov et al, 2005).
For transaldolase reaction, f6p+e4p ↔ s7p+g3p, in addition to the forward and reverse fluxes through the whole reaction cycle two additional isotope exchange fluxes should be evaluated: f6p ↔ g3p and s7p ↔ e4p.
For transketolase reactions (x5p+r5p ↔ s7p+g3p, f6p+g3p ↔ x5p+e4p, f6p+r5p ↔ s7p+e4p) in addition to the six forward and reverse fluxes through these three reaction cycles the isotope exchanges x5p ↔ g3p, f6p ↔ e4p and s7p ↔ r5p should be evaluated.
All the isotope exchange fluxes necessary for the calculation of isotopomer dynamics could be evaluated if the total metabolite concentrations are known after the execution of kinetic model. When the kinetic simulation is done, Isodyn calculates the necessary isotope exchange fluxes in one step, as explained in more details in (Selivanov et al, 2005) and Additional file 1.
Step 3
When all the preparations are done (i.e. total metabolite concentrations and necessary fluxes are computed), the system of equations of type (3) for all isotopomers could be solved. Isodyn calls a module, which performs the computation of isotopomer dynamics.
The differential equations for isotopomers of each of m considered metabolites after the decomposition of fluxes could be presented in the form similar to (2):
here w_{j }is individual reaction rate that changes the concentration of isotopomer x_{sj}, which depends on the decomposition u_{jh }of rate v_{j}, vector of total concentrations c and isotopomer concentrations x.
Equation (3) describes the evolution of concentration of isotopomer i of substance s. To simulate the isotopomers distribution in metabolites, which are described by kinetic model (1), the equations of type (3) must be written for all the isotopomers of considered metabolites. This procedure is automated in our software.
To calculate the derivatives of isotopomers, the functions, which are specific for each reaction type, simulate the reaction mechanism in order to define the products for each isotopomersubstrate and also calculate the reaction rate for each given isotopomer transformation.
The simulation of reaction mechanism is based on operations with binary numbers. Since only two carbon isotopes are considered, any combination of carbons in the skeleton of a molecule could be reflected by a respective binary number, where "1" states for ^{13}C and "0" states for ^{12}C. In this way all hexose isotopomers could be represented by numbers from 000000 (decimal 0) to 111111 (binary equivalent of decimal 63), triose isotopomers range from 000 to 111 (decimal 7), etc.
When Isodyn simulates reactions for isotopomers it splits/recombines the binary representation of isotopomers in the same way as enzymes split and reform molecules. For instance since aldolase splits hexose (fbp) into two trioses, the respective Isodyn function splits binary numbers of fbp representatives, taking them one by one subsequently, e.g. when it takes "010101", it splits it to "010" and "101". The number 010101 (binary equivalent of decimal 21) indicates the position of the value of this isotopomer concentration (C_{010101 }with binary index or C_{21 }with decimal index) in the array of concentrations of fbp isotopomers. Knowing thus the position (and hence the concentration) of given isotopomer the algorithm calculates the reaction rate for a each isotopomer in a given reaction based on the total reaction rate (u^{ald}_{t}) and total concentration (C^{fbp}_{t}) known from the kinetic model simulation:
this rate is subtracted from the derivative of 21^{st }isotopomer of fbp and added to the derivatives of trioses "010" (decimal 2) and "101" (decimal 5). The arrays of derivatives are organized in the same way as those for concentrations.
The same isotopomers could participate in various reactions. Isodyn simulates all of them adding the reaction rates (with correct sign) to the respective derivatives. The functions performing such simulations are described in Additional file 1. They constitute a library, which can be used selectively.
In general, large system for isotopomers (3) depends on and could be solved simultaneously with (1). However, if the dynamics of isotopomer distribution is simulated in the conditions of metabolic steady state, the procedure of numerical solution could be simplified so that the general kinetic equations (1) could be solved separately from the solution for isotopomers (3). This case is presented here as the steps 13.
Here is a small example of application of above methodology, where we organize the calculation of derivatives for all the isotopomers of system consisting of substrate and product of pdh reaction with constant pyruvate input: v_{0 }→ pyr → accoa →. Let it is at metabolic steady state with initial isotopomer distribution for pyruvate (C_{000}, C_{001}, C_{010}, C_{011}, C_{100}, C_{101}, C_{110}, C_{111}) and accoa (C_{00}, C_{01}, C_{10}, C_{11}). Let v_{0 }is constant input of uniformely labeled pyruvate (111). In this system at steady state all rates are v_{0}, and let the computed total concentrations for some given set of parameters are C_{pyr }and C_{accoa}. Simulating this process Isodyn calls three functions that simulate respectively three reactions of the system. First function simulates constant input, it simply gives value v_{0 }to the derivative of uniformely labeled pyruvate, not touching other derivatives (D_{111 }= v_{0}). Then, the function, which simulates pdh takes the arrays for pyr and accoa, calculates what accoa isotopomer is produced from each substrate simulating decarboxylation by removing the first digit from binary representation of pyr, calculates the rates for each isotopomer transition, and adds this rate to the value of respective derivative as it is demonstrated in Table 4.
Table 4. The algorithm for the automated calculation of reaction rates (column "rate") and time derivatives for all the isotopomers (shown in left column) of pyruvate (d_{pyr }/dt) and accoa (d_{accoa }/dt) exemplified for pdh reaction (pyr → accoa).
Then Isodyn calls a function that simulates efflux of accoa as it is demonstrated in Table 5.
Table 5. The algorithm for the automated calculation of reaction rates (column "rate") and time derivatives for all the isotopomers (shown in left column) of accoa (d_{accoa }/dt) exemplified for the reaction of efflux of accoa (accoa →).
After the simulation in such a way of all the reactions of considered pathway, the whole array of derivatives (3) for all isotopomers at a given time point is formed. The function that calculates derivatives as described above could be called by any ODE solver, which solves the ODE system thus constructed. Isodyn implements several methods for ODE solving provided for C++ by Press et al [30], including fourthorder RungeKutta, BulirschStoer and Rosenbrock method for stiff systems. Also is implemented implicit RungeKutta 5^{th }order method for stiff systems (Radau5), described in [31] and backward differentiation formulas as their implemented in the solver DASSL [32] written in Fortran but linked with the C++ code of Isodyn. The implementation of various methods allows to find the fastest solution for a given problem. In the problem considered here the combination of DASSL for kinetic model with RungeKutta method as it is implemented in [30] for isotopomer distribution provided the fastest way to obtain solution.
The fact that solving the equations of kinetic model is followed by solving the system for isotopomers gives a way of checking solutions. Isodyn checks if the sum of isotopomers obtained after solving large isotopomer system is the same as total concentrations of metabolites obtained by kinetic model. The large isotopomer system could easily became stiff so that the employed method of solution does not provide necessary accuracy. In this case the program indicates the inaccuracy of solution.
Optimization
Simulation with an initial set of kinetic parameters gives a set of fluxes and corresponding isotopomer distribution, which could be compared with the measured distribution. Mass isotopomers of lactate and ribose were measured; to compare them with model prediction, the sums corresponding to respective measured mass isotopomers were calculated and the fractions with respect to the total amount were found. The difference between experimental data and the prediction was characterized by normalized square deviations (χ^{2}=Σ_{i}((f^{e}_{i}f^{t}_{i})/σ^{e}_{i})^{2 }), where f^{e}_{i }is experimental mass isotopomer fraction, f^{t}_{i }is the predicted one, σ^{e}_{i }is experimental standard deviation. The objective was to find the set of parameters (Vmax for simple reactions described by Michaelis' equation with fixed Km, and rate constants for elementary steps of more complex reactions like transketoliase and transaldolase, described elsewhere [9], in total 29 parameters), which minimize χ^{2}. To this end we subsequently used modified Simulated Annealing algorithm [33] and genetic algorithm. Random change of parameter values in Simulated Annealing we combined with Powell's method of coordinate descent [30], modifying each parameter in the direction that provides decrease of χ^{2}. After each stochastic perturbation of parameters and descent to a local minimum of χ^{2}, Isodyn saved the set of parameters and respective fluxes corresponding to the local minimum. These sets were used as an initial population for genetic algorithm. It performed crossover of these sets, "mutations" of randomly selected parameters and selection of obtained sets using the same criterion of χ^{2}. After the finding the global minimum of χ^{2 }a range of fluxes corresponding to a specific level of confidence of χ^{2}, was taken as a confidence interval corresponding to the chosen level of confidence [30].
Experimental methods
Cell culture conditions
Jurkat (acute T cell leukaemia), obtained form the American Type Culture Collection, were grown in RPMI 1640 culture medium (SigmaAldrich Co, St Louis, MO) supplemented with 10% heatinactivated FCS (PAA Laboratories, Pasching, Austria), 2 mM Lglutamine, 100 U × mL^{1 }of penicillin and 100 μg × mL^{1 }of streptomycin (Invitrogen, Paisley, UK) at 37°C in a humidified atmosphere of 5% CO_{2}. At the time of treatment with edelfosine for isotopomeric distribution analysis, cells (250000 cells × mL^{1}) were incubated in 75cm^{2 }Petri dishes in the RPMI 1640 medium, and with 10 mmol × L^{1 }of [1,2^{13}C_{2}]glucose (SigmaAldrich Co, St Louis, MO) at 50% isotope enrichment for 48 hours. At the end of the incubations, cells were centrifuged (1,500 rpm for 5 minutes) and medium was collected for glucose and lactate analysis, whereas cell pellets were frozen for RNA ribose analysis. Cells were counted with a haemocytometer, and edelfosineinduced apoptosis was assessed by flow cytometry using a fluorescenceactivated cell sorter (FACS), as the percentage of cells in the subG0 region (hypodiploidy) in cell cycle analysis.
Glucose and lactate concentration
From culture medium, glucose and lactate concentration were determined as previously described [34,35] using a Cobas Mira Plus chemistry analyzer (HORIBA ABX, Montpellier, France) at the beginning and at the end of the treatments.
The calculation of fluxes of glucose consumption and lactate production normalized per intracellular volume
The glucose (C_{g}) consumption rate (v) normalized per number of cells can be defined as follows. Glucose consumption rate for the given number of cells n_{t }is
Assuming exponential cell growth n_{t }could be expressed as:
Separation of variables in (5) with substitution (6) gives:
Integration of (7) gives
k could be defined from (6):
, and substitution in (8) gives:
During the 48 hours of incubation without edelfosine (control) the number of cells in average increased from 0.25 × 10^{6 }cells × mL^{1 }to 0.36 × 10^{6 }cells × mL^{1}. Glucose concentration decreased in average from 8.66 mM to 4.76 mM.
This consumption of glucose induces changes in intracellular concentrations of metabolites. To compare this flux with intracellular metabolic fluxes the units were changed to mM/min in the intracellular volume by dividing the above value of v to the volume of 10^{6 }cells, which was 0.05 mL:
This value characterizes the change of intracellular concentration per minute and all the fluxes were expressed in the same units.
Lactate isotopomer distribution
Lactate from the cell culture medium was extracted by ethyl acetate after acidification with HCl. Lactate was derivatized to its propylamideheptafluorobutyric form and the m/z 328 (carbons 13 of lactate, chemical ionization) was monitored as described [36].
RNA Ribose isotopomer distribution
RNA ribose was isolated by acid hydrolysis of cellular RNA after Trizol (Invitrogen) purification of cell extracts. Ribose isolated from RNA was derivatized to its aldonitrile acetate form using hydroxylamine in pyridine and acetic anhydride. We monitored the ion cluster around the m/z 256 (carbons 15 of ribose, chemical ionization) to find the molar enrichment and distribution of ^{13}C labels in ribose^{43}.
Gas Chromatography/Mass Spectrometry
Mass spectral data were obtained on the HP5973 mass selective detector connected to an HP6890 gas chromatograph. The settings are as follows: GC inlet 230°C, transfer line 280°C, MS source 230°C, MS quad 150°C. An HP5 capillary column (30 m length, 250 μm diameter, 0.25 μm film thickness) was used for analysis of ribose and lactate.
In vitro experiments were carried out using duplicate cultures each time for each treatment regimen. Mass spectral analyses were carried out by three independent automated injections of 1 μl of each sample and were accepted only if the standard sample deviation was less than 1% of the normalized peak intensity.
ROS production
ROS production was monitored using the fluorescente probe 2',7'dichlordihydrofluorescein diacetate (H_{2}DCFDA) (Invitrogen, Carlsbad, CA). Jurkat cells (150,000 cells/well) were grown in RPMI 1640 medium (as described above) in 6well culture plates. Prior to Edelfosine treatment, cells were harvested by centrifugation and preloaded with 10 μM H_{2}DCFDA in PBS for 30 min at 37°C, followed by wash in PBS to remove unloaded probe, addition of fresh medium containing 02 μg × mL^{1 }edelfosine, and incubation at 37°C/5% CO_{2 }for 48 hr. After treatment, cells were harvested, washed in PBS twice, and resuspended in PBS before DCF fluorescence was read with excitation and emission wavelengths at 495 nm and 527 nm, respectively. All fuorescence readings were normalizad to cell counts. The same treatment was also performed for cells grown in cover slips and DCF fluorescence examined using a BX51 fluorescence microscope (Olympus, Melvilla, NY).
List of abbreviations
Metabolites: glc: glucose; g6p, glucose6phosphate; f6p: fructose6phosphate; g3p: glyceraldehydes3phosphate; dhap: dihydroxyacetone phosphate; pep: phosphoenolpyruvate; pyr: pyruvate; lac: lactate; xu5p: xylulose5phosphate; r5p: ribose5phosphate; s7p: sedoheptulose7phosphate; e4p: erythrose4phosphate; accoa: acetylCoA; cit: citrate; kg: αketoglutarate; glut: glutamate; oaa: oxaloacetate. Enzyme reactions: hk: hexokinase; pfk: phosphofructokinase; fbpase: fructose bisphosphatase; pk: pyruvate kinase; pepck: phosphoenolpyruvate caboxykinase; pdh: pyruvate dehydrogenase complex; pc: pyruvate carboxylase; CitSyn: citrate synthase.
Authors' contributions
VAS performed the theoretical analysis of data and wrote the paper, PV performed experiments and wrote the paper, FM analyzed data, TWMF analysed data and performed experiments, WNPL analysed data, MC analysed data and wrote the paper. All authors read and approved the final manuscript.
Acknowledgements
This study was supported by Spanish Government: Ministerio de Ciencia e Innovación (SAF200800164, SAF200504293 and SAF200802251); from Red Temática de Investigación Cooperativa en Cáncer, Instituto de Salud Carlos III, Spanish Ministry of Science and Innovation & European Regional Development Fund (ERDF) "Una manera de hacer Europa" (ISCIIIRTICC grants RD06/0020/0046 and RD06/0020/1037); government of Catalonia (2009SGR1308) and European Commission (FP6, BioBridge LSHGCT2006037939, FP7, Diaprepp HealthF22008202013, Etherpath KBBEgrant agreement 222639). Additional support from the US National Institute of Health (grant # 1R01CA10119901 and 1R01CA11843401A2), National Science Foundation EPSCoR (grant # EPS0447479) is acknowledged. Mass spectrometry facility was supported by grants to WNP Lee from UCLA Center of Excellence (PO1 AT00396001) and from HarborUCLA GCRC (MO1 RR0042533).
References

Vizan P, SanchezTena S, AlcarrazVizan G, Soler M, Messeguer R, Pujol MD, Lee WP, Cascante M: Characterization of the metabolic changes underlying growth factor angiogenic activation: identification of new potential therapeutic targets.
Carcinogenesis 2009, 30:946952. PubMed Abstract  Publisher Full Text

Vizan P, Boros LG, Figueras A, Capella G, Mangues R, Bassilian S, Lim S, Lee WP, Cascante M: Kras codonspecific mutations produce distinctive metabolic phenotypes in NIH3T3 mice [corrected] fibroblasts.
Cancer Res 2005, 65:55125515. PubMed Abstract  Publisher Full Text

Marin S, Lee WP, Bassilian S, Lim S, Boros LG, Centelles JJ, FernAndezNovell JM, Guinovart JJ, Cascante M: Dynamic profiling of the glucose metabolic network in fasted rat hepatocytes using [1,213C2]glucose.
Biochem J 2004, 381:287294. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Boros LG, Lapis K, Szende B, TomoskoziFarkas R, Balogh A, Boren J, Marin S, Cascante M, Hidvegi M: Wheat germ extract decreases glucose uptake and RNA ribose formation but increases fatty acid synthesis in MIA pancreatic adenocarcinoma cells.
Pancreas 2001, 23:141147. PubMed Abstract  Publisher Full Text

Boros LG, Cascante M, Lee WNP: Metabolic profiling of cell growth and death in cancer: applications in drug discovery.
Drug Discov Today 2002, 7:364372. PubMed Abstract  Publisher Full Text

Boren J, Cascante M, Marin S, CominAnduix B, Centelles JJ, Lim S, Bassilian S, Ahmed S, Lee WN, Boros LG: Gleevec (STI571) influences metabolic enzyme activities and glucose carbon flow toward nucleic acid and fatty acid synthesis in myeloid tumor cells.
J Biol Chem 2001, 276:3774737753. PubMed Abstract  Publisher Full Text

Wiechert W: 13C metabolic flux analysis.
Metab Eng 2001, 3:19520.
Review. PubMed PMID: 11461141
PubMed Abstract  Publisher Full Text 
Selivanov VA, Puigjaner J, Sillero A, Centelles JJ, RamosMontoya A, Lee PW, Cascante M: An optimized algorithm for flux estimation from isotopomer distribution in glucose metabolites.
Bioinformatics 2004, 20:33873397. PubMed Abstract  Publisher Full Text

Selivanov VA, Meshalkina LE, Solovjeva ON, Kuchel PW, RamosMontoya A, Kochetov GA, Lee PWN, Cascante M: Rapid simulation and analysis of isotopomer distributions using constraints based on enzyme mechanisms: an example from HT29 cancer cells.
Bioinformatics 2005, 21:35583564. PubMed Abstract  Publisher Full Text

Selivanov VA, Marin S, Lee PWN, Cascante M: Software for dynamic analysis of tracerbased metabolomic data: estimation of metabolic fluxes and their statistical analysis.
Bioinformatics 2006, 22:28062812. PubMed Abstract  Publisher Full Text

Wahl SA, Noh K, Wiechert W: 13C labeling experiments at metabolic nonstationary conditions: an exploratory study.
BMC Bioinformatics 2008, 9:152. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Hanahan D, Weinberg RA: The hallmarks of cancer.
Cell 2000, 100:5770. PubMed Abstract  Publisher Full Text

Johnstone RW, Ruefli AA, Lowe SW: Apoptosis: a link between cancer genetics and chemotherapy.
Cell 2002, 108:153164. PubMed Abstract  Publisher Full Text

Ferreira CG, Epping M, Kruyt FAE, Giaccone G: Apoptosis: target of cancer therapy.
Clin Cancer Res 2002, 8:20242034. PubMed Abstract  Publisher Full Text

Vrablic AS, Albright CD, Craciunescu CN, Salganik RI, Zeisel SH: Altered mitochondrial function and overgeneration of reactive oxygen species precede the induction of apoptosis by 1Ooctadecyl2methylracglycero3phosphocholine in p53defective hepatocytes.
FASEB J 2001, 15:17391744. PubMed Abstract  Publisher Full Text

Mates JM, SanchezJimenez FM: Role of reactive oxygen species in apoptosis: implications for cancer therapy.
Int J Biochem Cell Biol 2000, 32:157170. PubMed Abstract  Publisher Full Text

Engel RH, Evens AM: Oxidative stress and apoptosis: a new treatment paradigm in cancer.
Front Biosci 2006, 11:30012. PubMed Abstract  Publisher Full Text

Xu R, Pelicano H, Zhang H, Giles FJ, Keating MJ, Huang P: Synergistic effect of targeting mTOR by rapamycin and depleting ATP by inhibition of glycolysis in lymphoma and leukemia cells.
Leukemia 2005, 19:21532158. PubMed Abstract  Publisher Full Text

Gajate C, An F, Mollinedo F: Differential cytostatic and apoptotic effects of ecteinascidin743 in cancer cells. Transcriptiondependent cell cycle arrest and transcriptionindependent JNK and mitochondrial mediated apoptosis.
J Biol Chem 2002, 277:4158041589. PubMed Abstract  Publisher Full Text

Gajate C, Fonteriz RI, Cabaner C, AlvarezNoves G, AlvarezRodriguez Y, Modolell M, Mollinedo F: Intracellular triggering of Fas, independently of FasL, as a new mechanism of antitumor ether lipidinduced apoptosis.
Int J Cancer 2000, 85:674682. PubMed Abstract  Publisher Full Text

Gajate C, Del CantoJanez E, Acuna AU, AmatGuerri F, Geijo E, SantosBeneit AM, Veldman RJ, Mollinedo F: Intracellular triggering of Fas aggregation and recruitment of apoptotic molecules into Fasenriched rafts in selective tumor cell apoptosis.
J Exp Med 2004, 200:353365. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Mollinedo F, MartinezDalmau R, Modolell M: Early and selective induction of apoptosis in human leukemic cells by the alkyllysophospholipid ET18OCH3.
Biochem Biophys Res Commun 1993, 192:603609. PubMed Abstract  Publisher Full Text

Mollinedo F, FernandezLuna JL, Gajate C, MartinMartin B, Benito A, MartinezDalmau R, Modolell M: Selective induction of apoptosis in cancer cells by the ether lipid ET18OCH3 (Edelfosine): molecular structure requirements, cellular uptake, and protection by Bcl2 and BclX(L).
Cancer Res 1997, 57:13201328. PubMed Abstract  Publisher Full Text

Gajate C, SantosBeneit A, Modolell M, Mollinedo F: Involvement of cJun NH2terminal kinase activation and cJun in the induction of apoptosis by the ether phospholipid 1Ooctadecyl2Omethylracglycero3phosphocholine.
Mol Pharmacol 1998, 53:602612. PubMed Abstract  Publisher Full Text

Gajate C, Mollinedo F: Edelfosine and perifosine induce selective apoptosis in multiple myeloma by recruitment of death receptors and downstream signaling molecules into lipid rafts.
Blood 2007, 109:711719. PubMed Abstract  Publisher Full Text

Na H, Surh Y: The antitumor ether lipid edelfosine (ET18OCH3) induces apoptosis in Hras transformed human breast epithelial cells: by blocking ERK1/2 and p38 mitogenactivated protein kinases as potential targets.
Asia Pac J Clin Nutr 2008, 17(Suppl 1):204207. PubMed Abstract

Zhang H, Gajate C, Yu L, Fang Y, Mollinedo F: Mitochondrialderived ROS in edelfosineinduced apoptosis in yeasts and tumor cells.
Acta Pharmacol Sin 2007, 28:888894. PubMed Abstract  Publisher Full Text

Selivanov VA, de Atauri P, Centelles JJ, Cadefau J, Parra J, Cussó R, Carreras J, Cascante M: The changes in the energy metabolism of human muscle induced by training.
J Theor Biol 2008, 252:402410. PubMed Abstract  Publisher Full Text

Selivanov VA, Votyakova TV, Zeak JA, Trucco M, Roca J, Cascante M: Bistability of mitochondrial respiration underlies paradoxical reactive oxygen species generation induced by anoxia.
PLoS Comput Biol 2009, e1000619. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Press W, Flannery B, Teukolsky S, Vetterling W: Numerical Recipes in C: The Art of Scientific Computing. Cambridge University Press, New York, USA; 2002.

Hairer E, Wanner G: Solving Ordinary Differential Equations II. Stiff and DifferentialAlgebraic Problems.
Springer Series in Computational Mathematics 14, SpringerVerlag 1996.

Petzold LR: A description of DASSL: A differential/algebraic system solver.

Kirkpatrick S, Gelatt CD, Vecchi MP: Optimization by Simulated Annealing.
Science 1983, 220:671680. PubMed Abstract  Publisher Full Text

Kunst A, Draeger B, Ziegenhorn J: DGlucose; UVmethods with hexokinase and glucose6phosphate dehydrogenase. In Methods of Enzymatic Analysis. Verlag Chemie, Weinheim, Germany; 1984.

Passonneau JV, Lowry OH: Enzymatic analysis: a practical guide. The Humana Press Inc, Totowa, New Jersey, USA; 1993.

Lee WN, Boros LG, Puigjaner J, Bassilian S, Lim S, Cascante M: Mass isotopomer study of the nonoxidative pathways of the pentose cycle with [1,213C2]glucose.
Am J Physiol 1998, 274:E843E851. PubMed Abstract  Publisher Full Text