The HeartRateLab.jl: a toolkit for heart beats time series

Alberto Barradas

Overview

  • Existing efforts
  • Introduction to HRV
  • Inter-Beat-Intervals
  • HRV Domains and Metrics
  • The HeartRateLab
    • Preprocessing
    • Feature Domains
    • Modeling
  • Applications
  • Metascience and Agentic tools

HRV Analysis

  • Electro Cardiography (ECG)
  • photoplethysmography (PPG)
  • Series of inter-beat-intervals (IBIs)

Cardiac Cycle and Heart Rate

By DrJanaOfficial - Own work: Cardiac Cycle, Wiggers diagram derived from User:Xavax, CC BY-SA 4.0, https://commons.wikimedia.org/w/index.php?curid=54393545

Cardiac Cycle

Sinoatrial and Atrioventricular Nodes

By J. Heuser - self made, based upon Image:Heart anterior view coronal section.jpg by Patrick J. Lynch (Patrick J. Lynch; illustrator; C. Carl Jaffe; MD; cardiologist Yale University Center for Advanced Instructional Media ), CC BY 2.5, https://commons.wikimedia.org/w/index.php?curid=1686121

Heart Rate Variability (HRV)

By YitzhakNat - Own work + ECG wave from, CC BY-SA 4.0, https://commons.wikimedia.org/w/index.php?curid=121207798

Interbeat Intervals

Code 
# infile = joinpath(@__DIR__, "..", "..", "test", "testdata", "example.txt")
# infile = joinpath(@__DIR__, "..", "..", "test", "testdata", "e1304.txt")
infile = joinpath(@__DIR__, "example.xdf")
data = Float64.(HeartRateLab.read_xdf(infile))

println("Data loaded from: ", infile)
println("Data length: $(length(data)) heartbeats")
println("Type of data: ", typeof(data))
Data loaded from: /home/beto/.julia/dev/HeartRateLab/.worktrees/cl/docs/slides/example.xdf
Data length: 341 heartbeats
Type of data: Vector{Float64}

Generic Heart Rate Monitor

Polar H10

Optical (PPG) sensors

Polar Verity Sense - optical armband

COROS Heart Rate Monitor - optical armband

Data Output

Although some heart rate bands have the capability of recording ECG, communicating the IBIs is more effective for low energy devices.

Code 
f = Figure(size=(1200, 600));
ax = Axis(f[1, 1], 
    title="Raw Data",
    xlabel="Time (s)",
    ylabel="IBI (ms)")
n = data[1:50]
t = Observable(cumsum(n) ./ 1000) # to seconds
# x = Observable(60000 ./ n) # to bpm
x = Observable(n)
stem!(ax, t, x, color=:blue)
ylims!(ax, (minimum(n) - 20, maximum(n) + 20))
f

IBIS are inversely proportional to Heart Rate

Code 
f = Figure(size=(1400, 600));
ax1 = Axis(f[1, 1], 
    title="Heart Rate",
    xlabel="Time (s)",
    ylabel="Heart Rate (bpm)")
n = data[1:50]
t = Observable(cumsum(n) ./ 1000) # to seconds
x = Observable(60000 ./ n) # to bpm
lines!(ax1, t, x, color=:blue, label="Heart Rate")
ax2 = Axis(f[1, 2], 
    title="Inter-Beat Intervals",
    xlabel="Time (s)",
    ylabel="IBI (ms)")
lines!(ax2, t, n, color=:red, label="IBI")
f
Code 
f = Figure(size=(1600, 600));
ax = Axis(f[1, 1], 
    title="Heart Rate from Inter-Beat Intervals",
    xlabel="Heart Rate (bpm)",
    ylabel="Time (s)")
n = data[1:50]
t = Observable(cumsum(n) ./ 1000) # to seconds
x = Observable(60000 ./ n) # to bpm
lines!(ax, t, x, color=:teal, label="Continuous Heart Rate")
hlines!(ax, mean(x[]), label="Mean Heart Rate", linestyle=:dash, color=:coral, linewidth=2)

