Mathematical modeling of depressive disorders: Circadian driving, bistability and dynamical transitions.

The hypothalamus-pituitary-adrenal (HPA) axis is a key neuroendocrine system implicated in stress response, major depression disorder, and post-traumatic stress disorder. We present a new, compact dynamical systems model for the response of the HPA axis to external stimuli, representing stressors or therapeutic intervention, in the presence of a circadian input. Our work builds upon previous HPA axis models where hormonal dynamics are separated into slow and fast components. Several simplifications allow us to derive an effective model of two equations, similar to a multiplicative-input FitzHugh-Nagumo system, where two stable states, a healthy and a diseased one, arise. We analyze the effective model in the context of state transitions driven by external shocks to the hypothalamus, but also modulated by circadian rhythms. Our analyses provide mechanistic insight into the effects of the circadian cycle on input driven transitions of the HPA axis and suggest a circadian influence on exposure or cognitive behavioral therapy in depression, or post-traumatic stress disorder treatment.

Introduction

When presented with external stimuli and challenges, organisms activate a series of physiological and behavioral actions to minimize departure from homeostasis. These body-brain responses are mostly coordinated by the hypothalamic–pituitary–adrenal (HPA) axis which controls the expression of three major stress-related hormones: CRH (corticotropin-releasing hormone), ACTH (adrenocorticotropic hormone), and glucocorticoids, through a complex set of interactions, ligand-receptor binding events, and feedback loops [1], [2]. The most prevalent glucocorticoid in humans is cortisol. The above hormones are secreted on a circadian basis: under normal conditions levels are low at night, peak in the early morning hours and decline slowly throughout the day [3]. Non-stimulated basal CRH is released in a pulsatile manner with a frequency of about two to three episodes per hour; similarly ACTH and cortisol manifest burst-like releases that follow a 60–90 min periodicity [4], [5], [6], [7], [8], [9]. As a result, an ultradian rhythm is superimposed on the circadian cycle of each hormone. Animal studies suggest that pulsatile CRH is not the source of the ultradian ACTH and cortisol rhythms[10]. On the other hand, ACTH and cortisol pulses are highly correlated, with cortisol typically trailing ACTH by a 15-min delay [11].

The HPA axis regulates many physiological processes, including digestion, the immune system, mood and emotions, sexuality and metabolism. It is a highly conserved system that is present in many vertebrate, invertebrate and mono-cellular species. Due to its central role in the body’s response to environmental stimuli and demands, disruptions to the HPA axis and related neuroendocrine activity are associated to a wide variety of pathologies, including stress-related ones [12]. For example, over-production of CRH may be linked to major depressive disorder (MDD) [13], [14], [15], [16], [17] and to anorexia nervosa where hypercortisolemia is also observed [18], [19], [20]. Patients with acute MDD have also reported elevated cortisol levels [21]. Increased ACTH and/or cortisol levels are observed in patients with Nelson’s syndrome, Cushing syndrome or Cushing disease [22], [23], [24], [25]. Cortisol production irregularities are also linked to Sheehan’s syndrome and are often observed in patients with pituitary tumors [26], [27], [28]. Those at risk for depression have been found to exhibit greater waking cortisol and a larger cortisol awakening response [29], while abnormally low cortisol levels are observed in post-traumatic stress disorder (PTSD) patients as measured by urinary and/or salivary samples [30], [31], [32], [33], [34], [35], [36], [37], [38], [39]. Addison disease, or primary adrenal insufficiency, is marked by low levels of cortisol due to adrenal gland disorders; secondary adrenal insufficiency is marked by low levels of ACTH and cortisol due to dysregulation of the pituitary gland; tertiary adrenal insufficiency is marked by low levels of CRH, ACTH, and cortisol, due to hypothalamic diseases [40]. Addison disease may also be characterized by high levels of ACTH. Here, since the adrenal cortex is unable to synthesize and/or secrete cortisol, the pituitary gland increases production of ACTH in an effort to stimulate it. Some studies associate elevated levels of CRH and cortisol to the onset of dementia and Alzheimer’s disease, whereas others observe low CRH concentrations in patients with advanced neurodegenerative conditions [41]. These seemingly contrasting findings may be due to chronic HPA axis overstimulation, where high levels of CRH and cortisol may induce permanent damage to neuronal connectivity, so that eventually secretion of CRH is reduced [42]. Other contradictory findings may emerge in surveying the literature on alterations to HPA-regulated hormonal activity. For example not all studies report cortisol deficiencies [43] among PTSD patients. These discrepancies may be due to methodological variabilities in sampling and/or timing, limited number of participants, confounding effects such as patients of different ages or who are taking medications that affect hormonal expression, and other extraneous factors.

As described above, one of the most important functions of the HPA axis is to maintain homeostasis in response to mild to severe stress via the enhanced release of CRH from the hypothalamus, initiating the CRH–ACTH–cortisol cascade. The secretion of cortisol helps the body cope with the stressor, for example by facilitating the release of glucose to activate “fight or flight” responses. However, sometimes the trauma is so severe (or the individual experiencing it is particularly susceptible) that the normal HPA axis response is disrupted and abnormal levels of stress hormones are produced on a permanent basis. Thus, acute stressful events such as major accidents, combat, assaults, natural disasters, death of a loved one, can lead to neuroendocrine dysfunction. This happens in PTSD and in MDD: both can be triggered by episodes of acute stress, and are characterized by abnormally low cortisol levels (one biomarker of PTSD) and over-expression of CRH (a biomarker of MDD). The evidence that the HPA axis is hypo/hyper active in PTSD and/or MDD comes from clinical trials, biochemical studies, functional HPA axis tests, neuro-imaging and postmortem studies. Based on these observations, we formulate our main modeling assumption, that external stressors may cause long-lasting damage.

The neuroendocrine dynamics of the HPA axis has been well-studied using mathematical models that describe the abundances of CRH, ACTH, cortisol, and glucocorticoid receptors and their dynamics in response to external inputs [44], [45], [46], [47], [48], [49], [50], [51], [52], [53], [54], [55], [56]. These models include physiologically motivated feedback loops and delays that allow for the emergence of the ultradian oscillations. The external input, which may be time dependent, is interpreted as environmental stimulus, stressor or trauma, while the circadian rhythm is typically neglected. Under given parametric regimes and modeling conditions two stable states or limit cycles emerge: one is identified as the “healthy” or “healthy” state, the other as a psychiatrically “perturbed” or “diseased” condition [47], [49]. Transitions between the two states may arise from parameter changes, which we interpret as physical injury, or as a function of external input, which we interpret as psychological trauma. In the latter case, HPA stress-related disorders may be viewed as a consequence of the bistability of the system, with the diseased state emerging as an organism’s response to psychological stress, rather than to physiological damage. Mood disorders often appear as abnormalities in CRH expression, which also manifest downstream, in ACTH and/or cortisol expression. We thus seek for bistability at the origin of the hormonal sequence that is regulated by the HPA axis, that is at the CRH level.

Separating CRH dynamics into a fast CRH release and a slow CRH synthesis yields a larger parameter space and higher likelihood for the emergence of the healthy-diseased bistable configuration [55], [56]. Apart from the action of external stressors, the time-varying input into the hypothalamus can also model cognitive behavioral therapy (CBT) or exposure therapy (ET), where patients are subjected to psychological interventions or exposed to the stressor in a controlled manner to relieve them of depression, PTSD and other stress related symptoms. Within these dynamical systems models, external inputs can trigger transitions between the healthy and diseased states depending on input duration, amplitude and time of application relative to the intrinsic oscillations of the ultradian cortisol pulses [55], [56].

The circadian cycle, which is known to affect the regulation and metabolism of several hormones, has been neglected in previous mathematical analyses. Here, we incorporate it as an input to the dynamical system. Physiologically, the relevant oscillating pacemaker is located in the suprachiasmatic nucleus (SCN), a group of neurons in the hypothalamus that respond to photosensitive ganglion cells in the retina [57], [58], [59]. The circadian rhythm can modulate input from the hypothalamus to stimulate the secretion of CRH hormones in the portal vessel which connects the hypothalamus to the anterior pituitary [60]; it may also cause oscillations in the relevant tissues and organs that comprise the HPA axis and that respond to the SCN independently of the hypothalamus.

