The Hepato-Hypothalamic-Pituitary-Adrenal-Renal Axis: Mathematical Modeling of Cortisol’s Production,Metabolism,and Seasonal Variation

Abstract

Cortisol dynamics are governed by the integration of influences from the suprachiasmatic nucleus (SCN), the hypothalamic-pituitary-adrenal (HPA) axis, and metabolic enzymes, such as the 11β–hydroxysteroid dehydrogenase (HSD) family, which are highly expressed in hepatic and renal tissue. The coordinated regulation of cortisol dynamics is essential for the maintenance of a healthy state, and aberrant cortisol circadian rhythms are associated with various pathophysiological conditions. The duration of the light-dark cycle, or photoperiod, which regulates SCN activity, varies seasonally, and the shorter photoperiod winter season is associated with elevated cortisol levels, peak inflammatory disease incidence, and symptom exacerbation. Elevated expression and activity of 11β-HSD1 protein, assumed to also occur during the winter, have been allied with numerous inflammatory conditions. A comprehensive understanding of the communication between the underlying regulatory mechanisms of cortisol as well as how changes in their activity could lead to the development of disease is yet to be elucidated. In this work, we propose the use of a semimechanistic mathematical model to explore the impact of the hepato-hypothalamic-pituitary-adrenal-renal axis in modulating neuroendocrine-immune system dynamics. Our model predicts the predominance of a winter proinflammatory state and that genetic variations could alter 11β-HSD enzyme functionality, rendering certain subpopulations more susceptible to disease as a consequence of HPA axis dysregulation.

Glucocorticoids (GCs) are a class of steroid hormones implicated in the regulation of immune (Barnes, 1998), cardiovascular (Whitworth et al., 2005), metabolic (Shamoon et al., 1980), and behavioral (Oquendo et al., 2003) processes. Both cortisol (predominant GC in humans) and corticosterone (predominant GC in rodents) are downstream effectors of the hypothalamic-pituitary-adrenal (HPA) axis, which is a major neuroendocrine stress-response system. Physical, physiological, and environmental stressors activate the HPA axis and induce the release of corticotropin-releasing hormone (CRH) from the hypothalamus. Secreted CRH stimulates the anterior pituitary to produce and release adrenocorticotropic hormone (ACTH) into the systemic circulation that subsequently promotes the synthesis and secretion of GCs from the zona fasciculata of the adrenal cortex. The produced GCs ultimately mediate their regulatory effects by binding to 2 types of ubiquitously distributed intracellular receptors (Joëls, 2006): the glucocorticoid receptor (GR) and the mineralocorticoid receptor (MR).

Hormone production and other physiological phenomena such as the rest-activity cycle and body temperature display circadian rhythmicities that are in synchrony with the 24-h daily fluctuations of light and darkness (Chung et al., 2011). This entrainment to the external light-dark cycle is arbitrated by the suprachiasmatic nucleus (SCN) located in the anterior hypothalamus. Photosensitive retinal-ganglion cells located within the retina relay photic information to the SCN via the retinohypothalamic tract (Newman et al., 2003). The SCN, or master circadian pacemaker, translates the light signal into neuronal, endocrine, and/or metabolic signals that facilitate temporal coordination of various other biological activities. GCs are considered circadian hormones as both cortisol (Selmaoui and Touitou, 2003) and corticosterone (Malisch et al., 2008) concentrations vary predictably on a daily basis, with levels peaking just before the onset of activity. The disruption of corticosterone secretion following SCN lesioning underscores the contributing role of the central pacemaker in regulating GC production (Buijs et al., 1993).

Apart from being centrally regulated by the HPA axis and the SCN, GC availability in the periphery is controlled by the activity of an enzymatic system comprising 11β-hydroxysteroid dehydrogenase (HSD), 20β-oxoreductase, 6β-hydroxylase, 5β-reductase, 5α-reductase, and 3α-HSD (Gathercole et al., 2013). Of these, the 11β-HSD family, which consists of the 11β-HSD1 and 11β-HSD2 proteins, has been identified as a core modulator of GC bioavailability. In humans, these intracellular gatekeepers of GC action catalyze the reciprocal conversion between cortisol and its inactive metabolite, cortisone (Chapman et al., 2013). 11β-HSD2 is NAD⁺ dependent, catalyzes the oxidation of cortisol to cortisone, and is localized in aldosterone-sensitive tissues such as the kidneys, colon, and salivary glands. 11β-HSD1 is NADPH dependent; present in metabolically active tissues such as the central nervous system, liver, skeletal muscle, and adipose tissue; and catalyzes the regeneration of cortisol from cortisone. The deactivation of cortisol by 11β-HSD2 serves to prevent the illicit stimulation of the mineralocorticoid receptor by cortisol, whereas the regeneration by 11β-HSD1 amplifies cortisol action (Seckl and Walker, 2001). Altered 11β-HSD1 activity and expression have been associated with metabolic syndrome (Esteves et al., 2014), rheumatoid arthritis (RA; Hardy et al., 2008), systemic lupus erythematosus, pleuritis, acute interstitial pneumonitis (Ichikawa et al., 1977), and depression (Dekker et al., 2012). Although a local increase in 11β-HSD2 gene expression was observed in peripheral blood mononuclear cells obtained from RA patients (Olsen et al., 2004), defective 11β-HSD2 protein expression and activity are commonly associated with apparent mineralocorticoid excess, a disease characterized by symptoms such as hypertension, hypokalemia, and low plasma renin activity (White, 2001). Gene expression studies in mice emphasize the profound influence of 11β-HSD1 in regulating HPA axis dynamics as suppression of its expression results in the elevation of basal ACTH, corticosterone levels, and exacerbation of the acute restraint stress response (Harris et al., 2001; Carter et al., 2009). Restoration of the corticosterone profile and reversal of the aggravated stress response to control conditions is conferred by crossing the 11β-HSD1–deficient mice with transgenic 11β-HSD1–overexpressing mice (Paterson et al., 2007).