f
Code 
f = Figure(size=(1600, 600));
ax = Axis(f[1, 1], 
    title="Heart Rate from Inter-Beat Intervals",
    xlabel="Heart Rate (bpm)",
    ylabel="Time (s)")
n = data
t = Observable(cumsum(n) ./ 1000) # to seconds
x = Observable(60000 ./ n) # to bpm
lines!(ax, t, x, color=:teal, label="Continuous Heart Rate")
hlines!(ax, mean(x[]), label="Mean Heart Rate", linestyle=:dash, color=:coral, linewidth=2)

f

HRV Domains and Features

  • Time & statistics
  • Frequency
  • Geometric
  • Nonlinear

Time Domain

  • Mean
  • Standard Deviation (SDNN)
  • Median
  • Extrema and Range
  • BPM representations
  • SDSD
  • RMSSD
  • SDANN
  • PNN50, PNN20 and percentages
  • CVSD
  • rRR

An empirical random variable

Code 
# A histogram of the first 50 heartbeats
f = Figure(size=(1400, 600));
ax1 = Axis(f[1, 1], 
    title="Histogram of measured Heart Rate",
    xlabel="Heart Rate (bpm)",
    ylabel="Frequency")
n = data[1:50]
t = Observable(cumsum(n) ./ 1000) # to seconds
x = Observable(60000 ./ n) # to bpm
hist!(ax1, x, color=:teal, bins=40)
vlines!(ax1, mean(x[]), color=:red, linewidth=2, label="Mean")

ax2 = Axis(f[1, 2], 
    title="Estimated distribution of measured Heart Rate",
    xlabel="Heart Rate (bpm)",
    ylabel="Frequency")
# Makie.boxplot!(ax2, x, color=:teal, label="Heart Rate Boxplot"; orientation=:horizontal)
density!(ax2, x, color=:coral, label="Distribution Density")
vlines!(ax2, mean(x[]), color=:red, linewidth=2, label="Mean")

f

Numeric Derivative

Code 
f = Figure(size=(1200, 600));
ax = Axis(f[1, 1], 
    title="Difference of Sequential Inter-Beat Intervals",
    xlabel="Time (s)",
    ylabel="Inter-Beat Interval (ms)")
n = data[1:end]
t = Observable(cumsum(n) ./ 1000) # to seconds
x = Observable([0; diff(n)])
x_above_zero = Observable(clamp.(x[], 0, Inf))
x_below_zero = Observable(clamp.(x[], -Inf, 0))
zero_line = Observable(zeros(length(x[])))
lines!(ax, t, x, color=:teal)
# Band above zero
# band(x, ylower, yupper; kwargs...)
band!(ax, t, zero_line, x_above_zero,
    color=:teal, alpha=0.5)
band!(ax, t, zero_line, x_below_zero,
    color=:coral, alpha=0.5)
# The +/-50 ms threshold: successive differences beyond it are the ones counted
# by pNN50. Paint those beats blue, the rest red.
hlines!(ax, [50.0, -50.0], color=:gray, linestyle=:dash, linewidth=1)
pnn50_color = [abs(v) > 50 ? :blue : :red for v in x[]]
scatter!(ax, t, x, color=pnn50_color, markersize=8)
f

RMSSD

Code 
n = data[1:50]
t = Observable(cumsum(n) ./ 1000) # to seconds
x = Observable([0; diff(n)])
x_above_zero = Observable(clamp.(x[], 0, Inf))
x_below_zero = Observable(clamp.(x[], -Inf, 0))
zero_line = Observable(zeros(length(x[])))
f = Figure(size=(1600, 600));
ax1 = Axis(f[1, 1], 
    title="Difference of Sequential IBIs",
    xlabel="Time (s)",
    ylabel="ΔTime (ms)")
lines!(ax1, t, x, color=:teal, label="ΔRR[n, n+1]")
# Band above zero
# band(x, ylower, yupper; kwargs...)
band!(ax1, t, zero_line, x_above_zero, 
    color=:teal, alpha=0.5)
band!(ax1, t, zero_line, x_below_zero, 
    color=:coral, alpha=0.5)
scatter!(ax1, t, x, color=:red, markersize=5)

