Onset, timing, and exposure therapy of stress disorders: mechanistic insight from a mathematical model of oscillating neuroendocrine dynamics

Abstract

Background

The hypothalamic-pituitary-adrenal (HPA) axis is a neuroendocrine system that regulates numerous physiological processes. Disruptions in the activity of the HPA axis are correlated with stress-related diseases such as post-traumatic stress disorder (PTSD) and major depressive disorder. In this paper, we characterize “normal” and “diseased” states of the HPA axis as basins of attraction of a dynamical system describing the inhibition of peptide hormones such as corticotropin-releasing hormone (CRH) and adrenocorticotropic hormone (ACTH) by circulating glucocorticoids such as cortisol (CORT).

Results

In addition to including key physiological features such as ultradian oscillations in cortisol levels and self-upregulation of CRH neuron activity, our model distinguishes the relatively slow process of cortisol-mediated CRH biosynthesis from rapid trans-synaptic effects that regulate the CRH secretion process. We show that the slow component of the negative feedback allows external stress-induced reversible transitions between “normal” and “diseased” states in novel intensity-, duration-, and timing-dependent ways.

Conclusion

Our two-step negative feedback model suggests a mechanism whereby exposure therapy of stress disorders such as PTSD may act to normalize downstream dysregulation of the HPA axis. Our analysis provides a causative rationale for improving treatments and guiding the design of new protocols.

Reviewers

This article was reviewed by Dr. Daniel Coombs, Dr. Yang Kuang, and Dr. Ha Youn Lee.

Open peer review

This article was reviewed by Dr. Daniel Coombs, Dr. Yang Kuang, and Dr. Ha Youn Lee. For the full reviews, please go to the “Reviewers’ comments” section.

Background

Stress is an essential component of an organism’s attempt to adjust its internal state in response to environmental change. The experience, or even the perception of physical and/or environmental change, induces stress responses such as the secretion of glucocorticoids hormones (CORT) – cortisol in humans and corticosterone in rodents – by the adrenal gland. The adrenal gland is one component of the hypothalamic-pituitary-adrenal (HPA) axis, a collection of interacting neuroendocrine cells and endocrine glands that play a central role in stress response. The basic interactions involving the HPA axis are shown in Fig. 1. The paraventricular nucleus (PVN) of the hypothalamus receives synaptic inputs from various neural pathways via the central nervous system that are activated by both cognitive and physical stressors. Once stimulated, CRH neurons in the PVN secrete corticotropin-releasing hormone (CRH), which then stimulates the anterior pituitary gland to release adrenocorticotropin hormone (ACTH) into the bloodstream. ACTH then activates a complex signaling cascade in the adrenal cortex, which ultimately releases glucocorticoids (Fig. 1 b). In return, glucocorticoids exert a negative feedback on the hypothalamus and pituitary, suppressing CRH and ACTH release and synthesis in an effort to return them to baseline levels. Classic stress responses include transient increases in levels of CRH, ACTH, and cortisol. The basic components and organization of the vertebrate neuroendocrine stress axis arose early in evolution and the HPA axis, in particular, has been conserved across mammals [1].

Fig. 1

Schematic of HPA axis. a Stress is processed in the central nervous system (CNS) and a signal is relayed to the PVN in the hypothalamus to activate CRH secretion into the hypophyseal portal system. b CRH diffuses to the pituitary gland and activates ACTH secretion. ACTH travels down to the adrenal cortex to activate cortisol (CORT) release. Cortisol inhibits both CRH and ACTH secretion to down-regulate its own production, forming a closed loop. In the pituitary gland, cortisol binds to glucocorticoid receptors (GR) (yellow box) to inhibit ACTH and self-upregulate GR production. This part of the axis comprises the PA subsystem. c Negative feedback of cortisol affects the synthesis process in the hypothalamus, which indirectly suppresses the release of CRH. External inputs such as stressors and circadian inputs also directly affect the release rate of CRH

Dysregulation in the HPA axis is known to correlate with a number of stress-related disorders. Increased cortisol (hypercortisolism) is associated with major depressive disorder (MDD) [2, 3], while decreased cortisol (hypocortisolism) is a feature of post-traumatic stress disorder (PTSD), post infectious fatigue, and chronic fatigue syndrome (CFS) [4–7]. Since PTSD develops in the aftermath of extreme levels of stress experienced during traumatic incidents like combat, sexual abuse, or life-threatening accidents, its progression may be strongly correlated with disruption of the HPA axis caused by stress response. For example, lower peak and nadir cortisol levels were found in patients with combat-related PTSD [8].

Mathematical models of the HPA axis have been previously formulated in terms of dynamical systems of ordinary differential equations (ODEs) [9–12] or delay differential equations (DDEs) [13–15] that describe the time-evolution of the key regulating hormones of the HPA axis: CRH, ACTH, and cortisol. These models [13, 14, 16] incorporate positive self-regulation of glucocorticoid receptor expression in the pituitary, which may generate bistability in the dynamical structure of the model [17]. Of the two stable equilibrium states, one is characterized by higher levels of cortisol and is identified as the “normal” state. The other is characterized by lower levels of cortisol and can be interpreted as one of the “diseased” states associated with hypocortisolism. Stresses that affect the activity of neurons in the PVN are described as perturbations to endogenous CRH secretion activity. Depending on the length and magnitude of the stress input, the system may or may not shift from the basin of attraction of the normal steady state towards that of the diseased one. If such a transition does occur, it may be interpreted as the onset of disease. A later model [16] describes the effect of stress on the HPA axis as a gradual change in the parameter values representing the maximum rate of CRH production and the strength of the negative feedback activity of cortisol. In this model, cortisol secretion patterns are assumed to depend solely on physiological changes arising from e.g., anatomical or biochemical changes in cells or tissues. Such structural-level variations can be mathematically represented by changes in physiological parameter values.

These two classes of models imply qualitatively different time courses of disease progression [16, 17]. The former suggests that the abnormal state is a pre-existing basin of attraction of a dynamical model that stays dormant until a sudden transition is triggered by exposure to trauma [17]. In contrast, the latter assumes that the abnormal state is reached by the slow development of structural changes in physiology due to the traumatic experience [16]. Although both models [16, 17] describe changes in hormonal levels experienced by PTSD patients, they both fail to exhibit stable ultradian oscillations in cortisol, which is known to play a role in determining the responsiveness of the HPA axis to stressors [18].

In this study, we consider a number of distinctive physiological features of the HPA axis that give a more complete picture of the dynamics of stress disorders and that have not been considered in previous mathematical models. These include the effects of intrinsic ultradian oscillations on HPA dysregulation, distinct rapid and slow feedback actions of cortisol, and the correlation between HPA imbalance and disorders induced by external stress. As with the majority of hormones released by the body, cortisol levels undergo a circadian rhythm, starting low during night sleep, rapidly rising and reaching its peak in the early morning, then gradually falling throughout the day. Superposed on this slow diurnal cycle is an ultradian rhythm consisting of approximately hourly pulses. CRH, ACTH, and cortisol are all secreted episodically, with the pulses of ACTH slightly preceding those of cortisol [19].

As for many other hormones such as gonadotropin-releasing hormone (GnRH), insulin, and growth hormone (GH), the ultradian release pattern of glucocorticoids is important in sustaining normal physiological functions, such as regulating gene expression in the hippocampus [20]. It is unclear what role oscillations play in homeostasis, but the time of onset of a stressor in relation to the phase of the ultradian oscillation has been shown to influence the physiological response elicited by the stressor [21].

To distinguish the rapid and slow actions of cortisol, we separate the dynamics of biosynthesis of CRH from its secretion process, which operate over very different timescales [22]. While the two processes are mostly independent from each other, the rate of CRH secretion should depend on the synthesis process since CRH peptides must be synthesized first before being released (Fig. 1 c). On the other hand, the rate of CRH peptide synthesis is influenced by cortisol levels, which in turn, are regulated by released CRH levels. We will investigate how the separation and coupling of these two processes can allow stress-induced dysregulations of the HPA axis.