In Section 2, we derive a dynamical systems model of the HPA axis. The “fundamental” form we present in Eqs. (12), (13) is based on existing models but it is more mathematically compact, allowing us to better study the parameter regimes that yield healthy-diseased bistability, and to include circadian driving. In Section 3 we discuss how the circadian driving affects the dynamics and equilibria of the fundamental form, through analytical and numerical analysis. In Section 4 we include various forms of external inputs representing trauma and/or exposure therapy. The complete stimulus to the HPA axis is given by the superposition of basal, rhythmic and external/reactive terms. We show that the magnitude and duration of the external input, its timing relative to the phase of the circadian cycle, and the amplitude of the circadian rhythm, strongly affect the transition between healthy and diseased states. We offer conclusions and a brief discussion in Section 5.

Reduced model without circadian drive

The model we use in this paper is a compact, reduced description of the HPA neuroendocrine system that is derived from previous work [53], [49], [54], [48], [50], [51], [52], [55], [56]. The simplifications introduced allow us to incorporate the circadian drive and perform complete analytical analyses. We begin by illustrating the full, initial model which includes five variables that represent the most relevant hormones and receptors involved in stress response [55], [56]. Later we show how our new, compact two-variable system is obtained.

The CRH neurons are the first to be activated by physiological changes such as stress and the HPA axis model we consider includes two compartments for it. One is CRH storage, that obeys relatively slow dynamics, on the order of ten to twelve hours; the other is CRH secretion, that follows relatively fast dynamics, on the order of minutes. Secreted CRH initiates a cascade of events that includes ACTH and cortisol production. Feedback loops emerge from circulating cortisol binding to glucocorticoid receptors leading to inhibition of ACTH and self-upregulation of glucocorticoid receptor production. Depending on the parameters chosen, a maximum of two stable fixed points may arise, which are interpreted as the healthy and diseased states. Specific time-dependent perturbations of the external input into the hypothalamus can drive transitions from one stable point to the other, suggesting that onset of dysfunction in stress response and its treatment can be framed in terms of dynamical bistability. In non-dimensional terms, the full model [55], [56] is given by

Here, represents the slowly evolving concentration of stored CRH in the neurons of the paraventricular nucleus of the hypothalamus, which is not directly affected by the external input, but rather depends on circulating cortisol. Stored CRH may be quickly released into the hypophyseal portal vessels before being transported to the anterior pituitary. The concentration of CRH in this “circulating” pool is denoted . The stored CRH concentration evolves on the long time scale and relaxes towards an equilibrium value that depends on cortisol levels o through an indirect negative feedback process. From experimental observations on rats so that at steady state the stored decreases with , leveling at for large cortisol levels as modulated by the parameter b [61], [62]. Aside from natural degradation modeled by the term, the dynamics of is driven by two processes. The first is the stimulus that triggers release of stored CRH: this is modeled by the product between and . In this expression is the maximum possible secretion rate of the stored pool to the circulating CRH pool which is achieved in the limit. For , there is no stored CRH to release. The coefficient k modulates the response between these two limits. The second process is self-upregulation, whereby circulating CRH stimulates further release of the same hormone [62]. Since upregulation is mediated by CRH/receptor binding, we use a Hill-type increasing function that varies between 0 and . The coefficient is the inverse of the CRH concentration that produces half maximum self-upregulation, and the Hill coefficients used are , although other choices of n will not qualitatively change our results. The concentration of ACTH generated in the pituitary gland is denoted as a. As can be seen from Eqs. 1, its production is driven by circulating c levels but inhibited by cortisol o bound to glucocorticoid receptors r through the negative feedback term . ACTH is cleared at rate . Cortisol production in the adrenal gland is driven by and follows a natural decay. Finally, glucocorticoid receptors in the anterior pituitary self-upregulate through cortisol binding to the receptors themselves [63]. The dynamics of r thus assumes a cortisol-independent production rate and cortisol-mediated generation represented by the term; the clearance rate is . The non-dimensional parameters in Eqs. 1 are highly variable and may depend on genetic traits, age, gender, and other environmental factors. Parameter choices areimportant as they may (or may not) lead to bistable solutions, representing individuals who are less (or more) resistant to dysfunctions of the HPA axis. The parameters we use in this work for Eqs. 1 lead to bistability. They are set as in [56], and are listed in Table. 1.

Table 1

Non-dimensional parameter values of the model Eqs. 1 as taken from Refs. [55], [56]. We also include non-dimensional parameter values for the 24-h circadian rhythm.

Parameter	Value	Description
c¯_∞	0.2	minimal stored baseline CRH
b	0.6	stored CRH decay as a function of cortisol
t_c	69.3	CRH biosynthesis timescale
q_0	28.0	maximum release rate of CRH in basal state
I_0	1.0	basal level of the external stimuli
k	2.83	relates stored CRH to CRH release rate
g_c,max	42.0	maximum auto/paracrine effect of CRH in the pituitary
n	5	Hill coefficient describing the self-upregulation of CRH
q1-1	25.0	circulating CRH conc. at half-maximum self-upregulation
q_2	1.8	ratio of CRH and cortisol decay rates
p2-1	0.067	or-complex conc. for half-maximum negative feedback
p_3	7.2	ratio of ACTH and cortisol decay rates
p_4	0.05	(or-complex conc)2 at half-maximum positive feedback on r production
p_5	0.11	basal GR production rate by pituitary
p_6	2.9	ratio of GR and cortisol decay rates
ω	0.045	frequency of 24-h circadian rhythm

Note that the model presented in Eqs. 1 is non-dimensionalized using the inverse of the decay rate of cortisol as the reference time scale and various parameter combinations to obtain [54], [55], [56]. Furthermore, in the original formulation of the model, a delay was included in the ACTH activated synthesis of cortisol so that [55], [56]. The delayed response of the adrenal gland to ACTH is well established and estimated to be roughly 15 min [11]. Physiologically, the delay is due to the delivery of ACTH from the pituitary to the adrenal gland and to the subsequent synthesis of cortisol. Mathematically, the delay allows for the recovery of the ultradian rhythm [54], [55], [56]. We do not consider this delay in the current work, as its mathematical implications are well understood and our goal is to understand how the circadian rhythm, which occurs on a much longer time scale than the delay, affects the dynamics. For better insight into the model, and its non-dimensionalization, the reader is referred to the analyses in [54], [55], [56].

The separation between CRH synthesis and release, processes operating on two distinct time scales, yields interesting behavior, including bistability and the possibility of transitions between the respective basins of attraction. Eqs. 1 however are mathematically cumbersome and the inclusion of a circadian rhythm would add algebraic tedium without offering clear insight. We thus present a simplified version of Eqs. 1 that exhibits the same main features but that is much simpler to analyze. To proceed, we focus only on and c in Eqs. 1 and set the subsystem to equilibrium under a given c, effectively parameterizing . We use this steady-state approximation since the subsystem evolves on a much faster timescale than .

We keep the expression for in Eqs. 1 and simplify the dynamics of c with the goal of preserving bistability. Henceforth, unless specified otherwise, we will work in space. For a given value of , equilibrium is attained when the two nullclines obtained by setting and in Eqs. 1 intersect. Since they evolve on different time scales, the first, -nullcline is sometimes referred to as the slow nullcline, the latter, the c-nullcline, as the fast one. Finally, note that the equilibrium values of define single-valued, positive, real functions of c; hence, once c is specified they are uniquely defined and no bistability can emerge at this stage. Conversely, if bistability emerges in , space it will be reflected in space as well, as parametrized by the two bistable values of c.

We begin by simplifying the right-hand side of Eqs. 1 driving the dynamics. Since we must find an explicit expression for so that our new model is self-contained in space. Under the assumption that the subsystem is equilibrated, can be found through algebraic manipulation of Eqs. 1 with the derivatives set to zero so that

As written above, o is the root of a quartic polynomial where c appears as a parameter in some of the coefficients. Writing the full exact solution to Eq. 2, while possible, is cumbersome. We note however that using realistic parameters given in Refs.[55], [56] the termis relatively small compared to others. Thus, the relation between o and c can be approximated by adding the three above terms to Eq. 2, leading to a quadratic equation in :

This simplified quadratic equation can be solved for and the resulting can be inserted into in the first of Eqs. 1 so that the evolution of is completely described in space. The approximated and full-model are shown in Fig. 1. In Appendix 6 we discuss further approximations to the term in Eqs. 1.

Fig. 1

The curve expressing cortisol as a function of circulating CRH from the full model is plotted by solving the exact Eq. 2 and is shown as a dashed line. The curve from the simplified model derived from the approximate Eq. 3 is shown as a solid line. Henceforth we will use the latter as it allows for an analytical solution, while preserving the scale and salient features of the full model . In both cases parameters from Table 1 are used.