ax2 = Axis(f[1, 2], 
    title="IBIs and its Numeric Derivative",
    xlabel="Time (s)",
    ylabel="ΔTime (ms)")
lines!(ax2, t, x, color=:teal, alpha=0.5)
band!(ax2, t, zero_line, x_above_zero, color=:teal, alpha=0.25)
band!(ax2, t, zero_line, x_below_zero, color=:coral, alpha=0.25)
scatter!(ax2, t, x, color=:red, markersize=5, alpha=0.5)

ax3 = Axis(f[1, 2], 
    xlabel="Time (s)",
    ylabel="RR-interval (ms)",
    yticklabelcolor=:blue,
    yaxisposition=:right
    )
lines!(ax3, t, n, color=:blue, linewidth=3, label="Detrended RR-interval signal")

f
Code 
println("SDNN: ", round(StatsBase.std(n), digits=2))
println("RMSSD: ", round(StatsBase.std(x[]), digits=2))
SDNN: 85.4
RMSSD: 30.03

Frequency Domain

  • Power Spectral Density (PSD) over four frequency bands:
    • High Frequency (HF): 9-24 BPM (0.15-0.4 Hz)
    • Low Frequency (LF): 2.4-9 BPM (0.04-0.15 Hz)
    • Very Low Frequency (VLF): .18-2.4 BPM (0.003-0.04 Hz)
    • Ultra Low Frequency (ULF): < .18 BPM (< 0.003 Hz)
  • Absolute and relative power
  • Peak frequencies
  • Proportions and ratios (LF/HF, VLF/LF, etc.)

Geometric Domain

  • Poincaré Plot
    • SD1: Short-term variability
    • SD2: Long-term variability
    • SD1/SD2 ratio, area
    • cvi
    • ccsi
  • Histogram Based:
    • Triangular Index
    • Triangular Interpolation of NN intervals (TINN)
Code 
function getellipsepoints(cx, cy, rx, ry, θ)
    t = range(0, 2*pi, length=100)
    ellipse_x_r = @. rx * cos(t)
    ellipse_y_r = @. ry * sin(t)
    R = [cos(θ) sin(θ); -sin(θ) cos(θ)]
    r_ellipse = [ellipse_x_r ellipse_y_r] * R
    x = @. cx + r_ellipse[:,1]
    y = @. cy + r_ellipse[:,2]
    (x,y)
end
x = [data[i] for i in 1:length(data)-1]
y = [data[i+1] for i in 1:length(data)-1]
sd1 = sqrt(StatsBase.var((x - y) / sqrt(2)))
sd2 = sqrt(StatsBase.var((x + y) / sqrt(2)))
f = Figure(size=(1200, 600));
ax = Axis(f[1, 1], 
    title="Poincaré Plot",
    xlabel="ΔRR(n) (ms)",
    ylabel="ΔRR(n+1) (ms)")
scatter!(ax, x, y, markersize=5, label="ΔRR[i, i+1]")
cx = mean(x)
cy = mean(y)
ellipse_x, ellipse_y = getellipsepoints(cx, cy, sd2, sd1, π/4)
lines!(ax, ellipse_x, ellipse_y, linewidth=2, label="Variance Ellipse", color=:coral)
f[1,2] = Legend(f, ax, "Poincaré Plot",
    position=:bottomright, 
    title="Poincaré Plot Legend",
    fontsize=10,
    bgcolor=:white,
    framecolor=:black)

f

Nonlinear Domain

  • Approximate Entropy (ApEn)
  • Sample Entropy (SampEn)
  • Hurst Exponent
  • Renyi Entropy
  • Detrended Fluctuation Analysis (α1, α2)
Code 
# Congruent with HeartRateLab's DFA defaults (Features.config): an INTEGER box
# grid, order-1 detrending, non-overlapping windows - α1 over boxes 4-16 and
# α2 over 16-64.
ns   = collect(4:64)
fluc = [DFA.dfa(data, k; order=1, overlap=0.0) for k in ns]
logn = log10.(Float64.(ns))
logf = log10.(fluc)

