BAND: A Unified Probabilistic Framework for Non-Stationary Heart Rate Variability Analysis

1. Introduction

A major goal of contemporary biomedical science is to mechanistically connect the dynamic behavior of physiological systems to outcomes in human health and disease. The primary obstacle to this goal is the inferential gap between what can be measured with high-temporal resolution via wearable sensors and the underlying neurobiological processes that generate these measurements[1]. This inverse problem is particularly evident in the study of the autonomic nervous system (ANS), where heart rate variability (HRV) has become the dominant non-invasive marker[2,3,4]. The R-R interval (RRi) time-series is a highly informative signal, yet it presents a substantial analytical challenge, as it reflects the integrated and dynamic interplay of multiple regulatory systems[5]. Addressing this challenge requires formal methodologies capable of deconstructing the composite signal into quantitative, interpretable, and testable estimates of its underlying generative processes. The urgency of this task is highlighted by the central role of autonomic dysregulation in conditions ranging from cardiovascular diseases[6,7] to anxiety and post-traumatic stress[8,9,10], where precise characterization of system dynamics is essential for diagnosis and treatment.

The analytical tools most widely employed in this field remain poorly suited to the task, as the prevailing paradigm relies on calculating summary statistics within fixed temporal windows, thereby imposing a static and linear worldview on a system that is inherently dynamic and non-linear. These methods assume stationarity, an assumption routinely disrupted during physiological transitions of scientific or clinical interest[11,12]. As a result, they can produce misleading or uninterpretable outcomes. For example, by averaging across stress recovery, such methods conflate the initial sympathetic surge with the subsequent parasympathetic rebound, yielding a single opaque number that obscures the system’s adaptive capacity[11,12]. In this way, the unavoidable distortion of RRi signals under non-stationary conditions may mask critical differences between a resilient individual, characterized by rapid recovery, and a vulnerable one, marked by dampened recovery, thereby concealing an essential biomarker of pathology[13]. As a result, the “discretize-and-summarize” approach discards the temporal structure of the data and precludes direct estimation of core mechanistic properties such as response latencies, rates of change, and adaptive magnitudes, which represent the essential features that define the physiological process[14].

More sophisticated approaches have been proposed, but they too often fail to resolve the core inferential challenge, presenting instead a false dichotomy between phenomenological description and non-interpretable abstraction. On one hand, time–frequency decompositions such as the wavelet transform provide detailed visualizations of spectral dynamics, making them valuable for exploratory analysis[15]. Yet they remain fundamentally descriptive, unable to instantiate a generative model of the underlying biology. These methods can demonstrate that a change has occurred, yet they are unable to provide a parameterized and testable explanation of the underlying nature of that change[16]. On the other hand, modeling traditions like state-space models (SSMs) offer mathematical sophistication for tracking and forecasting system evolution[17,18]. While powerful, these models frequently sacrifice physiological fidelity for mathematical convenience, as their latent states are often defined by statistical rather than biological properties[18,19]. Even in tailored applications, SSMs excel at tracking continuous variation but are not explicitly structured to parameterize the discrete phases of physiological events, such as onset rate, recovery latency, or adaptive magnitude in response to stress[18,19]. The field thus remains methodologically constrained, caught between descriptive tools that cannot test mechanisms and abstract models whose parameters do not map onto physiological hypotheses.

Bridging this gap requires a conceptual shift from model-agnostic data description to the construction of hypothesis-driven frameworks grounded in generative modeling. Instead of simply characterizing what the data looks like, the focus should be on constructing explicit, testable, and quantitative models of the biological processes that generate the data. Here, we introduce and validate the Biphasic Autonomic Non-stationary Decomposition (BAND) framework, which reframes physiological analysis as a form of formal hypothesis testing. We hypothesize that autonomic responses to acute perturbations constitute structured, dynamic processes that can be effectively modeled as continuous-time stochastic process. Specifically, we employ a double-logistic function as a parsimonious hypothesis for biphasic dynamics, capturing the characteristic onset and recovery phases[20]. This representation translates raw signals into a small set of interpretable mechanistic parameters, including the latency, rate, and magnitude of each phase, thereby providing direct insight into autonomic reactivity and resilience.

2. Results

We propose a framework based on a generative Bayesian model that yields robust parameter estimates with full uncertainty quantification, supports hierarchical inference across individuals and groups, and enables formal model comparison. The model decomposes the RRi series into a time-varying mean, total variance, and a dynamically shifting oscillatory component controlled by a double-logistic latent dynamic. This formulation provides a more accurate and physiologically interpretable representation of autonomic dynamics than conventional fixed-window methods, offering quantitative indices of reactivity and recovery and a principled means to translate complex time series into mechanistic, falsifiable knowledge.

After its construction, we test the BAND model’s internal consistency by assessing its ability to recover its own parameters from simulated data. We then benchmark its descriptive fidelity against conventional methods to quantify its effectiveness in characterizing dynamic signals. Finally, we use the framework to generate a structured, phenomenological description of an empirical exercise-response signal, illustrating its potential for developing precise and data-driven hypotheses.

The core components of the model, related the underlying data generation process being modeled, are illustrated in Figure 1.

2.1. Model Validation and Parameter Recovery on Synthetic Data

To establish the model’s internal consistency and the identifiability of its parameters, we first assessed its ability to perform inference on synthetic data generated from its own structure. This in silico procedure serves as a necessary check of the model’s self-consistency and the estimation algorithm’s ability to navigate the complex posterior geometry. To test the model’s performance, three distinct scenarios were designed: (A) a canonical response with complete recovery; (B) an incomplete recovery with persistence in both time and frequency domains; and (C) a high-noise condition with a low structured variance fraction ($w$ = 0.6). The ground-truth dynamics are illustrated in Figure 2, with corresponding generative parameter values listed in Table 1. This procedure tests for inferential integrity under the assumption that the model is correctly specified; it does not test for robustness to model misspecification, a challenge inherent to all empirical data.

When fit to this synthetic data, the model demonstrated a high-fidelity reconstruction of the generative process across all three scenarios. Figure 3 presents a visual summary of these results, plotting the model’s posterior estimates against the known ground-truth trajectories. The posterior means for the baseline RRi, total SDNN, and the time-varying spectral proportions closely track the true underlying dynamics. The 95% credible intervals, which represent the range of plausible parameter values given the data, consistently envelop the ground-truth trajectories.

Beyond recovering the continuous dynamic trajectories, the model also recovered the discrete parameters that govern the shape and timing of these dynamics. Table 2 presents the posterior summaries for these key parameters, comparing the posterior median and 95% credible intervals against their true values. For all parameters across all scenarios, the 95% credible intervals successfully contained the ground-truth values. Successfully recovering all parameters confirms that they are practically identifiable from this type of data, at least when the data-generating process aligns with the model’s structural assumptions.

Benchmarking the Model’s Descriptive Fidelity

Having established the model’s internal consistency, we next benchmarked its descriptive performance against conventional, widely adopted analytical methods. The objective of this comparison was to quantify the core differences in how these analytical frameworks represent and summarize dynamic signals, and to highlight the specific conditions under which our proposed model offers a more advantageous description and performance. The same synthetic datasets from the previous section were analyzed using a 60-second sliding-window analysis for time-domain metrics (mean RRi and SDNN), and a Short-Time Fourier Transform (STFT) with a similar window length for characterizing spectral dynamics. The STFT, and similar Fourier methods, allows to decompose the signal into frequency bands that compose the complete signal, which helps in the identification of the relative contribution of major frequency bands (which BAND model mechanistically performs as part of its architecture).

The results of this comparative analysis highlight the descriptive distortions that are inherent to any method that relies on assumptions of local stationarity when applied to a continuously varying, non-stationary signal. Figure 4 illustrates the outputs from the windowed methods. Visually, these methods produced a delayed and morphologically blunted approximation of the underlying dynamics. The sharp, non-linear transitions that characterize the onset and recovery phases of the signal are smeared and flattened by the averaging process inherent in the windowed calculation. Consequently, the windowed approach failed to accurately capture both the shape and the precise timing of the recovery process. Instead of reflecting the smooth, logistic curve defined by the ground-truth process, the windowed analysis superimposed the square-wave artifacts of its own computational structure onto the data, a clear instance of the analytical method distorting the phenomenon it is intended to measure.

A formal quantitative analysis confirms the descriptive efficiency and accuracy of the generative approach for signals of this particular nature. Table 3 provides a summary of key performance metrics. Across all scenarios and all signal components, the generative model consistently provided a more accurate and parsimonious description of the data. It achieved substantially lower error (RMSE and Bias) and captured a significantly higher proportion of the signal’s variance, with $R^2$ values typically exceeding 0.90. It is important to interpret this result with caution. The model’s performance is, in one sense, expected, as the synthetic data were generated from its own structure. However, it illustrates a core theoretical argument of this work: for signals that are plausibly characterized by smooth, continuous transitions between distinct physiological states, a model whose parametric form explicitly encodes this structure will naturally provide a more efficient, less distorted, and more interpretable quantitative summary than a non-parametric, windowed approach that was designed for stationary or quasi-stationary signals.

Further model diagnostics testing the model’s overall stability and reliability on synthetic data for in silico validation can be found in the supplementary material. Prior predictive checks are shown in Figure S1, prior sensitivity analysis in Figure S2, full posterior distributions in Figure S3, and model fit to asymmetric recovery data from a different data generation process in Figure S4. Additionally, non-Gaussian distributions are presented in Figure S5.

2.2. A Phenomenological Description of an Empirical Exercise Response

To illustrate the framework’s intended application in a real-world context, we analyzed an RRi time series recorded from a single healthy subject undergoing a standardized two-minute submaximal step test. The express purpose of this N-of-1 analysis is to demonstrate how the model translates a complex, high-dimensional, and noisy time series into a concise, low-dimensional, and interpretable set of parameters that form the basis for generating new, quantitative, and ultimately falsifiable hypotheses.

Upon fitting the model to the empirical data, it yields a phenomenological decomposition of the entire exercise-recovery response, as depicted in Figure 5. The framework parses the signal into its constituent components, describing the physiological process as a rapid initial tachycardia (a decrease in the baseline RRi component, panel A and B in Figure 5) and a concurrent suppression of total heart rate variability (a decrease in the SDNN component, panel C in Figure 5). This initial response is followed by a slower, more gradual recovery phase after the cessation of the exercise stimulus. From a spectral perspective, the model attributes the signature of this response to a near-complete withdrawal of power in the high-frequency (HF) band, and a corresponding surge in the proportion of very-low-frequency (VLF) power (panel D in Figure 5). The model captures these dynamics as continuous, smooth trajectories, providing an integrated narrative of the entire physiological event.

The framework’s primary scientific utility, however, lies in generating quantitative, testable claims, which are encapsulated by the posterior distributions of its parameters. A key feature of the model is its separate parameterization of the recovery dynamics of the signal’s first moment (mean RRi, via the coefficient $c_r$) and its second moment (total variability, or SDNN, via the coefficient $c_s$). This explicit separation allows for a direct, quantitative description of their potential divergence. In this particular recording, the posterior distributions for these two recovery parameters were credibly distinct; the 95% credible interval for $c_r$ (median of the posterior distribution = 0.69, CI_95%[0.66, 0.73]) was clearly separated from that for $c_s$ (median = 0.97, CI_95%[0.91, 1.00]). This descriptive feature of the model’s output, which we term model-inferred dissonance, provides a data-driven and falsifiable hypothesis: that the recovery of mean heart rate and the recovery of its variability can follow distinct and dissociable time courses in response to an acute exercise stressor. This hypothesis is a quantitative statement, grounded in the posterior probability distributions of the model’s parameters.

Empirical and visual model diagnostics can be seen in the supplementary material. Complete traceplots, posterior distributions and MCMC diagnostics can be seen in Figure S6, S7 and Table S1, respectively. Posterior and distributional predictive checks can be seen in Figure S8 and S9.

3. Discussion

The quantitative analysis of HRV has long been a primary method for the non-invasive study of the autonomic nervous system[2,21]. A persistent challenge in this field is the inherent non-stationarity of physiological signals, which directly conflicts with the foundational assumptions of traditional analytical techniques[11,12,22], which were often adopted due to the computational constraints of their era. This discrepancy necessitates the use of methodological compromises, most notably time-windowing, which can introduce significant analytical artifacts and obscure the very dynamics they seek to quantify[13,23]. The present study was undertaken to develop and evaluate a generative probabilistic model for RRi time series, explicitly designed to represent and infer the parameters of transient physiological events. Our findings, based on a comprehensive suite of simulations and an empirical application, suggest that the BAND framework offers a more accurate and granular characterization of autonomic dynamics compared to what conventional windowed approaches can achieve. Our generative model provides a more detailed characterization of autonomic dynamics; however, its interpretation requires careful consideration of its underlying assumptions.

The primary contribution of this research is the synthesis and application of established Bayesian inference principles to construct a purpose-built, physiologically oriented measurement instrument[24]. The model formalizes a specific, falsifiable hypothesis about the data-generating process underlying dynamic RRi signals. We employ the term generative to signify this complete probabilistic specification, which enables a coherent inferential procedure regarding the model’s latent parameters. The internal validity and statistical integrity of this inferential process were examined through extensive simulation studies. In these controlled in silico experiments, where the ground truth was known, the model consistently recovered the generative parameters, and its posterior credible intervals were well-calibrated[25,26]. The model’s successful parameter recovery in these controlled simulations provides baseline confidence in its function as a reliable measurement device, but only under conditions where its numerous and restrictive structural assumptions are perfectly met. When benchmarked against conventional methods, its superior performance in tracking dynamic signals was observed. This outperformance is an expected consequence of its design, which replaces the crude filtering of windowing with a formal, continuous model of the underlying process. However, the apparent sophistication of this approach belies a series of deep-seated assumptions and conceptual challenges that demand scrupulous consideration.

3.1. Methodological Superiority and the Abandonment of Stationarity

A core assertion of this work is that the analytical shortcomings of traditional HRV analysis are conceptual, not merely technical. Windowed methods are, in essence, digital filters[27]. As such, they inevitably distort the signal they are intended to measure, introducing lag, attenuating amplitudes, and enforcing a false dichotomy between temporal and spectral resolution as dictated by the Heisenberg-Gabor uncertainty principle[28]. The inferential approach adopted here avoids these specific artifacts by modeling the entire time-course of an event holistically. Information from every data point is leveraged to inform the estimate at every other data point, conditioned on the assumed model structure[29].

However, this approach introduces its own set of abstractions. A potentially problematic assumption is the model’s imposition of orthogonality on a non-orthogonal system. The model specifies the observed RRi signal as a linear, additive combination of a baseline trend, several sinusoidal components, and a residual noise term. This mathematical decomposition is elegant and ensures parameter identifiability[30,31]; nonetheless, it is a convenient simplification of cardiovascular physiology.

Similarly, in physiology, the systems underlying these model components are deeply and non-linearly coupled. Baroreflex activity (a major contributor to the low-frequency component) and respiratory sinus arrhythmia (the high-frequency component) are not independent processes that sum together; they are neurologically intertwined within the brainstem and are engaged in a constant, dynamic feedback loop[32,33]. By construction, the model is designed to assign variance to its clean, orthogonal components in a way that may not map cleanly onto the underlying biology. The estimate for the amplitude of the high-frequency component, for instance, should be interpreted as the magnitude of the best-fitting sinusoidal component in that band, assuming it can be linearly and additively separated from all other ongoing dynamics, rather than a pure measure of respiratory sinus arrhythmia. The model’s simplifying assumption of linear, additive separation is an important caveat for physiological interpretation that requires careful consideration[34].

3.2. From Statistical Constructs to Biological Reality

The intended value of this framework is its ability to distill a complex signal into a set of “interpretable” parameters. This claimed interpretability, however, occupies a position between useful abstraction and the risk of equating a parameter with a physical process[35,36]. The very act of assigning a label like “recovery rate” to a parameter ($\varphi$) cognitively biases the researcher to perceive it as a distinct, tangible biological process[37]. While the parameter perfectly describes a feature of the fitted mathematical curve, it is an unsubstantiated inferential leap to equate it with a singular neurophysiological mechanism[36,38]. The observed recovery trajectory is the net result of an immensely complex cascade of interacting neural, endocrine, and humoral processes[13,33,39]. The model’s parameterization of this trajectory is a descriptive summary, rather than an elucidation of its constituent parts.

Such interpretive discipline, separating a mathematical parameter from a singular biological process, must be applied to all model outputs. The concept of “dissonant autonomic recovery”, for example, is a statistically identifiable pattern where heart rate returns to baseline far more quickly than measures of heart rate variability[40]. Yet, its proposed link to “lingering sympathetic activation” remains a speculative hypothesis requiring external validation. Could this pattern alternatively reflect differential downregulation rates of adrenergic and cholinergic receptors, or a mismatch in the resetting of central autonomic command versus peripheral feedback? These are precisely the kinds of hypotheses the model allows one to formulate, but it possesses no inherent capacity to adjudicate between them. It provides the tool to quantify the phenomenon with precision, thereby allowing for the examination of such hypotheses[41], but it cannot, on its own, provide definitive causal evidence.

Nowhere is this interpretive caution more warranted than with the structured variance fraction, $w$. Its proposed role as an index of “autonomic coherence” is a compelling narrative, but one that rests on a fragile foundation. At its core, $w$ is simply the proportion of signal variance captured by a fixed set of sinusoids. A reduction in $w$ could plausibly reflect a pathological dysregulation of autonomic rhythms[42]. However, it could equally reflect a plethora of other confounders: an increase in movement artifact, the presence of uncorrected ectopic beats, a change in breathing frequency outside the canonical high-frequency band, or activation of other physiological processes that contribute to RRi variance but are not sinusoidal[12,13,14]. To label this and other model parameters with a term as specific is, at best, premature and risks imbuing it with a degree of physiological certainty that it does not currently warrant.

3.3. Limitations

The model’s architecture is built upon a series of rigid assumptions. The choice of a symmetric double-logistic function is primarily mathematical[20], not a physiological imperative. Many biological recovery processes exhibiting asymmetry, hysteresis, or multi-timescale recovery would be poorly approximated by this functional form. Moreover, the model is a specialist instrument designed for a single, isolated perturbation-recovery event[20]. This makes it highly effective for analyzing clean, stimulus-response experimental paradigms but makes it unsuited for the exploratory analysis of complex, naturalistic data containing multiple, overlapping events of unknown form[12,22]. It is a tool for confirmatory, not exploratory science.