The regulation of GC dynamics is therefore multitiered, controlled at both the central and peripheral levels, with the central clock (i.e., the SCN) entrainable to external agents such as the light-dark cycle. These zeitgebers, or environmental cues, can therefore modulate physiological processes by way of the SCN. In temperate regions, the length of the light phase in the light-dark cycle predictably fluctuates throughout the year. The photoperiod, or light fraction of the day, which decreases toward winter and then increases from winter through summer, is a robust environmental entrainer that regulates various seasonal metabolic and endocrine changes (Hazlerigg and Wagner, 2006). Peak cortisol concentrations in the winter (Walker et al., 1997; Persson et al., 2008; Hadlow et al., 2014) reflect this seasonal plasticity of endocrine function. Furthermore, elevation of the ratio of cortisol to cortisone urinary metabolites and enhanced sensitivity to topical GC application in the winter suggests the amplification of GC action by augmented 11β-HSD1 activity (Walker et al., 1997).

A coordinated balance between stress-induced and basal regulatory mechanisms governing GC secretion is imperative for homeostasis (Dedovic et al., 2009), and abnormal cortisol dynamics are associated with RA (Perry et al., 2009), advanced stages of breast cancer (Sephton et al., 2000), type 2 diabetes (Lederbogen et al., 2011), cardiovascular disease mortality (Kumari et al., 2011), and behavioral disorders such as posttraumatic stress disorder and depression (Sriram et al., 2012). In addition, coronary heart disease (Douglas et al., 1995), inflammatory diseases such as RA (Schlesinger and Schlesinger, 2005), and seasonal affective disorder (winter depression; Lam and Levitan, 2000) display seasonal variation with elevated incidence and intensified symptoms occurring in the winter season. Experimental restrictions limit the comprehensive evaluation of the underlying mechanisms contributing to aberrant cortisol profiles in disease states, disease and cortisol production seasonality, and the regulatory influence of cortisol metabolism on HPA axis activity. Procedures to accurately quantify expression of 11β-HSD enzymes are invasive, requiring tissue explantation, while in vivo measurements of enzyme activity are extensive and typically involve radiation exposure (Andrew et al., 2002; Basu et al., 2004; Hughes et al., 2012). In this work, we employ a mathematical modeling approach to further develop the understanding of the underlying dynamics. To the authors’ knowledge, no mathematical model currently exists that explores the dynamics of 11β-HSD activities or its contributions to HPA axis activity. With cortisol as the representative effector molecule, and extending from our previous modeling efforts (Pierre et al., 2016), our integrated model explores the influence of the hepato-hypothalamic-pituitary-adrenal-renal axis (Edwards, 2012) in modulating both baseline and seasonal neuroendocrine-immune system dynamics.

Materials and Methods

A 2-compartment model was employed that resolved the system into central and peripheral compartments (Fig. 1). The hypothalamus and anterior pituitary gland of the brain form the central compartment, while the peripheral compartment consists of adrenal, renal, hepatic, and immune subsystems. Simplifying a large-scale complex system into a subset of interacting functional subsystems enables a fundamental understanding of the underlying structure and dynamics. Similar approaches have been employed by others to predict the impact of inflammatory stimuli on respiratory function (Reynolds et al., 2010) and to describe a glucose-insulin control system that considers interactions among liver, renal, pancreas, muscle, adipose tissue, and gastrointestinal subsystems to characterize postprandial physiological events (Man et al., 2007). The central compartment is entrained to an external light-dark cycle that was modeled as a step function with values of 1 and 0 during the light and dark phases, respectively. The seasonal light profiles are depicted in Supplementary Figure S1. The photoperiod was varied for each season, with the spring model’s (S1B) light schedule being the baseline condition with 12 h of light and darkness (12L/12D; Jackson et al., 2014). Autumn (S1D) was modeled with a photoperiod of 10L/14D (Nowakowska et al., 2011), winter (S1A) with 8L/16D, and summer (S1C) with a photoperiod of 16L/8D (Bodenstein et al., 2012).

Figure 1.

