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

Background The hypothalamic-pituitary-adrenal (HPA) axis is a neuroendocrine system that regulates numerous physiological processes. Disruptions in the activity of the HPA axis are correlated with stress-related diseases such as post-traumatic stress disorder (PTSD) and major depressive disorder. In this paper, we characterize “normal” and “diseased” states of the HPA axis as basins of attraction of a dynamical system describing the inhibition of peptide hormones such as corticotropin-releasing hormone (CRH) and adrenocorticotropic hormone (ACTH) by circulating glucocorticoids such as cortisol (CORT). Results In addition to including key physiological features such as ultradian oscillations in cortisol levels and self-upregulation of CRH neuron activity, our model distinguishes the relatively slow process of cortisol-mediated CRH biosynthesis from rapid trans-synaptic effects that regulate the CRH secretion process. We show that the slow component of the negative feedback allows external stress-induced reversible transitions between “normal” and “diseased” states in novel intensity-, duration-, and timing-dependent ways. Conclusion Our two-step negative feedback model suggests a mechanism whereby exposure therapy of stress disorders such as PTSD may act to normalize downstream dysregulation of the HPA axis. Our analysis provides a causative rationale for improving treatments and guiding the design of new protocols. Reviewers This article was reviewed by Dr. Daniel Coombs, Dr. Yang Kuang, and Dr. Ha Youn Lee. Electronic supplementary material The online version of this article (doi:10.1186/s13062-016-0117-6) contains supplementary material, which is available to authorized users.


Additional File 1
Below, we provide appendices describing the mathematical details of our analysis.

Nondimensionalization
Our equations are nondimensionalized in a manner similar to that used by Walker et al.: Here, c s , c, a, r, o are the dimensionless versions of the original concentrations C s , C, A, R, O, respectively. C s is normalized byC s , which denotes the typical maximum amount of releasable CRH in the physiological range. Upon using these variables, the dimensionless forms of Eqs. 9-13 are expressed in Eqs. 14-18. The parameters q i , p i are dimensionless combinations conveniently defined to be analogous to those used by Walker et al.: Using these scalings, we arrive at the dimensionless Eqs. 14-19.

Parameter estimates
Many of the numerous physiological parameters in our model can be estimated or constructed from previous studies on the HPA axis. For example, as shown in Fig. A1, the parameters forming the function c ∞ (o) are derived from fitting to data on adrenalectomized male rats. From the fitting, we estimate the baseline level  Figure A1: Fitting c∞(o) to rat data. Nondimensionalized data taken from Watts and fitted using the form for c∞(o) given in Eq. 19. Cortisol levels were arbitrarily rescaled according to 125ng/ml → 3. c ∞ 0.2, and the decay rate b 0.56. Furthermore, the dimensionless parameters p 2 , . . . , p 5 and t d will be fixed to those used in Walker et al.: p 2 = 15, p 3 = 7.2, p 4 = 0.05, p 5 = 0.11, p 6 = 2.9 and t d = 1.44 (T d = 15 min). Although it is not possible to determine all of the remaining parameters from data, we will use reasonable estimates. The half-life of cortisol was estimated to be about 7.2min while the halflife of CRH has been estimated to be about 4min. Therefore, q 2 = d C d −1 O ≈ 1.8. Of the remaining parameters (n, µ c , q 0 , q 1 , k), the dependence on n will turn out to be quantitative so we henceforth set n = 5. These estimated parameters are listed in Table A1.
Even though one expects the values of these effective parameters to be highly variable, we fix them in order to concretely investigate the mathematical structure and qualitative predictions of our model. The parameters µ c , q 0 , q 1 , and k remain undetermined; however, it is instructive to treat k as a control parameter and explore the nullcline structure in µ c , q 0 , q 1 parameter space.