We now explore approximations to the c dynamics in the second of Eqs. 1 which couples and c. The analytical work in Ref.[55], [56] reveals that the key to bistability is the sigmoid, “S” shape of the c-nullcline. Since the cubic is one of the simplest forms to yield a sigmoid shape, we tailor an ad hoc cubic function as a proxy for the c-nullcline. To do this in a consistent manner, certain conditions must be met. We first consider the case of a fixed input, . As can be seen in [55], [56] increasing will shift the c-nullcline to the left in space; also the c-loci of the turning points of the c-nullcline are mostly insensitive to . The c-nullcline must also pass through the origin so that if no CRH is synthesized, no CRH can be released. Similarly, when and , indicating that the release of CRH increases when stored, synthesized CRH is present. We also impose that when and so that when there is no stored CRH present, no CRH can be released and its concentration c will decrease due to degradation. We thus posit where is a cubic in c with and for . The proportionality constant should be positive so that values will yield . Furthermore, to preserve positivity along the c-nullcline, , we must also impose that ; positive values of also guarantee that if then , as discussed above. Typical values of c are one order of magnitude larger than those of [55], [56]; as a result must be much smaller than , so that the nullcline equation , which includes a cubic in c, will yield reasonable values of . Under the assumption , we finally writeso that increases in shift the nullcline to the left in space as imposed above. The quantities represent the turning points of the cubic defined by , so that . Without loss of generality we assume . Since must be positive at the turning points, we must also impose and . The form of the right-hand-side of Eq. 4 implies that the c-nullcline follows a sigmoidal “S” shape where the c-coordinate of the turning points are independent of as desired.

We can arrive at a similar expression for Eq. 4 through a different route. Performing a Taylor expansion in of the right-hand side of the c dynamics in Eqs. 1, the term yields to first order, whereas the expansion in of results in to second order so that setting yields a cubic term with a negative coefficient. Upon including a time-dependent form for , the substitution will shift the c-nullcline accordingly. The full model can thus be re-written aswhere is the approximated form given implicitly in Eq. 3, and where time-dependent inputs can be easily incorporated in lieu of . To guarantee bistability, the and the c nullclines must allow for multiple intersections. Let us denote the coordinates of the turning points of the cubic on the c-nullcline as and , where and with . Since c decreases along the slow nullcline as increases, multiple crossings of the two nullclines will arise if lies to the right of the slow, nullcline, and lies to its left. Mathematically this is translated into the following bistability condition

As discussed above, c is typically one order of magnitude larger than [55], [56]. We thus further rescale by a factor so that the two CRH components are comparable in size and define

Note that the rescaling will also affect the time scale t, as seen from Eq. 11. For simplicity, we drop the prime notation from and introduce . Since the dynamics of unfolds on the longer time scale compared to that of c, we also assume so that . Our reduced model is now complete. The ODE system written in terms of is given bywhere is the real, positive solution toand are the rescaled values according to Eq. 10 and subject to the following constraints

The form of our reduced model in Eqs. (12), (13) under the constraints given by Eqs. (15), (16) is reminiscent of that of the FitzHugh-Nagumo model for neuron spiking [64], which exhibits rich dynamics in response to external excitations such as graded responses, the appearance and disappearance of limit cycles, and large excursions in phase space depending on the amplitude of the external input. In Fig. 2 we show the area in parameter space that yields bistability (i.e. two stable equilibria) obtained by taking into account all constraints, Eqs. (15), (16) and the assumption that . For given values selected from this region, the corresponding values can be calculated through Eqs. (15), (16), leading to and coordinates for the turning points of the cubic term in the fast y-nullcline.

Fig. 2

Shaded in gray is the parameter region in space that yields bistability for . The constraints that lead to this region are detailed in the text. The three blue lines are and . The dotted curve is implicitly defined via ; the dashed one via as in Eqs. (15), (16) respectively, with . All other parameters are as listed in Table 1. We set as denoted by the red dot. These values fall in the bistable, gray region and allow for good qualitative agreement between Eqs. (12), (13), and the full model in Eqs. 1.

Once all the physiological parameters are chosen, transitions between the two stable equilibria can be induced by perturbations, even transient ones, to the stimulus . In principle, one can consider a more complex description of the system where instead of bistability, the nullclines intersect three or more times. Two of these equilibria would then represent healthy and diseased conditions, the others would be interpreted as prenosological states. Mental disorders in combatants for example can be preceded by milder neurotic conditions [65]. In this case, changes to could induce direct transitions between healthy and diseased states, but could also modulate a first passage from the healthy to the prenosological state, and then from the prenosological to the diseased state. Another scenario is that of a tailored that once the prenosological state is reached would lead to a reverse transition back to the healthy one, representing early intervention. For simplicity we only consider bistability in the remainder of this paper.

Analysis of the dynamical model without circadian drive

The steady-state solutions of the reduced model Eqs. (12), (13) are found by setting their right-hand sides to zero so that

If the turning point coordinates are chosen as in Fig. 2, two distinct solutions arise and or , labelled according to . In Fig. 3 we compare the nullclines derived from Eqs. (12), (13) with those arising from the full model in Eqs. 1, showing that qualitative features are preserved. Two different ways for the equilibrium values to relate to each other are shown in Fig. 4. To the left, the stable points are characterized by and ; to the right and . When identifying healthy and diseased states with either of the two equilibria, the relation between will be important since, as discussed in the Introduction, stored CRH (a decreasing function of cortisol) and circulating CRH may be both over- or under-expressed, or one may be deficient while the other is produced in excess. The left-hand panel represents the case where abnormalities are marked by high stored CRH and low circulating CRH (or vice versa); the right-hand panel represents the case where abnormalities manifest via both high (or both low) stored and circulating CRH. We will mostly focus on the left-hand side representation but our analysis and conclusions are similarly applicable to the right-hand panel.

Fig. 3

Fast CRH release (solid, red) and slow CRH synthesis (solid, blue) nullclines of the fundamental HPA axis model obtained by setting and respectively in Eqs. (12), (13). For comparison we also plot the corresponding c-nullcline (dashed, red) and -nullcline (dashed, blue) from Eqs. 1 as a function of and . Although the details of the nullclines from the model (Eqs. (12), (13)) and from the full model (Eqs. 1) differ, the main features persist. All parameters used are as in Table 1.

Fig. 4

Fast (red) and slow (blue) nullclines for two schematic realizations of the reduced model in Eqs. (12), (13). (a) The slow nullcline defines a negative slope resulting in the two nullclines intersecting at equilibrium points and , where and . (b). The slow nullcline defines a positive slope corresponding to equilibrium points and , where and . In panel (a) we show the intersections , in panel (b) we show the turning points and. .

We can further analyze Eqs. (12), (13) by noting that the dynamics contains two time scales: t, over which y evolves, and over which x evolves. To make analytical progress we use asymptotic expansion methods [66] by first calculating solutions in these two time scales respectively and then matching the solutions in the intermediate scale where are valid. We begin with the t time scale, and pose

Upon inserting Eqs. (19), (20) into Eqs. (12), (13), we find to leading order in t,which imply the fast dynamics will occur along the vertical, y axis while remains constant. The slower motion instead will arise from assuming Eqs. (12), (13) evolve over the time scale . We thus poseand use the chain rule to derive

By matching orders of we find

Eqs. (27), (28) imply that over long time scales the trajectory moves along the slow nullcline, and that are related via Eq. 28. The two sets of solutions, at the short and long time scales must coincide at intermediate times when and +. We thus impose

Solving Eqs. (21), (22) leads to , while can be found by factoring the right-hand side of Eq. 22 once has been inserted, and solving via separation of variables. On the other hand, solving Eqs. (27), (28) requires specifying the exact form of which can be in principle found through Eqs. 2 or 3. However, since all functions involved are analytic and since upon inspection of Eq. 28 is a function of , we can linearize as a function of and write , where . Under this assumption and using Eq. 29 we find

We can finally estimate by inserting Eq. 31 into Eq. 28 and by solving the resulting cubic equation for . At steady state, two stable solutions emerge for corresponding to . From Eq. 31 we can also estimate the typical time scale to reach steady state as , or . The fully linearized problem, where we set directly into Eqs. (12), (13), allows us to write the two nullclines as a vertical line intersecting a cubic, under fast y and slow x dynamics. We consider this problem in Appendix 6.

Reduced model with circadian drive

Eqs. (12), (13) describe the fundamental dynamics of the full model without circadian driving. One important physiological feature we now add is the circadian cycle, as modulated by the SCN in the hypothalamus. This rhythm is manifest as a small, periodic variation in the basal input . As a result, the full model presented in Eqs. 1 must be rewritten with a time-dependent term which we model as