The mathematical model we derive incorporates the above physiological features and reflects the basic physiology of the HPA axis associated with delays in signaling, fast and slow negative feedback mechanisms, and CRH self-upregulation [23]. Within an appropriate parameter regime, our model exhibits two distinct stable oscillating states, of which one is marked by a larger oscillation amplitude and a higher base cortisol level than the other. These two states will be referred to as normal and diseased states. Our interpretation is reminiscent of the two-state dynamical structure that arises in the classic Fitzhugh-Nagumo model of a single neuron, in which resting and spiking states emerge as bistable modes of the model [24], or in models of neuronal networks where an “epileptic brain” is described in terms of the distance between a normal and a seizure attractor in phase-space [25]. Our main result is that stress-driven transitions between normal and diseased states can arise when a two-stage negative feedback (of cortisol on CRH) mechanism is incorporated. The possibility of such transitions lead to a number of novel features in the overall system.

Methods

Models of HPA dynamics [13, 14, 16, 17, 26] are typically expressed in terms of ordinary differential equations (ODEs):

$$\begin{array}{*{20}l} \frac{\mathrm{d} C}{\mathrm{d} T} &= p_{C} I(T)f_{C}(O)-d_{C}(C), \end{array} $$

((1))

$$\begin{array}{*{20}l} \frac{\mathrm{d} A}{\mathrm{d} T} &= p_{A}Cf_{A}(OR, O)-d_{A}(A), \end{array} $$

((2))

$$\begin{array}{*{20}l} \frac{\mathrm{d} O}{\mathrm{d} T} &= p_{O}A(T)-d_{O}(O), \end{array} $$

((3))

$$\begin{array}{*{20}l} \frac{\mathrm{d} R}{\mathrm{d} T} &= p_{R}g_{R}(OR) -d_{R}(R), \end{array} $$

((4))

where C(T),A(T), and O(T) denote the plasma concentrations of CRH, ACTH, and cortisol at time T, respectively. R(T) represents the availability of glucocorticoid receptor (GR) in the anterior pituitary. The amount of cortisol bound GR is typically in quasi-equilibrium so concentration of the ligand-receptor complex is approximately proportional to the product O(T)R(T) [17]. The parameters p _α (α∈{C,A,O,R}) relate the production rate of each species α to specific factors that regulate the rate of release/synthesis. External stresses that drive CRH release by the PVN in the hypothalamus are represented by the input signal I(T). The function f _C (O) describes the negative feedback of cortisol on CRH levels in the PVN while f _A (O R,O) describes the negative feedback of cortisol or cortisol-GR complex (at concentration O(T)R(T)) in the pituitary. Both are mathematically characterized as being positive, decreasing functions so that f _A,C (·)≥0 and f A,C′(·)<0. On the other hand, the function g _R (O R) describes the self-upregulation effect of the cortisol-GR complex on GR production in the anterior pituitary [27]. In contrast to f _A,C (·), g _R (·) is a positive but increasing function of OR so that g _R (·)≥0 and g R′(·)>0. Finally, the degradation functions d _α (·) describe how each hormone and receptor is cleared and may be linear or nonlinear.

Without including the effects of the glucocorticoid receptor (neglecting Eq. 4 and assuming f _A (O R,O)=f _A (O) in Eq. 2), Eqs. 1–3 form a rudimentary “minimal” model of the HPA axis [9, 28]. If f _A,C (·) are Hill-type feedback functions dependent only on O(T) and d _α (·) are linear, a unique global stable point exists. This equilibrium point transitions to a limit cycle through a Hopf bifurcation but only within nonphysiological parameter regimes [9]. The inclusion of GR and its self-upregulation in the anterior pituitary [17] creates two stable equilibrium states of the system, but still does not generate oscillatory behavior. More recent studies extend the model (represented by Eqs. 1–4) to include nonlinear degradation [16] or constant delay to account for delivery of ACTH and synthesis of glucocorticoid in the adrenal gland [13]. These two extended models exhibit only one intrinsic circadian [16] or ultradian [13] oscillating cycle for any given set of parameter values, precluding the interpretation of normal and diseased states as bistable oscillating modes of the model.

Here, we develop a new model of the HPA axis by first adapting previous work [13] where a physiologically-motivated delay was introduced into Eq. 3, giving rise to the observed ultradian oscillations [13]. We then improve the model by distinguishing the relatively slow mechanism underlying the cortisol-mediated CRH biosynthesis from the rapid trans-synaptic effects that regulate CRH secretion. This allows us to decompose the dynamics into slow and fast components. Finally, self-upregulation of CRH release is introduced which allows for bistability. These ingredients can be realistically combined in a way that leads to novel, clinically identifiable features and are systematically developed below.

Ultradian rhythm and time delay

Experiments on rats show a 3–6 min inherent delay in the response of the adrenal gland to ACTH [29]. Moreover, in experiments performed on sheep [30], persistent ultradian oscillations were observed even after surgically removing the hypothalamus, implying that oscillations are inherent to the pituitary-adrenal (PA) subsystem. Since oscillations can be induced by delays, we assume, as in Walker et al. [13], a time delay T _d in the ACTH-mediated activation of cortisol production downstream of the hypothalamus. Equation 3 is thus modified to

$$ \frac{\mathrm{d} O}{\mathrm{d} T} = p_{O}A\left(T-T_{\mathrm{d}}\right)-d_{O}O. $$

((5))

Walker et al. [13] show that for fixed physiological levels of CRH, the solution to Eqs. 2, 4 and 5 leads to oscillatory A(T),O(T), and R(T). In order to describe the observed periodic cortisol levels in normal and diseased states, the model requires two oscillating stable states. We will see that dual oscillating states can arise within our model when the delay in ACTH-mediated activation of cortisol production is coupled with other known physiological processes that we describe below.

Synthesis of CRH

CRH synthesis involves various pathways, including CRH gene transcription and transport of packaged CRH from the cell body (soma) to their axonal terminals where they are stored prior to release. Changes in the steady state of the synthesis process typically occur on a timescale of minutes to hours. On the other hand, the secretory release process depends on changes in membrane potential at the axonal terminal of CRH neurons, which occur over millisecond to second timescales.

To model the synthesis and release process separately, we distinguish two compartments of CRH: the concentration of stored CRH within CRH neurons will be denoted C _s(T), while levels of released CRH in the portal vein outside the neurons will be labeled C(T) (Fig. 1 c). Newly synthesized CRH will first be stored, thus contributing to C _s. We assume that the stored CRH level C _s relaxes toward a target value set by the function C _∞ (O):

$$ \frac{\mathrm{d} C_{\mathrm{s}}}{\mathrm{d} T} = \frac{C_{\infty}(O)-C_{\mathrm{s}}}{T_{C}}. $$

((6))

Here, T _C is a characteristic time constant and C _∞ (O) is the cortisol-dependent target level of stored CRH. Equation 6 also assumes that the relatively small amounts of CRH released into the bloodstream do not significantly deplete the C _s pool. Note that the effects induced by changing cortisol levels are immediate as the production term C _∞ (O)/T _C is adjusted instantaneously to current cortisol levels. Our model thus does not exclude cortisol rapidly acting on the initial transcription activity, as suggested by CRH hnRNA (precursor mRNA) measurements [31]. On the other hand, the time required to reach the steady state for the completely synthesized CRH peptide will depend on the characteristic time scale constant T _C . Ideally, T _C should be estimated from measurements of the pool size of releasable CRH at the axonal terminals. To best of our knowledge, there are currently no such measurements available, so we base our estimation on mRNA level measurements. We believe this is a better representation of releasable CRH than hnRNA levels since mRNA synthesis is a further downstream process. Previous studies have shown that variations in CRH mRNA due to changes in cortisol levels take at least twelve hours to detect [32]. Therefore, we estimate $T_{C} \gtrsim 12~\text {hrs} = 720~\text {min}$. The negative feedback of cortisol on CRH levels thus acts through the production function C _∞ (O) on the relatively slow timescale T _C . To motivate the functional form of C _∞ (O), we invoke experiments on rats whose adrenal glands had been surgically removed and in which glucocorticoid levels were subsequently kept fixed (by injecting exogenous glucocorticoid) for 5–7 days [22, 33]. The measured CRH mRNA levels in the PVN were found to decrease exponentially with the level of administered glucocorticoid [22, 33]. Assuming the amount of releasable CRH is proportional to the amount of measured intracellular CRH mRNA, we can approximate C _∞ (O) as a decreasing exponential function of cortisol level O.

