BACKGROUND
sphingolipids  are categorized as lipids with a sphingoid base backbone  <cit>  that is often derivatized with an amide-linked fatty acid to make ceramides  and more structurally complex sl with diverse biological functions  <cit> .
sl in essentially every subcategory, from the lipid backbones  <cit>  to complex sl  <cit> , are highly bioactive and play important roles in diseases  <cit> ; hence, methods for "lipidomic" analysis of sl and sl metabolism are important for an in-depth understanding of these enigmatic compounds.
in recent years, a number of large-scale experimental and bioinformatics projects have begun to address the complexity of the lipidome.
examples include the lipid metabolites and pathways strategy  consortium  <cit> , the lipid library  <cit> , cyberlipid center  <cit>  and lipidbank  <cit> .
in particular, lipid maps has become a comprehensive resource for information on classification, structures and quantitative data on lipids and provides an opportunity for developing quantitative models of lipid synthesis and metabolism thus facilitating a mechanistic and systems-level understanding.
the de novo biosynthesis of sl begins with production of the sphingoid base, which utilizes serine and palmitoyl-coenzyme a  and various fatty acyl-coas to make n-acylsphinganines  that are desaturated to cer  and incorporated into more complex sl such as ceramide 1-phosphate , sphingomyelin , glucosyl- and galactosyl-ceramide  and more complex glycosphingolipids  <cit> .
ceramide can also be synthesized by recycling of sphingosine from turnover of sl such as sm  <cit> ; furthermore, sphingosine and sphinganine are phosphorylated to sphingosine 1-phosphate  and sphinganine 1-phosphate  which are intermediates of sphingoid base degradation  <cit>  and cell signaling molecules  <cit> .
due to the complexity of this pathway, and the paucity of data for its many metabolites, there are only a few models of sl metabolism available in the literature  <cit> .
the lipid maps consortium  <cit>  has quantified the global changes in lipid metabolites and genes in raw  <dig>  macrophage cells treated with kdo2-lipid a .
kla is the active component of inflammatory lipopolysaccharide which functions as a selective agonist of toll-like receptor  <dig>  and mimics bacterial infection.
the measurements are carried out over a 24-hour time period and the data is freely available via the lipid maps website  <cit> .
the goal of the work presented here is to develop a predictive kinetic model for sl metabolism using the lipidomics and transcriptomics data from the lipid maps studies.
this manuscript is organized as follows: we first briefly discuss the experimental data preprocessing and the methodology used for estimating the rate parameters, then we present the results of parameter estimation, followed by discussion and conclusions.
methods
network simplification
a detailed metabolic reaction network was developed using the information available in the literature and the kegg pathways database  <cit>  .
the c16-branch of cer biosynthesis  was selected for developing the model because this is a major subspecies for all categories of complex sl in the raw <dig>  cells.
vanted software was used to draw the reaction network  <cit> .
it is common in modeling studies for the network to contain several unmeasured nodes ; in our pathway , quantities are known for all of the metabolites and genes except dhgalcer and galcer  and are expressed as leaf-nodes  in the network.
one of the steps in our matrix-based fast algorithm for parameter-estimation requires experimental data on all metabolites except on the leaf-nodes/metabolites in the network.
a detailed procedure for simplifying the network, if the network contains unmeasured components, is described in our previous work  <cit> .
the leaf-nodes were exempted from the model described in this paper because the reactions leading to the unmeasured leaf metabolites were combined with the default degradation of their precursors.
experimental data and preprocessing
the lipid maps consortium has measured all the major lipids in mouse macrophage raw  <dig>  cells grown in 10% serum and treated with kla.
all metabolites were measured in pmol/μg-dna units.
the time-courses of lipids and genes were measured under two conditions:  normal condition designated as control and  inflamed condition .
time series comprising of  <dig> data points at  <dig>   <dig> ,  <dig>   <dig>   <dig>   <dig>   <dig> and  <dig> hr were measured with three biological replicates consisting of three technical replicates each.
the three technical replicate experiments were performed on the same day with a single batch of cells.
in addition, each time course was repeated three times on different days with different batches of cells .
for kinetic modeling purposes, outlier points were detected by a t-test and were excluded at each time-point.
the resulting data from all the replicates were averaged at each time-point.
data was processed for all metabolites under the treatment with kla and control conditions.
serine palmitoyltransferase activity in vitro was determined using the assay conditions described previously  <cit>  but with palmitoyl-coa as the labeled substrate.
after extraction of the products, the amounts of 3-ketosphinganine were determined by liquid chromatography, electrospray-ionization tandem mass spectrometry  <cit> .
the assay was conducted using  <dig> μg of protein obtained by sonication of raw <dig>  cells in the buffer used for the assay.
development of a kinetic model and parameter estimation
we have developed a kinetic model of the sl metabolism.
the procedure is similar to the one described in our previous work  <cit> , but it is presented here through an example of reactions from the sl network for completeness.
the reaction rates were described by linear or law of mass action kinetics with the assumption that for enzymatic reactions, the substrate concentrations are much smaller as compared to the corresponding michaelis constant, km.
for example, the following types of reaction schemes and rate expressions were used:
reactionrate expressionsa→bk*a+b→a.bk**a→bk**
the enzymes involved in sl metabolism can be regulated at multiple levels from mrna expression to posttranslational modification.
mrna data on the genes involved in the pathway is available from microarray experiments .
however, the corresponding proteomics data is not available.
hence, in our modeling approach, we have captured the effect of temporal changes in transcription and protein levels by utilizing the microarray data with a three hour time delay as an input to the model.
this assumes that the corresponding protein profile is similar to the mrna time-course with the three hour delay, i.e., p = g where p and g denote the level of protein and gene-mrna, respectively, at time t which is in hrs.
a three hour time-delay was chosen based on the general consensus on delay in the protein synthesis from mrna  <cit> .
here is an example of model formulation; the enzymes sphk <dig> and sphk <dig> catalyze the conversion of dhsph into dhsph1p .
this effect was functionally captured through the reactions:
dhsph→kf4dhsph1pfor the activity of sphk <dig>
dhsph→kf5dhsph1pfor the activity of sphk <dig>
the effective rate of dhsph1p production was written as ×.
the flux expressions obtained from this scheme were linear in rate parameters and nonlinear in metabolite concentrations.
the matrix-based approach to estimate the rate constants is described below in terms of the reaction numbers labeled in figure  <dig> and listed in table  <dig>  eq.
<dig> describes the rate of change of  and .
ddt=kf4+kf5− kf6ddt=kf7−kf <dig>
microarray data was used for protein concentration.
x → means default degradation of the metabolite x.
unit of first order reaction is 1/hr.
unit of second order reaction is 1/hr when it involves protein as a modifier as we have used fold change data with respect to control for these variables.
unit of third order reaction is μg dna/ when it involves protein as a modifier.
unit of third order reaction is μg dna/ when it involves dg and protein as modifiers.
where the rate constants ki  are as defined in table  <dig>
if the metabolite concentrations were known and the rate parameters were unknown, then the ordinary differential equations  can be rearranged in a matrix format as shown in eq.
<dig>
dtddt]=00−00 0−]ty=  x *   b
the coefficients in the matrix x are linear/nonlinear functions of metabolite and gene concentrations.
all the equations used in the simulation are listed in appendix a. x is completely defined.
the left hand side of the equations  was computed using discretization and the experimental data .
dxdt|t=tk=x|t=tk−x|t=tk−1tk−tk− <dig>
eq.
<dig> contains known matrices x and y, and the only unknown in this equation is the rate-constant vector b. the constrained least-squares approach  was used to solve b. lsqlin optimized the solution with objective function  with additional constraints that all parameter values have to be positive.
min‖y−xb‖2 where y and xb are:dydt|expand  dydt|pred
the estimated values of the parameters were further refined by using generalized constrained nonlinear optimization  where the objective  was to minimize the weighted fit-error between the experimental and predicted metabolite concentrations.
the algorithm of fmincon does not require a matrix form.
thus, numerical integration was used  to simulate the system to circumvent the discretization error.
the combined use of lsqlin and fmincon made the overall process computationally efficient.
the objective function for use with fmincon was:
mink,x0∑i=1nsp)2)k: prameters x0: initial conditions
where, nt is the number of time-points and nsp is the number of species.
the initial concentrations of the metabolites were also optimized in a narrow range around the experimental values.
when data on more than one condition was available, then all the data was used to compute the fit-error by simulating the model several times individually and minimizing the objective function collectively.
estimation of uncertainty in the optimized parameters
the variation among the different technical and biological replicates of lipid and gene data should be accounted for in the estimated values of the parameters.
hence, an uncertainty analysis was performed on the parameters.
their standard error of mean  was computed as follows:
<dig>  compute the sem in the lipid and gene data at each time point .
<dig>  create a candidate data set for parameter estimation by generating nsp x m random matrix utilizing the normal distribution, scale it with the corresponding sem and add the scaled-value matrix to the mean-value data on lipids and genes .
<dig>  estimate the parameters using the candidate data set thus producing a one parameter-value set.
<dig>  repeat steps 2- <dig> k times to generate k parameter-value sets .
<dig>  compute the sem for each parameter across the k sets.
RESULTS
fit to experimental data
odes were generated for all metabolites in the network and effective rate constants were estimated for the simplified sl reaction network using the approach described in the "methods" section.
table  <dig> lists the reactions and the corresponding estimated reaction-rate parameters included in the model.
lipid metabolism and signaling are complex processes and the mechanisms involved are only partially known.
kla treatment generated increases in most sphingolipids in raw cells.
the increase in sphinganine , which doubled in  <dig> hours, followed by increases in downstream metabolites, including n-acyl-sphinganines  after a lag of approximately  <dig> to  <dig> hours indicted the induction of de novo sphingolipid biosynthesis.
to account the effect of kla signaling, the time-courses of the concentrations of dhsph, coa <dig>  c <dig> dg, c <dig> gpcho and microarray data were used as input to the network .
linear interpolations of these time-courses were used in the integration.
the decline in sl control experiment data suggested that the control experiments were not at steady state  <cit> .
the assumption of steady-state was circumvented by including  data obtained in two experimental scenarios, namely, the treatment with kla and the corresponding control data set, during parameter-estimation.
the microarray data was used to represent the corresponding protein concentration with  <dig> hr time delay.
for microarray data, fold change with respect to the control was used in the simulation corresponding to the treatment experiment .
a fold-value of  <dig> was used in the simulation corresponding to the control experiment.
cers <dig> and degs <dig> did not show significant change with respect to the control in microarray data, thus the corresponding reactions  were written without these enzymes.
the omission of these enzymes indicated their constant activity in these reactions.
the simplified model is a reliable predictive model as evidenced by the good fit to experimental data for most metabolites .
for c <dig> dhcer , the fit is good up to  <dig> hrs after which it deteriorates.
one possible reason for this deterioration in the fit is the discrepancy between mrna and protein levels for the enzymes cers5/ <dig>
there were two problems in the optimization of parameters:  the wide range of metabolite concentrations, and  irregular time intervals .
the concentration of metabolites varied between  <dig>  pmol/μg dna  and  <dig> pmol/μg dna .
due to the orders of magnitude difference in the metabolite concentrations, the fits were biased towards metabolites with high concentration and resulted in a poor fit for the metabolites with lower concentrations.
to resolve this issue, the experimental values and predicted values for each metabolite were scaled/normalized by its maximum experimental value.
then, the sum of squares of the normalized fit-error on all metabolites was minimized.
this scaling, essentially, normalized the maximum experimental concentrations to  <dig> for all metabolites and resulted in approximately equal weight for all metabolites.
further, the data measured at irregular time intervals  also created a problem in optimization and led to relatively poor fit at later time points.
from parameter-estimation view-point, the measurements should be made at equal intervals so that equal weight is assigned to the entire time-course.
to account for this, the point-wise-error was scaled by the 1/ <dig> power of the length of the corresponding time-interval, resulting in a higher weight for later time points.
consequently, the quality of fit for the later time-points was improved.
for most time points, the difference between the predicted and experimental data was within the standard-error of mean .
the good fit was obtained for all metabolites under both treatment and control conditions.
parametric sensitivity and time-scale analysis
parametric sensitivity analysis: parametric sensitivity analysis was performed by varying all parameters  by two-fold up and down from its original  value.
the sensitivity of each metabolite was studied by plotting the fold-change at its maximum concentration as compared to the maximum concentration corresponding to the original value of the parameter vs. the ratio-change in the value of the perturbed parameter .
the numerical values of the sensitivity, i.e.
the slope of the plots at the optimized value of the parameter, are listed in additional file  <dig> table s <dig>  for each parameter and each metabolite, monotonic increase, decrease or no change was observed depending upon the respective location of the parameter and the metabolite chosen in the network.
for example, an increase in the parameter kf <dig>   produced a decrease in all upstream metabolites except dhsph1p .
an increase is observed for cer  and its downstream products such as cer-1p, glccer, galcer and sm.
this is meaningful mechanistically because the increase in kf <dig> increases the flux of the reaction c <dig> dhcer → c <dig> cer.
if the level of c <dig> dhcer were not to change much, then one would expect almost proportional increase in c <dig> cer .
however, this is not true because the level of c <dig> dhcer is reduced .
this results in a net sensitivity value of less than  <dig>  similarly, an increase in kf <dig> results in decreased dhcer-1p, dhglccer, and dhsm levels.
these metabolites are products of dhcer and hence the sensitivities are consistent with the structure of the biochemical reaction network .
small to moderate sensitivities  suggest that the biochemical system is robust with respect to parametric perturbations.
uncertainty analysis on the parameters: the sem of the parameters is calculated as described in the methods section and the results are reported in table  <dig>  overall, the parameters are well-behaved and the parameter-estimation procedure is reliable because the fractional sem for the lipid and gene data and the fractional sem for most of the parameters are of the same order .
time-scale analysis: time-scale is an important intrinsic property of dynamical biological systems.
while the time-scale of metabolites, at which they evolve, can be gleaned by analyzing several time-courses under different conditions, a more systemic picture can be obtained through eigenvalue and eigenvector analysis of the jacobian matrix of ordinary differential equations at the steady-state conditions.
in our computations, the steady-state was identified by simulating the system corresponding to the control condition  for a long time .
the jacobian matrix was computed through numerical differentiation of the expressions on the right hand sides of the odes with respect to the state variables.
the eigenvalues were split into three broad ranges.
for each eigenvalue, the metabolites with substantial contribution to the corresponding eigenvector were identified.
depending upon the eigenvalues and metabolites significantly contributing to the corresponding eigenvectors, these metabolites have been divided into three categories as listed in table  <dig>  when a metabolite contributed significantly in two eigenvectors spanning in two different eigen value ranges, the metabolite was assigned to the smaller eigen value because the fast manifold only determines its initial transients and the slow manifold governs the later response leading to steady state.
a comparison of table  <dig> and figure  <dig> shows that the time scale of the metabolite is dependent on its location  and its concentration.
the medium time-scale metabolites  are the leaf node metabolites having low concentration ; the slow time-scale metabolites  are the leaf node metabolite with high concentration .
dhcer and cer have fast time-scale because of their intermediate location in the network, moderate concentration  and high flux through these nodes in the production of sm.
discussion
this study has used the large data sets from mass spectrometric measurements of the sl and the microarray data of the mrnas of raw <dig>  cells generated by the lipid maps consortium to evaluate a model for sl metabolism.
importance of including the transcriptomics data and the data on fatty acyl-coa
in this study, we have included the microarray data for the lipid-related genes using a time-delay of three-hours to account for the time for mrna translation and protein translocation.
the use of mrna data for kinetic modeling was motivated by good correlations between specific genes and its metabolic products in the sphingolipid pathway  <cit> .
in general, protein levels follow the qualitative profile of mrna with appropriate time-delays.
however, we note the potential caveat that the protein levels for some proteins may not be even qualitatively similar to the transcriptional levels of their genes.
the discrepancy of mismatch between the mrna and protein profiles can arise due to several factors such as dependence of mrna translation on ribosome binding site  sequences, post translational modification of the protein, protein translocation and its stability.
this may introduce errors in the estimated values of the kinetic parameters.
to delineate the importance of the transcriptional data, we first developed the mathematical model without using the gene/protein data.
the rate parameters were estimated using the lipidomics data alone.
we visualized the fit to experimental data .
a reasonable fit was obtained for all metabolites except c16-dhcer.
to resolve disagreement in the shape of dhcer, we identified the reactions in which these lipids were consumed or produced.
we also analyzed the time-course of the mrna levels for the genes related to these reactions.
most of the genes exhibited differential regulation at later time points .
among several differentially regulated genes, the prominent were cers <dig> , smpd <dig> , degs <dig>  and ugcg .
the up-regulation and involvement of cers <dig> in the production of dhcer suggested that the gene data must be included in the network to capture the dhcer dynamics.
after adding the gene data with delay, we observed substantial improvement in the fit to the experimental data for dhcer .
the profile of palmitoyl-coa <dig> also increases monotonically during 0- <dig> hrs .
even though  increases till  <dig> hr and comes back to control level at  <dig> hr, the influx to dhcer is maintained throughout 0- <dig> hr because of increase in palmitoyl-coa and cers <dig> at later time points.
due to this reason, dhcer shows a monotonically increasing profile during 0- <dig> hr.
thus, the profiles of cers <dig> and palmitoyl-coa are important to obtain a good fit on dhcer.
rate parameters for the enzymes
we compared the combined rate constant for the enzymes cers5/cers <dig>  with an estimation of the maximum flux through the de novo biosynthetic pathway based on the activity of serine palmitoyltransferase assayed in vitro with optimal concentrations of substrates, which was  <dig> +  <dig> pmole/min/mg protein for control  raw <dig>  cells.
to convert the estimated values to the same units, we used estimations  for the relationships between dna, cell number and protein amount of ~ <dig> μg dna/ <dig> raw <dig>  cells, and  <dig> cells have ~ <dig>  mg protein.
thus, the computed value for c16-sl biosynthesis by cers5/cers <dig> is ~ <dig> pmol/min/mg protein, which is about one order of magnitude lower than the theoretical maximum rate of sphinganine production by the cells.
this might mean that the calculated value is low due to inaccuracies in some of the modeling approximations, including the use of linear kinetics instead of michaelis-menten kinetics; however, the differences might be real because serine palmitoyltransferase has usually not been found to be operating at vmax because its substrates are not saturating  <cit>  and other cers also contribute to the utilization of the sphinganine that is produced de novo.
in addition, sphinganine and sphinganine 1-phosphate are elevated to some extent in raw <dig>  cells, which implies that cerss are not trapping all of the sphinganine that is made.
similarity in the rate parameters for same gene/enzyme involved in different reactions
structurally, dhcer and cer differ by one double bond .
however, both are converted to their corresponding derivatives by four genes, namely ugcg, cerk and sms1/ <dig> .
to check the effects of this double bond on the rate parameter, we compared the rate parameters for the above genes/enzymes in their reactions.
the rate parameter for sms <dig>  and ugcg  was found approximately similar for the both reactions involving dhcer and cer.
for sms <dig> and cerk, the rate parameters for the two reactions differed by a factor of  <dig> .
to further test whether we can find the common rate constants for these enzymes in their reactions, we carried out the optimization using same rate constant for these enzymes in their reactions.
we were able to get reasonably good fit for all of the metabolites .
in the modified optimization, the values obtained for these rate-parameters were between the corresponding values for the two reactions obtained in the original optimization .
this result suggests that the affinities of these enzymes are similar for both the substrates cer and dhcer.
consideration of different chain-lengths
c <dig> sphingolipids are used in this model because  these are present in higher amounts than the sl of larger chain-lengths, and  they showed significant differences between the treatment and control experiments in response to kla.
we can extend our model to include the metabolites with higher chain-lengths .
the reactions in the existing model  will serve as a template that can be instantiated for higher chain-lengths.
to estimate the rate parameters for the entire model, the parameters for the c16-based model can provide bounds for the rate parameters in the overall model for the similar types of reactions.
for example, the rate constants for the reactions cn dhcer → cn dhglccer , can be constrained with a factor of  <dig>  to  <dig> of the rate constant for the reaction c <dig> dhcer → c <dig> dhglccer in the model developed here.
use of the model for in-silico perturbation experiments
the model can be used to perform knockdown  or other perturbation experiments in-silico.
such computational studies provide useful insights into the behavior of the system which can add in designing the actual perturbation experiments.
these simulations can assist in finding the propagation of the effects of kd in the network and lead to a better design of experiments .
simulation can also suggest whether a high level of kd is needed such that the differences between the kd and control scenarios would be statistically significant.
the data generated from new perturbation experiments can be used to further refine the model.
for the enzymes, the kd perturbations can be simulated by decreasing the corresponding rate parameter because amount of the  enzyme directly affects the rate parameters.
the kd perturbation simulation results will be similar to the sensitivity analysis results on the corresponding rate parameter.
for example, the effect of 50% kd of the enzyme for the reaction  <dig> , can be predicted from the simulation for sensitivity analysis in which the value of the parameter kf <dig> is reduced by 50% .
the case of knockdown of cerk  is also interesting due to its direct effect on ceramide-1-phosphate  which inhibits apoptosis and induces cell growth  <cit> .
the result of sensitivity analysis for the corresponding rate parameter, kf <dig> , is shown in additional file  <dig> figure s <dig>  decrease in kf <dig> results in a corresponding decrease in the cerp level.
these changes have also been observed in recent experiments validating our parametric sensitivity analysis and the predictive ability of the model  <cit> .
cerp has also been implicated in regulating the homeostasis of calcium  <cit>  thereby affecting the activity of several signaling pathways.
cerp and cerk mediate the effect of cytokines to activate cytosolic phospholipase a <dig>  and cyclooxygenase  <dig> , resulting in increased production of prostaglandin e <dig> , a mediator of inflammation  <cit> .
hence, it has been hypothesized that cerk could be a potential target for anti-inflammatory drugs  <cit> .
CONCLUSIONS
use of systems biology approaches is becoming more common in the study of lipids to elucidate their functions and roles in human health and diseases such as arthritis and cancer.
systems biology has already been recognized as an indispensable tool in pathway-based drug discovery.
here we have applied a matrix-based approach to develop a dynamic model of sl metabolism by integrating legacy information on the lipid pathways with novel experimental data.
the metabolic pathway was reconstructed using information from the kegg database and the existing literature.
based upon the network map reconstructed, we have developed an ordinary differential equations-based mathematical model.
parameter-estimation used a two-step approach.
in the first step, a matrix-based approach provided an initial guess.
the parameter-values were further refined in the second step.
the resulting model fitted the experimental data well for all species and demonstrated that the integrated metabolic and signaling network and the experimental data are consistent with each other.
the robustness of the model parameters was also validated through sensitivity analysis.
though we have used this two-step approach previously and applied it to eicosanoid lipid pathway, the major distinction lies in its application to the sl metabolic pathway and the integration of transcriptomic data with the metabolomic data along with legacy knowledge to develop the kinetic model.
previous computational models of sphingolipid metabolism were for non-mammalian systems in which only a few metabolites were measured as compared to the total number of metabolites in the reaction network.
in comparison, our model is based on a large amount of experimental time-course data where the concentrations of most metabolites and mrna levels of genes in the network are measured.
this provides a more context-specific model for raw cells in particular and for mammalian cell systems in general.
abbreviations used
cer: ceramide; cerp: ceramide-1-phosphate; cerk: ceramide kinase; degs1: degenerative spermatocyte homolog  <dig> ; degs2: degenerative spermatocyte homolog  <dig> , lipid desaturase; cers5: ceramide synthase 5; cers6: ceramide synthase 6; smpd1: sphingomyelin phosphodiesterase  <dig>  acid lysosomal; sphk1: sphingosine kinase 1; sphk2: sphingosine kinase 2; sms1: sphingomyelin synthase 1; sms2: sphingomyelin synthase 2; ugcg: udp-glucose ceramide glucosyltransferase; c <dig> cer: n-- sphing-4-enine ; c <dig> cerp: n-- sphing-4-enine-1-phosphate ; c <dig> glccer: n--1-β- sphing-4-enine ; c <dig> sm: n-- sphing-4-enine -1-phosphocholine ; c <dig> dhcer: n--sphinganine ; c <dig> dhcerp: n--sphinganine-1-phosphate; c <dig> dhglccer: n--1-β-glucosyl-sphinganine; c <dig> dhsm: n--sphinganine-1-phosphocholine; dhsph: sphinganine; dhsph1p: sphinganine-1-phosphate; dhgalcer: dihydro galactosylceramide; galcer: galactosylceramide;
authors' contributions
sg designed the simulations, wrote the computer program, analyzed the experimental data and the simulation results and wrote the manuscript.
ahm assisted in refining the reaction network.
mrm assisted in designing some of the simulations, wrote part of the computer program and contributed in writing the manuscript.
ahm, ckg and ss assisted in revising the manuscript.
sphingolipid and enzyme activity measurements were carried out in the laboratory of ahm.
transcriptomics experiments were carried out in the laboratory of ckg.
ss supervised the modeling study.
all authors have read and approved the final manuscript.
appendix a
the flux expressions for the reactions shown in figure  <dig> are as follows:
v <dig> = kf1
v <dig> = kf2
v <dig> = kf3
v <dig> = kf4
v <dig> = kf5
v <dig> = kf6
v <dig> = kf7
v <dig> = kf8
v <dig> = kf <dig> - kb9
v <dig> = kf <dig> - kb10
v <dig> = kf11
v <dig> = kf12
v <dig> = kf13
v <dig> = kf14
v <dig> = kf15
v <dig> = kf16
v <dig> = kf17
v <dig> = kf18
v <dig> = kf <dig> - kb19
v <dig> = kf <dig> - kb20
v <dig> = kf21
v <dig> = kf22
v <dig> = kf23
v <dig> = kf24
v <dig> = kf25
the differential equations describing the rate of change of metabolite concentrations are:
d/dt = v <dig> + v <dig> - v <dig> - v <dig> - v <dig> - v <dig> + v <dig> - v <dig> - v <dig> - v16
d/dt = v <dig> + v <dig> - v6
d/dt = v <dig> - v8
d/dt = v <dig> + v <dig> - v <dig> - v12
d/dt = v <dig> - v14
d/dt = v <dig> + v <dig> - v <dig> - v <dig> - v <dig> + v <dig> - v <dig> - v24
d/dt = v <dig> - v18
d/dt = v <dig> + v <dig> - v <dig> - v22
d/dt = v <dig> - v25
supplementary material
additional file 1
supporting material.
this file is in pdf format and contains additional figures s1-s <dig> and additional tables s1-s <dig>
click here for file
acknowledgements
we would like to acknowledge the assistance of eoin fahy for preparing the network map of sphingolipid metabolism, and elaine wang, kristin jones, samuel kelly, rebecca shaner and jeremy allegood for the analysis of the sphingolipids by mass spectrometry.
this research was supported by nih collaborative grant u <dig> gm69338- <dig> lipid maps  and niddk grant p01-dk <dig> .

 sphingolipids are categorized as lipids with a sphingoid base backbone that is often derivatized with an amide-linked fatty acid to make ceramides and more structurally complex sl with diverse biological functions. sl in essentially every subcategory, from the lipid backbones to complex sl, are highly bioactive and play important roles in diseases; hence, methods for "lipidomic" analysis of sl and sl metabolism are important for an in-depth understanding of these enigmatic compounds. in recent years, a number of large-scale experimental and bioinformatics projects have begun to address the complexity of the lipidome. examples include the lipid metabolites and pathways strategy consortium, the lipid library, cyberlipid center and lipidbank. in particular, lipid maps has become a comprehensive resource for information on classification, structures and quantitative data on lipids and provides an opportunity for developing quantitative models of lipid synthesis and metabolism thus facilitating a mechanistic and systems-level understanding. the de novo biosynthesis of sl begins with production of the sphingoid base, which utilizes serine and palmitoyl-coenzyme a and various fatty acyl-coas to make n-acylsphinganines that are desaturated to cer and incorporated into more complex sl such as ceramide 1-phosphate, sphingomyelin, glucosyl- and galactosyl-ceramide and more complex glycosphingolipids. ceramide can also be synthesized by recycling of sphingosine from turnover of sl such as sm; furthermore, sphingosine and sphinganine are phosphorylated to sphingosine 1-phosphate and sphinganine 1-phosphate which are intermediates of sphingoid base degradation and cell signaling molecules. due to the complexity of this pathway, and the paucity of data for its many metabolites, there are only a few models of sl metabolism available in the literature. the lipid maps consortium has quantified the global changes in lipid metabolites and genes in raw macrophage cells treated with kdo2-lipid a. kla is the active component of inflammatory lipopolysaccharide which functions as a selective agonist of toll-like receptor and mimics bacterial infection. the measurements are carried out over a 24-hour time period and the data is freely available via the lipid maps website. the goal of the work presented here is to develop a predictive kinetic model for sl metabolism using the lipidomics and transcriptomics data from the lipid maps studies. this manuscript is organized as follows: we first briefly discuss the experimental data preprocessing and the methodology used for estimating the rate parameters, then we present the results of parameter estimation, followed by discussion and conclusions. methods network simplification a detailed metabolic reaction network was developed using the information available in the literature and the kegg pathways database. the c16-branch of cer biosynthesis was selected for developing the model because this is a major subspecies for all categories of complex sl in the raw cells. vanted software was used to draw the reaction network. it is common in modeling studies for the network to contain several unmeasured nodes; in our pathway, quantities are known for all of the metabolites and genes except dhgalcer and galcer and are expressed as leaf-nodes in the network. one of the steps in our matrix-based fast algorithm for parameter-estimation requires experimental data on all metabolites except on the leaf-nodes/metabolites in the network. a detailed procedure for simplifying the network, if the network contains unmeasured components, is described in our previous work. the leaf-nodes were exempted from the model described in this paper because the reactions leading to the unmeasured leaf metabolites were combined with the default degradation of their precursors. experimental data and preprocessing the lipid maps consortium has measured all the major lipids in mouse macrophage raw cells grown in 10% serum and treated with kla. all metabolites were measured in pmol/μg-dna units. the time-courses of lipids and genes were measured under two conditions: normal condition designated as control and inflamed condition. time series comprising of data points at, and hr were measured with three biological replicates consisting of three technical replicates each. the three technical replicate experiments were performed on the same day with a single batch of cells. in addition, each time course was repeated three times on different days with different batches of cells. for kinetic modeling purposes, outlier points were detected by a t-test and were excluded at each time-point. the resulting data from all the replicates were averaged at each time-point. data was processed for all metabolites under the treatment with kla and control conditions. serine palmitoyltransferase activity in vitro was determined using the assay conditions described previously but with palmitoyl-coa as the labeled substrate. after extraction of the products, the amounts of 3-ketosphinganine were determined by liquid chromatography, electrospray-ionization tandem mass spectrometry. the assay was conducted using μg of protein obtained by sonication of raw cells in the buffer used for the assay. development of a kinetic model and parameter estimation we have developed a kinetic model of the sl metabolism. the procedure is similar to the one described in our previous work, but it is presented here through an example of reactions from the sl network for completeness. the reaction rates were described by linear or law of mass action kinetics with the assumption that for enzymatic reactions, the substrate concentrations are much smaller as compared to the corresponding michaelis constant, km. for example, the following types of reaction schemes and rate expressions were used: reactionrate expressionsa→bk*a+b→a.bk**a→bk** the enzymes involved in sl metabolism can be regulated at multiple levels from mrna expression to posttranslational modification. mrna data on the genes involved in the pathway is available from microarray experiments. however, the corresponding proteomics data is not available. hence, in our modeling approach, we have captured the effect of temporal changes in transcription and protein levels by utilizing the microarray data with a three hour time delay as an input to the model. this assumes that the corresponding protein profile is similar to the mrna time-course with the three hour delay, i.e., p = g where p and g denote the level of protein and gene-mrna, respectively, at time t which is in hrs. a three hour time-delay was chosen based on the general consensus on delay in the protein synthesis from mrna. here is an example of model formulation; the enzymes sphk and sphk catalyze the conversion of dhsph into dhsph1p. this effect was functionally captured through the reactions:  dhsph→kf4dhsph1pfor the activity of sphk  dhsph→kf5dhsph1pfor the activity of sphk the effective rate of dhsph1p production was written as ×. the flux expressions obtained from this scheme were linear in rate parameters and nonlinear in metabolite concentrations. the matrix-based approach to estimate the rate constants is described below in terms of the reaction numbers labeled in figure and listed in table eq. describes the rate of change of and. ddt=kf4+kf5− kf6ddt=kf7−kf microarray data was used for protein concentration. x → means default degradation of the metabolite x. unit of first order reaction is 1/hr. unit of second order reaction is 1/hr when it involves protein as a modifier as we have used fold change data with respect to control for these variables. unit of third order reaction is μg dna/ when it involves protein as a modifier. unit of third order reaction is μg dna/ when it involves dg and protein as modifiers. where the rate constants ki are as defined in table if the metabolite concentrations were known and the rate parameters were unknown, then the ordinary differential equations can be rearranged in a matrix format as shown in eq. dtddt]=00−00 0−]ty=  x *   b the coefficients in the matrix x are linear/nonlinear functions of metabolite and gene concentrations. all the equations used in the simulation are listed in appendix a. x is completely defined. the left hand side of the equations was computed using discretization and the experimental data. dxdt|t=tk=x|t=tk−x|t=tk−1tk−tk− eq. contains known matrices x and y, and the only unknown in this equation is the rate-constant vector b. the constrained least-squares approach was used to solve b. lsqlin optimized the solution with objective function with additional constraints that all parameter values have to be positive. min‖y−xb‖2 where y and xb are:dydt|expand  dydt|pred the estimated values of the parameters were further refined by using generalized constrained nonlinear optimization where the objective was to minimize the weighted fit-error between the experimental and predicted metabolite concentrations. the algorithm of fmincon does not require a matrix form. thus, numerical integration was used to simulate the system to circumvent the discretization error. the combined use of lsqlin and fmincon made the overall process computationally efficient. the objective function for use with fmincon was: mink,x0∑i=1nsp)2)k: prameters x0: initial conditions  where, nt is the number of time-points and nsp is the number of species. the initial concentrations of the metabolites were also optimized in a narrow range around the experimental values. when data on more than one condition was available, then all the data was used to compute the fit-error by simulating the model several times individually and minimizing the objective function collectively. estimation of uncertainty in the optimized parameters the variation among the different technical and biological replicates of lipid and gene data should be accounted for in the estimated values of the parameters. hence, an uncertainty analysis was performed on the parameters. their standard error of mean was computed as follows: compute the sem in the lipid and gene data at each time point. create a candidate data set for parameter estimation by generating nsp x m random matrix utilizing the normal distribution, scale it with the corresponding sem and add the scaled-value matrix to the mean-value data on lipids and genes. estimate the parameters using the candidate data set thus producing a one parameter-value set. repeat steps 2- k times to generate k parameter-value sets. compute the sem for each parameter across the k sets.  fit to experimental data odes were generated for all metabolites in the network and effective rate constants were estimated for the simplified sl reaction network using the approach described in the "methods" section. table lists the reactions and the corresponding estimated reaction-rate parameters included in the model. lipid metabolism and signaling are complex processes and the mechanisms involved are only partially known. kla treatment generated increases in most sphingolipids in raw cells. the increase in sphinganine, which doubled in hours, followed by increases in downstream metabolites, including n-acyl-sphinganines after a lag of approximately to hours indicted the induction of de novo sphingolipid biosynthesis. to account the effect of kla signaling, the time-courses of the concentrations of dhsph, coa c dg, c gpcho and microarray data were used as input to the network. linear interpolations of these time-courses were used in the integration. the decline in sl control experiment data suggested that the control experiments were not at steady state. the assumption of steady-state was circumvented by including data obtained in two experimental scenarios, namely, the treatment with kla and the corresponding control data set, during parameter-estimation. the microarray data was used to represent the corresponding protein concentration with hr time delay. for microarray data, fold change with respect to the control was used in the simulation corresponding to the treatment experiment. a fold-value of was used in the simulation corresponding to the control experiment. cers and degs did not show significant change with respect to the control in microarray data, thus the corresponding reactions were written without these enzymes. the omission of these enzymes indicated their constant activity in these reactions. the simplified model is a reliable predictive model as evidenced by the good fit to experimental data for most metabolites. for c dhcer, the fit is good up to hrs after which it deteriorates. one possible reason for this deterioration in the fit is the discrepancy between mrna and protein levels for the enzymes cers5/ there were two problems in the optimization of parameters: the wide range of metabolite concentrations, and irregular time intervals. the concentration of metabolites varied between pmol/μg dna and pmol/μg dna. due to the orders of magnitude difference in the metabolite concentrations, the fits were biased towards metabolites with high concentration and resulted in a poor fit for the metabolites with lower concentrations. to resolve this issue, the experimental values and predicted values for each metabolite were scaled/normalized by its maximum experimental value. then, the sum of squares of the normalized fit-error on all metabolites was minimized. this scaling, essentially, normalized the maximum experimental concentrations to for all metabolites and resulted in approximately equal weight for all metabolites. further, the data measured at irregular time intervals also created a problem in optimization and led to relatively poor fit at later time points. from parameter-estimation view-point, the measurements should be made at equal intervals so that equal weight is assigned to the entire time-course. to account for this, the point-wise-error was scaled by the 1/ power of the length of the corresponding time-interval, resulting in a higher weight for later time points. consequently, the quality of fit for the later time-points was improved. for most time points, the difference between the predicted and experimental data was within the standard-error of mean. the good fit was obtained for all metabolites under both treatment and control conditions. parametric sensitivity and time-scale analysis parametric sensitivity analysis: parametric sensitivity analysis was performed by varying all parameters by two-fold up and down from its original value. the sensitivity of each metabolite was studied by plotting the fold-change at its maximum concentration as compared to the maximum concentration corresponding to the original value of the parameter vs. the ratio-change in the value of the perturbed parameter. the numerical values of the sensitivity, i.e. the slope of the plots at the optimized value of the parameter, are listed in additional file table s for each parameter and each metabolite, monotonic increase, decrease or no change was observed depending upon the respective location of the parameter and the metabolite chosen in the network. for example, an increase in the parameter kf produced a decrease in all upstream metabolites except dhsph1p. an increase is observed for cer and its downstream products such as cer-1p, glccer, galcer and sm. this is meaningful mechanistically because the increase in kf increases the flux of the reaction c dhcer → c cer. if the level of c dhcer were not to change much, then one would expect almost proportional increase in c cer. however, this is not true because the level of c dhcer is reduced. this results in a net sensitivity value of less than similarly, an increase in kf results in decreased dhcer-1p, dhglccer, and dhsm levels. these metabolites are products of dhcer and hence the sensitivities are consistent with the structure of the biochemical reaction network. small to moderate sensitivities suggest that the biochemical system is robust with respect to parametric perturbations. uncertainty analysis on the parameters: the sem of the parameters is calculated as described in the methods section and the results are reported in table overall, the parameters are well-behaved and the parameter-estimation procedure is reliable because the fractional sem for the lipid and gene data and the fractional sem for most of the parameters are of the same order. time-scale analysis: time-scale is an important intrinsic property of dynamical biological systems. while the time-scale of metabolites, at which they evolve, can be gleaned by analyzing several time-courses under different conditions, a more systemic picture can be obtained through eigenvalue and eigenvector analysis of the jacobian matrix of ordinary differential equations at the steady-state conditions. in our computations, the steady-state was identified by simulating the system corresponding to the control condition for a long time. the jacobian matrix was computed through numerical differentiation of the expressions on the right hand sides of the odes with respect to the state variables. the eigenvalues were split into three broad ranges. for each eigenvalue, the metabolites with substantial contribution to the corresponding eigenvector were identified. depending upon the eigenvalues and metabolites significantly contributing to the corresponding eigenvectors, these metabolites have been divided into three categories as listed in table when a metabolite contributed significantly in two eigenvectors spanning in two different eigen value ranges, the metabolite was assigned to the smaller eigen value because the fast manifold only determines its initial transients and the slow manifold governs the later response leading to steady state. a comparison of table and figure shows that the time scale of the metabolite is dependent on its location and its concentration. the medium time-scale metabolites are the leaf node metabolites having low concentration; the slow time-scale metabolites are the leaf node metabolite with high concentration. dhcer and cer have fast time-scale because of their intermediate location in the network, moderate concentration and high flux through these nodes in the production of sm. discussion this study has used the large data sets from mass spectrometric measurements of the sl and the microarray data of the mrnas of raw cells generated by the lipid maps consortium to evaluate a model for sl metabolism. importance of including the transcriptomics data and the data on fatty acyl-coa in this study, we have included the microarray data for the lipid-related genes using a time-delay of three-hours to account for the time for mrna translation and protein translocation. the use of mrna data for kinetic modeling was motivated by good correlations between specific genes and its metabolic products in the sphingolipid pathway. in general, protein levels follow the qualitative profile of mrna with appropriate time-delays. however, we note the potential caveat that the protein levels for some proteins may not be even qualitatively similar to the transcriptional levels of their genes. the discrepancy of mismatch between the mrna and protein profiles can arise due to several factors such as dependence of mrna translation on ribosome binding site sequences, post translational modification of the protein, protein translocation and its stability. this may introduce errors in the estimated values of the kinetic parameters. to delineate the importance of the transcriptional data, we first developed the mathematical model without using the gene/protein data. the rate parameters were estimated using the lipidomics data alone. we visualized the fit to experimental data. a reasonable fit was obtained for all metabolites except c16-dhcer. to resolve disagreement in the shape of dhcer, we identified the reactions in which these lipids were consumed or produced. we also analyzed the time-course of the mrna levels for the genes related to these reactions. most of the genes exhibited differential regulation at later time points. among several differentially regulated genes, the prominent were cers, smpd, degs and ugcg. the up-regulation and involvement of cers in the production of dhcer suggested that the gene data must be included in the network to capture the dhcer dynamics. after adding the gene data with delay, we observed substantial improvement in the fit to the experimental data for dhcer. the profile of palmitoyl-coa also increases monotonically during 0- hrs. even though increases till hr and comes back to control level at hr, the influx to dhcer is maintained throughout 0- hr because of increase in palmitoyl-coa and cers at later time points. due to this reason, dhcer shows a monotonically increasing profile during 0- hr. thus, the profiles of cers and palmitoyl-coa are important to obtain a good fit on dhcer. rate parameters for the enzymes we compared the combined rate constant for the enzymes cers5/cers with an estimation of the maximum flux through the de novo biosynthetic pathway based on the activity of serine palmitoyltransferase assayed in vitro with optimal concentrations of substrates, which was + pmole/min/mg protein for control raw cells. to convert the estimated values to the same units, we used estimations for the relationships between dna, cell number and protein amount of ~ μg dna/ raw cells, and cells have ~ mg protein. thus, the computed value for c16-sl biosynthesis by cers5/cers is ~ pmol/min/mg protein, which is about one order of magnitude lower than the theoretical maximum rate of sphinganine production by the cells. this might mean that the calculated value is low due to inaccuracies in some of the modeling approximations, including the use of linear kinetics instead of michaelis-menten kinetics; however, the differences might be real because serine palmitoyltransferase has usually not been found to be operating at vmax because its substrates are not saturating and other cers also contribute to the utilization of the sphinganine that is produced de novo. in addition, sphinganine and sphinganine 1-phosphate are elevated to some extent in raw cells, which implies that cerss are not trapping all of the sphinganine that is made. similarity in the rate parameters for same gene/enzyme involved in different reactions structurally, dhcer and cer differ by one double bond. however, both are converted to their corresponding derivatives by four genes, namely ugcg, cerk and sms1/. to check the effects of this double bond on the rate parameter, we compared the rate parameters for the above genes/enzymes in their reactions. the rate parameter for sms and ugcg was found approximately similar for the both reactions involving dhcer and cer. for sms and cerk, the rate parameters for the two reactions differed by a factor of. to further test whether we can find the common rate constants for these enzymes in their reactions, we carried out the optimization using same rate constant for these enzymes in their reactions. we were able to get reasonably good fit for all of the metabolites. in the modified optimization, the values obtained for these rate-parameters were between the corresponding values for the two reactions obtained in the original optimization. this result suggests that the affinities of these enzymes are similar for both the substrates cer and dhcer. consideration of different chain-lengths c sphingolipids are used in this model because these are present in higher amounts than the sl of larger chain-lengths, and they showed significant differences between the treatment and control experiments in response to kla. we can extend our model to include the metabolites with higher chain-lengths. the reactions in the existing model will serve as a template that can be instantiated for higher chain-lengths. to estimate the rate parameters for the entire model, the parameters for the c16-based model can provide bounds for the rate parameters in the overall model for the similar types of reactions. for example, the rate constants for the reactions cn dhcer → cn dhglccer, can be constrained with a factor of to of the rate constant for the reaction c dhcer → c dhglccer in the model developed here. use of the model for in-silico perturbation experiments the model can be used to perform knockdown or other perturbation experiments in-silico. such computational studies provide useful insights into the behavior of the system which can add in designing the actual perturbation experiments. these simulations can assist in finding the propagation of the effects of kd in the network and lead to a better design of experiments. simulation can also suggest whether a high level of kd is needed such that the differences between the kd and control scenarios would be statistically significant. the data generated from new perturbation experiments can be used to further refine the model. for the enzymes, the kd perturbations can be simulated by decreasing the corresponding rate parameter because amount of the enzyme directly affects the rate parameters. the kd perturbation simulation results will be similar to the sensitivity analysis results on the corresponding rate parameter. for example, the effect of 50% kd of the enzyme for the reaction, can be predicted from the simulation for sensitivity analysis in which the value of the parameter kf is reduced by 50%. the case of knockdown of cerk is also interesting due to its direct effect on ceramide-1-phosphate which inhibits apoptosis and induces cell growth. the result of sensitivity analysis for the corresponding rate parameter, kf, is shown in additional file figure s decrease in kf results in a corresponding decrease in the cerp level. these changes have also been observed in recent experiments validating our parametric sensitivity analysis and the predictive ability of the model. cerp has also been implicated in regulating the homeostasis of calcium thereby affecting the activity of several signaling pathways. cerp and cerk mediate the effect of cytokines to activate cytosolic phospholipase a and cyclooxygenase, resulting in increased production of prostaglandin e, a mediator of inflammation. hence, it has been hypothesized that cerk could be a potential target for anti-inflammatory drugs.  use of systems biology approaches is becoming more common in the study of lipids to elucidate their functions and roles in human health and diseases such as arthritis and cancer. systems biology has already been recognized as an indispensable tool in pathway-based drug discovery. here we have applied a matrix-based approach to develop a dynamic model of sl metabolism by integrating legacy information on the lipid pathways with novel experimental data. the metabolic pathway was reconstructed using information from the kegg database and the existing literature. based upon the network map reconstructed, we have developed an ordinary differential equations-based mathematical model. parameter-estimation used a two-step approach. in the first step, a matrix-based approach provided an initial guess. the parameter-values were further refined in the second step. the resulting model fitted the experimental data well for all species and demonstrated that the integrated metabolic and signaling network and the experimental data are consistent with each other. the robustness of the model parameters was also validated through sensitivity analysis. though we have used this two-step approach previously and applied it to eicosanoid lipid pathway, the major distinction lies in its application to the sl metabolic pathway and the integration of transcriptomic data with the metabolomic data along with legacy knowledge to develop the kinetic model. previous computational models of sphingolipid metabolism were for non-mammalian systems in which only a few metabolites were measured as compared to the total number of metabolites in the reaction network. in comparison, our model is based on a large amount of experimental time-course data where the concentrations of most metabolites and mrna levels of genes in the network are measured. this provides a more context-specific model for raw cells in particular and for mammalian cell systems in general. abbreviations used cer: ceramide; cerp: ceramide-1-phosphate; cerk: ceramide kinase; degs1: degenerative spermatocyte homolog; degs2: degenerative spermatocyte homolog, lipid desaturase; cers5: ceramide synthase 5; cers6: ceramide synthase 6; smpd1: sphingomyelin phosphodiesterase acid lysosomal; sphk1: sphingosine kinase 1; sphk2: sphingosine kinase 2; sms1: sphingomyelin synthase 1; sms2: sphingomyelin synthase 2; ugcg: udp-glucose ceramide glucosyltransferase; c cer: n-- sphing-4-enine; c cerp: n-- sphing-4-enine-1-phosphate; c glccer: n--1-β- sphing-4-enine; c sm: n-- sphing-4-enine -1-phosphocholine; c dhcer: n--sphinganine; c dhcerp: n--sphinganine-1-phosphate; c dhglccer: n--1-β-glucosyl-sphinganine; c dhsm: n--sphinganine-1-phosphocholine; dhsph: sphinganine; dhsph1p: sphinganine-1-phosphate; dhgalcer: dihydro galactosylceramide; galcer: galactosylceramide; authors' contributions sg designed the simulations, wrote the computer program, analyzed the experimental data and the simulation results and wrote the manuscript. ahm assisted in refining the reaction network. mrm assisted in designing some of the simulations, wrote part of the computer program and contributed in writing the manuscript. ahm, ckg and ss assisted in revising the manuscript. sphingolipid and enzyme activity measurements were carried out in the laboratory of ahm. transcriptomics experiments were carried out in the laboratory of ckg. ss supervised the modeling study. all authors have read and approved the final manuscript. appendix a the flux expressions for the reactions shown in figure are as follows: v = kf1 v = kf2 v = kf3 v = kf4 v = kf5 v = kf6 v = kf7 v = kf8 v = kf - kb9 v = kf - kb10 v = kf11 v = kf12 v = kf13 v = kf14 v = kf15 v = kf16 v = kf17 v = kf18 v = kf - kb19 v = kf - kb20 v = kf21 v = kf22 v = kf23 v = kf24 v = kf25 the differential equations describing the rate of change of metabolite concentrations are: d/dt = v + v - v - v - v - v + v - v - v - v16 d/dt = v + v - v6 d/dt = v - v8 d/dt = v + v - v - v12 d/dt = v - v14 d/dt = v + v - v - v - v + v - v - v24 d/dt = v - v18 d/dt = v + v - v - v22 d/dt = v - v25 supplementary material additional file 1 supporting material. this file is in pdf format and contains additional figures s1-s and additional tables s1-s click here for file acknowledgements we would like to acknowledge the assistance of eoin fahy for preparing the network map of sphingolipid metabolism, and elaine wang, kristin jones, samuel kelly, rebecca shaner and jeremy allegood for the analysis of the sphingolipids by mass spectrometry. this research was supported by nih collaborative grant u gm69338- lipid maps and niddk grant p01-dk.