Here, is the amplitude of the circadian perturbation, and its frequency, defining the period . Since we assume the time-dependent term is small compared to the basal term, , the analysis performed in Section 2 remains valid. We can thus study our reduced problem, Eqs. (12), (13), by replacing with as expressed in Eq. 32 and with the rescaling shown in Eq. 11. Under this non-dimensionalization, , and the rescaled frequency is , corresponding to a period with given in Table 1. Henceforth we consider the rescaled version of Eq. 32, where is as in Eq. 11 and drop the prime notation.

Before numerically analyzing the circadian version of Eqs. (12), (13) driven by Eq. 32, we present some analytical approximations. Since both are assumed to be small, our analysis will depend on how the two relate to each other. If , the circadian rhythm defines the smallest perturbation and solutions to the model without the circadian rhythm presented in Section 2.1 (Eqs. (12), (13) with ) are valid zero-th order approximations to the circadian problem. If instead the analysis presented in Section 2.1 is no longer valid and the contribution of CRH hormonal storage and of the circadian rhythm must be jointly considered. We consider the two cases below.

Small circadian drive limit

We first consider the case, and expand solutions to the circadian problem in Eqs. (12), (13) in different orders of as followswhere are solutions to the reduced non-circadian model, that solve Eqs. (12), (13) with . These solutions equilibrate towards the steady-state values presented in Eqs. (17), (18), , under the bistable conditions shown in Fig. 2. Upon inserting Eqs. (33), (34) into the circadian model with given by Eq. 32 we find

To further simplify our analysis, without loss of generality we also assume that ) in Eqs. (35), (36) are set as one of the two stable, fixed-points under so that Eqs. (35), (36) become an inhomogeneous linear ODE system with respect to . We can now rewrite Eqs. (35), (36) as a matrix equationwhere , , and

Using standard methods to solve linear ODE systems driven by a periodic term we findwhere and are the eigenvalues of , corresponding to eigenvectors and , respectively. We now compare and contrast the analytical results in Eqs. (33), (34) truncated at first order in and with given in Eqs. (39), (40) with numerical evaluation for the two steady-states. For concreteness we specify and in addition to the other parameters given in Table 1. The above choices yield which fall in the bistability region shown in Fig. 2, and lead to . These choices ensure all constraints listed in Sect. 2 are met and that the reduced model is as close as possible to the original full model in Eqs. 1. Unless otherwise noted, the above parameters will be fixed at the above values for the remainder of this work.

In Fig. 5, we show numerical results for . In the top row we consider perturbations around the fixed-point . Similarly, those around the fixed point are shown in the bottom row. For , analytical results from Eqs. (33), (34) agree well with numerical ones for both fixed-points, as increases discrepancies between the two become more pronounced. This is to be expected, as the current analytical results are valid only insofar as . Not shown in Fig. 5 are numerical and analytical results for values of which are also in good agreement. Note that the analytical approximations define curves that are centered about both fixed points . This is because the first-order approximations in Eqs. (39), (40) contain superpositions of oscillatory terms with the same frequency , so that the dynamics in space is symmetric about the central, fixed point.

Fig. 5

Periodic orbits in space driven by a circadian rhythm of amplitude , from left to right. Analytical estimates are derived from Eqs. (33), (34); numerical results are obtained from Eqs. (12), (13). Parameters are listed in Sections 3.1, 3.2 and Table 1. For these values . Analytical estimates closely match numerical results for , as detailed in the text. Top row: periodic orbits for perturbations around the fixed point . Bottom row: periodic orbits for perturbations around the fixed point .

Circadian drive comparable to stored CRH dynamics

We now assume and pose . Since the and contributions in Eqs. (12), (13) and in Eq. 32 are of the same order of magnitude, we must derive our analytical results independently of the results found in Section 2.1. We again expand with respect to as follows

We expect the zero-th order solution, in Eqs. (41), (42) to be the same as the fast solution to the non-circadian problem in Eqs. (21), (22) given that now all perturbations arise to order . We also include second order terms because, as we shall see below, the first order term is constant. Upon inserting Eqs. (41), (42) into Eqs. (12), (13), and using Eq. 32 we find the following sets of identities stemming from the three orders of

Solutions to Eqs. (43), (44), (45), (46), (47), (48) will depend on the initial conditions. As expected, Eqs. (43), (44) are the same as Eqs. (21), (22), so the zeroth order dynamics is the same as in the non-circadian case. If we also assume that are the time-independent equilibrium solutions to the non-circadian model, we also find that the right-hand side of Eq. 45 vanishes, leading to a constant value of . We can now solve Eq. 46 using the above assumptions. We findwhere

Note that , since and and . Eq. 49 depends on the constant value for and on the initial condition . If we assume , the first order dynamics along the y axis is given by

We can also impose and assume are related viato find a purely oscillatory solution

To find the dynamics for and we assume for simplicity and set as in Eq. 52 leading to Eq. 53 so thatresulting in

The value of is arbitrary and can be set to zero; we can also choose so that the higher order solution remains oscillatory. This is shown in Appendix 7, where we also derive , the solution to Eq. 48, which follows from tedious but straightforward computations.

In Fig. 6 we plot as derived from Eqs. (41), (42), truncated to second order. For simplicity, we focus on purely periodic solutions, by setting the proper initial conditions to all the relevant orders of . Specifically we use as given by Eq. 53 and and as given in Eqs. (88), (86) respectively in Appendix 7. As done in Fig. 5, we consider both steady-state values and compare analytical and numerical evaluations for , and with . The results shown in Fig. 6 reveal good agreement between analytical and numerical curves. Upon comparing Fig. 5, Fig. 6 however we find that when , Eqs. (33), (34) are a closer approximation to the numerical results than Eqs. (41), (42), but the reverse is true for , as can be expected since we are here considering the limit. The analytical approximations define curves that are not centered about the fixed points , to the contrary of what observed in Fig. 5. This is because the second-order approximations to Eqs. (47), (48), shown in Eq. 55 and in Eq. 85 of Appendix 7, superimpose constants and oscillatory terms of frequency to the first order approximations centered about the fixed point and containing oscillations of frequency . These second-order terms break the symmetry about the fixed point, so that is no longer central as and in Eqs. (41), (42) oscillate. More details are shown in Appendix 8. We investigate the behavior of the system for larger values of in Appendix 9. As can be seen, for period-doubling, subharmonics, and chaotic behavior emerge. While mathematically interesting, this limit does not allow to clearly distinguish between healthy and diseased states. In the remainder of this work we thus study stress-induced transitions to the HPA axis by keeping the amplitude of the circadian rhythm .

Fig. 6

Periodic orbits in space driven by a circadian rhythm of amplitude , from left to right. Analytical estimates are derived from Eqs. (41), (42); numerical results are obtained from Eqs. (12), (13). Parameters are listed in Sections 3.1, 3.2, Table 1 and Fig. 5. Analytical estimates closely match numerical results for , as detailed in the text. Top row: periodic orbits for perturbations around . Bottom row: periodic orbits for perturbations around . Note the better agreement between analytic and numerical curves for than Fig. .6.

Inducing transitions between steady states

We have hitherto assumed that the two stable states in the absence of the circadian drive, for , represent a healthy and a diseased state. For concreteness and without loss of generality we identify with the healthy state, and with the diseased one, where and . This choice is simply made for illustrative purposes, since as discussed in the Introduction, over- or under-expression of any of the hormones regulated by the HPA axis induce stress-related pathologies. As shown in Sect. 3 including the circadian drive yields limit cycles around both for . Since at steady state limit cycles will orbit around them, we use these points as markers for healthy and/or diseased conditions even in the presence of the circadian drive. Finally, we limit our analysis to the regime, since larger values of may lead to chaotic behavior.

We now study how the system transitions between the two limit cycles, or equilibrium values, in response to external perturbations. Within the diseased/healthy context introduced above, transitions from to (or from and to the limit cycles orbiting around them) represent the onset of disease, and the opposite progression, healing. These transitions may arise through parameter changes, which we associate to physical injury or surgical intervention, or in response to external input such as psychological trauma, cognitive behavioral or exposure therapy. We only consider the latter scenario and study how an external input superimposed to the circadian drive may or may not induce transitions between the two limit cycles. Thus, the basal value in Eq. 13 is replaced by . While several shapes are possible we only consider illustrative examples where a simple pulse of fixed amplitude , duration , is applied at a given phase (modulo ) of the circadian rhythm. The complete form of includes basal, circadian, and reactive inputs to the HPA axis. In previous work, where the circadian drive was not included, the amplitude and duration of the external input were shown to greatly influence transitions between steady states [55], [56]. Here, the acute stressors may evoke strong responses, even surpassing the rhythmic ones; however the amplitude, duration and timing of the circadian drive contribute to the emergence of transitions in non-trivial ways.