Secretion of CRH

To describe CRH secretion, we consider the following three factors: synaptic inputs to CRH cells in the PVN, availability of releasable CRH peptide, and self-upregulation of CRH release.

CRH secretion activity is regulated by synaptic inputs received by the PVN from multiple brain regions including limbic structures such as the hippocampus and the amygdala, that are activated during stress. It has been reported that for certain types of stressors, these synaptic inputs are modulated by cortisol independent of, or parallel to, its regulatory function on CRH synthesis activity [34]. On the other hand, a series of studies [35–37] showed that cortisol did not affect the basal spiking activity of the PVN. We model the overall synaptic input, denoted by I(T) in Eq. 1, as follows

$$ I(T)=I_{\text{base}}+I_{\text{ext}}(T), $$

((7))

where I _base and I _ext(T) represent the basal and stress-dependent synaptic input of the PVN, respectively. As the effect of cortisol on the synaptic input during stress is specific to the type of stressor [38–40], we assume I _ext(T) to be independent of O for simplicity and generality. Possible implications of a cortisol-dependent input function I _ext(T,O) on model behavior will be discussed in the Additional file 1.

The secretion of CRH will also depend upon the amount of stored releasable CRH, C _s(T), within the neuron and inside the synaptic vesicles. Therefore, C _s can also be factored into Eq. 1 through a source term h(C _s) which describes the amount of CRH released per unit of action potential activity of CRH neurons. Finally, it has been hypothesized that CRH enhances its own release [23], especially when external stressors are present. The enhancement of CRH release by CRH is mediated by activation of the membrane-bound G-protein-coupled receptor CRHR-1 whose downstream signaling pathways operate on timescales from milliseconds to seconds [41, 42]. Thus, self-upregulation of CRH release can be modeled by including a positive and increasing function g _C (C) in the source term in Eq. 1.

Combining all these factors involved in regulating the secretion process, we can rewrite Eq. 1 by replacing f _C (O) with h(C _s)g _C (C) as follows

$$ \frac{\mathrm{d} C}{\mathrm{d} T} = p_{C}I(T)h(C_{\mathrm{s}})g_{C}(C)-d_{C}C. $$

((8))

In this model (represented by Eqs. 6, 8, 2, 5, and 4), cortisol no longer directly suppresses CRH levels, rather, it decreases stored CRH availability, C _s, through Eq. 6, which in turn decreases the secretion rate of CRH. The combination h(C _s)g _C (C) in Eq. 8 indicates the release rate of stored CRH decreases when either C _s or C decrease. We assume that synaptic inputs into the CRH neurons modulate the overall release process with weight p _C .

Complete delay-differential equation model

We are now ready to incorporate the mechanisms described above into a new, more comprehensive mathematical model of the HPA axis, which, in summary, includes

• A delayed response of the adrenal cortex to cortisol (Eq. 5).

• A slow time-scale negative feedback by cortisol on CRH synthesis (through the production term C _∞ (O) in Eq. 6).

• A fast-acting positive feedback of stored and circulating CRH on CRH release (through the factor h(C _s)g _C (C) of the production term in Eq. 8).

Our complete mathematical model thus consists of Eqs. 2, 4, 5, 6, and 8. We henceforth assume f _A (O R,O)=f _A (O R) depends on only the cortisol-GR complex and use Hill-type functions for f _A (O R) and g _R (O R) [13, 14, 16, 17]. Our full theory is characterized by the following system of delay differential equations:

$$\begin{array}{*{20}l} \frac{\mathrm{d}C_{\mathrm{s}}}{\mathrm{d} T} &= \frac{C_{\infty}(O)-C_{\mathrm{s}}}{T_{C}}, \end{array} $$

((9))

$$\begin{array}{*{20}l} \frac{\mathrm{d} C}{\mathrm{d} T} &= p_{C}I(T)h(C_{\mathrm{s}})g_{C}(C)- d_{C}C, \end{array} $$

((10))

$$\begin{array}{*{20}l} \frac{\mathrm{d} A}{\mathrm{d} T} &= p_{A}C\left(\frac{K_{A}}{K_{A}+OR}\right)-d_{A}A, \end{array} $$

((11))

$$\begin{array}{*{20}l} \frac{\mathrm{d} O}{\mathrm{d} T} &= p_{O}A(T-T_{d})-d_{O}O, \end{array} $$

((12))

$$\begin{array}{*{20}l} \frac{\mathrm{d} R}{\mathrm{d} T} &= p_{R}\left(1-\frac{\mu_{R}{K_{R}^{2}}}{{K_{R}^{2}}+(OR)^{2}}\right) -d_{R}R. \end{array} $$

((13))

The parameters K _A,R represent the level of A and R at which the negative or positive effect are at their half maximum and 1−μ _R represents the basal production rate for GR when O R=0.

Of all the processes modeled, we will see that the slow negative feedback described in Eq. 9 will be crucial in mediating transitions between stable states of the system. The slow dynamics will allow state variables to cross basins of attraction associated with each of the stable states.

Nondimensionalization

To simplify the further development and analysis of our model, we nondimensionalize Eqs. 9–13 by rescaling all variables and parameters in a manner similar to that of Walker et al. [13], as explicitly shown in the Additional file 1. We find

$$\begin{array}{*{20}l} \frac{\mathrm{d} c_{\mathrm{s}}}{\mathrm{d} t} & = \frac{c_{\infty}(o)-c_{\mathrm{s}}}{t_{c}}, \end{array} $$

((14))

$$\begin{array}{*{20}l} \frac{\mathrm{d} c}{\mathrm{d} t} & = q_{0}I(t)h(c_{\mathrm{s}})g_{c}(c)-q_{2}c, \end{array} $$

((15))

$$\begin{array}{*{20}l} \frac{\mathrm{d} a}{\mathrm{d} t} & = \frac{c}{1+p_{2}(or)} -p_{3}a, \end{array} $$

((16))

$$\begin{array}{*{20}l} \frac{\mathrm{d} o}{\mathrm{d} t} & = a(t-t_{\mathrm{d}}) -o, \end{array} $$

((17))

$$\begin{array}{*{20}l} \frac{\mathrm{d} r}{\mathrm{d} t} & = \frac{(or)^{2}}{p_{4}+(or)^{2}}+p_{5} -p_{6}r, \end{array} $$

((18))

where c _s,c,a,r,o are the dimensionless versions of the original concentrations C _s,C,A,R,O, respectively. The dimensionless delay in activation of cortisol production by ACTH is now denoted t _d. All dimensionless parameters q _i ,p _i ,t _d, and t _c are combinations of the physical parameters and are explicitly given in the Additional file 1. The functions c _∞ (o), h(c _s), and g _c (c) are dimensionless versions of C _∞ (O), h(C _s), and g _C (C), respectively, and will be chosen phenomenologically to be

$$\begin{array}{*{20}l} c_{\infty}(o) = & \bar{c}_{\infty}+ e^{-b o}, \\ h(c_{\mathrm{s}}) = & 1-e^{-kc_{\mathrm{s}}}, \\ g_{c}(c) = & 1- \frac{\mu_{\mathrm{c}}}{1+(q_{1}c)^{n}}. \end{array} $$