Schematic representation of our 2-compartment model. The central compartment is entrained to the external light/dark cycle. Cortisol in the periphery (F_periphery) entrains the molecular clock components and modulates the dynamics of the proinflammatory mediator (P) and its receptor (R_p). Peripheral cortisol availability is also regulated by the enzymatic activities of 11β-HSD1 and 11β-HSD2.

A detailed description of the developed mathematical model that characterizes the photoperiod-induced seasonal entrainment of cortisol dynamics and explores the consequential effects of interacting central and peripheral processes is contained in the Model Construction section of the Supplementary Online Material. Briefly, cortisol secretion, modulated by the photoperiod, subsequently synchronizes dynamics in the periphery to the seasonal photoperiod following binding to and activation of the MR and GR. The nuclear cortisol-glucocorticoid receptor complex (FGR[N]) entrains the dynamics of peripheral molecular clock components such as PER-CRY and BMAL1 mRNA and also suppresses transcription of the generic proinflammatory mediator mRNA (P_mRNA), whereas the cortisol-mineralocorticoid complex (FMR[N]) stimulates transcription of the proinflammatory receptor mRNA ( ${\mathrm{mRNA}}_{{\mathrm{R}}_{\mathrm{P}}}$ ). Transcription of 11β-HSD1 and 11β-HSD2 mRNA is also regulated by the cortisol-receptor complexes, and in turn, the translated enzymes modulate the bioavailability of peripheral cortisol (F_periphery) and its inactive metabolite cortisone (E).

Results

Seasonal Changes in Circadian Dynamics

Representative circadian profiles for all seasons are depicted in Figure 2. The oscillatory amplitudes of ACTH (Fig. 2A), F_periphery (Fig. 2B), E (Fig. 2C), 11β-HSD1 (Fig. 2E), and 11β-HSD2 (Fig. 2F) proteins; the proinflammatory mediator P (Fig. 2G); and its receptor R_p (Fig. 2H) were greatest for the winter model and lowest for the summer model. Alternatively, the amplitude of the molecular clock protein BMAL1 (Fig. 2D) was heightened in the summer and diminished in the winter model. 11β-HSD1 and P displayed similar acrophases (time of peak levels) near the minima of F_periphery while 11β-HSD2 and R_p rhythms were at a maxima shortly after the peak of F_periphery. The acrophases of all rhythms were delayed as the seasons progressed from summer through winter, following which they were advanced. The underlying mechanisms governing our predicted seasonal variations of F_periphery were then assessed. Figure 3 depicts the dynamics of ACTH-, 11β-HSD1–, 11β-HSD2–, and degradation enzyme–mediated rate of change of cortisol concentration for each season. $\frac{{\mathrm{k}}_{\mathrm{p}3}.{\mathrm{ACTH}}^{{\mathrm{n}}_{\mathrm{season}}}}{{\mathrm{K}}_{\mathrm{p}3}^{{\mathrm{n}}_{\mathrm{season}}}+{\mathrm{ACTH}}^{{\mathrm{n}}_{\mathrm{season}}}}$ regulates the rate of change of adrenal production, $\frac{{\mathrm{k}}_{\mathrm{cat}}.11\mathrm{\beta }-\mathrm{HSD}1.\mathrm{E}}{{\mathrm{K}}_{\mathrm{mF}}+\mathrm{E}}$ governs hepatic activation, $\frac{{\mathrm{k}}_{\mathrm{cat}2}.11\mathrm{\beta }-\mathrm{HSD}2.{\mathrm{F}}_{\mathrm{periphery}}}{{\mathrm{K}}_{\mathrm{mE}}+{\mathrm{F}}_{\mathrm{periphery}}}$ modulates the rate of renal inactivation, and $\frac{{\mathrm{V}}_{\mathrm{d}3}.{\mathrm{F}}_{\mathrm{periphery}}}{{\mathrm{K}}_{\mathrm{d}3}+{\mathrm{F}}_{\mathrm{periphery}}}$ controls the rate of enzymatic degradation of F_periphery. The upper panels show the rate of change versus substrate concentrations, while the lower panels depict the circadian variation of these rates. Increases in ACTH and F_periphery substrate concentrations enhance the rates of production (Fig. 3A) and degradation of cortisol (Fig. 3D), which eventually saturates at higher substrate concentrations. The reduced substrate concentrations for the summer model restrict the ranges of the rates of cortisol production, activation, inactivation, and degradation. Seasonal disparities were observed for the rate of adrenal cortisol production whereby lower ACTH concentrations resulted in higher rates of production in the summer model. At greater ACTH concentrations, however, the rate of cortisol production was greater in the winter model. At comparable F_periphery concentrations, the rate of degradation was seasonally independent. The enzyme-mediated rates of activation and inactivation of cortisol, arbitrated by 11β-HSD1 and 11β-HSD2, were predicted to display hysteretic behavior (Figs. 3B, 3C) with comparable substrate concentrations yielding different responses. The predicted hysteretic phenomena are attributed to the circadian variations of both the substrate and enzyme levels (Suppl. Figs. S2-S3). The slower descent of 11β-HSD1 in the early morning prolongs its availability. This, coupled with increasing E levels, results in greater changes in F_periphery activation rates when compared with later in the day. Similarly, higher 11β-HSD2 levels during the afternoon/night results in greater changes in F_periphery deactivation rates when compared with the early morning. For all mechanisms, the circadian rate of change of cortisol concentration for the winter model displayed the greatest dynamic range (Figs. 3E-3H).

