Home | About | Journals | Submit | Contact Us | Français |

**|**BMC Syst Biol**|**v.4; 2010**|**PMC2984393

Formats

Article sections

- Abstract
- Background
- Results and Discussion
- Conclusions
- Methods
- List of abbreviations
- Authors' contributions
- Supplementary Material
- References

Authors

Related links

BMC Syst Biol. 2010; 4: 135.

Published online 2010 October 6. doi: 10.1186/1752-0509-4-135

PMCID: PMC2984393

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}

Vitaly A Selivanov: ude.bu@vonaviles; Pedro Vizán: ku.gro.recnac@naziV.ordeP; Faustino Mollinedo: se.lasu@nillomf; Teresa WM Fan: ude.sivadcu@nafwt; Paul WN Lee: gro.demoibal@eel; Marta Cascante: ude.bu@etnacsacatram

Received 2010 January 18; Accepted 2010 October 6.

Copyright ©2010 Selivanov et al; licensee 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.

This article has been cited by other articles in PMC.

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.

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.

The application of Isodyn to the analysis of mechanism of edelfosine-induced apoptosis revealed primary drug-induced metabolic changes, which are important for the subsequent initiation of apoptotic program. Initiation of such metabolic changes could be exploited in anticancer therapy.

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 non-stationary 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 pre-steady state phase for all isotopomers. Of course, there is always a possibility of measuring the labeling of such stores and apply classical ^{13}C-MFA for the "fast" intermediates of central metabolism considering that they are in quasi-steady 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 non-stationary 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 non-stationary 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 (1-*O*-octadecyl-2-*O*-methyl*rac*-glycero-3-phosphocholine, ET-18-OCH_{3}) selectively induces apoptosis in cancer cells [19-25]. The cell-killing 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.

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 Table1.1. These values calculated directly from experimental data were used in the fitting of isotopic isomer distribution described below.

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 ribose-5-phosphate (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 edelfosine-treated cells are shown in Table Table22.

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 Table2.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 non-labeled 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 Table2,2, in edelfosine-treated cells dilution decreased, which indicates that a greater fraction of RNA was synthesized *de novo*.

Table Table33 show the fluxes corresponding to the best fit shown in Table Table22 and indicates the fluxes for which the difference between treated and non-treated 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 Table33 illustrates that a small change in the distribution of mass isotopomers shown in Table Table22 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 4-5% 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 edelfosine-treated cells. Subsequent fitting procedure decreases χ^{2 }to 5.75 (Table (Table2)2) just by increase of metabolic fluxes through the TCA cycle and changes in the pentose phosphate pathway as Table Table33 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 Figure1.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.

Another important distinction in the revealed flux profiles of control and edelfosine-treated Jurkat cells is the difference in the fluxes of pentose phosphate pathway. In control cells all pentose phosphate fluxes were higher than those in edelfosine-treated cells (Table (Table3).3). Moreover, the model predicted an increase in ribose consumption for RNA synthesis in edelfosine-treated cells, which led to a decrease in ribose 5-phosphate concentration, such that its synthesis became practically unidirectional. Thus, newly synthesized ribose 5-phosphate was consumed for RNA synthesis to a greater extent in edelfosine-treated cells than control (see flux r5p, Table Table3),3), in accordance with the value of ^{13}C dilution in r5p (see Table Table2).2). In contrast, in control cells, the newly synthesized ribose 5-phosphate is not massively used for RNA synthesis, but reutilized in the central energetic metabolism.

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.

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 isotope-exchange 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 1-3, 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 1-3 with the objective to decrease χ^{2}. The steps 1-3 and the procedure of optimization are described next.

The classical kinetic model is represented by a system of ordinary differential equations (ODE) simulating the evolution of concentrations of metabolites:

$$\text{d}c/\text{dt}=N\text{}v\left(c,p\right)$$

(1)

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 Figure2.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 Figure2.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).

According to (1) differential equation for the concentration of a metabolite s (c_{s}) is:

$${\text{dc}}_{\text{s}}/\text{dt}=\Sigma {\text{n}}_{\text{sj}}{\text{v}}_{\text{j}}\left(\text{c},\text{p}\right)$$

(2)