((19))

The form of c _∞ (o) is based on the above-mentioned exponential relation observed in adrenalectomized rats [22, 33]. The parameters $\bar {c}_{\infty }$ and b represent the minimum dimensionless level of stored CRH and the decay rate of the function, respectively. The function h(c _s) describes how the rate of CRH release increases with c _s. Since the amount of CRH packaged in release vesicles is likely regulated, we assume h(c _s) saturates at high c _s. The choice of a decreasing form for c _∞ (o) implies that increasing cortisol levels will decrease the target level (or production rate) of c _s in Eq. 14. The reduced production of c _s will then lead to a smaller h(c _s) and ultimately to a reduced release source for c (Eq. 15). As expected, the overall effect of increasing cortisol is a decrease in the release rate of CRH. Finally, since the upregulation of CRH release by circulating CRH is mediated by binding to CRH receptor, g _c (c) will be chosen to be a Hill-type function, with Hill-exponent n, similar in form to the function g _R (O R) used in Eqs. 13 and 18. The parameter 1−μ _c represents the basal release rate of CRH relative to the maximum release rate and $q_{1}^{-1}$ represents the normalized CRH level at which the positive effect is at half-maximum.

Fast-slow variable separation and bistability

Since we assume the negative feedback effect of cortisol on synthesis of CRH operates over the longest characteristic timescale t _c in the problem, the full model must be studied across two separate timescales, a fast timescale t, and a slow timescale τ=t/t _c≡ε t. The full model (Eqs. 14–18) can be succinctly written in the form

$$\begin{array}{*{20}l} \frac{\mathrm{d} c_{\mathrm{s}}}{\mathrm{d} t} & = \varepsilon(c_{\infty}(o)-c_{\mathrm{s}}), \end{array} $$

((20))

$$\begin{array}{*{20}l} \frac{\mathrm{d} \textbf{x}}{\mathrm{d} t} & = \textbf{F}(c_{\mathrm{s}}, \textbf{x}), \end{array} $$

((21))

where x=(c,a,o,r) is the vector of fast dynamical variables, and F(c _s,x) denotes the right-hand-sides of Eqs. 15–18. We refer to the fast dynamics described by dx/dt=F(c _s,x) as a fast flow. In the ε→0 limit, it is also easy to see that to lowest order c _s is constant across the fast timescale and is a function of only the slow variable τ.

Under this timescale separation, the first component of Eq. 21 (Eq. 15) can be written as

$$ \frac{\mathrm{d} c}{\mathrm{d} t} = q(c_{\mathrm{s}}(\tau),I) g_{c}(c) - q_{2}c, $$

((22))

where $q(c_{\mathrm {s}}(\tau),I) \equiv q_{0}I h(c_{\mathrm {s}}(\tau))=q_{0}I (1-e^{-kc_{\mathrm {s}}(\tau)})\phantom {\dot {i}\!}$ is a function of c _s(τ) and I. Since c _s is a function only of the slow timescale τ, q can be viewed as a bifurcation parameter controlling, over short timescales, the fast flow described by Eq. 22. Once c(t) quickly reaches its non-oscillating quasi-equilibrium value defined by dc/dt=q g _c (c)−q ₂ c=0, it can be viewed as a parametric term in Eq. 16 of the pituitary-adrenal (PA) subsystem.

Due to the nonlinearity of g _c (c), the equilibrium value c(q) satisfying q g _c (c)=q ₂ c may be multi-valued depending on q, as shown in Fig. 2 a and b. For certain values of the free parameters, such as n,1−μ _c, and q ₁, bistability can emerge through a saddle-node bifurcation with respect to the bifurcation parameter q. Figure 2 b shows the bifurcation diagram, i.e., the nullcline of c defined by q g _c (c)=q ₂ c.

Fig. 2

Nonlinear g _c (c) and bistability of fast variables. a The stable states of the decoupled system in Eq. 22 can be visualized as the intersection of the two functions q g _c (c) (dashed curve) and q ₂ c (gray line). For a given Hill-type function g _c (c), Eq. 22 can admit one or two stable states (solid circles), depending on function parameters. The unstable steady state is indicated by the open circle. b Bifurcation diagram of the decoupled system (Eq. 22) with q as the bifurcation parameter. Solid and dashed segments represent stable and unstable steady states of the fast variables, respectively. L and U label basins of attraction associated with the lower and upper stable branches of the c-nullcline. Left and right bifurcation points (q _L,c _L) and (q _R,c _R) are indicated. Fixed points of c appear and disappear through saddle node bifurcations as q is varied between q _L and q _R

For equilibrium values of c lying within a certain range, the PA-subsystem can exhibit a limit cycle in (a,o,r) [13] that we express as (a ^∗(θ;c),o ^∗(θ;c),r ^∗(θ;c)), where θ=2π t/t _p(c) is the phase along the limit cycle. The dynamics of the PA-subsystem depicted in Fig. 3 indicate the range of c values that admit limit cycle behavior for (a,o,r), while the fast c-nullcline depicted in Fig. 2 b restricts the range of bistable c values. Thus, bistable states that also support oscillating (a,o,r) are possible only for values of c that satisfy both criteria.

Fig. 3

Dynamics of the oscillating PA-subsystem as a function of fixed c. a Maximum/minimum and period-averaged values of ACTH, a(t), as a function of circulating CRH. b Maximum/minimum and period-averaged values of cortisol o(t). Within physiological CRH levels, ACTH, GR (not shown), and cortisol oscillate. The minima, maxima, and period-averaged cortisol levels typically increase with increasing c. The plot was generated using dimensionless variables c, a, and o with parameter values specified in [13] and t _d=1.44, corresponding to a delay of T _d=15 min

Since in the ε→0 limit, circulating CRH only feeds forward into a,o, and r, a complete description of all the fast variables can be constructed from just c which obeys Eq. 22. Therefore, to visualize and approximate the dynamics of the full five-dimensional model, we only need to consider the 2D projection onto the fast c and slow c _s variable. A summary of the time-separated dynamics of the variables in our model is given in Fig. 4.

Fig. 4

Classification of variables. Variables of the full five-dimensional model are grouped according to their dynamical behavior. c _s(τ) is a slow variable, while x(t)=(c,a,o,r) are fast variables. Of these, (a,o,r) form the typically oscillatory PA-subsystem that is recapitulated by c. In the ε=1/t _c≪1 limit, the variable c _s(τ) slowly relaxes towards a period-averaged value 〈c _∞ (o(c))〉. Therefore, the full model can be accurately described by its projection onto the 2D (c _s,c) phase space

To analyze the evolution of the slow variable c _s(τ), we write our equations in terms of τ=ε t:

$$\begin{array}{*{20}l} \frac{\mathrm{d} c_{\mathrm{s}}}{\mathrm{d} \tau} & = (c_{\infty}(o)-c_{\mathrm{s}}), \end{array} $$

((23))

$$\begin{array}{*{20}l} \varepsilon\frac{\mathrm{d} \textbf{x}}{\mathrm{d} \tau} & = \textbf{F} (c_{\mathrm{s}}, \textbf{x}). \end{array} $$

((24))

In the ε→0 limit, the “outer solution” F(c _s,x)≈0 simply constrains the system to be on the fast c-nullcline defined by q g _c (c)=q ₂ c. The slow evolution of c _s(τ) along the fast c-nullcline depends on the value of the fast variable o(t) through c _∞ (o). To close the slow flow subsystem for c _s(τ), we fix c to its equilibrium value as defined by the fast subsystem and approximate c _∞ (o(c)) in Eq. 23 by its period-averaged value

$$ \langle c_{\infty}(c) \rangle \equiv \int_{0}^{2\pi}c_{\infty}(o^{*}(\theta;c)) {\mathrm{d} \theta\over 2\pi} = \bar{c}_{\infty} + \int_{0}^{2\pi} e^{-bo^{*}(\theta;c)}{\mathrm{d} \theta \over 2\pi}. $$

