Non-Linear HRV Monitoring After LBBAP: An N-of-1 Framework for Physiological Pacing Quality Surveillance

Why ask the question this way

A common informal hypothesis among engaged patients and clinicians is that non-linear heart rate variability (HRV) might forecast LBBAP lead malfunction before it becomes electrically manifest. The premise needs unpacking before it can be operationalized.

LBBAP lead failure modes are mechanically and electrically analogous to those of any active-fixation transvenous lead: threshold rise, microdislodgement with progressive loss of conduction system capture (selective → non-selective → myocardial-only → loss), insulation breach, conductor fracture. None of these failure mechanisms communicates with the autonomic nervous system. There is no physiologic chain by which an impending capture-vector change at the helix-tissue interface would alter brainstem-mediated cardiac vagal tone before the electrical event itself occurs. HRV cannot forecast acute lead malfunction.

However, a related and tractable question exists: does sustained preservation of favorable non-linear HRV indices serve as a non-invasive surrogate for sustained physiological-pacing quality — that is, persistent non-selective LBB capture with preserved electromechanical synchrony? This reformulation inverts the causal arrow. We are no longer asking HRV to predict an electrical event upstream of itself. We are asking whether the autonomic-hemodynamic consequence of degraded resynchronization (e.g., a transition from non-selective LBB capture to myocardial-only capture, with re-introduction of dyssynchrony) leaves a detectable longitudinal signature in non-linear HRV metrics.

Mechanistic substrate

Non-selective LBB capture produces a paced QRS duration typically 100–120 ms with a right-bundle-like morphology in V1 (terminal R or r′) and a short V6 R-wave peak time (<80 ms). The hemodynamic correlate is near-physiological electrical and mechanical activation, with preserved myocardial work efficiency comparable to intrinsic conduction.

Loss of LBB capture with preservation of septal myocardial capture transitions the patient acutely to a wide, RV-pacing-equivalent QRS (typically 150–180 ms) with restored dyssynchrony. The hemodynamic consequences — altered baroreflex loading, transient stroke volume reduction, sympathetic activation — are precisely the inputs that drive non-linear HRV. This is the substrate the framework exploits.

The metrics

Phase-Rectified Signal Averaging Deceleration Capacity (PRSA-DC)

DC was introduced by Bauer et al. (Lancet 2006) to isolate the vagal component of heart rate fluctuation from a non-stationary RR series. In the original 1,455-patient post-MI cohort with blinded validation in London (n=656) and Oulu (n=600) cohorts, DC outperformed both LVEF and SDNN for mortality prediction. In heart failure populations specifically, DC achieves AUC 0.88 in men and 0.97 in women for classifying HF status, with cutoffs of approximately 4.55 ms and 4.85 ms respectively.

Algorithm

1. Anchor identification. Walk the RR series. An anchor is any RR(i) such that RR(i) > RR(i−T), with the fractional increase bounded (typically ≤5%) to exclude artifact and ectopic beats.
2. Segment extraction. Around each anchor, extract a window of 2L+1 beats (L beats before, L beats after). For short-term vagal DC, L = 2.
3. PRSA averaging. Align all extracted segments at their anchors and average beat-by-beat. Quasi-periodic oscillations that consistently align with deceleration anchors survive; uncorrelated fluctuations average toward zero. This is the phase-rectification step that gives PRSA its robustness to non-stationarity.
4. DC computation from the averaged PRSA curve:

where X(0) is the value at the anchor index in the averaged segment. Units are milliseconds; higher values indicate greater vagal deceleration capacity.

Population reference

For the post-LBBAP monitoring use case, population norms are largely irrelevant. The signal of interest is within-subject longitudinal trend against the patient’s own baseline.

Detrended Fluctuation Analysis short-term scaling exponent (DFA α1)

DFA quantifies the fractal correlation properties of the heartbeat time series. The short-term scaling exponent α1, calculated over windows of 4–16 beats, has shown independent prognostic value across post-MI, chronic HF, and HFpEF populations.

