Abstract
Background
The study of phenotype transitions is important to understand progressive diseases, e.g., diabetes mellitus, metabolic syndrome, and cardiovascular diseases. A challenge remains to explain phenotype transitions in terms of adaptations in molecular components and interactions in underlying biological systems.
Results
Here, mathematical modeling is used to describe the different phenotypes by integrating experimental data on metabolic pools and fluxes. Subsequently, trajectories of parameter adaptations are identified that are essential for the phenotypical changes. These changes in parameters reflect progressive adaptations at the transcriptome and proteome level, which occur at larger timescales. The approach was employed to study the metabolic processes underlying liver X receptor induced hepatic steatosis. Model analysis predicts which molecular processes adapt in time after pharmacological activation of the liver X receptor. Our results show that hepatic triglyceride fluxes are increased and triglycerides are especially stored in cytosolic fractions, rather than in endoplasmic reticulum fractions. Furthermore, the model reveals several possible scenarios for adaptations in cholesterol metabolism. According to the analysis, the additional quantification of one cholesterol flux is sufficient to exclude many of these hypotheses.
Conclusions
We propose a generic computational approach to analyze biological systems evolving through various phenotypes and to predict which molecular processes are responsible for the transition. For the case of liver X receptor induced hepatic steatosis the novel approach yields information about the redistribution of fluxes and pools of triglycerides and cholesterols that was not directly apparent from the experimental data. Model analysis provides guidance which specific molecular processes to study in more detail to obtain further understanding of the underlying biological system.
Background
Cardiovascular and metabolic diseases such as diabetes mellitus and metabolic syndrome are progressive in time [15]. Progressive diseases are often being studied by experimentally comparing different states: a control state representing a healthy phenotype, and one or more adapted states representing phenotypes of certain stages of the disease. Experimentally observed differences between phenotypes provide information about biological processes that are involved in the pathogenesis. Most research is carried out using mouse models, having many practical advantages such as short generation times, reduced genetic variation, and the possibility to apply gene manipulation technology [68]. For instance, the genetic leptindeficient (ob/ob) or leptinresistant (db/db) mouse are frequently used to study metabolic pathologies, e.g., obesity, insulin resistance, and diabetes [912]. A challenging task is to explain phenotypical characteristics and the progression of phenotype transitions in terms of adaptations in molecular components and interactions in underlying biological systems. This is especially the case for the study of progressive diseases in which multiple processes, operating on various length and timescales, are altered.
In systems biology mathematical modeling is applied to integrate different sources of experimental data of a phenotype and to investigate the complex interactions of underlying biological systems [1319]. However, several issues complicate the simulation and prediction of molecular adaptations associated with progressive diseases. One problem is to cover large differences in timescales. Computational models in molecular systems biology are typically constructed to simulate processes on a single timescale. These range from seconds in signal transduction and metabolic network models to hours for genetic networks [2024]. On the other hand, progressive diseases often comprise of a combination of these processes and typically develop over a time span of years in humans. Another issue is that mathematically describing progressive adaptations could become unfeasible when sufficient information of the underlying biological system, such as network structure, molecular concentrations and fluxes, as well as their interaction mechanisms, is lacking.
In the present work, we propose a novel computational approach to analyze molecular adaptations in a biological system to overcome these problems. We use mathematical modeling to quantitatively integrate metabolic data of different phenotypes and subsequently exploit this mathematical framework to analyze which molecular processes have changed and are collectively responsible for the shift between phenotypes. This information is obtained by identifying the progression of necessary parameter changes required for the model to be consistent with the experimental data of these phenotypes. These changes in parameters reflect progressive adaptations at the transcriptome and proteome level, which occur at larger timescales than the metabolic processes. The approach involves consecutive steps of data gathering, model development, and parameter estimation, which will be discussed in detail. An advantage of our approach is that mathematical models containing processes at any timescale of interest can be used, while their longterm adaptations are captured by identifying necessary parameter changes. This enables us to study longterm aspects of shortterm processes. Furthermore, in cases when the amount of information of the underlying biological system is limited, our approach could provide a means to describe adaptations in molecular processes without the necessity to develop detailed kinetic models of the modulating mechanisms. For instance, if one is interested in studying a metabolic pathway which is adapting due to activation of a signal transduction pathway, the modulating effects can be captured by identifying necessary changes in the metabolic pathway parameters rather than developing a mathematical model that includes an explicitly modeled signal transduction pathway. The approach, which is applicable to a multitude of biological systems, is demonstrated on the basis of a case involving the activation of the liver X receptor (LXR), a promising drug target for atherosclerotic therapies [25,26].
The family of liver X receptors (LXRα and LXRβ) is involved in the control of cellular lipid metabolism. LXRs, when ligandactivated by oxysterols, heterodimerize with the retinoid X receptor (RXR) and bind to LXR responsive elements on the DNA [27], where they induce the transcription of lipogenic genes such as SREBP1C, FAS, ABCA1, and ACC1. Hereby they modulate the control of cholesterol, fatty acid, triglyceride, and lipoprotein metabolism. As a consequence, LXRs have emerged as promising drug targets for pharmacological LXR agonists to treat metabolic diseases like atherosclerosis and type 2 diabetes [28]. In rodents it has been shown that synthetic LXR agonists (T0901317, GW3965, and WAY252623) promote cellular cholesterol efflux, transport, and excretion, herewith halting the progression of atherosclerosis. However, pharmacological LXR activation also induces hepatic steatosis and promotes the secretion of enlarged atherogenic verylowdensitylipoprotein (VLDL) particles by the liver, complicating the clinical application of LXR agonists [29,30]. In the present study, we applied our computational approach to determine which metabolic processes change upon LXR activation, and identify the progression of molecular adaptations that collectively result in a shift of phenotype (wildtype versus LXR activated state). Parameters that are critical to the phenotype transition are considered candidates as biomarkers for disease diagnosis, treatment, or even prevention.
Methods
Several theoretical sections are presented describing the methodology of the computational approach, which involves consecutive steps of data gathering, model development, and various parameter estimation steps.
Model development
The computational approach is developed to analyze progressively adapting biological systems that are modeled using a system of (non)linear ordinary differential equations (ODEs):
where is a vector of flux descriptions of molecular species x, which are given by a set of functions f, that in turn depend on kinetic parameters θ and model inputs u. The initial concentrations of molecular species x are given by x_{0}. The vector of model outputs y is given by:
which is described by a set of functions g depending on molecular species x, kinetic parameters θ, and model inputs u.
Simulating the biological system of phenotype A
Once a network topology of molecular species and corresponding flux descriptions are defined, values for the kinetic parameters θ have to be specified in order to perform simulations and make predictions. One way of determining parameter values is to directly measure them. However, this could become impractical when it is not possible to perform the necessary experiments, or model parameters do not have a welldefined physiological meaning, e.g., when multiple processes are lumped into a single model parameter. Another way to obtain parameter values, which was employed here, is to estimate them by minimizing the difference between experimental data and corresponding model simulations [31]. The amount of experimental data is usually limited compared to the number of parameters, resulting in nonunique solutions for the model parameters. Hence, multiple parameter sets exist that adequately describe the experimental data. Conversely, model predictions of unmeasured molecular species might potentially vary greatly depending on the chosen parameter set. To assess the uncertainty associated with model predictions, differences between feasible parameter sets must be examined [3234]. A largescale parameter estimation protocol was employed to capture multiple parameter sets describing the biological system of phenotype A. First, parameter regions were identified that are most likely to describe the experimental data. To this end, a collection of one hundred million parameter sets was sampled from a loguniform distribution, capturing a parameter range of twelve orders of magnitude (10^{6 }to 10^{6}). For each parameter set a simulation to steadystate was carried out. Subsequently, the weighted sum of squared errors X_{d }(θ) between the experimental data of phenotype A and corresponding steadystate model outputs were determined:
where N is the number of measurements, d^{A }and σ^{A }respectively the means and standard deviations of the experimental data of phenotype A, and y the corresponding model outputs. Furthermore, a Monte Carlo approach was employed to account for experimental uncertainties. Each simulation a different realization for d^{A }was used. It was assumed that the experimental data is Gaussian distributed with means μ^{A }and standard deviations . Subsequently, the ten thousand best parameter sets (lowest X_{d }values) were selected and optimized to describe the experimental data, by applying a weighted nonlinear least squares algorithm that minimizes X_{d }(θ):
where represents the optimized parameter set. An optimized parameter set was acceptable if corresponding model outputs were in the confidence intervals of the experimental data. A significance level of 0.05, adjusted by Bonferroni correction to account for the number of comparisons being performed (number of model outputs), was used [35].
Identification of molecular adaptations from phenotype A to phenotype B
Parameter estimation to describe phenotype B
The mathematical model together with the collection of acceptable parameter sets, represents the biological system of phenotype A. Molecular processes that are responsible for the transition of the biological system from phenotype A to phenotype B, are determined by identifying kinetic parameters that necessarily have to change in order to describe the biological system of phenotype B. A first approach could be to repeat the largescale parameter estimation protocol, employed on phenotype A, for phenotype B. However, apart from being computationally expensive, comparing parameter sets from different phenotypes with each other is problematic, as they are obtained independently from each other. For instance, in the case when multiple separate minima exist, it would not be possible to know which realization of phenotype A is the reference for a specific realization of phenotype B. However, the fact that phenotype B originates from phenotype A could be used to address latter problem. The acceptable parameter sets from phenotype A could be used as initial values and reoptimized by once more applying a weighted nonlinear least squares algorithm, minimizing X_{d}(θ) with respect to the experimental data of phenotype B. Subsequently, necessary parameter adaptations can be identified which are responsible for the change of phenotype.
Iterative data integration and parameter estimation
Parameter adaptations describing a phenotype transition are often not unique. For instance, in order to increase a specific molecular concentration, corresponding production and degradation parameters can be changed in infinitely many different ways to accomplish this. Here, we assume that adaptations are minimal and proceed progressively in time. Therefore, the concept described in the previous section was extended to study progressively adapting biological systems, by defining artificial intermediate phenotypes. Hereto, the experimental data is interpolated from phenotype A to phenotype B in a number of steps. For instance, for a linear interpolation scheme this would imply d^{q }= (1  q)d^{A }+ qd^{B}, where d^{A }and d^{B }respectively represent the experimental data of phenotype A and B, and q a coordinate ranging from zero (completely phenotype A) to one (completely phenotype B). At each interpolation step the parameters are reoptimized in order to describe the newly interpolated data. The final values of the model states and parameters of the current optimization step are used as initial values for the next optimization step. This procedure is repeated until the final state representing phenotype B is reached and a parameter adaptation trajectory is obtained. The new objective function becomes as follows:
Similar as in equation (3), for each parameter trajectory different realizations for d^{A }and d^{B }were used to account for experimental uncertainties.
Regularization of parameter adaptation trajectories
It is assumed that adaptations are minimal and proceed progressively in time. Therefore, the parameter estimation protocol was extended to avoid needless change of parameters, hereby identifying minimal parameter adaptations that are necessary to describe a phenotype transition. To this end, X_{d }could be combined with a regularization term X_{r }given by the sum of squared parameter changes. When changing a parameter is costly, it will be avoided if not necessary. The new objective function is given by:
where M is the number of parameters, the initial parameter set representing phenotype A, and λ a constant determining the strength of the regularization term.
Consistency of parameter adaptation trajectories
The identification of parameter adaptation trajectories was performed for each acceptable parameter set, which gives information about the possible dispersion of parameter trajectories due to kinetic variations between the different acceptable parameter sets. However, given the uncertainties arising from experimental data and parameter estimates, the reliability of individual parameter trajectories is also a relevant topic to explore; is an observed trajectory consistent or is its path just a coincidental result? Given a certain parameter trajectory, it is important to analyze how reliable and consistent its path is to eventually draw conclusions about potential molecular adaptations that could have taken place. Therefore, the protocol described above was extended by not only determining parameter trajectories from phenotype A to phenotype B, but also backwards from phenotype B to phenotype A. A backward trajectory is obtained by interpolating the data from phenotype B to phenotype A, whilst reoptimizing the parameters. The final values of the model states and parameters obtained from the forward trajectory are used as initial values to calculate the backward trajectory. Furthermore, the reference parameter set (equation 6) is exchanged in this case by (the initial parameter set representing phenotype B) in order to regularize the backward trajectory. This process can be repeated an arbitrary number of steps, each time using the newly obtained values for the model states and parameters as initial values. The obtained parameter trajectories have been analyzed for consistency, which gives information regarding how well these adaptations are constrained by the data and can be predicted by the model. It must be noted that the calculation of backward trajectories is mainly a mathematical technique to assess the robustness of a specific solution. Hence, these trajectories do not necessarily have to exist physiologically.
Results
We presented a computational approach to analyze molecular adaptations in a biological system evolving through various phenotypes, which is generically applicable to different biological systems. In this section, the computational approach is demonstrated by applying it to a case of liver X receptor induced hepatic steatosis.
Experimental data
The acquisition of quantitative experimental data of different phenotypes is essential to gain insight in the progression of molecular adaptations in underlying biological systems. The available experimental data determines to a large extend the development of a mathematical model. The level of detail and precision at which certain biological processes can be integrated in a mathematical model, is determined by the selection of molecular species, as well as the type and quality of the measurements. With respect to the LXR case, several datasets of wildtype and T0901317 LXR activated C57BL/6J mice were obtained. Data was included containing measurements of hepatic triglyceride, free cholesterol, and cholesterylester levels, as well as plasma triglyceride, highdensitylipoprotein (HDL) cholesterol, total cholesterol, and free fatty acid levels in overnightfasted mice [29]. Furthermore, data of nascent produced VLDL particles such as diameter, triglyceride/cholesterol composition ratio, and VLDL triglyceride production rate was used [29]. Data was included containing rate measurements of hepatic cholesterol production, hepatic cholesterol uptake via HDL, and cholesterol uptake by peripheral tissues [36]. Information about the deposition and production of hepatic triglycerides in cytoplasmic and microsomal fractions was included [37]. An overview of the obtained experimental data is included in Additional file 1.
Additional file 1. Supplementary material. Description of model equations, additional analyses, implementation details, and experimental data.
Format: PDF Size: 880KB Download file
This file can be viewed with: Adobe Acrobat Reader
Computational model of hepatic lipid and plasma lipoprotein metabolism
A mathematical multicompartment model was constructed, based on the available experimental data, which integrates metabolic processes involved in hepatic lipid metabolism, as well as plasma lipoprotein metabolism (Figure 1). The mathematical model contains three compartments representing the liver, blood plasma, and periphery. The liver compartment includes reactions representing the production, utilization and storage of triglycerides and cholesterols. Furthermore, the model includes the mobilization of these metabolites to the endoplasmic reticulum, where they are incorporated into nascent produced VLDL particles. These VLDL particles are subsequently secreted in the plasma compartment where they serve as nutrients for peripheral tissues, e.g., muscle, heart, and adipose tissue. Remnant particles are taken up and cleared by the liver. The model furthermore includes the hepatic uptake of free fatty acids and the reverse transport of cholesterol via HDL. The model size and complexity of the reaction equations was kept to a minimum to preserve feasibility of model analyses and parameter estimation. The model developed contains eight molecular species x and twentytwo kinetic parameters θ. The flux descriptions f are all based on mass action kinetics. A description of the mathematical model, including equations, is presented in Additional file 1. Furthermore, an implementation of the model is available in SBML format (Additional file 2).
Figure 1. Hepatic lipid and plasma lipoprotein metabolism. The mathematical model has three compartments representing the liver, blood plasma, and peripheral tissues. The liver compartment includes reactions representing the production, utilization and storage of triglycerides and cholesterols, and the mobilization of these metabolites to the endoplasmic reticulum, where they are incorporated into nascent produced VLDL particles. The VLDL particles are secreted in the plasma compartment where they serve as nutrients for peripheral tissues. Remnant particles are taken up and cleared by the liver. The model furthermore includes the hepatic uptake of free fatty acids and reverse cholesterol transport via HDL. ABCA1, ATPbinding cassette transporter 1; ACAT, acylCoA: cholesterol acyltransferase; ApoB, apolipoprotein B; CE, cholesterylester; DGAT, diglyceride acyltransferase; ER, endoplasmic reticulum; FFA, free fatty acid; FC, free cholesterol; HDL, highdensitylipoprotein; HSL, hormonesensitive lipase; IDL, intermediate density lipoprotein; LDL, low density lipoprotein; LDLr, low density lipoprotein receptor; LPL, lipoprotein lipase; MTP, microsomal triglyceride transfer protein; SRB1, scavenger receptor class B1; TG, triglyceride; TGH, triglyceride hydrolase; VLDL, very low density lipoprotein; VLDLr, very low density lipoprotein receptor.
Additional file 2. SBML file. Implementation of the mathematical model in SBML format.
Format: XML Size: 39KB Download file
Simulating the wildtype mouse
A largescale parameter estimation protocol was employed to capture multiple parameter sets that describe the experimental data of phenotype A (wildtype C57BL/6J mice). Mass isotopomer distribution analyses indicate that the metabolic fluxes are expected to be in the μM/h range [38,39]. Therefore, parameter sets corresponding to unphysiologically high fluxes for any of the reactions (> 100 mM/h) were removed from further analyses. Finally, a collection of 2909 acceptable parameter sets was obtained that describe the experimental data. With respect to the parameter values, it appeared that several are very constrained by the data and have a well defined value, whereas others show a larger spread of possible outcomes. Figure 2 shows an example of four parameter combinations, in which the black dots represent the initial sampled parameter sets, the red dots represent the ten thousand best parameter sets, and the green dots represent the optimized acceptable parameter sets that describe the experimental data. The observed variation in several parameters is reflected in specific model predictions. Figure 2 shows two examples of model predictions obtained for all acceptable parameter sets for the depositioning of hepatic triglycerides and cholesterylesters in cytoplasmic and endoplasmic reticulum fractions. Note that only the total pools of triglycerides and cholesterylesters were measured [29]. Nonetheless, the predictions for the triglyceride fractions are consistent, due to the data of triglyceride deposition and production rates in cytoplasmic and endoplasmic reticulum fractions [37]. However, the predictions for the cholesterylester fractions show a larger spread of possible outcomes. The latter case illustrates the importance of exploring differences between feasible parameter sets to assess the uncertainty associated with model predictions.
Figure 2. Parameter scatter plots and predictions. Several parameters are very constrained by the data and have a well defined value (A and B), whereas others show a larger spread of possible values (D and E). The black dots represent the 10^{8 }initial sampled parameter values (individual dots not visible), whereas the red dots represent the 10^{4 }best parameter sets which were optimized. The resulting 2909 acceptable parameter sets that describe the experimental data are shown in green. Model predictions for the depositioning of hepatic triglycerides (C) and cholesterylesters (F) in cytoplasmic and endoplasmic reticulum fractions were obtained for all acceptable parameter sets. Note that only the total pools of triglyceride and cholesterylester were measured. The predictions for the triglyceride fractions are consistent, whereas the predictions for the cholesterylester fractions show a larger spread of possible outcomes.
Parameter adaptations from the wildtype to the LXR activated phenotype
Using the previously described techniques, an analysis was carried out to study the metabolic consequences of T0901317 induced LXR activation. It was assumed that metabolic adaptations upon LXR activation proceed linearly in time [40]. Therefore, a linear interpolation scheme was used for the stepwise optimization to describe the transition between phenotypes. A beneficial consequence of the approach is that the stepwise optimization guides the parameter estimation algorithm and hereby could overcome potential practical problems, such as convergence to local unacceptable minima. Figure 3 shows an example of an acceptable parameter set describing the wildtype phenotype, which was not successfully reoptimized by singlestep optimization to describe the LXR activated phenotype, whereas this problem was circumvented by multistep optimization.
Figure 3. Singlestep versus multistep optimization. A) An example of an acceptable parameter set describing the wildtype phenotype, which was not successfully optimized by singlestep optimization to describe the LXR activated phenotype. This problem was circumvented by multistep optimization. In the latter case, the parameter estimation is carried out in a stepwise fashion. Hereto, the experimental data is interpolated from the wildtype phenotype to the LXR activated phenotype. At each interpolation step the parameters are reoptimized in order to describe the newly interpolated data. This procedure is repeated until the final state representing phenotype B is reached. B, C) Two examples of corresponding parameter trajectories from the wildtype phenotype to the LXR activated phenotype, normalized by the initial parameters of the wildtype phenotype.
The parameter trajectories were regularized according to equation (6) to avoid needless change of parameters. A potential risk of regularization, as always with multiobjective optimization, is that for a low λ the regularization term has no effect, whereas for a large λ the parameter estimation algorithm might focus on minimizing the regularization term while describing the experimental data inaccurately. Therefore, the effect of λ on the sum of squared model errors X_{d }and the sum of squared parameter differences X_{r }was investigated for a collection of acceptable parameter sets. Figure 4a shows X_{d }for increasing λ, where green indicates an acceptable data fit and red an unacceptable data fit. Figure 4b shows X_{r }for increasing λ. Note that a small λ is already sufficient to minimize parameter changes, while the experimental data is still described very well. It is preferred to bias the data fitting as little as possible and therefore a λ of 0.1 was selected (both for the forward and backward trajectories). Figure 5 shows three examples of parameter trajectories from the wildtype phenotype to the LXR activated phenotype obtained without regularization (blue dashed) and with regularization (red). Both the regularized and unregularized parameter trajectories are acceptable in terms of model error X_{d}. Note that the triglyceride production and metabolism parameters counteract each others effect and not necessarily have to change to describe the change in phenotype (Figure 5a,b). In some cases a less prominent change is sufficient to describe the change in phenotype (Figure 5c).
Figure 4. Effect of regularization strength. A) Model error X_{d }for all acceptable parameter sets for increasing λ, where green indicates an acceptable data fit and red an unacceptable data fit. B) Regularization error X_{r }for increasing λ. Note that a small λ is already sufficient to minimize parameter changes, while describing the experimental data still very well. A λ of 0.1 was selected (dashed line).
Figure 5. Regularization of parameter trajectories. Three examples of parameter trajectories from the wildtype phenotype to the LXR activated phenotype obtained without regularization (blue dashed) and with regularization using λ = 0.1 (red).
The uncertainty associated with parameter trajectories was investigated, among other things, by repeatedly calculating forward and backward trajectories. Figure 6 shows three examples of back and forward parameter trajectories from the wildtype phenotype to the LXR activated phenotype, using a hundred repetitions. Some parameters change consistently (Figure 6a,b), whereas others show a large spread in possible outcomes (Figure 6c).
Figure 6. (In)consistency of parameter adaptation trajectories. Three examples of back and forward parameter trajectories from the wildtype phenotype to the LXR activated phenotype, using a hundred repetitions.
The parameter adaptation trajectories were determined for all acceptable parameter sets and used to determine how the fluxes of triglycerides and cholesterols change in time from the wildtype phenotype to the LXR activated phenotype. The data interpolation was carried out in a hundred steps and the back and forward flux trajectories were determined using a hundred repetitions. Figure 7 shows pairs of flux trajectories of several metabolic processes included in the model, where the large green and red dots respectively represent the wildtype phenotype and the LXR activated phenotype. The small dots represent the artificial intermediate phenotypes. The majority of these flux trajectories are reproduced very consistently for the different parameter sets. The main findings are that both the VLDLTG and VLDLCE production are increased (Figure 7a), whereas the production of apolipoprotein B is slightly decreased (Figure 7b). The hepatic and wholebody uptake of triglycerides and cholesterols are increased (Figure 7c, e, and 7f). The increased hepatic triglyceride fluxes are especially stored in cytosolic fractions, rather than in endoplasmic reticulum fractions (Figure 7d). Furthermore, the net synthesis of cholesterylester from endogenous free cholesterol is decreased in the cytosol, yet increased in the endoplasmic reticulum (Figure 7h).
Figure 7. Flux adaptations upon T0901317 induced LXR activation. Flux trajectories from the wildtype phenotype (green) to the LXR activated phenotype (red). The data interpolation was carried out in a hundred steps and the back and forward flux trajectories were determined using a hundred repetitions. Molecular fluxes (AC, EH) are given in mM/h, whereas the triglyceride concentrations presented in (D) are given in mM.
As described in previous sections, several parameters are not constrained by the data and show a large spread of possible outcomes. This makes the calculation of consistent quantitative trajectories impossible. Nonetheless, relative changes with respect to the initial values of the wildtype phenotype could still provide useful information, e.g., to determine ranges of feasible fold inductions of molecular concentrations and fluxes, and to discriminate between different possible scenarios. The latter could be used to generate new hypotheses and to guide the design of new experiments. An example is depicted in Figure 8, showing adaptations in metabolic processes/components involved in hepatic cholesterol metabolism, normalized by the values of corresponding wildtype phenotype. The green dots represent the wildtype phenotype, whereas the blue and black dots represent the LXR activated phenotype. A positive correlation between HDLCE synthesis and HDLCE uptake by the liver was observed (Figure 8a). Both fluxes are either increased or decreased depending on the chosen parameter set. To investigate how these different scenarios are reflected in other related metabolic processes, solutions corresponding to an increased HDLCE synthesis/uptake rate are colored blue, whereas solutions corresponding to a decreased HDLCE synthesis/uptake rate are colored black. Different clusters of possible scenarios exist depending on how the HDLCE synthesis/uptake rate adapts. The ellipses were calculated by means of principal component analysis (PCA) and contain 95% of the corresponding solutions.
Figure 8. Adaptations in cholesterol metabolism upon T0901317 induced LXR activation. Adaptations in metabolic processes/components involved in hepatic cholesterol metabolism, normalized by the values of corresponding wildtype phenotype. The green dots represent the wildtype phenotype, whereas the blue and black dots represent the LXR activated phenotype. Solutions corresponding to an increased HDLCE synthesis/uptake rate are colored blue, whereas solutions corresponding to a decreased HDLCE synthesis/uptake rate are colored black. The ellipses were calculated by means of principal component analysis (PCA) and contain 95% of the corresponding solutions. ABCA1, ATPbinding cassette transporter 1; ABCG5, ATPbinding cassette transporter G5; SRB1, scavenger receptor class B1.
Discussion
To improve our understanding of progressive diseases such as diabetes mellitus and metabolic syndrome, the study of phenotype transitions is important. A challenging task is to explain the progression of phenotype transitions in terms of molecular adaptations in underlying biological systems. Here, we propose a novel computational approach to analyze biological systems evolving through various phenotypes and to predict which molecular processes are responsible for the transition. We presented a case involving the activation of the liver X receptor, a promising drug target for atherosclerotic therapies.
Parameter adaptation trajectories during phenotype transitions: strengths and limitations
A largescale parameter estimation protocol was employed to capture multiple parameter sets describing the biological system of phenotype A. A collection of 2909 acceptable parameter sets were obtained that describe the experimental data. A substantial fraction of the optimized parameter sets were not acceptable. These parameter sets did either not describe the experimental data or displayed unphysiologically high fluxes for any of the reactions. It appeared that the latter criterion contributed significantly to the rejection of parameter sets. The efficiency of sampling acceptable parameter sets could potentially be improved by including the rejection criteria in the optimization objective function. Note, that the computational approach is not restricted to the parameter estimation protocol presented here. Various parameter optimization methods exist that minimize the difference between experimental data and corresponding model simulations, e.g., trustregion optimization methods, simulated annealing, and genetic algorithms [31]. All these methods have their own merits and shortcomings, and therefore the preference for a certain protocol varies on a casebycase basis.
Parameter trajectories describing the phenotype transition were determined by interpolating the data between phenotypes. The data interpolation was carried out in a hundred steps. We have performed several tests by using different numbers of steps. Performing more than a hundred steps did not change the results significantly. The computational approach allows free choice of interpolation schemes. Hence, when information is available about the progressive nature of certain biological processes, this information could be incorporated in the interpolation scheme. Furthermore, the computational approach could be used to explore different possible transition scenarios by employing a variety of different interpolation schemes. The latter could be useful when sufficient information about the transition characteristics is lacking, e.g., to test hypotheses about the feasibility of specific transitions with respect to the available experimental data. In this work we focused on diseases that arise progressively, e.g., hepatic steatosis, diabetes type 2, and metabolic syndrome. However, some diseases arise abruptly like in case of diabetes type 1. In latter cases it could be relevant to explore switchlike interpolation schemes and investigate whether the computational model can exhibit bistable behavior [4145]. Here, it has been assumed that metabolic adaptations upon LXR activation proceed linearly in time. Although there is limited information about the dynamic response upon T0901317 induced LXR activation, this assumption is supported by experimental observations from Okazaki et al. showing a fairly linear response in plasma triglyceride and cholesterol levels in wildtype and Ldlr ^{/ }mice treated with T0901317 [40]. Although initial and final points of the trajectories were consistent with experimental data, the actual trajectories between phenotypes partly depended on the selected interpolation scheme. If more timecourse data of LXR activated C57BL/6J mice would be available, a more realistic interpolation scheme could be defined. Although the dynamic behavior of parameter trajectories depends on the selected interpolation scheme, the relation between parameter trajectories (as visualized in Figure 7) does not necessarily have to change for different interpolation schemes. Namely, in the case all measured metabolite concentrations/fluxes adapt in a similar way, i.e. it can be assumed that the interpolation scheme is identical for each measurement, the relation between parameter trajectories remains identical. The results depicted in Figure 7 were reproduced using a quadraticlike and inversequadraticlike interpolation scheme for the measurements (Additional file 1). To identify minimal parameter adaptations that are necessary to describe a phenotype transition, the parameter estimation protocol was extended by including a regularization term given by the sum of squared parameter changes. This prevents unnecessarily changes of parameter values. The strength of the regularization term, determined by λ, was chosen carefully. It is preferred to bias the data fitting as little as possible and therefore a minimal value for λ, while still being effective, was selected. From a physical point of view, the regularization term could be interpreted as a measure for the energy cost, e.g., in terms of ATP production, to achieve a certain system adaptation. In future research, the approach could be refined by introducing multiple λ parameters to take account for different energy costs for the various processes included in a model.
Metabolic adaptations upon T0901317 induced LXR activation
A computational model of hepatic lipid and plasma lipoprotein metabolism was developed to study the metabolic consequences of LXR activation. We were able to quantitatively integrate data of wildtype and LXR activated C57BL/6J mice into a consistent model and identified trajectories of parameter adaptations, describing the change in phenotype. The presented model predictions are in good agreement with experimental observations by other groups and contribute to the current understanding of LXR activation. The VLDLTG production rate increases about 2.6 times upon LXR activation, as predicted by the model and experimentally measured [29]. A novel finding is that model predictions indicate that the VLDLCE production rate increases as well (2.5 fold induction) and hereby contributes to the increase of plasma cholesterol levels. Model predictions indicate that the production of apolipoprotein B decreases slightly, which was also observed by [29,30,46]. This is reflected in an increase of VLDL particle diameter (94 ± 12 nm to 129 ± 9 nm). A novel model prediction is that the liver plays a major role in the reuptake of lipoproteins (2.5 fold induction) and hereby prevents plasma hyperlipidemia. This flux prediction was not directly measured, but is in good agreement with gene expression data showing increased hepatic levels of the VLDL and LDL receptor [29]. Interestingly, model predictions indicate that the uptake of lipoproteins at peripheral tissues is negligible. Model analysis reveals that the uptake of triglycerides through lipolysis by lipases is increased as well (2.6 fold increase), which is in correspondence with gene expression data showing a significant induction of the lipoprotein lipase gene [29,30,47]. A significant increase in hepatic triglyceride level (6.92 ± 2.65 versus 57.74 ± 16.61 nmol/mg liver) was observed by [29], which is partly caused by an induction of lipogenic genes [29,30,4650]. A novel model prediction is that the increased triglyceride fluxes are especially stored in cytosolic fractions, rather than in endoplasmic reticulum fractions which are predominantly used for incorporation into nascent produced VLDL particles. The increased level of ER triglycerides is partly caused by stimulation of the mobilization of triglycerides from the cytosol to the ER. This is confirmed by several studies indicating that a large part of secreted VLDL triglycerides are derived via lipolysis of cytosolic triglyceride storage pools [5155]. A relevant followup study would be to determine whether these differences are associated with alterations in diglyceride acyltransferase activities (DGAT1 and DGAT2), which play a crucial role in the biosynthesis and deposition of triglycerides [5659]. Another interesting example which could guide the design of new experiments is depicted in Figure 8, showing adaptations in metabolic processes/components involved in hepatic cholesterol metabolism. Different clusters of possible scenarios exist depending on how the HDLCE synthesis/uptake rate adapts. Hence, many solutions could potentially be excluded by measuring one of these fluxes/components. With respect to this, the 'blue' scenario is probably more plausible for several reasons. First, these solutions are associated with an increased level of the ATPbinding cassette transporter G5 (ABCG5), resulting in an increased biliary excretion of free cholesterol. Secondly, these solutions correspond to an increased level of the ABCA1 transporter, which is responsible for the efflux of cholesterol from peripheral tissues to HDL [30,47,49,50].
Mathematical modeling of progressively adapting biological systems
Mathematical modeling is well suited for integrating different sources of experimental data for a certain phenotype and allows investigating of the complex interactions of underlying biological systems. A mathematical model can be used to obtain thorough understanding of a biological system, e.g. by investigating its complex behavior in response to various stimuli. However, simulating and predicting longterm progressive adaptations is challenging. The multiscale nature of progressively adapting biological systems complicates the development of predictive computational models. As such, one of the central and formidable challenges in systems biology is to develop multiscale computational models and methods that can be used to study molecular mechanisms underlying progressive diseases [6065]. Furthermore, model parameters that determine the kinetics of molecular processes are often assumed to be constant in time and between phenotypes. This is most probably a valid assumption to study shortterm processes, e.g., initial response kinetics to perturbations of a biological system. In case of progressively adapting biological processes, it is questionable whether this assumption still holds. For instance, effects of aging, changes in body composition, or other developmental changes, influence the phenotypical characteristics and the transition between phenotypes.
The computational approach presented here was employed to study the metabolic consequences of LXR activation, which displays several of the aforementioned issues. For example, the underlying biological system contains processes at timescales ranging from seconds to hours, whereas the phenotypical characteristics develop at a timescale ranging from days to weeks in mice. Our approach has as advantage that it can readily deal with large differences in timescales. For instance, longterm changes in shortterm processes could be studied by explicitly modeling the shortterm processes, whereas the longterm modulation could be captured by identifying necessary parameter changes. This implies that molecular adaptations could be described without the necessity to develop detailed kinetic models of the modulating mechanisms. This is a major advantage, e.g., for the LXR case, as LXRs modulate a wide range of heavily interlinked complex metabolic processes and signal transduction pathways of which the kinetics and molecular mechanisms are not well understood.
Conclusions
The study of phenotype transitions is important to understand disease progression. We developed a novel computational approach to analyze molecular adaptations in a biological system evolving through various phenotypes, which is generically applicable to different biological systems. For the case of liver X receptor induced hepatic steatosis the novel approach yields information about the redistribution of fluxes and pools of triglycerides and cholesterols that was not directly apparent from the experimental data. The collection of parameter and corresponding flux trajectories give a broad overview of keyprocesses that are involved in the phenotype transition and how they potentially change in time. Model analysis provides guidance which specific molecular processes to study in more detail to obtain further understanding of the underlying biological system. The main findings are that both the VLDLTG and VLDLCE production rates are increased, as well as the uptake of triglycerides through lipolysis. The liver plays a major role in the reuptake of lipoproteins and hereby prevents plasma hyperlipidemia. The increased triglyceride fluxes are especially stored in cytosolic fractions, rather than in endoplasmic reticulum fractions.
Authors' contributions
CT developed the mathematical model, performed the computational analysis, and wrote the paper. JV contributed to the computational analysis, improved the computational efficiency of the software, and revised the paper. PH supervised the study and revised the paper. NR supervised the study, contributed to the computational analysis, and revised the paper. All authors read and approved the final manuscript.
Acknowledgements
Research was funded by the Netherlands Consortium for Systems Biology (NCSB). We gratefully thank Aldo Grefhorst, Maaike Oosterveer, Barbara Bakker, and Albert Groen for useful discussions.
References