Transitions between healthy and diseased states

In Figs. 7 to 11 we show trajectories under piece-wise constant perturbations of magnitude lasting for assuming the initial condition is along the normal or diseased limit cycles for , or, if , at the normal or diseased fixed point. Although in our computations is measured in non-dimensional units, according to the scaling provided in Eq. 11 and in Ref. [55], for concreteness, we describe following figures using dimensional values in the captions. Superimposing an external input to the circadian drive changes the shape of the fast y nullcline, but not of the slow x one. Following [55], we adopt the convention , since the majority of neural circuits that project to the PVN are excitatory. This setting also avoids the emergence of negative values in the driving stimulus, since when . As a result, the fast y nullcline is compressed or stretched horizontally when the circadian and/or external drive are added to the basal . If , the smallest compression factor arises when and the largest compression factor occurs when . If the nullcline is stretched by a factor when and similarly compressed by when . The fast nullcline oscillates between the bounds provided by the two limits above. Whether transitions arise or not depends on the intricate interplay between the external input and the oscillating nullclines. In certain regimes transitions will arise, in others, they will not. Unless otherwise specified, parameters for the figures shown in this subsection and the next are as listed in Sect. 3.1, 3.2 and Table 1; illustration conventions are as described in Fig. 7, Fig. 8.

Fig. 7

Dynamics in phase space for (upper left) , (upper right) , (lower left) , (lower right) , when is active from to hours and the system is initiated at the healthy state. All other parameters are as listed in Sect. 3.1, 3.2, and Table 1. For , the red solid curve is the fast nullcline when , the red dashed curve is the fast nullcline when , and the blue curve is the slow nullcline. The black curve is the trajectory. The three red solid curves that arise for are the fast nullclines when and the circadian drive defines a phase (modulo 2) from right to left, respectively. The fast nullcline oscillates between these three curves. The red dashed curves are the fast nullclines when . We also use three values of the circadian phase , however the corresponding fast nullclines are not sufficiently resolved and appear as a single dashed line. The blue curve is the slow nullcline. When , the external input is sufficient to induce a transition from the healthy state to the diseased one. Including a circadian drive with intermediate hinders the transition and the system remains in the healthy limit cycle. Further increases to restores the transition to the diseased state.

Fig. 11

Dynamics in phase-space when an external input lasting for variable is applied in the presence of the circadian rhythm. In the left panel is applied to the healthy state at time () for  min (solid black curve), and for  min (dash-dot black curve). The amplitude of the circadian rhythm is set at . When  min the external stressor induces a transition to the diseased state, when  min no transition arises. Similarly in the right panel is applied for a system initially at the diseased state, at time  min () for  min (solid black curve) and  min (dash-dot black curve) for . In the first case the trajectory returns to the diseased state, in the second case, the pulse allows the system to transition to the healthy state.

Fig. 8

Dynamics in phase-space for (upper left) and (upper right) when is active from to = 12 h, and for (lower left) and (lower right) when is active for hours, from  min () to  min. The initial condition is always the diseased state. The first row shows the circadian-induced, diseased-to-healthy transition. The second row shows the circadian drive impeding the same transition. In the left panels, the solid red curves are the fast nullcline in the basal state, . The dashed red curves are the fast nullcline when and . The blue curve is the slow nullcline. In the right panels, the three red solid curves to the outermost right are the fast nullcline when and the circadian drive defines a phase , from right to left, respectively. For , the fast nullcline oscillates between them. The three dashed curves to the outermost left are the fast nullcline when and similarly from right to left. For the fast nullcline oscillates between them. Due to the numerical values of some of the fast nullclines superimpose. From right to left, the driving terms of the fast nullclines are (superimposed), (superimposed), . We plot the healthy limit cycle as a dotted black curve in the lower right panel for reference.

We begin in Fig. 7 by showcasing the dynamics of the system when an external stressor is applied for hours with and without a circadian rhythm of varying amplitude . The initial state is the healthy state at . The four panels show that increasing the magnitude of affects the transitions from the healthy to the diseased state in a non-trivial way. As can be seen, in the absence of the circadian drive, induces a transition from to , a modest amplitude hinders the transition and the system remains in the healthy state, but larger values restore the transition to the diseased one. In the top row of Fig. 8 we initiate the system at the diseased state and apply a pulse for hours with and without a circadian rhythm of amplitude . For the system persists at the diseased state even after the pulse is applied, whereas for induces a transition to the healthy state. In the bottom row of Fig. 7 the opposite outcome is observed. Here, the system is also originally initiated in the diseased state and is applied for hours. In the absence of the circadian drive, the external input induces a transition to the healthy state, however a circadian rhythm with amplitude hinders the same transition, with unable to dislodge the system from the diseased state. These two examples show that the circadian drive may strongly influence how external pulses, representing therapeutic intervention, affect the system. Subjects may thus respond differently to cognitive behavior treatment depending on how sensitive they are to the circadian rhythm.

To understand the role of the magnitude , in Fig. 9 we consider the non-circadian case, where the external pulse is active from to hours and where the magnitude is increased from to , while keeping all other parameters fixed. In both cases, the external perturbation dislodges the system from the diseased state towards the perturbed equilibrium given by the intersection of the slow nullcline with the perturbed fast nullcline. Since is large enough, the system will settle in this new equilibrium, marked by the red cross in Fig. 9. However, once is terminated, the system must return to either of the original steady states. For , the perturbed equilibrium falls into the basin of attraction of the original healthy equilibrium, so the trajectory will settle into the healthy steady state and a transition is recorded, as shown in the left panel of Fig. 9. The value instead leads the system further away from the original nullclines so that here the new perturbed equilibrium, marked by the red cross, falls into the basin of attraction of the original diseased state. As a result, after a large excursion in phase space, the trajectory returns towards the diseased state. This is a slightly counter-intuitive finding, since one may expect that larger pulses may better facilitate transitions to the healthy state.

Fig. 9

Dynamics in phase-space for (left) and (right) active from to hours. In the left panel, causes a transition from the initial diseased state to the perturbed equilibrium, which falls into the basin of attraction of the non-perturbed healthy equilibrium. Once is terminated, the system transitions to this healthy state. A larger disrupts the dynamics further and the new perturbed equilibrium falls into the basin of attraction of the non-perturbed diseased state. After a long excursion, the termination of returns the system to the diseased state.

In Fig. 10 we show that another determinant of the transition between states is given by the starting time of relative to the phase of the circadian rhythm. In the left panel we apply two external inputs lasting hours. In the first case is applied from to  min. The trajectory in the black solid line, is perturbed but stays in the healthy state. In the second case, although similarly and hours, the external input is delayed, starting at  min and terminating at  min. As can be seen, the transition to the diseased state does occur. The same dynamics is observed in the right panel where the system is initiated at the diseased state and two external inputs lasting hours are applied. When is applied from to  min, a transition towards the healthy state arises, whereas when the external stressor is delayed and initiated at  min and terminated at  min, the transition does not occur and the trajectory returns to the diseased state.

Fig. 10

Dynamics in phase-space where an external input of magnitude lasting is applied at different phases relative to the circadian rhythm. In the left panel is applied for a system initially at the healthy state, for  min, with start times of (, solid black curve) and  min (, dash-dot black curve). The amplitude of the circadian rhythm is set at . When is initiated at the trajectory remains at the healthy state, when is initiated at  min the trajectory transitions to the diseased state. Similar dynamics are shown in the right panel when is applied for a system initially at the diseased state, for  min, with start times of (, solid black curve) and  min (, dash-dot black curve). When is initiated at the trajectory transitions to the healthy state, when is initiated at  min the trajectory returns to the diseased state.

