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) = \frac{\cos(t)+2}{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)) = \frac{1}{3} and max(s(t))=1\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:

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−γ||x-\gamma| serves as a suitable metric, where x represents an individualߣs age and y 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 y 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(\frac{|x-\gamma|}{\max(D)}\right),

where x ∈ [21, 70] = D is the age of the individual, y 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 = \frac{\sigma W(B(t)+\beta t)}{\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_LTable 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) = \frac{\delta z N}{\sigma W} - \beta t.

Observe that the BAC decreases linearly at a rate of β with respect to time t. The initial BAC level, B(ti)B(t_i), can be defined as:

5limt→tiB(t)=δzNσW−βti=B(ti), \lim_{t \to t_i} B(t) = \frac{\delta z N}{\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 = \frac{\delta z N}{\sigma W \beta}

For a single instance of drinking, a day d_j, we define the average BAC over the period [tij,tfj][t_i^j, t_f^j] as:

7dj=1tfj−tij∫tijtfjBj(t)dt=1tfj−tij∫tijtfj(δzNjσW−βt)dt, d_j = \frac{1}{t_f^j - t_i^j} \int_{t_i^j}^{t_f^j} B_j(t) dt = \frac{1}{t_f^j - t_i^j} \int_{t_i^j}^{t_f^j} \left( \frac{\delta z N_j}{\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][t_i^j, t_f^j] 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][t_i^j, t_f^j] meaning an individual will not consume more alcohol after, during the day d_j). We define tijt_i^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, \mu = \frac{1}{A} \sum_{j=1}^A d_j

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(\mu) = \frac{2}{1+e^{-\mu}} g_2(\mu) = \frac{2\mu}{\sqrt{1+\mu^2}} + 1, g_3(\mu) = \arctan(\mu) + 1.

Observe that these are constructed so that g1′≠g2′≠g3′g_1' \ne g_2' \ne g_3' and g1(0)=g2(0)=g3(0)=1g_1(0) = g_2(0) = g_3(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(μ), h_i(t,x,\mu) = s(t) \cdot a(x) \cdot g_i(\mu)

where hi(t,x,μ)>0h_i(t,x,\mu) > 0 for all t, x, and μ. This condition is satisfied since s(t):[0,2π]→[13,1]s(t): [0, 2\pi] \to [\frac{1}{3}, 1], a(x):[21,70]→[0.8411,1]a(x): [21, 70] \to [0.8411, 1] (since here, max(|21−γ|,|70−γ|)=70−γ\max(|21-\gamma|, |70-\gamma|) = 70-\gamma, with γ = 30 for this model, Section 2.2) with the domain representing ages in the scope of this model, and all g_i satisfy gi:[0,1.9516]→[1,2.7799]g_i: [0, 1.9516] \to [1, 2.7799], with 2.7799 inclusive only for g₂. The minimum of the domain for gi corresponds to no alcohol being consumed in the period of A = 30 days, whereas the maximum of the domain corresponds to a female consuming 10 drinks ((Nj=10 for 1≤j≤A)N_j = 10 of liquor (zLz_L) every day for A = 30 days. Therefore, the absolute minimum of h is 0.8411⋅13=0.28040.8411 \cdot \frac{1}{3} = 0.2804. It is also important to note that the absolute maximum of any h is 2.7799. Thus, formally, min(h)=0.2804\min(h) = 0.2804 and max(h)=2.7799\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 of 5 nonlinear ODE’s as follows,

10dcs(t)dt=c˜∞+e−bu(c)−cstc \frac{dc_s(t)}{dt} = \frac{\tilde{c}_\infty + e^{-bu(c)} - c_s}{t_c}

11dc(t)dt=q0I(t)(1−e−kcs)+gc,max(q1c)n1+(q1c)n−q2c, \frac{dc(t)}{dt} = q_0 I(t) (1 - e^{-kc_s}) + g_{c,\max} \frac{(q_1 c)^n}{1 + (q_1 c)^n} - q_2 c,

12da(t)dt=c1+p2(ur)−p3a, \frac{da(t)}{dt} = \frac{c}{1 + p_2(ur)} - p_3 a,

13du(t)dt=a−u, \frac{du(t)}{dt} = a - u

14dr(t)dt=(ur)2p4+(ur)2+p5−p6r. \frac{dr(t)}{dt} = \frac{(ur)^2}{p_4 + (ur)^2} + p_5 - p_6 r.

Here, csc_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\frac{du}{dt} < 0. With the addition of our function h, the difference a - u is modulated. Even if a<ua < u, we can have h⋅a>uh \cdot 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.

Table 2

Parameter	Value	Description	Source(s)
c¯∞\overline{c}_{\infty}	0.2	Minimal stored baseline CRH	[28-30]
b	0.6	Stored CRH decay rate as a function of cortisol	[28-30]
tct_c	69.3	CRH biosynthesis timescale	[28-30]
q0q_0	28.0	Maximum release rate of CRH in basal state	[28-30]
I0I_0	1.0	Basal level of the external stimuli	[28-30]
k	2.83	Relates stored CRH to CRH release rate	[28-30]
gc,maxg_{c,\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]
q2q_2	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]
p3p_3	7.2	Ratio of ACTH and cortisol decay rates	[28-30]
p4p_4	0.05	(or-complex conc.)2^2 at half-maximum positive feedback on r production	[28-30]
p5p_5	0.11	Basal GR production rate by pituitary	[28-30]
p6p_6	2.9	Ratio of GR and cortisol decay rates	[28-30]
ω\omega	0.045	Frequency of 24hr circadian rhythm	[28-30]
zBz_B	0.6	Fluid ounces of alcohol for 12oz. beer	[25]
zWz_W	0.75	Fluid ounces of alcohol for 6.25oz. wine	[25]
zLz_L	0.8	Fluid ounces of alcohol for 2oz. liquor	[25]
β\beta	0.00015	The alcohol elimination rate in kg/L/hr	[24,31]
σ\sigma	0.6	Volume of distribution (a constant relating the distribution of water in the body in L/kg)	[24]
δ\delta	0.8	The density of ethanol (0.8 oz. per fluid ounce)	[24]
WfW_f	2732.8	Average weight of a female in oz.	[32]
WmW_m	3196.8	Average weight of a male in oz.	[32]
A	30	Number of days in the period of concern
BACmaxBAC_{max}	0.0008	Legal limit of BAC in the United States in kg/L	[33]
λ\lambda	(varies) with λ∈[0,7]\lambda \in [0,7]	Number of days per week (7 days) where alcohol was consumed of any quantity and type.

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=30A=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\lambda=4) [34], everyday (λ=7\lambda=7), weekdays (λ=5\lambda=5), weekends (λ=2\lambda=2), and sporadic (λ∈[0,7]\lambda \in [0,7]) drinking schedules where λ\lambda 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, λ\lambda, is the main difference schedule to schedule. The sporadic schedule is A-many randomly generated numbers, denoted by α\alpha, such that α∈{0,1}\alpha \in \{0,1\}. (Note: 0=not drinking, 1=drinking on a particular day).

