Autonomic Dysregulation in Cardiac Arrhythmias Revealed by Reproducible Long-Term Heart Rate Variability Analysis

1. Introduction

The introduction should briefly place the study in a broad context and highlight why it is important. It should define the purpose of the work and its significance.

Heart rate variability (HRV) describes temporal fluctuations between consecutive R–R intervals of the electrocardiogram (ECG) and represents an indirect marker of the dynamic balance between the sympathetic and parasympathetic branches of the autonomic nervous system (ANS). Since the landmark Task Force report (1996) [1], HRV has been established as a noninvasive tool for evaluating cardiac autonomic modulation, physiological homeostasis, and the organism’s adaptive capacity to internal and external stimuli. Reduced HRV has been consistently associated with increased cardiovascular mortality, autonomic dysfunction, chronic stress, and elevated risk of sudden cardiac death [2,3,4].

Twenty-four-hour Holter monitoring enables continuous assessment of cardiac electrical activity under real-life conditions, facilitating the detection of transient arrhythmic episodes, silent ischemia, and rhythm disturbances that may not be evident in short-duration recordings. Unlike brief ECG acquisitions, Holter recordings capture full circadian cardiovascular dynamics, allowing HRV analysis across sleep–wake cycles and in response to physiological, emotional, and environmental variations [1,5,6].

Beyond its established clinical role, HRV has gained renewed relevance in recent years due to the availability of large open-access databases and the development of reproducible computational frameworks. Public repositories such as PhysioNet, together with robust environments like MATLAB, have facilitated standardized signal processing, transparent methodological reporting, and improved inter-study comparability in cardiovascular research [7,8,9].

Heart rate variability (HRV) describes temporal fluctuations between consecutive R–R intervals of the electrocardiogram (ECG) and represents an indirect marker of the dynamic balance between the sympathetic and parasympathetic branches of the autonomic nervous system (ANS). Since the landmark Task Force report [1], HRV has been established as a noninvasive tool for evaluating cardiac autonomic modulation, physiological homeostasis, and the organism’s adaptive capacity to internal and external stimuli. Reduced HRV has been consistently associated with increased cardiovascular mortality, autonomic dysfunction, chronic stress, and elevated risk of sudden cardiac death [2,3,4].

Twenty-four-hour Holter monitoring enables continuous assessment of cardiac electrical activity under real-life conditions, facilitating the detection of transient arrhythmic episodes, silent ischemia, and rhythm disturbances that may not be evident in short-duration recordings. Unlike brief ECG acquisitions, Holter recordings capture full circadian cardiovascular dynamics, allowing HRV analysis across sleep–wake cycles and in response to physiological, emotional, and environmental variations [1,5,6].

Over the past two decades, HRV analysis has evolved from an exploratory autonomic physiology tool into a widely used clinical biomarker in cardiology, sleep medicine, psychophysiology, and applied computational sciences [6]. This progress has been driven by the availability of large public repositories of physiological signals, particularly PhysioNet, and by robust computational environments such as MATLAB, which have promoted reproducibility and standardized analysis in cardiovascular research [10,11,12].

Numerous studies have demonstrated the clinical relevance of HRV for cardiovascular risk stratification, disease progression, and adverse event prediction. In post-myocardial infarction populations, HRV metrics derived from 24-hour Holter recordings—such as SDNN, spectral components, and nonlinear parameters—have shown independent prognostic value for total mortality and heart failure [8]. More recent investigations have confirmed the utility of HRV in emerging clinical contexts, including post-COVID-19 syndrome, type 2 diabetes, and chronic systemic diseases, revealing persistent sympathovagal alterations in long-term recordings [2,3,13].

Between 2015 and 2025, advanced methodological approaches have emerged for HRV analysis, incorporating nonlinear and multiscale metrics to better characterize cardiovascular system complexity. Techniques such as sample entropy, fractal fluctuation analysis, and multiscale entropy have demonstrated higher sensitivity for detecting subclinical autonomic alterations, even when traditional linear metrics remain unchanged [14,15,16].

In parallel, machine learning and deep learning techniques applied to long-term ECG signals have integrated HRV features as key physiological markers for predicting arrhythmias, mortality, and cardiovascular risk, reinforcing the value of prolonged Holter monitoring over isolated measurements [17,18,19]. These trends align with the emerging concept of digital biomarkers in cardiology, where HRV bridges physiology, computational analysis, and clinical decision-making [20].