The effects of duration of the external input are shown in Fig. 11 where is applied at but for different durations . In the left panel, the system starts in the healthy state and . When an external stressor is applied for hours (solid black curve), a transition to the diseased state is observed, however when is applied for hours (dash-dot black curve), the system remains at the healthy state. Similarly, in the right panel the system starts in the diseased state and . When an external pulse is applied for hours (solid black curve), the system does not transition to the healthy state, however when is applied for hours (dash-dot black curve), a transition to the healthy state is observed. Upon further increasing to hours, the system will remain in the diseased state. We do not plot this case, but the trajectory follows very closely the one observed for hours. These results show that the amplitude of the circadian drive , the magnitude , onset phase , and duration of the external input affect transitions between states in subtle ways. Finally, Fig. 9, Fig. 10, Fig. 11 show that under stressed conditions, trajectories may access perturbed equilibria that persist as long as the stress remains active. These perturbed equilibria can be considered meta-stable states that vanish once the external stressor/input is removed. The system will then return to its original healthy or diseased, non-stressed state or transition to the opposite one.

Transition diagrams between healthy and diseased states

In this section we study the interplay between and , by plotting transition diagrams in space for two possible initial conditions (healthy and diseased) and several values of . These diagrams will indicate whether a transition has occurred or not after a stressor or pulse characterized by applied at phase is terminated.

In Fig. 12 we start the system in the healthy (left panel) and the diseased (right panel) states and show the final equilibrium configuration in the absence of the circadian rhythm (). As can be seen in the left panel, transitions from the healthy to the diseased state arise for sufficiently large external stressors of sufficient duration . Interestingly, in the right-hand panel, long lived pulses are associated with transitions from the diseased to the healthy state only if is of intermediate value, but not too small or too large. Why would large pulses that last long enough not induce transitions? The answer is the asymmetry induced by selecting . This choice always compresses the fast nullcline and its corresponding limit cycle to the left. When starting in the diseased state, small values of are not large enough to sufficiently perturb the system, so no transitions are observed, as can be expected. Intermediate values may induce transitions, however if is too large, the compression to the left may cause the intersection between the fast and slow nullclines corresponding to the diseased state to vanish. The trajectory must leave the initial (diseased) equilibrium and start evolving towards the perturbed one arising under . For short durations , upon termination of the external pulse, the system may still be en route to the perturbed equilibrium and/or be in the basin of attraction of the non-perturbed healthy state and a transition will be observed. However, when is sufficiently large, the trajectory will be in the proximity, or at, the perturbed equilibrium while is still active. Once the pulse is terminated, if is large enough, the perturbed equilibrium may be sufficiently distant from the non-perturbed nullclines and fall into the basin of attraction of the original, diseased equilibrium. In this case, there will be a large excursion in space but the system will eventually return to the original diseased state. Representative dynamics are shown in Fig. 9 where an intermediate value of for sufficiently large induces a transition, but larger values of do not. These results show that large, long lived external inputs will always result in the system stabilizing at the diseased state, regardless of initial conditions, indicating that therapeutic intervention should be applied judiciously.

Fig. 12

Final equilibrium states in space in the absence of the circadian rhythm (), starting from the healthy (left) and diseased state (right). Durations are set as , with , magnitudes are set at , with , defining 1600 combinations, corresponding to 1600 rectangles centered at . Colors represent the final equilibrium state after cessation of the perturbation: yellow ones indicate a final healthy equilibrium at , purple ones indicate a final diseased equilibrium at . Note that large and long perturbations always lead to the diseased state, regardless of the initial configuration.

We present similar results in Figs. 13 to 18 where the amplitude of the circadian rhythm is set at , and where the system is initiated at the healthy state in Fig. 13, Fig. 14, Fig. 15, and at the diseased state in Fig. 16, Fig. 17, Fig. 18. Plotting conventions are the same as in Fig. 12; however, here we also consider the phase at which is superimposed on the initial, healthy or diseased, limit cycle relative to the circadian rhythm. Specifically, we set , with , resulting in eight panels for each choice of . see Fig. 17.

Fig. 13

Starting from a healthy-state limit cycle, we show the final configurations in space for and eight values of the start time of the external stress relative to the circadian drive .

Fig. 18

Starting from a diseased-state limit cycle, we show the final configurations in space for and eight values of the start time of the external stress relative to the circadian drive .

Fig. 14

Starting from a healthy-state limit cycle, we show the final configurations in space for and eight values of the start time of the external stress relative to the circadian drive .

Fig. 15

Starting from a healthy-state limit cycle, we show the final configurations in space for and eight values of the start time of the external stress relative to the circadian drive .

Fig. 16

Starting from a diseased-state limit cycle, we show the final configurations in space for and eight values of the start time of the external stress relative to the circadian drive .

Fig. 17

Starting from a diseased-state limit cycle, we show the final configurations in space for and eight values of the start time of the external stress relative to the circadian drive .

In Fig. 13, is applied when and the system is in a limit cycle about the healthy state. For relatively low and/or , no transitions to the diseased state occur and the system remains in the healthy state. However, external stressors with longer duration and larger magnitude do not necessarily imply a higher likelihood of transitioning to the diseased state, as observed in the non-circadian case. Rather, regions in space where transitions to the diseased state occur are separated from those where no transitions occur by an undulating parameter boundary exhibiting multiple minima in , a consequence of the pulsatile circadian rhythm imposed on the HPA dynamics. For a null onset phase of the external stressor relative to the circadian rhythm, , the nadir of the separatrix emerges at . As increases from to , this minimum shifts to shorter following an approximate trend. At , an additional nadir emerges at , which follows a trend as further increases. The time at which the stressor ends defines a phase relative to the circadian rhythm given by where we add the onset phase and the phase defined by the duration of the external stressor . For all recorded minima in Fig. 13, (modulo 2), indicating that the lowest that can be applied to induce a transition from the healthy to the diseased limit cycle should last long enough, and be terminated when the circadian rhythm is at .

Qualitatively, low values of and/or do not allow the system, initially at the healthy limit cycle, to reach the diseased one, so transitions are unlikely. Once and/or are sufficiently large however, the trajectory will reach the perturbed limit cycle and remain anchored to it. Further increasing and/or will still result in the trajectory oscillating about the perturbed healthy limit cycle so that beyond a certain threshold the magnitude of and/or become irrelevant. Thus, upon terminating the external input, whether the trajectory returns to the healthy limit cycle or is able to escape towards, and finally attracted into, the diseased one, is highly sensitive to the end phase , which determines the position of the oscillating fast nullcline, and much less to and/or . The optimal end phase for the trajectory to remain at the diseased state is the “neutral” configuration .

These features persist for , as shown in Fig. 14 where all minima in space are marked by an end phase (modulo 2), just as for . The region where the transition is observed (depicted in purple) however is larger than for . Here, the larger oscillations induced by the circadian drive allow lower values of and/or to induce transitions. Finally, in Fig. 15 where and circadian oscillations are even wider, transitions occur in most of phase space, except for very low values of .

In Fig. 16, Fig. 17, Fig. 18 we initiate the system at the diseased limit cycle around and set the amplitude of the circadian rhythm to , respectively. In Fig. 16 a small parameter region with low and a broad range of is observed where trajectories remain at the diseased state. For slightly larger values of , pulses with intermediate durations will induce a transition to the healthy limit cycle; however, larger may not be conducive to transitions. The dynamics here mirrors the non-circadian case, where large values of may cause the diseased steady state to vanish. Finally note that the boundaries between regions in space where transitions do and do not emerge shares some similarities to the ones observed in Figs. 13 to 15.

Discussion and conclusions