i1 = findall(n -> 4 <= n <= 16, ns)     # α1 (short-term) boxes
i2 = findall(n -> 16 <= n <= 64, ns)    # α2 (long-term) boxes
b1, α1 = DFA.polyfit(logn[i1], logf[i1])
b2, α2 = DFA.polyfit(logn[i2], logf[i2])

f = Figure(size=(1400, 650));
ax = Axis(f[1, 1],
    title="DFA on this recording:  α1 = $(round(α1, digits=3))   |   α2 = $(round(α2, digits=3))",
    xlabel="log10(box size n)",
    ylabel="log10 F(n)",
)
scatter!(ax, logn, logf, color=:black, markersize=8)
x1 = range(minimum(logn[i1]), maximum(logn[i1]), length=50)
lines!(ax, x1, b1 .+ α1 .* x1, color=:crimson, linewidth=3,
    label="α1 (4-16) = $(round(α1, digits=3))")
x2 = range(minimum(logn[i2]), maximum(logn[i2]), length=50)
lines!(ax, x2, b2 .+ α2 .* x2, color=:royalblue, linewidth=3,
    label="α2 (16-64) = $(round(α2, digits=3))")
axislegend(ax; position=:lt)
f

The HeartRateLab

Preprocessing

  • Input from TXT: EliteHRV, Kubios
  • Input as WFDB: PhysioNet
  • Input as XDF: sub ms precision, markers
  • Live input from LSL:
    • Synchronized signals and data
    • Combine with myriad of other devices
  • Replace Zeros
  • Replace biological outliers (200 to 3000 ms) and breaks
  • Replace statistical outliers (Z(0.025, 0.975))
  • Ectopic beat detection:
    • ACAR Acar et al. (2000)
    • Karlsson Karlsson et al. (2012)
    • Malik Malik (1996)
    • Kamath Kamath and Fallen (1995)
    • Custom absolute threshold (0.2)
  • Signal interpolation:
    • Constant
    • Linear
    • Quadratic
    • Cubic
  • Windowed analysis:
    • Sliding window
    • Overlapping windows
    • Non-overlapping windows

Raw data

Code 
f = Figure(size=(800, 600));
ax = Axis(f[1, 1], 
    title="Raw Data",
    xlabel="Time (s)",
    ylabel="Inter-Beat Interval (ms)",
    )
n = data[1:50]
x = Observable(n)
t = Observable(cumsum(n) ./ 1000) # to seconds
stem!(ax, t, x, color=:red)
ylims!(ax, (minimum(n) - 20, maximum(n) + 20))
f

Ectopic-beat removal

Code 
f = Figure(size=(800, 600));
ax = Axis(f[1, 1], 
    title="Raw Data",
    xlabel="Time (s)",
    ylabel="Inter-Beat Interval (ms)")
n = data[1:end]
x = Observable(n)
t = Observable(cumsum(n) ./ 1000) # to seconds
lines!(ax, t, x)
scatter!(ax, t, x, color=:red, markersize=5)

# Detect outliers
# xn = HeartRateLab.replace_ectopic_beats(n, method=:karlsson)
xn = HeartRateLab.replace_ectopic_beats(n, method=:acar, threshold=0.20)
x = Observable(xn)
scatter!(ax, t, x, color=:blue, markersize=5, label="Ectopic Beats Removed")

f

Feature Extraction

1

Open the full vertical graph

Rolling Window Features

Code 
ds_windowed = HeartRateLab.Features.windowed_feature_set(
    data, features=["vlf", "rmssd"], window_size=60, stride=10
)

f = Figure(size=(1200, 600));
ax = Axis(f[1, 1], 
    title="RMSSD through Time WS=60, S=10",
    xlabel="Time (s)",
    ylabel="RMSSD (ms²)")
x = ds_windowed.rmssd
t = 1:length(x)
lines!(ax, t, x, color=:blue)
f
Extracted 29 windows with 2 features each.

With different resolutions

Code 
ds_windowed = HeartRateLab.Features.windowed_feature_set(
    data, features=["vlf", "rmssd"], window_size=30, stride=1
)

ax = Axis(f[2, 1], 
    title="RMSSD through Time WS=30, S=1",
    xlabel="Time (s)",
    ylabel="RMSSD (ms²)")