Furthermore, we acknowledge that the empirical validation of this framework is presently confined to a single-subject case study. This was a deliberate choice reflecting the manuscript’s primary aim: to introduce and validate a new methodological framework, not to make a generalizable empirical claim. The N-of-1 analysis serves as a proof-of-concept, demonstrating the model’s ability to perform a deep phenotyping of an individual’s response and, most importantly, to generate precise, falsifiable hypotheses. We contend that establishing this hypothesis-generating capacity is a necessary prerequisite to the much larger, and computationally demanding, task of applying the model hierarchically to a full cohort. This foundational work is offered to empower the broader research community to test these newly generated hypotheses on larger, more diverse datasets.

The adoption of a Bayesian framework, while providing intuitive probabilistic outputs, comes with its own methodological assumptions. The posterior credible interval for a parameter is a statement of subjective belief about its value, conditional on the data, the chosen model structure, and the specified prior distributions[24]. It is not equivalent to a frequentist confidence interval and does not account for structural uncertainty, the possibility that the model itself is misspecified. This leads to what might be termed the “Rashomon effect” in modeling, where multiple, structurally different models can provide equally good fits to the data while telling different and contradictory stories about the underlying process[43]. Without a systematic exploration of this “model space”, the conclusions drawn from this single model must be regarded as highly provisional.

The model’s reliance on Hamiltonian Monte Carlo sampling makes it computationally expensive. The required runtime may not only be a practical barrier for large datasets but could also create an “analytical divide”, where such advanced methods are only accessible to research groups with substantial computational resources. This effectively sidelines the method from deployment in real-time or near-real-time clinical settings, such as intensive care monitoring, where rapid feedback is paramount. Furthermore, the model’s performance is closely tied to the fidelity of upstream data. Uncorrected premature ventricular contractions or excessive measurement noise, for instance, could drastically alter the local mean and variance of model parameters. While the model may be robust to small deviations, its behavior in the face of such gross, non-stationary artifacts is unknown and could lead to biased inference on the entire stress-recovery curve. Finally, the treatment of the residual term as unstructured “noise” is a convenient assumption. In biology, true randomness is rare; the residuals likely contain unmodeled, structured biological processes or deterministic chaos.

3.4. Future Directions

The limitations catalogued above define a necessary line for future research. This line of research should move beyond incremental model refinement and toward the systematic characterization of the model as a scientific instrument.

The most urgent task is to initiate the calibration process for this instrument. Just as a new telescope must have its optical aberrations mapped and its resolution quantified, this model requires a comprehensive characterization of its operational properties. An ideal future study would involve establishing a formal calibration protocol, creating a standardized open-source library of synthetic RRi signals with known complexities, asymmetric responses, overlapping events, various artifact types, and non-Gaussian noise, and then systematically testing the model against this library. The results, published as a detailed “specification sheet”, would delineate the model’s biases, sensitivities, and failure modes, defining its safe operating envelope for future users.

Furthermore, future work should focus on utilizing the model for quantitative hypothesis testing and falsification. The framework should be applied in experimental contexts where specific interventions are known to target certain physiological pathways. For instance, testing the effect of beta-blockade on the response amplitude and recovery rate parameters following a cold pressor or exercise-based test would provide crucial external validation, linking the model’s abstract parameters to tangible pharmacological effects. This shifts the focus from “fitting a model” to “using a model to test a theory”[41,43].

The generative framework presented here represents a conceptually powerful alternative to conventional methods for analyzing dynamic HRV. Ultimately, this model should be viewed as a specialized instrument, rather than as a perfect reflection of reality. Its value lies in its ability to translate a complex signal into a set of precise, falsifiable hypotheses. The difficult work of calibrating this instrument, understanding its distortions, and cautiously deploying it to test concrete physiological hypotheses will ultimately determine its lasting contribution to the science of autonomic control.

Methods

3.5. Model Formulation

This section details the formulation of the BAND model, a comprehensive generative framework for the RRi signal, observed at a discrete set of time points $\{t_i{\}}_{i=1}^N$. The model’s architecture is designed to deconstruct the signal into its constituent physiological components, providing a mechanistic account of the processes that shape heart rate dynamics. This is achieved by conceptualizing the signal as a probabilistic process controlled by a time-varying mean, $\mu \left(t_i\right)$, and a partitioned, time-varying standard deviation. This approach moves beyond simple curve-fitting to formalize a theory of autonomic control and its temporal evolution within a fully Bayesian framework.

A central component of the model is the partitioning of the total signal variance into two distinct, time-dependent components: a structured variance, ${\sigma }_{\mathrm{struct}}^2\left(t_i\right)$, which captures the physiologically meaningful, oscillatory patterns of HRV, and a residual variance, ${\sigma }_{\mathrm{resid}}^2\left(t_i\right)$, which accounts for unstructured, moment-to-moment fluctuations or measurement noise. The complete observation model, which integrates these components, is defined by the Normal likelihood in Equation 1.

$$\mathrm{R}\mathrm{R}\mathrm{i}\left(t_i\right)\sim \mathcal{N}\left(\mu \left(t_i\right),{\sigma }_{\mathrm{resid}}^2\left(t_i\right)\right)$$

(1)

The core of the model lies in the construction of the mean and variance components from a shared set of underlying dynamic functions. The mean trajectory, $\mu \left(t_i\right)$, is a superposition of a smoothly varying baseline trend and the synthesized structured signal itself, as shown in Equation 2.

$$\mu \left(t_i\right)=\underset{\begin{array}{c}\mathrm{Gross}\\ \mathrm{RRi}\mathrm{Trend}\end{array}}{\underbrace{\mathrm{R}{\mathrm{R}}_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}\left(t_i\right)}}+\underset{\begin{array}{c}\mathrm{Structured}\mathrm{Variability}\mathrm{Signal}\end{array}}{\underbrace{A\left(t_i\right)\cdot {\displaystyle \sum _{j=1}^J}{p}_j\left(t_i\right)\cdot S_j\left(t_i\right)}}$$

(2)

Here, $\mathrm{R}{\mathrm{R}}_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}\left(t_i\right)$ represents the gross, underlying heart period trajectory. The structured variability signal is synthesized from a set of dynamically evolving spectral oscillators, $S_j\left(t_i\right)$, which represent activity in different physiological frequency bands ($j=1,2,3$ for VLF, LF, and HF, respectively). These oscillators are weighted by time-varying proportions, $p_j\left(t_i\right)$, and their overall magnitude is controlled by a scaling amplitude, $A\left(t_i\right)$. This amplitude is deterministically calculated to ensure the signal’s variance exactly matches the target structured variance, ${\sigma }_{\mathrm{struct}}^2\left(t_i\right)$.