Parameter space and nullcline structure
To determine how the q-nullcline crosses the c-nullcline, we substitute c s by its equilibrium period-averaged value c ∞ (c) . If we assume a basal input level I = 1, the values of k that will position the basal q-nullcline to just pass through the left and right bifurcation points (q L , c L ) and (q R , c R ) can be found by solving q L,R = q 0 (1 − e −k c∞(cL,R) ): All possible ways in which the nullclines can cross each other as k is varied are illustrated in Fig. A2.
Type I Type II Figure A2: The possible number of equilibria of the reduced (q, c) system. (A) A Type 0 scenario in which kR < kL permits only one nullcline intersection, either one the lower stable branch, the unstable branch, or the upper stable branch. (B) In this Type I parameter regime, the c-nullcline is shaped and positioned such that kL < kR. Therefore, it is possible for the model to exhibit two oscillating stable states provided kL < k < kR. For k < kL (k > kR), the q-nullcline shifts to the left(right) and the intersection with the upper(lower) branch of the c-nullcline disappears, leading to only one stable point. (C) A Type II c-nullcline. For k < kL, there is only one intersection at the lower branch. For all k > kL there are two intersections.
The specific locations of the bifurcation points, as well as k L and k R , are complicated functions of all parameters. However, Eqs. A3 allows us to distinguish three qualitatively different regimes. The first possibility is k L > k R , where there can be at most only one intersection between the slow and fast nullclines. We denote this as a Type 0 scenario (Fig. A2A) characterized by having at most a single stable state towards which the system will always return upon cessation of external stress. In Type 0 situations with intermediate values of k, the intersection will arise in the unstable branch of the c-nullcline. In this case, we expect the system to oscillate between the two stable branches of the c-nullcline. Here, the fast variables a, o, and r will cycle periodically between two oscillating levels.
In order for the two nullclines to intersect three times (twice on stable branches of the c-nullcline), the q-nullcline must "fit" within the bistable region of the cnullcline. As shown in Fig. A2, there are two separate subcases of nullclines that intersect twice. If k L < k R , a value of k L < k < k R would imply that the q−nullcline can intersect both stable branches of the c-nullcline, leading to two stable solutions. We refer to this case as Type I (Fig. A2B).
Another possibility is that the right bifurcation point is beyond the maximum value q = q 0 dictated by the function h( c ∞ (c R ) ). As shown in Fig. A2C, the bistable c-nullclines exhibits only one bifurcation point within the domain of q. The lower branch of the c-nullclines in this set extends across the entire range of physiological values of q, ensuring that the q-nullcline will intersect with the lower branch for any value of k. Therefore, to determine if there are two intersections we only need to check that k L ≤ k is satisfied. In this Type II case, the system is either perpetually in the diseased low cortisol state, or is bistable between the diseased and normal states; the system will always be at least susceptible to low-cortisol disease. Summarizing, -Type 0: Exactly one solution (one nullcline intersection) exists for the reduced subsystem. Here, k R < k L and the intersection may occur on the lower or upper stable branches, or on the unstable branch of the c-nullcline. The system is either permanently diseased, permanently resistant, or oscillates between normal and diseased states.
-Type I: At least one solution exists. A stable diseased solution exists if k < k L , two stable solutions (diseased and normal) arise if k L ≤ k ≤ k R , and fully resistant state arises if k > k R .
-Type II: At least one solution exists. A stable diseased state arises if k < k L while both diseased and normal solutions arise if k > k L . A fully diseaseresistant state cannot arise. With the parameters fixed according to Table A1, we will treat k as a control parameter and exhaustively sweep the three-dimensional parameter space (q 0 , q 1 , µ c ) to determine the regions which lead to each of the nullcline structural types. In addition, we restrict the parameter domain to regions which admit oscillating solutions of the full problem. In other words, parts of both stable branches of the c-nullclines must fall within values of c which support oscillations in the PA-subsystem (Fig. 3). The regions in (µ c , q 0 , q 1 ) space that satisfy these conditions and that yield each of the types of nullcline crossings are indicated in Fig. A3. Based on measurements of self-upregulation of CRH secretion during stress, µ c = 0.6 is chosen to set the baseline level of the Hill function g c (c = 0) ≈ 0.4. q 1 is approximated by setting the inflection point of g c (c) to arise at c ≈ 25, the average value used by Walker et al. Assuming c ≈ 25 is a fixed point of Eq. 15 when I = 1 and c s ≈ c ∞ (25) , q 0 can be estimated as a root of the right-hand-side of Eq. 15. This choice for the remaining parameters puts our nullcline system into the Type I category that can exhibit one or two stable states with oscillating (a, o, r) subsystems. We restricted the analysis of our model to Type I systems.

