Edelfosineinduced metabolic changes in cancer cells that precede the overproduction of reactive oxygen species and apoptosis
 Vitaly A Selivanov†^{1, 2},
 Pedro Vizán†^{1, 7},
 Faustino Mollinedo^{3},
 Teresa WM Fan^{4, 5},
 Paul WN Lee^{6} and
 Marta Cascante^{1}Email author
DOI: 10.1186/175205094135
© Selivanov et al; licensee BioMed Central Ltd. 2010
Received: 18 January 2010
Accepted: 6 October 2010
Published: 6 October 2010
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 [1–6]. 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") [8–10] 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 [15–17] 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 [19–25]. 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 metabolic fluxes of glucose consumption (J_{Glc}) and lactate production (J_{lac}).
J_{Glc} (mM/min)  J_{Lac} (mM/min)  

Con  0.09  0.1225 
e0.5  0.1196  0.1565 
e1  0.1359  0.1279 
Measured isotopomer distributions
Isotopomer distribution in lactate secreted into the medium and RNA ribose.
Lactate:  control  sd  fit  0.5 μg/mL  sd  fit  1 μg/mL  sd  fit 

m0  0.7900  0.003  0.7900  0.7860  0.001  0.7860  0.7910  0.003  0.7920 
m1  0.0097  0.001  0.0096  0.0099  0.000  0.0099  0.0098  0.001  0.0098 
m2  0.1990  0.003  0.1990  0.2030  0.001  0.2030  0.1980  0.003  0.1970 
m3  0.0013  0.001  0.0013  0.0013  0.001  0.0013  0.0010  0.002  0.0011 
Rib5P  
m0  0.5480  0.011  0.5480  0.5570  0.003  0.5570  0.5790  0.011  0.5790 
m1  0.2310  0.004  0.2310  0.2210  0.007  0.2250  0.2130  0.004  0.2110 
m2  0.1490  0.005  0.1490  0.1410  0.003  0.1430  0.1340  0.005  0.1390 
m3  0.0423  0.001  0.0420  0.0483  0.005  0.0449  0.0426  0.001  0.0418 
m4  0.0299  0.001  0.0296  0.0325  0.002  0.0293  0.0312  0.001  0.0291 
dilution:  0.3450  0.2080  0.2920  
χ  0.1760  4.7400  5.7500 
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 [8–10] and Additional file 1.
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
Simulated metabolic fluxes corresponding to the best fit presented in Table 2.
flux  Control  Edelf: 0.5 μg/mL  Edelf:1 μg/mL  

hk***  0.0901118  1  0.12  1  0.135  1 
pfk  0.0968532  1.07481  0.120985  1.00821  0.135572  1.004237 
fbpase  0.00716729  0.07954  1.46E3  0.01219  0.0009715  0.007197 
pep > g3p  0.00694225  0.07704  0.00749499  0.06246  0.0025560  0.018933 
g3p > pep  0.186074  2.06492  0.246177  2.05148  0.271442  2.010681 
pk  0.181224  2.01110  0.253755  2.11463  0.28068  2.079111 
TCA > pyr  2.23E3  0.02470  0.0152166  0.12681  0.0118396  0.087701 
lac >  0.12219  1.35598  0.153661  1.28051  0.128901  0.954822 
pdh***  0.057984  0.64347  0.0815881  0.67990  0.140089  1.037696 
pc  1.05E3  0.01165  0.0185053  0.15421  0.0116895  0.086589 
CitSyn**  0.0351159  0.38969  0.0547404  0.45617  0.120038  0.889170 
cit > mal**  0.0366362  0.40656  0.0549581  0.45798  0.120969  0.896067 
PPP*  1.32E4  1.47E3  2.80E6  2.33E5  3.93E6  2.91E5 
r5p >  4.38E4  0.00486  6.84E4  0.00570  0.0005799  0.004296 
cit > glt  3.50E3  0.03884  2.13E3  0.01775  0.0014090  0.010437 
 > cit  0.00549634  0.06099  2.35E3  0.01959  0.0023398  0.017332 
