# Why slow breathing near 6 breaths/min enlarges heart-rate oscillations: a delayed negative-feedback
# loop (simplified model of the baroreflex, after the resonance account of Vaschillo et al. 2006).
# Breathing pushes heart rate (respiratory sinus arrhythmia); blood pressure follows heart rate with a
# vascular lag; the baroreflex corrects heart rate after a delay. When the loop's total lag equals half a
# breathing cycle, the correction arrives in step with the next push and the swings add up.
import cmath, math

G, DELAY, TAU_V = 0.8, 3.5, 2.0       # loop gain, reflex delay (s), vascular lag (s): illustrative values

def gain(f):
    w = 2 * math.pi * f
    loop = G * cmath.exp(-1j * w * DELAY) / (1 + 1j * w * TAU_V)
    return abs(1 / (1 + loop))

def simulate(f, seconds=300, dt=0.01):
    """Euler simulation of the same loop; returns peak-to-peak heart-rate swing in the last minute."""
    n_delay = int(DELAY / dt)
    hist = [0.0] * n_delay
    p, lo, hi = 0.0, 1e9, -1e9
    for k in range(int(seconds / dt)):
        t = k * dt
        h = math.sin(2 * math.pi * f * t) - G * hist[k % n_delay]   # pressure from DELAY seconds ago
        hist[k % n_delay] = p
        p += dt * (h - p) / TAU_V
        if t > seconds - 60:
            lo, hi = min(lo, h), max(hi, h)
    return hi - lo

print("breaths/min   frequency (Hz)   heart-rate swing (x of the breathing push alone)")
for bpm in (4, 5, 6, 7, 8, 10, 12, 15):
    f = bpm / 60
    print(f"{bpm:8d}        {f:6.3f}           {gain(f):4.2f}   (simulation {simulate(f) / 2:4.2f})")
best = max((b / 10 for b in range(30, 200)), key=lambda b: gain(b / 60))
print(f"peak of the loop: {best:.1f} breaths/min ({best / 60:.3f} Hz)")
# breaths/min   frequency (Hz)   heart-rate swing (x of the breathing push alone)
#        4         0.067           1.20   (simulation 1.20)
#        5         0.083           1.73   (simulation 1.73)
#        6         0.100           1.99   (simulation 1.99)
#        7         0.117           1.64   (simulation 1.64)
#        8         0.133           1.28   (simulation 1.28)
#       10         0.167           0.92   (simulation 0.92)
#       12         0.200           0.81   (simulation 0.80)
#       15         0.250           0.82   (simulation 0.82)
# peak of the loop: 5.8 breaths/min (0.097 Hz)

# RMSSD and SDNN from a simulated series of heartbeat (RR) intervals while breathing at 6 or 12 per minute.
# The respiratory swing is taken from the loop model above (about 2.0 x at 6/min, 0.8 x at 12/min).
import math, random, statistics

def rr_series(bpm, swing, minutes=5, mean_rr=850.0, base_amp=25.0, noise=12.0, seed=2):
    random.seed(seed)
    t, rr = 0.0, []
    while t < minutes * 60:
        x = mean_rr + base_amp * swing * math.sin(2 * math.pi * bpm / 60 * t) + random.gauss(0, noise)
        rr.append(x)
        t += x / 1000
    return rr

def rmssd(rr):
    d = [b - a for a, b in zip(rr, rr[1:])]
    return math.sqrt(sum(x * x for x in d) / len(d))

for bpm, swing in ((12, 0.81), (6, 1.99)):
    rr = rr_series(bpm, swing)
    print(f"breathing {bpm:2d}/min: {len(rr)} beats, SDNN {statistics.pstdev(rr):5.1f} ms, RMSSD {rmssd(rr):5.1f} ms")
print("Slow breathing nearly doubles the slow swing (SDNN) but changes RMSSD much less, because")
print("RMSSD looks at beat-to-beat differences, and a slow wave changes little from one beat to the next.")
# breathing 12/min: 354 beats, SDNN  19.0 ms, RMSSD  22.9 ms
# breathing  6/min: 354 beats, SDNN  36.2 ms, RMSSD  25.0 ms
# Slow breathing nearly doubles the slow swing (SDNN) but changes RMSSD much less, because
# RMSSD looks at beat-to-beat differences, and a slow wave changes little from one beat to the next.