((25))

Since period-averaged values of o ^∗ increases with c, 〈c _∞ (c)〉 is a decreasing function of c under physiological parameter regimes. This period-averaging approximation allows us to relate the evolution of c _s(τ) in the slow subsystem directly to c. The evolution of the slow subsystem is approximated by the closed (c _s,c) system of equations

$$\begin{array}{*{20}l} \frac{\mathrm{d} c_{\mathrm{s}}}{\mathrm{d} \tau} & = {\langle c_{\infty}(c) \rangle-c_{\mathrm{s}}}, \end{array} $$

((26))

$$\begin{array}{*{20}l} 0 & = q_{0}h(c_{\mathrm{s}})I(t)g_{c}(c) - q_{2}c. \end{array} $$

((27))

with 〈c _∞ (c)〉 evaluated in Eq. 25. By self-consistently solving Eqs. 26 and 27, we can estimate trajectories of the full model when they are near the c-nullcline in the 2D (c _s,c)-subsystem. We will verify this in the following section.

Nullcline structure and projected dynamics

The separation of timescales results in a natural description of the fast c-nullcline in terms of the parameter q (Fig. 2) and the slow c _s-nullcline (defined by the relation c _s=〈c _∞ (c)〉 relating c _s to c) in terms of c. However, the c-nullcline is plotted in the (q,c)-plane while the c _s-nullcline is defined in the (c,c _s)-plane. To plot the nullclines together, we relate the equilibrium value of c _s, 〈c _∞ (c)〉, to the q coordinate through the monotonic relationship $q(c_{\mathrm {s}}) = q_{0}I h(\langle c_{\infty }(c)\rangle)= q_{0}I(1-e^{-k\langle c_{\infty }(c)\rangle })$ and transform the c _s variable into the q parameter so that both nullclines can be plotted together in the (q,c)-plane. These transformed c _s-nullclines will be denoted “q-nullclines.”

We assume a fixed basal stress input I=1 and plot the q-nullclines in Fig. 5 a for increasing values of k, the parameter governing the sensitivity of CRH release to stored CRH. From the form $h(\langle c_{\infty }(c)\rangle)=(1-e^{-k\langle c_{\infty }(c)\rangle })$, both the position and the steepness of the q-nullcline in (q,c)-space depend strongly on k. Figure 5 b shows a fast c-nullcline and a slow q-nullcline (transformed c _s-nullcline) intersecting at both stable branches of the fast c-nullcline. Here, the flow field indicates that the 2D projected trajectory is governed by fast flow over most of the (q,c)-space.

Fig. 5

Slow and fast nullclines and overall flow field. a The nullcline of c _s in the ε→0 limit is defined by c _s=〈c _∞ (c)〉. To plot these slow nullclines together with the fast c-nullclines, we transform the variable c _s and represent it by q through the relation q=q ₀ h(c _s). These transformed nullclines then become a function of c and can be plotted together with the fast c-nullclines. For each fixed value of c, o(t;c) is computed by employing a built-in DDE solver dde23 in MATLAB. The numerical solution is then used to approximate 〈c _∞ (c)〉 in Eq. 25 by Euler’s method. The q-nullcline shifts to the right and gets steeper as k increases. b The fast c-nullcline defined by q g _c (c)=q ₂ c (black curve) is plotted together with the slow c _s-nullcline plotted in the (q,c) plane (“q-nullcline,” blue curve). Here, two intersections arise corresponding to a high-cortisol normal (N) stable state and a low-cortisol diseased (D) stable state. The flow vector field is predominantly aligned with the fast directions toward the c-nullcline

How the fast and slow nullclines intersect controls the long-term behavior of our model in the small ε limit. In general, the number of allowable nullcline intersections will depend on input level I and on parameters (q ₀,…,p ₆,b,k,n,μ _c,t _d).

Other parameters such as q ₀, q ₁, and μ _c appear directly in the fast equation for c and thus most strongly control the fast c-nullcline. Figure 6 a shows that for a basal stress input of I=1 and an intermediate value of k, the nullclines cross at both stable branches of the fast subsystem. As expected, numerical simulations of our full model show the fast variables (a,o,r) quickly reaching their oscillating states defined by the c-nullcline while the slow variable q=q ₀ I h(c _s) remains fairly constant. Independent of initial configurations that are not near the c-nullcline in (q,c)-space, trajectories quickly jump to one of the stable branches of the c-nullcline with little motion towards the q-nullcline, as indicated by ξ _f in Fig. 6 a.

Fig. 6

Equilibria at the intersections of nullclines. a For intermediate values of k, there are three intersections, two of them representing stable equilibria. Solid red lines are projections of two trajectories of the full model, with initial states indicated by red dots and final stable states shown by black dots. The full trajectories approach the intersections of the q-nullcline (blue) and c-nullcline (black). b For large k there is only one intersection at the upper branch of the c-nullcline. Two trajectories with initial states near different branches of the c-nullcline both approach the unique intersection (black dot) on the upper branch. The scenario shown here corresponds to a Type I nullcline structure as described in the Additional file 1

Once near the c-nullcline, say when |F(c _s,x)|≪ε, the trajectories vary slowly according to Eq. 23. Here, the slow variable c _s relaxes to its steady state value while satisfying the constraint F(c _s,x)≈0. In (q,c)−space, the system slowly slides along the c-nullcline towards the q-nullcline (the ξ _s paths in Fig. 6 a). This latter phase of the evolution continues until the system reaches an intersection of the two nullclines, indicated by the filled dot, at which the reduced subsystem in c _s and c reaches equilibrium.

For certain values of k and if the fast variable c is bistable, the two nullclines may intersect within each of the two stable branches of the c-nullcline and yield the two distinct stable solutions shown in Fig. 6 a. For large k, the two nullclines may only intersect on one stable branch of the c-nullcline as shown in Fig. 6 b. Trajectories that start within the basin of attraction of the lower stable branch of the c-nullcline (“initial state 2” in Fig. 6 b) will stay on this branch for a long time before eventually sliding off near the bifurcation point and jumping to the upper stable branch. Thus, the long-term behavior of the full model can be described in terms of the locations of the intersections of nullclines of the reduced system.

Results and discussion

The dual-nullcline structure and existence of multiple stable states discussed above results from the separation of slow CRH synthesis process and fast CRH secretion process. This natural physiological separation of time scales ultimately gives rise to slow dynamics along the fast c-nullcline during stress. The extent of this slow dynamics will ultimately determine whether a transition between stable states can be induced by stress. In this section, we explore how external stress-driven transitions mediated by the fast-slow negative feedback depend on system parameters.

Changes in parameters that accompany trauma can lead to shifts in the position of the nullclines. For example, if the stored CRH release process is sufficiently compromised by trauma (smaller k), the slow q-nullcline moves to the left, driving a bistable or fully resistant organism into a stable diseased state. Interventions that increase k would need to overcome hysteresis in order to restore normal HPA function. More permanent changes in parameters are likely to be caused by physical rather than by psychological traumas since such changes would imply altered physiology and biochemistry of the person. Traumatic brain injury (TBI) is an example in which parameters can be changed permanently by physical trauma. The injury may decrease the sensitivity of the pituitary to cortisol-GR complex, which can be described by decreasing p ₂ in our model. Such change in parameter would lead to a leftward shift of the q-nullcline and an increased likelihood of hypocortisolism.

In the remainder of this work, we focus on how external stress inputs can by themselves induce stable but reversible transitions in HPA dynamics without changes in physiological parameters. Specifically, we consider only temporary changes in I(t) and consider the time-autonomous problem. Since the majority of neural circuits that project to the PVN are excitatory [43], we assume external stress stimulates CRH neurons to release CRH above its unit basal rate and that I(t)=1+I _ext(t) (I _base=1) with I _ext≥0.