We partition the total signal variability, $\mathrm{S}\mathrm{D}\mathrm{N}\mathrm{N}{\left(t_i\right)}^2$, into two components. This partition is governed by a single parameter, $w\in \left[0,1\right]$, which represents the fraction of total variance that is considered physiologically structured (i.e., oscillatory). The remaining fraction, ($1-w$), is treated as unstructured residual noise. This relationship is formalized in in Equation 3.

$$\begin{array}{cc}{\sigma }_{\mathrm{struct}}^2\left(t_i\right)& =w\cdot \mathrm{S}\mathrm{D}\mathrm{N}\mathrm{N}{\left(t_i\right)}^2\hfill \\ {\sigma }_{\mathrm{resid}}^2\left(t_i\right)& =\left(1-w\right)\cdot \mathrm{S}\mathrm{D}\mathrm{N}\mathrm{N}{\left(t_i\right)}^2\end{array}$$

(3)

This formulation allows the model to simultaneously infer not only how the total amount of HRV changes over time (via $\mathrm{S}\mathrm{D}\mathrm{N}\mathrm{N}\left(t_i\right)$), but also how the nature of that variability evolves, that is the balance between predictable, oscillatory patterns and unpredictable noise (via $w$). This unified analysis facilitates a deeper understanding of autonomic regulation by capturing phenomena across both time and frequency domains within a single, coherent inferential framework.

3.5.1. Baseline Heart Period and Total Variability Trajectories

The model posits that the primary trends in both the mean heart period and its total variability are driven by a common underlying physiological response to a stimulus. To capture this, both the $\mathrm{R}{\mathrm{R}}_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}\left(t_i\right)$ and $\mathrm{S}\mathrm{D}\mathrm{N}\mathrm{N}\left(t_i\right)$ trajectories are parameterized using the same flexible double-logistic functional form.

The baseline heart period, $\mathrm{R}{\mathrm{R}}_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}\left(t_i\right)$, which quantifies the gross, underlying variations in the mean RRi, is defined in Equation 4.

$$\mathrm{R}{\mathrm{R}}_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}\left(t_i\right)=\underset{\begin{array}{c}\mathrm{Resting}\mathrm{RRi}\end{array}}{\underbrace{{\alpha }_r}}+\underset{\begin{array}{c}\mathrm{Perturbation}-\mathrm{induced}\\ \mathrm{RRi}\mathrm{Drop}\end{array}}{\underbrace{{\beta }_r\cdot {\mathcal{D}}_1\left(t_i\right)}}-\underset{\begin{array}{c}\mathrm{Post}-\mathrm{perturbation}\\ \mathrm{RRi}\mathrm{Recovery}\end{array}}{\underbrace{c_r{\beta }_r\cdot {\mathcal{D}}_2\left(t_i\right)}}$$

(4)

In this formulation, ${\alpha }_r$ represents the initial, stable heart period. The parameter ${\beta }_r$ signifies the magnitude of the decline or increase in RRi induced by the perturbation. The fractional recovery amplitude is denoted by $c_r$, where $c_r>1$ indicates an overshoot.

Concurrently, the total instantaneous variability, $\mathrm{S}\mathrm{D}\mathrm{N}\mathrm{N}\left(t_i\right)$, follows a parallel dynamic trajectory, as defined in Equation 5.

$$\mathrm{S}\mathrm{D}\mathrm{N}\mathrm{N}\left(t_i\right)=\underset{\begin{array}{c}\mathrm{Resting}\mathrm{SDNN}\end{array}}{\underbrace{{\alpha }_s}}+\underset{\begin{array}{c}\mathrm{Perturbation}-\mathrm{induced}\\ \mathrm{SDNN}\mathrm{Drop}\end{array}}{\underbrace{{\beta }_s\cdot {\mathcal{D}}_1\left(t_i\right)}}-\underset{\begin{array}{c}\mathrm{Post}-\mathrm{perturbation}\\ \mathrm{SDNN}\mathrm{Recovery}\end{array}}{\underbrace{c_s{\beta }_s\cdot {\mathcal{D}}_2\left(t_i\right)}}$$

(5)

Here, ${\alpha }_s$, ${\beta }_s$, and $c_s$ are analogous to their counterparts in the baseline model, representing the resting total SDNN, the magnitude of its suppression or increase, and its fractional recovery, respectively. The $\mathrm{S}\mathrm{D}\mathrm{N}\mathrm{N}\left(t_i\right)$ trajectory thus serves as a master controller for the total instantaneous variability in the signal. This total variance is then partitioned into a structured, oscillatory component and an unstructured, residual noise component, as defined in Equation 3.

The dynamics for both trajectories (i.e., $\mathrm{R}{\mathrm{R}}_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}$ and $\mathrm{S}\mathrm{D}\mathrm{N}\mathrm{N}$), are driven by a shared pair of logistic transition functions, ${\mathcal{D}}_1\left(t_i\right)$ and ${\mathcal{D}}_2\left(t_i\right)$, defined in Equation 6.

$$\begin{array}{cc}{\mathcal{D}}_1\left(t_i\right)& ={\left(1+e^{-\lambda \left(t_i-\tau \right)}\right)}^{-1}\hfill \\ {\mathcal{D}}_2\left(t_i\right)& ={\left(1+e^{-\varphi \left(t_i-\tau -\delta \right)}\right)}^{-1}\end{array}$$

(6)

The shared timing parameters enforce physiological coupling: $\tau$ is the inflection point of the initial response, with rate $\lambda$. The second transition, representing recovery, is offset by a delay $\delta$ and proceeds at a rate $\varphi$. This shared structure ensures that all dynamic components of the model are driven by a single, unified underlying process.

3.5.2. Generative Model for the Structured Signal

The model partitions the signal into structured components, composed of three standard physiological frequency bands: very-low-frequency (VLF) from 0.003 to 0.04 Hz, low-frequency (LF) from 0.04 to 0.15 Hz, and high-frequency (HF) from 0.15 to 0.4 Hz.

The structured component of the signal is where the model’s spectral properties are defined. Its construction involves three key elements; the dynamic spectral proportions, the latent spectral oscillators (whose amplitudes are governed by a Gaussian Process), and a deterministic amplitude inversion that ties them to the target variance.

Dynamic Frequency Band Proportions

