Onset, timing, and exposure therapy of stress disorders: mechanistic insight from a mathematical model of oscillating neuroendocrine dynamics
 Lae U. Kim^{1},
 Maria R. D’Orsogna^{2} and
 Tom Chou^{3}Email author
DOI: 10.1186/s1306201601176
© Kim et al. 2016
Received: 24 December 2015
Accepted: 17 March 2016
Published: 25 March 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.
Keywords
HPAaxis PTSD Stress disorders Dynamical systemOpen 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
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
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
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.
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.
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.
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

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).
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
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
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 τ.
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.
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.”
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}).
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.
Dimensionless parameter values of our full model
Parameter  Value  Source and Ref.  Description 

n  5  Assumed  Hill coefficient in upregulation function g _{ c }(c) 
\(\bar {c}_{\infty }\)  0.2  Estimated from [22]  Baseline stored CRH level 
b  0.6  Estimated from [22]  Relates cortisol to stored CRH level 
k  Undetermined  ·  Relates stored CRH to CRH release rate 
μ _{c}  Undetermined  ·  Basal CRH release rate 
q _{0}  Undetermined  ·  Maximum CRH release rate 
\(q_{1}^{1}\)  Undetermined  ·  Circulating CRH for halfmaximum selfupregulation 
q _{2}  1.8  Estimated from [21]  Ratio of CRH and cortisol decay rates 
\(p_{2}^{1}\)  0.067  \(p_{2}^{1}\) [13]  (o r)complex level for halfmaximum feedback 
p _{3}  7.2  p _{3} [13]  Ratio of ACTH and cortisol decay rates 
p _{4}  0.05  p _{4} [13]  (o r)complex level for halfmaximum upregulation 
p _{5}  0.11  p _{5} [13]  Basal GR production rate by pituitary 
p _{6}  2.9  p _{6} [13]  Ratio of GR and cortisol decay rates 
t _{c}  69.3  Assumed  CRH biosynthesis timescale 
t _{d}  1.44  “ τ” [13]  Delay in ACTHactivated cortisol release 
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
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
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
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.
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
Declarations
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.
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.
Authors’ Affiliations
References
 Denver R. Structural and functional evolution of vertebrate neuroendocrine stress systems. Ann N Y Acad Sci. 2009; 1163(1):1–16.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticleGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 Vinther F, Andersen M, Ottesen JT. The minimal model of the hypothalamicpituitaryadrenal axis. J Math Biol. 2011; 63:663–90.View ArticlePubMedGoogle Scholar
 Jelic S, Cupic Z, KolarAnic L. Mathematical modeling of the hypothalmicpituitaryadrenal system activity. Math Biosci. 2005; 197:173–87.View ArticlePubMedGoogle Scholar
 Kyrylov V, Severyanova L, Vieira A. Modeling robust oscillatory behavior of the hypothalamicpituitaryadrenal axis. IEEE Trans Biomed Eng. 2005; 52(12):1977–83.View ArticlePubMedGoogle Scholar
 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.Google Scholar
 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.View ArticleGoogle Scholar
 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.View ArticleGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticleGoogle Scholar
 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.View ArticlePubMedPubMed CentralGoogle Scholar
 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.PubMedGoogle Scholar
 Chrousos G. Editorial: ultradian, circadian, and stressrelated hypothalamicpituitaryadrenal axis activity — a dynamic digitaltoanalog modulation. Endocrinology. 1998; 139(2):437–40.PubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.PubMedGoogle Scholar
 Watts A. Glucocorticoid regulation of peptide genes in neuroendocrine CRH neurons: a complexity beyond negative feedback. Front Neuroendocrinol. 2005; 26(3):109–30.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedPubMed CentralGoogle Scholar
 FitzHugh R. Mathematical models of threshold phenomena in the nerve membrane. Bull Math Biophys. 1955; 17(4):257–78.View ArticleGoogle Scholar
 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.View ArticleGoogle Scholar
 BenZvi A, Vernon SD, Broderick G. Modelbased therapeutic correction of hypothalamicpituitaryadrenal axis dysfunction. PLoS Comput Biol. 2009; 5(1):1000273.View ArticleGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 Papaikonomou E. Rat adrenocortical dynamics. J Physiol. 1977; 265(1):119–31.View ArticlePubMedPubMed CentralGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedPubMed CentralGoogle Scholar
 Ma X, Aguilera G. Differential regulation of corticotropinreleasing hormone and vasopressin transcription by glucocorticoids. Endocrinology. 1999; 140(12):5642–50.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedPubMed CentralGoogle Scholar
 Tasker J, Di S, MalcherLopes R. Rapid glucocorticoid signaling via membraneassociated receptors. Endocrinology. 2006; 147(12):5549–56.View ArticlePubMedPubMed CentralGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedPubMed CentralGoogle Scholar
 Papadimitriou A, Priftis K. Regulation of the hypothalamicpituitaryadrenal axis. Neuroimmunomodulation. 2009; 16(5):265.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 McEwen BS, Stellar E. Stress and the individual: mechanisms leading to disease. Arch Intern Med. 1993; 153:2093–101.View ArticlePubMedGoogle Scholar
 McEwen BS. Stress, adaptation, and disease: Allostasis and allostatic load. Ann N Y Acad Sci. 1998; 840(1):33–44.View ArticlePubMedGoogle Scholar
 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.Google Scholar
 Hauger RL, Thrivikraman KV, Plotsky PM. Agerelated alterations of hypothalamicpituitaryadrenal axis function in male Fischer 344 rats. Endocrinology. 1994; 134(3):1528–36.PubMedGoogle Scholar
 Averill P, Beck J. Posttraumatic stress disorder in older adults: a conceptual review. J Anxiety Disord. 2000; 14(2):133–56.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.PubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 Yehuda R, LeDoux J. Response variation following trauma: a translational neuroscience approach to understanding PTSD. Neuron. 2007; 56:19–32.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.
 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.View ArticlePubMedGoogle Scholar
 Trouche S, Sasaki J, Tu T, Reijmers L. Fear extinction causes targetspecific remodeling of perisomatic inhibitory synapses. Neuron. 2013; 80(4):1054–65.View ArticlePubMedGoogle Scholar
 Lemieux A, Coe C. Abuserelated posttraumatic stress disorder: evidence for chronic neuroendocrine activation in women. Psychosom Med. 1995; 57(2):105–15.View ArticlePubMedGoogle Scholar
 Baum A. Implications of psychological research on stress and technological accidents. Am Psychol. 1993; 48(6):665.View ArticlePubMedGoogle Scholar
 Rosenmund C, Stevens C. Definition of the readily releasable pool of vesicles at hippocampal synapses. Neuron. 1996; 16(6):1197–207.View ArticlePubMedGoogle Scholar