Building on previous work, we incorporated a circadian cycle into a dynamical systems model of the HPA axis. We start by considering the dynamics and interplay of five main quantities involved in the HPA axis response to external perturbations (stored CRH, circulating CRH, ACTH, cortisol and glucocorticoid receptor concentrations), and introduce several simplifications to build a simpler two-dimensional reduced model. We first review the non-circadian limit and neglect the known delay in ACTH activation of adrenal gland secretion of cortisol. This delay of approximately fifteen minutes allows for the emergence of ultradian (hourly) oscillations. Since we are interested in the effects of the circadian (diurnal) cycle, which unfolds over a longer time frame than hourly oscillations, we do not include delays, and refer the interested reader to previous work where its impacts are discussed [55], [56]. We also note that the dynamics of ACTH, cortisol and glucocorticoid receptor concentrations, evolve on much shorter time scales than stored CRH, allowing us to consider the steady state values of these quantities, effectively projecting the original system of five equations to a set of two coupled ones for stored and circulating CRH. Of these, the first evolves on a longer time scale, of the order of hours, the other on a shorter one on the order of minutes. Once the equations for stored and circulating CRH are solved, values for the other quantities (ACTH, cortisol and glucocorticoid receptors) can be derived by substitution in the respective steady state expressions. Using order of magnitude estimates for the various terms involved we finally express the right-hand side of the dynamics of circulating CRH via a cubic expression. This form bears no physiological relationship to the actual evolution of circulating CRH, however it is a very good ad hoc substitute in that it allows for a thorough mathematical analysis while preserving the main features of the system, namely the emergence of bistability, i.e. of two steady states, one marked by low values of circulating CRH, and the other by higher values. These stable states emerge as intersections of a fast nullcline, when circulating CRH reaches equilibrium, and a slow nullcline, when stored CRH reaches equilibrium. We interpret them as diseased and healthy states, respectively. After analyzing the two coupled, self-contained equations using a cubic approximation (without the circadian rhythm, Eqs. (12), (13)), we included the circadian drive and considered several possibilities for the amplitude of the circadian component relative to the evolution of the slow nullcline. We show that low amplitudes of the circadian rhythm turn the two fixed points from the non-circadian system into limit cycles about them, allowing us to identify healthy and diseased limit cycles. We also obtain analytical approximations for them. For larger values of the circadian amplitude, the two limit cycles merge into one and chaotic behavior is observed, reminiscent of Duffing oscillator dynamics. Finally, we include an external, constant pulse input (box function) of amplitude and duration , and investigate the trajectories of the system and whether transitions from healthy to diseased states or vice versa can occur. We find that whether transitions arise or not depends in a complex way on and the phase at which the external input is applied or terminated relative to the circadian rhythm. Interestingly, we find that an important determinant for the existence of transitions depends on the phase at which the external stressor is terminated relative to the circadian rhythm. We also find that large perturbations applied for sufficiently long times on the diseased state may greatly disrupt the system leading to large excursions in the dynamics of stored and circulating CRH. This disruption may however be temporary and the system will eventually return to the original diseased state. In the context of therapeutic intervention, this result may signify that large, long lasting external inputs may not be as effective as more moderate ones, applied for shorter times. We also find that due to the nullcline structure of the system, remaining in/transitioning to the diseased state is more likely for large circadian amplitudes than for smaller ones.

In this work we have assumed that the diseased state corresponds to larger stored CRH and lower circulating CRH, exemplified by the steady state (or the limit cycle) about . This configuration corresponds to low circulating levels of cortisol since and, at steady state, is a decreasing function of cortisol. Conversely, the healthy equilibrium or limit cycle at is characterized by lower stored CRH (and higher cortisol) and larger circulating CRH. This is one of many choices, since other pathologies may exhibit different relative hormonal levels between healthy and diseased states. For example, acute manifestations of major depressive disorder (MDD) are marked by high levels of secreted CRH and of cortisol [17], [21], implying that the healthy state can be identified with and the diseased one with . In this case the external stressor applied to the healthy state can be identified as a triggering event that leads to MDD, or, if applied to the diseased state, as a pulse associated to exposure or cognitive behavioral therapy (CBT). Mutatis mutandis, our results would also indicate that transitions to MDD may be more easily induced by lower values of the circadian amplitude , consistent with experimental findings that find blunted amplitude of the circadian rhythm in depressed patients [67]. On the other hand, PTSD is often associated with low levels of cortisol but relatively high values of CRH. This would translate to large values for the diseased state, or equivalently large . In our specific model, the slow nullcline carries a negative slope as shown in the left-hand panel of Fig. 4. In order for the intersection between nullclines to yield a set of larger values in the diseased state, parameters would have to be chosen so that the slow nullcline has a positive slope has a positive slope, as in the right-hand panel of Fig. 4.

The SCN modulates not only the input signal , but the activity of the adrenal gland as well. The latter responds to circadian stimuli via the splanchnic nerves that relay synaptic signals from the central nervous system to the peripheral sympathetic neurons. These innervate, among other organs, the adrenal medulla, allowing input from the SCN to be relayed to the adrenal tissues. Thus, another possible way of including circadian responses is to include oscillations in the sensitivity parameter that regulates cortisol release in the adrenal gland, as stimulated by ACTH. Specifically, further analysis would include diurnal periodicity in the production of cortisol, driven by an ad hoc periodic form so that in Eqs. 1. The two oscillatory responses to the SCN, and , may be out of phase with each other and could lead to interesting “constructive” and “destructive” interference. Diurnal responses in the HPA axis are also known to depend on age, gender, neuroadaptation, exposure to light, jet-lag, and the use of medication [68], [69], [70], [71], [72], [73]. Interindividual variability and/or stochastic fluctuations may also be present. These influences can be incorporated in the analysis by introducing time-dependent parameters, by coupling to environmental stimuli that modulate and/or and introduce possible time-dependent phases, and/or by including random noise to the circadian amplitude . Different circadian clock speeds could also be studied following Aschoff’s rule, whereby the circadian period is shortened in diurnal mammals exposed to bright, constant light [74]. Conversely, individuals suffering from non-24-h sleep-wake disorder, display abnormally long circadian rhythms, with a period of 25 or 26 h [75]. These scenarios lead to circadian inputs of the type (bright, constant light) and (sleep-wake disorder) that could be compared to the standard circadian drive where , and to investigate how shorter or longer periods affect the dynamics, and the system’s response to acute external stressors.

Also of interest is the pineal gland which translates light/dark inputs from the retina into hormonal signals. Specifically, the pineal gland releases melatonin via the activation of beta–adrenergic receptors by norepinephrine to regulate the circadian rhythm. Stress-induced activation of the HPA appears to correlate with an increase in melatonin production [76]; other studies hypothesize that the pineal gland produces CRH-inhibiting factor concomitantly with melatonin, and that melatonin may inhibit production of cortisol in the adrenal gland [77]. Vice-versa, the beta-adregenic receptors of the pineal gland are believed to be sensitive to ACTH concentrations [78]. These findings suggest strong interplay between stress and the circadian system [79], so that a natural extension of our work would be to relate pineal gland and HPA axis dynamics through proper feedback and feedforward equations. Other interesting avenues of research would include analysis of neuroendocrine systems upstream of the HPA axis or more general forms of , in the context of control theory. Similarly, environmental stochasticity could also be explored. Preliminary results where is modeled as a random variable subject to white noise, show that when the noise amplitude is large enough, transitions not previously realized can arise, indicating that external fluctuations may strongly affect the dynamics. Our mechanistic model provides a framework through which therapeutic interventions can be explored, such as cognitive behavioral therapy, often recommended as a first intervention for many psychological disorders.

Appendix – Linearized model

To better study the dynamics of our specific problem we linearize Eq. 12 by expanding as a function of x. This allows us to also formulate a more general system made of two intersecting nullclines: a slow one, given by a vertical line, and a fast one, given by a cubic curve. In this Appendix we study this simplified system as it may serve as a new paradigm for slow-fast bistable dynamical systems. To begin, we write Eqs. (12), (13) aswhere . We also rescale and introduce as follows

We finally arrive at a fundamental formwhere we dropped the prime notation. The slow u, and fast v nullclines respectively are

The formulation in Eqs. (66), (67) allows us to better study the dynamics and to classify different behaviors. Along the fast v nullcline the following holdswhich implies that if , the v nullcline is an increasing function of v, and the intersection with the slow u nullcline will yield only one stable fixed-point. Conversely, for , the fast v nullcline will have two extrema, located at and corresponding to . The and points also represent the turning points of the cubic equation. Note that for , if , then , and that for is negative. If or only one stable fixed-point will arise from the intersection of the two nullclines, whereas if , three intersections exist that yield two stable, fixed-points. We denote these by and , with . Henceforth, we assume to be in the bistable regime, . The fast dynamics occurs along the following trajectoryand thus, initially at least, the system evolves only along the v axis, with u being fixed at its initial condition . By introducing , we can also write the equations that describe the slow dynamics

Since Eq. 71 can be solved analytically, the slow motion along the u axis is

Upon inserting Eq. 73 in Eq. 72 we can implicitly determine the slow dynamics of

Although trajectories usually move quickly along the v direction, towards the fast v nullcline, for some initial conditions the dynamics may evolve in more subtle ways and additional local analyses may be required. One interesting case is if the initial position is within a small neighborhood of the turning points of the fast v nullcline. Without loss of generality we focus on the turning point and assume , with a small perturbation from , and similarly . As per Eqs. (64), (65), the dynamics is such that if the trajectory will escape the initial neighborhood of the turning point and reach the upper stable fixed-point at . Similarly, if the trajectory will reach the lower stable fixed-point at . Determining the dynamics in other cases depends on the amplitude of . To be concrete, we assume that . In this case, if reaches before reaches , the trajectory will reach the upper stable fixed-point at , whereas if the opposite is true, the lower stable fixed-point at will be reached. This is depicted in Fig. 19, Fig. 20, where if the trajectory crosses the blue-dashed boundary it will reach , while if it crosses the red-dashed boundary it will reach. For a given initial condition we thus estimate the time it takes for , and the time it takes for , and compare the two to determine which basin of attraction the initial condition belongs to. We begin by studying the dynamics for the slow variable by writing the solution to Eq. 64 as followsso that at . Note that by writing Eq. 75 yields . We now pose and find the time evolution for . Upon inserting Eq. 75 into Eq. 65 and expanding the cubic around we find that obeys the following

