A new model for the HPA axis explains dysregulation of stress hormones on the timescale of weeks

A model for HPA axis dynamics that includes functional mass changes

We begin by showing that the classic HPA model cannot produce the observed dysregulation. We then add new equations for the gland masses and show that they are sufficient to explain the dysregulation.

The classical understanding of the HPA model is described by several minimal models (Gupta et al, 2007; Sriram et al, 2012; Andersen et al, 2013; Bangsgaard & Ottesen, 2017). These models are designed to address the timescale of hours to days and capture HPA dynamics on this timescale including circadian and ultradian rhythms. The input to these HPA models is the combined effects of physical and psychological stresses, including low blood glucose, low blood pressure, inflammation signals, psychological stressors, or effects of drugs such as alcohol. All inputs acting at a given time‐point are considered as a combined input signal which we denote as u. The concentration of the three hormones CRH, ACTH, and cortisol are x ₁, x ₂, x ₃. The three‐hormone cascade, with feedback by x ₃, is described by

where the hormone secretion parameters are k ₁, k ₂, k ₃, and hormone removal rates are w ₁, w ₂, w ₃. Hormone half‐lives, given by log 2/w _i, are 4 minutes for x ₁, 20 minutes for x ₂, and 80 min for x ₃ (Bingzheng et al, 1990; Table 1). The feedback functions g ₁(x ₃) and g ₂(x ₃) are Hill functions which describe the negative effect of cortisol on secretion of x ₁ and x ₂ (see Materials and Methods).

These equations have a single stable steady‐state solution (Andersen et al, 2013). Their response to prolonged stress, namely a pulse of input u(t) that lasts for a few weeks, shows elevated hormones during the stress, and a return to baseline within hours after the stress is over. Figure 2A shows a simulated CRH test, in which external CRH is added at a given time‐point. ACTH shows no blunted ACTH responses in either early, intermediate, or late withdrawal phases (Fig 2A). This behavior can also be shown analytically (Materials and Methods).

To describe HPA dynamics over weeks, we add to the classic model two interactions between the hormones and the total functional mass of the cells that secrete these hormones. We introduce two new variables, the functional mass of the corticotrophs, C(t), and the functional mass of the adrenal cells that secrete cortisol, A(t). To focus on the role of the mass, we separate the secretion parameters into a product of secretion per cell times the total cell mass. The secretion parameter of ACTH is thus k ₂ = b ₂C, where C is the corticotroph mass and b ₂ is the rate of ACTH secretion per unit corticotroph mass. The parameter b ₂ includes the metabolic capacity of the corticotrophs, the number of CRH receptors and the total blood volume which dilutes out ACTH. A similar equation describes the secretion parameter of cortisol, k ₃ = b ₃ A. To isolate the effects of mass changes, we assume for simplicity that the per‐unit‐biomass secretion rates b ₂ and b ₃ are constant, whereas A(t) and C(t) can vary with time.

This introduces two new equations for the functional masses, which have a slow timescale of weeks. Corticotrophs proliferate under control of x ₁, and adrenal cortex cells under control of x ₂, and thus

where the cell‐mass production rates are k _C x ₁ and k _A x ₂, and the cell‐mass removal rates are w _C and w _A. The removal rates are taken from experimental data that indicate cell half‐lives of days–weeks (Swann, 1940; Westlund et al, 1985; Gulyas et al, 1991). We use half‐life of 6 days for C and 12 days for A (parameters given in Table 1), but the results do not depend sensitively on these parameters, as shown below.

The new model thus has five equations, three on the fast timescale of hours and two on the slow timescale of weeks. One can prove that they have a single stable steady state. We non‐dimensionalized the equations to provide hormone and cell‐mass steady‐state levels of 1 (Materials and Methods; equations (1), (2), (6), (7), (8)).

Model shows HPA axis dysregulation after prolonged activation

The HPA model with gland‐mass changes captures the experimentally observed dysregulation (Fig 2A and B). In response to a prolonged stress input u of a few weeks or longer, the adrenal gland‐mass A and corticotroph mass C both grow. When the prolonged stress input is over, the glands are large (Fig 2B—early withdrawal). They gradually return to baseline, but the corticotroph mass returns with an undershoot, dropping below baseline mass and then returning over several weeks (Fig 2B). These qualitative properties of the dynamics can be shown analytically (Materials and Methods) and do not depend on model parameters. These transients of gland masses are at the core of the hormone dysregulation.