Despite these advances, relevant methodological challenges persist, including lack of standardization in R-peak detection, sensitivity of certain metrics to artifacts and ectopic beats, and variability across spectral methods, which hinder inter-study comparability [8,21,22]. These limitations have motivated initiatives promoting transparency and reproducibility through open databases, code sharing, and cross-validation of analysis pipelines [7,23,24].

Within this context, the objective of the present study is to analyze heart rate variability in 24-hour Holter recordings obtained from PhysioNet using a reproducible pipeline implemented in MATLAB R2025a, integrating time-domain, frequency-domain, and nonlinear metrics to characterize autonomic modulation and arrhythmic burden in a cohort of subjects with rhythm disorders.

2. Materials and Methods

The Materials and Methods should be described with sufficient details to allow others to replicate and build on the published results.

2.1. Data Sources

Long-term ECG recordings were obtained from publicly available databases hosted on the PhysioNet platform, a reference repository for reproducible research in complex physiological signals [7,9]. Databases widely used in HRV and long-term Holter studies were selected.

Specifically, the following datasets were included: (i) the Japanese Atrial Fibrillation Holter Database (SHDB-AF), comprising 98 recordings from 93 patients with atrial fibrillation and significant arrhythmias, sampled at 125 Hz using two leads; (ii) the Long-Term ST Database, including recordings from 85 subjects (healthy and with ischemic ST-segment episodes), acquired at 250 Hz with 12-bit resolution; and (iii) the MIT-BIH Normal Sinus Rhythm Database, consisting of 18 recordings from healthy subjects without clinically significant arrhythmias, sampled at 360 Hz. Records with excessive noise or incomplete annotations were excluded prior to analysis.

Long-term ECG recordings were obtained from publicly available databases hosted on the PhysioNet platform, a reference repository for reproducible research in complex physiological signals [7,9]. Databases widely used in HRV and long-term Holter studies were selected.

Specifically, the following datasets were included: (i) the Japanese Atrial Fibrillation Holter Database (SHDB-AF), comprising 24-hour Holter recordings from patients with atrial fibrillation and significant arrhythmias, sampled at 125 Hz using two leads; (ii) the Long-Term ST Database, including 21–24 hour recordings from healthy subjects and patients with ischemic ST-segment episodes, acquired at 250 Hz with 12-bit resolution; and (iii) the MIT-BIH Normal Sinus Rhythm Database, consisting of 24-hour ECG recordings from healthy subjects without clinically significant arrhythmias, sampled at 360 Hz.

2.2. Signal Preprocessing

Signal preprocessing was performed in MATLAB R2025a following widely accepted methodological recommendations for HRV analysis in long-term recordings [1,22,23]. A fourth-order Butterworth band-pass filter (0.5–40 Hz) was applied to preserve the main ECG spectral components while removing baseline wander and high-frequency noise. Additionally, a 60 Hz notch filter was used to suppress power-line interference. R-peak detection was performed using an optimized implementation of the Pan–Tompkins algorithm. Abnormal R–R intervals deviating more than 20% from the local mean were identified as artifacts or ectopic beats and corrected using linear interpolation, a threshold commonly adopted to balance sensitivity and preservation of physiological variability. Each recording was segmented into non-overlapping one-hour windows to enable circadian HRV analysis.

Signal preprocessing was performed in MATLAB R2025a following widely accepted methodological recommendations for HRV analysis in long-term recordings [1,22,23]. A fourth-order Butterworth band-pass filter (0.5–40 Hz) and a 60 Hz notch filter were applied to remove baseline wander and power-line interference. R-peak detection was performed using an optimized implementation of the Pan–Tompkins algorithm. Abnormal R–R intervals deviating more than 20% from the local mean were corrected using linear interpolation. Each recording was segmented into non-overlapping one-hour windows to enable circadian HRV analysis.

2.3. HRV Parameter Computation

HRV parameters were computed in time, frequency, and nonlinear domains according to international standards [1,11,16].

Table 1. Equations used to compute heart rate variability (HRV).

PROCESS	EQUATIONS
R-peak (QRS) detection [25]	If $x(n)$is the digitized ECG signal: Derivative:
	$y(n)={\displaystyle \frac{1}{8T}}\left[-x(n-2)-2x(n-1)+2x(n+1)+x(n+2)\right](Eq.1)$
	Squared signal:
	$z(n)=[y(n){]}^2(Eq.2)$
	R peaks are identified when $z(n)$exceeds a dynamic threshold. RR intervals are defined as:
	$R{R}_i=t_{R(i+1)}-t_{R(i)}(Eq.3)$