Figure 2.

Seasonal circadian profiles. The winter circadian rhythms (solid gray) of ACTH (A), F_periphery (B), E (C), 11β-HSD1 (E), 11β-HSD2 (F), P (G), and Rp (H) displayed the highest amplitudes and were the most phase delayed, whereas the summer rhythms (dotted black) of these variables were phase advanced with attenuated amplitudes. The amplitude of BMAL1 (D) was elevated in the summer and reduced in the winter model.

Figure 3.

Seasonal and circadian variation of mechanisms regulating F_periphery availability. ACTH (A, E), 11β-HSD1– (B, F), 11β-HSD2– (C, G), and degradation enzyme–mediated (D, H) rate of change of cortisol concentration for each season. The upper row depicts the rate of change versus substrate concentration, while the lower row is plotted with respect to time. Both seasonal and daily variations of cortisol’s production, activation, inactivation, and degradation rates are evident, with the winter model exhibiting the greatest dynamic range. Spring (Spr) and summer (Sum) are represented by open circles in (D).

Population Seasonal Changes in Circadian Dynamics

Expanding our assessment to account for subject variability, we used our in silico 741 subject data set (see the Parameter Sampling section of the Supplementary Online Material for details on data set generation) to analyze seasonal variations in the mechanisms regulating F_periphery availability at a population level. The ACTH-, 11β-HSD1–, and 11β-HSD1–mediated and degradation terms regulating cortisol’s rate of change were numerically integrated (from t = 0 to t = 24 h) to estimate the contribution of each on cortisol’s daily levels. The ratio of cortisol activated by 11β-HSD1 to cortisol inactivated by 11β-HSD2 was also determined. Figure 4 illustrates the seasonal whisker plots for total daily cortisol regulated by each mechanism with the exclusion of 11β-HSD2–inactivated cortisol, which displayed no significant differences between seasons. ACTH-induced cortisol (Fig. 4A) was highest in the winter and autumn and lowest during the summer. Pairwise comparison (see the Statistical Analysis section of the Supplementary Online Material) revealed that no significant difference existed between the winter and autumn values. 11β-HSD1–activated cortisol in the summer (Fig. 4B) was predicted to be significantly higher than winter and autumn values, with the ratio of activated to inactivated cortisol also highest in the summer season (Fig. 4C). Enzymatic degradation of cortisol was the most pronounced in the spring season (Fig. 4D). Furthermore, as is shown in Figure 4E, our model predicted the greatest F_periphery to E ratio for the winter model followed by autumn, spring, and summer.

Figure 4.

Seasonal whisker plots for F_periphery added or removed from the system by various mechanisms and the ratio of F_periphery to E. F_periphery produced following ACTH stimulation (A) was highest in the winter, whereas 11β-HSD1–activated F_periphery (B) was highest in the summer season, which is also characterized by elevated ratios of activated to deactivated F_periphery (C). The highest F_periphery to E ratio was predicted for the winter season. Winter, spring, summer, and autumn are represented by Win, Spr, Sum, and Aut, respectively.

Subpopulation Classification and In Silico 11β-HSD1 Knockout

HPA axis activity is characterized by pronounced interindividual variability under both basal and stressed conditions. Differences in HPA axis responses underlie differential vulnerability to a variety of disorders; however, the bases of these variations are poorly understood (Federenko et al., 2004). To investigate the influence of cortisol metabolism heterogeneity in the regulation of HPA axis activity, a virtual population of subjects was generated and subjected to cluster analysis as well as in silico 11β-HSD1 mRNA knockout (KO). Refer to the Parameter Sampling section of the Supplementary Online Material for more details on modeling subject variability, in silico KO, and cluster analysis. Figure 5 illustrates circadian amplitude scatter plots for the 11β-HSD2 (y axis) and 11β-HSD1 (x axis) enzymes. Each subject is color- and shape-coded by cluster. Four subpopulations are clearly identified for any given season: those with low 11β-HSD1 and 11β-HSD2 (cluster 1, $\overline{x}=0.0646,\overline{y}=28.0250$ ), intermediate 11β-HSD1 and low 11β-HSD2 (cluster 2, $\overline{x}=0.1216,\overline{y}=34.3723$ ), low 11β-HSD1 and high 11β-HSD2 (cluster 3, $\overline{x}=0.0920,\overline{y}=45.7049$ ), and high 11β-HSD1 and 11β-HSD2 (cluster 4, $\overline{x}=0.1702,\overline{y}=47.5001)$ amplitudes. Of the 16 parameters sampled, the distributions of 11 were determined to be significantly different (see the Statistical Analysis section of the Supplementary Online Material). These parameters as well as their median values are shown in Supplementary Table S1 with representative distributions of 11β-HSD1 transcription rate (k8_a), 11β-HSD2 transcription rate (k_10a), and the degradation rate of E (k_12d), arranged by cluster illustrated in Supplementary Figure S4. Lower transcription and translation rates for 11β-HSD1 and 11β-HSD2 mRNAs and proteins (k_8a, k_9a, k_10a, k_11a) were observed in the subjects of cluster 1, a group characterized by low-amplitude 11β-HSD1 and 11β-HSD2 protein oscillations. The highest production rates (transcription and translation) are associated with clusters 3 and 4. Alternatively, degradation rate constants were maximal in cluster 1 when compared with the other 3 clusters, with the exception of the degradation rate of cortisone (k_12d), which was elevated in cluster 4.