The changed masses of the glands affect the response of the hormones to a CRH test. In early withdrawal, cortisol is high and ACTH responses are blunted (Fig 2B). Then, for a period of several weeks in intermediate withdrawal, cortisol has returned to its original baseline but ACTH remains blunted. Finally, after several months, both cortisol and ACTH return to their original baselines. Thus, the model recapitulates the experimentally observed dynamics of Fig 1.

To understand these dynamics in detail, we plot in Fig 3 the full behavior of the functional masses and hormones during and after the prolonged pulse of input (Fig 3A). Importantly, the qualitative conclusions are insensitive to the precise values of the model parameters. Figure 3 shows the dynamics of corticotroph and adrenal masses (Fig 3B), as well as the hormone levels (Fig 3C). To compare to CRH tests, we also simulated a CRH test at each time‐point (Fig 3D). To visualize the result of the CRH tests, we plot the ratio of the peak hormone level after CRH administration relative to a control CRH test without the prolonged stress input (Fig 3E). Blunted responses correspond to values less than 1.

One can see two phases during the stress pulse. The initial phase occurs after the onset of the stressor and before adaptation to the stressor (Marked ONSET in Fig 3). It lasts several weeks. In this phase, the increase in input u causes elevated levels and responses of CRH, ACTH, and cortisol (Fig 3C). However, over weeks of stress input, the corticotroph and adrenal masses grow (Fig 3B). These gland masses thus effectively adjust to the stressor, as example of the more general phenomenon of dynamical compensation in physiological systems (Karin et al, 2016). The mass growth causes a return to baseline of hypothalamic CRH and ACTH, due to negative feedback by cortisol. More precisely, a larger adrenal functional mass means that less ACTH is needed to produce a concentration of cortisol that drives ACTH down to baseline.

Such a return to baseline is called exact adaptation. Exact adaptation is a robust feature of this circuit due to a mathematical principle in the functional mass equations, equations ((4) and (5)), called integral feedback (Karin et al, 2016; Materials and Methods): The only steady‐state solution of equations ((4) and (5)) is that the hormones x ₁ and x ₂ balance proliferation and growth parameters. Exact adaptation does not occur in models without the effects of functional mass changes—the hormones do not adapt to the stressor (Appendix Fig S1, Appendix Section 1).

The enlarged functional masses result in elevated cortisol levels during the stress period, but in adapted (that is, baseline) levels of CRH, ACTH, and blunted responses of CRH and ACTH to inputs (Fig 3E). During prolonged stress, there is thus a transition from an elevated to a blunted response of ACTH that occurs due to changes in functional masses, and results from cortisol negative feedback.

At the end of the prolonged stress pulse, the early withdrawal (or EW) phase, the functional masses are abnormal and take weeks to months to recover. This fundamental process is the reason for the hormonal dysregulation that is the subject of this study. In the first weeks after the stressor is removed, the adrenal and corticotroph functional masses shrink, accompanied by dropping cortisol levels. ACTH responses are blunted, and blunting may even worsen over time.

Then, cortisol and CRH levels and responses simultaneously normalize. This marks the beginning of the next phase, intermediate withdrawal (IW, Fig 3). In this phase, ACTH responses remain blunted, despite the fact that cortisol is back to baseline, because adrenal functional mass is enlarged, and corticotroph functional mass is deficient. Finally, over time, the entire dynamics of the HPA axis normalize (Late Withdrawal, LW in Fig 3), and the system has fully recovered.

These recovery phases are robust features of the HPA model with mass dynamics. After withdrawal of the stressor, cortisol and CRH levels and dynamics recover together, before the recovery of ACTH. This order occurs regardless of parameter values such as turnover times of the tissues (proof in Materials and Methods). The intuitive reason for this is that before CRH returns to baseline, the growth rate of the pituitary corticotrophs is negative, preventing ACTH from returning to baseline.

We conclude that the gland‐mass model is sufficient to explain the dynamics of recovery from chronic HPA activation in the conditions mentioned in the introduction—anorexia, alcohol addiction, and pregnancy (Fig 3). In order to explain the timescales of recovery, the only model parameters that matter are the tissue turnover times. Good agreement is found with turnover times on the scale of 1–3 weeks for the corticotrophs and adrenal cortex cells (see Appendix Fig S2 for comparison of different turnover times, Appendix Section 2). The model therefore explains how ACTH responses remain blunted despite normalization of cortisol baseline and dynamics.