Time-domain HRV analysis [11]	Let $R{R}_i$be the series of normal NN intervals, with $N$samples. Mean NN interval:
	$\overline{RR}={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N}R{R}_i(Eq.4)$
	SDNN (standard deviation of NN intervals):
	$SDNN=\sqrt{{\displaystyle \frac{1}{N-1}}{\displaystyle \sum _{i=1}^N}(R{R}_i-\overline{RR}{)}^2}(Eq.5)$
	Represents global HRV (widely used in 24-h Holter recordings). RMSSD:
	$RMSSD=\sqrt{{\displaystyle \frac{1}{N-1}}{\displaystyle \sum _{i=1}^{N-1}}(R{R}_{i+1}-R{R}_i{)}^2}(Eq.6)$
	Related to parasympathetic activity. pNN50:
	$pNN50={\displaystyle \frac{\mathrm{number}\mathrm{of}\mid R{R}_{i+1}-R{R}_i\mid >50\mathrm{ms}}{N-1}}\times 100(Eq.7)$
Frequency-domain HRV analysis [1]	Power spectral density:
	$PSD(f)=\mid \mathcal{F}\{RR(t)\}{\mid }^2(Eq.8)$
	Spectral bands: VLF: 0.0033–0.04 Hz (hormonal/thermoregulatory modulation) LF: 0.04–0.15 Hz (sympathetic + parasympathetic modulation) HF: 0.15–0.40 Hz (parasympathetic activity) Power in each band:
	$P_{band}={\int }_{f_1}^{f_2}PSD(f)df(Eq.9)$
	LF/HF ratio:
	$LF/HF(Eq.10)$
	Indicator of sympathovagal balance.
Nonlinear HRV analysis [26,27,28]	Particularly relevant for long-term ECG/Holter recordings. Poincaré plot:
	$(R{R}_i,R{R}_{i+1})(Eq.11)$
	SD1 (short-term variability):
	$SD1=\sqrt{{\displaystyle \frac{1}{2}}}\cdot RMSSD(Eq.12)$
	SD2 (long-term variability):
	$SD2=\sqrt{2SDN{N}^2-SD{1}^2}(Eq.13)$
	Approximate entropy (ApEn):
	$ApEn(m,r,N)={\mathsf{\Phi }}^m(r)-{\mathsf{\Phi }}^{m+1}(r)(Eq.14)$
	where $m$is the embedding dimension and $r$is the tolerance (typically $0.2\cdot SDNN$). Sample entropy (SampEn):
	$SampEn(m,r,N)=-\mathrm{l}\mathrm{n}({\displaystyle \frac{A}{B}})(Eq.15)$
	Lower values indicate pathological regularity. Detrended Fluctuation Analysis (DFA):
	$F(n)\propto n^{\alpha }(Eq.16)$
	where ${\alpha }_1$represents short-term and ${\alpha }_2$long-term scaling exponents.

Developed by the authors based on the cited references.

Table 1. Equations used to compute heart rate variability (HRV).

PROCESS	EQUATIONS
R-peak (QRS) detection [25]	If $x(n)$is the digitized ECG signal: Derivative:
	$y(n)={\displaystyle \frac{1}{8T}}\left[-x(n-2)-2x(n-1)+2x(n+1)+x(n+2)\right](Eq.1)$
	Squared signal:
	$z(n)=[y(n){]}^2(Eq.2)$
	R peaks are identified when $z(n)$exceeds a dynamic threshold. RR intervals are defined as:
	$R{R}_i=t_{R(i+1)}-t_{R(i)}(Eq.3)$
Time-domain HRV analysis [11]	Let $R{R}_i$be the series of normal NN intervals, with $N$samples. Mean NN interval:
	$\overline{RR}={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N}R{R}_i(Eq.4)$
	SDNN (standard deviation of NN intervals):
	$SDNN=\sqrt{{\displaystyle \frac{1}{N-1}}{\displaystyle \sum _{i=1}^N}(R{R}_i-\overline{RR}{)}^2}(Eq.5)$
	Represents global HRV (widely used in 24-h Holter recordings). RMSSD:
	$RMSSD=\sqrt{{\displaystyle \frac{1}{N-1}}{\displaystyle \sum _{i=1}^{N-1}}(R{R}_{i+1}-R{R}_i{)}^2}(Eq.6)$
	Related to parasympathetic activity. pNN50:
	$pNN50={\displaystyle \frac{\mathrm{number}\mathrm{of}\mid R{R}_{i+1}-R{R}_i\mid >50\mathrm{ms}}{N-1}}\times 100(Eq.7)$
Frequency-domain HRV analysis [1]	Power spectral density:
	$PSD(f)=\mid \mathcal{F}\{RR(t)\}{\mid }^2(Eq.8)$
	Spectral bands: VLF: 0.0033–0.04 Hz (hormonal/thermoregulatory modulation) LF: 0.04–0.15 Hz (sympathetic + parasympathetic modulation) HF: 0.15–0.40 Hz (parasympathetic activity) Power in each band:
	$P_{band}={\int }_{f_1}^{f_2}PSD(f)df(Eq.9)$
	LF/HF ratio:
	$LF/HF(Eq.10)$
	Indicator of sympathovagal balance.
Nonlinear HRV analysis [26,27,28]	Particularly relevant for long-term ECG/Holter recordings. Poincaré plot:
	$(R{R}_i,R{R}_{i+1})(Eq.11)$
	SD1 (short-term variability):
	$SD1=\sqrt{{\displaystyle \frac{1}{2}}}\cdot RMSSD(Eq.12)$
	SD2 (long-term variability):
	$SD2=\sqrt{2SDN{N}^2-SD{1}^2}(Eq.13)$
	Approximate entropy (ApEn):
	$ApEn(m,r,N)={\mathsf{\Phi }}^m(r)-{\mathsf{\Phi }}^{m+1}(r)(Eq.14)$
	where $m$is the embedding dimension and $r$is the tolerance (typically $0.2\cdot SDNN$). Sample entropy (SampEn):
	$SampEn(m,r,N)=-\mathrm{l}\mathrm{n}({\displaystyle \frac{A}{B}})(Eq.15)$
	Lower values indicate pathological regularity. Detrended Fluctuation Analysis (DFA):
	$F(n)\propto n^{\alpha }(Eq.16)$
	where ${\alpha }_1$represents short-term and ${\alpha }_2$long-term scaling exponents.

Developed by the authors based on the cited references.

Time-domain metrics included mean NN, SDNN, RMSSD, and pNN50. Frequency-domain analysis employed Welch’s method to estimate spectral power in the LF and HF bands and the LF/HF ratio. HRV metrics were computed for each one-hour segment and subsequently averaged across the full 24-hour recording to obtain subject-level descriptors. Nonlinear analysis included Poincare plot descriptors (SD1 and SD2) and sample entropy. All computations were implemented using custom MATLAB scripts and validated with PhysioNet and Signal Processing Toolbox functions.

HRV parameters were computed in the time, frequency, and nonlinear domains according to international standards [1,11,16]. Time-domain metrics included mean NN, SDNN, RMSSD, and pNN50. Frequency-domain analysis employed Welch’s method to estimate LF, HF, and LF/HF ratio. Nonlinear analysis included Poincare plot descriptors (SD1, SD2) and sample entropies. All computations were implemented using custom MATLAB scripts and validated with PhysioNet and Signal Processing Toolbox functions.

3. Results

Analysis of 24-hour Holter recordings from 85 subjects with arrhythmia revealed reduced global HRV and substantial inter-subject variability. Mean SDNN was 0.165 ± 0.057 s, RMSSD was 0.187 ± 0.082 s, and pNN50 reached 52.2 ± 27.8%, reflecting heterogeneous parasympathetic modulation across subjects. Spectral analysis showed low absolute power in both LF (0.0075 ± 0.0148 s²) and HF (0.0096 ± 0.0084 s²) bands, with an average LF/HF ratio of 0.78 ± 0.99. Automated arrhythmia detection identified a high arrhythmic burden, with mean counts of 2996 ± 2773 tachycardia episodes, 333 ± 723 bradycardia episodes, and 1923 ± 1060 premature beats per subject, highlighting pronounced rhythm instability.

Table 2. Summary of HRV analysis results in MATLAB and mean ± standard deviation.