College	Everyday	Weekdays	Weekend	Sporadic
0001111	1111111	1111100	0000011	ααααααα\alpha\alpha\alpha\alpha\alpha\alpha\alpha

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∈{zB,zW,zL}z \in \{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, vj∈{0,1}v_j \in \{0,1\}, the binary representation if a drink was consumed on a day, d_j,

15N=[f(v1)f(v2)⋮f(vA)],f(vj)={randint{1,10}if vj=10if vj=0for1≤j≤A. N = \begin{bmatrix} f(v_1) \\ f(v_2) \\ \vdots \\ f(v_A) \end{bmatrix}, \quad f(v_j) = \begin{cases} \text{randint}\{1,10\} & \text{if } v_j = 1 \\ 0 & \text{if } v_j = 0 \end{cases}

The N_j th 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)s(t), and age, a(x)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\frac{du}{dt}=h\cdot a-u. Thus, our full model is defined as,

16dcs(t)dt=c¯∞+e−bu(c)−cstc, \frac{dc_s(t)}{dt} = \frac{\tilde{c}_\infty + e^{-bu(c)} - c_s}{t_c}

17dc(t)dt=q0I0(1−e−kcs)+gc,max(q1c)n1+(q1c)n−q2c, \frac{dc(t)}{dt} = q_0 I_0 (1 - e^{-kc_s}) + g_{c,\max} \frac{(q_1 c)^n}{1 + (q_1 c)^n} - q_2 c,

18da(t)dt=c1+p2(ur)−p3a, \frac{da(t)}{dt} = \frac{c}{1 + p_2(ur)} - p_3 a,

19du(t)dt=h(t,x,μ)a−u, \frac{du(t)}{dt} = h(t, x, \mu)a - u,

20dr(t)dt=(ur)2p4+(ur)2+p5−p6r, \frac{dr(t)}{dt} = \frac{(ur)^2}{p_4 + (ur)^2} + p_5 - p_6 r,

where u≡uiu \equiv u_i and h:=h(t,x,μ)≡hi(t,x,μ)=s(t)⋅a(x)⋅gi(μ)h := h(t,x,\mu) \equiv h_i(t,x,\mu) = s(t) \cdot a(x) \cdot g_i(\mu) for i=1,2,3i=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)I(t) with I0I_0, 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 cs(t)c_s(t), stored CRH (SCRH), evolves on a slower timescale hence the division by tct_c [28-30]. The equilibrium of csc_s is c¯∞+e−bu(c)\overline{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 cs(t)c_s(t) approaches 0, u(t)u(t) tends to 0, and when u(t)u(t) goes to infinity, cs(t)c_s(t) levels out at c¯∞+1tc\frac{\overline{c}_{\infty}+1}{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 q2cq_2 c term represents the natural degradation of circulating CRH. The concentration of released CRH (rCRH) is driven by the q0I0(1−e−kcs(t))q_0 I_0 (1 - e^{-kc_s(t)}) and gc,max(q1c)n1+(q1c)ng_{c,\max} \frac{(q_1 c)^n}{1 + (q_1 c)^n} terms. Within these catalytic terms, we have the baseline external stimuli I0I_0 multiplying into q0(1−e−kcs(t))q_0(1-e^{-kc_s(t)}), where q0q_0 is the maximum secretion rate (conversion rate) of stored CRH to circulating CRH occurring when cs(t)→∞c_s(t) \to \infty When cs(t)→0c_s(t) \to 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)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 gc,maxg_{c,\max} Within this, q1q_1represents 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, p2p_2 is similar to q1q_1 in Eq.(17) as it is the inverse of the ur concentration that would result in half-maximum negative feedback. The -p3a term models the clearance of ACTH (naturally) and p3p_3 is the clearance rate.

The cortisol equation, Eq. (19), is a positive feedback equation. For the case of h=1h=1 if a<ua<u, then u(t)→∞u(t) \to \infty; if a<ua<u, then u(t)→0u(t) \to 0 and u(t)u(t) is stable if a=ua=u. Biologically, h≠1h \ne 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\frac{du}{dt} < 0) 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 (μ→∞\mu \to \infty) results in u(t)→∞u(t) \to \infty. This is since, in that limit, a>>ua >> u given a>u>0a > u > 0 initially. For our model, our timescale is small such that age, x, is fixed for a specific study (i.e. A<365A < 365). However, the model can be easily adapted to consider a larger timescale across many years where x varies (A≥365A \ge 365).