We also tested several alternative mechanisms with a slow timescale of weeks. We tested models with constant gland masses and the following processes to which we assigned time constants on the order of 1 month: GR resistance following HPA activation (Appendix Fig S1, purple lines), slow changes in input signal (Appendix Fig S1, blue lines), and slow changes in cortisol removal rate (Appendix Fig S1, gray lines). None of these putative models show the dysregulation that we consider. The reason is that these slow processes do not cause a mismatch between ACTH and cortisol needed to capture the blunting of ACTH despite normal cortisol responses during intermediate withdrawal.

Deficient GR feedback exacerbates HPA dysregulation following prolonged stress

We next considered the role of glucocorticoid receptor (GR) in recovery from prolonged HPA activation. GR mediates the negative feedback of cortisol on CRH and ACTH secretion. Impaired feedback by GR is observed in many cases of depression. For example, administration of dexamethasone, which binds the GR in the pituitary, fails to suppress cortisol secretion in the majority of individuals suffering from depression (Coppen et al, 1983). Reduced expression of GR and impaired GR function in people with depression was also demonstrated in post mortem brains (López et al, 1998; Webster et al, 2002; McGowan et al, 2009; Pandey et al, 2013) and in peripheral tissues (Pariante, 2004). The feedback strength of the GR is regulated epigenetically and is affected by early‐life adversity (Weaver et al, 2004; McGowan et al, 2009).

The relation between depression and impaired GR function seems paradoxical, since GR signaling mediates many of the detrimental effects associated with high cortisol levels such as hippocampal atrophy (Sapolsky et al, 1985). One explanation for the association between impaired GR feedback and depression is that impaired GR feedback leads to failure of the HPA axis to terminate the stress response on the timescale of hours, leading to excessive cortisol levels (De Kloet et al, 2005; Anacker et al, 2011). GR feedback also plays a role in ultradian and circadian rhythms in the HPA axis (Stavreva et al, 2009; Walker et al, 2010; Sriram et al, 2012).

Here, we present an additional explanation for the association between GR feedback and depression, based on the functional mass dynamics model. Model simulations show that negative feedback by GR protects the HPA axis against large changes in hormone levels and responses on the timescale of weeks after a prolonged stress input (Fig 4A). The impaired dynamics of cortisol and ACTH is reduced when GR affinity is high. In other words, dysregulation is more severe the weaker the feedback from GR.

The reason for this effect is as follows. After an increase in stress levels (u ₁ → u ₂), the adrenal gland mass increases (Fig 4B). When GR feedback is weak (K _GR ≫ u ₂), the adrenal increases by about a factor of u ₂/u ₁. However, if GR feedback is strong (K _GR ≪ u ₂), the adrenal increases to a smaller extent, because less cortisol is required to inhibit ACTH to the level required for precise adaptation. The smaller adrenal mass means a smaller dysregulation of cortisol and ACTH after stress. Strong GR feedback therefore provides resilience to the HPA axis against stress on the slow timescale.

We also tested the effects of variation in all other model parameters. We find that the hormone production and removal rates affect baseline levels but not the dynamics on the timescale of weeks. The parameters that do affect the timescale of weeks are the cell turnover times. We find that the present results are found for wide variations in these parameters (Appendix Fig S2, Appendix Section 2). We also tested other putative circuit designs in which the hormones control mass changes in alternative ways and find that the interactions considered here are among the very few designs that provide the observed HPA behavior (Appendix Section 3).

Mass changes provide robustness and dynamic compensation to the HPA axis

Finally, we discuss potential advantages provided by the functional mass changes in the HPA axis. The first advantage is organ size control for the corticotrophs and adrenal cortex. The circuit provides a way for the cell populations to maintain their steady‐state mass, by balancing growth and removal. Growth and removal must be precisely balanced in order to avoid aberrant growth or atrophy of the tissue. Thus, although the cells turn over, the feedback through the hormones couples with the hormone trophic effects to set the masses at a steady functional level.