Transketolase:  
xu5p > s7p***  0.00826654  0.09174  2.26E7  1.88E6  1.46E6  1.08E5 
s7p > xu5p  0.0084719  0.09402  4.24E4  0.00353  0.0003579  0.002651 
f6p > xu5p  5.66E4  0.00628  3.54E5  0.00030  3.03E5  0.000224 
xu5p > f6p***  5.38E4  0.00597  3.45E9  2.87E8  1.89E8  1.40E7 
f6p > s7p  0.00534401  0.05930  2.46E4  0.00205  0.0001979  0.001466 
s7p > f6p  0.0052113  0.05783  4.49E5  0.00037  3.08E5  0.000228 
xu5p <  > g3p***  8.75E4  9.71E3  3.26E8  2.71E7  2.21E7  1.64E6 
f6p <  > e4p  3.48E4  0.00386  3.75E6  3.12E5  2.63E6  1.95E5 
r5p <  > s7p  0.0800302  0.88812  2.94E3  0.02453  0.0023626  0.017501 
Transaldolase:  
f6p > s7p  1.54E4  0.00171  2.57E4  0.00214  0.0002003  0.001484 
s7p > f6p  2.48E5  0.00028  2.25E5  0.00019  4.54E6  3.36E5 
f6p <  > g3p  2.64E7  2.93E6  1.09E6  9.05E6  2.51E7  1.86E6 
s7p <  > e4p  0.0144573  0.16044  0.0053368  0.04447  0.0037384  0.027691 
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.
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
 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
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
 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.
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.
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.
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).
reaction  rate  d_{pyr}/dt  d_{accoa}/dt 

000 → 00  u_{000}= v_{0}C_{000}/C_{pyr}  D_{000}= u_{000}  D_{00}= u_{000} 
001 → 01  u_{001}= v_{0}C_{001}/C_{pyr}  D_{001}= u_{001}  D_{01}= u_{001} 
010 → 10  u_{010}= v_{0}C_{010}/C_{pyr}  D_{010}= u_{010}  D_{10}= u_{010} 
011 → 11  u_{011}= v_{0}C_{011}/C_{pyr}  D_{011}= u_{011}  D_{11}= u_{011} 
100 → 00  u_{100}= v_{0}C_{100}/C_{pyr}  D_{100}= u_{100}  D_{00}= u_{000}+u_{100} 
101 → 01  u_{101}= v_{0}C_{101}/C_{pyr}  D_{101}= u_{101}  D_{01}= u_{101}+u_{001} 
110 → 10  u_{110}= v_{0}C_{110}/C_{pyr}  D_{110}= u_{110}  D_{10}= u_{110}+u_{010} 
111 → 11  u_{111}= v_{0}C_{111}/C_{pyr}  D_{111}= v_{0}u_{111}  D_{11}= u_{111}+u_{011} 
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 →).
reaction  rate  d_{accoa} /dt 