To be more concrete in our analysis, we now choose our nullclines by specifying parameter values. We estimate many of the dimensionless parameters by using values from previous studies, as listed in Table 1. Of the four remaining parameters, μ _c,q ₀,q ₁, and k, we will study how our model depends on k while fixing μ _c,q ₀, and q ₁. Three possible nullcline configurations arise according to the values of μ _c,q ₀, and q ₁ and are delineated in the Additional file 1. We have also implicitly considered only parameter regimes that yield oscillations in the PA subsystem at the stable states defined by the nullcline intersections.

Table 1 Dimensionless parameter values of our full model

Given these considerations, we henceforth chose μ _c=0.6, q ₁=0.04, and q ₀=77.8 for the rest of our analysis. This choice of parameters is motivated in the Additional file 1 and corresponds to a so-called “Type I” nullcline structure. In this case, three possibilities arise: one intersection on the lower stable branch of the c-nullcline if k<k _L, two intersections if k _L<k<k _R (Fig. 6 a), and one intersection on the upper stable branch of the c-nullcline if k>k _R (Fig. 6 b). For our chosen set of parameters and a basal stress input I=1, the critical values k _L=2.50<k _R=2.54 are given by Eq. A3 in the Additional file 1.

Normal stress response

Activation of the HPA axis by acute stress culminates in an increased secretion of all three main hormones of the HPA axis. Persistent hypersecretion may lead to numerous metabolic, affective, and psychotic dysfunctions [44, 45]. Therefore, recovery after stress-induced perturbation is essential to normal HPA function. We explore the stability of the HPA axis by initiating the system in the upper of the two stable points shown in Fig. 7 a and imposing a 120 min external stress input of I _ext=0.1. The HPA axis responds with an increase in the peak level of cortisol before relaxing back to its original state after stress is terminated (Fig. 7 b). This transient process is depicted in the projected (q,c)-space in Fig. 7 a.

Fig. 7

Normal stress response. Numerical solution for the response to a 120 min external stress I _ext=0.1. a At the moment the external stress is turned on, the value of (q,c) increases from its initial stable solution at (64.4,27) to (71,27) after which the circulating CRH level c, quickly reaches the fast c-nullcline (black) before slowly evolving along it towards the slow q-nullcline (blue). After short durations of stress, the system returns to its starting point within the normal state basin. b The peaks of the cortisol level are increased during stress (red) but return to their original oscillating values after the stress is turned off

Upon turning on stress, the lumped parameter q and the slow nullcline shift to the right by 10 % since q=q ₀(1+I _ext)h(〈c _∞ (c)〉) (see Fig. 7 a). The trajectory will then move rapidly upward towards the new value of c on the c-nullcline; afterwards, it moves very slowly along the c-nullcline towards the shifted q-nullcline. After 120 min, the system arrives at the “ ×” on the c-nullcline (Fig. 7 a). Once the stress is shut off the q-nullcline returns to its original position defined by I=1. The trajectory also jumps horizontally back to near the initial q value and quickly returns to the original upper-branch stable point.

External stress induces transition from normal to diseased state

We now discuss how transitions from a normal to a diseased state can be induced by positive (excitatory) external stress of sufficient duration. In Fig. 8, we start the system in the normal high-c state.

Fig. 8

Stress-induced transitions into an oscillating low-cortisol diseased state. An excitatory external stress I _ext=0.2 is applied for 30 hrs. Here, the system reaches the new stable point set by I=1.2 before stress is terminated and the q-nullcline reverts to its original position set by I=1. a At intermediate values of 2.50<k<2.54, when two stable state arise, a transition from the normal high-cortisol state into the diseased low-cortisol state can be induced by chronic external stress. b Numerical solutions of cortisol level o(T) plotted against the original time variable T shows the transition to the low-cortisol diseased state shortly after cessation of stress. c and d If k>k _R=2.54, only the normal stable state exists. The system will recover and return to its original normal state after a transient period of low cortisol

Upon stimulation of the CRH neurons through I _ext>0, both CRH and average glucocorticoid levels are increased while the average value of c _∞ (o(t)) is decreased since c _∞ (o) is a decreasing function of o. As c _s(τ) slowly decays towards the decreased target value of 〈c _∞ (o(c))〉, h(c _s(τ)) and hence q(c _s), also decrease. As shown in Fig. 8 a, much of this decrease occurs along the high-c stable branch of the c-nullcline. Once the external stress is switched off, q will jump back down by a factor of 1/(1+I _ext). If the net decrease in q is sufficient to bring it below the bifurcation value q _L ≈64 at the leftmost point of the upper knee, the system crosses the separatrix and approaches the alternate, diseased state. Thus, the normal-to-diseased transition is more likely to occur if the external stress is maintained long enough to cause a large net decrease in q, which includes the decrease in q incurred during the slow relaxation phase, plus the drop in q associated with cessation of stress. The minimum duration required for normal-to-diseased transition should also depend on the magnitude of I _ext. The relation between the stressor magnitude and duration will be illustrated in the Additional file 1 (Fig. A4).

A numerical simulation of our model with a 30 hr I _ext=0.2 was performed, and the trajectory in (q,c)-space is shown in Fig. 8 a. The corresponding cortisol level along this trajectory is plotted in Fig. 8 b, showing that indeed a stable transition to the lower cortisol state occurred shortly after the cessation of stress.

In addition to a long-term external stress, the stable transition to a diseased state requires 2.50<k<2.54 and the existence of two stable points. On the other hand, when k>k _R=2.54, the enhanced CRH release stimulates enough cortisol production to drive the sole long term solution to the stable upper normal branch of the c-nullcline, rendering the HPA system resistant to stress-induced transitions.

The response to chronic stress initially follows the same pattern as described above for the two-stable-state case, as shown in Fig. 8 c. However, the system will continue to evolve along the lower branch towards the q-nullcline, eventually sliding off the lower branch near the right bifurcation point (indicated in Fig. 2b and Fig. A2 in the Additional file 1 by (q _R,c _R)) and returning to the single normal equilibrium state. Thus, when k is sufficiently high, the system may experience a transient period of lowered cortisol levels after chronic stress but will eventually recover and return to the normal cortisol state. The corresponding cortisol level shown in Fig. 8 d shows this recovery at T≈3400 min, which occurs approximately 1500 min after the cessation of stress.

Transition to diseased state depends on stress timing

We have shown how transitions between the oscillating normal and diseased states depend on the duration of the external stress I _ext. However, whether a transition occurs also depends on the time – relative to the phase of the intrinsic ultradian oscillations – at which a fixed-duration external stress is initiated. To illustrate this dependence on phase, we plot in Fig. 9 a and b two solutions for o(T) obtained with a 250 min I _ext=0.1 initiated at different phases of the underlying cortisol oscillation. If stress is initiated near the nadir of the oscillations, a transition to the low-cortisol diseased state occurs and is completed at approximately T=1000 min (Fig. 9 a, c). If, however, stress is initiated near the peak of the oscillations, the transition does not occur and the system returns to the normal stable state (Fig. 9 b, d). In this case, a longer stress duration would be required to push the trajectory past the low-q separatrix into the diseased state.

Fig. 9

Stress timing and transition to low-cortisol oscillating state. Cortisol levels in response to I _ext=0.1 applied over 250 min. a If stress is initiated at T=150 min, a transition to the low-cortisol diseased state is triggered. b If stress is initiated at T=120 min, the system returns to its normal high-cortisol state. Note that the first peak (marked by “$\blacktriangledown$”) during the stress in (a) is higher than the first peak in (b). c If stress is initiated at T=150 min, stress cessation and the slow relaxation along the c-nullcline during stress are sufficient to bring q just left of the separatrix, inducing the transition. d For initiation time T=120 min, q remains to the right of the separatrix, precluding the transition

As discussed earlier, an increase in period-averaged cortisol level during stress drives a normal-to-diseased state transition. We see that the period-averaged level of cortisol under increased stress is different for stress started at 120 min from stress started at 150 min.