In Eq.(20), r, the glucocorticoid receptor complexes (GR) concentration has a baseline production rate of p5p_5, which is independent of cortisol concentration, and a clearance rate p6p_6. The concentration of GR is cortisol dependent and modeled by (ur)2p4+(ur)2\frac{(ur)^2}{p_4 + (ur)^2} , similar to a Hill-type increasing function, where p4p_4 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 ((cs,c)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 csc_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 ((cs,c)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\frac{dc_s}{dt}=0 and dcdt=0\frac{dc}{dt}=0 Recall that csc_s evolves on a slower timescale and thus, later, the csc_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)u(c) which shows up in the term e−bu(c)e^{-bu(c)} previously mentioned, to do this we will set the {a,u,r} subsystem to equilibrium. Thus, setting dadt=dudt=drdt=0\frac{da}{dt}=\frac{du}{dt}=\frac{dr}{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. [p_2 p_3^2 + p_2^3 p_3^3 p_4 p_5 + p_2 p_3^3 p_5]u^4 - [2p_2 p_3^2 hc + 2p_2 p_3^2 p_5 hc - p_3^3 p_6]u^3 + [p_2 p_3 h^2 c^2 + p_2 p_3 p_5 h^2 c^2 - p_2^2 p_3^2 p_4 p_6 hc - 3 p_3^2 p_6 hc]u^2 + 3p_3 p_6 h^2 c^2 u - p_6 h^3 c^3 = 0.

During this derivation, we see that having h>0h>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 e 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, [p_2 p_3^2 + p_2^3 p_3^3 p_4 p_5 + p_2 p_3^3 p_5]u^4 + [2p_2 p_3^2 hc + 2p_2 p_3^2 p_5 hc - p_3^3 p_6]u^3 - 3p_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]u2−p6h3c3=0 2[p_2 p_3^2 + p_2^3 p_3^3 p_4 p_5 + p_2 p_3^3 p_5]u^4 + [p_2 p_3 h^2 c^2 + p_2 p_3 p_5 h^2 c^2 - p_2^2 p_3^2 p_4 p_6 hc]u^2 - p_6 h^3 c^3 = 0

which is a quadratic equation in u2u^2. Now, Eq.(23) can be solved for u2(c)u^2(c) and the resulting u(c)u(c) can be used in c¯∞+e−bu(c)\overline{c}_{\infty}+e^{-bu(c)} so that cs(t)c_s(t) completely lives in the (cs,c)(c_s, c) space. The solutions for u(c)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

T 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>0h>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.77990.2804 \le h \le 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)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.7799h=2.7799 in u(c)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⋅ch \cdot c, the solutions will be proportional to h. Hence, the behavior of the roots of the polynomial u(h⋅c)u(h \cdot c) is preserved in the extreme cases of h=0.2804h=0.2804 and h=2.7799h=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)f(c). We will continue to fix I(t)=I0I(t)=I_0 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\frac{dc}{dt}>0 when cs>0c_s>0 and c=0c=0 (basal). This lends itself to CRH release being increased when there exists stored CRH. Similarly, dcdt<0\frac{dc}{dt}<0 when cs=0c_s=0 and c>0c>0. Hence, when there is no sCRH, there will be a decrease in rCRH as the concentration of csc_s will approach 0 (below baseline). In the previous works [12], the authors posit dcdt=k2I0cs−k1f(c)\frac{dc}{dt}=k_2 I_0 c_s - k_1 f(c) where f(c)=c33−c1*+c2*2c2+c1*c2*cf(c)=\frac{c^3}{3}-\frac{c_1^*+c_2^*}{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)f(c), such that df(ci*)dc=0\frac{df(c_i^*)}{dc}=0 Also, f(c)f(c) follows that f(0)=0f(0)=0 and f(c)>0f(c)>0 for c>0c>0 We require that when csc_s is positive and c=0c=0, we have dcdt>0\frac{dc}{dt}>0 resulting in a positive value of k2k_2. To further preserve positivity along the c-nullcline, cs=k1f(c)k2I0c_s = \frac{k_1 f(c)}{k_2 I_0}, we will have k1>0k_1 > 0 to guarantee if cs=0c_s=0 and c>0c>0, then dcdt<0\frac{dc}{dt}<0. Since k1,k2k_1, act as proportionality constants, we will also let k1<<k2k_1 << k_2 due to levels of c being larger than csc_s. Thus,

24dcdt=k2I0cs−k1(c33−c1*+c2*2c2+c1*c2*c). \frac{dc}{dt} = k_2 I_0 c_s - k_1 \left( \frac{c^3}{3} - \frac{c_1^* + c_2^*}{2} c^2 + c_1^* c_2^* c \right) .

We let 0<c1*<c2*0<c_1^*<c_2^* in f(c)f(c). When c=ci*csc=c_i^*, must be positive, thus we impose the following relationships, c2*>c1*3c_2^* > \frac{c_1^*}{3} and c1*>c2*3c_1^* > \frac{c_2^*}{3} for the turning points of f(c)f(c). The full model in the (cs,c)(c_s, c) space can be defined as

25dcsdt=c¯∞+e−bu(c)−cstc \frac{dc_s}{dt} = \frac{\tilde{c}_\infty + e^{-bu(c)} - c_s}{t_c}