00 →  u_{00}= v_{0}C_{00}/C_{accoa}  D_{00}= u_{000}+u_{100}u_{00} 
01 →  u_{01}= v_{0}C_{01}/C_{accoa}  D_{01}= u_{101}+u_{001}u_{01} 
10 →  u_{10}= v_{0}C_{10}/C_{accoa}  D_{10}= u_{110}+u_{010}u_{10} 
11 →  u_{11}= v_{0}C_{11}/C_{accoa}  D_{11}= u_{111}+u_{011}u_{11} 
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
k could be defined from (6):
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).
Notes
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.
Declarations
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).
Authors’ Affiliations
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. 10.1093/carcin/bgp083View ArticlePubMedGoogle Scholar
 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. 10.1158/00085472.CAN050074View ArticlePubMedGoogle Scholar
 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. 10.1042/BJ20031737PubMed CentralView ArticlePubMedGoogle Scholar
 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. 10.1097/0000667620010800000004View ArticlePubMedGoogle Scholar
 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. 10.1016/S13596446(02)021797View ArticlePubMedGoogle Scholar
 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.PubMedGoogle Scholar
 Wiechert W: 13C metabolic flux analysis. Metab Eng. 2001, 3: 19520. Review. PubMed PMID: 11461141 10.1006/mben.2001.0187View ArticlePubMedGoogle Scholar
 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. 10.1093/bioinformatics/bth412View ArticlePubMedGoogle Scholar
 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. 10.1093/bioinformatics/bti573View ArticlePubMedGoogle Scholar
 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. 10.1093/bioinformatics/btl484View ArticlePubMedGoogle Scholar
 Wahl SA, Noh K, Wiechert W: 13C labeling experiments at metabolic nonstationary conditions: an exploratory study. BMC Bioinformatics. 2008, 9: 152 10.1186/147121059152PubMed CentralView ArticlePubMedGoogle Scholar
 Hanahan D, Weinberg RA: The hallmarks of cancer. Cell. 2000, 100: 5770. 10.1016/S00928674(00)816839View ArticlePubMedGoogle Scholar
 Johnstone RW, Ruefli AA, Lowe SW: Apoptosis: a link between cancer genetics and chemotherapy. Cell. 2002, 108: 153164. 10.1016/S00928674(02)006256View ArticlePubMedGoogle Scholar
 Ferreira CG, Epping M, Kruyt FAE, Giaccone G: Apoptosis: target of cancer therapy. Clin Cancer Res. 2002, 8: 20242034.PubMedGoogle Scholar
 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. 10.1096/fj.000300comView ArticlePubMedGoogle Scholar
 Mates JM, SanchezJimenez FM: Role of reactive oxygen species in apoptosis: implications for cancer therapy. Int J Biochem Cell Biol. 2000, 32: 157170. 10.1016/S13572725(99)000886View ArticlePubMedGoogle Scholar
 Engel RH, Evens AM: Oxidative stress and apoptosis: a new treatment paradigm in cancer. Front Biosci. 2006, 11: 30012. 10.2741/1798View ArticlePubMedGoogle Scholar
 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. 10.1038/sj.leu.2403968View ArticlePubMedGoogle Scholar
 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. 10.1074/jbc.M204644200View ArticlePubMedGoogle Scholar
 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. 10.1002/(SICI)10970215(20000301)85:5<674::AIDIJC13>3.0.CO;2ZView ArticlePubMedGoogle Scholar
 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. 10.1084/jem.20040213PubMed CentralView ArticlePubMedGoogle Scholar
 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. 10.1006/bbrc.1993.1458View ArticlePubMedGoogle Scholar
 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.PubMedGoogle Scholar
 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.PubMedGoogle Scholar
 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. 10.1182/blood200604016824View ArticlePubMedGoogle Scholar
 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.PubMedGoogle Scholar
 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. 10.1111/j.17457254.2007.00568.xView ArticlePubMedGoogle Scholar
 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. 10.1016/j.jtbi.2007.09.039View ArticlePubMedGoogle Scholar
 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, e1000619Google Scholar
 Press W, Flannery B, Teukolsky S, Vetterling W: Numerical Recipes in C: The Art of Scientific Computing. 2002, Cambridge University Press, New York, USAGoogle Scholar
 Hairer E, Wanner G: Solving Ordinary Differential Equations II. Stiff and DifferentialAlgebraic Problems. Springer Series in Computational Mathematics 14, SpringerVerlag. 1996Google Scholar
 Petzold LR: A description of DASSL: A differential/algebraic system solver. SAND828637, Sandia National Laboratories, Livermore. 1982Google Scholar
 Kirkpatrick S, Gelatt CD, Vecchi MP: Optimization by Simulated Annealing. Science. 1983, 220: 671680. 10.1126/science.220.4598.671View ArticlePubMedGoogle Scholar
 Kunst A, Draeger B, Ziegenhorn J: DGlucose; UVmethods with hexokinase and glucose6phosphate dehydrogenase. Methods of Enzymatic Analysis. 1984, Verlag Chemie, Weinheim, GermanyGoogle Scholar
 Passonneau JV, Lowry OH: Enzymatic analysis: a practical guide. 1993, The Humana Press Inc, Totowa, New Jersey, USAView ArticleGoogle Scholar
 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.PubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.