The reactions j, which change concentrations c_{s}, change also the concentrations of isotopomers **x _{s }**= {x

1. e + fbp ↔ e-fbp

2. e-fbp ↔ e-dhap + g3p

3. e-dhap ↔ e + dhap,

the fluxes through the whole reaction cycle (reactions 1-3) in the forward and reverse directions (*u _{f}*,

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.

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):

$${\text{dx}}_{\text{si}}/\text{dt}={\Sigma}_{\text{j}}{\Sigma}_{\text{h}}{\text{n}}_{\text{sj}}{\text{w}}_{\text{j}}\left({\text{u}}_{\text{jh}}\left(c,p\right),c,x\right)$$

(3)

here *w _{j }*is individual reaction rate that changes the concentration of isotopomer

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 isotopomer-substrate 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:

$${\text{u}}^{\text{ald}}{}_{21}={\text{u}}^{\text{ald}}{}_{\text{t}}{\text{*C}}_{\text{21}}{\text{/C}}^{\text{fbp}}{}_{\text{t}}$$

(4)

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 1-3.

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 Table44.

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 Table55.

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 fourth-order Runge-Kutta, Bulirsch-Stoer and Rosenbrock method for stiff systems. Also is implemented implicit Runge-Kutta 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 Runge-Kutta 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.

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].

Jurkat (acute T cell leukaemia), obtained form the American Type Culture Collection, were grown in RPMI 1640 culture medium (Sigma-Aldrich Co, St Louis, MO) supplemented with 10% heat-inactivated FCS (PAA Laboratories, Pasching, Austria), 2 mM L-glutamine, 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 (Sigma-Aldrich 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 edelfosine-induced apoptosis was assessed by flow cytometry using a fluorescence-activated cell sorter (FACS), as the percentage of cells in the sub-G0 region (hypodiploidy) in cell cycle analysis.

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 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

$$\frac{d{C}_{g}}{dt}=v\times {n}_{t}$$

(5)

Assuming exponential cell growth n_{t }could be expressed as:

$${n}_{t}={n}_{0}\times {e}^{k\times t}$$

(6)

Separation of variables in (5) with substitution (6) gives:

$$d{C}_{g}=\frac{v\times {n}_{0}}{k}\times {e}^{k\times t}d(k\times t)$$

(7)

Integration of (7) gives

$$\Delta {C}_{g}=\frac{v}{k}({n}_{t}-{n}_{0})\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}v=\frac{\Delta {C}_{g}\times k}{{n}_{t}-{n}_{0}}$$

(8)

k could be defined from (6):

$k=\frac{\mathrm{ln}({n}_{t}/{n}_{0})}{t}$, and substitution in (8) gives:

$$v=\frac{\Delta {C}_{g}\times \mathrm{ln}({n}_{t}/{n}_{0})}{({n}_{t}-{n}_{0})\times t}$$

(9)

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.

$$v=\frac{3.9\mu mol\times \mathrm{ln}(0.36/0.25)}{{10}^{6}\times (0.36-0.25)cells\times 2880\mathrm{min}}=0.0045\frac{\mu mol}{{10}^{6}cells\times \mathrm{min}}$$

(10)

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:

$$v=0.0045\frac{\mu mol}{{10}^{6}cells\times \mathrm{min}}=0.09\frac{mM}{\mathrm{min}}$$

(11)

This value characterizes the change of intracellular concentration per minute and all the fluxes were expressed in the same units.

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 1-3 of lactate, chemical ionization) was monitored as described [36].

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 hydroxyl-amine in pyridine and acetic anhydride. We monitored the ion cluster around the m/z 256 (carbons 1-5 of ribose, chemical ionization) to find the molar enrichment and distribution of ^{13}C labels in ribose^{43}.

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 HP-5 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 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 6-well 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 0-2 μ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).

Metabolites: glc: glucose; g6p, glucose-6-phosphate; f6p: fructose-6-phosphate; g3p: glyceraldehydes-3-phosphate; dhap: dihydroxyacetone phosphate; pep: phosphoenolpyruvate; pyr: pyruvate; lac: lactate; xu5p: xylulose-5-phosphate; r5p: ribose-5-phosphate; s7p: sedoheptulose-7-phosphate; e4p: erythrose-4-phosphate; accoa: acetyl-CoA; 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.

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.

**Supplementary materials**. PDF containing supplemental figures and other information.

Click here for file^{(478K, PDF)}

This study was supported by Spanish Government: Ministerio de Ciencia e Innovación (SAF2008-00164, SAF2005-04293 and SAF2008-02251); 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" (ISCIII-RTICC grants RD06/0020/0046 and RD06/0020/1037); government of Catalonia (2009-SGR1308) and European Commission (FP6, BioBridge LSHG-CT-2006-037939, FP7, Diaprepp Health-F2-2008-202013, Etherpath KBBE-grant agreement 222639). Additional support from the US National Institute of Health (grant # 1R01CA101199-01 and 1R01CA118434-01A2), National Science Foundation EPSCoR (grant # EPS-0447479) is acknowledged. Mass spectrometry facility was supported by grants to WNP Lee from UCLA Center of Excellence (PO1 AT003960-01) and from Harbor-UCLA GCRC (MO1 RR00425-33).