The proportions $p_j\left(t_i\right)$ dictate how the total structured variance, ${\sigma }_{\mathrm{struct}}^2\left(t_i\right)$, is allocated across the different frequency bands at each moment. The model captures the evolution of these proportions as a smooth transition between two distinct spectral states: a baseline state (${\overrightarrow{\pi }}_{\mathrm{base}}$) and a perturbed state (${\overrightarrow{\pi }}_{\mathrm{pert}}$). The transition is determined by a single master controller function, $C\left(t_i\right)$, as shown in Equation 7. Here, $\overrightarrow{p}\left(t_i\right)$ is the 3x1 vector of proportions at time $t_i$, whose j-th element is $p_j\left(t_i\right)$

$$\overrightarrow{p}\left(t_i\right)=\left(1-C\left(t_i\right)\right)\cdot {\overrightarrow{\pi }}_{\mathrm{base}}+C\left(t_i\right)\cdot {\overrightarrow{\pi }}_{\mathrm{pert}}$$

(7)

Here, ${\overrightarrow{\pi }}_{\mathrm{base}}$ and ${\overrightarrow{\pi }}_{\mathrm{pert}}$ are simplex vectors representing the characteristic spectral distributions at rest and during peak perturbation. The master controller, $C\left(t_i\right)$, is itself built from the same logistic building blocks, ensuring the spectral transition is synchronized with the primary physiological response, as defined in Equation 8.

$$C\left(t_i\right)={\mathcal{D}}_1\left(t_i\right)\cdot \left(1-c_c\cdot {\mathcal{D}}_2\left(t_i\right)\right)$$

(8)

This function naturally transitions from 0 towards 1, with the parameter $c_c\in \left[0,1\right]$ allowing for an incomplete spectral recovery.

Latent Spectral Oscillators via a Gaussian Process Prior

A central task of the model is to infer the power of oscillations across a range of frequencies from the observed data. A simple approach might impose a rigid mathematical structure on the spectrum (e.g., a $1/f^b$ power law) or treat the power at each frequency as independent. However, physiological spectra are rarely so simple; they typically exhibit smooth, continuous shapes with broad peaks.

To capture this realistic structure, we employ a flexible, non-parametric Gaussian Process (GP) prior. A GP is a statistical tool that allows us to place a prior on an unknown function, enforcing the assumption that it is smooth[44]. In this context, it formalizes the intuition that the amplitudes of nearby frequencies should be similar, allowing the model to infer the characteristic shape of the spectrum within each band directly from the data.

Each underlying oscillator, $S_j\left(t_i\right)$, which contributes to the structured variability, is modeled as a sum of simple sine and cosine waves with pre-specified frequencies ($f_{j,k}$) but unknown amplitudes ($u_{j,k}$), as shown in Equation 9.

$$S_j\left(t_i\right)={\displaystyle \sum _{k=1}^{K_j}}\left[u_{j,k}^{\left(\mathrm{s}\mathrm{i}\mathrm{n}\right)}\mathrm{s}\mathrm{i}\mathrm{n}\left(2\pi f_{j,k}{t}_i\right)+u_{j,k}^{\left(\mathrm{c}\mathrm{o}\mathrm{s}\right)}\mathrm{c}\mathrm{o}\mathrm{s}\left(2\pi f_{j,k}{t}_i\right)\right]$$

(9)

The key innovation is the GP prior placed on the log-amplitude envelope, a vector denoted ${\overrightarrow{v}}_j$. To make the prior on the GP’s smoothness parameter stable and interpretable, the GP is defined over a log-frequency axis that has been scaled to the unit interval for each band. Let this scaled log-frequency be denoted ${\tilde{f}}_{j,k}$. The log-amplitude envelope is modeled as a single draw from a GP with a zero-mean function and a squared exponential covariance kernel (${\mathbf{K}}_j$), also known as a radial basis function (RBF) kernel, defined in Equation 10. The term in the numerator represents the squared Euclidean distance between two points on the scaled log-frequency axis. This kernel is chosen because it formalizes the assumption that the spectral envelope is a smooth, continuous function, which is a strong but reasonable prior for physiological spectra.