Donath M, Shoelson S: Type 2 diabetes as an inflammatory disease.
Nat Rev Immunol 2011, 11:98107. PubMed Abstract  Publisher Full Text

Minamino T, Komuro I: Vascular aging: insights from studies on cellular senescence, stem cell aging, and progeroid syndromes.
Nat Clin Pract Card 2008, 5(10):637648. Publisher Full Text

Drüeke T, Massy Z: Atherosclerosis in CKD: differences from the general population.
Nat Rev Nephrol 2010, 6(12):723735. PubMed Abstract  Publisher Full Text

Yeh M, Brunt E: Pathology of nonalcoholic fatty liver disease.
Am J Clin Pathol 2007, 128(5):837847. PubMed Abstract  Publisher Full Text

Rutledge J, Ng K, Aung H, Wilson D: Role of triglyceriderich lipoproteins in diabetic nephropathy.
Nat Rev Nephrol 2010, 6(6):361370. PubMed Abstract  Publisher Full Text

Speakman J, Hambly C, Mitchell S, Krol E: The contribution of animal models to the study of obesity.
Lab Anim 2008, 42(4):413432. PubMed Abstract  Publisher Full Text

Manis J: Knock out, knock in, knock downgenetically manipulated mice and the Nobel Prize.
N Engl J Med 2007, 357(24):24262429. PubMed Abstract  Publisher Full Text