Figure 5.

Scatter plots of 741 virtual subjects’ 11β-HSD2 versus 11β-HSD1 enzyme amplitudes for all 4 seasons with each point color- and shape-coded by cluster association. Cluster 1 (black-edged circles) represents those with low 11β-HSD1 and 11β-HSD2 amplitudes, whereas cluster 4 (gray-edged circles) samples exhibit high-amplitude oscillations of both enzymes. Cluster 2 (asterisks) is characterized by intermediate 11β-HSD1 and low-amplitude 11β-HSD2 amplitudes, while cluster 3 (triangles) represents those displaying high 11β-HSD2 but low 11β-HSD1 amplitudes.

Cortisol’s (F_periphery) zenith (maxima) and nadir (minima) values following in silico 11β-HSD1 mRNA KO were compared with their control levels. Figure 6 illustrates scatter plots of the zenith and nadir fold changes of cortisol. Fold-change distributions for each cluster are shown in Supplementary Figures S5 and S6. Fold change was computed as shown below as the ratio of cortisol’s maxima or minima following KO to the control case.

$$\mathrm{Fold\; Change}=\frac{{\mathrm{F}}_{{\mathrm{periphery}}_{\mathrm{Zenith}{\left(\mathrm{or\; Nadir}\right)}_{{\mathrm{HSD}11\mathrm{B}1}_{\mathrm{KO}}}}}}{{\mathrm{F}}_{{\mathrm{periphery}}_{\mathrm{Zenith}{\left(\mathrm{or\; Nadir}\right)}_{\mathrm{Control}}}}}.$$

Figure 6.

Cortisol fold-change scatter plots. Fold change was computed as the ratio of zenith and nadir F_periphery levels following in silico 11β-HSD1 mRNA knockout (KO) to control. The zeniths of the in silico subjects in clusters 2, 3, and 4 were drastically decreased, while those in cluster 1 were relatively unchanged following KO. Alternatively, nadir cortisol levels were generally enhanced in clusters 2, 3, and 4 and decreased in cluster 1. The winter model (A) predicted the greatest increase in cortisol nadir following KO with the summer model (C) being the least affected.

The maxima levels of in silico subjects in cluster 1 were generally unaffected by KO with fold changes (x-axis) close to 1. Alternatively, the maxima of profiles in clusters 2, 3 and 4 were diminished and exhibited fold change values less than 1. Augmented minima values following $11\mathrm{\beta }-\mathrm{HSD}1$ mRNA KO were observed in clusters 2, 3, and 4, whereas diminished levels were predicted for the simulated subjects of cluster 1. Overall, the winter model (Fig. 6A; Suppl. Figs. S5A and S6A) displayed the most marked changes, while the summer model (Fig. 6C; Suppl. Figs. S5C and S6C) exhibited the least. Figure 7 depicts the mean circadian profiles of F_periphery under control (solid lines) and in silico 11β-HSD1 mRNA KO conditions (dashed lines) for the spring model. Simulated subjects of cluster 1 (Fig. 7A) were less affected by 11β-HSD1 mRNA KO versus the other clusters as the amplitude of the rhythms was only slightly decreased. Despite comparable baseline F_periphery oscillations among clusters, the activation of F_periphery by 11β-HSD1 is lowest in cluster 1 (Fig. 8A), highlighting a lower dependency on 11β-HSD1 enzyme activity for the maintenance of normal F_periphery circadian rhythms.

Figure 7.

Comparison of F_periphery circadian profiles of each cluster for the spring model. Profiles shown represent the cluster means for control conditions (solid lines) and following simulated 11β-HSD1 mRNA knockout (dashed lines). The in silico knockout resulted in minimal change of the circadian profile in cluster 1 (A) and pronounced rhythm attenuation in clusters 2 (B), 3 (C), and 4 (D). The gray area represents ±1 standard deviation.

Figure 8.

Cluster variation in F_periphery activation for the spring (baseline) model under normal conditions. Simulated subjects of cluster 1 (A) are characterized by low rates of 11β-HSD1 activity, whereas those of cluster 4 (D) exhibit the highest. These profiles indicate that subjects of cluster 1 are less dependent on 11β-HSD1 enzyme activity for the maintenance of normal F_periphery circadian rhythms. The gray area represents ±1 standard deviation.

Discussion

