Journal of Integrative Bioinformatics, 7(3):129, 2010

http://journal.imbio.de

A Model of Phospholipid Biosynthesis in Tumor in Response to an Anticancer Agent in Vivo

Mahsa Behzadi1,* ,Aicha Demidem2 ,Daniel Morvan2,4 ,Laurent Schwartz1,3 ,Georges Stepien2 ,Jean-Marc Steyaert1
1

LIX, Ecole Polytechnique, Palaiseau Cedex 91128, France.

INRA UMR1019, Centre de Clermont-Ferrand Theix, 63122 Saint-Gens Champanelle, France.

` Service de Radiotherapie, Hopital de la Pitie-Salpetriere, Boulevard de lHopital 75013 Paris.

4

Universite Clermont-1, Boulvard Mitterand, 63001 Clermont-Ferrand, France. Summary

We study, in this paper, a model for the core of the system of the Glycerophospholipid metabolism in the murine cells. It comprises the simple and enzymatic reactions of PhosphatidylEthanolamine and the PhosphatidylCholine . The models general structure is taken from a number of books and articles. We translate this model into a set of ordinary differential equations (ODEs), to propose a quantitative explanation of the experimental experiences and the observed results. In order to make it usable as a basis for simulations and mathematical analysis we need to make precise the various constants present in the equations but which are usually not directly accessible in the literature. In a rst step we considered experimental data of rats liver cells obtained by NMR spectroscopy: given the values of metabolite concentrations we nd appropriate parameter values which allow us to describe the system with ODEs. We have then performed several analyses using the developed model such as stability analysis. A rst interesting result is the global stability of the system which was observed by simulation and then proved by mathematical arguments. A second important result is that we observe on the diagrams that the steady state for normal cells is precisely a singular point of order two, whereas tumoral cells present different characteristics; this fact has been proved for PhosphatidylEthanolamine N-Methyl transferase (PEMT), an enzyme which seems to be identied for the rst time as a crucial element in the tumoral process. In a second step we applied our model to experimental data of proton HRMAS NMR spectroscopy for solid B16 melanoma and Lewis lung (3LL) 3LL carcinoma cells treated by Chloroethyl Nitrosourea (CENU). We performed a complete comparative analysis of parameters in order to learn the predictive statements to explain increases and decreases which one can observe in concentrations.

1 Introduction
Phospholipids are a major component of biological membranes. They are a class of lipids formed from four components: fatty acids, a negatively-charged phosphate group, alcoholamine and a backbone. PhosphatidylCholine (PtdCho) and PhosphatidylEthanolamine (PtdEth) are two of the most abundant Phospholipids. In most eukaryotic cells, PtdCho is synthesized through two different pathways [1]; in the cytidine diphosphate-choline (CDP-choline)
* mahsa.behzadi@polytechnique.edu

doi:10.2390/biecoll-jib-2010-129

Journal of Integrative Bioinformatics, 7(3):129, 2010

http://journal.imbio.de

pathway (Kennedy pathway) and via the transmethylation of PtdEth catalysed by PhosphatidylEthanolamine N-Methyl transferase (PEMT). Choline, supplied by food, is principally in the form of PtdCho but also exists as free Choline [2, 3]. The Kennedy pathway for producing PtdCho, involves the activation of Choline (Cho) to CDPcholine through an intermediate product, PhosphoCholine (P-Cho). The second pathway to produce PtdCho consists of three sequential methylations of PtdEth. Cho derived from turnover of PtdCho produced by the methylation pathway is used for PtdCho synthesis through the Kennedy pathway. Therefore the activity of the Kennedy pathway does not reduce even in the absence of Cho in the growth medium [4]. In 1975, Sundler et al. used radioisotope methods to examine the rates of synthesis for PtdCho and PtdEt of liver [5, 6]. In their study, there are still questions about metabolic pathways to be answered. However evidence of the two different pathways of PtdCho synthesis and the relative activities of these pathways was provided by Vance et al. [7, 8]. The Nuclear Magnetic Resonance (NMR) spectroscopy method has been used to study the biosynthesis of PtdCho and PtdEth [9, 10]. The NMR technique can also provide a detailed examination of the specic metabolic pathways. Reo et al. performed kinetic analyses of liver PtdCho and PtdEth biosynthesis using 13 C NMR spectroscopy [11]. The development of methods for pathway-specic analyses of phospholipid biosynthesis in intact tissue can help in our understanding of numerous cellular processes, and may be important for cancer studies. This is why the Phospholipid metabolism has attracted the attention in cancer research. It is of interest to biologists to be able to follow the phospholipid metabolism in circumstances in which cell survival and cell proliferation are of concern, e.g. neurological disorders and cancer [12, 13]. Thus there is a need to develop a model for their biosynthesis and turnover. This is why we tried to nd a model for the GlyceroPhospholilid metabolism in the murine liver cell. Our goal is to build a model with which one could simulate the behavior of Phospholipids interactions. Due to the complexity of this system, mathematical modeling and numerical simulation is necessary to enable a compact representation of the current knowledge and to make meaningful quantitative predictions guiding future experimental studies. Once the model is developed to represent the metabolism of the murine liver cell, one can apply it to study the Phospholipid metabolism of mouse melanoma and 3LL carcinoma cells. In the recent years several studies have been carried out to perform chloroethylnitrosourea (CENU) chemotherapy for the treatment of B16 melanoma and Lewis lung (3LL) carcinoma tumors in vivo [14, 15, 16, 17, 18]. We apply our model to study the effects of such treatments. For each of these two tumors we have experimental data for three different phases: Control(CTL), Inhibition(INH) and Recovery(REC)[27].The results of our comparative analysis, based on our simulation, show good agreements with experimental data [16]. In our study, we have directly translated the biochemical reactions into ODEs, following the Michaelis-Menten chemical paradigm. Other authors have used probabilistic models or computer science models based on -calcul; under very reasonable assumption all these models yield ODEs as stated by L. Cardelli in his series on Articial Biochemistry [19].

2 Methods
In this section, we rst describe our model for the Phospholipids metabolism which is supplied from bibliographical references(e.g. M.Israel and L.Schwartz [20, 30, 31, 32, 33, 34]). Then we introduce our methods to obtain experimental data. Next, we develop an ODEs-based model. Finally, we study different phase spaces of the system.
doi:10.2390/biecoll-jib-2010-129 2

Journal of Integrative Bioinformatics, 7(3):129, 2010

http://journal.imbio.de

Figure 1: Schematic representation of the model. Arrows with V Mi and KMi parameters refer to enzymatic reactions while the rest represent simple reactions. Reactants: Cho(Choline), PCho(Phospho-Choline), PtdCho(Phosphatidyle-Choline), GPC(GlyceroPhosphoCholine), Eth(Ethanolamine), PEth(Phospho-Ethanolamine), PtdEth(PhosphatidyleEthanolamine), GPE(Glycero-PhosphoEthanolamine). Enzymes: CK(Choline-Kinase), EK(Ethanolamine-Kinase),CPT(PC-transferase),CCT(PhosphoCholine-Cytidyl-Transferase), ECT/EPT(PhosphoEthanolamine-Cytidyl-Transferase), PEMT(PhosphatidyleEthanolamine-Nmethyl-Transferase),PlpA2 (PhosphoLipase A2), PlpC(PhosphoLipase C), PlpD(PhosphoLipase D). Parameters: VM(Michaelis maximum reaction rate), KM(Michaelis concentration constant), k1 -k12 ( Rate constants for external reactions).

2.1 Biochemistry of the phospholipid metabolism Our analysis concerns twenty-four biochemical reactions (Fig. 1). In this system there are two main sub-systems with similar reaction structures; the rst one is the Choline (Cho) cycle and the second one is the Ethanolamine (Eth) cycle. In order to have a more complete model several reactions involving external reactants are also considered in the model (Fig. 1). Choline cycle: Cho is phosphorylated in a reaction catalyzed by Choline-Kinase (CK), resulting in the formation of PhosphoCholine (PCho)[30]. PCho is converted to PtdCho in a two step reaction, rst catalyzed by PhosphoCholine-Cytidyl-transferase (CCT), then by PCtransferase (CPT)[30, 31]. PtdCho is converted to Glycero-PhosphoCholine (GPC) in the reaction catalyzed by Phospholipase A2 (PlpA2)[31]. In addition, PCho and Cho can be synthesized from hydrolysis of PtdCho through the reactions catalyzed by Phospholipase C (PlpC) and Phospholipase D(PlpD) respectively[31, 32]. Cho can be also synthesized from GPC[32]. Ethanolamine cycle: Eth is phosphorylated in an enzymatic reaction catalyzed by Ethanolamine-Kinase (EK), resulting in the formation of PhosphoEthanolamine (PEth)[33]. PEth is converted to PtdEth in a two steps reaction, rst catalyzed by PhosphoEthanolamine-Cytidyltransferase (ECT), then by PE-transferase (EPT)[33, 34]. PtdEth is converted to GlyceroPhosphoEthanolamine (GPE) in the reaction catalyzed by Phospholipase A2 (PlpA2)[8, 31]. Eth is synthesized from GPE[32]. The above two sub-systems are related through the reaction between PtdEth and PtdCho where PhosphatidylEthanolamine N-Methyl Transferase (PEMT) plays the role of catalyst[4, 7, 8]. This reaction seems to be an important reaction in this system, and the basis of main analysis in our study, since homeostasis of PtdCho is essential to maintain cell survival.

doi:10.2390/biecoll-jib-2010-129

Journal of Integrative Bioinformatics, 7(3):129, 2010

http://journal.imbio.de

External reactions: In addition to the reactions described so far, most of the reactants in phospholipids metabolism model have external reactions. For example there is a reversible reaction in which PhosphatidylSerine (PtdSer) releases CO2 and PtdEth as products [21]. In the same way there are several external reactions in which Cho, Eth, PCho, PEth, PtdCho and PtdEth have the role of substrate or product. We present these external reactions by input or output arrows in the model(Fig. 1)[9, 21, 30]. 2.2 Mathematical simulation The mathematical simulation of the model is performed via a free mathematical software, Scilab. For the simulation of the biological system we use traditional reaction-rate approach by dening the equations describing the system and setting the initial parameters required for the calculation. In this approach, the chemical reactions are modelled by ordinary differential equations (ODEs) representing the concentrations of the substances. We formulate the basic model for the given chemical reactions in terms of a system of differential equations , which consists of one differential equation for the kinetics of each of the reactants. In each equation, [X] represents the concentration of a given reactant X, whose values are expressed in mol.g 1 . The ki is the reaction rate coefcient or rate constant in simple reactions. For enzymatic reactions the parameters V Mi and KMi are the maximum rate and Michaelis constant respectively. The model writes as follows:
V M2 .[P tdCho] [Cho] = k2 .e[Cho] + t KM2 + [P tdCho] + V M1 .[Cho] V M3 .[GP C] KM3 + [GP C] KM1 + [Cho] (1) V M4 .[P Cho] [P tdCho] = k12 .e[P tdCho] + t KM4 + [P Cho] + V M6 .[P tdCho] V M7 .[P tdEth] KM7 + [P tdEth] KM6 + [P tdCho] V M5 .[P tdCho] V M2 .[P tdCho] KM5 + [P tdCho] KM2 + [P tdCho] V M12 .[P tdCho] k5 .[P tdCho] KM12 + [P tdCho] (5)

k1 .[Cho]

V M10 .[GP E] [Eth] = k10 .e[Eth] + t KM10 + [GP E] V M11 .[Eth] k11 .[Eth] KM11 + [Eth] (2)

V M12 .[P tdCho] [P tdEth] = k7 .e[P tdCho] + t KM12 + [P tdCho] + (3) V M6 .[P tdCho] V M3 .[GP C] [GP C] = t KM6 + [P tdCho] KM3 + [GP C] (7) V M9 .[P Eth] V M8 .[P tdEth] KM9 + [P Eth] KM8 + [P tdEth] V M7 .[P tdEth] k6 .[P tdEth] KM7 + [P tdEth] (6)

V M1 .[Cho] [P Cho] = k4 .e[P Cho] + t KM1 + [Cho] + V M5 .[P tdCho] KM5 + [P tdCho] V M4 .[P Cho] k3 .[P Cho] KM4 + [P Cho]

V M11 .[Eth] [P Eth] = k9 .e[P Eth] + t KM11 + [Eth] V M9 .[P Eth] k8 .[P Eth] KM9 + [P Eth] (4)

V M8 .[P tdEth] V M10 .[GP E] [GP E] = t KM8 + [P tdEth] KM10 + [GP E]

(8)

In this model, each of these differential equations expresses the rate of change of one reactant as a sum of fractional terms for enzymatic reactions and non-fractional terms for simple reactions. Furthermore, we proposed an exponential formula with concentrations of studied reactants in our system as variables, to explain the kinetics of reactions with reactants from the external environment. Diffusion phenomena are the reason of this exponential form. Molecular

doi:10.2390/biecoll-jib-2010-129

Journal of Integrative Bioinformatics, 7(3):129, 2010

http://journal.imbio.de

diffusion, often called simply diffusion, is a net transport of molecules from a region of higher concentration to one of lower concentration by random molecular motion. 2.3 Model analysis for healthy liver cells In this section, we use the described mathematical model in two steps: rst we try to nd the different required parameters in system of equations, such as the rate constants for each reaction, using the experimental values of the concentrations in healthy rat liver metabolism. The second step is to study the phase spaces diagrams and also the different stability analyses using the parameters obtained for healthy liver cells. 2.3.1 Concentrations and Parameter estimations The described kinetic equations require different parameters, such as the rate constants for each reaction. In the rst application of our mathematical model, the experimental values that we use, are derived from the concentrations of rat liver metabolism measured at several instants during infusion with Choline and Ethanolamine [11]. Concentration of [PtdCho] and [PtdEth] were measured from the 31 P NMR spectra of the lipid extracts [11]. A description of these experimental analysis is given in [23, 24]. The parameter values have an important effect on the precision of the model which is representing this biological system. However these values, ki , V Mi and KMi , are difcult to estimate experimentally and many are unknown; that is why we estimate them by means of a numerical method. There are several ways of doing so. Since the changes in the concentration values in the different instants are small, we predict that the system is initially close to its steady state. We use this assumption to treat our system of ODEs as a nonlinear algebric problem (since all derivatives are zero).This simplication now leaves a system of equations with a number of variables greater than the number of equations. Therefore mathematically we have a set of possible solutions of this ODE system. Furthermore we know that the vector of these rate constants needs to insure the behavior of the model in such a way that the cell is viable. For example there exists a specic limited range of concentration, for some metabolites, in which the cell can stay alive. These ranges can give upper and lower bounds for parameters.[25]. Once we take into account all these biological constraints, the possible rate constant vectors fall into a subset of the parameter space. Characterizing this subset would be a prediction of the model, and so would be characterizing the set of all the dynamics of the model consistent with the parameter vectors in this subset. The vector of parameters which are shown in Table 2 (cf. Annexe 1) is one of these possible solutions. To obtain this, we rst dene the sum of squares of rate equations of section 2.2 as a function. Then we nd a solution which minimizes the value of this function, in the viable range for parameters. For this aim, since we have two inequality constraints as a lower bound and an upper bound for each of parameters, we use the Karush-Kuhn-Tucker(KKT) theorem [35]. The KKT conditions are necessary for a solution in nonlinear programming to be optimal, provided some regularity conditions are satised. It is a generalization of the method of Lagrange multipliers to inequality constraints[35]. This method gives a solution that best ts with the biological measurements. Table 1 (cf. Annexe 1) represents the average of concentrations which are measured experimentally, for 8 reactants of our system. We rst try to use the average values presented in Table 1 (cf. Annexe 1), to obtain a possible vector of the constants of reactions ki , the maximum velocity and kinetic constant of Michaelis-Menten model V Mi and KMi (Table 2 (cf. Annexe 1)). We also try to nd a possible vector of the constants, for each set of concentrations measured in each of 6 time points [11]. Comparing the vectors obtained in each of these two cases, we do
doi:10.2390/biecoll-jib-2010-129 5

Journal of Integrative Bioinformatics, 7(3):129, 2010

http://journal.imbio.de

not observe big variations in the parameters, which means this is a robust solution. The results and variations are shown in Table 2 (cf. Annexe 1). We also observe the similar results for up to 100 other vectors in the subset of parameter space.

Figure 2: Changes in initial concentrations. A: Changes of initial concentration of PhosphatidylCholine. B: Phase space for Choline and PhosphoCholine. C: An Example of different initial concentrations traversing into steady state. In A and B the color change from red to blue refers to approaching the steady state. In C each color associates to the concentration of one of the reactants. Concentrations values are given in mol.g 1 .

2.3.2 Phase spaces and Stability Analyses Phase spaces: Here we calculate and study the phase spaces diagrams using the parameters obtained in our model. The goal is to obtain the behaviour of the system with respect to time. The simulation results showed that for the obtained parameter values there exists only one steady state point in the range near the studied initial concentration values. When we try to change one or several initial concentrations at time t0 , we see that after passing a period of time the concentrations of all the reactants converge nally toward the concentrations of the steady state point. Therefore one can conclude that the change in the initial concentration of each of the reactants does not modify the behaviour of the system at innity (+). (Fig. 2). Stability Analyses: As a rst set of stability analyses we study the changes in steady states by modifying the concentrations of enzymes in enzimatic reactions. To determine the maximum rate of an enzymatic reaction (like most of the reactions in our model) we used the MichaelisMenten model. In this model, the maximum initial velocity (a kinetic constant of the enzymatic reactions) reects the activity of an enzyme and is proportional to its concentration. Therefore, in our simulation and in order to represent changes in enzyme concentrations, we simply modied the value of the maximum initial velocities of the reactions (V Mi ). For example in the case of the reaction between PtdEth and GPE, the diagrams of changes of steady state resulting from the change of maximum velocity are shown in Fig. 3(a). On each of these diagrams, each point (VM, Xi) corresponds to a concentration of the reactant Xi at the steady state.
doi:10.2390/biecoll-jib-2010-129 6

Journal of Integrative Bioinformatics, 7(3):129, 2010

http://journal.imbio.de

Figure 3: Steady State Concentration vs. Enzyme Concentration. A: Change of steady state point for the reaction of PtdEth and GPE. On each of these diagrams, each point (VM, Xi) corresponds to a concentration of the reactant Xi at the steady state. B: Change of steady state point for the reaction of PtdEth and PtdCho. The red point in each of these diagrams is associated to the concentration at the steady state of the experimental values. Concentration values are given in mol.g 1 .

The diagrams were as expected. When the velocity of reaction from PtdEth to GPE increases, the concentrations of GPE, PEth and Eth increase and the concentrations of PtdEth, Cho, PtdCho, PCho and GPC decrease. Here it is worth to recall the special role of PtdEth N-methyltransferase enzyme(PEMT) in phospholipids biosynthesis. The PEMT pathway is especially functional in the liver. It is a minor pathway for PtdCho synthesis from PtdEth, the major pathway being the Kennedy pathway which involves the metabolism of Choline taken from the blood. The PEMT pathway is implicated in the biosynthesis of lipoproteins in the liver. In contrast, its role is poorly known in tumors, excepted in hepatocarcinoma [26]. It was recently shown that some tumor cell types could compensate for the deciency of the Kennedy pathway by upregulating the PEMT pathway, thus surviving. Therefore a complementary aim for our modelling was to get further insight into the possible role of the PEMT pathway in the regulation of tumor phospholipid metabolism. To our knowledge, the implication of PEMT in response to an anticancer agent, in melanoma and 3LL carcinoma has never been investigated. Now we get back to our model and take the reaction of rst order, related to PtdEth and PtdCho for which PtdEth N-methyltransferase (PEMT) plays the role of enzyme. Fig. 3(b) represents the changes of the steady state point associated to different reactants. On each of its diagrams, every point ([PEMT], Xi) corresponds to a concentration of reactant Xi at steady state. So if we change the enzyme concentration, for an arbitrary value of the initial concentration like Xi0 the concentration [Xi ] will tend towards its steady state concentration value. In Fig. 3(b), the red point (also shown by a circle) in each of the diagrams is associated to the concentration of the steady state among the experimental values. In the second set of stability analysis, we studied the effects of the changes in reaction rates on concentrations and their relation to the steady state. For this aim, we performed several different simulations for each reaction. In the diagrams presented in the example of Fig. 4(a), each colour is associated to one of simulations. Each of these diagrams represents the changes of concentration of a reactant vs. its reaction rate for that reactant. Let us recall that the point where the rate of reaction reaches zero, is called the steady state point for a given reactant. To make it clearer, let us explain the corresponding diagram for PtdEth. This is also shown in Fig. 4(b) with a better resolution. The steady state is obtained when the concentration is around 11 mol.g 1 . This is obtained by taking some random concentrations and measuring the reaction rate for each of them. Then we connected the resulting points to see when the zero rate is obtained. For any other points in this diagram, which are not the steady
doi:10.2390/biecoll-jib-2010-129 7

Journal of Integrative Bioinformatics, 7(3):129, 2010

http://journal.imbio.de

Figure 4: Rate vs. Concentration. A: Rate vs. Concentration for PhtdEth: Each experiment is shown by a different colour in this diagram. The steady state is obtained when the concentration of PtdEth is around 11 mol.g 1 . B: Rate vs. Concentration: Each diagram corresponds to the change in concentration of one reactant vs. its reaction rate. We dene the slope (k) as rate coefcient which is used as a parameter in analysis concerning the speed of reaching the steady state. Values are given in mol.g 1 for Concentrations and in mol.g 1 .s1 for Rates .

state, concentration tends towards the steady state concentration. In other words the normal cell behaviour corresponds to a superstable steady state. As shown in Fig. 4(a) and Fig. 4(b), for all of the reactants, these diagrams have shown to be linear and they coincide for any xed reactant. In the third set of stability analyses on the concentrations of rat liver metabolism, we tried to study the speed to reach the steady state point. In Fig. 4(b) for the points with the concentrations far from the steady state point (which is shown by an arrow), the absolute value of this rate is bigger than for the points which have a concentration close to the steady state. This means that, the speed to reach the steady state point increases when we try to change the concentration of reactants. One of the parameters which could inuence this speed is the concentration of enzyme. To study the effect of the change of concentration or activity of an enzyme on this speed, as it is shown in Fig. 4(b), we considered the slopes of the diagrams as an indicating coefcient for the reaction speed. Let us call this slope the rate coefcient (k). If the concentration of one of the enzymes changes, the rate coefcient does so. When the rate coefcient is small, the rate changes more slowly than when the rate coefcient is larger. The diagrams of rate coefcients vs. concentration of enzyme PEMT are shown in Fig. 5.

3 Results and Discussion

In this Section we will rst study some interesting results obtained from the stability analyses of section 3.2. We will also discuss the complexity of the algorithm applied to Rate vs. Concentration stability analysis. Then in the next step we give a mathematical proof for the stability of this model of equations. Finally as the most interesting application of our proposed mathematical model, we provide a complete comparative analysis using the experimental data which come from chloroethylnitrosourea (CENU) chemotherapy for the treatment of B16 melanoma and Lewis lung (3LL) carcinoma tumoral mice cells in vivo [19].

doi:10.2390/biecoll-jib-2010-129

Journal of Integrative Bioinformatics, 7(3):129, 2010

http://journal.imbio.de

Figure 5: Speed Analysis. -k (Rate Coefcients) vs. Concentration of enzyme PEMT(real values). The red point on each diagram is associated to rate coefcient at the steady state of the experimental values. Concentrations values are given in mol.g 1 .

Figure 6: Treatment Phases. Growth curves of untreated (white circles) and CENU-treated tumors during the growth inhibition phase to treatment (black circles) and the growth recovery phase (gray circles). CENU was given intratumorally at days 11, 14, and 18. Bars, SD.

3.1 Stability analyses results and Complexity study The analysis of reactants concentrations vs. enzymes concentrations in section 3.2. shows that: The PEMT Enzyme is found to provoke a reciprocal trend between PEth and PCho. This means that there is a balance between the change of concentration of PEth and the change of concentration of PCho. This result is in good agreement with experimental biological results reported in B16 melanoma cell cultures responding to an anticancer agent [14]. One can note from Fig.3(b) that when the concentration of PEMT increases the concentration of Cho and PC(PCho) saturate after certain values of [PEMT], while the concentration of PtdCho still increases. This ts with the fact that the methylation of PtdEth may relay PtdCho biosynthesis when the choline pathway is saturated or blocked[7, 8, 26]. The interesting result about the red points (also shown by a circle) in Fig.3(b), is that, these positions are the places where the behaviour of diagrams changes (either inexion point in mathematical term which means where the second derivative is zero or in some situations it is the third derivative which vanishes).

doi:10.2390/biecoll-jib-2010-129

Journal of Integrative Bioinformatics, 7(3):129, 2010

http://journal.imbio.de

In the stability analyses of Rate vs. Concentration in section 3.2., given the values of concentrations {Ci } which were observed experimentally, we managed to nd appropriate parameter values {Pi }. With these parameter values, the ODEs system has a stable solution, and the resulting concentration values are equal to the initial ones ({Ci }). However, one could ask whether a change in parameter values could give an unstable or oscillating solution. For that purpose we studied the eigenvalues of the Jabobian matrix in a number of points close to {Ci }. As we had 41 parameters we just could not try all variations of them at the same time. For instance, even to try 10 values for each parameter would take O(1041 ) operations. So we used random numbers to change all parameter values simultaneously. In each of 10 000 experiments we made, each parameter took a value, where rand factor was an uniformly distributed random value in a range from 0 to 20. In all the experiments we observed the stability of the system. 3.2 Mathematical proof of stability (sketch) We tried to nd a xed point for the proposed ODEs system using the parameter values that we found in section 3.1, such as it would be the closest xed point to {Ci }. The eigenvalues for Jacobian matrix of the system in such xed points are always real and negative. This indicates that the solution is always stable and without oscillations regardless of parameters values of the system. More generally if we pay attention to our ODEs system where most of the equations have the general form of Michaelis-Menten, we can prove that all solutions for this system are stable. The proof is that we have a Jacobian matrix with a dominant diagonal. This means that, in every row of this matrix, the magnitude of the diagonal entry in the row is larger than the sum of the magnitudes of all the other (non-diagonal) entries in that row. More precisely, the matrix A is diagonally dominant if |aii | j=i |aij | for all i, where aij denotes the entry in the ith row and jth column. The Jacobi method for solving a linear system converges if the matrix is diagonally dominant. The eigenvalues for the Jacobian matrix in such systems are always real and negative. This indicates that the solution is always stable and without oscillations regardless the particular parameters of the system. In the other word all the evolutions of the cell metabolism are stable in the proposed mathematical model which is based on MichaelisMenten kinetics. 3.3 Application of the model on Tumors treated by CENU 3.3.1 In vivo experiments: Treatment protocol Six-to-eight-weeks-old C57BL6/6J male mice were purchased from IFFA CREDO, (LArbresle, France). Mice were shaved before s.c. injections into their ank of 5105 tumor cells(B16 melanoma or 3LL cells). B16 melanoma or 3LL tumors became palpable at days 8-10 after cell inoculation. Mice were divided into two groups, a Control (CT) group, which received sham injections of saline solution, and a Treated (TR) group. The TR group received intratumor chloroethylnitrosourea (CENU) injections at a dose of 15 g/g body weight. CENU was injected at days 11, 14, and 18 from B16 cell inoculation. At dened times of tumor evolution (days 10, 12, 15, 20, 24, and 29 after B16 cell inoculation for CT and TR tumors and prolonged to days 35, 43, and 54 for TR tumors), three mice of each group were sacriced according to institutional guidelines for animal welfare and experimental conduct. Tumors were dissected and weighed. The dissection of the s.c. tumor took <2 min. A piece of the tumor <50 mg of the whole tumor was immediately prepared for NMR Spectroscopy as described below or frozen at 80 C in case of delayed examination.

doi:10.2390/biecoll-jib-2010-129

10

Journal of Integrative Bioinformatics, 7(3):129, 2010

http://journal.imbio.de

Tumor growth curves were tted to the Gompertz function and modied for the TR group to include a growth delay period. Ranges for maximum attainable weights and other model parameters obtained in CENU-treated melanoma models have been published previously [22]. 3.3.2 Results and discussion The aim of this section is to provide insights into metabolic pathways from biochemical data derived from 1H-NMR spectroscopy-based metabolite proling of tumors [27]. Proton twodimensional NMR spectroscopy analysis has been shown to cover very well the subset of phospholipid derivatives [14, 16, 29], including the most concentrated phospholipids (PtdCho and PtdEth), water-soluble precursors (choline, phosphocholine, cytidyl-diphosphocholine, ethanolamine, phosphoethanolamine, cytidyl-diphosphoethanolamine), phospholipid hydrolysis products (glycerophosphocholine, glycerophosphoethanolamine), and oxidization products (betaine). Besides these technical conveniences, phospholipid metabolism is crucial for the build-up of cellular membranes thus for tumor cell proliferation, a major phenotypic feature of tumors. Recently, as an anticancer treatment strategy, it was proposed to inhibit key-enzymes of phospholipid metabolism (choline-kinase) to slow down tumor cell proliferation [28]. We hypothetized that, by modelling phospholipid derivative content variations between two conditions at steady state, we could give insight, through the used set of parameters, into the induced regulations of phospholipid metabolism. We thus compared phospholipids metabolism alterations in murine tumors between baseline and the stable phase of their response to an anticancer agent. Based on the classical hypothesis that pathways of phospholipids metabolism are very similar in liver cells and tumor cells [9], we applied our mathematical model to study the effects of such treatments. For each of these two tumors we have experimental data for three different phases: Control(CTL), Inhibition(INH) and Recovery(REC) [27](See Fig. 6). The average concentrations measured experimentally at steady state for each of these phases are shown in Table 3 and 4 (cf. Annexe 1). At a rst step we tried to obtain a possible vector for appropriate parameter values for each phase of treatment, applying the same methodology we used for the liver cell metabolites. The results are shown in Table 5 (cf. Annexe 1). Afterwards a complete comparative parameter analysis is performed in order to understand the background of the observed increaments or decrements of concentrations(Table 3 and 4 (cf. Annexe 1). Let us explain an example in more details: As it is shown in Table 3 (cf. Annexe 1), for B16 melanoma, we observe an increase of 5.33 mol.g 1 melanoma for Phosphoethanolamine(PEth). The reason is easily explainable by our parameter analysis for change of rate of all reactions in which PEth plays the role of a substrate or product. Here, PEth concentration increases since V M11 is increasing rapidly and V M9 is decreasing, despite PEMT activity is decreasing. By a similar analysis we obtain the explanation for the decreases in GPC and GPE in 3LL carcinoma [16].

4 Conclusion
Understanding cell metabolism evolution and changes is for many scientists more than a challenge; it is the key to a thorough understanding of cell dysfunction and very likely a step toward the elucidation of carcinogenesis along the lines of Warburgs seminal papers. In this paper we presented mathematical analysis of the metabolic pathways which control and command the production of Glycerophospholipids through the enzymatic reactions of PtdEth and
doi:10.2390/biecoll-jib-2010-129 11

Journal of Integrative Bioinformatics, 7(3):129, 2010

http://journal.imbio.de

PtdCho. Analysis shows that the normal cell stands at very special points of equilibrium. We also checked our model against a series of experiments and gave evidence for the crucial role of PhosphatidylEthanolamine N-Methyl transferase (PEMT) in tumor cells under CENU treatment. Our results show that: -The model ts in vivo observations and experiments with CENU tumor inhibitor, and provides new hypotheses on metabolic pathway activity from metabolite proling of phospholipids derivatives. and provides: -All the evolutions of the cell metabolism are stable in the Michaelis-Menten formula. -The normal cell behaviour corresponds to a superstable steady state.

Acknowledgements
Authors would like to thank M. Regnier, researcher in INRIA Rocquencourt, France. We also thank E.A. Vassilieva, researcher in the French National Center for Scientic Research (CNRS) who is working at the Laboratory of Computer Science (LIX) in Ecole Polytechnique, Paris.

References
[1] Henry S.A., Patton-Vogt J.L., Genetic regulation of phospholipid metabolism in yeast. Progress in Nucleic Acid Research and Molecular Biology 61,133-179,1998. [2] Blusztajn JK., Choline a vital amine., Science ,281,794-795,1998. [3] Zeisel SH., Choline A nutrient that is involved in the regulation of cell proliferation, cell death, and cell transformation., Adv Exp Med Biol,399,131-141,1996. [4] Patton-Vogt J.L., Griac P., Sreenivas A., Bruno V., Dowd S., Swede M.J.,Henry S.A., Role of the yeast phosphatidylinositol/phosphatidylcholine transfer protein (Sec14p) in phosphatidylcholine turnover and INO1 regulation, J. Biol. Chem., 272, 20873-20883,1997. [5] Sundler R., Akesson B.,Regulation of phospholipid biosynthesis in isolated rat hepatocytes,Biochem.J.,146,309-315,1975. [6] Sundler R.,The phospholipid substrates in the Ca2+ -stimulated incorporation of nitrogen bases into microsomal phospholipids,Biochim.Biophys.Acta,306,218-226,1973. [7] Vance D.E., Walkey C.J., Roles for the methylation of phosphatidylethanolamine,Curr. Opin. Lipidol.,9,125-130,1998. [8] Vance D.E., Houweling M., Lee M., Cui Z., Phosphatidylethanolamine methylation and hepatoma cell growth, Anticancer Res.,16,1413-1416, 1996. [9] Podo F., Tumor phospholipid metabolism NMR, Biomed,12,413-439,1999. [10] Ackerstaff E., Pug B.R., Nelson J.B., Bhujwalla Z.M., Detection of increased choline compounds with proton nuclear magnetic resonance spectroscopy subsequent to malignant transformation of human prostatic epithelial cells,Cancer Res,61,3599-3603,2001.

doi:10.2390/biecoll-jib-2010-129

12

Journal of Integrative Bioinformatics, 7(3):129, 2010

http://journal.imbio.de

[11] Reo N.V., Adinehzadeh M., Foy B.D., Kinetic analyses of liver phosphatidylcholine and phosphatidylethanolamine biosynthesis using (13)C NMR spectroscopy, Biochim. Biophys. Acta,1580,171-188,2002. [12] Farooqui A.A., Horrocks L.A., Glycerophospholipids in brain:their metabolism,incorporation into membranes,functions, and involvement in neurological disorders.Chem Phys Lipids,106,129,2000. [13] Gupta C.M. Phospholipids in disease, In: Cevc G, editor. Phospholipids handbook. New York: Marcel Dekker,895-908,1993. [14] Morvan D., Demidem A., Papon J., Madelmont J.C., Quantitative HRMS Proton Total Correlation Spectroscopy Applied to Cultured Melanoma Cells Treated by Chloroethyl Nitrosourea: demonstration of phospholipid metabolism alternation, Mgan Reson Med,49,241-248,2003. [15] Morvan D., Demidem A., Guenin S., Madelmont J.C., Methionine-Dependence Phenotype of Tumors: Metabolite proling in a Melanoma model using L-[methyl13C]methionine and high-resolution magic angle spinning 1H-13C nuclear magnetic resonance spectroscopy. Mgan Reson Med, 55,984-996,2006. [16] Demidem A., Morvan D., Madelmont J.C.,Bystander effects are induced by CENU treatment and associated with altered protein secretory activity of treated tumor cells. A relay for chemotherapy? Int J of Cancer,119,992-1004,2006. [17] Merle P., Morvan D., Madelmont J.C., Caillaud D., Demidem A., CENU treatment induced bystander effects which are effective on parental and non-parental tumors and have a phospholipid metabolism proton NMR spectroscopy signature, Int J Oncol,29,637642,2006. [18] Merle P.,Morvan D., Caillaud D., Demidem A., Chemotherapy-induced Bystander Effect in Response to Several Chloroethylnitrosoureas: An origin independent of DNA damage? Anticancer Research,28(1A),21-7,2008. [19] Cardelli L., Articial Biochemistry.Technical Report TR-08-2006, The Microsoft Research - University of Trento Centre for Computational and Systems Biology,2006. [20] ISRAEL M., SCHWARTZ L., English Book: drome.Publisher:John Libbey,2-7420-0597-8,2006. Cancer,a dysmethylation syn-

[21] Kanfer J., Kennedy E.P., Metabolism and function of bacterial lipids.II. Biosynthesis of phospholipids in escherchia coli.J Biol Chem. Jun,239,1720-6,1964. [22] Demidem A., Morvan D., Papon J., De Latour M.,Madelmont J.C., Cystemustine induces redifferentiation of primary tumors and confers protection against secondary tumor growth in a melanoma murine model.Cancer Res,61,2294-2300,2001. [23] Adinehzadeh M., Reo N.V., Effects of Peroxisome Proliferators on Rat Liver Phospholipids: Sphingomyelin Degradation May Be Involved in Hepatotoxic Mechanism of Peruorodecanoic Acid, Chem Res. Toxicol, 11,428-440,1998.

doi:10.2390/biecoll-jib-2010-129

13

Journal of Integrative Bioinformatics, 7(3):129, 2010

http://journal.imbio.de

[24] Adinehzadeh M., Reo N.V., NMR analysis of liver phospholipids: temperature dependence of 31P chemical shifts and absolute quantitation, Q. Magn. Reson. Biol. Med,3,171-176,1996. [25] Cuia Z., Houwelingb M., Phosphatidylcholine and cell death. Biochimica et Biophysica Acta,2003. [26] Tessitor L. et al. Expression of phosphatidylethanolamine N-methyltransferase in human hepatocellular carcinomas, Oncology, 65,152-8,2003. [27] Morvan D., Demidem A., Metabolomics by proton nuclear magnetic resonance spectroscopy of the response to chloroethylnitrosourea reveals drug efcacy and tumor adaptive metabolic pathways. Cancer Res, 67,2150-2159,2007. [28] Hernando E., Sarmentero-Estrada J., Koppie T., Belda-Iniesta C., Ramrez de Molina V. et al., Lacal JC.A critical role for choline kinase-alpha in the aggressiveness of bladder carcinomas.Oncogene,28,2425-35,2009. [29] Morvan D., Demidem A., Papon J., De Latour M., Madelmont J.C.,Melanoma tumors acquire a new phospholipid metabolism phenotype under cystemustine as revealed by high-resolution magic angle spinning proton nuclear magnetic resonance spectroscopy of intact tumor samples.Cancer research.62(6),1890-7,2002. [30] Kennedy E.,Roelofsen B. The biosynthesis of phospholipids. In Lipids and Membranes: Past, Present and Future,Elsevier Science,171-206,1986. [31] Vance D.E., Phosphatidylcholine Metabolism,CRC Press,22,33-46,1989. [32] Kanoh H.,Kai M. ,Wada I., Phosphatidic acid phosphatase from mammalian tissues: discovery of channel-like proteins with unexpected functions. ,Biochim. Biophys. Acta.,1348,56-62,1997. [33] Jelsema C.L. , More D.J.,Distribution of phospholipid biosythetic enzymes among cell components of rat liver. J Biol Chem,253,7960-7971,1978. [34] Miller M., Kent C.,Characterization of the pathways for phosphatidylethanolamine biosynthesis in Chinese hamster ovary mutant and parental cell lines. J Biol Chem,261,9753-9761,1986. [35] Karush W., Kuhn W., Tucker W.,Characterization of the pathways for Nonlinear programming, The Karush-Kuhn-Tucker Theorem. University of California Press.2,481492,1951.

doi:10.2390/biecoll-jib-2010-129

14

Journal of Integrative Bioinformatics, 7(3):129, 2010

http://journal.imbio.de

Annexe 1
Table 1: Initial concentrations of rat liver metabolites
Choline Ethanolamine PhosphoCholine PhosphoEthanolamine PhosphatidylCholine PhosphatidylEthanolamine GlyceroPhosphoCholine GlycerophosphoEthanolamine 0.510.11 0.110.01 0.880.18 1.020.15 29.572.76 9.851.15 0.410.09 0.350.05

Values are given in mol.g 1 liver, measured by P-NMR. doi:10.1016/S1388-1981(01)00202-5

Table 2: Estimated parameter values for rat liver metabolites

Param k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 Value 1.4650.091 1.8820.110 0.2810.013 2.9810.118 0.0640.007 0.0540.003 1.0010.081 0.7860.041 2.2060.109 1.2380.074 0.9800.079 1.0000.042 Param V M1 V M2 V M3 V M4 V M5 V M6 V M7 V M8 V M9 V M10 V M11 V M12 Value 2.2890.105 0.6240.076 0.8140.040 4.8980.122 0.5750.032 0.6960.009 10.4510.094 0.5050.008 3.6290.086 0.5770.002 2.8250.075 1.6570.009 Param KM1 KM2 KM3 KM4 KM5 KM6 KM7 KM8 KM9 KM10 KM11 KM12 Value 0.5670.015 29.7070.195 0.5490.033 0.9500.010 29.8120.252 29.6340.310 4.8450.046 10.0600.096 1.8880.015 0.4570.002 0.1210.011 29.8440.155

Values are given in mol.g 1 liver for KMi and in mol.g 1 .s1 for V Mi and s1 for ki .

Table 3: Average concentrations of B16 melanoma tumor model metabolites

B16: Cho Eth PCho PEth PtdCho PtdEth GPC GPE CTL 0.395 0.1 1.091 4.001 15.030 5.010 0.367 0.703 INH 0.636 0.102 1.847 9.337 25.782 8.594 1.561 1.401 REC 0.518 0.101 1.842 6.304 21.915 7.305 0.131 0.665

Values are given in mol.g 1 melanoma, measured by 1 H NMR spectroscopy.

Table 4: Average concentrations of 3LL carcinoma tumor metabolites

3LL: Cho Eth PCho PEth PtdCho PtdEth GPC GPE CTL 0.946 0.1 0.996 3.533 10.456 3.485 4.643 1.356 INH 1.746 0.991 1.610 6.458 11.470 3.823 1.370 1.025 REC 1.596 0.101 1.542 5.468 13.071 4.357 2.434 0.873

Values are given in mol.g 1 3LL carcinoma,measured by 1 H NMR spectroscopy.

doi:10.1002/ijc.21761

doi:10.2390/biecoll-jib-2010-129

15

Journal of Integrative Bioinformatics, 7(3):129, 2010

http://journal.imbio.de

Table 5: Estimated parameter values for mouse melanoma and 3LL carcinoma metabolites which is used in our comparative analysis
Cell: Parameter: k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 V M1 KM1 V M2 KM2 V M3 KM3 V M4 KM4 V M5 KM5 V M6 KM6 V M7 KM7 V M8 KM8 V M9 KM9 V M10 KM10 V M11 KM11 V M12 KM12 CTL 1.5800.28 1.5440.19 0.6640.07 2.4690.12 0.0690.04 0.0410.05 1.0550.04 0.0090.04 1.1150.04 0.9210.09 1.0160.02 1.0000.01 3.2530.08 1.3960.15 0.2900.02 15.190.03 0.8180.07 1.5590.03 3.1580.35 2.3980.20 0.3170.07 15.120.02 0.3130.06 15.150.01 2.4260.01 6.5390.50 0.4100.06 5.2950.01 1.8630.12 4.1300.16 0.4900.01 1.0270.07 1.9800.05 0.1120.04 1.0720.18 15.130.10 B16-Treated INH 0.9750.11 1.7870.04 0.2080.08 1.9070.01 0.0050.08 0.0640.02 1.0020.002 0.0290.01 1.0000.001 0.9130.19 1.0980.09 1.0000.001 2.8940.51 1.3450.23 0.4600.10 25.880.01 0.8710.14 2.1130.31 3.0380.05 2.9000.09 0.6690.15 25.850.003 0.7420.17 25.870.01 1.5970.26 9.4190.06 0.2620.01 8.7550.001 1.1590.15 9.5500.17 0.3350.03 2.2080.29 3.2840.26 0.2870.05 1.7540.09 25.820.005 REC 0.7040.28 2.4020.42 0.0110.01 2.2540.14 0.1130.01 0.0070.06 1.0070.04 0.0820.03 1.0120.001 1.9330.036 1.0280.09 1.0000.001 3.1480.29 0.6330.11 0.4020.036 22.040.02 0.7520.151 0.5440.06 4.5700.29 1.9200.31 0.9590.05 21.970.006 0.2930.05 22.010.31 3.6230.82 8.0420.45 0.4520.08 7.4820.03 2.7560.27 6.5510.14 0.7030.25 1.4270.09 3.6760.69 0.0960.02 1.2990.04 22.100.11 CTL 0.5030.19 1.3220.05 0.1970.11 1.8260.07 0.0930.01 0.0590.01 1.1480.07 0.0640.01 1.1000.04 1.0130.11 1.0090.01 1.0000.001 1.6130.18 1.1440.004 0.5850.06 10.610.006 0.8040.02 4.6810.01 3.2740.13 1.1300.06 0.6460.02 10.590.001 0.8060.14 10.600.001 2.1430.08 3.6090.06 0.4480.08 3.6980.01 1.6750.04 3.5170.08 0.4890.002 1.6910.05 2.3480.17 0.1270.017 1.2070.08 10.540.01 3LL-Treated INH 0.2040.01 1.3800.03 0.0660.006 1.7200.07 0.0210.014 0.1390.04 1.1270.007 0.0690.03 1.0080.04 1.2910.14 0.9740.012 1.0000.001 0.8920.18 1.8360.04 0.4400.06 11.590.21 0.7320.09 1.6650.13 2.1380.06 1.7610.18 0.6940.06 11.570.12 0.6620.05 11.530.08 1.5390.15 4.1650.25 0.2320.007 4.0980.32 1.5060.23 6.7640.27 0.2720.09 1.4660.2 2.7720.05 0.1340.017 1.2500.08 11.670.1 REC 0.3180.05 1.4200.08 0.0850.008 1.7390.19 0.0140.01 0.0960.05 1.0720.14 0.0890.01 1.0190.15 1.2660.21 0.9840.05 1.0000.001 0.9460.07 1.7220.16 0.5840.08 13.170.11 0.7970.01 2.6560.17 2.3750.15 1.8690.15 0.7510.07 13.140.01 0.7640.07 13.130.02 1.6990.01 4.6740.26 0.2880.06 4.5990.35 1.4510.19 5.7940.14 0.3030.09 1.0120.21 2.4530.27 0.1060.008 1.3340.07 13.280.14

Values (means estimateSD estimate) are given in mol.g 1 tumor for KMi and in mol.g 1 .s1 for V Mi and s1 for ki .

doi:10.2390/biecoll-jib-2010-129

16