26dcdt=k2I0cs−k1(c33−(c1*+c2*)2c2+c1*c2*c). \frac{dc}{dt} = k_2 I_0 c_s - k_1 \left( \frac{c^3}{3} - \frac{(c_1^* + c_2^*)}{2} c^2 + c_1^* c_2^* c \right) .

Recall that u(c)u(c) is determined using Eq.(23) by finding the roots of the quartic polynomial in u2(c)u^2(c). At this stage, we re-parametrized the (cs,c) model again, following the guidelines of [12]. In doing so, we need to guarantee the intersection of the cSc_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)f(c) and solving for csc_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}^* = \frac{k_1 f(c_i^*)}{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 (csc_s nullcline) as cs→∞c_s \to \infty Therefore, multiple intersections of the nullclines will occur, if the following conditions hold:

27cs,1*>c¯∞+e−bu(c1*) c_{s,1}^* > \overline{c}_\infty + e^{-bu(c_1^*)}

28cs,2*<c¯∞+e−bu(c2*) c_{s,2}^* < \overline{c}_\infty + e^{-bu(c_2^*)}

In other words, if (cs,1*,c1*)(c_{s,1}^*, c_1^*) is to the right of the csc_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 csc_s, we rescale csc_s by κ=k2k1>>1\kappa = \frac{k_2}{k_1} >> 1 which will increase its magnitude in order for comparability. Also, we define the following as part of re-parametrization,

29x=cs x = c_s

30y=cκ1/3 y = \frac{c}{\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 csc_s evolves, we will disregard the prime notation on t and define ε=1t′=1(tck11/3k22/3)\epsilon = \frac{1}{t'} = \frac{1}{(t_c k_1^{1/3} k_2^{2/3})} We also will assume that tc>>1k1k223t_c >> \sqrt[3]{\frac{1}{k_1 k_2^2}} which results in ε<<1\epsilon << 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) \frac{dx}{dt} = \epsilon (\overline{c}_\infty + e^{-bu(y)} - x)

33dydt=I0x−f(y), \frac{dy}{dt} = I_0 x - f(y) ,

where

f(y)=(y33−(y1*+y2*)2y2+y1*y2*y) f(y) = \left( \frac{y^3}{3} - \frac{(y_1^* + y_2^*)}{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)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, 2[p_2 p_3^2 + p_2^3 p_3^3 p_4 p_5 + p_2 p_3^3 p_5]u^4 + [p_2 p_3 \kappa^{2/3} h^2 y^2 + p_2 p_3 p_5 \kappa^{2/3} h^2 y^2 - p_2^2 p_3^2 p_4 p_6 \kappa^{1/3} hy - 3p_3^2 p_6 \kappa^{1/3} hy]u^2 - p_6 \kappa h^3 y^3 = 0,

where yis apositive real number. Furthermore, utilizing the inequalities in Eq.(27) and Eq.(28), we constrain y_i^* > 0 as follows

35y1*>y2*3>0x1*=3y1*2y2*-y1*36I0>c¯∞+e-bu(y1*) y_1^* > \frac{y_2^*}{3} > 0 \quad x_1^* = \frac{3y_1^{*^2} y_2^* - y_1^{*^3}}{6I_0} > \bar{c}_\infty + e^{-bu(y_1^*)}

36y2*>y1*3>0x2*=3y2*2y1*-y2*36I0<c¯∞+e-bu(y2*), y_2^* > \frac{y_1^*}{3} > 0 \quad x_2^* = \frac{3y_2^{*^2} y_1^* - y_2^{*^3}}{6I_0} < \bar{c}_\infty + e^{-bu(y_2^*)}

where h shifts the constraints in Eq.(35) and Eq.(36), through the solutions to 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.

As previously mentioned, with the addition of our function h, we had to ensure that the main properties of the system were intact. Note that for all domain values, regardless of the function g, the function h is always positive. More specifically, at the end of Section 2.3, we found that 0.2804 ≤ h ≤ 2.7799 for any function g. We observed that the shape of the solutions to u(c) are maintained at the bounds of h for the full and simplified model, the space (c_s, c), in Figure 1.

Fig. 2

The grey region represents bistability in the (y1*,y2*)(y_1^*, y_2^*) plane for K=13×104K = \frac{1}{3} \times 10^4. The three solid blue lines are defined as y2*=3y1*y_2^* = 3y_1^*, y2*=y1*y_2^* = y_1^*, and y2*=13y1*y_2^* = \frac{1}{3}y_1^*. The red dotted and dashed curves are defined by Eq.(35) and Eq.(36). Parameters used are from Table 2, and the red dot indicates the intersection of Eq.(35) and Eq.(36) when h = 1 for reference.

In Figure 2, the region enclosed by the solid blue and the dotted/dashed red lines is the area of bistability, the gray shaded region, following the conditions outlined in Eqs.(27), (28), (35) and (36). The intersection point of the red curves when h = 1 is defined by the red dot at the point, y1*=1.04y_1^* = 1.04 and y2*=1.87y_2^* = 1.87. Observe in Figure 2 (a) and Figure 2 (b), that the shape and behavior of the graphs is maintained at the bounds of h. We notice Figure 2 (a), the graph corresponding to min(h), has a larger region of bistability compared to the graph with max(h). Recall that the maximum value of h is achieved when an individual consumes high volumes of alcohol, while the minimum value is achieved in the absence of alcohol. In Figure 2, increasing an individualߣs alcohol consumption will decrease the bistability of the HPA axis system, making it less predictable and volatile. The reduced stability results in hypercortisolism in the system/individual due to disruption in the negative feedback loop. The case of min(h) (minimal to no significant volumes of alcohol consumption) will be the most stable value of the system. In the general case, as h → 2.7799, the area of bistability decreases, whereas when h → 0.2804, the area of bistability increases.

Fig. 3

Dotted red - Fast c nullcline, full model; Solid red - Fast y nullcline, simplified model; Dashed blue - Slow c_s nullcline, full model; Solid blue - Slow x nullcline, simplified model.

In Figure 3, the shape of the c_s, c-nullclines is maintained at the bounds of the function h. Figure 3 (a) and Figure 3 (b). Comparing the nullclines at the bounds h = 0.2804 and h = 2.7799, we see that as h → 2.7799 the slow nullclines shift to the left. This shift to the left alludes to the fact that sCRH evolves on a slower timescale, meaning less available sCRH within the main timescale of the model. This is the case when an individual consumes large amounts of alcohol, h = 2.7799, disrupting future stress responses by delaying sCRH replenishment in the hypothalamus. The fast nullclines shift slightly rightward when h → 2.7799, alludes to rCRH evolving faster within the model timescale resulting in more available rCRH. At this limit of h, the increase in rCRH concentration is directly related to hypercortisolism, as rCRH acts as a catalyst for cortisol production. The higher concentration of rCRH also creates a higher basal level of sCRH because of its role as a regulatory signal. Further, if sCRH is evolving slower while rCRH is evolving faster, we see that this phenomenon exhibits more evidence of decreased bistability within the HPA axis in the presence of alcohol consumption. An individual who consumes alcohol frequently and in high quantities will experience continued stress responses due to their BAC consistently being nonzero. This leads to the need for higher baseline hormone levels in an individual as stress tolerance builds up, similar to alcohol tolerance in chronic drinkers. However, in Figure 3, there is not enough sCRH in the system to maintain the continued stress responses needed when alcohol is consumed, causing a delayed decay of cortisol and glucocorticoid complexes due to the inhibition of negative feedback provided by the adrenals. The delay in negative feedback occurs because the stressor, alcohol (as indicated by the individual's BAC), remains present, prompting the adrenals to produce cortisol and therefore pulling more signals from the hypothalamus and pituitary gland to maintain appropriate cortisol levels through the catalysts CRH and ACTH, respectively. However, the concentration of sCRH, the initiating signal, is insufficient, while circulating CRH (rCRH) remains elevated - leading to sustained cortisol production despite inadequate CRH reserves. This imbalance will cause an individual to experience increased stress for a longer period of time until the negative feedback response manifests. The negative feedback process will be activated when sufficient sCRH is present in an individual and therefore when B(t) → 0.

Fig. 4

h2(t,33,μ)h_2(t, 33, \mu): These plots represent the variation in the rate h₂, using g₂, when fixing an age to x = 33 and applying different drinking schedules and alcohol types/concentrations (See Appendix 7). Recall t = 0 corresponds to 8:00 AM. (g₂ and x = 33 were chosen for simulation purposes. The plots utilizing g₁ can be found in the GitHub within Appendix 7).

In Figure 4, we observe the differences of the function h₂ depending on the frequency/quantity of alcohol consumption, the type of alcohol, and gender. These graphs take into account a period of 30 days, i.e. A = 30 in Eq.(8), with a fixed alcohol type, same number of drinks for both males and females for the each day according to the drinking schedule, and an individual does not deviate from the chosen schedule. For all drink types (beer, wine, liquor), the frequency of alcohol consumption, i.e. different drinking schedules, has a positive correlation with h₂. This is expected as frequent alcohol consumption leads to elevated basal cortisol levels and increased basal concentrations of catalytic hormones. Due to this fact, observe that graphs are higher on the h-axis in the everyday and weekday schedules, as an individual will drink more frequently when following those schedules. In general, the number of days drinking during the period of A days for each schedule can be calculated via ⌊A7⌋⋅λ\lfloor \frac{A}{7} \rfloor \cdot \lambda, where λ follows the definition from Table 2. For example, in a period of 30 days, someone following the college schedule (λ = 4) will drink ⌊307⌋⋅4=16\lfloor \frac{30}{7} \rfloor \cdot 4 = 16 times whereas someone following the weekend schedule (λ = 2) will drink ⌊307⌋⋅2=8\lfloor \frac{30}{7} \rfloor \cdot 2 = 8 times. However, it is also important to note that as h₂ → ∞, then c_s, c, a → 0 faster, causing an individual to continue to drink more as a way to “reduce” stress due to social norms [9]. Thus, as the subsystem {c_s, c, a} decays to zero, the only way to reduce cortisol levels, u, is by time. The time it takes is relative to the decay factor of B(t), β. In Figure 4, all the graphs demonstrate that the nadir value is higher than in the absence of alcohol (h₂ = 0.2804), and the peak values are also elevated. Due to the function h accounting for the circadian rhythm, the increase in nadir levels of cortisol causes a disruption in the circadian rhythm, which is in part why many individuals report cases of insomnia after a day/night of drinking. Furthermore, we see that in every graph, females will have higher values of h₂, i.e. cortisol, due to the lower average weight of females than males, differences in biological structure, and biochemical differences (i.e. hormones) that lead to females absorbing more alcohol with a longer period of metabolization. Studying these plots, in correspondence to the results found in the bistability (Figure 2) and nullcline (Figure 3) figures, we can readily observe that there will be a higher concentration of cortisol when consuming alcohol. This is an indicator of Cushing's Syndrome and hypercortisolism, which can lead to obesity, osteoporosis, and various other ailments [42,43].

4.1

Instantaneous circadian impact

4.1.1

Preliminary enhanced circadian function

Thus far, we have examined the cumulative effects of alcohol consumption on both diurnal and nocturnal cortisol levels. This study specifically investigates the impact of alcohol on basal and peak cortisol concentrations in individuals who consume alcohol at varying quantities and frequencies. In addition to these cumulative effects, alcohol consumption has been shown to induce a phase shift (phase delay) in the circadian rhythm, which typically exhibits a sinusoidal pattern [44]. This phase shift is analogous to that observed in cases of jet lag, where the circadian clock becomes misaligned with the astronomical clock [8]. For the purpose of this study, we define the standard circadian rhythm as the uninterrupted rhythm in the absence of alcohol consumption, and the affected circadian rhythm as the disrupted rhythm resulting from alcohol intake. Accordingly, any quantity or type of alcohol consumed results in a phase shift from the standard to the affected circadian rhythm. It is important to emphasize that an individual possesses only one circadian rhythm; the use of two rhythm types in this description is intended solely for illustrative purposes. In Section 2.1, we defined the circadian rhythm function, s, as:

s(t)=cos(t)+23, s(t) = \frac{\cos(t)+2}{3}

where s: [0, 2π] → [1/3, 1]. Clearly, the function is not impacted by an individualߣs BAC, however, we can introduce the Widmark equation which is given by Eq.(4) to formulate a phase shift φ(t˜)\varphi(\tilde{t}). We define, φ(t˜)=B(t˜)=δzNσW−βt˜\varphi(\tilde{t}) = B(\tilde{t}) = \frac{\delta z N}{\sigma W} - \beta \tilde{t}. The phase shift, φ(t˜)\varphi(\tilde{t}), is dependent upon an individuals BAC at a time φ(t˜)\varphi(\tilde{t}). Observe that φ(t˜)\varphi(\tilde{t}) is measured in kg/L; therefore, we require a dimensionless equivalent for compatibility within the cosine function, as t is in hours. To achieve this, φ(t˜)\varphi(\tilde{t}) is normalized to make it dimensionless. Naturally, the norm is chosen to be the maximum norm, which corresponds to the legal limit of BAC in the United States (0.0008 kg/L) [33]. Therefore, it follows that:

37s(t,t~)=s(t−φ~(t~))=s(t−1BACmaxφ~(t~))=cos(t−10.0008(δzNδW−βt2))+23, \overline{\varphi}(\tilde{t}) = \frac{1}{BAC_{\max}}\varphi(\tilde{t}) = \frac{\cos\left(1-\frac{1}{0.0008}\left(\frac{\delta z N}{\sigma W} - \beta \tilde{t}\right)\right)+2}{3}

where φ¯(t˜)=1BACmaxφ(t˜)\overline{\varphi}(\tilde{t}) = \frac{1}{BAC_{max}}\varphi(\tilde{t}). Note that φ¯(t˜)\overline{\varphi}(\tilde{t}), the dimensionless phase shift, is calculated with t˜\tilde{t} in hours, while the phase time t is computed in the radian timescale. Both t˜\tilde{t} and t represent the same time but in different forms through the relation π ≡ 12 hours, as described in Section 2.1. We further impose that if φ(t˜)<0\varphi(\tilde{t}) < 0 then φ(t˜)=0\varphi(\tilde{t}) = 0, since a negative blood alcohol content is not biologically possible. Hence, with that constraint, s(t−φ¯(t˜))=s(t)s(t - \overline{\varphi}(\tilde{t})) = s(t) when N = 0, i.e. when an individual has a BAC of 0kg/L.

Fig. 5

Plots for consuming ten drinks following the everyday schedule for beer, wine, and liquor. The time of drinking is 9:00 PM, 13 hours after peak cortisol levels (8:00 AM). The top subplot for each plot shows the linear decay of the BAC level for each occurrence of drinking. In the bottom plots, the solid purple line represents the individual's circadian rhythm while the black dotted line represents the circadian rhythm if no alcohol was consumed. The vertical red dotted lines (t˜=24,48,…,168\tilde{t} = 24, 48, \ldots, 168 hours) represent a 24 hour block of time (8:00AM - 8:00AM (next day)).

In Figure 5, a seven-day period is modeled in which a male (left column) and a female (right column) each consume 10 drinks of different types of alcohol (one type per row), following the everyday drinking schedule described in Section 3.2. Recall that, in our model, the time required to consume N drinks is not considered; thus, the effects on an individualߣs BAC are assumed to be instantaneous. This assumption is evident in the upper plots, where the leftmost edge of each triangular BAC curve is vertical. In some of these plots, we observe peak BAC levels reaching approximately 0.004kg/L. This unrealistically high value results from the assumption that all N drinks are consumed instantaneously. In practice, an individualߣs BAC would not reach such a potentially fatal level [45] because alcohol is continuously metabolized at a rate β while the drinks are consumed over time. A notable contrast between the male and female plots is that the female's circadian rhythm is disrupted more significantly by alcohol consumption. This is primarily due to physiological differences; on average, males weigh more than females [46], and body weight directly affects BAC through the Widmark equation given by Eq.(4). A higher body weight implies a larger blood volume, resulting in a lower BAC for the same amount of alcohol consumed. That is, the greater the weight, the more diluted the alcohol concentration (kg/L) becomes. Furthermore, the type of alcohol consumed, determined by the alcohol content per fluid ounce, plays a crucial role in the magnitude of circadian rhythm disruption, with liquor causing the most significant phase shifts. Recall that cortisol is a key hormone involved in the circadian rhythm and plays a central role in the sleep-wake cycle [47]. This relationship is illustrated in the baseline case, where φ(t˜)=0\varphi(\tilde{t}) = 0 (indicated by the black dotted line). During the seven-day period, the maximum of s(t), our modeled circadian rhythm curve, consistently aligns with the vertical red dotted lines at 8:00 AM, the typical time of peak cortisol levels, facilitating the transition from sleep to wakefulness. The gradual increase in s(t) approaching 8:00 AM reflects the natural rise in cortisol that helps to initiate alertness upon waking. However, when φ(t˜)>0\varphi(\tilde{t}) > 0, this alignment is disrupted. The peak of s(t) is delayed and occurs several hours after the expected time. The magnitude of this delay depends on the value of φ(t˜)\varphi(\tilde{t}), which in turn is influenced by the number of drinks consumed (N), the type of alcohol, the individual's gender, and body weight. This explains why females generally experience more pronounced circadian disruptions than males, and why certain alcohol types (e.g. liquor) induce larger phase shifts. Table 3 summarizes the number of hours by which the cortisol peak is shifted in each plot presented in Figure 5.

Table 3

	Beer	Wine	Liquor
Male	max(φ¯(t˜))\overline{\varphi}(\tilde{t}) = 3.182 hours	max(φ¯(t˜))\overline{\varphi}(\tilde{t}) = 3.910 hours	max(φ¯(t˜))\overline{\varphi}(\tilde{t}) = 4.171 hours
Female	max(φ¯(t˜))\overline{\varphi}(\tilde{t}) = 3.659 hours	max(φ¯(t˜))\overline{\varphi}(\tilde{t}) = 4.574 hours	max(φ¯(t˜))\overline{\varphi}(\tilde{t}) = 4.879 hours

In Table 3, the values of φ¯(t˜)\overline{\varphi}(\tilde{t}) are inherently dimensionless. However, in the context of inducing a shift in the circadian rhythm, it is evident that these values can be interpreted as having units of hours, as illustrated in Figure 5 and Figure 6. At the time alcohol is consumed, we observe an immediate stress response, characterized by an increase in cortisol, which induces a shift in the timing of peak cortisol levels. Practically, this shift would occur more gradually (smoother) and would be proportional to the rate of alcohol consumption. The rightward shift further contributes to post-drinking symptoms such as veisalgia, commonly referred to as hangover in everyday vernacular, which is caused by elevated cortisol levels during the day following alcohol consumption [48]. The plots in Figure 5 represent individuals with chronic drinking habits who consume ten drinks everyday. Clearly, the severity of the disruption to an individualߣs circadian rhythm is pronounced in those figures. Below, in Figure 6, we again followed the everyday schedule; however, instead of fixing N = 10 daily, we generated a random number (between 1 and 10) of drinks for each day. This represents potentially a more practical representation of drinking habits/volumes for an individual consuming alcohol daily. The same observations drawn from Figure 5 can be made for Figure 6, however, there will be minimal rightward shifts if N is sufficiently small for a given day; nevertheless, a disruption in circadian rhythm will still be present.

Fig. 6

Plots for consuming a random number of drinks each day (Nj=randint{1,10})(N_j = \text{randint}\{1,10\}), following the everyday schedule for beer, wine, and liquor and varying W as needed. The value of N_j for each day can be found in Appendix 7.

Table 4

	Beer	Wine	Liquor
Male	0.3128 ≤ φ¯(t˜)\overline{\varphi}(\tilde{t}) ≤ 2.5025 (hours)	0.3910 ≤ φ¯(t˜)\overline{\varphi}(\tilde{t}) ≤ 3.5191 (hours)	0.4171 ≤ φ¯(t˜)\overline{\varphi}(\tilde{t}) ≤ 3.7538 (hours)
Female	0.3659 ≤ φ¯(t˜)\overline{\varphi}(\tilde{t}) ≤ 2.9274 (hours)	0.9148 ≤ φ¯(t˜)\overline{\varphi}(\tilde{t}) ≤ 4.1167 (hours)	0.4879 ≤ φ¯(t˜)\overline{\varphi}(\tilde{t}) ≤ 3.4153 (hours)

In Table 4, we present the closed interval in which φ¯(t˜)\overline{\varphi}(\tilde{t}) lies, based on the randomly determined number of drinks N_j. For all types of drinks, the lower bounds for males and females will result in minimal horizontal shifts in the circadian rhythm, potentially leading only to a sharp increase in cortisol, as observed on certain days in Figure 6. On the other hand, we observe significant horizontal shifts in the bottom plots of Figure 6, s, when φ¯(t˜)→max(φ¯(t˜))\overline{\varphi}(\tilde{t}) \to \max(\overline{\varphi}(\tilde{t})) for each range.

4.1.2

Enhanced circadian function

To further refine our circadian shift analysis, we constructed a piecewise function, s, based on the principles established in Section 4.1.1. This function is designed to reach its minimum at t=3π2t = \frac{3\pi}{2} corresponding to 2:00 AM (see Section 2.1), thereby aligning more accurately with the circadian rhythm than prior iterations of the function s. We define the piecewise function, s, as follows:

38s(t,t~)=s(t−φ~(t~))={cos(23θ−φ~(t~))+23,if θ∈[0,3π2),cos(2θ−φ~(t~))+23,if θ∈[3π2,2π), s(t,\tilde{t}) = s(t-\overline{\varphi}(\tilde{t})) = \begin{cases} \frac{\cos(\theta - \overline{\varphi}(\tilde{t}))+2}{3}, & \text{if } \theta \in [0,\pi], \\ \frac{\cos(2\theta - \overline{\varphi}(\tilde{t}))+2}{3}, & \text{if } \theta \in (\pi, 2\pi), \end{cases}

where θ := t (mod 2π). Reducing t modulo 2π ensures continuity by mapping all instances of t = 2πk (for any integer k) back to θ = 0. This is significant because ≡ does not define an equivalence relation in the standard sense ([0] ≠ [2πk]) for all positive integers k. As a result, the second case in the piecewise function excludes 2π, and values of θ that are integer multiples of 2π are instead handled by the first case, where θ = 0. The motivation for this piecewise function analysis is to model the antisymmetric nature of circadian rhythm about 12 hours from wake (8:00 AM) [15]. Recall that t˜\tilde{t}, in hours, is not reduced modulo 2π since uses a non-radian timescale, unlike t in radians, with both related through π ≡ 12 hours, as noted in Section 4.1.1. For example, if t=22π7t = \frac{22\pi}{7}, then θ≡22π7≡8π7 (mod 2π)\theta \equiv \frac{22\pi}{7} \equiv \frac{8\pi}{7} \text{ (mod } 2\pi) while t˜=22π7=22(12)7=37.714\tilde{t} = \frac{22\pi}{7} = \frac{22(12)}{7} = 37.714 hours. With this piecewise function, the drinking time of 9:00 PM is now associated with a decrease in cortisol (circadian rhythm approaching the nadir) as supposed to during an increase in cortisol concentration (circadian rhythm approaching the peak) in Section 4.1.1. Additionally, this section aims to substantiate findings from previous iterations of s, demonstrating that the prior symmetric/non-piecewise form retains accurate results.

Fig. 7

Plots for consuming ten drinks following the everyday schedule for beer, wine, and liquor. The time of drinking is 9:00 PM, 13 hours after peak cortisol levels (8:00 AM). The top subplot for each plot shows the linear decay of the BAC level for each occurrence of drinking. In the bottom plots, the solid purple line represents the individual's circadian rhythm while the black dotted line represents the circadian rhythm if no alcohol was consumed. The vertical red dotted lines (t = 24, 48, ..., 168 hours) represent a 24 hour block of time (8:00AM - 8:00AM (next day)).

Table 5

	Beer	Wine	Liquor
Male	max(φ¯(t˜))\overline{\varphi}(\tilde{t}) = 3.182 hours	max(φ¯(t˜))\overline{\varphi}(\tilde{t}) = 3.910 hours	max(φ¯(t˜))\overline{\varphi}(\tilde{t}) = 4.171 hours
Female	max(φ¯(t˜))\overline{\varphi}(\tilde{t}) = 3.659 hours	max(φ¯(t˜))\overline{\varphi}(\tilde{t}) = 4.574 hours	max(φ¯(t˜))\overline{\varphi}(\tilde{t}) = 4.879 hours

In Figure 7, we observe that the horizontal shift now occurs while s is decreasing, compared to Figure 5, where the horizontal shift occurs during an increase in s. Despite this notable difference between the two enhanced circadian rhythm functions, Eq. (37) and Eq.(38), we still observe similar phase shifts when alcohol is consumed at 9:00 PM. Furthermore, the values of φ¯(t˜)\overline{\varphi}(\tilde{t}) are equivalent in Table 3 and Table 5, as φ¯(t˜)\overline{\varphi}(\tilde{t}) is solely dependent on an individualߣs BAC when consuming ten drinks ((Nj=10)N_j = 10) of differing alcoholic content. This supports the results obtained in Section 4.1.1; however, the piecewise function introduced in this section is more biologically accurate. In Figure 8, we observe the induced phase shifts to the circadian rhythm when imposed with a random N_j, as done in Section 4.1.1. The conclusions drawn from Figure 8 are consistent with those of Figure 6. Again, these results reinforce the previous findings in Section 4.1.1, with the primary distinction between the resulting figures being the timing of the nadir cortisol concentration and the point in the circadian rhythm at which the phase shift occurs.

Fig. 8

Plots for consuming a random number of drinks each day, following the everyday schedule for beer, wine, and liquor. The value of N_j for each day can be found in Appendix 7.

Table 6

	Beer	Wine	Liquor
Male	0.3128 ≤ φ¯(t˜)\overline{\varphi}(\tilde{t}) ≤ 2.8153 (hours)	0.3910 ≤ φ¯(t˜)\overline{\varphi}(\tilde{t}) ≤ 3.5191 (hours)	0.8342 ≤ φ¯(t˜)\overline{\varphi}(\tilde{t}) ≤ 2.9196 (hours)
Female	0.3659 ≤ φ¯(t˜)\overline{\varphi}(\tilde{t}) ≤ 3.2933 (hours)	1.3722 ≤ φ¯(t˜)\overline{\varphi}(\tilde{t}) ≤ 4.1167 (hours)	0.9758 ≤ φ¯(t˜)\overline{\varphi}(\tilde{t}) ≤ 2.9274 (hours)

Our results indicate hypercortisolism associated with increased alcohol consumption. The severity of the dysregulation within the HPA axis depends on the amount of alcohol consumed (floz. of a drink and its alcohol concentration), gender, weight, alcohol elimination rate, and frequency of alcohol consumption which our model considers. Furthermore, we observed severe impacts to the circadian rhythm. These impacts were observed in vertical and lateral shifts from the standard circadian rhythm when alcohol was consumed. Again, the severity of the vertical and lateral shifts is due to factors that determine an individualߣs BAC (Section 2.3). Our model highlighted the basic differences between the activation of the HPA axis between males and females under alcohol consumption, as well as considered the basic age-related differences. Alcohol consumption can significantly disrupt physiological processes, potentially leading to hypercortisolism and circadian rhythm disturbances, as demonstrated by our model.