As detailed in the Additional file 1 (Fig. A5), the amplitude of the first cortisol peak after the start of stress is significantly lower when the applied stress is started during the falling phase of the intrinsic cortisol oscillations. The difference between initial responses in o(t) affects the period-averaging in 〈c _∞ (o)〉 during external stress, ultimately influencing c _s and consequently determining whether a transition occurs. Note that this phase dependence is appreciable only when stress duration is near the threshold value that brings the system close to the separatrix between normal and diseased basins of attraction. Trajectories that pass near separatrices are sensitive to small changes in the overall negative feedback of cortisol on CRH synthesis, which depend on the start time of the stress signal.

Stress of intermediate duration can induce “reverse” transitions

We can now use our theory to study how positive stressors I _ext may be used to induce “reverse” transitions from the diseased to the normal state. Understanding these reverse transitions may be very useful in the context of exposure therapy (ET), where PTSD patients are subjected to stressors in a controlled and safe manner, using for example, computer-simulated “virtual reality exposure.” Within our model we can describe ET as external stress (I _ext>0) applied to a system in the stable low-c diseased state. The resulting horizontal shift in q causes the system to move rightward across the separatrix and suggests a transition to the high-c normal state can occur upon termination of stress.

As shown in Fig. 10 a, if stressor of sufficient duration is applied, the trajectory reaches a point above the unstable branch of the c-nullcline upon termination leading to the normal, high-cortisol state (Fig. 10 b). Since the initial motion is governed by fast flow, the minimum stress duration needed to incite the diseased-to-normal transition is short, on the timescale of minutes. However, if the stressor is applied for too long, a large reduction in q is experienced along the upper stable branch. Cessation of stress might then lower q back into the basin of attraction of the low-cortisol diseased state (Fig. 10 c). Figure 10 d shows the cortisol level transiently increasing to a normal level before reverting back to low levels after approximately 1400 min.

Fig. 10

Stress-induced transitions to high-cortisol oscillating state. a Projected 2D system dynamics when a stressor of amplitude I _ext=0.2 is applied for 9 min starting at T=120 min. c is increased just above the unstable branch (c≈20) to allow the unstressed system to cross the separatrix and transition to the normal high-c stable state. b The plot of o(T) shows the transition to the high-cortisol, high-oscillation amplitude state shortly after the 9 min stress. c A stressor turned off after 780 min (13 hrs) leaves the system in the basin of attraction of the diseased state. d Cortisol levels are pushed up but after about 1400 min relax back to levels of the original diseased state

Within our dynamical model, stresses need to be of intermediate duration in order to induce a stable transition from the diseased to the normal state. The occurrence of a reverse transition may also depend on the phase (relative to the intrinsic oscillations of the fast PA subsystem) over which stress was applied, especially when the stress duration is near its transition thresholds. For a reverse diseased-to-normal transition to occur, the decrease in c _s cannot be so large that it brings the trajectory past the left separatrix, as shown in Fig. 10 c. Therefore, near the maximum duration, stress initiated near the nadir of cortisol oscillation will be more effective at triggering the transition to a normal high-cortisol state. Overall, these results imply that exposure therapy may be tuned to drive the dynamics of the HPA axis to a normal state in patients with hypocortisolism-associated stress disorders.

Conclusions

We developed a theory of HPA dynamics that includes stored CRH, circulating CRH, ACTH, cortisol and glucocorticoid receptor. Our model incorporates a fast self-upregulation of CRH release, a slow negative feedback effect of cortisol on CRH synthesis, and a delay in ACTH-activated cortisol synthesis. These ingredients allow our model to be separated into slow and fast components and projected on a 2D subspace for analysis.

Depending on physiological parameter values, there may exist zero, one, or two stable simultaneous solutions of both fast and slow variables. For small k, the parameter that relates the amount of stored releasable CRH vesicles to its secretion rate, CRH release is weak and only the low-CRH equilibrium point arises; an individual with such k is trapped in the low-cortisol “diseased” state. For large k, only the high-CRH normal state arises, rendering the individual resistant to acquiring the long-term, low-cortisol side-effect of certain stress disorders. When only one stable solution arises, HPA dysregulation must depend on changes in parameters resulting from permanent physiological modifications due to e.g., aging, physical trauma, or stress itself [45, 46]. For example, it has been observed that older rats exhibit increased CRH secretion while maintaining normal levels of CRH mRNA in the PVN [47]. Such a change could be interpreted as an age-dependent increase in k, which, in our model, implies that aging makes the organism more resistant to stress-induced hypocortisolism. Indeed, it has been suggested that prevalence of PTSD declines with age [48, 49].

Other regulatory systems that interact with or regulate the HPA axis can also affect parameter values in our model. Gonadal steroids, which are regulated by another neuroendocrine system called the hypothalamic-pituitary-gonadal (HPG) axis, activate the preoptic area (POA) of the hypothalamus [50, 51], which in turn attenuates the excitatory effects of medial amygdala stimulation of the HPA axis [52]. Thus, low testosterone levels associated with hypogonadism would effectively increase I(t) within our model, shift the q-nullcline in the (q,c)-space, and in turn increase cortisol levels. One might consider this as a possible explanation for chronically elevated cortisol levels observed in major depressive disorder patients who suffers from hypogonadism. Although it is beyond the scope of this paper, one may further investigate role of gonadal hormones, or role of any other interacting systems, in mediating stress response by considering which parameters would be affected in our model.

Within certain parameter regimes and for intermediate k, our theory can also exhibit bistability. When two stable solutions arise, we identify the states with low oscillating levels of cortisol as the diseased state associated with hypocortisolism. Transitions between different stable states can be induced by temporary external stress inputs, implying that HPA dysregulation may develop without permanent “structural” or physiological changes. Stresses that affect secretion of CRH by the PVN are shown to be capable of inducing transitions from normal to diseased states provided they are of sufficient duration (Fig. 8).

Our model offers a mechanistic explanation to the seemingly counter-intuitive phenomenon of lower cortisol levels after stress-induced activation of cortisol production. Solutions to our model demonstrate that the negative-feedback effect of a temporary increase in cortisol on the synthesis process of CRH can slowly accumulate during the stress response and eventually shift the system into a different basin of attraction. Such a mechanism provides an alternative to the hypothesis that hypocortisolism in PTSD patients results from permanent changes in physiological parameters associated with negative-feedback of cortisol [53, 54].

We also find that external stress can induce the “reverse” transition from a diseased low-cortisol state to the normal high-cortisol state. Our results imply that re-exposure to stresses of intermediate duration can drive the system back to normal HPA function, possibly “decoupling” stress disorders from hypocortisolism.

Interestingly, we show that the minimum duration required for either type of transition depends on the time at which the stress is initiated relative to the phase of the intrinsic oscillations in (a,o,r). Due to subtle differences in cortisol levels immediately following stress initiation at different phases of the intrinsic cortisol oscillation, the different cumulative negative-feedback effect on CRH can determine whether or not a trajectory crosses the separatrix (Fig. 9). When the duration of external stress is near its threshold, normal-to-diseased state transitions are easier to induce when stress is initiated near the nadir of cortisol oscillations. Reverse diseased-to-normal transitions are more easily induced when stress is initiated near the peak of the oscillations.

In summary, our theory provides a mechanistic picture that connects ortisol dysregulation with stress disorders and a mathematical framework one can use to study the downstream effects of therapies such as brief eclectic psychotherapy (BEP) and exposure therapy (ET). Both therapies involve re-experiencing stressful situations directly or through imagination, and have been consistently proven effective as first-line treatments for PTSD symptoms [55–57]. Our results suggest that ET can directly alter and “decouple” the expression of cortisol from an underlying upstream disorder. Changes in neuronal wiring that typically occur over slower times scales is also expected after ET [58]. In our model, such changes would lead to slow variations in the basal input I(t). Thus, cortisol level may not be tightly correlated with PTSD, particularly in the context of ET.