Algorithm

1. Integration. Compute the cumulative sum y(k) = Σ[RR(i) − mean(RR)] for i = 1 to k, producing an integrated profile.
2. Windowing and detrending. For each window size n in [4, 16], partition y(k) into non-overlapping windows of length n. In each window, fit a linear trend by least squares and subtract it.
3. Fluctuation function. F(n) = √[(1/N) Σ (y(k) − y_n(k))²], where y_n(k) is the locally detrended profile.
4. Scaling exponent. α1 is the slope of log F(n) vs log n in the 4–16 beat range.

Interpretation

Exercise physiology layer

DFA α1 has a useful dynamic range across exercise intensities. In healthy and trained populations, α1 ≈ 1.0 at rest, drops to ≈0.75 at the aerobic threshold (HRV-derived ventilatory threshold, "HRVT"), and reaches ≈0.5 at the anaerobic threshold. This creates a confounder for the monitoring use case — training state itself shifts α1 — but also creates an additional metric: the workload at which α1 crosses 0.75 (DFA-HRVT) should be stable in a patient with sustained physiological pacing quality at constant training state. An unexplained downward drift of DFA-HRVT, controlled for detraining, would be a candidate signal.

Post-exercise heart rate deceleration gradient (HRD)

An adjacent metric with direct CRT-response validation. Korosoglou and colleagues showed that in 37 patients receiving biventricular pacemakers, the gradient of heart rate decline 30 seconds after cessation of cardiopulmonary exercise testing predicted both functional response (AUC 0.87) and echocardiographic response (AUC 0.82) at 3 months. A cut-off of HRD30 ≥ 3 (≥5% HR drop from peak at 30 seconds) optimized the sensitivity/specificity trade-off.

The physiologic substrate is vagal reactivation capacity, which depends both on intrinsic autonomic competence and on the hemodynamic state being restored at exercise cessation. In the LBBAP monitoring context, this is one of the few metrics with a published cut-off in a CRT population and is easily extractable from any chest-strap RR recording of a structured exercise session.

Conceptual map: HRV behavior across capture states

The N-of-1 protocol

Hardware and software stack

Recording schedule

Ground-truth anchoring

Every scheduled device interrogation becomes a ground-truth anchor in the longitudinal dataset: paced QRS duration, V6 R-wave peak time, V1 morphology, capture threshold, lead impedance, sensed R-wave amplitude. The retrospective analysis at each interrogation asks two questions:

1. Concordance. Did any HRV metric track interrogation findings? Stable thresholds and morphology with stable HRV trends is the boring-but-desired state. Threshold drift coupled with deteriorating HRV is concordance worth investigating.
2. Lead/lag. When concordance occurs, does HRV change precede electrical drift? This is the properly framed forecasting question.

Decision logic

Triggers for action (12-lead ECG with V1 morphology check, earlier device interrogation) are defined as deviation from the patient's own 30-day rolling baseline:

• Resting morning DC dropping >2 SD below 30-day rolling mean, sustained for ≥7 consecutive days
• Resting α1 deviation >0.2 from baseline mean in either direction, sustained for ≥7 days
• HRD30 post-rowing-cessation falling below patient's established baseline by >1 SD over 4 consecutive sessions
• DFA-HRVT downward drift >10 W (or equivalent intensity unit) at constant training state

These thresholds are conservative starting points. After 60–90 days of personal baseline data, they should be tightened or loosened based on observed within-subject variability.

Reference implementation

A self-contained Python implementation of the core algorithms, suitable for batch analysis of Polar H10 exports or any RR series in milliseconds. Dependencies: NumPy, Pandas. The full module (including longitudinal deviation detection and a demonstration block) is available on www.abcfarma.net.