x = ds_windowed.rmssd
t = 1:length(x)
lines!(ax, t, x, color=:blue)
f
Extracted 312 windows with 2 features each.
Code 
ds_windowed = HeartRateLab.Features.windowed_feature_set(
    data, features=["vlf", "rmssd"], window_size=120, stride=5
)

ax = Axis(f[3, 1], 
    title="RMSSD through Time WS=120, S=5",
    xlabel="Time (s)",
    ylabel="RMSSD (ms²)")
x = ds_windowed.rmssd
t = 1:length(x)
lines!(ax, t, x, color=:blue)
f
Extracted 45 windows with 2 features each.

Windowed Frequency Analysis

Code 
using HeartRateLab.Frequency: lomb_scargle
ds_windowed = HeartRateLab.windowed(data, window_size=60, stride=5)
freq = lomb_scargle(ds_windowed[1]).freq
y_values = 1:length(ds_windowed)
X = []
min_size = Inf
for i in 1:length(ds_windowed)
    w = ds_windowed[i]
    y = y_values[i]
    pgram = lomb_scargle(w)
    x = pgram.freq
    z = pgram.power
    push!(X, z)
    min_size = min(min_size, length(x))
end
# Crop every element to the minimum size
println("Minimum size: ", min_size)
println("Minimum in X: ", minimum([length(x) for x in X]))
for i in 1:length(X)
    X[i] = X[i][1:Int(min_size)]  # Ensure all have the same length
end
X = hcat(X...) # Windows by Frequency
y_values = 1:size(X, 2) # Window Index
freq = freq[1:Int(min_size)] # Crop frequency to the minimum size
f = Figure(size=(1600, 800));
ax = Axis3(f[1, 1], 
    title="Windowed Power Spectral Density",
    xlabel="Window Index",
    ylabel="Frequency (Hz)",
    zlabel="Power (ms²/Hz)",
    # limits=((0, length(ds_windowed)), (0, 0.5), (0, nothing)),
)
for i in 1:size(X, 2)
    x = freq # Frequency
    y = y_values[i] # Window Index
    z = X[:, i]
    points = Point3.(x, y, z)  # Data to Point3
    base   = Point3.(x, y, minimum(z))
    band!(ax, base, points, alpha=0.5, color=Makie.wong_colors()[i % length(Makie.wong_colors()) + 1])
    lines!(ax, points, linewidth=2, color=Makie.wong_colors()[i % length(Makie.wong_colors()) + 1])
end
f
Minimum size: 86.0
Minimum in X: 86

Statistical Analysis of single recordings

Code 
feature_registry = HeartRateLab.Features.feature_registry
feature_set = String.(keys(feature_registry))
println("Feature set: ")
for f in feature_set
    println(" - ", f)
end
Feature set: 
 - sd2_sd1
 - range
 - hf_percentage
 - mean
 - sdsd
 - lf
 - mean_hr
 - tp
 - renyi2
 - hf_peak
 - tinn
 - sdnn
 - max
 - sd1_sd2_area
 - lf_relative
 - rRR
 - hurst
 - renyi1
 - pnn50
 - pnn20
 - dfa2
 - cvsd
 - sdann
 - lf_hf_ratio
 - apen
 - sampen
 - lf_percentage
 - median
 - sd1
 - hf
 - ccsi
 - min_hr
 - renyi0
 - rmssd
 - std_hr
 - triangular_index
 - hf_relative
 - min
 - cvi
 - max_hr
 - lf_peak
 - sd2
 - vlf
 - ulf