Voshol P, Rensen P, van Dijk K, Romijn J, Havekes L: Effect of plasma triglyceride metabolism on lipid storage in adipose tissue: studies using genetically engineered mouse models.
Biochim Biophys Acta 2009, 1791(6):479485. PubMed Abstract  Publisher Full Text

Ilan Y, Maron R, Tukpah A, Maioli T, Murugaiyan G, Yang K, Wu H, Weiner H: Induction of regulatory T cells decreases adipose inflammation and alleviates insulin resistance in ob/ob mice.
P Natl Acad Sci USA 2010, 107(21):97659770. Publisher Full Text

Zhang C, Baffy G, Perret P, Krauss S, Peroni O, Grujic D, Hagen T, VidalPuig A, Boss O, Kim Y, et al.: Uncoupling Protein2 Negatively Regulates Insulin Secretion and Is a Major Link between Obesity, β Cell Dysfunction, and Type 2 Diabetes.
Cell 2001, 105(6):745755. PubMed Abstract  Publisher Full Text

De Luca C, Kowalski T, Zhang Y, Elmquist J, Lee C, Kilimann M, Ludwig T, Liu S, Chua Jr S: Complete rescue of obesity, diabetes, and infertility in db/db mice by neuronspecific LEPRB transgenes.
J Clin Invest 2005, 115(12):34843493. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Banks A, Kon N, Knight C, Matsumoto M, GutiérrezJuárez R, Rossetti L, Gu W, Accili D: SirT1 gain of function increases energy efficiency and prevents diabetes in mice.
Cell Metab 2008, 8(4):333341. PubMed Abstract  Publisher Full Text

