A mathematical model of HPA axis dynamics and impacts of alcohol consumption

The HPA axis plays a crucial role in regulating an individual’s circadian rhythm, largely driven by cortisol and glucocorticoids. An insufficiency in the concentration of these hormone molecules can result in significant physiological consequences [14]. Typically, cortisol levels peak in the early morning around 7:00 - 8:00 AM, which aids the body in waking up and transitioning out of sleep. Afterward, cortisol levels gradually decline within 2-3 hours following sleep onset, with the lowest concentration, or nadir, occurring around 2:00 - 3:00 AM [15]. Of course, variation between individuals exists. Given this cyclical nature, the circadian rhythm can be modelled using a periodic function. For simplicity, at this point in our model, we will define t = 0 as 8:00 AM and assume the nadir of the circadian rhythm occurs at t = 12, corresponding to 8:00 PM. In Section 4.1.2, we will not assume that the circadian rhythm is symmetric about t = 12 (8:00 PM) and explore an antisymmetric circadian rhythm function with a minium occurring at 2:00 AM [15]. Importantly, the circadian rhythm exerts a positive influence on the HPA axis throughout the day, with varying magnitudes of influence depending on the time and proximity to nadir levels. Thus, our periodic function must always remain positive for any given time t. Based on this rationale, we define the circadian drive function in its most basic form as: 1s(t)=cos(t)+23.s(t) = {{\cos (t) + 2} \over 3}.

For this study, we define the circadian drive function, s(t), as a positive periodic function with t, measured in hours past 8:00 AM. The function is constrained such that min(s(t))=13\min (s(t)) = {1 \over 3} and max(s(t)) = 1. The minimum value occurs around t = 12, corresponding to 8:00 PM, while the maximum value occurs at t = 0, or 8:00 AM. Again, these times are not concrete and may vary between individuals along with day by day differences.

To represent the circadian rhythm over a single 24-hour period, we define s(t) within the domain t ∈ [0,2π], where π is equivalent to 12 hours. Thus, the periodic function s(t) captures the natural rise and fall of cortisol levels over a day, with π marking the mid-point (8:00 PM). The following Table 1 illustrates the correspondence between t values and specific times of day for clarity:

Table 1

t (hours past 8:00 AM)	Time of Day
0	8:00 AM
π2{\pi \over 2}	2:00 PM
π	8:00 PM
3π2{{3\pi } \over 2}\,	2:00 AM (next day)
2π	8:00 AM (next day)

This definition ensures that s(t) effectively models the expected oscillatory behavior of cortisol levels throughout the day, with peak levels at 8:00 AM and the lowest levels around 8:00 PM. The function s(t) will maintain a smooth transition between these points, preserving the physiological characteristics of the circadian rhythm.

As humans age, physiological systems within the body gradually lose their functionality due to the natural process of atrophy [16]. This atrophic decline can be accelerated by various external factors, including environmental conditions, health status, and genetic predispositions, all of which significantly influence an individual’s long-term health. Another external factor that can exacerbate the rate of atrophy is alcohol consumption, which introduces toxins that impair the body’s normal functioning. The specific impact of alcohol and its byproducts on stress response and other systems is discussed in Section 2.3.

Stress response varies significantly across different life stages. For instance, the stress regulation mechanisms in a newborn differ markedly from those in an adult, and there are further distinctions when comparing an adult with an elderly individual. A study by the University of Geneva, showed that older adults (65-84) had less subjective stress than younger adults (21-30) [17]. Consequently, in our model, we assume age 30 to be the “optimal” age for stress response, as we consider a periodic function and want the higher values to result from the domain’s minimum values, which occur between ages 21 and 30. Further, this choice of optimal age represents a balance between the availability of life experiences, emotional stability, and physical youthfulness.

However, there are known hormones that can influence activation of the HPA axis and potentially change this optimal age of stress response depending on the gender of an individual. Estrogen, the main female gonadal hormone, has been shown to increase cortisol levels when its concentration within an individual increases [18]. High cortisol levels in females can inhibit progesterone (P), cortisol modulator, production causing anxiety and menstrual irregularities. On the other hand, high P levels during certain phases of the menstrual cycle can lead to greater free cortisol levels in response to stress [19]. However, if progesterone levels are sufficient, it can adequately modulate cortisol, in acute and chronic stress scenarios, which is a crucial function especially during pregnancy to protect the fetus. Hence, the relationship between cortisol and progesterone is not as dichotomous as estrogen and cortisol’s relationship, but rather characterized by phases of the menstrual cycle or if a female is in a menopausal state. Thus, if we consider estrogen levels in females, an age in the range of 20-29 would be more accurate as an optimal age for stress responses as this is the point in a females life in which estrogen levels are the highest [20]. It is important to note that estrogen is also dependent on the phase of the menstrual cycle and menopause in females. Estrogen is also present in males, of course serving a different function. The primary gonadal hormone in males, testosterone, has effects similar to estrogen, as increased testosterone levels lead to elevated cortisol levels. Typically, the highest levels of testosterone in males occur in the age range of 17-29 [21, 22]. Thus, if we consider testosterone levels in males, an age in the range of 17-29 would be more accurate as an optimal age for stress response. There are even more factors that can influence the choice of an optimal stress response age, independent of gender, such as work conditions, genetics, disease, grief, etc. Thus, without loss of generality, we assume once again that 30 years old is the optimal age for stress response in both males and females.