Numerous seasonal profiles predicted by our model and shown in Figure 2 are in agreement with our previous work (Pierre et al., 2016) and experimental findings. Cortisol levels are generally increased during the winter season (Walker et al., 1997; Persson et al., 2008; Hadlow et al., 2014) with an approximate 5% increase in median cortisol concentration for each hour delay in the time of sunrise (Hadlow et al., 2014). In addition, the phase-delayed F_periphery winter rhythm predicted by our model (Fig. 2B) has been observed experimentally (Vondrasova et al., 1997), and the amplitude increase of the proinflammatory mediator (P) for the winter model (Fig. 2G) is in accordance with the human interleukin (IL)-6 seasonal study (Maes et al., 1994). Increased expression of IL-6 mRNA and sIL-6R during the winter season and peak expression of the BMAL1 encoding gene, ARNTL, during the summer season has been previously observed (Dopico et al., 2015), and so our seasonal circadian profiles of R_p (Fig. 2H) and BMAL1 (Fig. 2D) are in agreement with these findings. Human ACTH (Cauter et al., 1981) and E (Walker et al., 1997) levels are higher in winter, and elevated concentrations of ACTH have been reported in rats following short day exposure (Otsuka et al., 2012), which validates our predictions shown in Figure 2A and 2C. No experimental studies report on the seasonal variation of 11β-HSD1 and 11β-HSD2 protein expression, but it has been theorized that enhanced 11β-HSD1 enzyme activity promotes elevated cortisol levels during the winter (Walker et al., 1997). Given that the 11β-HSD2 mRNA transcription is stimulated by cortisol, the augmentation of the enzyme’s oscillatory amplitude during the winter season (Fig. 2F) seems reasonable. The predicted elevation of 11β-HSD1 protein (Fig. 2E) levels promote increased F_periphery, which drastically reduces cortisol’s nadir during the winter season via negative feedback. This attenuation diminishes the anti-inflammatory influence of F_periphery on P and consequently results in the predicted elevated P rhythms for our winter model. Our results therefore reflect a general shift from an anti- to a proinflammatory state as the seasons change from summer through winter. This seasonal plasticity may play a critical role in the seasonality of diseases such as RA. Both the onset of RA and aggravation of symptoms display peak occurrence during the winter season (Fleming et al., 1976; Hawley et al., 2001).

Figure 3 depicts the seasonal and circadian differences between various mechanisms regulating cortisol’s availability. The seasonal variations of the ACTH-mediated rate of cortisol production versus ACTH levels are plotted in Figure 3A. The differences among seasons are due to the differing adrenal sensitivities to ACTH stimulation, which is heightened during the short photoperiod, winter season (Amirat and Brudieux, 1993; Otsuka et al., 2012). No mechanistic seasonal differences for the enzyme-mediated degradation of F_periphery as a function of substrate concentration were predicted (Fig. 3D). Our model predicted hysteretic behavior for both the 11β-HSD1 and 11β-HSD2 enzymes. Figures 3B and 3C show that while the rate of change of F_periphery activation and inactivation changed with varying E and F_periphery substrate concentrations, different reaction rates were also observed for the same concentrations as well. Mechanisms such as time-varying concentrations of substrate, product, activator, or inhibitor (Purich, 2010), ligand-induced slow isomerization of an enzyme from its inactive to active conformational form, and the slow displacement of an allosteric effector (N.Appaji Rao et al., 1987) have been attributed to modulate enzyme hysteresis. The activity of the phosphofructokinase enzyme is described as hysteretic as it is enhanced by the accumulation of its product, fructose diphosphate (Tornheim and Lowenstein, 1975). Recently, positive feedback loops have been shown to control enzyme hysteresis both experimentally (Kramer and Fussenegger, 2005) and in silico (Novak and Tyson, 2008). For our model, the circadian variation of enzyme and substrate concentrations resulted in the observed responses. The physiological function of enzyme hysteresis is not entirely understood, but the phenomenon has been shown to control sudden metabolite fluxes (Lalanne and Henderson, 1975) and modulate the amplitude of oscillations within pathways (Tornheim and Lowenstein, 1975). This behavior exhibited by the 11β-HSD enzymes may therefore be critical for precisely regulating the circadian bioavailability of F_periphery. The corresponding circadian variations of the rate of change of F_periphery regulated by the various mechanisms are depicted in the lower panels of Figure 3. Overall, a greater dynamic range of oscillations was observed in the winter model, corresponding to the high-amplitude substrate oscillations shown in Figure 2.