Fang X, Wallqvist A, Reifman J: A systems biology framework for modeling metabolic enzyme inhibition of Mycobacterium tuberculosis.
BMC Syst Biol 2009, 3:92. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Wu X, Jiang R, Zhang M, Li S: Networkbased global inference of human disease genes.
Mol Syst Biol 2008, 4:189. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Young D, Stark J, Kirschner D: Systems biology of persistent infection: tuberculosis as a case study.
Nat Rev Microbiol 2008, 6(7):520528. PubMed Abstract  Publisher Full Text

Lusis A, Attie A, Reue K: Metabolic syndrome: from epidemiology to systems biology.
Nat Rev Genet 2008, 9(11):819830. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Stefanini M, Wu F, Mac Gabhann F, Popel A: A compartment model of VEGF distribution in blood, healthy and diseased tissues.
BMC Syst Biol 2008, 2:77. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Hu P, Bader G, Wigle D, Emili A: Computational prediction of cancergene function.
Nat Rev Cancer 2006, 7:2334. PubMed Abstract  Publisher Full Text

Kitano H: Computational systems biology.
Nature 2002, 420(6912):206210. PubMed Abstract  Publisher Full Text

Schoeberl B, EichlerJonsson C, Gilles E, Müller G: Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors.
Nat Biotechnol 2002, 20(4):370375. PubMed Abstract  Publisher Full Text

Koschorreck M, Gilles E: Mathematical modeling and analysis of insulin clearance in vivo.
BMC Syst Biol 2008, 2:43. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Schilling M, Maiwald T, Hengl S, Winter D, Kreutz C, Kolch W, Lehmann W, Timmer J, Klingmüller U: Theoretical and experimental analysis links isoformspecific ERK signalling to cell fate decisions.
Mol Syst Biol 2009, 5:334. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Borisov N, Aksamitiene E, Kiyatkin A, Legewie S, Berkhout J, Maiwald T, Kaimachnikov N, Timmer J, Hoek J, Kholodenko B: Systemslevel interactions between insulinEGF networks amplify mitogenic signaling.
Mol Syst Biol 2009, 5:256. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Cedersund G, Roll J, Ulfhielm E, Danielsson A, Tidefelt H, Strålfors P: Modelbased hypothesis testing of key mechanisms in initial phase of insulin signaling.

Repa J, Mangelsdorf D: The liver X receptor gene team: potential new players in atherosclerosis.
Nat med 2002, 8(11):12431248. PubMed Abstract  Publisher Full Text

Joseph S, Castrillo A, Laffitte B, Mangelsdorf D, Tontonoz P: Reciprocal regulation of inflammation and lipid metabolism by liver X receptors.
Nat med 2003, 9(2):213219. PubMed Abstract  Publisher Full Text

