|Home | About | Journals | Submit | Contact Us | Français|
The hypothalamus–pituitary–adrenal (HPA) system is closely related to stress and the restoration of homeostasis. This system is stimulated in the second half of the night, decreases its activity in the daytime, and reaches the homeostatic level during the late evening. In this paper, we derive and discuss a novel model for the HPA system. It is based on three simple rules that constitute a principle of homeostasis and include only the most substantive physiological elements. In contrast to other models, its main components include, apart from the conventional negative feedback ingredient, a positive feedback loop. To validate the model, we present a parameter estimation procedure that enables one to adapt the model to clinical observations. Using this methodology, we are able to show that the novel model is capable of simulating clinical trials. Furthermore, the stationary state of the system is investigated. We show that, under mild conditions, the system always has a well-defined set-point, which reflects the clinical situation to be modeled. Finally, the computed parameters may be interpreted from a physiological point of view, thereby leading to insights about diseases like depression, obesity, or diabetes.
The hypothalamus–pituitary–adrenal (HPA) system is a neuroendocrine system that is closely linked to stress in humans. This system is responsible for a rapid response to stressful stimuli and for the return to homeostasis through complex feedback mechanisms. Cortical brain regions, e.g., hippocampus and amygdala, are connected by glutamatergic pathways to the hypothalamus . An activation of the hypothalamic neurons of the paraventricular nucleus by glutamate causes a release of corticotropin releasing hormone (CRH). By this, CRH is secreted into the hypophyseal portal circulation to reach the anterior pituitary, where it subsequently stimulates the release of adrenocorticotropic hormone (ACTH) into the circulation. ACTH stimulates the release of cortisol in the adrenals (see Fig. 1). Serum cortisol concentration has to be sufficiently regulated within a physiological range. Hypercortisolism can cause depression, diabetes, visceral obesity, or osteoporosis. Therefore, inhibition of cortisol secretion is an essential component of regulation within this system. The inhibition is partly achieved by cortisol binding to glucocorticoid receptors in the amygdala, hippocampus, hypothalamus, pituitary, and adrenals. However, it is also important to maintain cortisol concentrations above a critical threshold since low cortisol may result in a disturbed memory formation or a life-threatening adrenal crisis [2, 3].
Two types of corticoid receptors have been described based on biochemical and functional characteristics . The mineralocorticoid receptor (MR) has a specificity in selectively binding cortisol. In the brain, MR is most densely localized in hippocampal and septal neurons. While MR has a high affinity for cortisol, its efficacy in the periphery and lower brain regions is limited by 11 β hydroxysteroid dehydrogenase 2, which converts intracellularly active cortisol into inactive cortisone. Conversely, the glucocorticoid receptor (GR) is widely distributed and represents the predominant binding site for cortisol in hypothalamus, pituitary, and adrenals, as well as in organs and tissues in the periphery. However, GR binds cortisol with a lower affinity as compared to MR. These receptor characteristics complement each other and put the MR and GR in a position to modulate HPA responses. MR appears to be sensitive to low, and saturates at high, cortisol concentrations. On the other hand, GR generates its dynamics at high, while it appears to be noneffective at low cortisol concentrations. It is known from the literature that cortisol binding to GR leads to an inhibitory effect on cortisol secretion . However, it is unclear how cortisol-bound MR affects the HPA system.
The HPA system is a dynamical closed loop system, homeostatically regulated, and subject to a daily rhythm. In the second half of the night, the HPA system is stimulated during REM sleep phases . A maximum of cortisol and ACTH concentrations is attained in the early morning hours. The hormones undergo a constant decay in the daytime. However, cortisol and ACTH concentrations rise after meals or in a physical or psychological stress situation. During the first half of the night, ACTH and cortisol concentrations reach a homeostatic level (for illustration, see Fig. 2).
The question remains as to in which way the HPA system is regulated and reaches its homeostatic level. The answer to this fundamental question is of foremost interest to the treatment of various diseases. This includes, for example, defined hypo- or hypercortisolism, such as Addison’s and Cushing’s disease. Also, its close connection to energy and weight regulation (e.g., obesity, diabetes mellitus, and metabolic syndrome) and to psychological illnesses such as depression is a profound motivation to investigate the mechanisms of HPA regulation . From an evolutionary point of view, it is widely accepted that a simple and durable mechanism is necessary to provide a basis for a homeostatic regulated system. This led us to the following principle of homeostasis for the HPA system. Its interactions eventually result in a homeostatic state of cortisol .
Note that rule 3 constitutes a positive feedback of MR on cortisol. In this way, we introduce a novel aspect in the concept of homeostasis.
Here, we present for the first time the derivation, the theoretical analysis, and a parameter estimation procedure of our mathematical HPA model. It should be noted that our findings are supported by clinical investigations. More precisely, we study various clinical interventions on different subjects and analyze the physiological impact of the estimated parameters in the work of Peters et al. .
The remaining part of this paper is organized as follows. Based on the three rules, we develop a “homeostatic” mathematical model and prove in a strict sense that this system reaches a stable state over time, that is, in the present context, the system reaches a homeostatic state or set point. Next, we establish a parameter estimation procedure, which we bring forward to adapt the new model to clinical data. Finally, we discuss its clinical relevance and conclude.
While quite a number of mathematical models of glucose metabolism have been developed and published [9–11], only a small number of models of the HPA system can be found with rather different aims. The focus of these HPA models varies from the influence of the inner clock to the internal dynamics in this system [12–14]. Despite several dissimilarities, certain characteristics of the models are modeled along the same lines. The stimulation of CRH via ACTH on cortisol is constructed identically in previous HPA models. Additionally, the degradation rates of CRH, ACTH, and cortisol are assumed to be linear. Interestingly, almost all previous HPA models use only pure negative feedback elements [15–17]. The absence of positive feedback elements in these models turns out to be a shortcoming, since there exist data indicating a positive stimulus in the HPA system . We will now present a new model, which includes positive as well as negative feedback elements and which is based on the three rules of homeostasis. Following the first rule of homeostasis (see box), cortisol (Z) binds to the high affinity MR (R1) and to the low affinity GR (R2). Cortisol and MR form a ligand-receptor complex C1, and cortisol and GR a complex C2, respectively. Denoting the rates of reaction by k1,k−1,k2 and k−2, the reaction equations can be written as
Under the common assumption of biochemical reaction kinetics and the law of mass action of competitive bindings (see ), we arrive at the following dose-response relations
for C1 and C2, respectively. Here, z denotes the concentration of Z in the chemical equilibrium. The coefficients e1 and e2 represent the integrated maximal efficacies of all MRs and GRs localized in different brain regions (e.g., hippocampus, amygdala, hypothalamus, and pituitary), while
reflect the binding affinity of MR and GR, respectively (note Fig. 3). While the affinity of MR is higher than the affinity of GR, the inequality K1<K2 holds. Following the second rule of homeostasis (see box), we add the stimulation via the cortisol–MR complex and the inhibition of the cortisol–GR complex
We consider h as a feedback of cortisol in the HPA system (see Fig. 3). We pool the integrated influences of glutamate, CRH, and ACTH from different brain regions in one molecular cue of the brain/pituitary compartment named Y (with a physiological interpretation as plasma ACTH) and suppose that the concentration y gives a positive stimulus on Z. This pooling of the CRH and ACTH compartments can be physiologically justified by the strong and fast synchrony between these hormones . Additionally, Drolet et al. showed that effects on CRH may be directly observed through the concentration of ACTH . Accordingly, several mathematical models of the HPA system either treat ACTH as a transfer hormone [16, 22, 23] or even neglect the CRH compartment . By denoting the compartment transition of the adrenal cortex to the brain/pituitary compartment by a constant , we obtain the time-dependent differential equation
where we assume Y to have a linear degradation rate . The function models an external input, e.g., infusion of ACTH. With a linear degradation rate of Z and a linear stimulation rate of Y on Z, we arrive at the closed system of two ordinary differential equations (ODE)
with some user-prescribed initial conditions . Figure 4 displays the system of differential equations 3 schematically in a two-compartment model (brain/pituitary and adrenals). By construction, the mathematical model satisfies the three physiological rules of homeostasis (see box). In addition, each parameter in our model has a clear physiological interpretation. The physiological term of homeostasis can be mathematically interpreted as an equilibrium point that is asymptotically stable. A point x∞ is called an equilibrium point for the differential equation if f (x∞)=0 for all t[t0,∞). Moreover, x∞ is called Lyapunov stable if for every neighborhood U(x∞) there is a neighborhood such that every solution x starting in x(t0)V remains in U for all t≥t0; otherwise, it is called unstable. An equilibrium x∞ is asymptotically stable if it is Lyapunov stable and, in addition, V can be chosen such that || x(t)−x∞|| →0 as t →∞ for all x(t0)V (see ). From a physiological, as well as mathematical, point of view, it is interesting to investigate the stability of the derived differential equations 3. The following theorem gives an answer in this direction.
Theorem 1 If the system of differential equations (3), with p0, satisfies the inequality
then there exists a unique point (y∞, z∞), with y∞, z∞>0, which is asymptotically stable. Moreover, the point (0,0) is unstable.
In other words, provided the inequality (4) is satisfied, the solution of the system of differential equations settles over time to a set-point. How can the inequality (4) be interpreted? In this inequality, the addend is positive and the two other addends are negative. In comparison, the positive feedback element has to “dominate” the negative feedback element and the relation of degradation to forward stimulation . From a physiological point of view, the inequality holds if and only if a positive stimulus on ACTH and cortisol for arbitrarily small concentrations of cortisol is provided.
Proof The equilibrium points of the system of differential equations (3) (with p0) are precisely the zeros of the function
and the equation . Thus, we either have one equilibrium point or three equilibrium points given by
Now, inequality (4) states that g has a positive gradient at zero
and since g(0)=0, there exists an ε>0 so that g(z)>0 for all z(0,ε]. Since the leading coefficient of g is negative, we deduce that limz→∞g(z)=−∞. Considering the curvature behavior of g, the system of differential equations (3) (with p0) has a unique equilibrium point with y(2),z(2)>0 and
Hence, (3) (with p0) has precisely one positive equilibrium point . Next, we verify that is an asymptotic stable point of the system of differential equations (3) (with p0). Here, it is sufficient to show that the real parts of the eigenvalues of the Jacobian of the right-hand side of (3) with p0 are genuinely negative in (see ). The eigenvalues (dependent only on z) of the Jacobian
are given by
For any z, the real part of λ2 is genuinely negative . Thus, it is sufficient to show that for z=z(2). From
it is easy to deduce that
Theorem 1 guarantees—under certain conditions that are consistent with the principle of homeostasis—the uniqueness of the asymptotic stable point (y∞,z∞). Apart from this important statement, we are able to calculate the asymptotic stable point of the system of differential equations (3) analytically. Furthermore, it should be mentioned that Theorem 1 holds even in the presence of a smooth external input function p with sufficient decay, e.g., p from (6).
When unpaired, the HPA system is know to be homeostatic . However, two kinds of oscillation in this system are known from the literature : on the one hand, an ultradian rhythm with a frequency of one to three times per hour and, on the other hand, a circadian rhythm as illustrated in Fig. 2. Presumably, because of small amplitudes, the ultradian rhythm does not show any relevant impact on the concentrations of ACTH and cortisol in standard CRH tests . Savić and Jelić provide a modulation for the mathematical models of the HPA system to integrate the ultradian rhythm . The circadian rhythm of the HPA system is known to be controlled by the suprachiasmatic nucleus (SCN) . Therefore, the 24-h rhythm can be modeled by an oscillating function p(t)=b6u(t) with , , representing the modulation of the SCN . The Jacobian of this extended model is given by
The eigenvalues of are determined by λ1,2 and . Assuming inequality (4) holds, we obtain: similar results as in Theorem 1. Since the eigenvalues λ1,2 have negative real parts and λ3,4 are purely imaginary, the solution trajectories in the yz phase plane will tend toward a closed curve around the stationary point (y∞,z∞); see Fig. 5.
It should come as no surprise that the solutions y(t) and z(t) greatly vary with respect to the parameters b1,b2, b3, b4, e1, e2, K1, K2 and with respect to the external input p(t). It is the goal of this section to show that these parameters may be chosen such that the resulting mathematical model closely approximates data stemming from clinical trials. While the constants b1,b2,K1,K2 in the system of differential equations (3) are known from the literature [8, 27], the parameters b3,b4,e1,e2 are unknown. The unknown external input p(t) will be treated separately. We combine the unknown parameters e1,e2,b3, and b4 with some arbitrary initial conditions of the differential equation y0 and z0 into a parameter vector θ=(e1,e2,b3,b4,y0,z0). The variation of y and z is therefore dependent on θ, i.e., y( · ;θ) and z( · ;θ). Our goal is to solve the inverse problem: Find a vector θ min so that resulting concentrations y( · ;θ min ) and z( · ;θ min ) possess a best fit with respect to given concentrations of ACTH and cortisol given by clinical trials, among all possible vectors θ. Let and , with , be clinical data of the blood plasma concentration of ACTH and blood serum concentration of cortisol at the times ti and τj with i=1,...,k and j=1,...,, respectively. We assume that the data values and are associated with measurement errors, which are independent and normally distributed. To define a proper distance measure χ2:6→+, we make use of the standard maximum-likelihood approach (see ), that is,
Here, and denote the standard deviations, which are assumed to be given (see Fig. 8 below).
At each iteration step, the system of differential equations (3) has to be solved in order to evaluate the corresponding distance measure χ2. The required approximations of y and z are computed by means of a numerical ODE solver with respect to approximate initial conditions y(0,θ)=y0, z(0,θ)=z0. Due to the nature of the underlying equations, we make use of the Matlab function , which is based on a modified Rosenbrock formula of order 2 . The aim is now to minimize χ2 with respect to θ
To solve this unconstrained, nonlinear optimization problem, we employ a direct search method (a modified Nelder–Mead algorithm). The selected stop criterion ensures that θ min is a local minimum of χ2. In order to find the global minimum and avoid local minima, we started the process with several randomly chosen initial parameter sets and chose the best set afterwards (Fig. 6).
We performed a CRH challenge test for an evaluation of our mathematical model (see ). In this clinical trial, we injected 20 healthy subjects with a dose of g CRH per kilogram body weight at the time t=0 (4 p.m.). The blood plasma concentration of ACTH and the blood serum concentration of cortisol were measured during a time period of 4 h (t[0,240]). We define to with i=1,...,k (k=17) as the mean values of the blood plasma concentration of ACTH and to with j=1,..., (=29) as the mean values of the blood serum concentration of cortisol, while and specify the respective standard deviations (see Fig. 8).
Before we perform the process of parameter estimation, the external input p of CRH on ACTH (i.e., y) has to be identified (Fig. 7). In the considerations for estimating the external input, we only take three assumptions into account. First, the separation of external input and feedback, and secondly, the assumption that feedback on ACTH is only caused by cortisol. Thirdly, we assume that the external input p can be represented with some and α+ by the exponential function
which is a natural impulse (CRH) response (ACTH) curve (see, for instance, ). Since the effect of CRH on ACTH is unknown, the parameters b5 and α are to be estimated first. We consider y from the system of differential equations (3) without feedback, i.e., e1=e2=0. Since b1 is known from the literature, y can be calculated with some given values of y0, b5, and α. The feedback of cortisol has to be responsible for the difference of the calculated y and the measured data , namely, (separation of external input and feedback and the feedback only caused by cortisol). Therefore, the parameters b5 and α can be determined in such a way that cortisol concentrations have to correspond to the difference of y and at the best possible rate. The degradation rate b1 of ACTH is given by [min−1] . For the given data (mean values measured from the 20 subjects, as shown in Fig. 8), we computed b5=25.12 [pmol l−1] and α=0.0243 [min−1], which corresponds to a degradation time of approximately 28 min.
After these preparations, we are able to apply the process of parameter estimation to the system of differential equations (3) and the data of the CRH challenge test. As initial values y0 and z0 of the system of differential equations (3), we picked the data values , , and t1=τ1=0. The parameters b2, K1, and K2 are well understood and we took b2=28, K1=0.5 [nmol l−1], and K2=5.0 [nmol l−1], from [5, 27].
The outlined parameter estimation procedure reveals θ min =(0.1290, 0.1633, 0.0336, 2.2234, 2.5345, 148.0701) as the parameter vector with a distance . The dimension units are given by [pmol l−1 min−1] for e1,e2 and [nmol l−1 min−1] for b3,b4. It should be pointed out that it is not clear whether θ min constitutes a global minimum of the distance function χ2. However, since we started the parameter estimation procedure with arbitrary several times and is small, it is likely that θ min is the global minimum. The assumptions of Theorem 1 are fulfilled for the parameter vector θ min . The left-hand side of the inequality (4) is equal to 7.52·10−3 and, therefore, greater than zero. The asymptotic stable point is calculated to be (y∞,z∞)=(1.37, 90.67), where the units of measurement for y∞ and z∞ are [pmol l−1] and [nmol l−1], respectively.
We developed a system of differential equations (3) for plasma ACTH and serum cortisol that fulfills the postulated rules of homeostasis. We proved that, under mild assumptions, the system reaches a stable state over time. The fact that the system stabilizes in a well-defined state resembles the existence of a set-point for physiological systems that obey the principles of homeostasis. Despite its appealing simplicity, we were able to tune the parameter of the model such that it fits remarkably well to clinical data.
We deliberately abstained from adding further compartments to our mathematical model. However, it is imaginable to split the brain/pituitary compartments into two (e.g., hippocampus/amygdala/hypothalamus and pituitary) or three compartments (e.g., hippocampus/amygdala, hypothalamus, and pituitary) without changing the substantive implications of Theorem 1. Since, in clinical trials, brain glutamate and CRH signals cannot be assessed in living humans, and as three or four compartments result in more unknown parameters to estimate, we used the two-dimensional system of differential equations.
One notes that no additional signals (e.g., glutamatergic stimulus from different brain regions) are needed to generate a stable HPA system. With a loss of such stimuli from other systems or external influences, our modeled system does not collapse, as purely negative feedback systems (e.g., in [15, 16]) would. From the physiological point of view, the HPA system cannot be regarded as an isolated system. However, on the assumption of the homeostatic rules, related systems within the organism may serve as modulators and not as indispensable regulators. The positive feedback via MR is, therefore, the crucial component to self-stabilize this dynamic system. One should notice that the asymptotic stable point (y∞,z∞) could be used as an indicator for disorders such as depression, obesity, and diabetes.
The domain of initial concentration of ACTH and cortisol resulting in a homeostatic point is still unknown. Therefore, further investigations should address the question of global stability of the ODE system (3).
In , we tested a slightly modified system of differential equations referring to various clinical trials. Despite the simple structure of that system of differential equations, the process of parameter estimation performs strongly to various even pathological–clinical trials. Especially in pathological cases, our model is in high agreement with clinical data, while purely negative models, e.g., [8, 18, 30], result in a poor approximation. Thus, the presented model has a wider scope of validity than purely negative feedback systems. The approach to compute the parameter vector θ min enables us to classify the feedback parameters, to calculate the asymptotic stable point (y∞,z∞), and to identify potential pathological states in any subject individually.
Receptors with high and low affinities acting on the same ligand can be found throughout the entire organism . It therefore appears plausible that other control systems in the organism are regulated by a similar mechanism. There is physiological evidence that the rules of homeostasis and the structure of our mathematical modeling might be transferred to other homeostatic biological systems [7, 8].
This work was supported by grants (Clinical Research Group KFO-126) from the German Research Foundation.
Christian Hubold, Email: firstname.lastname@example.org.
Bernd Fischer, Email: ed.kcebeul-inu.htam@rehcsif.
Achim Peters, Email: email@example.com.