It has been previously hypothesized that the seasonal increase in cortisol during the winter season is a consequence of elevated 11β-HSD1 activity and a paradoxical reduction in cortisol secretion rate (Walker et al., 1997). In the study by Walker et al., the tetrahydrocortisol/tetrahydrocortisone ratio was used to quantify 11β-HSD1 activity, and the total excretion of cortisol metabolites was used as an index of cortisol’s secretion rate. Our model predicted slightly lower cortisol clearance (Fig. 4D) but greater adrenal secretion of cortisol in the winter versus the summer (Fig. 4A). Despite the prediction of higher-amplitude 11β-HSD1 protein oscillations in the winter model, 11β-HSD1 activation of F_periphery was elevated in the summer (Fig. 4B, 4C). For our model, this increase in 11β-HSD1 activity for the summer compensates for the reduced sensitivity to ACTH stimulation. The differences between our model’s predictions and the seasonal experimental study may be due to confounding enzymatic mechanisms. While urinary analysis of cortisol metabolites may provide insight into adrenal function and peripheral metabolism, it is important to note that cortisol metabolites are also regulated by other enzymes such as 20β-oxoreductase, 5β-reductase, and 5α-reductase, which may lead to conflicting findings. A recent study found no relationships between HSD11B1 (the 11β-HSD1 encoding gene) genotypes and urinary cortisol/cortisone metabolite ratios in a cohort of women who expressed higher 11β-HSD1 mRNA levels and exhibited elevated 11β-HSD1 activity (Gambineri et al., 2011). Walker et al. (1997) found no seasonal differences in urinary cortisol to cortisone ratios, whereas our model predicted an elevation of this ratio for the winter model (Fig. 4E). While this ratio is frequently used to characterize 11β-HSD1 activity, we again emphasize the involvement of other mechanisms, such as adrenal production, which also regulate F_periphery levels. These ratios, as metrics for 11β-HSD1 activity, should therefore be interpreted with caution.

Interindividual heterogeneity was modeled by sampling parameters describing cortisol metabolism dynamics in the periphery (see the Parameter Sampling section of the Supplementary Online Material). In so doing, we created a virtual population of subjects, with each simulated subject being representative of a unique individual with normal baseline cortisol circadian rhythms. The use of virtual patients is quite common in various phases of drug development as it allows investigators to assess the underlying factors that contribute to differential patient responses to treatment (Rostami-Hodjegan and Tucker, 2007; Teutonico et al., 2015). In this work, the ensemble of in silico subjects allowed us to capture the biological variability existent within a population such that the influences of diverse metabolic dynamics on HPA axis functionality could be evaluated. The seasonal scatter plots of 11β-HSD2 versus 11β-HSD1 enzyme amplitudes reveal both intersubject and seasonal variability (Fig. 5). The in silico subjects of cluster 1 (black-edged circles) are characterized by low-amplitude oscillations for both enzymes, while those in cluster 4 (gray-edged circles) display pronounced amplitude oscillations. Cluster 2’s simulated subjects (asterisks) exhibited intermediate and low-amplitude oscillations, with cluster 3 (triangles) showing low and elevated amplitudes for 11β-HSD1 and 11β-HSD2, respectively. The higher-amplitude oscillations of both enzymes during the winter (Fig. 5A) agree with the representative profiles shown in Figure 2. As shown in Supplementary Figure S4 and Table S1, transcriptional and translational rate constants governing enzyme dynamics varied among the subpopulations. Despite the possession of equivalent phenotypes (comparable circadian F_periphery profiles), the parameters regulating enzyme dynamics were substantially different. The lower production and higher enzyme degradation rate constants of cluster 1 promote lower-amplitude 11β-HSD1 and 11β-HSD2 protein oscillations, while those of cluster 4 result in elevated amplitudes. Normal F_periphery circadian rhythms are maintained as low enzyme levels and activity (cluster 1) are countered by lower cortisone degradation rates (k_12d) and high enzyme oscillations with higher degradation rates (cluster 4). In the case of cluster 1, for instance, a lower degradation rate prolongs the availability of cortisone, which allows lower levels of the 11β-HSD1 enzyme to significantly contribute to cortisol’s profile. Our multivariate multiple linear regression model (Suppl. Table S2) indicated that the relationships between the predictor variables (rate constants) and the enzyme’s amplitudes were cluster dependent, with changes in the transcription and translation rates of 11β-HSD1 predicted to have a greater impact on the amplitude of 11β-HSD1 protein oscillations in cluster 4 versus cluster 1. The differing regression models generated for each cluster underscores the fact that similar phenotypic representations under baseline (unperturbed) conditions are possible despite the variability in parameters regulating enzyme dynamics.