a quantitative model of the sphigolipid pathway is developed by integrating metabolomics and transcriptomics data with legacy knowledge. in this direction, the lipid maps consortium is developing methods to quantitate the sphingolipid metabolites in mammalian cells and is investigating their application to studies of the activation of the raw macrophage cell by a chemically defined endotoxin, kdo2-lipid a. herein, we describe a model for the c16-branch of sphingolipid metabolism integrating lipidomics and transcriptomics data and using a two-step matrix-based approach to estimate the rate constants from experimental data. the rate constants obtained from the first step are further refined using generalized constrained nonlinear optimization. sphingolipids play important roles in cell structure and function as well as in the pathophysiology of many diseases. the resulting model fits the experimental data for all species. the robustness of the model is validated through parametric sensitivity analysis.
 sphingolipids play important roles in cell structure and function as well as in the pathophysiology of many diseases. many of the intermediates of sphingolipid biosynthesis are highly bioactive and sometimes have antagonistic activities, for example, ceramide promotes apoptosis whereas sphingosine-1-phosphate can inhibit apoptosis and induce cell growth; therefore, quantification of the metabolites and modeling of the sphingolipid network is imperative for an understanding of sphingolipid biology.  in this direction, the lipid maps consortium is developing methods to quantitate the sphingolipid metabolites in mammalian cells and is investigating their application to studies of the activation of the raw macrophage cell by a chemically defined endotoxin, kdo2-lipid a. herein, we describe a model for the c16-branch of sphingolipid metabolism integrating lipidomics and transcriptomics data and using a two-step matrix-based approach to estimate the rate constants from experimental data. the rate constants obtained from the first step are further refined using generalized constrained nonlinear optimization. the resulting model fits the experimental data for all species. the robustness of the model is validated through parametric sensitivity analysis.  a quantitative model of the sphigolipid pathway is developed by integrating metabolomics and transcriptomics data with legacy knowledge. the model could be used to design experimental studies of how genetic and pharmacological perturbations alter the flux through this important lipid biosynthetic pathway. 