- Open Access
Differentiation between genomic and non-genomic feedback controls yields an HPA axis model featuring Hypercortisolism as an irreversible bistable switch
Theoretical Biology and Medical Modellingvolume 10, Article number: 65 (2013)
The hypothalamic-pituitary-adrenal axis (HPA axis) is a major part of the neuroendocrine system responsible for the regulation of the response to physical or mental stress and for the control of the synthesis of the stress hormone cortisol. Dysfunctions of the HPA axis characterized by either low (hypocortisolism) or increased (hypercortisolism) cortisol levels are implicated in various pathological conditions. Their understanding and therapeutic correction may be supported by mathematical modeling and simulation of the HPA axis.
Mass action and Michaelis Menten enzyme kinetics were used to provide a mechanistic description of the feedback mechanisms within the pituitary gland cells by which cortisol inhibits its own production. A separation of the nucleus from the cytoplasm by compartments enabled a differentiation between slow genomic and fast non-genomic processes. The model in parts was trained against time resolved ACTH stress response data from an in vitro cell culture of murine AtT-20 pituitary tumor cells and analyzed by bifurcation discovery tools.
A recently found pituitary gland cell membrane receptor that mediates rapid non-genomic actions of glucocorticoids has been incorporated into our model of the HPA axis. As a consequence of the distinction between genomic and non-genomic feedback processes our model possesses an extended dynamic repertoire in comparison to existing HPA models. In particular, our model exhibits limit cycle oscillations and bistable behavior associated to hypocortisolism but also features a (second) bistable switch which captures irreversible transitions in hypercortisolism to elevated cortisol levels.
Model predictive control and inverse bifurcation analysis have been previously applied in the simulation-based design of therapeutic strategies for the correction of hypocortisolism. Given the HPA model extension presented in this paper, these techniques may also be used in the study of hypercortisolism. As an example, we show how sparsity enforcing penalization may suggest network interventions that allow the return from elevated cortisol levels back to nominal ones.
The HPA axis is a major part of the neuroendocrine system in mammals and particularly in humans. A main task of this hormonal network is the regulation of the response to physical or mental stress that threatens to disrupt the homeostatic balance of the organism. If a stressor is sensed by the nervous system the hypothalamus is stimulated to produce and secret the corticotropin-releasing hormone (CRH). The secretion of CRH causes the anterior pituitary to synthesize adrenocorticotropin (ACTH). ACTH then stimulates the adrenal glands to release cortisol, which down regulates the blood concentration of CRH and ACTH via different negative feedback mechanisms[1–4].
The HPA axis is subject of intensive research in endocrinology as HPA malfunctions are implicated in various pathological conditions. These are often characterized by either elevated or insufficient blood cortisol levels compared to the average healthy human. For instance, hypocortisolism (insufficient cortisol level) is reported in patients suffering from the chronic fatigue syndrome and post traumatic stress disorders (cf.[5–9]), whereas hypercortisolism (elevated cortisol level) is observed in depression, dementia or postoperative delirium (cf.[10–15]). Especially in the context of personalized medicine the use of modeling and simulation of biological systems for the rational design of treatments and drug intervention strategies is more and more recognized[16–20]. For such endeavors the integration of biological information of different type into computational, hence predictive, models is a prerequisite. The emphasis of earlier HPA modeling approaches with ordinary and delay differential equations has been put on self regulatory ultradian and circadian oscillatory behavior in[21–27], oscillations in response to an independent external pacemaker from the suprachiazmatic nucleus have been described in[28, 29]. In comparison, the article stands out as it incorporates intracellular glucocorticoid receptor kinetics which mediate bistable behavior of the HPA axis. Despite its parsimonious character the four state ODE model of offers an explanation for hypocortisolism as an irreversible biological switch and served in as a basis for the design of a therapeutic corrections of the HPA dysfunction. In it is shown that the model of also exhibits stable limit cycle oscillations, in the four state rate equations of were modified in order to fit oscillatory data of patients suffering from post traumatic stress disorder.
A parsimonious HPA model featuring hypocortisolism
The HPA axis model of captures the basic feedback mechanisms and includes an intracellular glucocorticoid receptor GR in the anterior pituitary gland as one of the four state variables, see Figure1. The dynamics of the model are described by the ODE system
where the ODE variables [CRH], [ACTH], [GR] and [COR] denote the concentrations of CRH, ACTH, GR and cortisol. If GR is bound to cortisol, the resulting complex is denoted by [COR-GR]. As this binding is very fast compared to the other dynamics it is assumed to be in equilibrium, i.e. [COR-GR] = [COR] · [GR]. The GR dimerization may generate an irreversible bistable switch in the model by which the dysfunctional abidance of the system at an alternative stable steady state characterized by low cortisol level can be explained. The therapeutic strategy of for the correction of hypocortisolism based on (1)-(4) is to force the HPA axis by a further suppression of cortisol to a point where the only stable condition in proximity is the original healthy state. An alternative for supporting the search for pharmaceutical therapies is the manipulation of the feedback network by means of inverse bifurcation techniques in order to transform the irreversible switch into a reversible one, see[19, 33]. Though the model of sheds new light on possible causes for hypocortisolism, its dynamic repertoire could not be shown to feature the opposite scenario of hypercortisolism. Furthermore, while the model takes the intracellular glucocorticoid receptor GR into account, it does not consider feedback mechanisms known to act on faster time scales.
Motivation of model extensions
The above mentioned GR mediated negative feedback on the production of ACTH acts at the genomic level. Specifically, cortisol COR binds to the glucocorticoid receptor GR and forms a homodimer (cf.[34, 35]). This homodimer is known to work as a transcription factor, enhances the expression of the GR gene and down regulates the expression of the ACTH related gene. A consequence of the transcriptional regulation is that the GR feedback control works rather slowly compared to other cellular processes such as a vesicle transport, cell signaling or extracellular events such as changes in the CRH blood concentration. For mammalian cells one can expect at least a delay of the down regulation in the range of 15 minutes up to 2 hours (c.f.). However, it has been frequently reported that cortisol induced inhibitory effects on anterior pituitary gland cells can already occur after a few seconds (c.f.[36, 37]), a phenomenon that cannot be explained by means of the genomic feedback mechanism.
The common hypothesis for this observation is that besides the GR receptor a second receptor in the membrane may directly inhibit the release of ACTH when bound to cortisol (cf.[38–42]). In the authors were able to provide evidence of a glucocorticoid receptor (GPCR) in the anterior pituitary cell membrane, possibly enabling the fast response related to cortisol induced inhibition. In the following we incorporate this fast and non-genomic feedback mechanism and for this purpose integrate a detailed pituitary gland cell model into the HPA network of. As a consequence, the extended dynamic repertoire of the resulting HPA axis feedback network also features bistable behavior compatible with the dysfunctional state of hypercortisolism and the ability to simulate time resolved ACTH stress response data from in vitro AtT-20 pituitary tumor cells.
Using mass action and Michaelis Menten enzyme kinetics we developed a mechanistic ODE system model of the glucocorticoid feedback mechanisms within the anterior pituitary gland cell. In particular, we incorporated the glucocorticoid membrane receptor GPCR found in and distinguished between slow genomic and fast non-genomic processes by using different compartments for the nucleus and the cytoplasm. Figure2 displays the reaction network graph of the resulting HPA axis model while the details are discussed in the Methods section.
Model comparison with experimental data
The experimentally observed fast inhibitory effect of cortisol on the anterior pituitary gland cell motivated the model incorporation of the GPCR receptor in order to allow for non-genomic, hence fast, signal processing. In order to test the capability of our model to capture these fast effects against actual data, we conducted experiments in cultures of the murine AtT-20 pituitary cell line (cf.). After administering doses of CRH and cortisol samples were taken from the medium and analyzed with respect to the total ACTH concentration in the cell medium, see Figure4A for the data obtained for two different experimental setups. In the first setup (indicated by the triangle markers in Figure4A) we considered the situation of a single stress stimulus by 10nM CRH at time point zero in the absence of cortisol. We observed an immediate increase of the ACTH secretion, followed by a slight stagnation and then a constantly increasing total ACTH concentration. In the second experiment (indicated by the dot markers in Figure4A) 10nM CRH and in addition 100nM cortisol were administered. As depicted, this immediately activated the fast non-genomic pathway inhibiting the release of ACTH. This led to an increased amount of ACTH stored in the vesicles within the cell and thus to a delayed but strong release of ACTH after approximately 15 minutes. However, after 30 to 60 minutes the genomic feedback mechanisms become effective resulting in an overall reduced ACTH concentration in the medium.
For training the model against the experimental data by adaptation of the model parameters we followed a variational regularization approach (cf.[33, 45]), see Methods for details. Sparsity enforcing penalization was used to eliminate model components less relevant for reproducing the observed kinetics. Figure4B shows that the experimental data can be fitted by the HPA model of Figure2 even if the autocatalytic loop for the GR gene regulation would be neglected and if the production of the TF (transcription factors for CRH) receptors would be considered to be independent of CRH.
In order to assess the relevance of the featured GPCR receptor for the model we knocked out the corresponding pathway by elimination of the respective ODE equations from our model and addressed the data fitting task with the reduced model. The result in the data space is shown in Figure4C and indicates that the consideration of the GPCR feedback pathway in the model is in fact vital for reproducing the experimentally observed dynamics. For comparison we finally run the data fitting procedure with the parsimonious model (1) - (4) from. The model lacks a GPCR mediated feedback pathway for fast response and fails to reproduce the corresponding experimental observation, see Figure4D.
Exploration of the dynamic repertoire of the model
One important aspect of is the inclusion of the glucocorticoid receptor GR which generates bistability for a model based explanation of hypocortisolism as an irreversible biological switch. As the GR mediated feedback pathway is also a constituent of the extended HPA network model of Figure2, the capability of our model to reflect the bistable behavior of (1)-(4) is an imperative requirement. Tools[46, 47] for studying the dynamic repertoire of ODE systems with respect to bistability or oscillations solve optimization problems that involve the eigenvalues of the Jacobian matrix f x (x,q) (with f denoting the right hand side of an ODE model) and search for regions in the parameter space that give rise to limit-point or Hopf bifurcations. Figure5B displays a resulting bifurcation diagram for the extended HPA model from Figure2 with the stress level of Equation (8) as the bifurcation parameter. It gives an explanation of the dysfunctional alteration from a steady state of normal cortisol level (upper branch of the equilibrium curve) to a second steady state with low cortisol level (lower branch of the equilibrium curve) in response to a temporary stressor. If the abscissa of the limit point LP2 lies below the basal stress level, say for instance 0.1, the return to the upper branch is hampered and the HPA axis is trapped at low cortisol as observed in hypocortisolism. Figure4A shows a corresponding cortisol concentration profile obtained by simulation of the model from Figure2.
Further explorations of the dynamic repertoire of the extended HPA model revealed a region in the parameter space which allows for an alternative bistable behavior that can be linked to the pathological condition of hypercortisolism. Figure5D provides an example of a corresponding bifurcation diagram in which the qualitative scenario of Figure5B is inverted. If a stressor drives the HPA axis from the lower cortisol branch to the upper one, it abides at elevated cortisol levels even after removal of the stressor if the abscissa of the limit point LP2 lies below the basal stress level, again say for instance 0.1. Figure5C shows a corresponding cortisol concentration profile obtained by simulation of the model depicted in Figure2. While the positive feedback mediated by the glucocorticoid receptor GR dimerization has been already suggested in as the cause for the bistable behavior associated to hypocortisolism, our in-silico pathway knock out tests on the extended HPA model indicate that the positive feedback on the genomic level mediated by the CRH hormone is vital for the bistability associated to hypercortisolism. Without this component, we could not find a parameter region giving rise to a bifurcation diagram as displayed in Figure5D. Similarly, the results of our numerical analysis suggest that the dynamic repertoire of the parsimonious model of excludes the bifurcation behavior of Figure5D, possibly due to its simplistic description of the CRH pathway.
Finally, we explored the parameter space of the extended HPA network model with objective functions for the detection of Hopf bifurcations and found parameter regions allowing for stable limit cycle oscillations, see Figure5E. The displayed temporal sequence in which the core species CRH, ACTH and cortisol reach their maximum value is in accordance with the established cascade in release of the HPA hormones. While oscillatory behavior for the parsimonious model was obtained in in response to repeated stress pulses, limit cycle oscillation for (1)-(4) were found in and for a modified version also in. However, to our knowledge the HPA model shown in Figure2 is the first one able to capture two distinct types of bistability (as reflected in Figures5B and5D) as well as autonomous oscillatory behavior.
Model based correction of Hypercortisolism
Based on the parsimonious model of Figure1 from a therapeutic strategy for the correction of hypocortisolism has been suggested in. This therapy aims at perturbing the HPA axis trapped at the low cortisol branch by a further suppression of cortisol in order to enable the return to a healthy cortisol level. An alternative idea for the treatment of hypocortisolism has been presented in which aims at an intervention with the reaction network in order to turn the dysfunctional irreversible switch into a reversible one by shifting the critical limit point in the bifurcation diagram. As the extended HPA model of Figure2 also captures a switching scenario associated to hypercortisolism, the inverse bifurcation approach of[33, 48, 48, 49] can also be applied to study which parts of the HPA network need to be targeted in order to allow the return from the abnormal hypercortisol steady state back to a healthy cortisol level.
For the purpose of illustration we suppose that a representative q0 out of the parameter region compatible with the switching behavior of Figure5D has been chosen (for instance by fitting of corresponding patient data) in which the abscissa S2(q0) of the limit point LP2(q0) lies below the basal stress level S = 0.1 (on a fictive scale). Then, a short but strong stress signal may drive the system from a normal cortisol level to an abnormal elevated one where it abides even if the stress drops back to basal level, see Figure5C. A therapeutic qualitative change from an irreversible to a reversible steady state transition can be achieved by interfering with the HPA model such that the abscissa of the limit point LP2(q0) is shifted to a value with while, e.g., maintaining the abscissa value of the second limit point LP1(q0), see Figure5A.
A corresponding interference strategy can be derived by solving the nonlinear inverse problem
where F maps the correction x onto the abscissas S1(q0 + x) and S2(q0 + x) of the limit points LP1(q0 + x) and LP2(q0 + x) defined by the HPA model with altered parameter vector q0 + x. In order to keep the number of necessary interventions with the HPA network as low as possible a sparse solution of (5) is sought for by solving the optimization problem
with φ as in (35) and the sparsity promoting choice p < 1.
Figure6C displays the solution of (6) (in relation to q0) giving rise to the re-engineered bifurcation diagram of Figure6A (green curve) and allowing the return to normal cortisol levels, see Figure6B. The parameters K6, v3 and d7 to be altered concern the feedback controls related to the hormone CRH. A decrease of v3 and an increase of d7 correspond to a reduction of the sensitivity of the cells with respect to extracellular CRH. Moreover, an increase of the dissociation/Michaelis-Menten constant K6 reduces the sensitivity of the TFs transcription factor regulation by the CRH-CRHR related pathways. Furthermore, the increase of ktl2 calls for an increase of the translation rate of GR mRNA in order to further diminish the effect of extracellular CRH due to higher availability of GR. The intervention strategy shown in Figure6C also suggests to alter the parameters K3, d3 and d2 in order to reduce the production of CRH and to increase the sensitivity with respect to the feedback via cortisol. As these parameters concern processes outside of the anterior pituitary gland cells, an extension of HPA model by a mechanistic description of the hypothalamus is necessary in order to enable a more detailed interpretation of those parameter changes.
Motivated by the experimental detection of a glucocorticoid membrane receptor GPCR we have built a mathematical model of the anterior pituitary cell that differentiates between genomic and GPCR mediated non-genomic signaling pathways. After the integration of our model into a literature framework of the hypothalamus pituitary adrenal axis and the training against stress induced ACTH concentration profiles of murine AtT-20 pituitary tumor cells we obtained a dynamic repertoire capturing two types of bistable switches as well as stable limit cycle oscillations. While the first bistable switch has been previously associated to the diseased state of hypocortisolism and utilized in the development of therapeutic strategies, to our knowledge the second bistable switch is novel in the context of HPA axis modelling. Along similar lines it can be linked to another dysfunctionality of the HPA axis referred to as hypercortisolism in which a continued abidance of the axis at elevated cortisol levels is observed even after complete removal of the stressing factor. Our in-silico knockout studies in particular indicate that the positive feedback on the genomic level mediated by the CRH hormone is responsible for the generation of the hypercortisolism switch. Using inverse bifurcation analysis with sparsity enforcing penalization we finally detected key model parameters that need to be altered in order to turn the irreversible switch into a reversible one then featuring a healthy stress response characterized by a return to normal cortisol levels.
While systems biology strongly advertises the use of modeling and simulation in the design of future drugs and therapies, one current challenge is to appropriately picture the complexity of the biological system at hand in the defining mathematical equations and their parameters. Though with respect to HPA axis modeling we reached a certain level of detail by introducing both genomic and non-genomic signaling pathways in the pituitary cell, we only adopted the parsimonious description of the link between the pituitary and the hypothalamus as well as the adrenal from the literature. Furthermore, the information content of the available ACTH concentration did not allow to sufficiently constrain the model parameters such that we so far have to content ourselves with parameter space regions corresponding to different qualitative model behaviour. Nevertheless, even at its current stage our model offers to endocronologists a different view on hypercortisolism as an irreversible bistable switch while the outlined mathematical methodologies will also apply to more realistic models and improved data sets. Finally, the model may also serve as a tool for the design of a model predictive control based treatment of hypercortisolism in which the HPA axis is driven by a - though at first glance counter-intuitive - additional administration of cortisol to a point where the only stable condition in proximity is the original healthy state, compare with.
In order to distinguish between slow genomic and fast non-genomic processes we use the modeling concept of compartments and separate the nucleus from the cytoplasm of the pituitary gland cell. Moreover, we base our model on the principles of Michaelis-Menten enzyme kinetics ([50–52]) in order to deal with physical meaningful parameters.
We refrain from a detailed description of gene regulation, RNA processing or cleaving and rather concentrate on the rate limiting processes such as transport/diffusion in and out of the nucleus. We apply this modeling strategy not only to the glucocorticoid mediated feedbacks but also to the stimulating pathways involving the hypothalamic hormone CRH. CRH does not enter the pituitary gland cells but interacts with a specific membrane receptor (CRHR) only. Still, also the binding of CRH induces both a genomic and a non-genomic feedback (cf.[53–55]). As this positive feedback via CRH counteracts the negative feedback via cortisol on the genomic as well as on the non-genomic level, also the catalyzing mechanisms related to CRH have to be modeled in more detail than in for the purpose of differentiating between genomic and non-genomic effects in the anterior pituitary gland cells. Concerning the molecular mechanisms of the anterior pituitary cells, our model involves three different pathways which are associated with the three main receptors, namely the intracellular glucocorticoid receptor GR and the two membrane receptors (GPCR and CRHR) with respect to cortisol and CRH, see Figure2 and Figure3.
Extracellular species and the controls related to the adrenals and the hypothalamus
Regarding the feedback controls related to the adrenals and the hypothalamus we follow the approach of in order to insure consistency and comparability as far as possible. In particular, we consider the following ODE describing the extracellular concentration of cortisol, with an additional term reflecting the diffusion (diffusion constant kdiff1) of cortisol into the anterior pituitary cells:
Here Vex and Vin are the volumes of the extracellular compartment and the intracellular compartment. Analogously to we model the temporal development of the CRH concentration by
The GR mediated feedback control
One significant extension over in our modeling approach concerns the regulatory pathway related to the intracellular glucocorticoid receptor GR and in particular to the transcription factor, being the homodimer of the GR-COR complex, see Figure3A. The following equation describes the intracellular concentration of cortisol
The change in the amount of cortisol in the cytoplasm (CORin) is subject to the transport of the dimerized GR-COR complex into the nucleus. This is described by
which follows from the assumption of the fast binding of cortisol and GR and the quick formation of the homodimer. A detailed derivation is provided in the Additional file1.
The dimerized GR-COR complex is transported to the nucleus, where it acts as a transcription factor and binds to the respective DNA sites. In particular it binds in the promoter region of the GR gene and enhances the transcription of the GR RNA sequences
which then are transported to the cytoplasm and translated to the GR protein.
The second gene regulated by the GR-COR complex is the Proopiomelanocortin gene POMC, which yields the ACTH peptide after post-translational cleavage. While the POMC gene is up-regulated by certain transcription factors (TFs) related to the CRH stimulus, it is down-regulated by the GR-COR complex. It has been observed that these two transcription factors counteract via competitive inhibition (cf.[53, 56]), i.e. the binding of the GR-COR complex prevents the catalyzation of the POMC expression and thus down regulates the transcription. This can be described by the following equation for the transcribed RNA
The GPCR mediated feedback control
The second pathway we consider in our model is the negative feedback on the release of ACTH via the glucocorticoid membrane receptor (GPCR), see Figure3B. Concerning the production of the GPCR receptor we assume no regulation by means of other compounds and hence disregard the processes on the genomic level. With respect to the complex formation (COR-GPCR) we suppose a (quasi) equilibrium as it has been experimentally observed that two ligands bind to the receptor and exhibit positive cooperativity (cf.).
The GPCR mediated feedback control only regulates the release of ACTH from intracellular vesicles. Again this mechanism counteracts the positive stimulus by CRH which is also detected via a membrane bound receptor (CRHR). It is well known that the anterior pituitary gland cells show a fast response to exposure with cortisol and CRH and that this is achieved by signaling cascades, directly relating the stimuli to the fusing of the vesicles (c.f.[36, 37]). Consequently, we assume that the release of ACTH directly depends on the amount of formed receptor ligand complexes, where the catalysis of the CRH-receptor complex and the inhibitory effect of the cortisol-receptor complex are modeled as competitive inhibition. The following ODE describes the extracellular ACTH concentration
The CRHR mediated feedback control
The last pathway we consider concerns the feedback mechanism related to CRH and the respective membrane receptor CRHR, see Figure3C. The main effects are the stimulated release of ACTH and the enhanced expression of the POMC gene. The genomic effects of the CRH mediated feedback involve several transcription factors (the Nur factors, Tpit/Tbx19, Pitx1, cf.), which catalyze the expression of the POMC gene. However, for keeping the complexity moderate we only incorporate a single CRH related transcription pathway into our model and denote the representative factor by TFs, see Equation (12). The alternative would be to model each of the transcription pathways separately at the cost of additional equations and parameters to be identified, for details on transcription modeling at different levels see, e.g.,. Consequently, we also model the catalyzed production of TFs only in a coarse manner, where we account for the fast signaling of the binding of CRH onto the production of TFs by a hill factor, leading to
Concerning the fast response to CRH stimulation it is the common hypothesis that the release of ACTH is regulated via a signaling cascade which is induced by the binding of CRH to the CRHR receptor in the cell membrane. Signaling cascades are very fast compared to other cellular processes. In fact the response to cortisol or CRH has been observed already after a few micro seconds (cf.[40, 54, 58, 59]). Thus we assume an immediate effect of the membrane receptor complex formation on the release of ACTH, as modeled in the ODE (13).
Figure2 gives a graphical outline of the complete model which is mathematically described by means of the following non-linear ODE system with 46 modeling parameters in 15 equations
We now comment on the role of the parameters in our ODE model. We have three subsidiary parameters describing the average volumes of the three compartments Vin, Vex and Vnu. These parameters are phenomenological and are fixed in the subsequent analysis.
Another class of parameters are the constants related to diffusion and transport of the compounds. In our model kdiff1, k1ex, k2in, k2ex, k3ex and k3in comprise the list of these parameters. We do not explicitly model active transport, as this would involve additional regulatory species and mechanisms and consequently would lead to an even more complex model. Subsequently, we assume that the rate of diffusion/transport is constant on average.
The largest portion of model parameters is related to degradation processes whose rates are assumed to be proportional to the concentration of the respective species. The associated rate constants are denoted by d1, d2, …, d14 and d15.
As a consequence of the considered enzymatic reactions and ligand-receptor complexes the model includes several dissociation constants. Moreover, we have some constants arising due to the Michaelis-Menten approximations, i.e. the maximum reaction rates and the Michaelis-Menten constants. They are closely related to the dissociation constant of the respective complex and approximate them if the rate of the product formation dominates the degradation of the enzyme-substrate complex. kGRdim, kGC2, kCRC and kGRC denote the dissociation constants of the considered ligand-receptor complexes. K1, K3, …, K9 denote the Michaelis-Menten constants of the respective enzyme-substrate complexes. v1, ks, ksksb and v4 denote the corresponding limiting rates of the Michaelis-Menten approximations.
Our description of the genomic mechanisms involves rates for transcription (v5, v6, v7 and v8), translation (ktl1 and ktl2) and transport (ktrs1 and ktrs2) of mRNA. We neglect the post-processing of mRNA as our focus is on the respective timescales and as the time delay is mainly due to the translocation. The constants v2 and v3 denote the production rate of the two membrane receptors. As we have not considered their regulatory mechanism we lumped the transcription and translation rate together.
Model training against data
The central idea of variational regularization is to augment the least-squares term that measures the mismatch between the data and the model with a penalizing term φ that guarantees the continuous dependence of the parameter solution on the data
The ODE system (32) represents the HPA axis model of Figure2 along with its parameter vector q of length 46. ρ1 and ρ2 in (33) denote the initial conditions corresponding to the two different experimental settings as described above and the two resulting data sets and displayed in Figure5A. While the initial concentration of some species are known from the experimental setup and can be taken into account by the constraints (34), the remaining components of ρ1, ρ2 have to be estimated from the data as well. As the available amount of data is not sufficient to uniquely determine the model parameters, we have chosen a two step approach for finding the model components most relevant for reproducing the data. First, we applied classical Tikhonov regularization and solved (31) with p = 2 in the penalty term
using a combination of global and local optimization routines. The resulting parameter vector qtik then served as an initial guess for solving (31) with a p-value less than one. Sparsity enforcing ℓ p functionals with the choice 0 < p ≤ 1 aim at solutions whose less important components are driven to zero[60–62] and have been used in the context of parameter identification in[63, 64]. For nonlinear inverse problems, the regularization properties of (35) with 0 < p ≤ 1 have been analyzed in[65, 66]. Figure7 displays the parameter vector q∗ obtained for p = 0.9 in relation to qtik. The near zero components of q∗ correspond to reactions and pathways that can rather be neglected in fitting the experimental data, see Figure5B. In particular, the autocatalytic loop for the GR gene regulation, represented by the parameters ktrs, v7, v8, v9, d15, d11, seems negligible with respect to the fast dynamics. Furthermore, the near zero value of K6 indicates that the production of the TFs receptors is independent of CRH, i.e. the genomic part of the CRH related feedback.
Origin of experimental data
Murine anterior pituitary AtT-20 cells (passage 12-24) were seeded and maintained as previously described. After 92 hours of cell growth AtT-20 cells were exposed to doses of 10nM CRH and in addition 100nM cortisol. The supernatant was carefully removed from the cell layer after 1 minute to 1 hour of administration and centrifuged for 10 minutes with 800xg at 37°C. The concentration of ACTH-molecules in the cell-free supernatant was detected using Fluorescence Correlation Spectroscopy Immunoassay.
The AtT-20 cells (ATCC no. CCL-89) were purchased from the American Type Culture Collection (ATCC, Manassas, USA). Cell culture media and reagents were from Sigma (Sigma-Aldrich Inc., St. Louis, MI).
Alberts B, Johnson A, Lewis J, Raff M, Roberts K, Walter P: Molecular Biology of the Cell: Reference Edition. 2007, Garland Science
Melmed S, Kleinberg D, Ho K: Pituitary Physiology and Diagnostic Evaluation. Williams Textbook of Endocrinology: Expert Consult-Online and Print. Edited by: Melmed S, Polonsky KS, Larsen PR, Kronenberg HM. 2011, Philadelphia: Elsevier/Saunders, 175-228.
Habib KE, Gold PW, Chrousos GP: Neuroendocrinology of stress. Endocrinol Metab Clin North Am. 2001, 30: 695-728. 10.1016/S0889-8529(05)70208-5.
Papadimitriou A, Priftis KN: Regulation of the hypothalamic-pituitary-adrenal axis. Neuroimmunomodulation. 2009, 16: 265-271. 10.1159/000216184.
Demitrack MA, Dale JK, Straus SE, Laue L, Listwak SJ, Kruesi MJ, Chrousos GP, Gold PW: Evidence for impaired activation of the hypothalamic-pituitary-adrenal axis in patients with chronic fatigue syndrome. J Clin Endocrinol Metab. 1991, 73: 1224-1234. 10.1210/jcem-73-6-1224.
Cleare AJ: The HPA axis and the genesis of chronic fatigue syndrome. Trends Endocrinol Metab. 2004, 15: 55-59. 10.1016/j.tem.2003.12.002.
Di Giorgio A, Hudson M, Jerjes W, Cleare AJ: 24-hour pituitary and adrenal hormone profiles in chronic fatigue syndrome. Psychosom Med. 2005, 67: 433-440. 10.1097/01.psy.0000161206.55324.8a.
Jerjes WK, Peters TJ, Taylor NF, Wood PJ, Wessely S, Cleare AJ: Diurnal excretion of urinary cortisol, cortisone, and cortisol metabolites in chronic fatigue syndrome. J Psychosom Res. 2006, 60: 145-153. 10.1016/j.jpsychores.2005.07.008.
Yehuda R: Advances in understanding neuroendocrine alterations in PTSD and their therapeutic implications. Ann New York Acad Sci. 2006, 1071: 137-166. 10.1196/annals.1364.012. [http://dx.doi.org/10.1196/annals.1364.012]
Juruena MF, Cleare AJ, Pariante CM: The hypothalamic pituitary adrenal axis, glucocorticoid receptor function and relevance to depression. Rev Bras Psiquiatr. 2004, 26: 189-201. 10.1590/S1516-44462004000300009.
Gold PW, Chrousos GP: Organization of the stress system and its dysregulation in melancholic and atypical depression: high vs low CRH/NE states. Mol Psychiatry. 2002, 7: 254-275. 10.1038/sj.mp.4001032.
Varghese FP, Brown ES: The Hypothalamic-Pituitary-Adrenal Axis in major depressive disorder: a brief primer for primary care physicians. Prim Care Companion J Clin Psychiatry. 2001, 3: 151-155. 10.4088/PCC.v03n0401.
Lupien SJ, Nair NPV, Briére S, Maheu F, Tu MT, Lemay M, McEwen BS, Meaney MJ: Increased cortisol levels and impaired cognition in human aging Implication for depression and dementia in later life. Rev Neurosci. 2011, 10 (2): 91-173.
Plaschke K, Kopitz J, Mattern J, Martin E, Teschendorf P: Increased cortisol levels and anticholinergic activity in cognitively unimpaired patients. J Neuropsychiatry Clinic Neurosci. 2010, 22 (4): 433-441. 10.1176/appi.neuropsych.22.4.433.
Mu DL, Wang DX, Li LH, Shan GJ, Li J, Yu QJ, Shi CX: High serum cortisol level is associated with increased risk of delirium after coronary artery bypass graft surgery: a prospective cohort study. Critic Care. 2010, 14 (6): R238-10.1186/cc9393. [http://ccforum.com/content/14/6/R238]. [See related commentary by Kazmierski and Kloszewska,http://ccforum.com/content/15/1/102],
Hood L, Heath JR, Phelps ME, Lin B: Systems biology and new technologies enable predictive and preventative medicine. Science. 2004, 306 (5696): 640-643. 10.1126/science.1104635. [http://www.sciencemag.org/content/306/5696/640.abstract]
Butcher EC, Berg EL, Kunkel EJ: Systems biology in drug discovery. Nature Biotechnol. 2004, 22: 1253-1259. 10.1038/nbt1017.
Kuepfer L, Lippert J, Eissing T: Multiscale mechanistic modeling in pharmaceutical research and development. Advances in Systems Biology, Volume 736 of Advances in Experimental Medicine and Biology. Edited by: Goryanin II, Goryachev AB. 2012, New York: Springer, 543-561.
Lu J, August E, Koeppl H: Inverse problems from biomedicine. J Math Biol. 2012, 1-26. [http://dx.doi.org/10.1007/s00285-012-0523-z]
Nielsen K, Pociot F, Ottesen J: Bifurcation analysis of an existing mathematical model reveals novel treatment strategies and suggests potential cure for type 1 diabetes. Math Med Biol (Online). 2013
Li G, Liu B, Liu Y: A dynamical model of the pulsatile secretion of the hypothalamo-pituitary-thyroid axis. Biosystems. 1995, 35: 83-92. 10.1016/0303-2647(94)01484-O.
Liu BZ, Peng JH, Sun YC, Liu YW: A comprehensive dynamical model of pulsatile secretion of the hypothalamo-pituitary-gonadal axis in man. Comput Biol Med. 1997, 27 (6): 507-13. 10.1016/S0010-4825(97)00026-7.
Kyrylov V, Severyanova LA, Vieira A: Modeling robust oscillatory behavior of the hypothalamic-pituitary-adrenal axis. IEEE Trans Biomed Eng. 2005, 52 (12): 1977-1983. 10.1109/TBME.2005.857671.
Bairagi N, Chatterjee S, Chattopadhyay J: Variability in the secretion of corticotropin-releasing hormone, adrenocorticotropic hormone and cortisol and understandability of the hypothalamic-pituitaryadrenal axis dynamics–a mathematical study based on clinical evidence. Math Med Biol. 2008, 25: 37-63. 10.1093/imammb/dqn003.
Lenbury Y, Pornsawad P: A delay-differential equation model of the feedback-controlled hypothalamus-pituitary-adrenal axis in humans. Math Med Biol. 2005, 22: 15-33. 10.1093/imammb/dqh020. [http://imammb.oxfordjournals.org/content/22/1/15.abstract]
Vinther F, Andersen M, Ottesen J: The minimal model of the hypothalamic–pituitary–adrenal axis. J Math Biol. 2011, 63: 663-690. 10.1007/s00285-010-0384-2. [http://dx.doi.org/10.1007/s00285-010-0384-2]
Ottesen J: Mathematical modelling of the Hypothalamic-Pituritary-Adrenal glad (HPA)-axis; Including Hippocampal mechanicsms. Math Biosci. 2013, 246: 122-138. 10.1016/j.mbs.2013.08.010.
Savic D, Jelic S: A mathematical model of the hypothalamo-pituitary-adrenocortical system and its stability analysis. Chaos, Solitons & Fractals. 2005, 26 (2): 427-436. 10.1016/j.chaos.2005.01.013. [http://www.sciencedirect.com/science/article/pii/S0960077905000810]
Jelic S, Cupic Z, Kolar-Anic L: Mathematical modeling of the hypothalamic-pituitary-adrenal system activity. Math Biosci. 2005, 197 (2): 173-187. 10.1016/j.mbs.2005.06.006.
Gupta S, Aslakson E, Gurbaxani BM, Vernon SD: Inclusion of the glucocorticoid receptor in a hypothalamic pituitary adrenal axis model reveals bistability. Theor Biol Med Model. 2007, 4: 8-10.1186/1742-4682-4-8.
Ben-Zvi A, Vernon SD, Broderick G: Model-based therapeutic correction of hypothalamic-pituitary-adrenal axis dysfunction. PLoS Comput Biol. 2009, 5: e1000273-10.1371/journal.pcbi.1000273. [http://dx.doi.org/10.1371%2Fjournal.pcbi.1000273]
Sriram K, Rodriguez-Fernandez M, Doyle FJ: Modeling cortisol dynamics in the neuro-endocrine axis distinguishes normal, depression, and post-traumatic stress disorder (PTSD) in humans. PLoS Comput Biol. 2012, 8 (2): e1002379-10.1371/journal.pcbi.1002379.
Engl HW, Flamm C, Kügler P, Lu J, Müller S, Schuster P: Inverse problems in systems biology. Inverse Problems. 2009, 25 (12): 123014-10.1088/0266-5611/25/12/123014. [http://stacks.iop.org/0266-5611/25/i=12/a=123014]
Tsai SY, Carlstedt-Duke J, Weigel NL, Dahlman K, Gustafsson JA, Tsai MJ, O’Malley BW: Molecular interactions of steroid hormone receptor with its enhancer element: evidence for receptor dimer formation. Cell. 1988, 55: 361-369. 10.1016/0092-8674(88)90059-1.
Drouin J, Sun YL, Tremblay S, Lavender P, Schmidt TJ, de Lean A, Nemer M: Homodimer formation is rate-limiting for high affinity DNA binding by glucocorticoid receptor. Mol Endocrinol. 1992, 6: 1299-1309. 10.1210/me.6.8.1299.
Norman AW, Mizwicki MT, Norman DP: Steroid-hormone rapid actions membrane receptors and a conformational ensemble model. Nat Rev Drug Discov. 2004, 3: 27-41. 10.1038/nrd1283.
Losel R, Wehling M: Nongenomic actions of steroid hormones. Nat Rev Mol Cell Biol. 2003, 4: 46-56. 10.1038/nrm1009.
Dayanithi G, Antoni FA: Rapid as well as delayed inhibitory effects of glucocorticoid hormones on pituitary adrenocorticotropic hormone release are mediated by type II glucocorticoid receptors and require newly synthesized messenger ribonucleic acid as well as protein. Endocrinol. 1989, 125: 308-313. 10.1210/endo-125-1-308.
Redei E, Li L, Halasz I, McGivern RF, Aird F: Fast glucocorticoid feedback inhibition of ACTH secretion in the ovariectomized rat: effect of chronic estrogen and progesterone. Neuroendocrinology. 1994, 60: 113-23. 10.1159/000126741. [PMID: 7969768]
Tasker JG, Di S, Malcher-Lopes R: Minireview: rapid glucocorticoid signaling via membrane-associated receptors. Endocrinology. 2006, 147 (12): 5549-56. 10.1210/en.2006-0981.
Lightman SL: Patterns of exposure to glucocorticoid receptor ligand. Biochemical Society Transactions. 2006, 34: 1117-8. [PMID: 17073764]
Stellato C: Post-transcriptional and nongenomic effects of glucocorticoids. Proc Am Thorac Soc. 2004, 1 (3): 255-263. 10.1513/pats.200402-015MS.
Maier C, Runzler D, Schindelar J, Grabner G, Waldhausl W, Kohler G, Luger A: G-protein-coupled glucocorticoid receptors on the pituitary cell membrane. J Cell Sci. 2005, 118 (Pt 15): 3353-3361.
Puchinger M, Zarzer C, Kügler P, Gaubitzer E, Kohler G: In vitro detection of adrenocorticotropic hormone levels by fluorescence correlation spectroscopy immunoassay for mathematical modeling of glucocorticoid-mediated feedback mechanisms. EURASIP J Bioinform Sys Biol. 2012, 2012: 17-10.1186/1687-4153-2012-17. [http://bsb.eurasipjournals.com/content/2012/1/17]
Engl HW, Hanke M, Neubauer A: Regularization of Inverse Problems, Volume 375 of Mathematics and its Applications. 1996, Dordrecht: Kluwer Academic Publishers Group
Chickarmane V, Paladugu SR, Bergmann F, Sauro HM: Bifurcation discovery tool. Bioinformatics. 2005, 21 (18): 3688-3690. 10.1093/bioinformatics/bti603. [http://bioinformatics.oxfordjournals.org/content/21/18/3688.abstract]
Lu J: Inverse eigenvalue problems for exploring the dynamics of systems biology models. J Biol Eng. 2009, 1: 711-728.
Lu J, Engl HW, Schuster P: Inverse bifurcation analysis: application to simple gene systems. Algo Mol Biol. 2006, 1 (11): 3353-3361.
Lu J, Engl HW, Machné R, Schuster P: Inverse Bifurcation Analysis of a Model for the Mammalian G1/S Regulatory Module. Bioinformatics Research and Development Volume 4414 of Lecture Notes in Computer Science. Edited by: Hochreiter S, Wagner R. 2007, Heidelberg: Springer Berlin, 168-184. [http://dx.doi.org/10.1007/978-3-540-71233-614]
Bisswanger H: Enzyme Kinetics: Principles and Methods. 2002, Weinheim: Wiley-VCH
Leskovac V: Comprehensive Enzyme Kinetics. 2004, New York, Boston, Dordrecht, London, Moscow: Kluwer Academic, Publishers
Cornish-Bowden A: Fundamentals of Enzyme Kinetics. 2004, Weinheim: Portland press
Lu L, Suzuki T, Yoshikawa Y, Murakami O, Miki Y, Moriya T, Bassett MH, Rainey WE, Hayashi Y, Sasano H: Nur-related factor 1 and nerve growth factor-induced clone B in human adrenal cortex and its disorders. J Clin Endocrinol Metab. 2004, 89: 4113-4118. 10.1210/jc.2004-0069.
Aguilera G: Corticotropin releasing hormone, receptor regulation and the stress response. Trends Endocrinol Metab. 1998, 9 (8): 329-336. 10.1016/S1043-2760(98)00079-4.
Ma XM, Aguilera G: Differential regulation of corticotropin-releasing hormone and vasopressin transcription by glucocorticoids. Endocrinology. 1999, 140 (12): 5642-5650. 10.1210/en.140.12.5642.
Murakami I, Takeuchi S, Kudo T, Sutou S, Takahashi S: Corticotropin-releasing hormone or dexamethasone regulates rat proopiomelanocortin transcription through Tpit/Pitx-responsive element in its promoter. J Endocrinol. 2007, 193 (2): 279-290. 10.1677/JOE-06-0143.
Wilkinson DJ: Stochastic Modelling for Systems Biology. 2012, Boca Raton, London, New York: CRC Press
Hashimoto K, Yunoki S, Takahara J, Ofuji T: ACTH release in pituitary cell cultures. Effect of neurogenic peptides and neurotransmitter substances on ACTH release induced by hypothalamic corticotropin releasing factor (CRF). Endocrinol Jpn. 1979, 26: 103-109. 10.1507/endocrj1954.26.103.
Lowry PJ, Estivariz FE, Gillies GE, Kruseman AC, Linton EA: CRF its regulation of ACTH and pro-opiomelanocortin peptide release and its extra hypothalamic occurrence. Acta Endocrinol Suppl (Copenh). 1986, 276: 56-62.
Bruckstein AM, Donoho DL, Elad M: From sparse solutions of systems of equations to sparse modeling of signals and images. SIAM Rev. 2009, 51 (11): 34-81.
Lai MJ: On sparse solution of underdetermined linear systems. J Concrete and Appl Math. 2010, 8: 296-327.
Daubechies I, Defrise M, De-Mol C: An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Commun on Pure Appl Math. 2004, 57 (11): 1413-1457. 10.1002/cpa.20042.
Bonneau R, Reiss D, Shannon P, Facciotti M, Hood L, Baliga N, Thorsson V: The Inferelator: an algorithm for learning parsimonious regulatory networks from systems-biology data sets de novo. Genome Biol. 2006, 7 (5): R36-10.1186/gb-2006-7-5-r36. [http://genomebiology.com/2006/7/5/R36]
Kügler P, Gaubitzer E, Müller S: Parameter identification for chemical reaction systems using sparsity enforcing regularization: a case study for the chlorite-iodide reaction. J Phys Chemis A. 2009, 113 (12): 2775-2785. 10.1021/jp808792u.
Zarzer CA: On Tikhonov regularization with non-convex sparsity constraints. Inverse Problems. 2009, 25: 025006-10.1088/0266-5611/25/2/025006.
Grasmair M, Haltmeier M, Scherzer O: Sparse regularization with ℓ q penalty term. Inverse Problems. 2008, 24 (5): 055020-10.1088/0266-5611/24/5/055020. [http://stacks.iop.org/0266-5611/24/i=5/a=055020]
Financial support by the Viennese Science and Technology Funds WWTF, project grant MA07–30, is kindly acknowledged.
The authors declare that they have no competing interests.
PK and CZ wrote the manuscript, discussed and simulated the HPA model, CZ developed the HPA model, MP conducted the lab experiments, MP and GK designed the lab experiments. All authors read and approved the final manuscript.
Electronic supplementary material
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.