- Vizan P, Sanchez-Tena S, Alcarraz-Vizan 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:946–952. doi: 10.1093/carcin/bgp083. [PubMed] [Cross Ref]
- Vizan P, Boros LG, Figueras A, Capella G, Mangues R, Bassilian S, Lim S, Lee WP, Cascante M. K-ras codon-specific mutations produce distinctive metabolic phenotypes in NIH3T3 mice [corrected] fibroblasts. Cancer Res. 2005;65:5512–5515. doi: 10.1158/0008-5472.CAN-05-0074. [PubMed] [Cross Ref]
- Marin S, Lee WP, Bassilian S, Lim S, Boros LG, Centelles JJ, FernAndez-Novell JM, Guinovart JJ, Cascante M. Dynamic profiling of the glucose metabolic network in fasted rat hepatocytes using [1,2-13C2]glucose. Biochem J. 2004;381:287–294. doi: 10.1042/BJ20031737. [PubMed] [Cross Ref]
- Boros LG, Lapis K, Szende B, Tomoskozi-Farkas 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:141–147. doi: 10.1097/00006676-200108000-00004. [PubMed] [Cross Ref]
- Boros LG, Cascante M, Lee WNP. Metabolic profiling of cell growth and death in cancer: applications in drug discovery. Drug Discov Today. 2002;7:364–372. doi: 10.1016/S1359-6446(02)02179-7. [PubMed] [Cross Ref]
- Boren J, Cascante M, Marin S, Comin-Anduix 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:37747–37753. [PubMed]
- Wiechert W. 13C metabolic flux analysis. Metab Eng. 2001;3:195–20. doi: 10.1006/mben.2001.0187. Review. PubMed PMID: 11461141. [PubMed] [Cross Ref]
- Selivanov VA, Puigjaner J, Sillero A, Centelles JJ, Ramos-Montoya A, Lee PW, Cascante M. An optimized algorithm for flux estimation from isotopomer distribution in glucose metabolites. Bioinformatics. 2004;20:3387–3397. doi: 10.1093/bioinformatics/bth412. [PubMed] [Cross Ref]
- Selivanov VA, Meshalkina LE, Solovjeva ON, Kuchel PW, Ramos-Montoya 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:3558–3564. doi: 10.1093/bioinformatics/bti573. [PubMed] [Cross Ref]
- Selivanov VA, Marin S, Lee PWN, Cascante M. Software for dynamic analysis of tracer-based metabolomic data: estimation of metabolic fluxes and their statistical analysis. Bioinformatics. 2006;22:2806–2812. doi: 10.1093/bioinformatics/btl484. [PubMed] [Cross Ref]
- Wahl SA, Noh K, Wiechert W. 13C labeling experiments at metabolic nonstationary conditions: an exploratory study. BMC Bioinformatics. 2008;9:152. doi: 10.1186/1471-2105-9-152. [PMC free article] [PubMed] [Cross Ref]
- Hanahan D, Weinberg RA. The hallmarks of cancer. Cell. 2000;100:57–70. doi: 10.1016/S0092-8674(00)81683-9. [PubMed] [Cross Ref]
- Johnstone RW, Ruefli AA, Lowe SW. Apoptosis: a link between cancer genetics and chemotherapy. Cell. 2002;108:153–164. doi: 10.1016/S0092-8674(02)00625-6. [PubMed] [Cross Ref]
- Ferreira CG, Epping M, Kruyt FAE, Giaccone G. Apoptosis: target of cancer therapy. Clin Cancer Res. 2002;8:2024–2034. [PubMed]
- 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 1-O-octadecyl-2-methyl-rac-glycero-3-phosphocholine in p53-defective hepatocytes. FASEB J. 2001;15:1739–1744. doi: 10.1096/fj.00-0300com. [PubMed] [Cross Ref]
- Mates JM, Sanchez-Jimenez FM. Role of reactive oxygen species in apoptosis: implications for cancer therapy. Int J Biochem Cell Biol. 2000;32:157–170. doi: 10.1016/S1357-2725(99)00088-6. [PubMed] [Cross Ref]
- Engel RH, Evens AM. Oxidative stress and apoptosis: a new treatment paradigm in cancer. Front Biosci. 2006;11:300–12. doi: 10.2741/1798. [PubMed] [Cross Ref]
- 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:2153–2158. doi: 10.1038/sj.leu.2403968. [PubMed] [Cross Ref]
- Gajate C, An F, Mollinedo F. Differential cytostatic and apoptotic effects of ecteinascidin-743 in cancer cells. Transcription-dependent cell cycle arrest and transcription-independent JNK and mitochondrial mediated apoptosis. J Biol Chem. 2002;277:41580–41589. doi: 10.1074/jbc.M204644200. [PubMed] [Cross Ref]
- Gajate C, Fonteriz RI, Cabaner C, Alvarez-Noves G, Alvarez-Rodriguez Y, Modolell M, Mollinedo F. Intracellular triggering of Fas, independently of FasL, as a new mechanism of antitumor ether lipid-induced apoptosis. Int J Cancer. 2000;85:674–682. doi: 10.1002/(SICI)1097-0215(20000301)85:5<674::AID-IJC13>3.0.CO;2-Z. [PubMed] [Cross Ref]
- Gajate C, Del Canto-Janez E, Acuna AU, Amat-Guerri F, Geijo E, Santos-Beneit AM, Veldman RJ, Mollinedo F. Intracellular triggering of Fas aggregation and recruitment of apoptotic molecules into Fas-enriched rafts in selective tumor cell apoptosis. J Exp Med. 2004;200:353–365. doi: 10.1084/jem.20040213. [PMC free article] [PubMed] [Cross Ref]
- Mollinedo F, Martinez-Dalmau R, Modolell M. Early and selective induction of apoptosis in human leukemic cells by the alkyl-lysophospholipid ET-18-OCH3. Biochem Biophys Res Commun. 1993;192:603–609. doi: 10.1006/bbrc.1993.1458. [PubMed] [Cross Ref]
- Mollinedo F, Fernandez-Luna JL, Gajate C, Martin-Martin B, Benito A, Martinez-Dalmau R, Modolell M. Selective induction of apoptosis in cancer cells by the ether lipid ET-18-OCH3 (Edelfosine): molecular structure requirements, cellular uptake, and protection by Bcl-2 and Bcl-X(L) Cancer Res. 1997;57:1320–1328. [PubMed]
- Gajate C, Santos-Beneit A, Modolell M, Mollinedo F. Involvement of c-Jun NH2-terminal kinase activation and c-Jun in the induction of apoptosis by the ether phospholipid 1-O-octadecyl-2-O-methyl-rac-glycero-3-phosphocholine. Mol Pharmacol. 1998;53:602–612. [PubMed]
- 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:711–719. doi: 10.1182/blood-2006-04-016824. [PubMed] [Cross Ref]
- Na H, Surh Y. The antitumor ether lipid edelfosine (ET-18-O-CH3) induces apoptosis in H-ras transformed human breast epithelial cells: by blocking ERK1/2 and p38 mitogen-activated protein kinases as potential targets. Asia Pac J Clin Nutr. 2008;17(Suppl 1):204–207. [PubMed]
- Zhang H, Gajate C, Yu L, Fang Y, Mollinedo F. Mitochondrial-derived ROS in edelfosine-induced apoptosis in yeasts and tumor cells. Acta Pharmacol Sin. 2007;28:888–894. doi: 10.1111/j.1745-7254.2007.00568.x. [PubMed] [Cross Ref]
- 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:402–410. doi: 10.1016/j.jtbi.2007.09.039. [PubMed] [Cross Ref]
- 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. p. e1000619. [PMC free article] [PubMed] [Cross Ref]
- 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 Differential-Algebraic Problems. Springer Series in Computational Mathematics 14, Springer-Verlag. 1996.
- Petzold LR. A description of DASSL: A differential/algebraic system solver. SAND82-8637, Sandia National Laboratories, Livermore. 1982.
- Kirkpatrick S, Gelatt CD, Vecchi MP. Optimization by Simulated Annealing. Science. 1983;220:671–680. doi: 10.1126/science.220.4598.671. [PubMed] [Cross Ref]
- Kunst A, Draeger B, Ziegenhorn J. Methods of Enzymatic Analysis. Verlag Chemie, Weinheim, Germany; 1984. D-Glucose; UV-methods with hexokinase and glucose-6-phosphate dehydrogenase.
- 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,2-13C2]glucose. Am J Physiol. 1998;274:E843–E851. [PubMed]

Articles from BMC Systems Biology are provided here courtesy of **BioMed Central**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |