Parameter adaptations during phenotype transitions in progressive diseases
© Tiemann et al; licensee BioMed Central Ltd. 2011
Received: 15 June 2011
Accepted: 26 October 2011
Published: 26 October 2011
Skip to main content
© Tiemann et al; licensee BioMed Central Ltd. 2011
Received: 15 June 2011
Accepted: 26 October 2011
Published: 26 October 2011
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.
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.
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.
Cardiovascular and metabolic diseases such as diabetes mellitus and metabolic syndrome are progressive in time [1–5]. 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 [6–8]. For instance, the genetic leptin-deficient (ob/ob) or leptin-resistant (db/db) mouse are frequently used to study metabolic pathologies, e.g., obesity, insulin resistance, and diabetes [9–12]. 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 [13–19]. 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 [20–24]. 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 long-term adaptations are captured by identifying necessary parameter changes. This enables us to study long-term aspects of short-term 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 ligand-activated by oxysterols, heterodimerize with the retinoid X receptor (RXR) and bind to LXR responsive elements on the DNA , where they induce the transcription of lipogenic genes such as SREBP-1C, 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 . 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 very-low-density-lipoprotein (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 (wild-type versus LXR activated state). Parameters that are critical to the phenotype transition are considered candidates as biomarkers for disease diagnosis, treatment, or even prevention.
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.
which is described by a set of functions g depending on molecular species x, kinetic parameters θ, and model inputs u.
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 .
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 large-scale 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 non-linear 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.
Similar as in equation (3), for each parameter trajectory different realizations for d A and d B were used to account for experimental uncertainties.
where M is the number of parameters, the initial parameter set representing phenotype A, and λ a constant determining the strength of the regularization term.
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.
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.
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 wild-type 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, high-density-lipoprotein (HDL) cholesterol, total cholesterol, and free fatty acid levels in overnight-fasted mice . Furthermore, data of nascent produced VLDL particles such as diameter, triglyceride/cholesterol composition ratio, and VLDL triglyceride production rate was used . Data was included containing rate measurements of hepatic cholesterol production, hepatic cholesterol uptake via HDL, and cholesterol uptake by peripheral tissues . Information about the deposition and production of hepatic triglycerides in cytoplasmic and microsomal fractions was included . An overview of the obtained experimental data is included in Additional file 1.
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.
A large-scale 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., trust-region optimization methods, simulated annealing, and genetic algorithms . All these methods have their own merits and shortcomings, and therefore the preference for a certain protocol varies on a case-by-case 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 switch-like interpolation schemes and investigate whether the computational model can exhibit bistable behavior [41–45]. 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 wild-type and Ldlr -/- mice treated with T0901317 . 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 time-course 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 quadratic-like and inverse-quadratic-like 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.
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 wild-type 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 VLDL-TG production rate increases about 2.6 times upon LXR activation, as predicted by the model and experimentally measured . A novel finding is that model predictions indicate that the VLDL-CE 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 re-uptake 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 . 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 , which is partly caused by an induction of lipogenic genes [29, 30, 46–50]. 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 [51–55]. A relevant follow-up 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 [56–59]. 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 HDL-CE 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 ATP-binding 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 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 long-term 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 [60–65]. 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 short-term 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, long-term changes in short-term processes could be studied by explicitly modeling the short-term processes, whereas the long-term 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.
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 key-processes 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 VLDL-TG and VLDL-CE production rates are increased, as well as the uptake of triglycerides through lipolysis. The liver plays a major role in the re-uptake of lipoproteins and hereby prevents plasma hyperlipidemia. The increased triglyceride fluxes are especially stored in cytosolic fractions, rather than in endoplasmic reticulum fractions.
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.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.