$${\overrightarrow{v}}_j\sim \mathcal{G}\mathcal{P}\left(\overrightarrow{0},{\mathbf{K}}_j\right), \mathrm{where}{\mathbf{K}}_j=\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{{\left({\tilde{f}}_j-{\tilde{f}}_j^{\mathrm{\text{'}}}\right)}^2}{2\hspace{0.17em}{\rho }_{gp,j}^2}}\right)$$

(10)

The covariance matrix ${\mathbf{K}}_j$ is controlled by a single hyperparameter, the length-scale ${\rho }_{gp,j}$ (a parameter controlling smoothness, not to be confused with a correlation coefficient), which determines the smoothness of the spectral envelope. The marginal standard deviation of the GP is fixed to 1. This is because the GP’s role is strictly to define the shape of the spectrum, not its absolute magnitude; the magnitude is handled by downstream normalization and scaling. The length-scale ${\rho }_{gp,j}$ is constrained to the [0,1] interval and controls the smoothness of the function on the unit-scaled input axis. A value near 1 enforces a very smooth, slowly changing shape, while a value near 0 allows for more rapid variations and narrower spectral peaks.

To ensure this model can be estimated efficiently, we use a non-centered parameterization. This is a computational technique that improves the performance of the sampler by reducing correlations between parameters[45]. Instead of estimating the highly correlated vector ${\overrightarrow{v}}_j$ directly, the model estimates a vector of independent standard normal deviates, ${\overrightarrow{z}}_{gp,j}$. The target log-amplitude envelope is then deterministically constructed via the transformation ${\overrightarrow{v}}_j={\mathbf{L}}_j{\overrightarrow{z}}_{gp,j}$, where ${\mathbf{L}}_j$ is the Cholesky factor of the covariance matrix (${\mathbf{K}}_j={\mathbf{L}}_j{\mathbf{L}}_j^{\mathrm{\top }}$). This approach makes the estimation problem much simpler for the algorithm, leading to more robust and reliable inference.

Normalization for Identifiability

A critical step in the model’s construction is to ensure that its parameters are identifiable, meaning there is a unique set of parameters that explains the data. A potential ambiguity arises because two components of the model control amplitude: the GP’s marginal standard deviation (${\alpha }_{gp,j}$) and the overall signal amplitude, $A\left(t_i\right)$. Without a constraint, the model could achieve the same result by increasing the GP’s amplitude while decreasing $A\left(t_i\right)$, making it difficult to interpret either parameter.

This normalization procedure is a critical step to ensure that the model’s parameters are identifiable. Its purpose is to assign a clear division of labor. The Gaussian Process exclusively determines the relative shape of the spectrum, while the time-varying amplitude $A\left(t_i\right)$ exclusively determines its absolute magnitude.

This is achieved by normalizing the amplitude envelope derived from the GP. First, let ${\overrightarrow{a}}_j=\mathrm{e}\mathrm{x}\mathrm{p}\left({\overrightarrow{v}}_j\right)$ be the vector of positive amplitudes determined by the GP. We then create a normalized scaling vector, ${\overrightarrow{s}}_j$, such that the expected variance of the synthesized oscillator $S_j\left(t_i\right)$ is precisely 1.

Because the final coefficients $u_{j,k}$ are generated by scaling standard normal deviates ($z_{j,k}\sim \mathcal{N}\left(0,1\right)$), the appropriate normalization must be based on the expected variance of this random process. A key property of this process is that the expected value of cross-products between independent random coefficients is zero. Consequently, the total expected variance simplifies to a sum that depends only on the diagonal elements of the pre-computed Gram matrices (${\mathbf{G}}_{\mathrm{s}\mathrm{i}\mathrm{n},j}$ and ${\mathbf{G}}_{\mathrm{c}\mathrm{o}\mathrm{s},j}$). This procedure is formalized in Equation 11.

$${\overrightarrow{s}}_j={\displaystyle \frac{{\overrightarrow{a}}_j}{\sqrt{{\displaystyle \sum _{k=1}^{K_j}}{a}_{j,k}^2\left({\left({\mathbf{G}}_{\mathrm{s}\mathrm{i}\mathrm{n},j}\right)}_{kk}+{\left({\mathbf{G}}_{\mathrm{c}\mathrm{o}\mathrm{s},j}\right)}_{kk}\right)}}}$$

(11)

The final oscillator coefficients are then constructed as $u_{j,k}^{\left(\mathrm{s}\mathrm{i}\mathrm{n}\right)}=z_{j,k}^{\left(\mathrm{s}\mathrm{i}\mathrm{n}\right)}\cdot s_{j,k}$ and $u_{j,k}^{\left(\mathrm{c}\mathrm{o}\mathrm{s}\right)}=z_{j,k}^{\left(\mathrm{c}\mathrm{o}\mathrm{s}\right)}\cdot s_{j,k}$. This procedure ensures that the GP component only informs the spectral shape, normalized to have a unit expected variance, thereby removing the ambiguity and making the model’s parameters clearly interpretable.

Deterministic Amplitude Inversion

The final step is to ensure that the synthesized structured signal, $X\left(t_i\right)=A\left(t_i\right){\displaystyle {\sum }_{j=1}^J}{p}_j\left(t_i\right)S_j\left(t_i\right)$, has a variance equal to the target structured variance, ${\sigma }_{\mathrm{struct}}^2\left(t_i\right)$, at every time point. A key feature of this model is that it does not assume the sinusoids form an orthogonal basis over the discrete, finite time domain. Instead, it computes the exact sample variance for each synthesized oscillator, $\mathrm{V}\mathrm{a}\mathrm{r}\left[S_j\left(t_i\right)\right]$, given the current values of its coefficients. This exact calculation is critical to prevent the misestimation of spectral power that arises from the non-orthogonality of sinusoidal basis functions over finite, discretely sampled time intervals.

Let ${\mathsf{\Sigma }}_S$ be the diagonal $3\times 3$ matrix whose entries are these computed variances, ${\mathsf{\Sigma }}_S\left[j,j\right]=\mathrm{V}\mathrm{a}\mathrm{r}\left[S_j\left(t_i\right)\right]$. Assuming independence between the bands, the variance of the weighted sum is given by the quadratic form in Equation 12.

$$\mathrm{V}\mathrm{a}\mathrm{r}\left[{\displaystyle \sum _{j=1}^J}{p}_j\left(t_i\right)S_j\left(t_i\right)\right]=\overrightarrow{p}{\left(t_i\right)}^{\mathrm{\top }}{\mathsf{\Sigma }}_S\hspace{0.33em}\overrightarrow{p}\left(t_i\right)$$

(12)

The variance of the complete structured signal is $\mathrm{V}\mathrm{a}\mathrm{r}\left[X\left(t_i\right)\right]=A{\left(t_i\right)}^2\cdot \left(\overrightarrow{p}{\left(t_i\right)}^{\mathrm{\top }}{\mathsf{\Sigma }}_S\hspace{0.33em}\overrightarrow{p}\left(t_i\right)\right)$. To match our target, we set this equal to ${\sigma }_{\mathrm{struct}}^2\left(t_i\right)=w\cdot \mathrm{S}\mathrm{D}\mathrm{N}\mathrm{N}{\left(t_i\right)}^2$. Solving for $A\left(t_i\right)$ yields the inversion formula in Equation 13.

$$A\left(t_i\right)={\displaystyle \frac{\sqrt{w}\cdot \mathrm{S}\mathrm{D}\mathrm{N}\mathrm{N}\left(t_i\right)}{\sqrt{\overrightarrow{p}{\left(t_i\right)}^{\mathrm{\top }}{\mathsf{\Sigma }}_S\hspace{0.33em}\overrightarrow{p}\left(t_i\right)}}}$$

(13)

This inversion is essential. It dynamically adjusts the amplitude of the synthesized spectral signal to ensure its contribution to the total variance is exactly as prescribed by the model’s high-level parameters ($\mathrm{S}\mathrm{D}\mathrm{N}\mathrm{N}\left(t_i\right)$ and $w$), thereby connecting the time-domain and frequency-domain components of the model.

3.5.3. Model Parameterization and Priors

To ensure numerical stability and efficient sampling, all model parameters are estimated on an unconstrained real-valued scale ($\mathbb{R}$). Priors are placed on these unconstrained parameters, which are then transformed back to their constrained, physically meaningful scales within the model.

Timing, Rate, and Magnitude Parameters

Parameters constrained to a specific interval are parameterized on the logit scale, while positive parameters are on the log scale. For instance, the timing parameters $\tau$ and $\delta$ are mapped to the observed time interval ($t_{\mathrm{range}}$) and the remaining time, respectively, using the inverse logit transformation, such that $\tau ={\mathrm{logit}}^{-1}\left({\tau }_{\mathrm{logit}}\right)\cdot t_{\mathrm{range}}+t_{\min }$ and $\delta ={\mathrm{logit}}^{-1}\left({\delta }_{\mathrm{logit}}\right)\cdot \left(t_{\mathrm{range}}-\tau \right)$.

The positive rate parameters $\lambda$ and $\varphi$ are simply log-transformed, i.e., $\lambda =\mathrm{e}\mathrm{x}\mathrm{p}\left({\lambda }_{\log }\right)$. The recovery coefficients for the time-domain components ($c_r,c_s$) are mapped to the interval $\left[0,2\right]$, while the recovery parameter for spectral components ($c_c$) is mapped to the $\left[0,1\right]$ interval, using a scaled logit function. Similarly, the fraction of structured variance, $w$, is constrained between 0 and 1 via $w={\mathrm{logit}}^{-1}\left(w_{\mathrm{logit}}\right)$.

The magnitude parameters (${\alpha }_r,{\beta }_r,{\alpha }_s,{\beta }_s$) are scaled relative to data-derived quantities to create dimensionless parameters whose priors are easier to specify. These transformations ensure the sampler explores a valid and well-behaved parameter space.

Spectral and Oscillator Parameters

The spectral proportions, ${\overrightarrow{\pi }}_{\mathrm{base}}$ and ${\overrightarrow{\pi }}_{\mathrm{pert}}$, are mapped from the 3-dimensional simplex to 2-dimensional real vectors using the additive log-ratio (ALR) transformation. The model estimates the unconstrained vectors, which are mapped back to the simplex via the softmax function.

The GP length-scale, ${\rho }_{gp,j}$, is constrained to the [0,1] interval to correspond with the scaled log-frequency input. This is achieved by parameterizing it on the unconstrained logit scale in the sampler. Concretely, we sample ${\rho }_{gp,j}^{\left(\mathrm{logit}\right)}\sim \mathcal{N}\left(0,1\right)$, and then deterministically transform it via the inverse logit function, ${\rho }_{gp,j}={\mathrm{logit}}^{-1}\left({\rho }_{gp,j}^{\left(\mathrm{logit}\right)}\right)$. This prior is weakly informative, allowing the data to determine the appropriate level of smoothness for the spectral envelope in each band.

Finally, the non-centered parameterization is completed by placing standard Normal priors on all latent “z” variables: the GP deviates ($z_{gp}$), the sine coefficients ($z_{\mathrm{s}\mathrm{i}\mathrm{n}}$), and the cosine coefficients ($z_{\mathrm{c}\mathrm{o}\mathrm{s}}$). This technique is critical for efficient sampling, as it decouples the hierarchical dependencies and removes pathological posterior geometries.

3.6. Simulation Studies

To evaluate the model’s performance, we conducted a series of simulation studies. The objectives were to assess the model’s ability to recover known ground-truth parameters from synthetic data and to quantitatively benchmark its reconstruction of key HRV dynamics against conventional, windowed analysis techniques. Synthetic RRi time series were generated using the model’s own generative structure with pre-specified parameter values. This process provided datasets with perfectly known underlying dynamics for $\mathrm{R}\mathrm{R}\left(t_i\right)$, $\mathrm{S}\mathrm{D}\mathrm{N}\mathrm{N}\left(t_i\right)$, and $p_j\left(t_i\right)$, serving as an objective gold standard for evaluation.

We designed three distinct scenarios to probe the model’s capabilities under diverse and challenging conditions. These scenarios: (A) a classic sympatho-vagal response, (B) an incomplete recovery with time-domain and spectral persistence, and (C) an incomplete spectral but complete time-domain recovery with high noise, were chosen specifically to test the limitations inherent in standard analysis methods. Each simulated dataset was analyzed with our complete Bayesian model and, for comparison, with conventional techniques: 1) a sliding-window analysis for time-domain metrics and 2) a Short-Time Fourier Transform (STFT) for spectral dynamics.

3.6.1. Time-Domain Trajectory Reconstruction

The capacity to accurately capture time-domain dynamics was evaluated by comparing the model’s continuous estimates of $\mathrm{R}\mathrm{R}\left(t_i\right)$ and $\mathrm{S}\mathrm{D}\mathrm{N}\mathrm{N}\left(t_i\right)$ against the stepwise approximations derived from a standard sliding-window mean and standard deviation. To provide a comprehensive evaluation of reconstruction fidelity, model performance was quantified against the ground-truth trajectories using a suite of complementary metrics. The root mean squared error (RMSE) and mean absolute error (MAE) were used to measure the average magnitude of the estimation error, with MAE offering a more robust assessment against sporadic outliers. The model’s bias (mean error) was calculated to identify any systematic tendency for over- or underestimation. Finally, the coefficient of determination ($R^2$) was used to assess the proportion of variance in the true signal captured by the estimate.

3.6.2. Spectral Dynamics Reconstruction

To evaluate the reconstruction of spectral dynamics, our model’s continuous estimates of the spectral proportions, $p_j\left(t_i\right)$, were compared against those derived from an STFT. For a direct comparison, the spectrogram generated by the STFT was normalized at each time step to yield proportional power within the VLF, LF, and HF bands. The performance of both methods was then assessed by comparing their estimated trajectories to the known ground-truth proportions using the same comprehensive suite of metrics (RMSE, MAE, Bias, and $R^2$). This comparative framework was designed to objectively quantify the trade-offs in temporal and frequency resolution inherent to windowed methods versus the continuous estimation provided by our model, especially in the presence of rapid transitions or low signal-to-noise ratios. As part of the model’s internal validation, we also confirmed that the 95% posterior credible intervals for all generative parameters consistently contained the known ground-truth values, ensuring both inferential accuracy and appropriate uncertainty quantification.

3.7. Proof of Concept

For a practical demonstration of the framework, we analyzed heart rate variability data (RRi time series) collected from one healthy individual completing a two-minute submaximal step test. This individual-level (N-of-1) analysis was designed explicitly to show how the model simplifies a complex, information-rich, and noisy time series into a few interpretable, low-dimensional parameters. Ethical approval regarding experimental protocols were obtained from the Scientific Ethics Committee of the University of Magallanes, CEC-UMAG (Nº053/SH/2023) and all methods were carried out in accordance with relevant guidelines and regulations. The subject received detailed information regarding the study objectives, procedures, and potential implications. Informed consent was obtained to ensure ethical compliance and participant autonomy.