It is important to emphasize that we modeled neuroendocrine dynamics downstream of the stress input I _ext. Understanding how the form of the stress function I _ext depends on the type of stress experienced requires a more detailed study of more upstream processes, including how hormones might feed back on these higher-brain processes. Since higher cortisol levels are found among female PTSD patients with a history of childhood abuse [59] and among PTSD patients who have experienced a nuclear accident [60], future studies of such divergent, experience-dependent dysregulation will rely on more complex input functions I _ext(t). For example, under periodic driving, complex resonant behavior should arise depending on the amplitude and frequency of the external stress I _ext(t) and the nullcline structure of the specific system. Moreover, effects of other regulatory networks that interacts with the HPA axis can be included in our model through appropriate forms of I _ext(t). For example, the effects of gonadal steroids in the stress response mentioned above can be further investigated by considering a form of I _ext(t) that is dependent on gonadal steroids level. Many other interesting properties, such as response to dexamethasone administration, can be readily investigated within our model under different system parameters.

Reviewers’ comments

Reviewer’s report 1: Daniel Coombs, University of British Columbia, Canada

Summary: I find this to be a very interesting extension of previous modeling efforts on an important neuroendocrine system. This is a great example of dynamical systems application to physiology and the interpretation of normal and disease states as attractors of the system. The findings are based squarely on the modeling and the authors related their model and findings to experimental data and propose extensions to include additional biological knowledge in the model. Potential caveats are uncertainties in parameter estimates, omission of important interactions with other physiological systems, and some question of how external forcing (positive and negative stressors) should be input to the dynamical system. These possible shortcomings are acknowledged in the manuscript and provide motivation for future work within the existing framework.

Authors’ response: We thank Dr. Daniel Coombs for his appreciation of the work.

Recommendations: Many of the results appear to depend sensitively on the parameter k. This parameter describes the regulation of CRH secretion by stored CRH (mathematically, the link between d C/d t and C _s). For example, on page 16, it is stated that you need 2.5<k<2.54 to find a certain transition to the disease state of the model. The need for precision of a few % in a parameter value might suggest that the system is not terribly robust. Can you comment on this? Is there any basis on which we might estimate k for a given individual? Could k be manipulated somehow? Also, for readability: I think it would be good to remind the reader of the precise biological meaning of this important parameter before describing its effects in the model (in the “Conclusions” section).

Authors’ response: We thank Dr. Coombs for raising the important issue of robustness of dynamical structure of our model. The estimated range of 0.04 of k (2.5<k<2.54) for bistability may be interpreted too narrow or wide depending on the range of values of k observed from experimental data. To best of our knowledge, such data does not yet exist. Although direct measurements on humans are difficult, there have been in vitro studies on rats that estimated the releasable pool size of synaptic vesicles from hippocampal neurons [ 61 ]. We believe one could design a similar experiment to estimate the size of stored releasable CRH in the PVN. We are also very much interested in knowing if k can be modulated and manipulated, as it will have important implications on our model prediction of disease onset and development.

The range of k used for analysis in the manuscript was chosen arbitrarily for the purpose of concreteness. Depending on other parameters q ₀,q ₁ , and μ _c , the range of k could be significantly greater. In particular, the range of possible k that give rise to bistability is infinite in parameter regimes that correspond to “Type II" nullcline structure. We have chosen a k value that correspond to a “Type I" nullcline for for the richness of the corresponding mathematical structure. Details on critical values k _L,R and possible nullcline configurations are discussed and illustrated (Fig. A2) in the Additional file 1.

We also thank Dr. Coombs for his suggestion to recall the biological meaning of the parameter in the Summary and Conclusion section. We have edited the section accordingly.

Minor issues: No minor issues, the paper is written well and the figures are clear.

Reviewer’s report 2: Yang Kuang, Arizona State University, United States of America

Summary: The authors carefully extended some existing models of oscillating neuroendocrine dynamics by explicitly incorporating, in the cortisol equation, a discrete time delay T _d representing the time needed for ACTH-mediated activation of cortisol production downstream of the hypothalamus. Based on a systematical analysis of this more plausible model, they developed a compelling theory of the HPA dynamics that includes stored CRH, circulating CRH, ACTH, cortisol and glucocorticoid receptor. This is a significant enhancement over the existing theories.

Authors’ response: We thank Dr. Yang Kuang for a positive summary of our manuscript.

Recommendation: To facilitate a deeper biological appreciation of the fast and slow dynamics described by the model, I suggest the authors also present a parameter table for the original parameters, including their values, units and references. Some details on estimation of these parameters are also helpful.

Authors’ response: We thank Dr. Kuang for his suggestions. The list of parameters used in our analysis and numerical simulations are now listed in Table 1 in the manuscript, with details on their estimations included in the Additional file 1.

Minor issues: None.

Reviewer’s report 3: Ha Youn Lee, Keck School of Medicine University of Southern California, United States of America

Summary: Kim et al. modified the existing cortisol dynamics model in the HPA axis, mainly to replicate stable ultradian oscillations in cortisol. By introducing the slow variable, the stored CRH, authors were able to observe ultradian oscillations in cortisol level and stress-induced transitions into a low-cortisol diseased state. The manuscript is clearly written and I recommend a publication in Biology Direct.

Authors’ response: We thank Dr. Ha Youn Lee for her positive recommendation.

Recommendations: The validity of this model can be more clearly demonstrated by comparing the model solution with the previously published data of cortisol dynamics in normal and post-traumatic stress disorder (PTSD) in references ([53] and [16]).

Authors’ response: We thank Dr. Lee for the suggestion. It would indeed be interesting to fit our model to the data in [ 53 ] and compare the estimated parameter values to those of Siriam et al. [ 16 ], as their estimations were also based on the data from [ 53 ]. As noted in the manuscript, Siriam et al. have estimated K _i , the parameter analogous to K _A in our model that represents the strength of negative feedback imparted by cortisol in the pituitary. Their estimations of K _i were consistent with an enhanced negative feedback action of cortisol (i.e. decreased K _i ) in PTSD patients, as hypothesized by Yehuda et al. [ 53 ].

Instead of fitting our model to the limited number of data (time series data of three individuals, one from each of control, depressed, and PTSD group) used in [ 16 ], we can predict how our model solution will compare to the results of [ 16 ] by analyzing the effect of varying the analogous parameter ( K _A ) on the nucllcine structure of our model. When K _A is decreased (enhanced negative feedback in the pituitary), the q-nullcline in our model shifts toward the normal, high cortisol state (on the upper branch of the c-nullcline), away from the diseased state. Thus, contrary to the hypothesis of Yehuda et al. [ 53 ], enhanced negative feedback action of cortisol does not characterize low cortisol levels observed in PTSD patients within our model. We plan to provide further discussion on this matter in our future work.

Minor issues: The observation that transition to diseased state depends on stress timing is interesting but it can be addressed whether or not this is a biologically relevant phenomenon.

Authors’ response: In a previous study [ 21 ], it was shown that changes in corticosterone levels induced by acute auditory stressor in rats were dependent on the timing of stress onset relative to the phase of underlying corticosterone oscillations. Based on this observation, we believe that the timing of stress onset could be relevant in the transitions to diseased states. Please refer to the Additional file 1 for details on how the experimental observation can be explained within our model.

To address the timing issue in a more detailed manner, we need a better description of the synaptic input function of the PVN, I(t ), that models how and when stress response is initiated and terminated. We have shown in the manuscript that timing may be crucial in inducing transitions, but only when the stress duration is near its transition thresholds. A more realistic form of I(t ) will thus allow us to understand under what circumstances stress duration may be near such transition thresholds. We are are currently investigating the endocannabinoid system that is known to regulate the initiation and termination of stress response, which is subject to fast nongenomic actions of cortisol.