# A minimal CRH -> ACTH -> cortisol loop with glucocorticoid-receptor (GR) feedback, arbitrary units
# (structure after Vinther, Andersen & Ottesen 2011 and Gupta et al. 2007; no daily or ultradian rhythm).
# Feedback strength G = GR availability (1 = intact). Cortisol above 1.4 slowly lowers G
# (down-regulation); G recovers with a time constant of a week.
import random
H1, H2, H3, DT = 0.10, 0.07, 0.012, 0.5      # clearance per minute; cortisol half-life about 1 h

def run(starts, minutes, state=(1.0, 1.0, 1.0, 1.0)):
    """starts: minutes at which a 30-minute stressor begins. Returns final state and cortisol trace."""
    R, A, C, G = state
    on = set()
    for s in starts:
        on.update(range(int(s / DT), int((s + 30) / DT)))
    trace = []
    for k in range(int(minutes / DT)):
        drive = 5.0 if k in on else 1.0
        fb = 1 - 0.9 * G * C ** 2 / (1 + C ** 2)          # Hill-type GR-mediated feedback
        R += DT * H1 * (2 * drive * fb - R)
        A += DT * H2 * (2 * R * fb - A)
        C += DT * H3 * (A - C)
        G += DT * ((1 - G) / (7 * 1440) - 0.00025 * max(C - 1.4, 0) * G)
        trace.append(C)
    return (R, A, C, G), trace

random.seed(4)
print("30 days, stressors per day         GR left  resting cortisol  test peak  recovery (min)")
for label, n, rnd in (("0 (rest)", 0, False), ("1, same time, predictable", 1, False),
                      ("3, random times", 3, True), ("6, random times", 6, True), ("9, random times", 9, True)):
    starts = [d * 1440 + (random.uniform(0, 1410) if rnd else 600) for d in range(30) for _ in range(n)]
    state, _ = run(starts, 30 * 1440)
    state, _ = run([], 3 * 1440, state)                   # three quiet days: the new resting state
    rest = state[2]
    _, tr = run([0], 1440, state)                          # one standard 30-minute stressor
    back = next((i * DT for i, c in enumerate(tr) if i * DT > 30 and c < 1.1 * rest), None)
    rec = f"{back:8.0f}" if back else "  > 1440"
    print(f"{label:32s} {state[3]:6.2f}  {rest:15.2f}  {max(tr):9.2f}  {rec}")
# 30 days, stressors per day         GR left  resting cortisol  test peak  recovery (min)
# 0 (rest)                           1.00             1.07       1.71       121
# 1, same time, predictable          0.98             1.09       1.73       122
# 3, random times                    0.93             1.13       1.83       128
# 6, random times                    0.19             2.88       5.51       262
# 9, random times                    0.18             2.92       5.58       262

# How "allostatic load" is usually scored in studies (after Seeman et al. 1997): for each biomarker, count
# 1 if a person falls in the high-risk quartile of the sample; the score is the count. It is an index for
# research on groups, not a diagnosis and not a direct measurement of "load". Simulated sample.
import random
random.seed(7)
MARKERS = {   # name: (mean, sd, risk direction)
    "systolic BP (mmHg)": (125, 15, +1), "diastolic BP (mmHg)": (78, 10, +1), "waist-hip ratio": (0.90, 0.07, +1),
    "HDL cholesterol (mmol/L)": (1.35, 0.35, -1), "total/HDL ratio": (4.2, 1.1, +1), "HbA1c (%)": (5.5, 0.5, +1),
    "overnight urinary cortisol": (1.0, 0.35, +1), "noradrenaline": (1.0, 0.35, +1), "adrenaline": (1.0, 0.35, +1),
    "DHEA-S": (1.0, 0.35, -1),
}
people = [{m: random.gauss(mu, sd) for m, (mu, sd, _) in MARKERS.items()} for _ in range(1000)]
cut = {}
for m, (_, _, sign) in MARKERS.items():
    v = sorted(p[m] for p in people)
    cut[m] = v[749] if sign > 0 else v[250]           # top quartile, or bottom quartile when low is risky
def score(p):
    return sum(1 for m, (_, _, s) in MARKERS.items() if (p[m] >= cut[m] if s > 0 else p[m] <= cut[m]))
scores = [score(p) for p in people]
for k in range(0, 8):
    n = scores.count(k)
    print(f"score {k}: {n:4d} people {'#' * (n // 10)}")
print(f"mean score {sum(scores) / len(scores):.2f} of 10 (with independent markers the expected mean is 2.5)")
print("The same person can score higher simply because the comparison sample is healthier: the cut-offs are relative.")
# score 0:   46 people ####
# score 1:  186 people ##################
# score 2:  292 people #############################
# score 3:  265 people ##########################
# score 4:  134 people #############
# score 5:   60 people ######
# score 6:   13 people #
# score 7:    1 people
# mean score 2.51 of 10 (with independent markers the expected mean is 2.5)
# The same person can score higher simply because the comparison sample is healthier: the cut-offs are relative.