Dynamic Hormone Control of Stress and Fertility

Neuroendocrine axes display a remarkable diversity of dynamic signaling processes relaying information between the brain, endocrine glands, and peripheral target tissues. These dynamic processes include oscillations, elastic responses to perturbations, and plastic long term changes observed from the cellular to the systems level. While small transient dynamic changes can be considered physiological, larger and longer disruptions are common in pathological scenarios involving more than one neuroendocrine axes, suggesting that a robust control of hormone dynamics would require the coordination of multiple neuroendocrine clocks. The idea of apparently different axes being in fact exquisitely intertwined through neuroendocrine signals can be investigated in the regulation of stress and fertility. The stress response and the reproductive cycle are controlled by the Hypothalamic-Pituitary-Adrenal (HPA) axis and the Hypothalamic-Pituitary-Gonadal (HPG) axis, respectively. Despite the evidence surrounding the effects of stress on fertility, as well as of the reproductive cycle on stress hormone dynamics, there is a limited understanding on how perturbations in one neuroendocrine axis propagate to the other. We hypothesize that the links between stress and fertility can be better understood by considering the HPA and HPG axes as coupled systems. In this manuscript, we investigate neuroendocrine rhythms associated to the stress response and reproduction by mathematically modeling the HPA and HPG axes as a network of interlocked oscillators. We postulate a network architecture based on physiological data and use the model to predict responses to stress perturbations under different hormonal contexts: normal physiological, gonadectomy, hormone replacement with estradiol or corticosterone (CORT), and high excess CORT (hiCORT) similar to hypercortisolism in humans. We validate our model predictions against experiments in rodents, and show how the dynamic responses of these endocrine axes are consistent with our postulated network architecture. Importantly, our model also predicts the conditions that ensure robustness of fertility to stress perturbations, and how chronodisruptions in glucocorticoid hormones can affect the reproductive axis' ability to withstand stress. This insight is key to understand how chronodisruption leads to disease, and to design interventions to restore normal rhythmicity and health.

1. Introduction

A robust dynamic interplay between body rhythms is essential to sustain healthy states. This requires the coordination of several regulatory systems spanning multiple levels of organization, from molecular, to cellular, to the whole organism. Network physiology approaches employ analytical tools, such as mathematical modeling to investigate the interactions between organs and their integration into physiological systems. Neuroendocrine axes are the perfect example of such interlocked-regulatory systems controlling body rhythms, with the brain decoding circadian and stress inputs as well as integrating feedback signals from endocrine organs. The hypothalamic-pituitary-adrenal (HPA) axis and the hypothalamic-pituitary-gonadal (HPG) axis are the major neuroendocrine systems underpinning stress and fertility, respectively. These axes control a range of hormonal and neural activity rhythms exhibiting ultradian (<24 h), circadian (~24 h) and infradian (>24 h) periodicity (Walker J. et al., 2010), as well as responses to environmental, biological and behavioral perturbations. For example, the HPA axis uses feedback control to regulate stress responses while sustaining ultradian and circadian glucocorticoid (CORT) rhythms (Walker J. J. et al., 2010; Spiga et al., 2017). On the other hand, the HPG axis controls infradian oscillations of reproductive hormones secreted in response to changes in the ultradian frequency of gonadotropin-releasing hormone (GnRH). GnRH secretion is controlled by a hypothalamic pulse generator (PG) (Voliotis et al., 2019), which is in turn modulated by gonadal hormones (Figure 1A). Mathematical modeling has significantly contributed to our understanding of this rhythmic behavior (Zavala et al., 2019; Clément et al., 2020), as well as the ability of these systems to respond to perturbations (Spiga et al., 2017). For instance, a mathematical model of rhythmic HPA axis activity has been introduced in Walker J. J. et al. (2010). In this model, ultradian CORT pulsatility is generated by a pituitary-adrenal feedback loop. The model predicts that, in contrast to the reproductive axis, in the stress axis the hypothalamus only needs to provide circadian amplitude modulation of ultradian CORT pulses to explain experimental observations (Walker et al., 2012). Regarding the HPG axis, a combination of mathematical modeling and experimental physiology has shown how the hypothalamic kisspeptin neuronal network generates and sustains pulsatile LH secretion (Voliotis et al., 2019). More recently, a generalized integrate and fire model has been postulated as a simple mechanism to generate a range of rhythmic neuroendocrine signals (Churilov et al., 2020).

FIGURE 1

Figure 1. Pictorial representation of the model. (A) Physiological model of the stress and reproductive neuroendocrine axes controlling ultradian, circadian, and infradian hormone oscillations. Includes the KNDy neuronal network controlling the GnRH pulse generator. Adapted from Zavala et al. (2019). (B) Network model of the systems-level cross-regulation between glucocorticoid (CORT) rhythms, the hypothalamic GnRH pulse generator and the estrous cycle, subject to stress and circadian inputs. Includes a mean-field model of the KNDy network from Voliotis et al. (2019).

Most of the evidence on the dynamic interactions between the HPA and HPG axes comes from animal studies (Acevedo-Rodriguez et al., 2018; McCosh et al., 2019). For example, experiments in ovariectomized rats show a reduction in circadian levels of CORT, which is restored to physiological levels with estradiol administration (Seale et al., 2004a,b, 2005a,b). Data from rodents also shows how physiological and psychosocial stressors can temporarily disrupt GnRH pulse generator activity. These stressors range from isolation and restrain, to insulin induced hypoglycemia and high exogenous CORT, with evidence suggesting the involvement of kisspeptin neuron activity (Li X. F. et al., 2004; Luo et al., 2016; Yang et al., 2017; Ayrout et al., 2019; Kreisman et al., 2019). Further studies in macaque have shed light on the sensitivity and resilience of the reproductive axis to stress signals (Herod et al., 2011a,b). Human studies have also highlighted the profound effect of glucocorticoid excess on the menstrual cycle (Ding et al., 1988; Suh et al., 1988; Saketos et al., 1993; Crofford et al., 1999). However, there is still a limited understanding of whether and how the HPA and HPG axes coordinate their hormone rhythms, how perturbations to one axis impact upon the other, what makes their dynamics robust to such perturbations, and in what circumstances chrono-disruptions can lead to disease.