The second beneficial feature is the robustness of the steady‐state hormone levels with respect to physiological parameters. The simple form of the equations for the masses C(t) and A(t) (equations (4) and (5)) locks the hormones CRH and ACTH into a unique steady state determined only by the growth and removal parameters of the tissues. Thus ${\text{CRH}}_{\mathit{st}}=\frac{w_C}{k_C}$, and ${\text{ACTH}}_{\mathit{st}}=\frac{w_A}{k_A}$. This is remarkable because all other model parameters do not affect these steady states. Similarly, the cortisol baseline level, ${\text{CORT}}_{\mathit{st}}\approx \frac{k_1{k}_C}{w_1{w}_C}u$, depends on relatively few parameters. Cortisol baseline is proportional to the input u to the hypothalamus (averaged over weeks), which corresponds to physiological and psychological stresses. Cortisol in the present analysis is the only HPA hormone that does not return to baseline after prolonged stress input. This may correspond to its physiological role as a stress response hormone, with multiple physiological endpoints, including modulation of immunity, metabolism, and behavior (McEwen, 1998). Cortisol therefore adjusts these endpoints according to averaged stress levels, allowing ongoing response as long as stress persists. Together with the persistent responsiveness of cortisol to stress input u, the circuit allows cortisol level to be independent on most other parameters, including the hormone production rates per unit biomass $(b_2,b_3)$ or removal rates $(w_2,w_3)$, as well as the rates of the proliferation and removal of the adrenal cortex cells, k _A and w _A. This robustness is due to the ability of the functional masses to grow and shrink to buffer changes in these parameters. It allows cortisol to respond precisely to chronic stress input despite physiological variations in many of the circuit parameters.

This robustness also makes the steady‐state level of cortisol and ACTH independent of total blood volume. This is because blood volume only enters through the production parameters b ₂, b ₃. These parameters describe the secretion rate per cell, and since hormone concentration is distributed throughout the circulation, these parameters are inversely proportional to blood volume [k ₁ relates to the portal vein system that directly connects the hypothalamus and pituitary, and not to total blood volume (Owens & Nemeroff, 1991)]. The functional masses therefore adjust to stay proportional to blood volume.

The third feature of the present model is dynamical compensation, in which masses C and A change to make the fast‐timescale (minute‐hour) response of the system to acute stress inputs u invariant to changes in the production rates b ₂, b ₃ as well as several other model parameters. A full analysis of dynamical compensation in the HPA axis is provided in the SI (Appendix Section 4).

HPA axis model

We model HPA axis using the following non‐dimensionalized equations, in which hormone and cell‐mass steady states are equal to 1:

The cortisol feedback on CRH is through the GR and MR receptors, and on ACTH only through the GR receptor. We model the effect of the high‐affinity MR receptor, $M(x_3)$, as a Michaelis–Menten function in the saturated regime, because cortisol is thought to saturate MR at physiological levels (Andersen et al, 2013). Thus, $M(x_3)=1/x_3$, where the Michaelis–Menten constant is absorbed in the production term parameters. We model the cooperative, low‐affinity GR receptor using a Hill function $G(x_3)=1/\left(1+{\left(\frac{x_3}{K_{\mathit{GR}}}\right)}^n\right)$ with n = 3, where K _GR is GR receptor halfway‐effect parameter (K _GR>>1). We model the combined receptor effects multiplicatively so that $g_1(x_3)=M(x_3)G(x_3)$, and $g_2(x_3)=G(x_3)$. All parameters are provided in Table 1. The quasi‐steady state of equations (6), (7), (8) can be solved by setting time derivatives equal to zero, yielding at basal input u = 1:

Alternative models

The model without functional mass dynamics is provided by equations (1), (2), (3). To this model, we added, instead of gland‐mass dynamics, several alternative biological processes that have a slow timescale, potentially on the order of a month. The first process is glucocorticoid resistance, where chronically elevated cortisol levels cause weaker feedback from the GR (Schaaf & Cidlowski, 2002; Cohen et al, 2012). One possible mechanism involves epigenetic effects such as DNA methylation (Watkeys et al, 2018). To model GR resistance, we added a variable R that modifies the effective binding coefficient, yielding:

For R, we use the following equation, based on a model of leptin resistance by Jacquier et al (2015). It describes the decline of R with cortisol levels:

For high cortisol levels to induce strong resistance, we use the simple forms: f(x) = 1, h(x) = 1 + λx ².

The second alternative slow process is a slow decrease in input, u. We considered an exponentially decreasing input signal (Appendix Fig S1). The third alternative process is a putative slow decrease in cortisol clearance rate. To model this, we added a term C _R to the removal term in equation (8):

Because there is no well‐characterized biological process that governs removal on the scale of weeks, we use a putative description in which cortisol reduces its own removal rate by increasing C _R:

All putative processes were provided with a month timescale, by setting $w_{C_R}=w_R=\frac{\log (2)}{30}{\text{day}}^{-1}$. Simulation results are shown in the SI.

CRH test

We modeled the CRH test by adding the following equation for the concentration of externally administrated CRH, denoted x _1E:

where $\mathrm{\delta }(t)$ describes a dose D of external CRH injected at T _inj and lasting a time W:

and the removal rate is $w_{{\mathrm{CRH}}_E}$ as described in (Saphier et al, 1992). Extrinsically administrated CRH causes the pituitary to secrete ACTH, as described by adding CRH_E to the intrinsic CRH (units of CRH_E are set to have equal biological effect to CRH). Thus, equation (7) was modified to:

The dose D and pulse width W were calibrated to provide the observed mean CRH test results in non‐stressed control subjects, providing D = 20 and W = 30 min.

Proof for dysregulation of ACTH after cortisol normalization

We briefly show how the model equations (1), (2), (6), (7), (8) provides ACTH dysregulation even after cortisol normalizes. Consider a prolonged stressor, like the one presented in Fig 3A (that is, a pulse increase in the input u). Since ACTH and CRH are adapted to the stressor, within hours after withdrawal of the stressor, CRH and ACTH levels drop to below their pre‐stressor baseline, whereas cortisol remains above its pre‐stressor baseline. CRH and cortisol recover simultaneously over a few weeks because CRH dynamics depend only on cortisol and the input u. It therefore remains to be shown that ACTH does not recover before cortisol. Since ACTH recovers from below baseline x ₂ < 1, it recovers with a positive temporal derivative $\frac{\mathrm{d}{x}_2}{\mathrm{d}t}>0$. If (by negation) this happens before CRH and cortisol recover, we come to a contradiction: Since CRH is below its baseline x ₁ < 1, C has a negative derivative due to equation (1); the derivative of A is zero due to equation (8). Because ACTH levels are approximately proportional to $C^{\frac{1}{2}}{A}^{-\frac{1}{2}}$ (equation (11)), we find that ACTH has a negative temporal derivative when it crosses its baseline $\frac{\mathrm{d}{x}_2}{\mathrm{d}t}<0$ which is a contradiction. We conclude that CRH and cortisol return to baseline before ACTH does.

Cortisol can show normal fast‐timescale response despite ACTH blunting

We now use the HPA model equations (1), (2), (6), (7), (8) to show that when CRH and cortisol levels return to baseline after stress (point IW in Fig 3), their entire dynamics in response to any fast‐timescale input (such as a CRH test) normalize, even if ACTH responses have not yet normalized. The mass of the adrenal cortex and pituitary corticotrophs at baseline is denoted $A_0,{\mathrm{C}}_0$, and their size at the time‐point IW where cortisol and CRH first normalize (x ₃ = x ₁ = 1) is λ_A A ₀, λ_C C ₀. Because cortisol and CRH are at baseline, and both are functions of the product of gland masses AC (equations (11), (12), (13)) one obtains ${\mathrm{\lambda }}_A{\mathrm{\lambda }}_C=1$. Replacing $\widetilde{x_2}={\mathrm{\lambda }}_C{x}_2$ in equations (6), (7), (8) yields revised equations for the fast‐timescale dynamics:

Note that equations ((20) and (1)) are the same as equations (6), (7), (8) when A = A ₀ and C = C ₀, that is, they are independent of λ_A, λ_C. In addition, the initial conditions for the hormones are also independent of λ_A, λ_C, because x ₃ = x ₁ = 1 are at baseline, while x ₂, which is proportional to $C^{\frac{1}{2}}{A}^{-\frac{n_1}{1+n_1}}=C_0^{\frac{1}{2}}{A}_0^{-\frac{1}{2}}·{\mathrm{\lambda }}_C$, scales with λ_C, so that ${\widetilde{x_2}}_{\mathit{ST}}$ is also independent of λ_A, λ_C. We conclude that after cortisol and CRH return to baseline, the fast‐timescale dynamics of cortisol and CRH to any input (such as a CRH test) are equal to the dynamics before the stressor. In the simulations, after CRH and cortisol have normalized for the first time, they only deviate from baseline slightly (at most 9%), and so this consideration holds approximately for further times in the scenario of Fig 3.

Software

All simulations were performed using Python 3.7.3.