Our in silico 11β-HSD1 mRNA KO study yielded differing F_periphery responses that varied by cluster association and season (Fig. 6; Suppl. Figs. S5 and S6). Suppression of transcription of 11β-HSD1 mRNA resulted in the general elevation of cortisol’s minima (nadir fold change >1) and attenuation of its maxima (zenith fold change <1) for simulated subjects of clusters 2, 3, and 4, while cluster 1 exhibited slight decreases in minima and maxima. The simulated KO resulted in a loss of cortisol regeneration that reduced negative feedback to the central compartment. ACTH production subsequently increased as a compensatory response (Suppl. Fig. S7), which promoted elevation of cortisol’s minima. This response was greatest for the simulated subjects of cluster 4. The greater fold change predicted for the winter model is due to seasonal augmentation of adrenal sensitivity to ACTH. Figure 7 depicts the mean F_periphery profiles for both the control and 11β-HSD1 mRNA KO case. These profiles further echo our previous observations whereby cortisol’s circadian rhythm is only mildly affected for in silico subjects of cluster 1 (Fig. 7A) and vastly attenuated in the others. Those of cluster 1 are less responsive to the loss of 11β-HSD1 expression and activity as their F_periphery dynamics are less dependent on peripheral cortisol activation (Fig. 8), and so the loss of 11β-HSD1 mRNA resulted in minimal change. Differing strain-dependent responses following Hsd11b1 KO have been observed with reports of no change (Carter et al., 2009), elevated basal corticosterone and ACTH levels (Harris et al., 2001), as well as attenuated corticosterone levels (Paterson et al., 2007). In humans, polymorphisms in HSD11B1 have direct functional effects on 11β-HSD1 enzyme activity (Draper et al., 2003; Malavasi et al., 2010) and are associated with increased nadir cortisol levels and depression (Dekker et al., 2012), diabetes, and metabolic syndrome (Gambineri et al., 2011). In select ethnic groups, conflicting associations have been reported between HSD11B1 single nucleotide polymorphisms (SNPs) and diseases such as metabolic syndrome and type 2 diabetes (Devang et al., 2016). Furthermore, the A allele of rs13306421, an allelic variation known to enhance expression and activity of 11β-HSD1, is found in African American and Japanese populations but not in Chinese, European, or Nigerian (Malavasi et al., 2010). Given our model’s predictions, the strain-dependent HPA axis responses in mice Hsd11b1 KO studies as well as varying subpopulation incidence of HSD11B1 SNPs with ethnic group–dependent associations of these SNPs with disease may be a result of genetic variations that alter the dynamics of the 11β-HSDs. These metabolic variations could render select cohorts more susceptible to HPA axis dysfunction following perturbation, potentially leading to the development of a diseased state. As identified in Supplementary Tables S1 and S2, the in silico subjects of cluster 4 exhibit the highest 11β-HSD1 mRNA transcription and protein translation rates, are stimulated the strongest by the proinflammatory mediator, and respond the greatest to changes in production rates and proinflammatory stimulation. These simulated subjects also display the greatest decrease in cortisol maxima following 11β-HSD1 mRNA KO. Based on our analyses, subjects with a cluster 4–type physiology would be more susceptible to systemic perturbation following the presentation of inflammatory triggers but may also be more responsive to therapeutic, 11β-HSD1 inhibitory agents. These agents are currently being developed for the treatment of metabolic syndrome and obesity-related disorders (Anagnostis et al., 2013). Both the local and global sensitivity analyses (see the Sensitivity Analysis section of the Supplementary Online Material for more details on methods) indicated that our model output (cortisone) was most sensitive to changes in parameters regulating the production and degradation of the transcript and protein levels of the enzymes as well as cortisone’s degradation rate (Suppl. Figs. S8-S11). Although there has been a recent surge in development of 11β-HSD1 inhibitors, sufficient knowledge of the inhibitors’ impact on the HPA axis is lacking (Harno and White, 2010). Based on our sensitivity analyses and our in silico KO study, it is apparent that inhibition of 11β-HSD1 activity would indubitably affect the HPA axis by way of altered cortisol dynamics in the periphery, which would subsequently feed back to the central compartment. It is expected that this response will vary in a subject-dependent manner.

Our detailed model explores the emergent properties that arise out of complex interactions between the central and peripheral compartments. It is possible, however, to describe our model in a simplified form. Time-scale analysis (see Time-Scale Analysis section of the Supplementary Online Material) revealed that numerous variables, which shared similar reaction pathways, evolved with comparable time scales (Suppl. Fig. S12). Groups of species acting through the same pathways with comparable time scales share similar dynamics and so can be lumped together without compromising the system’s description (Whitehouse et al., 2004). While the implementation of a formal model reduction technique is outside the scope of this article, our system can potentially be described using a central compartment consisting of 3 core variables: ACTH, CRH, and F_periphery. The nuclear PER-CRY complex, BMAL1, and CLOCK-BMAL1 protein variables could characterize the peripheral clock dynamics, and only the activated FMR(N) and FGR(N) complexes are needed to mediate cortisol signaling in the periphery. Proinflammatory dynamics can be described primarily with P, while the active enzymes 11β-HSD1 and 11β-HSD2 with E (in conjunction with F_periphery) would be sufficient to describe the bioavailability of cortisol in the periphery.

The model presented in this work integrates the effects of photoperiod, HPA axis regulation, and 11β-HSD activity to predict how these influences interact to modulate cortisol dynamics. The photoperiod-induced seasonal predictions were qualitatively similar to experimental data. By sampling parameters regulating 11β-HSD dynamics, we were able to predict and hypothesize about the underlying mechanisms contributing to the experimentally observed, differential HPA axis responses following 11β-HSD1 mRNA KO in mice. Human subpopulation differences in the regulation of 11β-HSD enzymes may alter enzyme functionality and render select groups more susceptible to HPA axis dysfunction following system perturbation. This may lead to the irreparable transition from a normal to abnormal state. Our future work will involve the investigation of how these multitiered regulatory influences of cortisol dynamics respond to, and are affected by, chronic inflammatory triggers.