The hypothalamic–pituitary–adrenal (HPA) axis is a neuroendocrine system that regulates the secretion of cortisol by the adrenal cortex in response to stressors. To understand how the HPA axis functions, especially in the context of posttraumatic stress disorder (PTSD), a number of challenge tests have been developed. These tests typically involve measuring changes in key endogenous hormone levels after the administration of their synthetic analogues in both PTSD and normal subjects. In this article, we reinterpret the outcomes and dynamics of some of the challenge tests through mechanistic models of HPA axis dynamics.
Cortisol (a type of glucocorticoid steroid hormone) is a “stress hormone” that regulates or supports a variety of important cardiovascular, metabolic, immunologic, and homeostatic functions (Watson, Brüne, & Bradley, ). As is typical for hormones, maintaining cortisol concentrations within appropriate ranges during both stress response and in the basal state is essential for normal physiological function. A stress response is typically initiated when neurons in the paraventricular nucleus (PVN) of the hypothalamus receive increased synaptic inputs from various regions of the brain, each containing information about certain types of stressors. These synaptic inputs induce the PVN neurons to release corticotropin-releasing peptide hormone (CRH) into the portal blood vessel connecting the hypothalamus to the anterior pituitary. Released CRH travels to the anterior pituitary and activates the secretion of adrenocorticotropin hormone (ACTH). ACTH travels via the bloodstream to the adrenal cortex, located above the kidneys, where it stimulates cortisol secretion. Finally, cortisol travels back to both the pituitary and the hypothalamus to suppress their activities, completing the negative feedback loop and returning cortisol to a basal level. Both ACTH and cortisol exhibit ultradian (hourly) and circadian (daily) oscillations. The basic interactions regulating the HPA axis dynamics are summarized in Figure 1. If the HPA axis is dysregulated, once stimulated, cortisol may fail to return to basal levels, disrupting other functions and causing comorbidities. For example, excessive cortisol is associated with major depressive disorder (MDD; Gold & Chrousos, ), while low cortisol is generally reported among PTSD patients. Clini cal observations confirm lower than normal cortisol levels in the urine of PTSD patients col lected over a 24-hour period (Mason, Giller, Kosten, Ostroff, & Podd, ; Yehuda, Kahana, Binder-Brynes, Southwick, et al., ; Yehuda et al., ) and in blood samples drawn at 15-min intervals over 24 hours (Bremner, Vermetten, & Kelley, ). In one measurement, cortisol levels collected in urine of PTSD subjects were 40.9 ± 12.3 μg/day and appreciably lower than the 62.8 ± 22.2 μg/day collected in urine of normal subjects (Yehuda et al., ). Blood samples also showed that plasma cortisol levels were consistently lower among PTSD patients over the 24-hour period and significantly lower in the afternoon (Bremner et al., ).
Despite the general trend of lower cortisol under PTSD, a few contradicting reports may be found in the literature. Most of them are surveyed in a recent meta-analysis (Meewisse, Reitsma, Vries, Gersons, & Olff, ). For example, urine measurements in Rasmusson et al. () showed a marginally significant increase of cortisol in premenopausal women with PTSD; the basal cortisol level in Yehuda, Golier, Yang, and Tischler () showed no difference between PTSD and control groups. The overall pooled estimate of the mean cortisol level from the meta-analysis (Meewisse et al., ) indicated that cortisol was lower within the PTSD group (standardized mean difference =−0.12 with p value 0.24) but concluded that the difference was not significant. On the other hand, the difference between PTSD and control subjects was shown to be significant in some subgroups, such as females and those exposed to certain types of trauma (Meewisse et al., ). Most importantly, the meta-analysis found that cortisol levels were significantly lower among PTSD patients when compared to subjects who had not been previously exposed to trauma (standardized mean difference =−0.35 with p value 0.007).
Another, more recent, meta-analysis (Morris, Hellman, Abelson, & Rao, ) studied the relationship between cortisol levels measured within 72 hours from trauma exposure and subsequent PTSD, symptoms. Overall, lower cortisol levels were shown to be associated with PTSD but not in a statistically significant way (overall correlation r = −0.07 with p value 0.449). However, as with the Meewisse meta-analysis (Meewisse et al., ), cortisol levels for subjects aged 30 years or older showed a significant negative correlation to PTSD (correlation r = −0.19 with p value 0.053). Lower cortisol levels may be associated with PTSD for other subgroups. One of the goals of our study is to propose a theory to explain how trauma exposure can lead to altered cortisol dynamics and understand why certain subgroups are likely to express lower posttraumatic cortisol levels.
Previous mathematical models (Andersen, Vinther, & Ottesen, ; S. Gupta, Aslakson, Gurbaxani, & Vernon, ) also attributed alterations in cortisol levels to bistability but did not consider ultradian oscillations known to arise in HPA dynamics. Alternatively, Gudmand- Hoeyer, Timmermann, and Ottesen (), Sriram, Rodriguez-Fernandez, and Doyle (), and Bangsgaard and Ottesen () attributed altered levels in observed cortisol levels in PTSD and depression to changes in parameter values. In particular, Bangsgaard and Ottesen () recently found that their mechanistic model yielded good fits to measured cortisol and ACTH levels of depressed and control subjects, both individually and as a group, and identified parameters relevant in diagnosing depression. Other models have focused on showing whether endogenous ultradian oscillations are allowed in a class of models of the HPA axis (Savić, Jelić, & Burić, ; Vinther, Andersen, & Ottesen, ) or have investigated the feedback mechanism between the HPA axis and the memory system (Savić, Knežević, & Opačić, ). Further discussion on the structure and features of previous mathematical models can be found in Kim, D’Orsogna, and Chou () and Bertram ().
There is no consensus on which mathematical model is mechanistically the most accurate. Hosseinichimeh, Rahmandad, and Wittenborn () selected five publications (Andersen et al., ; Conrad, Hubold, Fischer, & Peters, ; Jelić, Smiljana, Čupić, & Kolar-Anić, ; Sriram et al., ; Vinther et al., ) and compared model predictions with cortisol and ACTH levels obtained from 17 healthy individuals. Measured ACTH levels were substituted into the equation for cortisol for each of the models, and cortisol data were generated. Later, the role of cortisol and ACTH was reversed: Measured cortisol levels were used to numerically determine ACTH. Simulated data and actual measurements were then compared. According to Hosseinichimeh et al. (), none of the five models analyzed provided a good fit with data: The best fitting model, by Andersen et al. (), still yielded a difference of an order of magnitude in cortisol between simulations and data. Moreover, Vinther et al. () and Andersen et al. () proved that oscillations could not arise within a class of mathematical models, including their own, without the inclusion of additional mechanisms.
Attempts at better calibrating the model for a more accurate description of the circadian rhythm yielded some improvement, with the mean percentage error decreasing from 94% to 66%. It is important to note that when computing the goodness of fit on the data of 17 subjects, parameters were fixed at the original values provided by the authors. Individual differences and/or variations in environments were not taken into account. These attempts at matching models with actual data show how challenging the task actually is. Many variables must be taken into account, including personalized parameter sets, circadian and ultradian rhythms, external inputs, and other factors, some of which may still be unidentified.
The focus of our work is not to exactly match data to a model but rather to provide a mechanistic framework to better understand the interplay between various components of the HPA axis and its overall qualitative behaviors. Our goal is to gain insight into how variations in parameters may affect the HPA equilibrium for each individual and how transitions may be induced between multiple equilibria if they exist. Moreover, truly validating our model would require extensive data measured both before and after exposure to trauma. Such data do not exist to the best of our knowledge, and we hope that our work could motivate new, specific measurements that may yield a complete picture of the HPA axis and its dynamics.
Previously, we developed a dynamical model of the HPA axis to describe the interactions among the key hormones mentioned earlier and the glucocorticoid receptor (GR) that mediates feedback activity of cortisol (Kim et al., ). In the current work, we use and adapt this model to specifically consider perturbations representing externally or internally induced changes as well as currently employed challenge test protocols and methodologies. For appropriate sets of parameters, our model exhibits bistability with two attracting limit cycles over which cortisol and ACTH oscillate with an ultradian (hourly) rhythm. Of the two distinct oscillating states, the one with lower averaged cortisol level was characterized as a diseased state, and the one with higher cortisol level was associated with the normal state. Within our model, prolonged stress-induced secretion of CRH can trigger transitions between normal and diseased states, suggesting a possible mechanism leading to the emergence of a low cortisol state after a traumatic experience. Within our previous analysis, all parameters were fixed, except for the one representing the net synaptic input experienced by the PVN neurons. In reality, many of the parameters may vary due to aging (D. Gupta & Morley, ), life experience (Dinces, Romeo, McEwen, & Tang, ), gender differences (Seeman, Singer, Wilkinson, & McEwen, ; Uhart, Chong, Oswald, Lin, & Wand, ), the use of medication (Simunkova et al., ), and especially in response to injuries. For example, selective serotonin reuptake inhibitor (SSRI) treatment was shown to increase hippocampus volume (Thomaes et al., ). This may in turn modulate the tonic inhibition exerted by the hippocampus on the PVN and subsequently change the baseline synaptic input to CRH neurons in the PVN. On the other hand, traumatic brain injury (TBI) can produce diffuse axonal injuries (DAI; Sharp, Scott, & Leech, ) that can affect the synaptic strength or activation pattern of the hypothalamus (McCullers, Sullivan, Scheff, & Herman, ). Recent studies on the influence of maternal HPA function in neonatal rats have suggested that early life experience can have long-lasting impacts on responsivity of HPA axis activation (Dinces et al., ). Finally, the strength of the negative feedback of cortisol in the pituitary is hypothesized to be enhanced under PTSD (Yehuda, Halligan, Grossman, Golier, & Wong, ). These examples are only a few of many conditions that may be associated with anatomical/physiological changes that can be interpreted, and thus incorporated, as alterations in the parameters of our model (Kim et al., ).
For a more accurate description of the physiology in the current work, we include the self-upregulation property of CRH through the auto/paracrine function of pituitary cells. This modification of the original model does not change its overall qualitative features. A detailed comparison between the original model (Kim et al., ) and the current model is included in Kim et al. (, Appendix B). Here we briefly analyze the nullcline structure of the new model as a function of parameters and analyze the conditions under which normal–diseased bistability arises, whether or not stress-induced transitions are possible between the bistable states and whether or not the experimentally observed ultradian oscillations can be reproduced. We also systematically investigate how the dynamics depend on relevant parameter changes associated with, for instance, cerebral ischemia. Understanding how changes in parameters affect transient and long-term behavior of the HPA axis could guide therapeutic strategies.
We specifically include perturbations that represent existing pharmacological challenge tests used to assess HPA axis function. For example, dexamethasone (DEX) has been used as an exogenous steroid that suppresses ACTH release in the pituitary. DEX suppression tests have shown that the relative decrease in cortisol is greater in PTSD patients than in those in a control group. The current viewpoint is that the greater decrease in cortisol seen in PTSD subjects is due to an enhanced negative feedback of cortisol. Within a mathematical framework, this effect would typically be modeled by changing the physiological parameters describing the negative feedback. Our model allows for a novel mechanism: Cortisol suppression may be due to transitions between bistable steady states without necessarily invoking parameter changes. In this picture, disruptions to the HPA axis could be alleviated by externally controllable inputs that lead to transitions between the bistable states rather than by permanently altering specific physiological parameters.
Finally, we demonstrate how our mathematical analysis can help address previously unresolved observations of challenge experiments. ACTH stimulation tests performed to assess adrenal sensitivity in PTSD subjects showed slightly increased cortisol response to ACTH administration among PTSD subjects (Radant et al., ). This observation contradicts the authors’ speculation of decreased adrenal sensitivity in PTSD subjects. We show that this experimental result may be in fact consistent with the decreased adrenal sensitivity hypothesis and emphasize that the interplay between all components of the HPA axis must be taken into account to fully understand its overall dynamics.
Model and Methods
Our nondimensionalized model of the HPA axis based on the interactions shown in Figure 1 consists of a system of delay differential equations, as follows:
Parameter | Description | Effects on nullclines (when increased) |
---|---|---|
q0 | maximum release rate of CRH in basal state | an upward shift of the upper branch and a leftward shift of both knees of the c-nullcline |
q1 | circulating CRH for half-maximum self-upregulation | an upward shift of the upper branch and a leftward shift of both knees of the c-nullcline |
q2 | ratio of CRH and cortisol decay rates | a downward shift of the upper branch and a rightward shift of both knees of the c-nullcline |
gc,max | maximum auto/paracrine effect of CRH in the pituitary | a rightward shift of the lower and upper knees of the c-nullcline and an upward shift of the upper branch |
n | Hill coefficient of gc(c) describing the self-upregulation of CRH | a leftward shift of the left knee and a rightward shift of the right knee of the c-nullcline |
k | relates stored CRH to CRH release rate | a leftward shift of the middle branch of the c-nullcline and an upward shift of the lower branch of the c-nullcline |
b | relates cortisol to stored CRH level | a leftward shift of the middle branch of the cs-nullcline |
p2 | (or)-complex level for half-maximum negative feedback | a rightward shift and elongation of the oscillatory regime of the cs-nullcline |
p3 | ratio of ACTH and cortisol decay rates | a rightward shift and elongation of the oscillatory regime of the cs-nullcline |
p4 | (or)-complex level for half-maximum positive feedback on r production | elongation of the oscillatory branch of the cs-nullcline |
p5 | basal GR production rate by pituitary | shortening of the oscillatory branch of the cs-nullcline |
p6 | ratio of GR and cortisol decay rates | a rightward shift and shortening of the oscillatory branch of the cs-nullcline |
td | delay in the adrenal cortex in response to ACTH | elongation of the oscillatory branch of the cs-nullcline (Walker, Terry, & Lightman, ). |
The introduction of cs is the most distinctive feature of our model in comparison to others (Andersen et al., ; S. Gupta et al., ; Sriram et al., ; Walker et al., ) and allows us to more realistically model aspects of CRH dynamics that occur on different timescales. The two variables, cs and c, distinguish the two stages of the CRH secretion process: (a) the “slow” synthesis and packaging process of CRH peptides as regulated by cortisol and (b) the “fast” release process of CRH into the median eminence governed by synaptic activities and nongenomic effects of cortisol. The constant tc reflects the slow timescale (minutes to hours) over which the amount of stored CRH, cs, relaxes toward a target value relative to the timescale of CRH release (milliseconds). The target value is set by circulating cortisol levels o(t) and embodies its negative feedback on CRH synthesis.
The two state variables, cs and c, are coupled by Equation 2 through a saturating and monotonically increasing function h(cs) = 1 − e−kcs so that the average release rate of CRH increases with more stored CRH available. Self-upregulation of CRH release by the hypothalamic PVN neurons was included in our previous model (Kim et al., ) and was based on an experiment (Ono, Castro, & McCann, ) that showed that when CRH is injected into the third ventricle, PVN neurons increase their rate of CRH release into the median eminence. Thus this experiment indicates self-upregulation only if the CRH released in the median eminence directly and immediately increases the CRH level in the cerebrospinal fluid (CSF) filling the third ventricle. This is certainly a plausible mechanism. Yet, because a direct connection between the two pools of CRH is not yet well established, we improve our model by considering another relevant source of CRH activity.
Other measurements (Giraldi & Cavagnini, ) have demonstrated that CRH is also produced and secreted by the pituitary itself and that ACTH secretion in the anterior pituitary is upregulated in an auto/paracrine fashion by CRH secretion. We correspondingly revise the model presented in Kim et al. () and in Equation 2 of the current work to include both CRH secretion due to the stimulation of the PVN neurons, as modulated by the synaptic input I(t) and the amount of stored CRH cs, and CRH secretion due to the auto/paracrine activity of the anterior pituitary, as described by the increasing Hill-type function gc(c). Its amplitude gc,max can be estimated by isolating the self-upregulated CRH activity of the pituitary (Giraldi & Cavagnini, ). Together with the decaying term, our new Equation 2 describes the time rate of change of the total CRH concentration that can influence the anterior pituitary.
To include changes in synaptic input under stress, I(t) is modeled as a time-dependent parameter I(t) = I0 + Iext(t), where I0 is the minimum basal input level and Iext(t) is the increase in the synaptic input induced by external stress.
The concentrations of ACTH a, cortisol in circulation o, and the availability of glucocorticoid receptor GR in the pituitary r obey Equations 3, 4, and 5, respectively, and comprise the pituitary–adrenal (PA) subsystem, in which c can be viewed as a control parameter. The nonlinear multistate dynamics are defined by the decreasing and saturating regulation terms fa(or) and gr(or), respectively. The characteristic timescale constant td is normalized by the decay rate of cortisol and is proportional to the time required for cortisol synthesis and release by the adrenal cortex after stimulation by ACTH. More details of the model, its dimensional form, and the choice of the functional forms used in Equations 1–dcdt=q0I(t)(1−e−kcs)︸h(cs)+gc,max(q1c)n1+(q1c)n︸gc(c)−q2c,(2)dadt=c11+p2(or)︸fa(or)−p3a,(3)drdt=(or)2p4+(or)2+p5︸gr(or)−p6r,(4)5 can be found in Appendix A of Kim et al. (). A comprehensive list of all parameters and their descriptions is included in Table 1.
Walker et al. () showed that the PA subsystem (Equations 3–5) exhibits a limit cycle for a range of fixed time delay, td. For concreteness, we consider a time delay of 15 min (in dimensionless units, td = 1.44) for the rest of this article. Moreover, the amplitude and the frequency of the limit cycle were shown to depend continuously on c, so that the limiting behavior of the PA subsystem could be unequivocally determined by the value of c. The separation of the two timescales—the faster timescale of the (a,r,o) limit cycle in the PA subsystem and the slower timescale governing the CRH synthesis process—allows us to study the dynamics of the entire system (Equations 3–drdt=(or)2p4+(or)2+p5︸gr(or)−p6r,(4)5) by confining our analysis to the reduced system on the (cs,c)-phase plane. This means that the long-term behavior of the entire system can be characterized by the structures of the nullclines of cs and c projected onto the (cs,c)-plane. In particular, for certain sets of parameters, the nullcline structure exhibits bistable fixed points on the (cs,c)-plane that can be characterized as the diseased and normal modes of the PA subsystem, each marked by ultradian oscillations and distinct mean cortisol levels.
In the following sections, we outline the analysis of the present model described by Equations 1–dcdt=q0I(t)(1−e−kcs)︸h(cs)+gc,max(q1c)n1+(q1c)n︸gc(c)−q2c,(2)dadt=c11+p2(or)︸fa(or)−p3a,(3)drdt=(or)2p4+(or)2+p5︸gr(or)−p6r,(4)5.
Parameter Dependencies
In this section, we investigate how the model (Equations 1–dcdt=q0I(t)(1−e−kcs)︸h(cs)+gc,max(q1c)n1+(q1c)n︸gc(c)−q2c,(2)dadt=c11+p2(or)︸fa(or)−p3a,(3)drdt=(or)2p4+(or)2+p5︸gr(or)−p6r,(4)5) behaves as their parameters are varied. This is important because some of the ones used in previous studies (Kim et al., ; Walker et al., ) were estimated from insufficient data or arbitrarily selected. We first examine the robustness of the long-term behavior of our model to changes in individual parameters. Because the long-term behavior of the model can be characterized by the nullcline structures of Equations 1–dcdt=q0I(t)(1−e−kcs)︸h(cs)+gc,max(q1c)n1+(q1c)n︸gc(c)−q2c,(2)dadt=c11+p2(or)︸fa(or)−p3a,(3)drdt=(or)2p4+(or)2+p5︸gr(or)−p6r,(4)5 projected onto the (cs,c)-plane, in the following subsections, we study how parameter changes affect the cs- and c-nullclines. Descriptions of each parameter and their effects on the nullclines are listed in Table 1.
Nullcline Analysis: c-Nullcline
The c-nullcline is defined as the set of (cs,c) that satisfies the following relation (obtained by setting Equation 2 equal to zero):
Increasing q0 narrows and shifts the bistable region in the (cs,c)-plane toward lower values of cs, while extending the nullclines to larger values of c (Figure 2A). When q1 is increased, the bistable region also narrows and shifts toward lower cs values, but the range of the corresponding c values on the nullclines does not change significantly (Figure 2B). Increasing q2 (Figure 2C) appears to have the opposite effect: The bistable region is widened and the range of corresponding c values decreases. This behavior is expected, because it can be shown that the roots of Equation 6 depend on the ratio q0/q2. When gc, max is increased, the upper branch of c-nullcline shifts toward higher values of c, and the bistable regime moves toward the lower cs values (Figure 2D). The Hill coefficient n of gc(c) (in Equation 2) also exhibits a saddle-node bifurcation at a critical point 4 < n* < 5 (Figure 2E). Once bistability emerges, the upper and lower knees of the c-nullcline shift toward lower and higher values of cs as n increases. The bistable regime in cs is elongated as the two knees shift in the opposite direction from each other. Lastly, increasing k (Figure 2E) shifts the bistable region toward lower cs values without appreciably changing the values of c over which bistability exists.
By understanding how each parameter affects the c-nullcline, we can identify and predict which parameter changes disrupt HPA axis function. For example, we predict that increasing k will make the lower cortisol state less accessible because the range of cs covered by the lower branch of the c-nullcline narrows as k increases. Furthermore, we can use our results to interpret experimental reports of perturbations to the HPA axis. Any physiological disruption observed or associated with changes in long-term HPA function can be mapped to relevant parameter changes in our model.
Nullcline Analysis: cs-Nullcline
The cs-nullcline is defined as the set of (cs,c) that satisfies the relation (obtained by setting Equation 1 equal to zero)
Figure 3A shows that increasing b shifts the cs-nullcline to the left in the (cs,c)-plane. This shift toward lower values of stored CRH, cs, is expected, because higher b corresponds to stronger cortisol-induced suppression of the CRH synthesis. On the other hand, as shown in Figures 3B and3C, increasing p2 or p3 lengthens the oscillatory regime of the cs-nullcline and shifts it toward the right. Increasing p4 elongates the oscillatory regime, while increasing p5 shrinks it. In both cases, a relatively small horizontal shift of the cs-nullcline is observed (Figures 3D and3E). Finally, increasing p6 increases the upper limit of c of the oscillatory branch and shifts it to smaller values of cs, as shown in Figure 3F.
How the long-term behavior of the HPA axis varies as physiological parameters change is now mapped out. We can use these results to predict what long-term changes would emerge under various pathological conditions, as illustrated in the next section. Note that many physiological disruptions may alter more than one parameter and that changing multiple parameters simultaneously can lead to a qualitatively different deformation of the cs- and c-nullcline. The effect of multiple parameter changes is beyond the scope of this article. However, we can easily extend our analysis to address the issue by altering all the parameters of interest when generating the nullclines.
Results and Discussion
In the previous section, we have shown that the long-term dynamics of the system are determined by the crossing of the c- and cs-nullclines defined by Equations 6 and 9. Thus we only need to understand how the intersections change upon varying model parameters to study how alterations in parameters will affect the system. As a demonstration, we will predict how the parameter change observed in cerebral ischemia, a condition in which there is insufficient blood flow to the brain, would affect cortisol levels using the nullcline analysis from the previous section. Our predictions will be then compared to experimental observations, showing good agreement between them.
Finally, having studied how the HPA axis responds to parameter changes, we use our model to better understand a series of pharmacological challenge tests used to assess HPA function in PTSD. We reexamine some current interpretations by comparing them to our predictions and uncover unforeseen intricacies underlying the response of the HPA axis. Our analysis allows us to present novel, alternative interpretations of the observed challenge test experimental data. Indeed, this is one of the main virtues of mathematical models: They allow for probing complex, nonlinear dependencies and for crafting nontrivial predictions. Although the model presented in Equations 1–dcdt=q0I(t)(1−e−kcs)︸h(cs)+gc,max(q1c)n1+(q1c)n︸gc(c)−q2c,(2)dadt=c11+p2(or)︸fa(or)−p3a,(3)drdt=(or)2p4+(or)2+p5︸gr(or)−p6r,(4)5 can be used to study any condition or experimental protocol that involves parameter alterations in the HPA axis, we focus, in this article, on experiments related to PTSD.
Physiological Changes
In a previous study on male rats (Patricia, Raymond, Milot, Merali, & Plamondon, ), a long-lasting (30 days) increase in the expression of CRH receptor of Type 1 (CRHR-1) was observed in the PVN after global cerebral ischemia. The receptor is a G-protein coupled receptor that stimulates CRH secretion of the PVN neurons by increasing intracellular Ca +2 (Hazell et al., ). Consequently, the maximum of the CRH secretion rate is increased, which corresponds to higher values of q0 in our model. Figure 2A implies that an increase in q0 will shift the c-nullcline leftward in the (cs,c)-plane, which in turn will shift the intersection of the cs- and c-nullclines toward higher c values. This upward shift predicts elevated CRH and cortisol level after ischemia, consistent with reports in Patricia et al. (). This example demonstrates how the nullcline structure analysis can provide insights in analyzing experimental observations. We will look at experimental results on PTSD subjects in the following subsections in a similar manner.
In particular, we use results from the previous section to revisit the hypothesis (Yehuda et al., ; Yehuda, Golier, Yang, & Tischler, ) of cortisol exerting an enhanced negative feedback on the pituitary under PTSD, which is commonly invoked to explain the low cortisol levels observed in PTSD patients. We can mathematically model the enhanced negative feedback of cortisol by increasing p2 in Equation 3, which controls the sensitivity of the feedback function fa(or) = 1/(1 + p2(or)). Our results are shown in Figure 3B and indicate that increasing p2 shifts the cs-nullcline toward the right, away from the lower branch of the c-nullcline. We thus predict that enhanced inhibition of the pituitary by cortisol will typically increase cortisol levels. Parameter changes reflecting enhanced inhibition in this class of models cannot explain the lower cortisol levels often associated with PTSD. However, earlier models (S. Gupta et al., ; Kim et al., ) and our current work allow for bistable steady states to arise, providing an alternative mechanism by which PTSD subjects may manifest low baseline cortisol. In particular, low cortisol diseased states arise as a transition from one stable state to another rather than from permanent physiological changes. To further investigate if this perspective is consistent with other experimental observations, we revisit a study in which dexamethasone (DEX) suppression was tested on PTSD patients.
Dexamethasone Suppression Test (DST)
The dexamethasone suppression test (DST) is a pharmacological challenge test typically used to identify the cause of abnormal cortisol levels observed in diseases such as Cushing syndrome. In DST, a cortisol analogue (dexamethasone, DEX) is used to probe the negative feedback of cortisol on ACTH secretion in the pituitary. In particular, some studies have used DST to test whether lower basal cortisol levels in PTSD result from an enhanced negative feedback of cortisol on pituitary activity. DEX is a synthetic glucocorticoid compound that suppresses ACTH secretion, and subsequently cortisol secretion, when it binds to glucocorticoid receptors (GR) in the pituitary (Cole, Kim, Kalman, & Spencer, ). In DEX suppression studies, cortisol levels are measured pre-DEX and post-DEX and used to calculate the percentage suppression of cortisol, defined as
The mean percentage suppression of cortisol s was shown to be greater in PTSD subjects (s = 83%) than in the control group without PTSD (s = 74%Yehuda et al., ). The difference in s was interpreted as due to heightened sensitivity of the negative feedback of cortisol in the pituitary (Stein, Yehuda, Koverola, & Hanna, ; Yehuda et al., ). This interpretation implicitly assumes that the suppression effect of cortisol (and dexamethasone) is directly proportional to its concentration. The linear relationship is convenient but neglects the contribution of other components of the system that can interact with the suppression activity to yield a more complex, possibly nonlinear relationship. We use our mathematical model to examine how exogenous DST affects the dynamics of the HPA axis. We model DEX administration by replacing o(t) with o(t) + oexo(t) in Equations 3 and 4 to yield
The numerical solutions of cortisol levels during DEX suppression test as modeled by Equations 12–dcdt=q0I(t)h(cs)+gc(c)−q2c,(13)dadt=c1+p2((o+oexo(t))r)−p3a,(14)drdt=((o+oexo(t))r)2p4+((o+oexo(t))r)2+p5−p6r,(15)16 are shown in Figures 4A and4B for normal and PTSD subjects, respectively. The parameters used in both figures are identical: Normal and PTSD subjects are characterized solely by which one of the two stable states they initially reside in (cs,c) space. Initial (cs,c) values are plotted in Figure 4C and labeled Npre and Dpre for normal and PTSD subjects, respectively. We take the average of o over a full cycle of oscillation to estimate pre-DEX cortisol values 〈o〉pre,N = 63% and 〈o〉pre,D = 73% for normal and PTSD subjects. These predicted values are in qualitative agreement with the experimentally reported percentage suppression sN = 74% and sD = 83% for normal and PTSD subjects (Yehuda et al., ). Owing to the oscillatory behavior of cortisol in the basal sate, the percentage suppression depends on the phase of the oscillation at the time of pre-DEX measurement. For example, if the pre-DEX measurement time is set near the peak of the oscillations, our model predicts sN = 76% for normal subjects and sD = 81% for PTSD subjects. For a more accurate comparison between our predictions and data, the timing of DEX administration with respect to the phase of the oscillating cortisol levels should be carefully controlled in experiments.
To understand the behavior of our HPA axis model under DEX administration, we observe how oexo(t) affects the (cs,c)-nullcline structure. Because oexo(t) contributes only to the PA subsystem (Equations 3–drdt=(or)2p4+(or)2+p5︸gr(or)−p6r,(4)5), it only affects the cs-nullcline. The latter is shown in Figure 4C (dark blue) near the time of measurement, 9 hours from DEX administration. Normal and diseased states initially resting at the intersection of the cs- and c-nullclines slowly evolve toward the shifted intersections (labeled Npost and Dpost for normal and PTSD subjects, respectively) defined by the new cs-nullcline under DEX administration.
We use the decrease in fa(or) = 1/(1 + p2(or)) in Equation 3 due to DEX administra tion as a measure of suppression of pituitary activity under the challenge test. Recall that fa(or) is a modulating factor of the ACTH production rate and represents the negative feedback of the cortisol on the pituitary. For the pre-DEX value of fa(or), we take the period-averaged value of or, denoted 〈or〉N and 〈or〉D for normal and diseased initial states, respectively. Note that the new shifted cs-nullcline under DEX administration is not associated with oscillatory behavior and that the new intersections represent nonoscillating equilibria. The cortisol levels depicted in Figures 4A and4B confirm that cortisol oscillations cease under DEX suppression, which has not yet been experimentally tested. The nonoscillating equilibrium values for o and r are denoted and , where α = N,D indicate normal and PTSD states. In Figure 4C, pre-DEX values of fa(〈or〉α) and post-DEX values of are plotted for normal (α = N) and PTSD (α = D) subjects. The decreases in fa(or) denoted as Δfa,α in Figure 4D indicate that suppression is greater for PTSD subjects than for normal subjects, despite the comparable change in or values. This is because of the form of fa(or): Decreases in fa(or) due to increases in o are greater for lower initial values of or. In other words, the state with low initial or value will experience greater suppression as o increases to o + oexo because fa(or) is convex (fa(or)″ > 0) and decreasing (fa(or)′ < 0) for o ≥ 0 and r ≥ 0.
Our model provides an additional interpretation of DEX administration results in PTSD patients. The enhanced suppression effect on cortisol may be due to the intrinsic dynamics rather than parametric changes; within our model, the negative feedback effect is dependent on the state of the system at the time of DEX administration, which in turn depends on dynamics and history of the system. Thus increased suppression of cortisol levels in PTSD subjects during DEX administration may simply indicate that PTSD subjects were in the low cortisol state before the test rather than implying that permanent changes had occurred in their system in the course of developing PTSD.
This alternate explanation has direct implications for how one should design therapeutic protocols for PTSD or other stress-related disorders associated with cortisol disruption. Under the enhanced negative feedback hypothesis, therapies or medications would attempt to lower the pituitary sensitivity to cortisol to bring it back to normal levels. Our model shows that this “straightforward” approach would fail and suggests that therapies should instead focus on devising appropriate perturbations to the system. For example, applying externally controlled inputs I(t) could induce transitions back to the normal stable state, as shown in Kim et al. (). This framework is consistent with current recommendations for treatment of stress disorders via exposure therapy and cognitive-behavioral techniques in which an individual is reexposed to trauma or stress in a controlled way.
ACTH Stimulation Test
In this section, we consider the ACTH stimulation test typically used to diagnose conditions associated with insufficient adrenal activity. Cortisol levels are measured after the administration of cosyntropin, a synthetic derivative of ACTH. Exogenous ACTH stimulates cortisol secretion to the same extent of endogenous ACTH and thus effectively increases a(t) in our model. As in the analysis of the DEX suppression test, we can model cosyntropin administration by replacing a(t) with a(t) + aexo(t), where aexo(t) denotes the concentration of the cosyntropin in circulation.
In a previous study (Radant et al., ), an ACTH stimulation test was administered to PTSD and normal subjects to measure potential differences in adrenal gland response between the two groups. It was hypothesized that a bolus of cosyntropin would lead to a smaller pulse of cortisol secretion in PTSD patients due to hyporeactivity of their adrenal glands. Adrenal hyporeactivity would also suggest lower baseline cortisol levels, consistent with a number of observations (Yehuda et al., ,s). Surprisingly, the main experimental finding was that cortisol response to cosyntropin was not significantly altered in PTSD subjects. Moreover, the baseline cortisol levels observed under PTSD were actually slightly higher than normal (Radant et al., ). It was thus concluded (Radant et al., ) that either adrenal reactivity is not altered in PTSD patients or that potential alterations do not affect HPA dynamics.
This interpretation relies on the intuition that a proportional relationship exists between adrenal reactivity and cortisol response so that upon stimulation, in this case by cosyntropin, any adrenal hyporeactivity in PTSD subjects would lead to smaller cortisol increases. This picture would be valid if the stimulating activity of ACTH in the adrenal gland were isolated from other ACTH interactions within the HPA axis. However, cortisol suppresses endogenous ACTH secretion in the pituitary, which in turn indirectly reduces cortisol secretion. In addition, glucocorticoid receptor (GR) concentration is coupled to cortisol through the dependence on or in Equation 4 and thus influences the negative feedback in the pituitary. As a result of these various nonlinear interactions, particularly those embodied by gr(or) and fa(or), cortisol response to ACTH stimulation is more complex than a simple proportionality relationship. We now use our model to explore these nonlinear interactions and develop a more nuanced interpretation of the experimental observations.
In particular, we contrast and compare our results to the previously described experimental data and show that, indeed, upon taking into consideration the full dynamics of the HPA axis, laboratory observations (Radant et al., ) may be consistent with the hypothesis of reduced adrenal reactivity in PTSD subjects, in contrast to the original interpretation. Within our model, adrenal gland reactivity to ACTH determines the parameters p2 and p4 in Equations 1–dcdt=q0I(t)(1−e−kcs)︸h(cs)+gc,max(q1c)n1+(q1c)n︸gc(c)−q2c,(2)dadt=c11+p2(or)︸fa(or)−p3a,(3)drdt=(or)2p4+(or)2+p5︸gr(or)−p6r,(4)5 (refer to Kim et al., , Appendix A, for details). To model hyporeactivity, we adjust both parameters to reflect a 10% decrease in the adrenal gland reactivity in PTSD subjects. The resulting numerical solutions are plotted in Figure 5A and show that basal cortisol levels are indeed increased in the basal state of PTSD subjects, consistent with the baseline cortisol measurements in Radant et al. (). This result emphasizes that the simple intuition of a direct relationship between adrenal reactivity and cortisol response can be misleading.
To now describe the response of cortisol under an ACTH stimulation test, we modify our model by rewriting Equation 5 as
Our model predicts that reduced adrenal reactivity to ACTH will increase baseline cortisol levels. This implies that the experimental result reported in Radant et al. () may indeed reflect a reduced adrenal reactivity in PTSD subjects, amending previous interpretations. The increase in cortisol response seen among PTSD patients can also vary depending on the timing of cosyntropin administration. For example, the relative increase in cortisol is greatest upon administering cosyntropin at the nadir of the cortisol oscillation, as seen in Figures 5B and5C. Note that if experiments on normal and PTSD individuals are not phase matched, the maximum response of a normal subject may be greater than the minimal response of a PTSD subject. Thus phase matching would be required to remove confounding effects arising from the timing of ACTH/cosyntropin administration. For example, the relative increase in cortisol is greatest upon administering cosyntropin during the increasing phase near the nadir of cortisol’s intrinsic oscillation, as seen in Figure 5C. Note that if experiments on normal and PTSD individuals are not phase matched, the response of a normal subject may be less than the response of a PTSD subject, as observed from the experiment. Because the sample size used in the experiment (Radant et al., ) was small (n = 8 subjects for PTSD and n = 9 subjects for control), the increased cortisol response among PTSD subjects may be explained as a confounding effect arising from the timing of ACTH/cosyntropin administration. Phase matching would be required for a proper interpretation of the measurement.
Predictions for a New Two-Stage Challenge Test
In the previous section, we reexamined the hypothesis that changes in physiological parameters (specifically p2) can result in a greater suppression of pituitary activity by dexamethasone or cortisol. This enhanced negative feedback has been proposed as the cause for lower basal cortisol levels in PTSD subjects. Using our model and analyses, we presented an alter native hypothesis based on the nonlinear interactions in our dynamical system. Here we outline new experiments that could be used to further evaluate and distinguish these two hypotheses. If the negative feedback of cortisol on the pituitary is enhanced in PTSD subjects, their response to stressors should be reduced, especially when the pharmacological suppression is in effect. One way to further probe the system is to combine a nonpharmacological, psychological stressor with the DST. The ensuing cortisol response can be used to probe the purported enhanced negative feedback on the pituitary. The new challenge test consists of two stages: (a) Dexamethasone is administered to suppress the pituitary activity and (b) a psychological stressor applied during the suppression is in effect.
Physiological stressors are processed in various regions of the brain, and the information is transferred to the PVN neurons via synaptic inputs. Changes in the synaptic input can be represented by a perturbation in the I term in Equation 2 in the form I(t) = I0 + Iext(t), where Iext(t) represents the change in synaptic input received by the PVN neurons.
During the second stage of our proposed test, we assume Iext(t) to be a positive constant when the external (psychological) stressor is on, and zero while off. In our thought experiment, we turn Iext on for 60 min starting 9 hours after DEX administration, when post-DEX cortisol levels are usually measured. As can be seen in Figures 6A and6B, our model predicts cortisol response to be generally greater in the low-cortisol PTSD state than in the normal state, in contrast to the enhanced negative feedback hypothesis. Moreover, the peak of cortisol under PTSD can surpass that of cortisol under normal conditions if Iext is sufficiently large.
To understand this unexpected “reversed” phenomenon, consider the nullcline structure during DEX administration and after the stressor is applied (Figures 6C and6D). Upon turning the stressor on, the Iext perturbation increases the secretion rate of CRH of the PVN neurons, effectively changing q0 in Equation 2. We have previously shown that increasing q0 shifts the upper branch of the c-nullcline and moves the bistable regime toward the left in the (cs,c)-plane, as shown in Figure 2A. For a stressor with Iext = 0.5, the c-nullcline is shifted so that the PTSD state on the lower branch is no longer in the bistable regime in the (cs,c)-plane and, consequently, the PTSD state jumps to the upper branch and relaxes toward the only available steady state, approximated by the intersection of the two nullclines in Figure 6D. Meanwhile, the normal state residing on the upper branch of the c-nullcline also jumps to the shifted nullcline, but the size of the jump is significantly smaller compared to the jump from the lower branch (as seen in Figure 6D).
The proposed two-step challenge protocol shows discrepancies between the results of our model and the enhanced negative feedback hypothesis. The increase in cortisol response in a subject with lower basal cortisol level cannot be explained by an altered negative feedback strength, while it can be understood as a natural consequence of the changing dynamical structure of the system owing to perturbations in the parameter. The experiment also eliminates the confounding effect of measurements taken from hourly oscillating cortisol levels.
Summary and conclusions
The HPA axis is a dynamical system continuously evolving to meet changing physiological needs and environmental stimuli. Even at equilibrium, key hormones, such as cortisol and ACTH, exhibit ultradian oscillations. To accurately interpret the response of the HPA axis to internal, anatomical changes or external inputs, such as the injection of exogenous hormones, we need to understand how the response depends on the state of the system itself and the interplay between its different components. To this end, we developed a mathematical model of HPA dynamics where changes in parameter values and exogenous sources of key hormones could be explicitly included. Of particular importance is the understanding of how pharmacological intervention might affect the long-term dynamics of the HPA axis. Our model is useful in this respect because the effect of medication, trauma, or disease can be mapped onto changes in given parameters. We thus performed a parameter sweep and were able to predict possible modifications to long-term behaviors induced by corresponding pharmacological challenge tests.
Measurements taken during an ACTH stimulation test have shown higher baseline cortisol levels in PTSD subjects and unchanged cortisol responses to exogenous ACTH administration (Radant et al., ). This result was interpreted as evidence against altered adrenal reactivity in PTSD subjects. Upon incorporating the altered adrenal gland reactivity into our model as a parameter change, we found that data from Radant et al. () can be explained by adrenal hyporeactivity in PTSD.
Our simulations show that cortisol levels in subjects with hyporeactive adrenal glands can increase by amounts that are comparable to the increases observed in normal subjects upon the administration of exogenous ACTH. We believe that the phase of the intrinsic oscillations in cortisol at the time of exogenous ACTH administration should be controlled for a more accurate interpretation of experiments.
The most well known feature of HPA dysfunction in PTSD patients is the low secretion of cortisol. The current viewpoint is that low cortisol levels are a consequence of an enhanced negative feedback in the HPA axis (Bremner et al., ; Yehuda et al., ). This conclusion is based on DEX suppression tests that show a greater percentage suppression of cortisol among PTSD subjects. Extending our model to include DEX administration, we have provided an additional mechanism to explain low cortisol levels in PTSD, namely, that it arises as a feature of an alternate stable state of the dynamical system. The diseased low-cortisol stable state exhibits relatively greater suppression of cortisol without enhancing the sensitivity of the negative feedback of cortisol in the pituitary and without the need to invoke externally induced parameter changes. Thus enhanced cortisol suppression is an inherent feature of a stable state in our bistable model.
Our model can help us understand the many ways trauma might dysregulate cortisol dynamics and identify subgroups of PTSD patients that may require different treatment approaches. For instance, one could measure specific parameters for subjects within one of the subgroups that showed a significant difference in cortisol level (Meewisse et al., ) and verify whether the associated nullcline structure would allow for bistability. If so, our model supports possible treatments in the form of appropriate external inputs to the system to induce the transition between the stable states. To the contrary, if bistability does not arise, treatment protocols should focus on adjusting the parameter values to correct the dynamics.
Although our model has provided a mechanistic description of the HPA axis behavior under two distinct pharmacological challenge tests, it does not provide a direct explanation for the variability observed in baseline cortisol levels in PTSD patients. One possibility is that the discrepancy merely reflects different stages or effects of disruptions in the HPA axis induced by exposure to trauma.
The literature on PTSD has been driven by diagnostic criteria that rely heavily on nonquantitative and subjective self-reports. This is also true with research on neuroendocrine alterations in PTSD, where the definition of PTSD often failed to take into account important factors, such as gender and type of trauma. Such issues have likely contributed to conflicting reports on how basal cortisol levels are affected by PTSD. Our mathematical model provides a framework to help interpret how external perturbations, such as pharmacological challenge tests, lead to abnormal dynamics and long-term behavior. Analysis of our model allows us to characterize, more mechanistically, PTSD by identifying specific components of the HPA axis that can be affected by trauma or medication.
AUTHOR CONTRIBUTIONS
Lae Un Kim, Maria R. D’Orsogna, and Tom Chou designed the research. Lae Un Kim per formed the analysis and simulations. Lae Un Kim, Maria R. D’Orsogna, and Tom Chou wrote and edited the manuscript.
FUNDING INFORMATION
This work was supported by the Army Research Office via grant W911NF-14-1-0472 and the NSF through grant DMS-1516675.
ACKNOWLEDGMENT
The authors also thank the late Professor T. Minor for insightful discussions.
Additional File
The additional file for this article can be found as follows:
Supplementary FileSupplementary Material. DOI: https://doi.org/10.1162/CPSY_a_00013.s1