Code 
ds_windowed = HeartRateLab.Features.windowed_feature_set(
    data, features=["rmssd", "vlf"], window_size=60, stride=10
)
f = Figure(size=(600, 600));
ax = Axis(f[1, 1],
    title="Windowed RMSSD histogram",
    xlabel="Time (s)",
    ylabel="Feature Value",
)
x = ds_windowed.rmssd
# Confidence intervals
mean_x = mean(x)
std_x = std(x)
ci_lower = StatsBase.quantile(x, 0.025)
ci_upper = StatsBase.quantile(x, 0.975)
density!(ax, x, label="Density", color=(:coral, 0.5))
vlines!(ax, mean_x, color=:blue, linewidth=2, label="Mean")
vlines!(ax, ci_lower, color=:green, linewidth=2, linestyle=:dash, label="95% CI Lower")
vlines!(ax, ci_upper, color=:green, linewidth=2, linestyle=:dash, label="95% CI Upper")
# hist!(ax, x, bins=30, label="RMSSD")
f[1, 2] = Legend(f, ax, "Windowed RMSSD Histogram",
    position=:bottomright, 
    title="Windowed RMSSD Histogram Legend",
    fontsize=10,
    bgcolor=:white,
    framecolor=:black)  
f
Extracted 29 windows with 2 features each.

Dynamic Systems

Generative models of the IBI series.

  • Dynamic Mode Decomposition - data-driven modal reconstruction
  • VanDerPol - nonlinear relaxation oscillator
  • Lorenz - deterministic chaos
  • LIF - leaky integrate-and-fire pacemaker

Dynamic Mode Decomposition

DMD on our IBI series - mode spectrum (top) and reconstruction against the real data (bottom). A low-rank linear model recovers the mean and the dominant low-frequency oscillation.

Van der Pol

Van der Pol phase portrait - vector field with trajectories converging to the limit cycle (red); mu = 0.63 fitted to our IBI data.

Van der Pol - animated

Van der Pol relaxation oscillator - the trajectory (red dot) spirals onto the limit cycle in phase space, with X(t), Y(t) and the resulting inter-beat intervals read out over time.

Lorenz

Lorenz system (x-z): the trajectory wanders the two wings; every peak in z (red) is read out as a heartbeat, and the peak-to-peak gaps become inter-beat intervals (right). σ=10, ρ=28, β=8/3.

Leaky Integrate-and-Fire Model (LIF)

Artificial LIF example - sub-threshold membrane integration to threshold, then reset (plot_lif, our Visualization stack).

Code 
using DifferentialEquations
using Random

# Parameters for the HRV model
params = (
    τ=0.6,        # Integration time constant
    I_base=1.1,   # Baseline input to the integrator
    threshold=0.9, # Firing threshold
    noise_amp=0.15,  # Amplitude of stochastic noise
)

# Define the integrate-and-fire dynamics
function hrv_dynamics(dV, V, p, t)
    dV[1] = (p.I_base - V[1]) / p.τ + p.noise_amp * randn()
end

# Time span for simulation
tspan = (0.0, 1000.0)  # Simulate for 100 seconds
V0 = [0.0]            # Initial integrator value

# Define the callback to reset the integrator and log spike times
function hrv_fire!(integrator)
    # Log the firing time
    global spike_times
    push!(spike_times, integrator.t)
    # Reset the integrator value
    integrator.u[1] = 0.0
end

# Condition for firing
function fire_condition(u, t, integrator)
    u[1] >= params.threshold
end

# Initialize a global spike time array
global spike_times = Float64[]

# Create a DiscreteCallback for firing events
fire_callback = DiscreteCallback(fire_condition, hrv_fire!)

# Solve the integrate-and-fire HRV model
hrv_problem = ODEProblem(hrv_dynamics, V0, tspan, params)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1000.0)
u0: 1-element Vector{Float64}:
 0.0

LIF Input Current

The LIF driven by our IBI series: the input current reconstructed from the recording feeds the leaky integrator, whose threshold crossings become beats.

Solving ODEs

Code 
hrv_solution = solve(hrv_problem, Tsit5(); callback=fire_callback)

# Interbeat intervals (IBIs) in milliseconds
IBIs = diff(spike_times) * 1000.0
IBIs = HeartRateLab.replace_bio_outliers(IBIs)

# Plot the results
f = Figure(size=(1600, 430));
ax = Axis(f[1, 1],
    title="Simulated Heart Rate",
    xlabel="Time (s)",
    ylabel="Heart Rate (bpm)",
)
# stem!(ax, spike_times[1:(end - 1)], IBIs,
#     marker=:circle, color=:blue, markersize=5)
t = cumsum(IBIs) ./ 1000 # Convert to seconds
x = 60000 ./ IBIs # Convert to bpm
lines!(ax, t, x, label="Simulated Heart Rate")
f