File	SDNN	RMSSD	pNN50	LF	HF	LFHF	Tachycardia episodes	Bradycardia episodes	Premature beats
00m.mat	0,128459456	0,063537573	17,50612288	0,000272532	0,000575891	0,473235382	1081	6	1099
01m.mat	0,282687272	0,390257347	82,69946259	0,004966746	0,033526281	0,148144841	3479	960	3998
03m.mat	0,268105594	0,175670813	30,96622541	0,005131433	0,007379025	0,695408078	1375	1429	1741
74m.mat	0,220109852	0,227641068	77,60969746	0,010686637	0,015663641	0,682257531	2952	464	2697
75m.mat	0,141914335	0,198288907	78,81355932	0,003758839	0,008210399	0,457814352	7491	23	3055
Mean values	0,165029097	0,186757171	52,21541698	0,007474031	0,00959514	0,78390094	2995,988095	333,1190476	1923,309524
Standard deviation	0,056790095	0,081669569	27,82885138	0,014794079	0,008438391	0,996784183	2772,85534	723,002398	1060,198674

Prepared by the authors.

Table 2. Summary of HRV analysis results in MATLAB and mean ± standard deviation.

File	SDNN	RMSSD	pNN50	LF	HF	LFHF	Tachycardia episodes	Bradycardia episodes	Premature beats
00m.mat	0,128459456	0,063537573	17,50612288	0,000272532	0,000575891	0,473235382	1081	6	1099
01m.mat	0,282687272	0,390257347	82,69946259	0,004966746	0,033526281	0,148144841	3479	960	3998
03m.mat	0,268105594	0,175670813	30,96622541	0,005131433	0,007379025	0,695408078	1375	1429	1741
74m.mat	0,220109852	0,227641068	77,60969746	0,010686637	0,015663641	0,682257531	2952	464	2697
75m.mat	0,141914335	0,198288907	78,81355932	0,003758839	0,008210399	0,457814352	7491	23	3055
Mean values	0,165029097	0,186757171	52,21541698	0,007474031	0,00959514	0,78390094	2995,988095	333,1190476	1923,309524
Standard deviation	0,056790095	0,081669569	27,82885138	0,014794079	0,008438391	0,996784183	2772,85534	723,002398	1060,198674

Prepared by the authors.

Analysis of 24-hour Holter recordings from 85 subjects with arrhythmia revealed reduced global HRV, low spectral power in both LF and HF bands, and a high arrhythmic burden characterized by frequent tachycardia, bradycardia, and premature beats. The wide dispersion of HRV indices reflected marked inter-subject heterogeneity, consistent with diverse arrhythmic patterns.

4. Discussion

The present analysis of HRV derived from 24-hour Holter recordings demonstrates that arrhythmias are associated with marked autonomic dysregulation, reflected by reduced global variability, low spectral power, and altered sympathovagal balance. These findings are consistent with previous PhysioNet-based studies employing reproducible computational pipelines, which reported similar reductions in SDNN, RMSSD, and HF power in populations with high arrhythmic burden [9,29].

The pronounced inter-subject heterogeneity observed across HRV metrics underscores the physiological diversity inherent to arrhythmic populations and highlights the importance of integrating nonlinear metrics and standardized preprocessing. The use of a fully reproducible MATLAB-based pipeline and open-access data represents a methodological strength, facilitating replication, benchmarking, and comparison across future studies.

The findings demonstrate that arrhythmias are associated with significant autonomic dysregulation, reflected by reduced HRV and altered sympathovagal balance. The combination of linear and nonlinear metrics provided complementary insight into autonomic impairment. The reproducible MATLAB-based pipeline and use of open PhysioNet data represent key methodological strengths, facilitating replication and comparison across studies.

5. Conclusions

Heart rate variability analysis from long-term Holter recordings enables robust characterization of autonomic modulation in patients with arrhythmias, revealing reduced variability, altered sympathovagal balance, and high arrhythmic burden. The substantial heterogeneity observed across subjects emphasizes the need for standardized and reproducible analytical approaches.

The proposed computational pipeline, implemented in MATLAB and based on open PhysioNet data, ensures transparency, replicability, and methodological consistency. These features support potential clinical translation, including integration into decision-support systems, cardiovascular risk stratification frameworks, and the development of digital biomarkers for autonomic dysfunction.

Heart rate variability analysis from long-term Holter recordings provides robust characterization of autonomic modulation in patients with arrhythmias. The proposed reproducible computational pipeline enables standardized HRV assessment and supports future applications in cardiovascular risk stratification and digital biomarker development.