Fig. 19

Detail of the fast nullcline structure in space near the lower turning point . The dynamics of any starting point in the quadrant defined by is determined by time and trajectories take to the horizontal or vertical lines at or , respectively. For the red initial condition at so the trajectory escapes the basin of attraction of the lower equilibrium point at and eventually reaches the higher equilibrium point at . For the blue initial condition at =, so the trajectory remains in the basin of attraction of the lower equilibrium point and settles at .

Fig. 20

Periodic orbits in space driven by a circadian rhythm of amplitude , from left to right. The green curves represent the exact, numerical solutions obtained from Eqs. (12), (13); since , the blue curves represent analytical approximations up to second order in with in Eq. 55 and first order in with given by Eq. 53; the yellow curves represent analytical approximations up to second order in where given by Eq. 85 is further included. The second order curves typically approximate the exact, numerical solutions with higher accuracy.

The linear term in in Eq. 76 does not contribute to the dynamics, since we are expanding around an extremum of the cubic curve. Also note that and that . We also neglect higher order terms. Eq. 76 defines a Riccati equation that can be solved by imposing and deriving a second order differential equation for . After tedious but straightforward algebra we find that can be written as a linear combination of Bessel functions and that obeyswhere and the argument of the and Bessel functions is

Finally, the prefactor H is

The time it takes for the trajectory to reach the vertical line at , from is given by imposing . From Eq. 75 we find

Similarly, the time to reach the horizontal line at , from is given by imposing . From Eq. 77 we find the implicit expression for

Thus, if the initial condition will belong to the basin of attraction of the lower fixed-point at , whereas if the initial condition will equilibrate at the higher fixed-point at . Finally, the dynamics under initial conditions and depend on how compares with , where is given by

If , and , the trajectory will equilibrate at the lower stable fixed-point at . If , the explicit time dependence of the trajectory must be evaluated.

Appendix - Higher order solutions to the circadian problem

In this Appendix we find , the second order solution in to the circadian problem presented in Eqs. (12), (13) as driven by Eq. 32. We thus solve Eq. 48 and, for simplicity, specifically seek periodic solutions. To this end, we set the proper initial conditions so that only oscillatory terms arise from Eq. 48. Upon inserting Eqs. (53), (55) into Eq. 48, with , we obtainwhereand where . Thus,

To obtain a purely oscillatory solution we impose that the prefactor be zero, so that

We also determine the value of that should be used as the proper initial condition to Eq. 55 if one is seeking a pure oscillatory solution to order . We thus derive the time dynamics for

After inserting Eqs. (53), (55) and 85 with into Eq. 87, we set the resulting constant term in the right-hand side to zero, so that upon integration of Eq. 87 there are no linear terms and is purely periodic. After some algebraic manipulations we find

Appendix - Numerical comparisons between first and second order approximations

In this section we show the difference between lower (first) and higher (second) order approximations to the limit cycles in space driven by the circadian rhythm where . In Fig. 20, the exact, numerical solutions obtained from Eqs. (12), (13) are shown in green; since , second order approximations of upto from Eq. 55 and first order approximations of upto from Eq. 53 are depicted in blue, and second order approximations both in and , with from Eq. 85 further included are plotted in yellow. Second order curves typically approximate the exact, numerical solutions with higher accuracy.

Appendix - Stored CRH dynamics as the smallest perturbation

We now explore numerical solutions for values of where period-doubling, subharmonics, and chaotic behavior emerge. Unless specified, all parameters are chosen as in Section 3.1, 3.2 and Table 1 with initial conditions set at . Henceforth all of our numerical integration will be performed using ode45 in MATLAB®. Before illustrating our results, we note that the circadian version of Eqs. (12), (13), where , is reminiscent of the well known Duffing oscillator, characterized by a periodic forcing term and nonlinear elasticity [80]. The Duffing oscillator displays period doubling, other subharmonic responses, and chaos, as the amplitude of the forcing term is varied. These features arise from the nonlinear elasticity, so that a forcing term of period T may yield a response of augmented period , where m is an integer. The case represents period-doubling. Other values of m are subharmonic responses of order , whereas very large values of m are associated to chaos. Our circadian system is more complicated than the Duffing oscillator, however, inertia, nonlinear terms and periodic forcing are similarly present. One key parameter is , the prefactor of the circadian term and the equivalent of the amplitude of the forcing term in the Duffing oscillator. We find that increasing leads to the emergence of period-doubling bifurcations, subharmonics, and chaos in the circadian version of Eqs. (12), (13), just as for the Duffing oscillator.

In Fig. 21, Fig. 22, Fig. 23, Fig. 24 we plot the time series and and the phase portraits as derived from Eqs. (12), (13) with , for increasing values of . All trajectories are initiated at the steady state of the reduced model Eqs. (12), (13) without the circadian term. We let the system run for and then reset the clock, following the dynamics for a further period. Since we are interested in the long term dynamics, we only plot the trajectories during the final period and do not show the initial transient. The sole exceptions are the chaotic cases in Fig. 21, Fig. 22 where the transient is set at . The red dots in Fig. 21, Fig. 22, Fig. 23, Fig. 24 represent the Poincaré section, the discrete set of points in phase space at times . If the response is a simple periodic orbit the Poincaré section is a single point; after a period-doubling bifurcation, the Poincaré section consists of two points, representing doubling of the periodicity. For a subharmonic response of order , the Poincaré section consists of m points, indicating that the period is increased by a factor m with respect to the circadian term. Chaos leads to a richer set of points in the Poincaré section. Chaotic behavior arises for as can be seen in Fig. 21, Fig. 22. In Fig. 23 we find a subharmonic with period for . The limit cycle passes in the vicinity of both and and crosses itself. A slight increase to leads to a subharmonic with period , yet further increases to result in period-doubling with the emergence of limit cycle of period . Finally, in Fig. 24 we observe the limit cycle acquiring simpler periodic forms from to . In Appendix 10 we show chaotic phase portraits for over a much longer time frame, up to . We also tested both lower accuracy (ode23) and higher accuracy (smaller time step-sizes and stricter tolerances for absolute and relative errors) methods to verify that chaotic behavior is an intrinsic feature of the dynamics and not an artifact of our numerical computation.

Fig. 21

Time series and phase portraits from Eqs. (12), (13) driven by for . Other parameters are listed in Sect. 3.1, 3.2 and Table 1. The red dots in the phase portraits indicate times that are integer multiples of . Periods are T, and (). Chaos emerges at .

Fig. 22

Time series and phase portraits from Eqs. (12), (13) driven by for . Other parameters are listed in Sect. 3.1, 3.2 and Table 1. The red dots in the phase portraits indicate times that are integer multiples of . Chaotic trajectories arise for . For , the period is .

Fig. 23

Time series and phase portraits from Eqs. (12), (13) driven by for . Other parameters are listed in Sect.,3.1, 3.2 and Table 1. The red dots in the phase portraits indicate times that are integer multiples of . Periods are (), (), and ().

Fig. 24

Time series and phase portraits from Eqs. (12), (13) driven by for . Other parameters are listed in Sect.,3.1, 3.2 and Table 1. The red dots in the phase portraits indicate times that are integer multiples of . Periods are (), (), and T ().

Appendix – Chaotic phase portraits

In Fig. 25 we show the chaotic phase portraits and the Poincaré maps arising from Eqs. (12), (13) driven by the circadian term for . All other parameters are chosen as in Sect. 3.2, Appendix 9 and Table 1. The initial condition is set at of the reduced problem Eqs. (12), (13) without the circadian drive. As described in Appendix 9 we let the system run for and then reset the clock, following the dynamics for an additional . We do not show the initial transient in any of the panels in Fig. 25 but follow the dynamics over the subsequent . There is no limiting period over the timeframe.

Fig. 25

Phase portraits in space arising from Eqs. (12), (13) driven by for . The initial transient is not shown. Numerical evaluations are carried out up to 10,000T. Each red dot represents times that are integer multiples of indicating chaotic behavior.

Declaration of Competing Interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.