"""
LBBAP Longitudinal HRV Monitoring
==================================
Reference implementation of:
  - Phase-Rectified Signal Averaging Deceleration Capacity (PRSA-DC)
    Bauer et al., Lancet 2006; 367:1674-1681
  - Detrended Fluctuation Analysis short-term scaling exponent (DFA alpha1)
    Peng et al., Chaos 1995; 5:82-87

Intended use: within-subject longitudinal monitoring of autonomic /
fractal HRV indices as a non-invasive surrogate for sustained
physiological pacing quality in a DDD-LBBAP system. NOT a substitute
for direct device interrogation.
"""

import numpy as np
import pandas as pd


# -----------------------------------------------------------
# RR preprocessing
# -----------------------------------------------------------

def clean_rr(rr_ms, low=300, high=2000, max_pct_change=0.20):
    """Artifact rejection: physiologic bounds + successive %-change filter."""
    rr = np.asarray(rr_ms, dtype=float)
    rr = rr[(rr >= low) & (rr <= high)]
    diffs = np.abs(np.diff(rr) / rr[:-1])
    keep = np.concatenate([[True], diffs < max_pct_change])
    return rr[keep]


# -----------------------------------------------------------
# PRSA Deceleration Capacity
# -----------------------------------------------------------

def deceleration_capacity(rr_ms, L=2, T=1, max_pct=0.05):
    """
    Bauer's PRSA Deceleration Capacity.

    Returns
    -------
    dc : float (ms)
        Population HF cutoffs: ~4.55 (M), ~4.85 (F)
        Severe post-MI risk: <2.5
        Trained endurance athletes typically 8-15
    n_anchors : int
    """
    rr = np.asarray(rr_ms, dtype=float)
    N = len(rr)

    # Identify deceleration anchors with bounded fractional increase
    anchors = []
    for i in range(T + L, N - L):
        if rr[i] > rr[i - T] and (rr[i] / rr[i - T] - 1) <= max_pct:
            anchors.append(i)

    if len(anchors) < 10:
        return np.nan, len(anchors)

    # Build PRSA signal: average aligned 2L+1 windows
    window_len = 2 * L + 1
    aligned = np.zeros(window_len)
    for a in anchors:
        aligned += rr[a - L : a + L + 1]
    prsa = aligned / len(anchors)

    # DC formula: ([X(0) + X(1)] - [X(-1) + X(-2)]) / 4
    X0, X1, Xm1, Xm2 = prsa[L], prsa[L + 1], prsa[L - 1], prsa[L - 2]
    dc = ((X0 + X1) - (Xm1 + Xm2)) / 4.0
    return dc, len(anchors)


# -----------------------------------------------------------
# Detrended Fluctuation Analysis
# -----------------------------------------------------------

def dfa(rr_ms, scales=None, order=1):
    """DFA. Returns (alpha, scales, F)."""
    rr = np.asarray(rr_ms, dtype=float)
    if scales is None:
        scales = np.arange(4, 17)
    scales = np.asarray(scales, dtype=int)

    # Step 1: integrate
    y = np.cumsum(rr - np.mean(rr))
    N = len(y)

    F = np.zeros(len(scales))
    for idx, n in enumerate(scales):
        n_windows = N // n
        if n_windows < 2:
            F[idx] = np.nan
            continue
        segments = y[: n_windows * n].reshape(n_windows, n)
        x_axis = np.arange(n)
        rms_vals = []
        for seg in segments:
            coeffs = np.polyfit(x_axis, seg, order)
            trend = np.polyval(coeffs, x_axis)
            rms_vals.append(np.sqrt(np.mean((seg - trend) ** 2)))
        F[idx] = np.mean(rms_vals)

    valid = ~np.isnan(F) & (F > 0)
    if valid.sum() < 3:
        return np.nan, scales, F
    log_n = np.log10(scales[valid])
    log_F = np.log10(F[valid])
    alpha, _ = np.polyfit(log_n, log_F, 1)
    return alpha, scales, F