Bayesian Inference and Statistical Modeling

Pictorial proof of Bayes’ theorem By Cmglee - Own work, CC BY-SA 4.0, https://commons.wikimedia.org/w/index.php?curid=143581610

Where we actually use it: fitting a Bayesian AR model to the RR-interval series for forecasting. Unlike the LIF (which has a closed-form inter-spike interval), the AR coefficients have no analytical solution, so we estimate a posterior over the mean level, the innovation scale, and the AR weights.

using Turing

# Bayesian AR(p) model of the RR-interval series, used for RR forecasting.
# (from test/tools/forecasting/transfer_and_bayes.jl)
@model function bayes_ar(y, p)
    n = length(y)
    mu    ~ Normal(mean(y), 200)         # mean RR level (ms)
    sigma ~ Exponential(50)              # innovation scale (ms)
    phi   ~ filldist(Normal(0, 0.5), p)  # AR coefficients
    for t in (p + 1):n
        m = mu
        for j in 1:p
            m += phi[j] * (y[t - j] - mu)
        end
        y[t] ~ Normal(m, sigma)
    end
end

# Fit to a recording's RR series, then draw the posterior predictive band
chain = sample(bayes_ar(rr, 8), NUTS(0.65), 250)

AR forecast of the RR series

AR(8) forecasts of the same RR series at beats 133, 230, 326 (120-beat train, 10-beat forecast); the 90% predictive band is what we score for calibration.

Visualization

All visualizations in this notebook can be generated live.

Code 
N = 10000
x = Observable(rand(N))
fig = Figure(size=(1800, 600));
ax = Axis(fig[1, 1], subtitle="Random Variable")
scatter!(ax, x)
ylims!(ax, -1, 2)
ax = Axis(fig[1, 2], subtitle="Distribution")
# density!(ax, x)
hist!(ax, x, bins=100)
xlims!(ax, -1, 2)
button = Button(fig[2, 1:2], label="New")
l = on(button.clicks) do b
    x[] = rand(N)
end
fig

DEMO

Problem Solving with the HeartRateLab

What have we used the HeartRateLab for?

Teaching, priors, and modeling

Teaching and Learning

  • Communication and Visualization
  • Test-driven feature definition
  • Ranges and normal values of features are part of our tests
  • Open source encourages student engagement
  • Online tools help us debug live during our experiments
  • Online demos aid visual understanding
  • Personal Biofeedback practice

Managing Expectations and Priors

  • Open Data Sets become the basis for our Normative Datasets
  • For every study, a graph of related studies can be selected from a superset to compare and contrast the results.
  • Given similar studies, we can design experiments with quantified expectations (Bayesian priors).
  • We can identify outliers in different populations, and under different conditions.

Modeling and Predicting

  • Models can be tested against normative datasets.
  • Their accuracy and precision can be quantified.
  • Artificial models encourage model-thinking: no true model, but useful models.

What do we do with the HeartRateLab?

  • Arousal and Relaxation
  • Cognitive Mental Workload
  • Attention and Dual-Tasking
  • Attentional Networks
  • Facial Expression Observation
  • Ocular Rivalry and Conscious Perception
  • Cognitive Flexibility
  • Heart Rate Variability Biofeedback
  • Mindfulness and Contemplative Practices
  • Neurodiversity, ADHD, and Autism

Scientific Infrastructure behind HeartRateLab

Heterogenous Distributed Computing

The HeartRateLab was developed and is currently used under a self-designed research infrastructure that enables metascientific evaluation of the research field:

  • The Wanderers Science-OS: A custom made NixOS system for a fleet of heterogenous devices that includes diferent forms of computing in diverse environments.
  • The Science Agents agentic science framework: A sistem that includes research field mapping, bibliography and bibilometry, researcher and topic tracikng, conference and CfP monitor, project management, and general scientific supervision.
  • Explicit Knowledge Management systems over three shared memory layers on memvault, running on the Wanderers fleet.

