 Research
 Open Access
 Published:
Onset, timing, and exposure therapy of stress disorders: mechanistic insight from a mathematical model of oscillating neuroendocrine dynamics
Biology Directvolume 11, Article number: 13 (2016)
Abstract
Background
The hypothalamicpituitaryadrenal (HPA) axis is a neuroendocrine system that regulates numerous physiological processes. Disruptions in the activity of the HPA axis are correlated with stressrelated diseases such as posttraumatic 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 corticotropinreleasing 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 selfupregulation of CRH neuron activity, our model distinguishes the relatively slow process of cortisolmediated CRH biosynthesis from rapid transsynaptic effects that regulate the CRH secretion process. We show that the slow component of the negative feedback allows external stressinduced reversible transitions between “normal” and “diseased” states in novel intensity, duration, and timingdependent ways.
Conclusion
Our twostep 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 hypothalamicpituitaryadrenal (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 corticotropinreleasing 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].
Dysregulation in the HPA axis is known to correlate with a number of stressrelated disorders. Increased cortisol (hypercortisolism) is associated with major depressive disorder (MDD) [2, 3], while decreased cortisol (hypocortisolism) is a feature of posttraumatic 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 lifethreatening 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 combatrelated 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 timeevolution of the key regulating hormones of the HPA axis: CRH, ACTH, and cortisol. These models [13, 14, 16] incorporate positive selfregulation 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 structurallevel 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 preexisting 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 gonadotropinreleasing 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 stressinduced 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 selfupregulation [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 twostate dynamical structure that arises in the classic FitzhughNagumo 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 phasespace [25]. Our main result is that stressdriven transitions between normal and diseased states can arise when a twostage 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):
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 quasiequilibrium so concentration of the ligandreceptor 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 cortisolGR 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 selfupregulation effect of the cortisolGR 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 Hilltype 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 selfupregulation 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 physiologicallymotivated 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 cortisolmediated CRH biosynthesis from the rapid transsynaptic effects that regulate CRH secretion. This allows us to decompose the dynamics into slow and fast components. Finally, selfupregulation 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 pituitaryadrenal (PA) subsystem. Since oscillations can be induced by delays, we assume, as in Walker et al. [13], a time delay T _{d} in the ACTHmediated activation of cortisol production downstream of the hypothalamus. Equation 3 is thus modified to
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 ACTHmediated 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):
Here, T _{ C } is a characteristic time constant and C _{ ∞ }(O) is the cortisoldependent 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 selfupregulation 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
where I _{base} and I _{ext}(T) represent the basal and stressdependent 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 cortisoldependent 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 membranebound Gproteincoupled receptor CRHR1 whose downstream signaling pathways operate on timescales from milliseconds to seconds [41, 42]. Thus, selfupregulation 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
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 delaydifferential 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 timescale negative feedback by cortisol on CRH synthesis (through the production term C _{ ∞ }(O) in Eq. 6).

A fastacting 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 cortisolGR complex and use Hilltype 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:
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
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
The form of c _{ ∞ }(o) is based on the abovementioned 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 Hilltype function, with Hillexponent 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 halfmaximum.
Fastslow 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
where x=(c,a,o,r) is the vector of fast dynamical variables, and F(c _{s},x) denotes the righthandsides 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
where \(q(c_{\mathrm {s}}(\tau),I) \equiv q_{0}I h(c_{\mathrm {s}}(\tau))=q_{0}I (1e^{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 nonoscillating quasiequilibrium value defined by dc/dt=q g _{ c }(c)−q _{2} c=0, it can be viewed as a parametric term in Eq. 16 of the pituitaryadrenal (PA) subsystem.
Due to the nonlinearity of g _{ c }(c), the equilibrium value c(q) satisfying q g _{ c }(c)=q _{2} c may be multivalued 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 _{1}, bistability can emerge through a saddlenode 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 _{2} c.
For equilibrium values of c lying within a certain range, the PAsubsystem 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 PAsubsystem depicted in Fig. 3 indicate the range of c values that admit limit cycle behavior for (a,o,r), while the fast cnullcline 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.
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 fivedimensional model, we only need to consider the 2D projection onto the fast c and slow c _{s} variable. A summary of the timeseparated dynamics of the variables in our model is given in Fig. 4.
To analyze the evolution of the slow variable c _{s}(τ), we write our equations in terms of τ=ε t:
In the ε→0 limit, the “outer solution” F(c _{s},x)≈0 simply constrains the system to be on the fast cnullcline defined by q g _{ c }(c)=q _{2} c. The slow evolution of c _{s}(τ) along the fast cnullcline 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 periodaveraged value
Since periodaveraged values of o ^{∗} increases with c, 〈c _{ ∞ }(c)〉 is a decreasing function of c under physiological parameter regimes. This periodaveraging 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
with 〈c _{ ∞ }(c)〉 evaluated in Eq. 25. By selfconsistently solving Eqs. 26 and 27, we can estimate trajectories of the full model when they are near the cnullcline 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 cnullcline 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 cnullcline 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(1e^{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 “qnullclines.”
We assume a fixed basal stress input I=1 and plot the qnullclines 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)=(1e^{k\langle c_{\infty }(c)\rangle })\), both the position and the steepness of the qnullcline in (q,c)space depend strongly on k. Figure 5 b shows a fast cnullcline and a slow qnullcline (transformed c _{s}nullcline) intersecting at both stable branches of the fast cnullcline. Here, the flow field indicates that the 2D projected trajectory is governed by fast flow over most of the (q,c)space.
How the fast and slow nullclines intersect controls the longterm 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 _{0},…,p _{6},b,k,n,μ _{c},t _{d}).
Other parameters such as q _{0}, q _{1}, and μ _{c} appear directly in the fast equation for c and thus most strongly control the fast cnullcline. 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 cnullcline while the slow variable q=q _{0} I h(c _{s}) remains fairly constant. Independent of initial configurations that are not near the cnullcline in (q,c)space, trajectories quickly jump to one of the stable branches of the cnullcline with little motion towards the qnullcline, as indicated by ξ _{f} in Fig. 6 a.
Once near the cnullcline, 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 cnullcline towards the qnullcline (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 cnullcline 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 cnullcline as shown in Fig. 6 b. Trajectories that start within the basin of attraction of the lower stable branch of the cnullcline (“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 longterm 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 dualnullcline 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 cnullcline 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 stressdriven transitions mediated by the fastslow 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 qnullcline 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 cortisolGR complex, which can be described by decreasing p _{2} in our model. Such change in parameter would lead to a leftward shift of the qnullcline 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 timeautonomous 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 _{0},q _{1}, and k, we will study how our model depends on k while fixing μ _{c},q _{0}, and q _{1}. Three possible nullcline configurations arise according to the values of μ _{c},q _{0}, and q _{1} 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.
Given these considerations, we henceforth chose μ _{c}=0.6, q _{1}=0.04, and q _{0}=77.8 for the rest of our analysis. This choice of parameters is motivated in the Additional file 1 and corresponds to a socalled “Type I” nullcline structure. In this case, three possibilities arise: one intersection on the lower stable branch of the cnullcline 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 cnullcline 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 stressinduced 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.
Upon turning on stress, the lumped parameter q and the slow nullcline shift to the right by 10 % since q=q _{0}(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 cnullcline; afterwards, it moves very slowly along the cnullcline towards the shifted qnullcline. After 120 min, the system arrives at the “ ×” on the cnullcline (Fig. 7 a). Once the stress is shut off the qnullcline 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 upperbranch 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 highc state.
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 highc stable branch of the cnullcline. 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 normaltodiseased 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 normaltodiseased 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 longterm 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 cnullcline, rendering the HPA system resistant to stressinduced transitions.
The response to chronic stress initially follows the same pattern as described above for the twostablestate case, as shown in Fig. 8 c. However, the system will continue to evolve along the lower branch towards the qnullcline, 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 fixedduration 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 lowcortisol 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 lowq separatrix into the diseased state.
As discussed earlier, an increase in periodaveraged cortisol level during stress drives a normaltodiseased state transition. We see that the periodaveraged 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 periodaveraging 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, computersimulated “virtual reality exposure.” Within our model we can describe ET as external stress (I _{ext}>0) applied to a system in the stable lowc diseased state. The resulting horizontal shift in q causes the system to move rightward across the separatrix and suggests a transition to the highc 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 cnullcline upon termination leading to the normal, highcortisol state (Fig. 10 b). Since the initial motion is governed by fast flow, the minimum stress duration needed to incite the diseasedtonormal 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 lowcortisol 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.
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 diseasedtonormal 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 highcortisol 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 hypocortisolismassociated 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 selfupregulation of CRH release, a slow negative feedback effect of cortisol on CRH synthesis, and a delay in ACTHactivated 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 lowCRH equilibrium point arises; an individual with such k is trapped in the lowcortisol “diseased” state. For large k, only the highCRH normal state arises, rendering the individual resistant to acquiring the longterm, lowcortisol sideeffect 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 agedependent increase in k, which, in our model, implies that aging makes the organism more resistant to stressinduced 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 hypothalamicpituitarygonadal (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 qnullcline 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 counterintuitive phenomenon of lower cortisol levels after stressinduced activation of cortisol production. Solutions to our model demonstrate that the negativefeedback 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 negativefeedback of cortisol [53, 54].
We also find that external stress can induce the “reverse” transition from a diseased lowcortisol state to the normal highcortisol state. Our results imply that reexposure 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 negativefeedback 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, normaltodiseased state transitions are easier to induce when stress is initiated near the nadir of cortisol oscillations. Reverse diseasedtonormal 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 reexperiencing stressful situations directly or through imagination, and have been consistently proven effective as firstline 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 higherbrain 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, experiencedependent 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 _{0},q _{1} , 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 ACTHmediated 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 stressinduced transitions into a lowcortisol 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 posttraumatic 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 qnullcline in our model shifts toward the normal, high cortisol state (on the upper branch of the cnullcline), 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.
Abbreviations
 ACTH:

adrenocorticotropic hormone
 BEP:

brief eclectic psychotherapy
 CFS:

chronic fatigue syndrome
 CORT:

glucocorticoids
 CRH:

corticotropinreleasing hormone
 CRHR1:

Corticotropinreleasing hormone receptor 1
 DDE:

delay differential equation
 ET:

exposure therapy
 GH:

growth hormone
 GnRH:

gonadotropinreleasing hormone
 GR:

glucocorticoid receptors
 MDD:

major depressive disorder
 hnRNA:

heterogeneous nuclear RNA
 HPA:

hypothalamicpituitaryadrenal
 HPG:

hypothalamicpituitarygonadal
 mRNA:

messenger ribonucleic acid
 ODE:

ordinary differential equation
 PA:

pituitaryadrenal
 POA:

preoptic area
 PTSD:

posttraumatic stress disorder
 PVN:

paraventricular nucleus
 TBI:

traumatic brain injury
References
 1
Denver R. Structural and functional evolution of vertebrate neuroendocrine stress systems. Ann N Y Acad Sci. 2009; 1163(1):1–16.
 2
Gold P, Chrousos G. Organization of the stress system and its dysregulation in melancholic and atypical depression: high vs low CRH/NE states. Mol Psychiatry. 2002; 7(3):254–75.
 3
Juruena M, Cleare A, Pariante C. The hypothalamic pituitary adrenal axis, glucocorticoid receptor function and relevance to depression. Rev Bras Psiquiatr. 2004; 26(3):189–201.
 4
Rohleder N, Joksimovic L, Wolf J, Kirschbaum C. Hypocortisolism and increased glucocorticoid sensitivity of proinflammatory cytokine production in bosnian war refugees with posttraumatic stress disorder. Biol Psychiatry. 2004; 55(7):745–51.
 5
Giorgio AD, Hudson M, Jerjes W, Cleare A. 24hour pituitary and adrenal hormone profiles in chronic fatigue syndrome. Psychosom Med. 2005; 67(3):433–40.
 6
Jerjesnd W, Peters T, Taylor N, Wood P, Wessely S, Cleare A. Diurnal excretion of urinary cortisol, cortisone, and cortisol metabolites in chronic fatigue syndrome. J Psychosom Res. 2006; 60(2):145–53.
 7
Crofford L, Young E, Cary NEK, Korszun A, Brucksch C, McClure L, Brown M, Demitrack M. Basal circadian and pulsatile ACTH and cortisol secretion in patients with fibromyalgia and/or chronic fatigue syndrome. Brain Behav Immun. 2004; 18(4):314–25.
 8
Yehuda R, Teicher M, Levengood R, Trestman R, Siever L. Circadian regulation of basal cortisol levels in posttraumatic stress disorder. Ann N Y Acad Sci. 1994; 746(1):378–80.
 9
Vinther F, Andersen M, Ottesen JT. The minimal model of the hypothalamicpituitaryadrenal axis. J Math Biol. 2011; 63:663–90.
 10
Jelic S, Cupic Z, KolarAnic L. Mathematical modeling of the hypothalmicpituitaryadrenal system activity. Math Biosci. 2005; 197:173–87.
 11
Kyrylov V, Severyanova L, Vieira A. Modeling robust oscillatory behavior of the hypothalamicpituitaryadrenal axis. IEEE Trans Biomed Eng. 2005; 52(12):1977–83.
 12
Savić D, Knežević G, Opačić G. A mathematical model of stress reaction: Individual differences in threshold and duration. Psychobiology. 2000; 28(4):581–92.
 13
Walker JJ, Terry JR, Lightman SL. Origin of ultradian pulsatility in the hypothalamic–pituitary–adrenal axis. Proc R Soc Lond B Biol Sci. 2010; 277(1688):1627–33.
 14
Rankin J, Walker J, Windle R, Lightman S, Terry J. Characterizing dynamic interactions between ultradian glucocorticoid rhythmicity and acute stress using the phase response curve. PloS ONE. 2012; 7(2):30978.
 15
Bairagi N, Chatterjee S, Chattopadhyay J. Variability in the secretion of corticotropinreleasing hormone, adrenocorticotropic hormone and cortisol and understandability of the hypothalamicpituitaryadrenal axis dynamics — a mathematical study based on clinical evidence. Math Med Biol. 2008; 25:37–63.
 16
Sriram K, RodriguezFernandez M, Doyle III FJ. Modeling cortisol dynamics in the neuroendocrine axis distinguishes normal, depression, and posttraumatic stress disorder (PTSD) in humans. PLoS Comput Biol. 2012; 8:1002379.
 17
Gupta S, Aslakson E, Gurbaxani BM, Vernon SD. Inclusion of the glucocorticoid receptor in a hypothalamic pituitary adrenal axis model reveals bistability. Theor Biol Med Model. 2007; 4:8.
 18
Windle R, Wood S, Lightman S, Ingram C. The pulsatile characteristics of hypothalamopituitaryadrenal activity in female Lewis and Fischer 344 rats and its relationship to differential stress responses. Endocrinology. 1998; 139(10):4044–52.
 19
Chrousos G. Editorial: ultradian, circadian, and stressrelated hypothalamicpituitaryadrenal axis activity — a dynamic digitaltoanalog modulation. Endocrinology. 1998; 139(2):437–40.
 20
ConwayCampbell B, Sarabdjitsingh R, McKenna M, Pooley J, Kershaw Y, Meijer O, Kloet ED, Lightman S. Glucocorticoid ultradian rhythmicity directs cyclical gene pulsing of the clock gene period 1 in rat hippocampus. J Neuroendocrinol. 2010; 22(10):1093–100.
 21
Windle R, Wood S, Shanks N, Lightman S, Ingram C. Ultradian rhythm of basal corticosterone release in the female rat: Dynamic interaction with the response to acute stress. Endocrinology. 1998; 139(2):443–50.
 22
Watts A. Glucocorticoid regulation of peptide genes in neuroendocrine CRH neurons: a complexity beyond negative feedback. Front Neuroendocrinol. 2005; 26(3):109–30.
 23
Ono N, Castro JD, McCann S. Ultrashortloop positive feedback of corticotropin (ACTH)releasing factor to enhance ACTH release in stress. Proc Natl Acad Sci. 1985; 82(10):3528–31.
 24
FitzHugh R. Mathematical models of threshold phenomena in the nerve membrane. Bull Math Biophys. 1955; 17(4):257–78.
 25
Silva FLD, Blanes W, Kalitzin S, Parra J, Suffczynski P, Velis D. Epilepsies as dynamical diseases of brain systems: basic models of the transition between normal and epileptic activity. Epilepsia. 2003; 44(s12):72–83.
 26
BenZvi A, Vernon SD, Broderick G. Modelbased therapeutic correction of hypothalamicpituitaryadrenal axis dysfunction. PLoS Comput Biol. 2009; 5(1):1000273.
 27
Tsai SY, CarlstedtDuke J, Weigel NL, Dahlman K, Gustafsson J, Tsai MJ, O’Malley BW. Molecular interactions of steroid hormone receptor with its enhancer element: evidence for receptor dimer formation. Cell. 1988; 55(2):361–9.
 28
Andersen M, Vinther F, Ottesen J. Mathematical modeling of the hypothalamic–pituitary–adrenal gland (HPA) axis, including hippocampal mechanisms. Math Biosci. 2013; 246(1):122–38.
 29
Papaikonomou E. Rat adrenocortical dynamics. J Physiol. 1977; 265(1):119–31.
 30
Engler D, Pham T, Liu J, Fullerton M, Clarke I, Funder J. Studies of the regulation of the hypothalamicpituitaryadrenal axis in sheep with hypothalamicpituitary disconnection. II, evidence for in vivo ultradian hypersecretion of proopiomelanocortin peptides by the isolated anterior and intermediate pituitary. Endocrinology. 1990; 127(4):1956–66.
 31
Weiser M, Osterlund C, Spencer R. Inhibitory effects of corticosterone in the hypothalamic paraventricular nucleus (pvn) on stressinduced adrenocorticotrophic hormone secretion and gene expression in the pvn and anterior pituitary. J Neuroendocrinol. 2011; 23(12):1231–40.
 32
Ma X, Aguilera G. Differential regulation of corticotropinreleasing hormone and vasopressin transcription by glucocorticoids. Endocrinology. 1999; 140(12):5642–50.
 33
Watts A, SanchezWatts G. Regionspecific regulation of neuropeptide mRNAs in rat limbic forebrain neurones by aldosterone and corticosterone. J Physiol. 1995; 484(3):721–36.
 34
Tasker J, Di S, MalcherLopes R. Rapid glucocorticoid signaling via membraneassociated receptors. Endocrinology. 2006; 147(12):5549–56.
 35
Kasai M, Yamashita H. Inhibition by cortisol of neurons in the paraventricular nucleus of the hypothalamus in adrenalectomized rats; an in vitro study. Neurosci Lett. 1988; 91(1):59–64.
 36
Kasai M, Yamashita H. Cortisol suppresses noradrenalineinduced excitatory responses of neurons in the paraventricular nucleus; an in vitro study. Neurosci Lett. 1988; 91(1):65–70.
 37
Jones M, Hillhouse E, Burden J. Dynamics and mechanics of corticosteroid feedback at the hypothalamus and anterior pituitary gland. J Endocrinol. 1977; 73(3):405–17.
 38
Ginsberg A, Campeau S, Day H, Spencer R. Acute glucocorticoid pretreatment suppresses stressinduced hypothalamicpituitaryadrenal axis hormone secretion and expression of corticotropinreleasing hormone hnRNA but does not affect cfos mRNA or fos protein expression in the paraventricular nucleus of the hypothalamus. J Neuroendocrinol. 2003; 15(11):1075–83.
 39
Chen Y, Hua S, Wang C, Wu L, Gu Q, Xing B. An electrophysiological study on the membrane receptormediated action of glucocorticoids in mammalian neurons. Neuroendocrinology. 1991; 53(Suppl. 1):25–30.
 40
Imaki T, XiaoQuan W, Shibasaki T, Yamada K, Harada S, Chikada N, Naruse M, Demura H. Stressinduced activation of neuronal activity and corticotropinreleasing factor gene expression in the paraventricular nucleus is modulated by glucocorticoids in rats. J Clin Investig. 1995; 96(1):231.
 41
Papadimitriou A, Priftis K. Regulation of the hypothalamicpituitaryadrenal axis. Neuroimmunomodulation. 2009; 16(5):265.
 42
Makino S, Hashimoto K, Gold P. Multiple feedback mechanisms activating corticotropinreleasing hormone system in the brain during stress. Pharmacol Biochem Behav. 2002; 73(1):147–58.
 43
Herman JP, Figueiredo H, Mueller NK, UlrichLai Y, Ostrander M, Choi D, Cullinan W. Central mechanisms of stress integration: hierarchical circuitry controlling hypothalamo–pituitary–adrenocortical responsiveness. Front Neuroendocrinol. 2003; 24(3):151–80.
 44
McEwen BS, Stellar E. Stress and the individual: mechanisms leading to disease. Arch Intern Med. 1993; 153:2093–101.
 45
McEwen BS. Stress, adaptation, and disease: Allostasis and allostatic load. Ann N Y Acad Sci. 1998; 840(1):33–44.
 46
Dince SM, Rome RD, McEwen BS, Tang AC. Enhancing offspring hypothalamicpituitaryadrenal (hpa) regulation via systematic novelty exposure: the influence of maternal HPA function. Front Behav Neurosci. 2014; 8:204. doi:10.3389/fnbeh.2014.00204.
 47
Hauger RL, Thrivikraman KV, Plotsky PM. Agerelated alterations of hypothalamicpituitaryadrenal axis function in male Fischer 344 rats. Endocrinology. 1994; 134(3):1528–36.
 48
Averill P, Beck J. Posttraumatic stress disorder in older adults: a conceptual review. J Anxiety Disord. 2000; 14(2):133–56.
 49
Regier D, Boyd J, Burke J, Rae D, Myers J, Kramer M, Robins L, George L, Karno M, Locke B. Onemonth prevalence of mental disorders in the united states: based on five epidemiologic catchment area sites. Arch Gen Psychiatr. 1988; 45(11):977–86.
 50
Simerly RB, Swanson LW, Chang C, Muramatsu M. Distribution of androgen and estrogen receptor mrnacontaining cells in the rat brain: An in situ hybridization study. J Comp Neurol. 1990; 294(1):76–95.
 51
Gréco B, Allegretto E, Tetel M, Blaustein J. Coexpression of ER β with ER α and progestin receptor proteins in the female rat forebrain: effects of estradiol treatment. Endocrinology. 2001; 142(12):5172–81.
 52
Feldman S, Conforti N, Saphier D. The preoptic area and bed nucleus of the stria terminalis are involved in the effects of the amygdala on adrenocortical secretion. Neuroscience. 1990; 37(3):775–9.
 53
Yehuda R, Teicher M, Levengood R, Trestman R, Levengood R, Siever L. Cortisol regulation in posttraumatic stress disorder and major depression: a chronobiological analysis. Biol Psychiatry. 1996; 40(2):79–88.
 54
Yehuda R, LeDoux J. Response variation following trauma: a translational neuroscience approach to understanding PTSD. Neuron. 2007; 56:19–32.
 55
Olff M, de Vries G, Güzelcan Y, Assies J, Gersons B. Changes in cortisol and DHEA plasma levels after psychotherapy for PTSD. Psychoneuroendocrinology. 2007; 32(6):619–26.
 56
Foa E, Keane T, Friedman M, Cohen J. Effective treatments for PTSD: practice guidelines from the international society for traumatic stress studies. New York: Guilford Press: 2008.
 57
Rauch S, Eftekhari A, Ruzek J. Review of exposure therapy: a gold standard for PTSD treatment. J Rehabil Res Dev. 2012; 49(5):679–88.
 58
Trouche S, Sasaki J, Tu T, Reijmers L. Fear extinction causes targetspecific remodeling of perisomatic inhibitory synapses. Neuron. 2013; 80(4):1054–65.
 59
Lemieux A, Coe C. Abuserelated posttraumatic stress disorder: evidence for chronic neuroendocrine activation in women. Psychosom Med. 1995; 57(2):105–15.
 60
Baum A. Implications of psychological research on stress and technological accidents. Am Psychol. 1993; 48(6):665.
 61
Rosenmund C, Stevens C. Definition of the readily releasable pool of vesicles at hippocampal synapses. Neuron. 1996; 16(6):1197–207.
Acknowledgments
This work was supported by the Army Research Office via grant W911NF1410472 and the NSF through grant BCS1348123. The authors also thank professors T. Minor and M. Wechselberger for insightful discussions.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
LK performed the mathematical analysis and computational simulations. MRD and TC formulated and designed the study. All authors contributed significantly to the analysis of the mathematical simulations and the writing of the manuscript. All authors read and approved the final manuscript.
Additional file
Additional file 1
The file provides appendices describing the mathematical details of our analysis. It includes discussions on nondimensionalization and parameter estimations used in our model. Analyses of the nullcline structure and the parameter space are also provided and the effect of stress onset timing on initial stress response is described in detail. Finally, a possible form of cortisol dependent I _{ext} is illustrated with its implications on model behavior. (PDF 2099 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 HPAaxis
 PTSD
 Stress disorders
 Dynamical system