Janowski B, Willy P, Devi T, Falck J, Mangelsdorf D: An oxysterol signalling pathway mediated by the nuclear receptor LXRα.
Nature 1996, 383(6602):728731. PubMed Abstract  Publisher Full Text

Oosterveer M, Grefhorst A, Groen A, Kuipers F: The liver X receptor: Control of cellular lipid homeostasis and beyond: Implications for drug design.
Prog Lipid Res 2010, 49(4):343352. PubMed Abstract  Publisher Full Text

Grefhorst A, Elzinga B, Voshol P, Plösch T, Kok T, Bloks V, van der Sluijs F, Havekes L, Romijn J, Verkade H, et al.: Stimulation of lipogenesis by pharmacological activation of the liver X receptor leads to production of large, triglyceriderich very low density lipoprotein particles.
J Biol Chem 2002, 277(37):3418234190. PubMed Abstract  Publisher Full Text

Grefhorst A, Parks E: Reduced insulinmediated inhibition of VLDL secretion upon pharmacological activation of the liver X receptor in mice.
J Lipid Res 2009, 50(7):13741383. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Ashyraliyev M, FomekongNanfack Y, Kaandorp J, Blom J: Systems biology: parameter estimation for biochemical models.
FEBS J 2009, 276(4):886902. PubMed Abstract  Publisher Full Text