def dfa_alpha1(rr_ms):
    """Short-term DFA scaling exponent (windows 4-16 beats)."""
    alpha, _, _ = dfa(rr_ms, scales=np.arange(4, 17), order=1)
    return alpha


# -----------------------------------------------------------
# Sliding-window alpha1 (exercise threshold detection)
# -----------------------------------------------------------

def sliding_dfa_alpha1(rr_ms, window_beats=120, step_beats=30):
    """Sliding alpha1 across a session. Useful for DFA-HRVT detection."""
    rr = np.asarray(rr_ms, dtype=float)
    cumulative_s = np.cumsum(rr) / 1000.0
    rows = []
    for start in range(0, len(rr) - window_beats, step_beats):
        window = rr[start : start + window_beats]
        alpha = dfa_alpha1(window)
        center = start + window_beats // 2
        rows.append({
            "center_beat": center,
            "center_time_s": cumulative_s[center],
            "alpha1": alpha,
            "mean_hr": 60000.0 / np.mean(window),
        })
    return pd.DataFrame(rows)


# -----------------------------------------------------------
# Post-exercise HR deceleration gradient (CRT-response metric)
# -----------------------------------------------------------

def hr_deceleration_gradient(ex_rr, rec_rr, peak_hr=None,
                              sample_points_s=(30, 60, 90, 120)):
    """
    HRD30/60/90/120 -- Korosoglou CRT-response metric.
    Cutoff: HRD30 >= 5% associated with positive CRT response (AUC 0.87).
    """
    ex_rr = np.asarray(ex_rr, dtype=float)
    rec_rr = np.asarray(rec_rr, dtype=float)
    if peak_hr is None:
        last_beats = ex_rr[-10:] if len(ex_rr) >= 10 else ex_rr
        peak_hr = 60000.0 / np.min(last_beats)

    cumulative_s = np.cumsum(rec_rr) / 1000.0
    results = {"peak_hr": peak_hr}
    for tgt in sample_points_s:
        idx = np.searchsorted(cumulative_s, tgt)
        if idx >= len(rec_rr):
            results[f"HRD_{tgt}"] = np.nan
            continue
        lo, hi = max(0, idx - 2), min(len(rec_rr), idx + 3)
        hr_at_t = 60000.0 / np.mean(rec_rr[lo:hi])
        results[f"hr_at_{tgt}s"] = hr_at_t
        results[f"HRD_{tgt}"] = (peak_hr - hr_at_t) / peak_hr * 100.0
    return results


# -----------------------------------------------------------
# Longitudinal deviation detection
# -----------------------------------------------------------

def detect_deviation(history_df, current_value, metric,
                     window=30, sd_threshold=2.0):
    """Flag deviation from rolling baseline."""
    recent = history_df.tail(window)[metric].dropna()
    if len(recent) < 10:
        return {"flag": False, "reason": "insufficient baseline"}
    baseline_mean = recent.mean()
    baseline_sd = recent.std(ddof=1)
    z = (current_value - baseline_mean) / baseline_sd if baseline_sd > 0 else 0
    return {
        "metric": metric,
        "current": current_value,
        "baseline_mean": round(baseline_mean, 3),
        "baseline_sd": round(baseline_sd, 3),
        "z_score": round(z, 2),
        "flag": abs(z) >= sd_threshold,
        "direction": "below" if z < 0 else "above",
    }

Limits of evidence

This framework is hypothesis-generating. The component metrics — PRSA-DC and DFA α1 — have robust prognostic validation in post-MI and chronic HF populations, and HRD30 has been validated as a CRT-response marker (Korosoglou et al.). However, there is no prospective clinical trial of non-linear HRV as a longitudinal surveillance tool specifically in LBBAP or conduction system pacing populations. Clinicians and engaged patients using this framework are operating at the frontier, not following established guidelines. Findings should be interpreted as triggers for direct electrical verification, never as definitive evidence of capture state.

FAQ

References