In this manuscript, we investigate the dynamic control of stress and fertility by means of a mathematical model that accounts for the complex interactions between the HPA and HPG axes. First, we postulate that these neuroendocrine systems behave as a network of coupled oscillators that coordinate ultradian, circadian and infradian rhythms, and validate the model predictions against published physiological observations in female rodents. Second, we consider the evidence on stress-induced suppression of GnRH pulse generator activity dependent of estradiol, and use the model to understand the role of estradiol-dependent effects on the HPA axis (Seale et al., 2004b; Phumsatitpong and Moenter, 2018). We also simulate the simultaneous effect of exogenous estradiol and glucocorticoids on the dynamics of the GnRH pulse generator (Kreisman et al., 2019). Third, we use the model to explore how perturbations, such as stressors and chronic changes in gonadal steroids and glucocorticoid levels can disrupt normal rhythmicity and lead to dysregulations that propagate from one neuroendocrine system to the other. To do so, we consider typical restraint stress signals (Li X. et al., 2004) to brain regions that are connected to the hypothalamus, thus affecting both the HPA and HPG axes (Li et al., 2005; Herod et al., 2011b). Importantly, our model considers the signaling role of regulatory neuropeptides (e.g., Neurokinin-B and Dynorphin) within the KNDy neural network in stress-induced suppression of the GnRH pulse generator (Lehman et al., 2010; Grachev et al., 2014; Voliotis et al., 2019), which has implications for our understanding of how stress signals are decoded by the reproductive axis. Lastly, we predict an increase of the estrous cycle length under hiCORT and discuss how our model can help understand the mechanisms allowing robust control of ovulation despite the effect of stressors.

2. Model and Methods

2.1. Mathematical Modeling

We focus on the systems-level outputs and cross-regulation of the stress and reproductive axes, which in turn we model as a network of coupled oscillators (Figure 1B). We modeled this through a system of Ordinary Differential Equations (ODEs), where each oscillator represents an aspect of neuroendocrine rhythmic activity that can be characterized by a phase φ, a frequency ω, and an amplitude A. Our model consists of a master circadian oscillator in the hypothalamus, a glucocorticoid (CORT) oscillator with ultradian rhythmicity driven by the circadian oscillator, a pulse generator oscillator governed by the Kisspeptin, Neurokinin B, Dynorphin (KNDy) network regulating pulses of GnRH secretion, and an oscillator representing the estrous cycle. The equations for these oscillators are listed below, with coupling functions, parameter values, and further details of the model development described in the Supplementary Material.

2.1.1. Circadian Cycle

A fixed period hypothalamic oscillator to control the circadian rhythm of CORT:

$\begin{array}{ll}\frac{d}{dt}{\phi }_H={\omega }_{H_0},& (1)\end{array}$

where φ_H is the hypothalamic phase and ω_H0 is the natural frequency of the hypothalamic circadian drive.

2.1.2. CORT Oscillator

Accounts for CORT ultradian oscillations originating from the pituitary-adrenal feedback loop (Walker J. J. et al., 2010). Its dynamics can be affected by stressors, exogenous CORT, and the estrous cycle. The phase φ_C is given by:

$\begin{array}{ll}\frac{d}{dt}{\phi }_C={\omega }_{C_0}-\alpha s({\phi }_H),& (2)\end{array}$

where ω_C0 is the natural frequency of CORT ultradian oscillations, s(φ_H) is a function accounting for a transient acute stressor (equal to zero in the absence of stress), and α is a scaling factor accounting for how strongly such stressor temporarily disrupts CORT ultradian rhythmicity. The amplitude A_C is given by:

$\begin{array}{ll}\frac{d}{dt}{A}_C=f_H({\phi }_H)\frac{{A_E}^n}{{A_E}^n+{K_E}^n}-A_C+s({\phi }_H),& (3)\end{array}$

where f_H(φ_H) is a function representing hypothalamic circadian modulation. A_E is the amplitude of the estrous cycle (representative of the level of sex steroids) which modulates A_C through a Hill type function with coefficient n and half-maximum constant K_E.

2.1.3. Pulse Generator

Accounts for the activity of the GnRH pulse generator. Its frequency is modulated by stressors, CORT levels, and the activity of the KNDy network (Voliotis et al., 2019), which is in turn influenced by the phase of the estrous cycle. The phase φ_PG is given by:

$\begin{array}{ll}\frac{d}{dt}{\phi }_{PG}={\omega }_{PG},& (4)\end{array}$

where ω_PG denotes the varying frequency of the pulse generator. This is given by:

$\begin{array}{ll}\frac{d}{dt}{\omega }_{PG}={\omega }_{P{G}_m}{f}_K(N,D,\tilde{C})-{\omega }_{PG},& (5)\end{array}$

where ω_PGm is the maximum frequency of the pulse generator and $f_K(N,D,\tilde{C})$ is a function accounting for the regulation from the KNDy network and CORT. Equations for the excitatory (N; e.g., Neurokinin B and glutamate) and inhibitory (D; e.g., Dynorphin) signals regulating the frequency of the KNDy network, and the slow genomic CORT effects ($\tilde{C}$) are given in the Supplementary Material.

2.1.4. Estrous Cycle

Accounts for the activity of the reproductive cycle. The phase φ_E is given by:

$\begin{array}{ll}\frac{d}{dt}{\phi }_E={\omega }_{E_m}{f}_{PG}({\omega }_{PG}),& (6)\end{array}$

where f_PG(ω_PG) is a function accounting for the effects of the pulse generator and ω_Em is the maximum frequency of the estrous cycle. The amplitude A_E is given by:

$\begin{array}{ll}\frac{d}{dt}{A}_E=\epsilon +\beta f_E({\phi }_E)-A_E,& (7)\end{array}$

where ε is the basal activity of the estrous cycle, f_E(φ_E) is a function representing the effects of the estrous cycle, and β is a scaling factor accounting for the strength of such effects.

2.2. Computer Simulations and Parameter Estimation

To simplify our analysis, CORT oscillations were normalized to the maximum levels observed in physiological conditions. That is, the CORT amplitude, which is modulated by the circadian drive, spans the range between 0 and 1 unless stressors or exogenous CORT act upon it. Similarly, the activity of the PG was represented by normalized oscillations, with a frequency that changes periodically according to the different stages of the estrous cycle.

For the scenarios in sections 3.2 and 3.3, we model the estrous cycle regulation of its amplitude and its effects on the KNDy network through a skewed sinusoidal function of the phase φ_E. This is given by:

$\begin{array}{ll}{f}_E({\phi }_E)=\{\begin{array}{ll}{\sin }^2\left({\phi }_E-\sigma {\sin }^2({\phi }_E)\right),\hfill & \text{Normal physiological}\hfill \\ 0.99,\hfill & \text{OVX}\hfill \\ 0.2,\hfill & {\text{OVX}+\text{E}}_2\hfill \end{array}& (8)\end{array}$

where σ is the skewness of the estrous cycle. We fix f_E(φ_E) to a constant value to simulate the OVX and OVX + E₂ scenarios (Seale et al., 2004b; Kreisman et al., 2019). Note that in those cases, parameters ε and β in the equation for A_E also need to change as indicated in Supplementary Table 3 to reflect estrous activity expected in the diestrus and proestrus phase.

The model equations were numerically solved and analyzed in MATLAB R2020a using ode45 routines. Details of the mathematical model development and parameter values are described in the Supplementary Material. The model parameters were estimated from the literature where available and manually calibrated to reproduce experimental observations of CORT and reproductive rhythms in rodents.

No new data involving animal or human subjects is presented in this paper.

3. Results

3.1. Normal Physiological HPA and HPG Rhythms

We calibrate the model parameters to reproduce physiological HPA and HPG rhythms observed in rats (Walker J. et al., 2010). Accordingly, our model simulates CORT oscillations with a 75 min period, while the amplitude of these ultradian pulses is modulated in a circadian manner, reaching a maximum at the start of the dark period (Figure 2A). Furthermore, one full estrous cycle lasts ~4 days, matching the average cycle length measured in rats (McClintock, 1984). A recent study using fiber photometry calcium imaging from arcuate kisspeptin neurons in mice revealed the dynamic modulation of GnRH pulse frequency along the estrous cycle (McQuillan et al., 2019). Following these findings, the activity of the PG in the model remains inhibited (below 1 pulse/h) during the post-ovulatory, estrous phase, rises steeply at the start of metestrus, and levels off at 2 pulses/h for the rest of the cycle (Figure 2B).

FIGURE 2

Figure 2. The model reproduces physiological rhythms in the HPA and HPG axis. (A) Normalized CORT levels as a function of time. The light-dark cycle is represented with intermittent black bars on the top. (B) Normalized pulse generator activity (blue) and pulse generator frequency (red) as a function of time. The phases of the estrous cycle are marked on the top: estrus (E); metestrus (M); diestrus (D); and proestrus (P).

3.2. Recovery of CORT Dynamics Following Ovariectomy

Previous findings suggest that gonadal steroids are integral to the increased CORT levels seen in females compared to males. This has been demonstrated by showing the effects of estrogen replacement in recovering physiological CORT levels following ovariectomy in rats (Seale et al., 2004b). We investigate the dynamic effects of these hormones by simulating the inhibition of HPA axis activity resulting from ovariectomy (OVX) and its restitution following 17β-estradiol (E₂) replacement. In the model, this is achieved by replacing the influx term in the right hand side of Equation (7) by a constant term representing a drop in E₂ levels following OVX (causing A_E to drop down to a constant level of 2% from the estrous peak) and by replacing the periodic sensitivity of the KNDy network to the estrous phase by a constant low value (Supplementary Material). The model predicts a drop in CORT levels down to ~30% from its physiological value without loss of circadian or ultradian CORT rhythmicity while keeping the PG frequency at a high constant value of 2 pulses/h (Figure 3A). We then simulated the effects of an E₂ pellet on OVX rats by increasing the constant value of the influx term in the right hand side of Equation (7) (98% from physiological A_E) and increasing the sensitivity of the KNDy network to the estrous phase (φ_E) by a constant value (Supplementary Material). In agreement with Seale et al. (2004b), the model predicts recovery of physiological CORT levels without loss of circadian or ultradian CORT rhythmicity while marginally reducing the PG frequency just below 2 pulses/h (Figure 3B).

FIGURE 3

Figure 3. The model explains how E₂ replacement recovers physiological CORT levels in OVX rats. (A) Simulated OVX reduced CORT oscillations down to ~30% of the maximum physiological levels while keeping a constant high PG activity. (B) Simulated OVX + E₂ recovered CORT oscillations to physiological levels while keeping a constant high PG activity.

3.3. Estradiol-Mediated Inhibition of HPG Dynamics by High CORT Doses

In a recent study, Kreisman et al. (2019) investigated the effect of chronic CORT administration on LH pulsatility and demonstrated the importance of gonadal steroid hormones in mediating the inhibitory effect of CORT on the HPG axis. The study showed that a pellet delivering a high dose of CORT over 48 h in OVX mice has no effect on LH pulsatility, whereas a significant reduction of LH pulse frequency is observed in OVX animals treated with a 17β-estradiol silastic implant (OVX + E₂). In our model, we accounted for the OVX and OVX + E₂ scenarios as described in the previous section, while the constantly high CORT levels were achieved by replacing the effective CORT levels modulating the KNDy network by a constant high value estimated from Kreisman et al. (2019) (see Supplementary Material).

Figure 4 illustrates the differential effect of chronically elevated CORT levels on the GnRH pulse generator frequency in OVX vs. OVX + E₂ animals. In the case of OVX animals, elevated CORT levels do not alter the frequency of the pulse generator, whereas in OVX animals treated with estradiol the frequency is halved for as long as CORT levels are elevated. This effect is linked to the modulation of the GnRH pulse generator by gonadal steroids, which sensitize the system to inhibitory signals, such as CORT or acute stressors as we show below (Figure 4B).

FIGURE 4

Figure 4. The model reproduces estradiol-mediated inhibition of PG activity following high doses of CORT. (A) High exogenous CORT over 48 h does not affect the PG dynamics in OVX mice. (B) In the presence of estradiol, high CORT doses temporarily reduce PG activity in OVX mice.

3.4. Acute Stress Effects on the HPA and HPG Axes Depend on the Estrous Cycle Phase

To study the effect of acute stress on the dynamics of the HPA and HPG axes, we extend the model to include transient stress-related neuronal inputs affecting both axes (Yang et al., 2017). In our model, we account for these transient inputs by simulating a 2 h square pulse of amplitude 1, equivalent to a restraint stressor causing a CORT increase from its circadian nadir up to its circadian peak (Kitchener et al., 2004). The stressor affects the phase and amplitude of the CORT rhythm (function s(φ_H) in Equations 2 and 3) as well as the frequency of the GnRH pulse generator (function $f_K(N,D,\tilde{C})$ in Equation 5 and Supplementary Material).

Figures 5A,B illustrate the effect that 2 h of stress activation has on the dynamics of the HPA and HPG axes when applied at different times along the cycle. Both CORT and GnRH pulse generator responses are dependent on the timing of the input pulse (Figure 5C). The amplitude of the CORT response shows a circadian dependency with stressors delivered during the circadian peak eliciting a stronger response. The GnRH pulse generator frequency response to acute stressors depends on the phase of the estrous cycle. In particular, the frequency of the pulse generator appears most sensitive to stressors during estrus to early diestrus phases, with little or no effect during the mid-cycle phase. This differential effect of acute stress on the frequency of the GnRH pulse generator activity highlights the cycle dependent modulation of the pulse generator dynamics, which makes the pulse generator more robust to perturbations in the diestrus and proestrus phases (Figure 5D).

FIGURE 5

Figure 5. The effect of acute stress on the dynamics of the HPA and HPG axes. (A,B) CORT levels and PG activity in response to a transient (2 h long) stressor initiated at two different times. (C) Peak CORT levels (black line) and mean PG frequency (continuous red line) elicited by a 2 h long acute stressor as a function of the time at which the stressor arrives during the estrous cycle. The PG frequency without any stress perturbation is shown for comparison (dashed red line). (D) State space diagram describing the effect of acute stress on the dynamics of the pulse generator. Points mark different stages along the estrus cycle: estrus midpoint (E); metestrus midpoint (M); and diestrus midpoint (D). The shaded gray area denotes the region of the state space corresponding to frequencies above 1 pulse/h under normal physiological conditions. Acute stress shrinks this region (red shaded area), but the dynamics of the pulse generator maintains robustness to perturbations during the diestrus phase.

3.5. CORT Excess Increases the Length of the Estrous Cycle and Modulates Responses to Acute Stressors

Last, we used the model to predict the effects of high excess CORT (hiCORT) —mimicking levels expected to be observed in people with hypercortisolism—on the estrous cycle. To do this, we considered the increase in baseline and maximum CORT amplitude with respect to physiological levels in humans (Vagnucci, 1979) and implemented the equivalent increase ratios for our simulations of CORT dynamics in rodents (Supplementary Material). Evidence from high frequency sampling in humans shows hypercortisolism is associated with a reduction in the ultradian period of CORT oscillations (Van Aken et al., 2005). Accordingly, we also adjusted this parameter when modeling hiCORT, while keeping circadian oscillations and all other parameters unchanged. Our simulations predict an increase in the period of the estrous cycle from a physiological value of T_phys = 4.2 days up to T_hiC = 5.1 days under hiCORT, which is equivalent to a ~21% increase in the estrous cycle length (Figures 6A,C).

FIGURE 6

Figure 6. Enduring and transient dynamic changes under hiCORT. (A) CORT and PG rhythms in physiological conditions. (B) Mean PG frequency and maximum CORT levels elicited by 2 h long hypothalamic stressors arriving at different times across the estrous cycle. (C) Under hiCORT, the range of CORT levels is increased and the estrous cycle peak is delayed by ~21% compared to physiological values. (D) Mean PG frequency and maximum CORT elicited by 2 h long stressors under hiCORT. In this scenario, the period in which PG activity remains unchanged by stressors starts about half a day later and is shortened compared to normal physiological conditions.

We then used the model to investigate the transient changes in the GnRH pulse generator frequency and CORT amplitude elicited by exogenous acute stressors under physiological conditions and hiCORT. In particular, we look at the effects of the timing of stressors within the estrous cycle. To do this, we calculated frequency and amplitude response curves by simulating a 2 h long stressor elicited at different stages across the estrous cycle using 30 min time steps. In the physiological scenario (Figure 6B), the model predicts that acute stressors suppress PG activity during most of the estrous cycle except during the diestrus and early proestrus phases. These stressors also elicit an increase in peak CORT levels to a range between 1 and 2. While the model under the hiCORT scenario predicts a similar behavior, the region where PG activity remains unaffected by acute stressors is reduced and delayed by about half a day compared to its physiological counterpart (Figure 6D and Supplementary Figure 2). This is not surprising considering that the model also predicts that hiCORT prolongs the estrous cycle. Regarding the CORT response to stressors under hiCORT, our model predicts a ~2 to ~3.5 increase in CORT levels compared to the normal physiological scenario. This is due to a compounded effect of CORT surges over an excess CORT baseline.

4. Discussion

We developed and studied a mathematical model that integrates components of the stress and reproductive axes at different spatial and temporal scales, from the molecular intricacies of the KNDy network, to GnRH and CORT oscillations, up to the estrous cycle (Figure 1A). Previous mathematical models of the HPA and HPG axes either focus on a specific process within an axis, or consider them as a whole, but isolated from each other (Walker J. J. et al., 2010; Spiga et al., 2017; Voliotis et al., 2019; Clément et al., 2020). In contrast, our model integrates these neuroendocrine axes by considering the complex interactions between them as a network of interlocked oscillators, hence enabling us to integrate different physiological observations and experiments into a single coherent theoretical framework and study the effect of transient perturbations on the overall dynamics. In particular, our model postulated a network architecture (Figure 1B) that reflects physiological observations of ultradian and circadian CORT rhythms, as well as ultradian and infradian rhythms of the GnRH pulse generator (Figure 2). The model reproduced the effects of ovarian hormone removal (OVX) and restitution (E₂) on the HPA and HPG axes dynamics, both under physiological conditions (Seale et al., 2004b) (Figure 3) and under exogenous CORT excess (Kreisman et al., 2019) (Figure 4).

In addition to these slow timescale perturbations, we also investigated the fast timescale perturbations elicited by acute stressors. Our model predicted that exogenous stress perturbations not only cause transient increases in CORT levels, but also transiently inhibit GnRH pulse generator activity with the magnitude of this inhibition being dependent on the estrous and circadian phases (Figure 5). This has important implications about understanding how the timing of a stressor affects its ability to temporarily suppress the GnRH/LH ovulatory surge. According to our model, the pulse generator activity is robust to stress perturbations arriving between the diestrus and early proestrus stage, but is fragile to stressors arriving at estrus and metestrus stages (Figure 5C). Uncovering the origin of this robustness is beyond the scope of our phenomenological model, but we can speculate that molecular mechanisms ensure the resilience of the reproductive cycle during the key stages leading to ovulation. While our model suggests that the GnRH/LH surge should be delayed under frequent exposure to stressors, if the exposure occurs too close to the proestrus stage then these resilience molecular mechanisms ensure the surge continues as normal and triggers ovulation (Wagenmaker et al., 2010; Wagenmaker and Moenter, 2017).

We also used the model to investigate the potential detrimental effects on fertility elicited by chronic hypercortisolism. Our model predicted that hiCORT delays the increase in activity of the GnRH pulse generator, effectively prolonging the estrous cycle (Figures 6A,C). While evidence suggests that HPA axis hyperactivity—and specifically, increased circulating glucocorticoids—are unlikely to be the sole mechanism behind stress-induced reproductive dysfunction (Herod et al., 2011a), our simulations show the cycle length depends on the GnRH pulse generator's sensitivity to CORT (Supplementary Figure 1A). Thus, our model provides insight into how for example a hyper-sensitized HPG axis may explain amenorrhea secondary to high serum cortisol levels (Ding et al., 1988; Suh et al., 1988; Saketos et al., 1993). Interestingly, our model predicted that a period of robustness of the GnRH pulse generator in the presence of stressors is preserved under hiCORT, albeit the robust period occurs about half a day later in the cycle and is shorter in duration. Our model simulations of pulse generator activity suggest that prolonging the estrous cycle as predicted under hiCORT arises from a combination of longer estrus and metestrus stages while diestrus and proestrus stages are shortened (Figures 6B,D).

Our model considers essential features of HPA and HPG axes oscillators in a phenomenological way. This approach facilitates the simulation of a range of physio-pathological scenarios, but inevitably imposes certain limitations. In contrast to mechanistic models where parameters are often linked to chemical kinetic rates, the parameters in our model represent natural and maximum frequencies, phase relationships, as well as the coupling strengths between oscillators and sensitivities to perturbations. While our phenomenological approach limits the ability of the model to support discovery of specific molecular mechanisms, it can be used to suggest experiments that explore systems level properties involving both neuroendocrine axes. For example, evidence shows that in addition to exhibiting circadian and ultradian fluctuations, CORT levels also change across the estrous cycle, with maximum levels around the diestrus and proestrus phase (Carey et al., 1995; Atkinson and Waddell, 1997; Pilorz et al., 2009). While our model lacks the level of detail to describe the molecular mechanisms that underpin estrous changes on CORT, it does suggest this is mediated by a regulatory link from the estrous oscillator to the CORT oscillator, thus inferring that ovarian steroids may be the culprit of estrous regulation of CORT instead of the hypothalamic GnRH pulse generator (Figure 1). In our model, we only explored the scenario where the strength of this regulatory link allows for strong perturbations (e.g., stressors, OVX, E₂) in the estrous oscillator to have an impact on the CORT dynamics, but not from milder estrous regulation of CORT levels (Supplementary Figure 1). We speculate that combining mechanistic modeling with experimental physiology to investigate the effects of estradiol and progesterone on CORT may uncover the origins of its estrous cycle modulation. The experiments could test the dosing effect, timing, and combined sensitivity of gonadal steroids on circadian CORT levels across the estrous cycle. The mechanistic model could in turn help understand the robustness of such regulatory mechanism to perturbations (Wagenmaker et al., 2010), and predict the scenarios in which chronodisruptions would lead to disease.

We believe that the first generation mathematical model presented here could be used to inform further investigations into the timing of stress perturbations in reproductive health, including dysregulations induced by strenuous exercise (Ding et al., 1988), mood disorders (Young and Korszun, 2002), as well as clinical interventions, such as in vitro fertilization (Massey et al., 2016). Our model is the latest of a class of mathematical models that can support or replace animal studies in endocrinology (Zavala et al., 2019). It can also help design new studies that reduce the number of experiments necessary to refine our understanding of the HPA and HPG axis. Furthermore, computational models like ours can be used to contextualize the results of clinical studies where experimentation is not possible. This can be done in combination with a range of tools from network physiology and machine learning that consider the dynamic links between coupled body rhythms, such as body temperature and sleep (Bashan et al., 2012; Bartsch et al., 2015; Ivanov et al., 2016). We anticipate that healthcare technologies, such as wearable devices and smartphone apps collecting vast amounts of data on body rhythms, together with computer algorithms characterizing inter-individual variability, will help refine and personalize neuroendocrinological models (Kim et al., 2020; Li et al., 2020; Wang et al., 2020).

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

EZ, MV, XL, JRT, SL, KO'B, and KT-A conceived and designed the study. EZ, MV, TZ, JT, JW, and KT-A developed the mathematical model. EZ, MV, and TZ performed the analytical calculations and computer simulations, and wrote the manuscript with support of all authors.

Funding

This work was funded by the Medical Research Council (MRC) fellowship MR/P014747/1 (to EZ), MRC fellowship MR/N008936/1 (to JW), MRC grant MR/N022637/1 (to SL and KO'B), Biotechnology and Biological Sciences Research Council (BBSRC) grant BB/S001255/1 (to MV, XL, KO'B, and KT-A), Engineering and Physical Sciences Research Council (EPSRC) grant EP/N014391/2 (to EZ, MV, TZ, JT, JW, JRT, SL, and KT-A), and the Wellcome Trust Grant WT105618MA (to JW, JRT, and KT-A).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2020.598845/full#supplementary-material

References

Acevedo-Rodriguez, A., Kauffman, A., Cherrington, B., Borges, C., Roepke, T., and Laconi, M. (2018). Emerging insights into hypothalamic-pituitary-gonadal axis regulation and interaction with stress signalling. J. Neuroendocrinol. 30:e12590. doi: 10.1111/jne.12590

CrossRef Full Text | Google Scholar

Atkinson, H. C., and Waddell, B. J. (1997). Circadian variation in basal plasma corticosterone and adrenocorticotropin in the rat: sexual dimorphism and changes across the estrous cycle. Endocrinology 138, 3842–3848. doi: 10.1210/endo.138.9.5395

CrossRef Full Text | Google Scholar

Ayrout, M., Le Billan, F., Grange-Messent, V., Mhaouty-Kodja, S., Lombès, M., and Chauvin, S. (2019). Glucocorticoids stimulate hypothalamic dynorphin expression accounting for stress-induced impairment of GnRH secretion during preovulatory period. Psychoneuroendocrinology 99, 47–56. doi: 10.1016/j.psyneuen.2018.08.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartsch, R. P., Liu, K. K., Bashan, A., and Ivanov, P. C. (2015). Network physiology: how organ systems dynamically interact. PLoS ONE 10:e0142143. doi: 10.1371/journal.pone.0142143

PubMed Abstract | CrossRef Full Text | Google Scholar

Bashan, A., Bartsch, R. P., Kantelhardt, J. W., Havlin, S., and Ivanov, P. C. (2012). Network physiology reveals relations between network topology and physiological function. Nat. Commun. 3, 1–9. doi: 10.1038/ncomms1705

PubMed Abstract | CrossRef Full Text | Google Scholar

Carey, M., Deterd, C., De Koning, J., Helmerhorst, F., and De Kloet, E. (1995). The influence of ovarian steroids on hypothalamic-pituitary-adrenal regulation in the female rat. J. Endocrinol. 144, 311–321. doi: 10.1677/joe.0.1440311

PubMed Abstract | CrossRef Full Text | Google Scholar

Churilov, A. N., Milton, J., and Salakhova, E. R. (2020). An integrate-and-fire model for pulsatility in the neuroendocrine system. Chaos 30:083132. doi: 10.1063/5.0010553

PubMed Abstract | CrossRef Full Text | Google Scholar

Clément, F., Crépieux, P., Yvinec, R., and Monniaux, D. (2020). Mathematical modeling approaches of cellular endocrinology within the hypothalamo-pituitary-gonadal axis. Mol. Cell. Endocrinol. 518:110877. doi: 10.1016/j.mce.2020.110877

PubMed Abstract | CrossRef Full Text | Google Scholar

Crofford, L. J., Jacobson, J., and Young, E. (1999). Modeling the involvement of the hypothalamicpituitary-adrenal and hypothalamic-pituitary-gonadal axes in autoimmune and stress-related rheumatic syndromes in women. J. Womens Health 8, 203–215. doi: 10.1089/jwh.1999.8.203

PubMed Abstract | CrossRef Full Text | Google Scholar

Ding, J.-H., Sheckter, C. B., Drinkwater, B. L., Soules, M. R., and Bremner, W. J. (1988). High serum cortisol levels in exercise-associated amenorrhea. Ann. Intern. Med. 108, 530–534. doi: 10.7326/0003-4819-108-4-530

PubMed Abstract | CrossRef Full Text | Google Scholar

Grachev, P., Li, X., Hu, M., Li, S., Millar, R. P., Lightman, S. L., et al. (2014). Neurokinin b signaling in the female rat: a novel link between stress and reproduction. Endocrinology 155, 2589–2601. doi: 10.1210/en.2013-2038

PubMed Abstract | CrossRef Full Text | Google Scholar

Herod, S. M., Dettmer, A. M., Novak, M. A., Meyer, J. S., and Cameron, J. L. (2011a). Sensitivity to stress-induced reproductive dysfunction is associated with a selective but not a generalized increase in activity of the adrenal axis. Am. J. Physiol. Endocrinol. Metab. 300, E28–E36. doi: 10.1152/ajpendo.00223.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Herod, S. M., Pohl, C. R., and Cameron, J. L. (2011b). Treatment with a CRH-R1 antagonist prevents stress-induced suppression of the central neural drive to the reproductive axis in female macaques. Am. J. Physiol. Endocrinol. Metab. 300, E19–E27. doi: 10.1152/ajpendo.00224.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanov, P. C., Liu, K. K., and Bartsch, R. P. (2016). Focus on the emerging new fields of network physiology and network medicine. New J. Phys. 18:100201. doi: 10.1088/1367-2630/18/10/100201

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, D. W., Zavala, E., and Kim, J. K. (2020). Wearable technology and systems modeling for personalized chronotherapy. Curr. Opin. Syst. Biol. 21, 9–15. doi: 10.1016/j.coisb.2020.07.007

CrossRef Full Text | Google Scholar

Kitchener, P., Di Blasi, F., Borrelli, E., and Piazza, P. V. (2004). Differences between brain structures in nuclear translocation and dna binding of the glucocorticoid receptor during stress and the circadian cycle. Eur. J. Neurosci. 19, 1837–1846. doi: 10.1111/j.1460-9568.2004.03267.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kreisman, M., McCosh, R., Tian, K., Song, C., and Breen, K. (2019). Estradiol enables chronic corticosterone to inhibit pulsatile LH secretion and suppress Kiss1 neuronal activation in female mice. Neuroendocrinology 110, 501–516. doi: 10.1159/000502978

PubMed Abstract | CrossRef Full Text | Google Scholar

Lehman, M. N., Coolen, L. M., and Goodman, R. L. (2010). Minireview: Kisspeptin/Neurokinin B/Dynorphin (KNDy) cells of the arcuate nucleus: a central node in the control of gonadotropin-releasing hormone secretion. Endocrinology 151, 3479–3489. doi: 10.1210/en.2010-0022

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, K., Urteaga, I., Wiggins, C. H., Druet, A., Shea, A., Vitzthum, V. J., et al. (2020). Characterizing physiological and symptomatic variation in menstrual cycles using self-tracked mobile-health data. NPJ Digit. Med. 3, 1–13. doi: 10.1038/s41746-020-0269-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Edward, J., Mitchell, J., Shao, B., Bowes, J., Coen, C., et al. (2004). Differential effects of repeated restraint stress on pulsatile lutenizing hormone secretion in female fischer, lewis and wistar rats. J. Neuroendocrinol. 16, 620–627. doi: 10.1111/j.1365-2826.2004.01209.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X. F., Bowe, J. E., Lightman, S. L., and O'Byrne, K. T. (2005). Role of corticotropin-releasing factor receptor-2 in stress-induced suppression of pulsatile luteinizing hormone secretion in the rat. Endocrinology 146, 318–322. doi: 10.1210/en.2004-0950

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X. F., Bowe, J. E., Mitchell, J. C., Brain, S. D., Lightman, S. L., and O'Byrne, K. T. (2004). Stress-induced suppression of the gonadotropin-releasing hormone pulse generator in the female rat: a novel neural action for calcitonin gene-related peptide. Endocrinology 145, 1556–1563. doi: 10.1210/en.2003-1609

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, E., Stephens, S. B., Chaing, S., Munaganuru, N., Kauffman, A. S., and Breen, K. M. (2016). Corticosterone blocks ovarian cyclicity and the LH surge via decreased kisspeptin neuron activation in female mice. Endocrinology 157, 1187–1199. doi: 10.1210/en.2015-1711

PubMed Abstract | CrossRef Full Text | Google Scholar

Massey, A. J., Campbell, B. K., Raine-Fenning, N., Pincott-Allen, C., Perry, J., and Vedhara, K. (2016). Relationship between hair and salivary cortisol and pregnancy in women undergoing IVF. Psychoneuroendocrinology 74, 397–405. doi: 10.1016/j.psyneuen.2016.08.027

PubMed Abstract | CrossRef Full Text | Google Scholar

McClintock, M. K. (1984). Estrous synchrony: modulation of ovarian cycle length by female pheromones. Physiol. Behav. 32, 701–705. doi: 10.1016/0031-9384(84)90181-1

PubMed Abstract | CrossRef Full Text | Google Scholar

McCosh, R. B., Breen, K. M., and Kauffman, A. S. (2019). Neural and endocrine mechanisms underlying stress-induced suppression of pulsatile LH secretion. Mol. Cell. Endocrinol. 498:110579. doi: 10.1016/j.mce.2019.110579

PubMed Abstract | CrossRef Full Text | Google Scholar

McQuillan, H. J., Han, S. Y., Cheong, I., and Herbison, A. E. (2019). GnRH pulse generator activity across the estrous cycle of female mice. Endocrinology 160, 1480–1491. doi: 10.1210/en.2019-00193

PubMed Abstract | CrossRef Full Text | Google Scholar

Phumsatitpong, C., and Moenter, S. M. (2018). Estradiol-dependent stimulation and suppression of gonadotropin-releasing hormone neuron firing activity by corticotropin-releasing hormone in female mice. Endocrinology 159, 414–425. doi: 10.1210/en.2017-00747

PubMed Abstract | CrossRef Full Text | Google Scholar

Pilorz, V., Steinlechner, S., and Oster, H. (2009). Age and oestrus cycle-related changes in glucocorticoid excretion and wheel-running activity in female mice carrying mutations in the circadian clock genes per1 and per2. Physiol. Behav. 96, 57–63. doi: 10.1016/j.physbeh.2008.08.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Saketos, M., Sharma, N., and Santoro, N. F. (1993). Suppression of the hypothalamic-pituitary-ovarian axis in normal women by glucocorticoids. Biol. Reprod. 49, 1270–1276. doi: 10.1095/biolreprod49.6.1270

PubMed Abstract | CrossRef Full Text | Google Scholar

Seale, J., Wood, S., Atkinson, H., Bate, E., Lightman, S., Ingram, C., et al. (2004a). Gonadectomy reverses the sexually diergic patterns of circadian and stress-induced hypothalamic-pituitary-adrenal axis activity in male and female rats. J. Neuroendocrinol. 16, 516–524. doi: 10.1111/j.1365-2826.2004.01195.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Seale, J., Wood, S., Atkinson, H., Harbuz, M., and Lightman, S. (2004b). Gonadal steroid replacement reverses gonadectomy-induced changes in the corticosterone pulse profile and stress-induced hypothalamic-pituitary-adrenal axis activity of male and female rats. J. Neuroendocrinol. 16, 989–998. doi: 10.1111/j.1365-2826.2004.01258.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Seale, J., Wood, S., Atkinson, H., Harbuz, M., and Lightman, S. (2005a). Postnatal masculinization alters the hpa axis phenotype in the adult female rat. J. Physiol. 563, 265–274. doi: 10.1113/jphysiol.2004.078212

PubMed Abstract | CrossRef Full Text | Google Scholar

Seale, J., Wood, S., Atkinson, H., Lightman, S., and Harbuz, M. (2005b). Organizational role for testosterone and estrogen on adult hypothalamic-pituitary-adrenal axis activity in the male rat. Endocrinology 146, 1973–1982. doi: 10.1210/en.2004-1201

PubMed Abstract | CrossRef Full Text | Google Scholar

Spiga, F., Zavala, E., Walker, J. J., Zhao, Z., Terry, J. R., and Lightman, S. L. (2017). Dynamic responses of the adrenal steroidogenic regulatory network. Proc. Natl. Acad. Sci. U.S.A. 114, E6466–E6474. doi: 10.1073/pnas.1703779114

PubMed Abstract | CrossRef Full Text | Google Scholar

Suh, B., Liu, J., Berga, S., Quigley, M., Laughlin, G., and Yen, S. (1988). Hypercortisolism in patients with functional hypothalamic-amenorrhea. J. Clin. Endocrinol. Metab. 66, 733–739. doi: 10.1210/jcem-66-4-733

PubMed Abstract | CrossRef Full Text | Google Scholar

Vagnucci, A. H. (1979). Analysis of circadian periodicity of plasma cortisol in normal man and in Cushing's syndrome. Am. J. Physiol. Regul. Integr. Compar. Physiol. 236, R268–R281. doi: 10.1152/ajpregu.1979.236.5.R268

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Aken, M., Pereira, A., Van Thiel, S., Van Den Berg, G., Frolich, M., Veldhuis, J. D., et al. (2005). Irregular and frequent cortisol secretory episodes with preserved diurnal rhythmicity in primary adrenal Cushing's syndrome. J. Clin. Endocrinol. Metab. 90, 1570–1577. doi: 10.1210/jc.2004-1281

PubMed Abstract | CrossRef Full Text | Google Scholar

Voliotis, M., Li, X. F., De Burgh, R., Lass, G., Lightman, S. L., O'Byrne, K. T., et al. (2019). The origin of GnRH pulse generation: an integrative mathematical-experimental approach. J. Neurosci. 39, 9738–9747. doi: 10.1523/JNEUROSCI.0828-19.2019

PubMed Abstract | CrossRef Full Text | Google Scholar

Wagenmaker, E. R., Breen, K. M., Oakley, A. E., Tilbrook, A. J., and Karsch, F. J. (2010). The estrous cycle of the ewe is resistant to disruption by repeated, acute psychosocial stress. Biol. Reprod. 82, 1206–1215. doi: 10.1095/biolreprod.109.078774

PubMed Abstract | CrossRef Full Text | Google Scholar

Wagenmaker, E. R., and Moenter, S. M. (2017). Exposure to acute psychosocial stress disrupts the luteinizing hormone surge independent of estrous cycle alterations in female mice. Endocrinology 158, 2593–2602. doi: 10.1210/en.2017-00341

PubMed Abstract | CrossRef Full Text | Google Scholar

Walker, J., Terry, J., Tsaneva-Atanasova, K., Armstrong, S., McArdle, C., and Lightman, S. (2010). Encoding and decoding mechanisms of pulsatile hormone secretion. J. Neuroendocrinol. 22, 1226–1238. doi: 10.1111/j.1365-2826.2010.02087.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Walker, J. J., Spiga, F., Waite, E., Zhao, Z., Kershaw, Y., Terry, J. R., et al. (2012). The origin of glucocorticoid hormone oscillations. PLoS Biol. 10:e1001341. doi: 10.1371/joural.pbio.1001341

PubMed Abstract | CrossRef Full Text | Google Scholar

Walker, J. J., Terry, J. R., and Lightman, S. L. (2010). Origin of ultradian pulsatility in the hypothalamic-pituitary-adrenal axis. Proc. R. Soc. B Biol. Sci. 277, 1627–1633. doi: 10.1098/rspb.2009.2148

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y. X., Arvizu, M., Rich-Edwards, J. W., Stuart, J. J., Manson, J. E., Missmer, S. A., et al. (2020). Menstrual cycle regularity and length across the reproductive lifespan and risk of premature mortality: prospective cohort study. BMJ 371:m3464. doi: 10.1136/bmj.m3464

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J. A., Song, C. I., Hughes, J. K., Kreisman, M. J., Parra, R. A., Haisenleder, D. J., et al. (2017). Acute psychosocial stress inhibits LH pulsatility and Kiss1 neuronal activation in female mice. Endocrinology 158, 3716–3723. doi: 10.1210/en.2017-00301

PubMed Abstract | CrossRef Full Text | Google Scholar

Young, E., and Korszun, A. (2002). The hypothalamic-pituitary-gonadal axis in mood disorders. Endocrinol. Metab. Clin. N. Am. 31, 63–78. doi: 10.1016/S0889-8529(01)00002-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zavala, E., Wedgwood, K. C., Voliotis, M., Tabak, J., Spiga, F., Lightman, S. L., et al. (2019). Mathematical modelling of endocrine systems. Trends Endocrinol. Metab. 30, 244–257. doi: 10.1016/j.tem.2019.01.008

PubMed Abstract | CrossRef Full Text | Google Scholar