Cedersund G, Roll J: Systems biology: model based evaluation and comparison of potential explanations for given biological data.
FEBS J 2009, 276(4):903922. PubMed Abstract  Publisher Full Text

Raue A, Kreutz C, Maiwald T, Bachmann J, Schilling M, Klingmüller U, Timmer J: Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood.
Bioinformatics 2009, 25(15):1923. PubMed Abstract  Publisher Full Text

Gutenkunst R, Waterfall J, Casey F, Brown K, Myers C, Sethna J: Universally sloppy parameter sensitivities in systems biology models.
PLoS Comput Biol 2007, 3(10):18711878. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Osono Y, Woollett L, Herz J, Dietschy J: Role of the low density lipoprotein receptor in the flux of cholesterol through the plasma and across the tissues of the mouse.
J Clin Invest 1995, 95(3):11241132. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Millar J, Stone S, Tietge U, Tow B, Billheimer J, Wong J, Hamilton R, Farese R Jr, Rader D: Shortterm overexpression of DGAT1 or DGAT2 increases hepatic triglyceride but not VLDL triglyceride or apoB production.
J Lipid Res 2006, 47(10):22972305. PubMed Abstract  Publisher Full Text

Oosterveer M, Van Dijk T, Tietge U, Boer T, Havinga R, Stellaard F, Groen A, Kuipers F, Reijngoud D: High fat feeding induces hepatic fatty acid elongation in mice.
PLoS One 2009, 4(6):e6066. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Wiegman C, Bandsma R, Ouwens M, van der Sluijs F, Havinga R, Boer T, Reijngoud D, Romijn J, Kuipers F: Hepatic VLDL production in ob/ob mice is not stimulated by massive de novo lipogenesis but is less sensitive to the suppressive effects of insulin.
Diabetes 2003, 52(5):10811089. PubMed Abstract  Publisher Full Text