Science Agents on the Wanderers

The Wanderers Research-OS: four repos (ScienceAgents, Supervision, ScienceMap, Bibliometry) over three memory layers (Identity in git; Episodic and Semantic-graph on memvault)

The four systems

  • ScienceAgents - agent infrastructure: shared identity and methodology.
  • Supervision - the AI research supervisor: critical-path tracking and weekly / Sunday / CFP / SWOT reports, currently driving the TU Graz PhD.
  • ScienceMap - the world model: a Scientific Entity Graph (researchers, venues, topics, opportunities) on semantic memory, with crawlers, tracking and visualization.
  • Bibliometry - Zotero ingest + Scholar enrichment; emits entity/edge tuples into ScienceMap. Uses AutoResearchClaw for literature search and 4-layer citation verification.

Three memory layers on memvault - Identity (git), Episodic, and the Semantic graph - replicated across the fleet.

On the Wanderers fleet: the agents run on Saturn; Moon and Phobos host the services (git, memvault, entity graph); Neptune and Jupiter provide GPU inference and fast builds.

AutoResearchClaw (ARC)

AutoResearchClaw’s 23-stage pipeline, which Bibliometry uses for literature search (OpenAlex, Semantic Scholar, arXiv) and 4-layer citation verification. Its broad search pulls thousands of candidates that the screening gate narrows to a confirmed shortlist; runs as Slurm jobs against a local models on available GPUs.

The HRV field over time

The HeartRateLab knowledge base as a field: 554 HRV papers (154 confirmed references plus the HRV papers they cite), by publication year and coloured by research field - clinical, sports & peak-performance, contemplative practice, and methods & foundations. Seminal signal-processing and physiology work (Shannon 1948, Welch 1967, Akselrod 1981) anchors the early years.

Research field networks

554 papers and 3,841 citations - 89% sit in one connected component, so the field is a single citation web, not scattered groups. Colour = research field; node size = how often a paper is cited within the set (the largest nodes are the shared foundations - Task Force 1996, Akselrod 1981, Peng/DFA 1995).

Tools and Resources

The cynefin framework Snowden and Boone (2007)

References

Acar, Burak, Irina Savelieva, Harry Hemingway, and Marek Malik. 2000. “Automatic Ectopic Beat Elimination in Short-Term Heart Rate Variability Measurement.” Computer Methods and Programs in Biomedicine 63 (2): 123–31.
Kamath, MV, and EL Fallen. 1995. “Correction of the Heart Rate Variability Signal for Ectopics and Missing Beats’ Heart Rate Variability. Malik m, Camm AJ, Eds.” Boston: Blackwell Publishing.
Karlsson, Marcus, Rolf Hörnsten, Annika Rydberg, and Urban Wiklund. 2012. “Automatic Filtering of Outliers in RR Intervals Before Analysis of Heart Rate Variability in Holter Recordings: A Comparison with Carefully Edited Data.” Biomedical Engineering Online 11: 1–12.
MacDonald, Eilidh A, Robert A Rose, and T Alexander Quinn. 2020. “Neurohumoral Control of Sinoatrial Node Activity and Heart Rate: Insight from Experimental Models and Findings from Humans.” Frontiers in Physiology 11: 170.
Malik, Marek. 1996. “Heart Rate Variability: Standards of Measurement, Physiological Interpretation, and Clinical Use: Task Force of the European Society of Cardiology and the North American Society for Pacing and Electrophysiology.” Annals of Noninvasive Electrocardiology 1 (2): 151–81.
Netter, F. H., and Ciba Pharmaceutical Products. 1972. The Ciba Collection of Medical Illustrations: Anatomy, Physiology, and Metabolic Disorders. Nervous System. Ciba Pharmaceutical Products. https://books.google.at/books?id=2qARAQAAMAAJ.
Snowden, David J, and Mary E Boone. 2007. “A Leader’s Framework for Decision Making.” Harvard Business Review 85 (11): 68.

Thank you

HeartRateLab ring logo

github.com/abcsds/HeartRateLab.jl

Alberto Barradas, JuliaCon 2026, Mainz, Germany