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., 2016; Walker et al., 2010) 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):

where There are a total of seven parameters in the constraint defining the c-nullcline in Equation 6: q₀,q₁,q₂,I,g_c,max,n, and k. Our first approach is to study how time-dependent changes in the synaptic input I affect the dynamics of the system. Because I and q₀ form a product, changes in I and q₀ are equivalent; therefore we only consider changes q₀,q₁,q₂,μ_c,n, and k. In Figure 2, we vary one parameter at a time, while all others are kept fixed at nondimensional reference values q₀ = 28.0,q₁ = 0.04,q₂ = 1.8,I = 1,g_c,max = 42,n = 5, and k = 2.83, as detailed in Kim et al. (2016) and in Appendix A of Kim et al. (2017).

Increasing q₀ 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 q₁ 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 q₂ (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 q₀/q₂. When g_{c, 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 g_c(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)

Equation 8 does not directly relate csto c but couples the two through o, which exhibits oscillatory behavior for physiological values of c. To directly relate c to cs, we average o over one full cycle of its oscillation, the values of which are fully determined by c through the PA subsystem, as illustrated in the previous section. The nullcline relation Equation 8 can thus be approximated and rewritten as where The term o^*(θ;c) represents the trajectory of o(t) with phase θ along the limit cycle of the PA subsystem defined by a given c. Therefore the cs-nullcline depends on the dynamics of the PA subsystem (Equations 3–4) and its parameters b,p₂,p₃,p₄,p₅, and p₆. We can now vary one of these parameters while fixing the others to (Kim et al., 2016; Walker et al., 2010) b = 0.6,p₂ = 15,p₃ = 7.2,p₄ = 0.05,p₅ = 0.11, and p₆ = 2.9, as illustrated in Figure 3. The part of the cs-nullcline that is approximated by averaging $c_{\infty }(o)$ over a full period of the limit cycle is indicated by dashed segments.

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 p₂ or p₃ lengthens the oscillatory regime of the cs-nullcline and shifts it toward the right. Increasing p₄ elongates the oscillatory regime, while increasing p₅ shrinks it. In both cases, a relatively small horizontal shift of the cs-nullcline is observed (Figures 3D and3E). Finally, increasing p₆ 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.

Physiological Changes

In a previous study on male rats (Patricia, Raymond, Milot, Merali, & Plamondon, 2014), 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 ⁺² (Hazell et al., 2012). Consequently, the maximum of the CRH secretion rate is increased, which corresponds to higher values of q₀ in our model. Figure 2A implies that an increase in q₀ 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. (2014). 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., 2004a; Yehuda, Golier, Yang, & Tischler, 2004b) 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 p₂ in Equation 3, which controls the sensitivity of the feedback function f_a(or) = 1/(1 + p₂(or)). Our results are shown in Figure 3B and indicate that increasing p₂ 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., 2007; Kim et al., 2016) 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, 2000). 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., 2004a). The difference in s was interpreted as due to heightened sensitivity of the negative feedback of cortisol in the pituitary (Stein, Yehuda, Koverola, & Hanna, 1997; Yehuda et al., 2004b). 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

Here oexo(t) denotes the concentration of circulating DEX converted to equivalents in cortisol concentrations based on the relative potency of DEX from Steven (1997). Note that we do not include DEX in Equation 12 because it was shown that DEX retention is much lower in the brain (Steckle, Kalin, Reul, & Hans, 2005) and that DEX does not affect GR in brain tissue (Cole et al., 2000). In a typical DEX challenge test, cortisol levels are measured in the morning (8:00 a.m.) to obtain basal pre-DEX cortisol levels. DEX is then orally administered, typically at night (11:00 p.m.), and post-DEX cortisol levels are measured again the following morning (8:00 a.m.). We assume circulating DEX levels follow a simple pharmacokinetic law: where p₇ represents the decay rate of DEX relative to the decay rate of cortisol. For simplicity, we use a rectangle function for Π_o(t). The width (30 min) of Π_o(t) is estimated based on the at-peak concentration of DEX (Perez, Rogers, Smith, & Weisman, 1998) and the height (2 in nondimensionalized units of o) is set to match the dosage used in the experiment (Yehuda et al., 2004a). Based on the half-life of DEX (∼ 240 min; Perez et al., 1998) and cortisol (∼ 7.2 min; Lightman et al., 2008), we set p₇ = 0.03.

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 N_pre and D_pre 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 s_N = 74% and s_D = 83% for normal and PTSD subjects (Yehuda et al., 2004a). 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 s_N = 76% for normal subjects and s_D = 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 N_post and D_post for normal and PTSD subjects, respectively) defined by the new cs-nullcline under DEX administration.

We use the decrease in f_a(or) = 1/(1 + p₂(or)) in Equation 3 due to DEX administra tion as a measure of suppression of pituitary activity under the challenge test. Recall that f_a(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 f_a(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 $o_{\alpha }^{*}$ and $r_{\alpha }^{*}$, where α = N,D indicate normal and PTSD states. In Figure 4C, pre-DEX values of f_a(〈or〉_α) and post-DEX values of $f_a((o_{\alpha }^{*}+o_{\text{exo}})r_{\alpha }^{*})$ are plotted for normal (α = N) and PTSD (α = D) subjects. The decreases in f_a(or) denoted as Δf_a,α 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 f_a(or): Decreases in f_a(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 + o_exo because f_a(or) is convex (f_a(or)″ > 0) and decreasing (f_a(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. (2016). 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., 2009), 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., 2004a,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., 2009). It was thus concluded (Radant et al., 2009) 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 g_r(or) and f_a(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., 2009) 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 p₂ and p₄ 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., 2017, 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. (2009). 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

where aexo(t) denotes the concentration of exogenous cosyntropin in circulation. We model the dynamics of aexo(t) using a pharmacokinetic description similar to that used for DEX: Here p₈ represents the clearance rate of cosyntropin scaled by the clearance rate of cortisol and Π_a(t) is a rectangle function. We set p₈ = 1.7 based on the 4.1-min half-life of exogenous ACTH (Matsuyama, Ruhmann-Wennhold, Johnson, & Nelson, 1972). We compare cortisol’s response to cosyntropin administered near the nadir of its ultradian oscillation for hyporeactive and near the peak for normal adrenal glands in Figure 5B. In this case, the response of a hyporeactive subject is higher than that of a normal subject. To study the dependence of cortisol response on the timing of cosyntropin administration, in Figure 5C, we plot the peak levels of cortisol under cosyntropin administration against the phase of the oscillation at which the administration took place. Overall, cortisol response is predicted to be slightly greater in normal subjects, with the maximum peak response of a hypo-sensitive subject (dashed red in Figure 5C) comparable to the minimum peak response of a normal subject (solid red in Figure 5C). This prediction implies that, depending on the phase of the intrinsic oscillation at the time of cosyntropin administration, the response of a hyporeactive adrenal gland could be greater than normal. However, the corresponding mean cortisol response should be lower if a sufficiently large sample size is used without controlling for the phase of cortisol at the time of cosyntropin administration.

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. (2009) 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., 2009) 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 p₂) 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) = I₀ + I_ext(t), where I_ext(t) represents the change in synaptic input received by the PVN neurons.

During the second stage of our proposed test, we assume I_ext(t) to be a positive constant when the external (psychological) stressor is on, and zero while off. In our thought experiment, we turn I_ext 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 I_ext 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 I_ext perturbation increases the secretion rate of CRH of the PVN neurons, effectively changing q₀ in Equation 2. We have previously shown that increasing q₀ 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 I_ext = 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.