To incorporate age variation into our model, we employ a function, a(x), that measures the absolute distance from the chosen optimal age. Given that deviations from the optimal age are important, whether younger or older, the absolute difference |x – γ| serves as a suitable metric, where x represents an individual’s age and γ denotes the optimal age. We further refine this metric by normalizing it via the maximum age within the scope of the model, such that the metric lies between 0 and 1 and a(x) attains its maximum value at the optimal age, γ. Finally, we will enclose this normalized metric within a cosine function to ensure that a(x) is symmetric about γ and that no individual, x, receives a value of zero. The proposed function is defined as: 2a(x)=cos(|x−γ|max(D)),a(x) = \cos \left( {{{|x - \gamma |} \over {\max (D)}}} \right), where x ∈ [21, 70] = D is the age of the individual, γ is the “optimal” age for stress response (in our model, γ = 30), and max(D) = 70 is the upper limit of the domain. The upper limit of 70 is chosen to exclude individuals older than 70 years old, as they fall outside the scope of this model. The lower bound of 21 is set based on the legal drinking age in the United States, considering that alcohol consumption plays a role in this study [23]. Thus, a: [21, 70] → [0.8411,1], reaching its maximum value, a(x) = 1, when x = γ and minimum value at x = max (D) = 70, as intended. Note that the function relies on radian calculations to maintain the periodic nature of the cosine function and ensure smooth transitions across the domain.

Our primary focus is to understand how alcohol consumption affects the dynamics of the HPA axis. A common biomarker used to gauge an individual’s alcohol consumption is the blood alcohol concentration (BAC). More specifically, an individual’s BAC over a certain time period can be modeled with the Widmark equation [24] defined as: 3N=σW(B(t)+βt)δz,N = {{\sigma W(B(t) + \beta t)} \over {\delta z}}, where N is the number of standard drinks consumed (for beer, 12 fl. oz with 5% alcohol; for wine, 6.25 fl. oz with 12% alcohol; for liquor, 2 fl. oz with 40% alcohol, Table 2), W is body weight in ounces (oz) (W varies whether we are considering a male or female via the values in Table 2, i.e. W_m and W_f), σ is the volume distribution constant (relating to water distribution in the body in L/kg), B(t) is the BAC in kg/L at time t, β is the alcohol elimination rate in kg/L/hr, t is the time since the first drink in hours, z is the fluid ounces (fl. oz) of alcohol per drink, and 𝒮 is the density of ethanol (0.8 oz per fl. oz) [24]. Traditionally, a standard drink refers to equivalent alcohol content for different types of alcohol (i.e. varying fluid ounces), resulting in the same value z for all types of alcohol. Hence, for this model to study the impact different types of alcohol have on stress response dynamics, we increase the number of fluid ounces for what we define as a standard drink (unique values of z) [25]. Therefore, we constrain z_B < z_W < z_L Table 2. The alcohol elimination rate, β, varies based off of the concentration of alcohol dehydrogenase (ADH) and aldehyde dehydrogenase (ALDH2) within an individual, however, for this model we will fix β [26]. Solving Eq.(3) for B(t), we obtain: 4B(t)=δzNσW−βt.B(t) = {{\delta zN} \over {\sigma W}} - \beta t.

Observe that the BAC decreases linearly at a rate of β with respect to time t. The initial BAC level, B(t_i), can be defined as: 5limt→tiB(t)=δzNσW−βti=B(ti),\mathop {\lim }\limits_{t \to {t_i}} B(t) = {{\delta zN} \over {\sigma W}} - \beta {t_i} = B({t_i}), where t_i = 0 for our model. However, you can choose any initial time t_i within the interval of a non-zero BAC. We can determine the time when the BAC reaches 0kg/L, denoted as t_f, by setting B(t) = 0, yielding: 6tf=δzNσWβ.{t_f} = {{\delta zN} \over {\sigma W\beta }}.