Minimum duration and magnitude of stress
We plot the minimum duration required for normal-to-diseased transition against stress magnitude (Fig. A4). Higher magnitude of I ext generally requires a shorter duration of stress, as expected. Note that the minimum duration is also dependent on the phase of intrinsic oscillations of the system at stress onset.

Timing of stress onset and transient response
Here, we show how the dynamics of the system changes after the onset and cessation of stress. In previous studies, changes in corticosterone levels in rats weremeasured in response to stress induced by noise applied at different phases of the animals oscillating cortisol cycle. It was observed that the timing of the stress onset relative to the ultradian phase was crucial in determining the magnitude of corticosterone response. Increases in corticosterone levels were markedly higher when noise was initiated during the rising phase than when initiated during the falling phase.
We can frame these experimental observations mechanistically within our theory. Following the experimental protocol, we simulate the stress response using a brief stressor with a duration of 30min. As shown in Fig. A5A, an external stress that is applied mostly over the falling phase of the cortisol oscillation results in a higher subsequent nadir in o(t) than one that is applied predominantly during a rising phase. However, as shown in Fig. A5B, stress applied mainly during the rising phase leads to a higher subsequent peak level. This observation is consistent with the results of the experiment on rats and can be explained by the dynamics inherent in our model.
The immediate increase in q = q 0 Ih(c s ) associated with the increase in I leads to a rapid increase in c, as shown in Fig. 7. This higher level of circulating CRH shifts the stable limit cycle of the PA subsystem to a new one with higher minimum and maximum values of ACTH and cortisol (as shown in Fig. 3). This new limit cycle is shown by the blue curve in Figs. A5C,D. Under external stress, a trajectory of the system quickly deviates and approaches the new limit cycle, but quickly returns to the original limit cycle after cessation of stress. Thus, depending on the position of the trajectory relative to that of the new stressed limit cycle, the initial deviation may try to reach the new limit cycle in the falling or rising cortisol phases as shown in Figs. A5C,D. Moreover, if the duration of the stress is shorter than the period of the inherent oscillation, the trajectory will return to its original limit cycle before completing a full cycle of the new limit cycle. These properties of the limit cycle dynamics explain the difference in the level of subsequent peak following the stress onset depending on the timing of the stress onset.
Cortisol-dependent I ext Since it has been shown that synaptic input into PVN cells is modulated by cortisol for certain stressor types, we briefly discuss how cortisol-dependent I ext (T, O) may affect results from our model. We assume a simple form of  Figure A5: Stress timing and cortisol response. (A) A stressor of duration of 30min with magnitude Iext = 0.2 was applied mainly over the falling phase of the underlying cortisol oscillation. The first peak after the stress onset was almost unchanged, but the first nadir was elevated. (B) The same stressor used in (A) applied during the rising phase led to a significantly increased subsequent peak while the first nadir was unaffected. (C) The trajectory of the system (red) is projected onto the cortisol-ACTH plane. The new limit cycle of the PA-subsystem corresponding to fixed I(t) = 1.2 is indicated by the blue curve. During stress, the trajectory of the system is attracted towards the new limit cycle. The system recovers after making a smaller cycle within the normal limit cycle, reaching a higher nadir. (D) The trajectory of the system deviates then recovers back through a trajectory above the normal limit cycle, reaching a higher peak.
will make the timing of stress onset become more relevant in predicting whether a stressor can induce transitions between normal and diseased states. Driven by the intrinsic oscillations in O(T ), I ext (T, O) will also oscillate during stress, driving the q-nullcline back and forth in the (q, c)-plane as shown in Fig. A6C. Stress-driven oscillations in the q-nullcline affect the net decrease in q, changing the position of the system on the (q, c)-plane relative to the separatrix between the normal and the diseased basins of attraction upon stress termination. Since transitions are sensitive to the position of q at stress termination, including a cortisol-dependent I ext (T, O) will make transitions more strongly dependent on the timing of stress onset. In turn, theses shifts in q-nullcline will affect the net decrease in q during stress.