Okazaki H, Goldstein J, Brown M, Liang G: LXRSREBP1cphospholipid transfer protein axis controls very low density lipoprotein (VLDL) particle size.
J Biol Chem 2010, 285(9):6801. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Csete M, Doyle J: Bow ties, metabolism and disease.
TRENDS in Biotechnology 2004, 22(9):446450. PubMed Abstract  Publisher Full Text

Kitano H, Oda K, Kimura T, Matsuoka Y, Csete M, Doyle J, Muramatsu M: Metabolic syndrome and robustness tradeoffs.
Diabetes 2004, 53(suppl 3):S6. PubMed Abstract  Publisher Full Text

Conradi C, Flockerzi D, Raisch J: Multistationarity in the activation of a MAPK: Parametrizing the relevant region in parameter space.
Mathematical biosciences 2008, 211:105131. PubMed Abstract  Publisher Full Text

SaezRodriguez J, HammerleFickinger A, Dalal O, Klamt S, Gilles E, Conradi C: Multistability of signal transduction motifs.
Systems Biology, IET 2008, 2(2):8093. Publisher Full Text

Wilhelm T: The smallest chemical reaction system with bistability.
BMC Sys Biol 2009, 3:90. BioMed Central Full Text

Schultz J, Tu H, Luk A, Repa J, Medina J, Li L, Schwendner S, Wang S, Thoolen M, Mangelsdorf D, et al.: Role of LXRs in control of lipogenesis.
Genes Dev 2000, 14(22):28312838. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Chisholm J, Hong J, Mills S, Lawn R: The LXR ligand T0901317 induces severe lipogenesis in the db/db diabetic mouse.
J Lipid Res 2003, 44(11):20392048. PubMed Abstract  Publisher Full Text

Ou J, Tu H, Shan B, Luk A, DeBoseBoyd R, Bashmakov Y, Goldstein J, Brown M: Unsaturated fatty acids inhibit transcription of the sterol regulatory elementbinding protein1c (SREBP1c) gene by antagonizing liganddependent activation of the LXR.
P Natl Acad Sci USA 2001, 98(11):60276032. Publisher Full Text

Miao B, Zondlo S, Gibbs S, Cromley D, Hosagrahara V, Kirchgessner T, Billheimer J, Mukherjee R: Raising HDL cholesterol without inducing hepatic steatosis and hypertriglyceridemia by a selective LXR modulator.
J Lipid Res 2004, 45(8):14101417. PubMed Abstract  Publisher Full Text

Quinet E, Savio D, Halpern A, Chen L, Miller C, Nambi P: Geneselective modulation by a synthetic oxysterol ligand of the liver X receptor.
J Lipid Res 2004, 45(10):19291942. PubMed Abstract  Publisher Full Text

Gibbons G, Islam K, Pease R: Mobilisation of triacylglycerol stores.
Biochim Biophys Acta 2000, 1483:3757. PubMed Abstract  Publisher Full Text

Gibbons G, Bartlett S, Sparks C, Sparks J: Extracellular fatty acids are not utilized directly for the synthesis of verylowdensity lipoprotein in primary cultures of rat hepatocytes.
Biochem J 1992, 287(Pt 3):749753. PubMed Abstract  PubMed Central Full Text

Wiggins D, Gibbons G: The lipolysis/esterification cycle of hepatic triacylglycerol. Its role in the secretion of verylowdensity lipoprotein and its response to hormones and sulphonylureas.
Biochem J 1992, 284(Pt 2):457462. PubMed Abstract  PubMed Central Full Text

Yang L, Kuksis A, Myher J, Steiner G: Origin of triacylglycerol moiety of plasma very low density lipoproteins in the rat: structural studies.
J Lipid Res 1995, 36:125136. PubMed Abstract  Publisher Full Text

Lankester D, Brown A, Zammit V: Use of cytosolic triacylglycerol hydrolysis products and of exogenous fatty acid for the synthesis of triacylglycerol secreted by cultured rat hepatocytes.
J Lipid Res 1998, 39(9):18891895. PubMed Abstract  Publisher Full Text

Yen C, Stone S, Koliwad S, Harris C, Farese R Jr: Thematic review series: glycerolipids. DGAT enzymes and triacylglycerol biosynthesis.
J Lipid Res 2008, 49(11):22832301. PubMed Abstract  Publisher Full Text

Monetti M, Levin M, Watt M, Sajan M, Marmor S, Hubbard B, Stevens R, Bain J, Newgard C, Farese R Sr, et al.: Dissociation of hepatic steatosis and insulin resistance in mice overexpressing DGAT in the liver.
Cell Metab 2007, 6:6978. PubMed Abstract  Publisher Full Text

Smith S, Cases S, Jensen D, Chen H, Sande E, Tow B, Sanan D, Raber J, Eckel R, Farese R: Obesity resistance and multiple mechanisms of triglyceride synthesis in mice lacking Dgat.
Nat Genet 2000, 25:8790. PubMed Abstract  Publisher Full Text

Cases S, Smith S, Zheng Y, Myers H, Lear S, Sande E, Novak S, Collins C, Welch C, Lusis A, et al.: Identification of a gene encoding an acyl CoA: diacylglycerol acyltransferase, a key enzyme in triacylglycerol synthesis.
P Natl Acad Sci USA 1998, 95(22):1301813023. Publisher Full Text

Noble D: Modeling the heartfrom genes to cells to the whole organ.
Science 2002, 295(5560):16781682. PubMed Abstract  Publisher Full Text

Kohl P, Noble D: Systems biology and the virtual physiological human.
Mol Syst Biol 2009, 5:292. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hunter P, Borg T: Integration from proteins to organs: the Physiome Project.
Nat Rev Mol Cell Bio 2003, 4(3):237243. Publisher Full Text

Kuepfer L: Towards wholebody systems physiology.
Mol Syst Biol 2010, 6:409. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

De Graaf A, Freidig A, De Roos B, Jamshidi N, Heinemann M, Rullmann J, Hall K, Adiels M, Van Ommen B: Nutritional systems biology modeling: from molecular mechanisms to physiology.
PLoS Comput Biol 2009, 5(11):e1000554. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Vincent N, de Chassey Benoit C, Vincent L: When the human viral infectome and diseasome networks collide: towards a systems biology platform for the aetiology of human diseases.
BMC Syst Biol 2011, 5:13. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text