For a single instance of drinking, a day d_j, we define the average BAC over the period [ tij,tfj ]\left[ {t_i^j,t_f^j} \right] as: 7dj=1tfj−tij∫tijtfjBj(t)dt=1tfj−tij∫tijtfj(δzNjσW−βt)dt,{d_j} = {1 \over {t_f^j - t_i^j}}\int_{t_i^j}^{t_f^j} {{B_j}} (t)dt = {1 \over {t_f^j - t_i^j}}\int_{t_i^j}^{t_f^j} {\left( {{{\delta z{N_j}} \over {\sigma W}} - \beta t} \right)} dt, where d_j represents the average BAC during this period (a day), dependent on the number of drinks, N_j. Note the superscript j is an indexing variable corresponding to the j^th day. The interval [ tij,tfj ]\left[ {t_i^j,t_f^j} \right], in its entirety, is the only time range of a nonzero BAC for a day d_j. (For every d_j, there exists only one interval [ tij,tfj ]\left[ {t_i^j,t_f^j} \right], meaning an individual will not consume more alcohol after tijt_i^j, during the day d_j). We define tij\mu = {1 \over A}\sum\limits_{j = 1}^A {{d_j}}, as the time the individual consumed their last drink (BAC is at its maximum). We also restrict d_j ≥ 0, since a negative BAC is not possible, so if d_j < 0 (when t → ∞ in Eq.(4)), we require d_j = 0. For simplicity, the time of the first drink is considered effectively simultaneous with the time of the last drink, as the focus of this study is on post-consumption effects rather than the drinking period itself. We focus on how varying volumes of alcohol intake influence the HPA axis. For our study, we assume that an individual will exclusively consume beer, wine, or liquor on any given day, without mixing different types of alcoholic beverages. However, information on a single day of drinking does not provide enough information about an individual’s alcohol consumption patterns. Thus, we consider the average BAC over multiple days. The average BAC over a time period of A days is defined as: 8μ=1A∑j=1Adj,\matrix{ {{g_1}(\mu ) = {2 \over {1 + {e^{ - \mu }}}},} \cr {{g_2}(\mu ) = {{2\mu } \over {\sqrt {1 + {\mu ^2}} }} + 1,} \cr {{g_3}(\mu ) = \arctan (\mu ) + 1.} \cr } where μ is the average BAC over A days, expressed in kg/L. The relationship between alcohol consumption and cortisol levels can be approximated using a logarithmic model. To capture this relationship, we employ a sigmoid function [3, 27], as the body can only produce cortisol at a finite rate regardless of alcohol intake. The three sigmoid functions that we have constructed for this study to model the alcohol factor are as follows: g1(μ)=21+e−μ,g2(μ)=2μ1+μ2+1,g3(μ)=arctan(μ)+1.{{g'}_1} \ne {{g'}_2} \ne {{g'}_3}

Table 2

Parameter	Value	Description	Source(s)
c¯∞{{\bar c}_\infty }	0.2	Minimal stored baseline CRH	[28–30]
b	0.6	Stored CRH decay rate as a function of cortisol	[28–30]
tc	69.3	CRH biosynthesis timescale	[28–30]
q0	28.0	Maximum release rate of CRH in basal state	[28–30]
I0	1.0	Basal level of the external stimuli	[28–30]
k	2.83	Relates stored CRH to CRH release rate	[28–30]
gc, max	42.0	Maximum auto/paracrine effect of CRH in the pituitary	[28–30]
n	5	Hill coefficient describing the self-up-regulation of CRH	[28–30]
q1−1q_1^{ - 1}	25.0	Circulating CRH conc. at half-maximum self-up-regulation	[28–30]
q2	1.8	Ratio of CRH and cortisol decay rates	[28–30]
p2−1p_2^{ - 1}	0.067	or-complex conc. for half-maximum negative feedback	[28–30]
P3	7.2	Ratio of ACTH and cortisol decay rates	[28–30]
p4	0.05	(or-complex conc.)2 at half-maximum positive feedback on r production	[28–30]
P5	0.11	Basal GR production rate by pituitary	[28–30]
P6	2.9	Ratio of GR and cortisol decay rates	[28–30]
ω	0.045	Frequency of 24hr circadian rhythm	[28–30]
zB	0.6	Fluid ounces of alcohol for 12oz. beer	[25]
zW	0.75	Fluid ounces of alcohol for 6.25oz. wine	[25]
zL	0.8	Fluid ounces of alcohol for 2oz. liquor	[25]
β	0.00015	The alcohol elimination rate in kg/L/hr	[24, 31]
σ	0.6	Volume of distribution (a constant relating the distribution of water in the body in L/kg)	[24]
δ	0.8	The density of ethanol (0.8 oz. per fluid ounce)	[24]
Wf	2732.8	Average weight of a female in oz.	[32]
Wm	3196.8	Average weight of a male in oz.	[32]
A	30	Number of days in the period of concern
BACmax	0.0008	Legal limit of BAC in the United States in kg/L	[33]
λ	(varies) with λ ∈ [0,7]	Number of days per week (7 days) where alcohol was consumed of any quantity and type.

Observe that these are constructed so that g′1≠g′2≠g′3{h_i}(t,x,\mu ) = s(t)\cdota(x)\cdot{g_i}(\mu ), and g₁(0) = g₂(0) = g₃(0) = 1. In Eq.(8), μ is in kg/L. However, when using the functions g_i, we convert μ to grams/L. Each function g_i models the alcohol factor at different rates, and we explore how each affects the behavior of the entire model. Thus, we define our function h as follows: 9hi(t,x,μ)=s(t)·a(x)·gi(μ),s(t):\left[ {0,2\pi } \right] \to \left[ {{1 \over 3},1} \right] where h_i(t, x, μ) > 0 for all t, x, and μ. This condition is satisfied since s(t):[ 0,2π ]→[ 13,1 ]0.8411\;\cdot\;{1 \over 3} = 0.2804, a(x): [21, 70] → [0.8411,1] (since here, max(|21 – γ|, |70 – γ|) = 70 – γ, with γ = 30 for this model, Section 2.2) with the domain representing ages in the scope of this model, and all g_i’s satisfy g_i: [0,1.9516] → [1,2.7799], with 2.7799 inclusive only for g₂. The minimum of the domain for g_i corresponds to no alcohol being consumed in the period ofA = 30 days, whereas the maximum of the domain corresponds to a female consuming 10 drinks (N_j = 10 for 1 ≤ j ≤ A) of liquor (z_L) every day for A = 30 days. Therefore, the absolute minimum of h is 0.8411·13=0.2804{{d{c_s}(t)} \over {dt}} = {{{{\bar c}_\infty } + {e^{ - bu(c)}} - {c_s}} \over {{t_c}}},. It is also important to note that the absolute maximum of any h is 2.7799. Thus, formally, min(h) = 0.2804 and max(h) = 2.7799.

We built our model in spirit of the work in [12] that studied circadian drive and depressive disorders. The model [12] is defined as a system of5 nonlinear ODE’s as follows, 10dcs(t)dt=c¯∞+e−bu(c)−cstc,{{d{c_s}(t)} \over {dt}} = {{{{\bar c}_\infty } + {e^{ - bu(c)}} - {c_s}} \over {{t_c}}}, 11dc(t)dt=q0I(t)(1−e−kcs)+gc,max(q1c)n1+(q1c)n−q2c,{{dc(t)} \over {dt}} = {q_0}I(t)(1 - {e^{ - k{c_s}}}) + {g_{c,max}}{{{{({q_1}c)}^n}} \over {1 + {{({q_1}c)}^n}}} - {q_2}c, 12da(t)dt=c1+p2(ur)−p3a,{{da(t)} \over {dt}} = {c \over {1 + {p_2}(ur)}} - {p_3}a, 13du(t)dt=a−u,{{du(t)} \over {dt}} = a - u, 14dr(t)dt=(ur)2p4+(ur)2+p5−p6r.{{dr(t)} \over {dt}} = {{{{(ur)}^2}} \over {{p_4} + {{(ur)}^2}}} + {p_5} - {p_6}r.

Here, c_s is the concentration of stored CRH (sCRH, basal-state), c is the concentration of released CRH (rCRH), a is the concentration of ACTH, u is the concentration of CORT, and r is the concentration of glucocorticoid receptor complexes (R). The formal definitions and values for all variables can be found in Section 3.1 (Table 2). Observe here that cortisol production only depends on the concentration of ACTH and CORT, the negative feedback loop of the model, since if no ACTH is produced, then dudt<0{{du} \over {dt}} < 0. With the addition of our function h, the difference a – u is modulated. Even if a < u, we can have h · a > u due to alcohol consumption increasing the stress response and, consequently, the ACTH concentration. This further supports the necessity of this modulating function in capturing biologically relevant phenomena.

In this subsection, we define the parameters and their values.

Following Section 2.3, we wish to look at a period of A-days of various drinking patterns. In this model, we will set A = 30 in Eq.(8), although A can be changed to focus on a longer or shorter time period as desired. The patterns (schedules) of drinking we will examine are the college (λ = 4) [34], everyday (λ = 7), weekdays (λ = 5), weekends (λ = 2), and sporadic (λ ∈ [0,7]) drinking schedules where λ is defined in Table 2. Below are the binary schedules (Monday-Sunday) for a 7 day period and is repeated for A-days total. The frequency of drinking, λ, is the main difference schedule to schedule. The sporadic schedule is A-many randomly generated numbers, denoted by α, such that α ∈ {0,1}. (Note: 0=not drinking, 1=drinking on a particular day).

To determine the number of drinks consumed on a day, N_j, a random integer between 1 and 10 was generated with a fixed alcohol type, Z ∈ {z_B, z_W, z_L}. The range 1 to 10 was chosen based off the guidelines that binge drinking is defined as 5 or more drinks for males and 4 or more drinks for females within 2 hours [35]. We will cap the upper limit at 10 drinks. Formally, v_j ∈ {0, 1}, the binary representation if a drink was consumed on a day, d_j, 15N=[ f(v1)f(v2)⋮f(vA) ], where f(vj)={ vj·randint{1,10}if vj=10if vj=0 , for 1≤j≤A.\matrix{ {N = \left[ {\matrix{ {f({v_1})} \cr {f({v_2})} \cr \vdots \cr {f({v_A})} \cr } } \right]{\rm{, where }}f({v_j}) = \left\{ {\matrix{ {{v_j}\cdot{\rm{randint}}\{ 1,10\} } \hfill & {{\rm{if }}{v_j} = 1} \hfill \cr 0 \hfill & {{\rm{if }}{v_j} = 0} \hfill \cr } } \right.,{\rm{ for }}1 \le j \le A.} \hfill \cr }

The N_jth entry of N corresponds to the number of drinks consumed on day d_j. Note, unless otherwise specified, i.e. for our ten drink simulation results, we use the random Python package with these constraints. The random values for said results can be found in Appendix 7. This process is followed in the calculation of d_j in Eq.(7).

A method to include circadian responses is to include oscillations in the sensitivity parameter that regulates cortisol release from stimulation of ACTH [12]. This is due to the adrenal glands response to circadian stimuli through splanchnic nerves which relay synaptic signals from the CNS (Central Nervous System) to the Peripheral SNS’s (Sympathetic Nervous System) neurons. This neural chain and the adrenal medulla allow input from the SCN (Suprachiasmatic Nucleus), the primary structure involved in the circadian rhythm, in the anterior hypothalamus to be received by the adrenal tissues. Further, the paper [12] specifies that the further analysis would include diurnal periodicity, s(t), and age, a(x), which impacts HPA axis response [36–39]. Thus, we considered to model the production of cortisol by a periodic function, h, such that dudt=h·a−u{{du} \over {dt}} = h\,\cdot\,a - u. Thus, our full model is defined as, 16dcs(t)dt=c¯∞+e−bu(c)−cstc,{{d{c_s}(t)} \over {dt}} = {{{{\bar c}_\infty } + {e^{ - bu(c)}} - {c_s}} \over {{t_c}}}, 17dc(t)dt=q0I0(1−e−kcs)+gc,max(q1c)n1+(q1c)n−q2c,{{dc(t)} \over {dt}} = {q_0}{I_0}(1 - {e^{ - k{c_s}}}) + {g_{c,max}}{{{{({q_1}c)}^n}} \over {1 + {{({q_1}c)}^n}}} - {q_2}c, 18da(t)dt=c1+p2(ur)−p3a,{{da(t)} \over {dt}} = {c \over {1 + {p_2}(ur)}} - {p_3}a, 19du(t)dt=h(t,x,μ)a−u,{{du(t)} \over {dt}} = h(t,x,\mu )a - u, 20dr(t)dt=(ur)2p4+(ur)2+p5−p6r,{{dr(t)} \over {dt}} = {{{{(ur)}^2}} \over {{p_4} + {{(ur)}^2}}} + {p_5} - {p_6}r, where u ≡ u_i and h := h(t, x, μ) = h_i(t, x, μ) = s(t) · a(x) · g_i(μ) for i = 1,2,3, as outlined in Eq.(9) depending on which function g we are studying the model with. Note that in Eq.(17), we replace I(t) with I₀, a constant, as our circadian drive will be explored in the form of the rate function, h. In Eq.(16), an important distinction is that c_s(t), stored CRH (sCRH), evolves on a slower timescale hence the division by t_c [28–30]. The equilibrium of c_s is c¯∞+e−bu(c){{\bar c}_\infty } + {e^{ - bu(c)}} and thus we see that it depends on cortisol concentration, u, through the indirect negative feedback loop. Observe that as c_s(t) approaches 0, u(t) tends to 0, and when u(t) goes to infinity, cs(t) levels out at c¯∞+1tc{{{{\bar c}_\infty } + 1} \over {{t_c}}}. In the exponential, b is a modulator of the decay based on cortisol levels.

In Eq.(17), we explore the circulating levels of CRH. The q_2^c term represents the natural degradation of circulating CRH. The concentration of released CRH (rCRH) is driven by the q₀I₀ (1 – e^−kc_s(t)) and gc,max(q1c)n1+(q1c)n{g_{c,max}}{{{{({q_1}c)}^n}} \over {1 + {{({q_1}c)}^n}}} terms. Within these catalytic terms, we have the baseline external stimuli I₀ multiplying into q₀(1 – e^–kc_s(t)), where q₀ is the maximum secretion rate (conversion rate) of stored CRH to circulating CRH occurring when c_s(t) → ∞. When c_s(t) → 0, there is no stored CRH to convert to rCRH. Similar to Eq.(16), the constant k in the exponential function is a modulator of the response between the limits just explored. The second term of this equation models the self-up-regulation property of c(t) which will further increase the concentration of rCRH. To model self-up-regulation, we utilize a Hill-type increasing function [40] which lies between 0 and g_c,max. Within this, q₁ represents the inverse of the rCRH concentration that would result in half-maximum up-regulation.

Next, Eq.(18) models ACTH which is released by the pituitary gland. ACTH will increase with rCRH levels however the rate of increase will be slowed down based upon the concentration of glucocorticoid receptor complexes, hence the ur term in the denominator. This ur term also functions as a negative feedback term due to it depending on later hormonal concentrations in the stress response process. Furthermore, p₂ is similar to q₁ in Eq.(17) as it is the inverse of the ur concentration that would result in half-maximum negative feedback. The –p₃a term models the clearance of ACTH (naturally) and p₃ is the clearance rate.

The cortisol equation, Eq.(19), is a positive feedback equation. For the case of h = 1, if a > u, then u(t) → ∞; if a < u, then u(t) → 0; and u(t) is stable if a = u. Biologically, h ≠ 1 (consistently) and varies due to factors outlined in the construction of h, Section 2. Our model aims to explore how these factors influence stress response dynamics. Cortisol production is inhibited based on the time of day due to the decrease in alertness (dudt<0)\left( {{{du} \over {dt}} < 0} \right) as the day goes on [41]. Alcohol consumption also inhibits normal cortisol production since it results in a cascading stress response as discussed in Section 1. In this case, an individual consuming alcohol (μ → ∞), results in u(t) → ∞. This is since, in that limit, a > > u given a ≥ u > 0 initially. For our model, our timescale is small such that age, x, is fixed for a specific study (i.e. A < 365). However, the model can be easily adapted to consider a larger timescale across many years where x varies (A ≥ 365).

In Eq.(20), r, the glucocorticoid receptor complexes (GR) concentration has a baseline production rate of p₅, which is independent of cortisol concentration, and a clearance rate p₆. The concentration of GR is cortisol dependent and modeled by (ur)2p4+(ur)2{{{{(ur)}^2}} \over {{p_4} + {{(ur)}^2}}}, similar to a Hill-type increasing function, where p₄ is the half-maximum positive feedback relating ur to r production.

At this point, we will follow a model derivation described in [12]. We will simplify the five equation system into two equations in the (c_s,c) space, the reactants of a stress response. We can set the {a, u, r} subsystem to its equilibrium in order to compute solutions to u in Eq.(16). We will keep c_s unchanged, Eq.(16), and simplify c, Eq.(17), with the goal of preserving bistability. Furthermore, the equilibrium solutions for a, u, r are defined as single-valued, positive, real functions for a fixed c. They are uniquely defined and no bistability could emerge at this stage. However, if bistability was to emerge in the (c_s, c) space, it would be reflected in the (a, u, r) space as well. Equilibrium in this space occurs at the intersection of dcsdt=0{{d{c_s}} \over {dt}} = 0 and dcdt=0{{dc} \over {dt}} = 0. Recall that c_s evolves on a slower timescale and thus, later, the c_s-nullcline will be refereed to as the slow nullcline and the the c-nullcline as the fast nullcline.

To begin, we need an explicit expression for u(c) which shows up in the term e^−bu(c) in Eq.(16). As previously mentioned, to do this we will set the {a, u, r} subsystem to equilibrium. Thus, setting dadt=dudt=drdt=0{{da} \over {dt}} = {{du} \over {dt}} = {{dr} \over {dt}} = 0, and performing algebraic manipulation, we obtain, 21[ p2p32+p23p33p4p5+p2p33p5 ]u4−[ 2p2p32hc+2p2p32p5hc−p33p6 ]u3+[ p2p3h2c2+p2p3p5h2c2−p22p32p4p6hc−3p32p6hc ]u2+3p3p6h2c2u−p6h3c3=0.\left[ {{p_2}p_3^2 + p_2^3p_3^3{p_4}{p_5} + {p_2}p_3^3{p_5}} \right]{u^4} - \left[ {2{p_2}p_3^2hc + 2{p_2}p_3^2{p_5}hc - p_3^3{p_6}} \right]{u^3} + \left[ {{p_2}{p_3}{h^2}{c^2} + {p_2}{p_3}{p_5}{h^2}{c^2} - p_2^2p_3^2{p_4}{p_6}hc - 3p_3^2{p_6}hc} \right]{u^2} + 3{p_3}{p_6}{h^2}{c^2}u - {p_6}{h^3}{c^3} = 0.

During this derivation, we see that having h > 0 is necessary for the division involved in the process as well as for this polynomial to maintain shape. Also, for this model x and μ are fixed for all h, hence h is ultimately dependent on time, t, and therefore multiplication of h and c is defined. In Eq.(21), u is the root of the quartic polynomial where h, c are parameters that appear in some of the coefficients. It would be inconvenient to write the full exact solution to Eq.(21) and thus using the realistic parameters from Table 2, we observe that the term, 22[ p2p32+p23p33p4p5+p2p33p5 ]u4+[ 2p2p32hc+2p2p32p5hc−p33p6 ]u3−3p3p6h2c2u,\left[ {{p_2}p_3^2 + p_2^3p_3^3{p_4}{p_5} + {p_2}p_3^3{p_5}} \right]{u^4} + \left[ {2{p_2}p_3^2hc + 2{p_2}p_3^2{p_5}hc - p_3^3{p_6}} \right]{u^3} - 3{p_3}{p_6}{h^2}{c^2}u, is relatively small compared to the other terms. Therefore, we can simplify the relation between u and c approximately by adding Eq.(21) and Eq.(22) and obtain 232[ p2p32+p23p33p4p5+p2p33p5 ]u4+[ p2p3h2c2+p2p3p5h2c2−p22p32p4p6hc−3p32p6hc ]u2−p6h3c3=02\left[ {{p_2}p_3^2 + p_2^3p_3^3{p_4}{p_5} + {p_2}p_3^3{p_5}} \right]{u^4} + \left[ {{p_2}{p_3}{h^2}{c^2} + {p_2}{p_3}{p_5}{h^2}{c^2} - p_2^2p_3^2{p_4}{p_6}hc - 3p_3^2{p_6}hc} \right]{u^2} - {p_6}{h^3}{c^3} = 0 which is a quadratic equation in u². Now, Eq.(23) can be solved for u²(c) and the resulting u(c) can be used in c¯∞+e−bu(c){{\bar c}_\infty } + {e^{ - bu(c)}} so that c_s(t) completely lives in the (c_s, c) space. The solutions for u(c) in the full Eq.(21) and approximated Eq.(23) model at the bounds of h can be seen below in Figure 1.

Fig. 1

The dotted curves represent the solutions to u(c) in the full model, Eq.(21). The solid curve represents the solutions to u(c) in the simplified model, Eq.(23). Parameters from Table 2.

As previously mentioned, with our addition of the function h, we have to make sure the main properties of the system are intact. Mainly, the fact that h > 0 no matter the function g or domain values we use. At the end of Section 2.3, we found that 0.2804 ≤ h ≤ 2.7799, hence h is positive for any function g and domain value. Again, we ensure via Figure 1 that the shape of the solutions to u(c) is maintained at the bounds of h for the full and simplified model.

We observe in Figure 1 that the shape of the roots is maintained in the simplified and full model. Comparing Figure 1(a) and Figure 1(b), we see that the graphs behave similarly with u in Figure 1(b) increasing more rapidly due to the multiplication of h = 2.7799 in u(c). This outcome is expected because as h increases, the cortisol concentration, u, also increases. Consequently, when solving for the roots of u in h · c, the solutions will be proportional to h. Hence, the behavior of the roots of the polynomial u(h · c) is preserved in the extreme cases ofh = 0.2804 and h = 2.7799, as desired, with anticipated differences at those bounds.

Next, we discuss the process used in [12] to form the c(t) equation based off of its dynamics. In the works [28, 29], the bistability is related to a sigmoid function. Thus, we will consider the simplest sigmoid function, a cubic function, f(c). We will continue to fix I(t) = I₀ since we are including circadian drive in the function h. Biologically, we expect that if no sCRH is being produced, then no rCRH can be present. Hence, we need the c-nullcline to pass through the origin. Thus, we need dcdt>0{{dc} \over {dt}} > 0 when c_s > 0 and c = 0 (basal). This lends itself to CRH release being increased when there exists stored CRH. Similarly, dcdt>0{{dc} \over {dt}} < 0 when c_s = 0 and c > 0. Hence, when there is no sCRH, there will be a decrease in rCRH as the concentration of c_s will approach 0 (below baseline). In the previous works [12], the authors posit dcdt=k2I0cs−k1f(c){{dc} \over {dt}} = {k_2}{I_0}{c_s} - {k_1}f(c) where f(c)=c33−c1∗+c2∗2c2+c1∗c2∗cf(c) = {{{c^3}} \over 3} - {{c_1^ * + c_2^ * } \over 2}{c^2} + c_1^ * c_2^ * c is the cubic we desired and the ci∗c_i^ * represent the turning points for f(c), such that df(ci∗)dc=0{{df(c_i^ * )} \over {dc}} = 0. Also, f(c) follows that f(0) = 0 and f(c) > 0 for c > 0. We require that when c_s is positive and c = 0, we have dcdt>0{{dc} \over {dt}} > 0 resulting in a positive value of k₂. To further preserve positivity along the c-nullcline, cs=k1f(c)k2I0{c_s} = {{{k_1}f(c)} \over {{k_2}{I_0}}}, we will have k₁ > 0 to guarantee if c_s = 0 and c > 0, then dcdt<0{{dc} \over {dt}} < 0. Since k₁, k₂ act as proportionality constants, we will also let k₁ << k₂ due to levels of c being larger than c_s. Thus, 24dcdt=k2I0cs−k1(c33−c1∗+c2∗2c2+c1∗c2∗c).\matrix{ {{{dc} \over {dt}} = {k_2}{I_0}{c_s} - {k_1}\left( {{{{c^3}} \over 3} - {{c_1^ * + c_2^ * } \over 2}{c^2} + c_1^ * c_2^ * c} \right).} \hfill \cr }

We let 0<c1*<c2*0 < c_1^* < c_2^* in f(c). When c=ci*c = c_i^*, c_s must be positive, thus we impose the following relationships, c2*>c1*3c_2^* > {{c_1^*} \over 3} and c1*>c2*3c_1^* > {{c_2^*} \over 3} for the turning points of f(c). The full model in the (c_s, c) space can be defined as 25dcsdt=c¯∞+e−bu(c)−cstc{{d{c_s}} \over {dt}} = {{{{\bar c}_\infty } + {e^{ - bu(c)}} - {c_s}} \over {{t_c}}} 26dcdt=k2I0cs−k1(c33−(c1*+c2*)2c2+c1*c2*c).{{dc} \over {dt}} = {k_2}{I_0}{c_s} - {k_1}\left( {{{{c^3}} \over 3} - {{(c_1^* + c_2^*)} \over 2}{c^2} + c_1^*c_2^*c} \right).

Recall that u(c) is determined using Eq.(23) by finding the roots of the quartic polynomial in u²(c). At this stage, we re-parametrized the (c_s, c) model again, following the guidelines of [12]. In doing so, we need to guarantee the intersection of the c_s and c nullclines on multiple occasions in order to maintain bistability. The coordinates of the turning points of the c-nullcline can be found by setting Eq.(24) to equilibrium with ci*c_i^* as input in f(c) and solving for c_s. Therefore, if we define the turning points as (cs,1*,c1*)(c_{s,1}^*,c_1^*) and (cs,2*,c2*)(c_{s,2}^*,c_2^*), then cs,i*=k1f(ci*)k2I0c_{s,i}^* = {{{k_1}f(c_i^*)} \over {{k_2}{I_0}}} by the above logic, where c1*<c2*c_1^* < c_2^*. Further, we know that c decreases along the slow nullcline (c_s-nullcline) as c_s → ∞. Therefore, multiple intersections of the nullclines will occur, if the following conditions hold: 27cs,1*>c¯∞+e−bu(c1*)c_{s,1}^* > {{\bar c}_\infty } + {e^{ - bu(c_1^*)}} 28cs,2*<c¯∞+e−bu(c2*).c_{s,2}^* < {{\bar c}_\infty } + {e^{ - bu(c_2^*)}}.

In other words, if (cs,1*,c1*)(c_{s,1}^*,c_1^*) is to the right of the c_s-nullcline and (cs,2*,c2*)(c_{s,2}^*,c_2^*) is to the left, then bistability will emerge (multiple intersections of the two nullclines). Due to the fact that c is an order of magnitude large than c_s, we rescale c_s by κ=k2k1>>1\kappa = {{{k_2}} \over {{k_1}}} > > 1 which will increase its magnitude in order for comparability. Also, we define the following as part of re-parametrization, 29x=csx = {c_s} 30y=cκ1/3y = {c \over {{\kappa ^{1/3}}}} 31t′=tck11/3k22/3.t' = {t_c}k_1^{1/3}k_2^{2/3}.

Since the rescaling will impact the timescale in which c_s evolves, we will disregard the prime notation on t and define ε=1t′=1(tck11/3k22/3)\varepsilon = {1 \over {t'}} = {1 \over {\left( {{t_c}k_1^{1/3}k_2^{2/3}} \right)}}. We also will assume that tc>>1k1k223{t_c} > > \root 3 \of {{1 \over {{k_1}k_2^2}}} which results in ε << 1. Therefore, the fully reduced model follows by making these substitutions into Eq.(25) and Eq.(26) to obtain 32dxdt=ε(c¯∞+e−bu(y)−x){{dx} \over {dt}} = \varepsilon ({{\bar c}_\infty } + {e^{ - bu(y)}} - x) 33dydt=I0x−f(y),{{dy} \over {dt}} = {I_0}x - f(y), where f(y)=(y33−(y1*+y2*)2y2+y1*y2*y)f(y) = \left( {{{{y^3}} \over 3} - {{(y_1^* + y_2^*)} \over 2}{y^2} + y_1^*y_2^*y} \right), yi*y_i^* are the rescaled values ci*c_i^* via Eq.(30), and u(y) is the solution to: 342[ p2p32+p23p33p4p5+p2p33p5 ]u4+ [ p2p3κ2/3h2y2+p2p3p5κ2/3h2y2 −p22p32p4p6κ1/3hy−3p32p6κ1/3hy ]u2−p6κh3y3=0,\matrix{ {2\left[ {{p_2}p_3^2 + p_2^3p_3^3{p_4}{p_5} + {p_2}p_3^3{p_5}} \right]{u^4}} \hfill \cr { + \left[ {{p_2}{p_3}{\kappa ^{2/3}}{h^2}{y^2} + {p_2}{p_3}{p_5}{\kappa ^{2/3}}{h^2}{y^2}} \right.} \hfill \cr {\left. { - p_2^2p_3^2{p_4}{p_6}{\kappa ^{1/3}}hy - 3p_3^2{p_6}{\kappa ^{1/3}}hy} \right]{u^2}} \hfill \cr { - {p_6}\kappa {h^3}{y^3} = 0,} \hfill \cr } where y is a positive real number. Furthermore, utilizing the inequalities in Eq.(27) and Eq.(28), we constrain yi*>0y_i^* > 0 as follows 35y1*>y2*3>0x1*=3y1*2y2*−y1*36I0>c¯∞+e−bu(y1*)\matrix{ {y_1^* > {{y_2^*} \over 3} > 0} \hfill & {x_1^* = {{3y_1^{{*^2}}y_2^* - y_1^{{*^3}}} \over {6{I_0}}} > {{\bar c}_\infty } + {e^{ - bu(y_1^*)}}} \hfill \cr } 36y2*>y1*3>0x2*=3y2*2y1*−y2*36I0<c¯∞+e−bu(y2*),\matrix{ {y_2^* > {{y_1^*} \over 3} > 0} \hfill & {x_2^* = {{3y_2^{{*^2}}y_1^* - y_2^{{*^3}}} \over {6{I_0}}} < {{\bar c}_\infty } + {e^{ - bu(y_2^*)}},} \hfill \cr } where h shifts the constraints in Eq.(35) and Eq.(36), through the solutions to u(yi*)u(y_i^*), i = 1,2. The simplified model in the (x, y) space resembles the neuron spiking model, the FitzHugh-Nagumo model [13], which has increased reactions to external differences. This similarity is important as the function h will lends itself to those external differences that we wish to see within the model.

Our model can be very informative for alcohol intervention clinicians given its role as a visual representation of the impact that alcohol consumption has on key parts of a patient’s life, such as sleep and high/low stress. This model can be fitted to any individual, provided that you have basic medical information, along with drinking habits and patterns. The model can also be used predictively to plan potential detoxification strategies and show the patient what their future circadian patterns and cortisol levels could be with treatment. Of course, our model displays the dangers of alcohol consumption in high volumes and frequency, thus in some intervention strategies, it might be useful to highlight the potential medical consequences and dangers. Such consequences can be accelerated aging (consistently high levels of cortisol cause high rates of telomere shortening), Cushing’s Syndrome, liver cirrhosis, inter alia, based on stress hormone levels outputted by this model.

Our model employs a simplified circadian drive function, s(t), up to Section 4.1. Improving its complexity and practicality in earlier sections of this study is a key objective. The function s(t) throughout the paper assumes uniform light exposure and no deviations from the sleep schedule. To enhance s(t), one could focus on incorporating the pineal gland’s light-dependent melatonin regulation [49]. When the retina receives light and dark inputs, it sends the feedback to the optic nerve of the brain, which then sends the signal to the SCN. At this point, the SCN will signal the pineal gland to produce or not produce melatonin. In conjunction with cortisol secretion, the pineal and adrenal glands control the circadian rhythm. Future models can explore a system of equations including melatonin production as one equation modeled by pineal gland activity. Including the pineal gland as another system/subsystem of this model would lead one to study the interplay between cortisol and melatonin in the presence of alcohol consumption. Considering varying amounts of sleep per day would also enhance the model. Additionally, in our model, there are many simplifying assumptions regarding age despite biological variation. We assume that our age factor, a(x), is solely dependent on an individual’s age and the normalized absolute distance from an optimal age. Considering different optimal ages for male and female, activity level, and health of an individual is of interest since these are important factors that determine the rate of atrophy for bodily systems in an individual. Including variations in γ related to BMI / weight could also potentially enhance the function.

Currently, our drinking schedules vary the frequency of drinking and do not represent the drinking tendencies for each schedule. A further improvement for our model would be to consider more accurate drinking amounts for each schedule along with varying the time of day alcohol is consumed for each schedule. The social context of alcohol consumption is significant; for example, in college, drinking typically occurs at parties or gatherings that can occur at late night hours. Thus, the college drinking schedule may exhibit a correlation with increased alcohol consumption (larger N_j and longer drinking intervals) due to peer pressure and social influences [50]. Whereas someone drinking everyday might have less to drink and is probably not in a social setting each time. Data supporting various drinking distributions such as probability density functions for the schedules would be of interest. Furthermore, incorporating a step function N(t) to describe the number of drinks consumed at time t would eliminate the assumption of instantaneous alcohol consumption that we observe in the top plots of Figure 5, Figure 6, Figure 7, and Figure 8. In addition, there is a statistical error that is not accounted for in the calculation of B(t), BAC, and accounting for this would further enhance the results [24]. The omission of varying weights and different values of β depending on age and alcohol consumption history limits the model as well. This can be refined through further simulation as well as inclusion of stochastic and statistical error modeling.

This model builds on previous work [12] and explores the impacts of varying degrees of alcohol consumption on the HPA axis. Furthermore, we included an age factor to explore how cortisol production within the HPA axis changes with age. The constructed functions, g, relate alcohol consumption to cortisol production via a logarithmic relationship. Hence, g models the biological limit for basal sCRH and CRH production, inhibiting a constant increase in cortisol. For any function g, there will not be any increase when you approach the case of μ→μ′i\mu \to {{\mu '}_i}, where μ′i≤10{{\mu '}_i} \le 10 is some value of reunitized BAC for which dgidμ|μ=μ′i=0{{d{g_i}} \over {d\mu }}{|_{\mu = {{\mu '}_i}}} = 0. In other words, for any function g, the tangent line at g(μ′i)g({{\mu '}_i}) is the horizontal asymptote for each logarithmic function; each function g has a unique asymptote (i.e. for our model, g1(μ′i)<g3(μ′3)<g2(μ′2){g_1}({{\mu '}_i}) < {g_3}({{\mu '}_3}) < {g_2}({{\mu '}_2})). Furthermore, we explore three different functions, g_i, where dg1dμ≠dg2dμ≠dg3dμ{{d{g_1}} \over {d\mu }} \ne {{d{g_2}} \over {d\mu }} \ne {{d{g_3}} \over {d\mu }}. However, these functions are not experimentally fitted, and thus potentially have unaccounted error to actual biological results. These functions do however follow general trends of the expected behavior for cortisol production driven by different levels of alcohol consumption. Substituting B(t) for μ in the functions g, would allow for a study of instantaneous results.

Our model explored the impacts alcohol consumption has on HPA axis function by looking at BAC via the Widmark equation [24]. An increase in alcohol consumption resulted in a decrease of bistability within the HPA axis Figure 2 and lower availability of sCRH for stress responses Figure 3. We were able to explore impacts on cortisol levels in the function h depending on long-term alcohol consumption habits through various drinking schedules/frequencies. In Figure 4, we see that basal and peak levels of cortisol are vastly higher when μ is larger. This increase becomes taxing on the adrenals and other bodily systems which can result in higher blood pressure and even hypertension due to cortisol’s role as a vasoconstrictor. Other medical complications can develop as previously discussed. Our model’s circadian shift analysis, with an enhanced s(t) function, yields promising results. Various φ¯(t˜)\bar \varphi (\tilde t) values, due to an individual’s alcohol consumption habits, show disruptions to sleep/wake cycles, Figure 5, Figure 6, Figure 7, and Figure 8. Considering both Eq.(37) and Eq.(38) integrated into the main model (Figure 2, Figure 3, Figure 4) would yield interesting results. One would have to check that the solutions to the quartic polynomial, u(c), are still defined at the limits of the piecewise form, s. Combining results from the subplots and circadian analysis, while having g(φ¯(t˜))g\left( {\bar \varphi (\tilde t)} \right) in the function h along with the shifted circadian function, s(t−φ¯(t˜))s\left( {t - \bar \varphi (\tilde t)} \right), could be studied. This function h would be defined as h(t,t˜,x)=s(t−ϕ¯(t˜))·a(x)·g(ϕ(t˜))h(t,\tilde t,x) = s(t - \bar \phi (\tilde t))\cdota(x)\cdotg(\phi (\tilde t)). Also, the time it takes for negative feedback to initiate in the presence of alcohol is of interest. The time is most likely proportional to the time it takes to reach a BAC of 0kg/L in an individual via Eq.(6) in the presence of alcohol, but biological tests would need to be completed. Lastly, our model is easily adapted to other stress catalysts such as caffeine and nicotine, where both substances are known to drive HPA axis activity [51, 52].