Heart Rate Fragmentation as a Novel Biomarker of Adverse Cardiovascular Events: The Multi-Ethnic Study of Atherosclerosis

Background: A major objective of precision medicine is the elucidation of non-invasive biomarkers of cardiovascular (CV) risk. Recently, we introduced a new dynamical marker of sino-atrial instability, termed heart rate fragmentation (HRF), which outperformed traditional and nonlinear heart rate variability metrics in separating ostensibly healthy subjects from patients with coronary artery disease. Accordingly, we hypothesized that HRF may be a dynamical biomarker of adverse cardiovascular events (CVEs).

Methods: This study employed data from a cohort of participants in the Multi-Ethnic Study of Atherosclerosis (MESA), a prospective study of sub-clinical heart disease. Interbeat interval time series (n = 1963), derived from the electrocardiographic channel of the polysomnogram study, were analyzed using the newly introduced metrics of fragmentation, as well as traditional heart rate variability (HRV) indices and the short-term detrended fluctuation analysis exponent. Cox regression analysis was used to assess the association between HR dynamic indices and CV outcomes in unadjusted and adjusted models.

Results: The mean (± SD) follow-up time was 2.97 ± 0.63 years. In adjusted models, higher fragmentation was significantly associated with incident CVEs (number of events; hazard ratio [95% confidence interval]: n = 72, 1.43 [1.16–1.76]) and CV death (n = 21; 1.65 [1.15–2.36]). The traditional HRV and the fractal indices were not associated with CVEs or CV death. The most discriminatory fragmentation indices added significant value to Framingham and MESA CV risk indices in all analyses.

Conclusion: Our findings show that HRF has promise as a non-invasive, automatable biomarker of CV risk. The basic mechanisms underlying fragmentation remain to be delineated. Its association with incident outcomes raises the possibility of connections to degenerative changes in the multisystem network controlling SAN function.

1. Introduction

This study describes a novel noninvasive biomarker of cardiovascular (CV) risk based on heart rate dynamics. In healthy adults at rest and during sleep, the highest frequency at which the sino-atrial node (SAN) rate fluctuates varies between ~0.15 and 0.40 Hz (Figures 1A₁–A₄). These oscillations, referred to as respiratory sinus arrhythmia, are due to vagally-mediated coupling between the SAN and breathing. However, not all fluctuations in heart rate (HR) at or above the respiratory frequency are attributable to vagal tone modulation. Under pathologic conditions, an increased density of reversals in HR acceleration sign, not consistent with short-term parasympathetic control, can be observed (Figures 1B₁–B₄). This dynamical biomarker of electrophysiologic instability has recently been identified and termed heart rate fragmentation (HRF) (Costa et al., 2017a). A set of metrics (computational probes) for its quantification was also introduced (Costa et al., 2017a,b).

FIGURE 1

Figure 1. Heart rate dynamics from two MESA participants, A with respiratory sinus arrhythmia (black tracings) and B with fragmented sinus rhythm (red tracings). Normal-to-normal (NN) sinus interval time series for the entire sleep period (A₁,B₁) and for a 70-second window (A₂,B₂). Twelve-second ECG recordings (A₃,B₃). Poincare plots for the entire sleep period (A₄,B₄) and the 70-second NN time series shown in A₂ and B₂ (C). The green circles highlight the “inflection points,” where the changes in heart rate acceleration sign occur. Histogram of the percentage of inflection points (PIP) calculated with a moving window of 1000 NN intervals (D). Neither participant had prevalent or incident CVEs. However, they were in two different risk categories: the Framingham CVD risk index was 2.4 (1st percentile of MESA participants) for A and 15.6 (55th percentile) for B. The time series from participant B was 30% more fragmented (average PIP = 65%) than the one from participant A (average PIP = 50%). Traditional HRV and DFA α₁ indices were comparable: mean NN interval, 906 and 957 ms; rMSSD, 25.3 and 27.6 ms; pNN50, 5.0 and 6.5%; HF power, 348 and 310 ms²; DFA α₁, 1.05 and 1.18, for participants A and B, respectively.

Perhaps the most explicit example of HRF is the subtle supraventricular arrhythmia termed sinus node alternans (Binkley et al., 1995), in which the time between consecutive sinus beats oscillates between two values, short (S) and long (L) following an SLSL pattern. However, HRF includes not only pure (2:1) sinus node alternans but also quasi-periodic and more irregular variants of normal-to normal (NN) alternation. As Figures 1A₃,B₃ illustrates, clinical recognition of such patterns is difficult from standard electrocardiograms (ECGs). Traditional heart rate variability (HRV) indices and other metrics such as those derived from Poincaré plots (Figures 1A₄,B₄,C) may be also of limited use in the identification of HRF patterns. The basic mechanisms of fragmentation, involving either anomalous sinus beats (Lewis, 1920) or supraventricular ones originating near the SAN, are unresolved (Costa et al., 2017a,b).

The potential importance of HRF is several-fold. First, it produces a high degree of short-term variability that may be mistaken as a marker of healthy vagal control when standard measures of short-term HRV are used. Second, its presence supports the delineation of a new class of biomarkers of cardiac risk. The latter is premised on the conjectured link between HRF and the breakdown in one or more components of the control system (and/or in their interactions) regulating SAN function. Notably, earlier reports of what were first termed “sinus extrasystoles” as well as sinus alternans (Lewis, 1920; Binkley et al., 1995) were from patients who were older or had organic heart disease. Third, the investigation of fragmentation may yield new insights into SAN functionality in health, aging and disease.

In recent studies (Costa et al., 2017a,b) we analyzed annotated Holter recordings [University of Rochester Telemetric Holter ECG Warehouse (THEW)] from healthy subjects and patients with advanced coronary artery disease (CAD) using the newly devised HRF metrics. Fragmentation was found to significantly increase as a function of the participants' age in both the healthy population and the one with CAD. In contrast, most short-term HRV indices did not significantly change with the participants' age in the CAD group. Furthermore, fragmentation was higher in patients with CAD than in healthy subjects, during both estimated awake and sleep periods, while traditional HRV metrics did not discriminate the two groups.

The general motivation for the present study was to assess the potential utility of the novel indices of HRF as predictors of adverse cardiovascular events (CVEs) and CV mortality, using the large Multi-Ethnic Study of Atherosclerosis (MESA). This ongoing prospective cohort study (Bild et al., 2002) was designed to investigate the prevalence, correlates and progression of subclinical cardiovascular disease (CVD) in a multi-ethnic population free of overt clinical CVD at study entry. We specifically hypothesized that HRF would: (1) be positively associated with cross-sectional age; (2) be positively associated with incident CVEs and CV death; and (3) outperform traditional HR dynamical measures. In addition, we also sought to determine if fragmentation metrics added value to prediction tools computed from clinical measures, namely the Framingham Heart Study (D'Agostino et al., 2008) and MESA CV risk indices (McClelland et al., 2015).

2. Methods

2.1. Study Population and Data Collection

The MESA study has been described in detail previously (Bild et al., 2002). Briefly, over a period of approximately 2-years, starting in July 2000, 6,814 persons between the ages of 45 and 84 years of age without evident clinical CVD were recruited at six field centers in the US. Institutional review boards from each study site approved the conduct of this study, and written informed consent was obtained from all participants.

A sleep ancillary study was conducted in conjunction with MESA's fifth examination (2010–2013). The study enrolled 2,060 participants who underwent unattended, in-home polysomnography (PSG) following a standardized protocol (Redline et al., 1998). The data obtained using the 15-channel Compumedics Somte System (Compumedics LTd., Abbottsville, Australia) were scored at the Brigham and Women's Hospital centralized reading center by trained technicians using published guidelines (Redline et al., 2007). The apnea-hypopnea index (AHI) was calculated based on the average number per hour of sleep of all apneas plus hypopneas associated with ≥3% oxygen desaturation or arousal.

The ECG channels, sampled at 256 Hz, were processed using Compumedics Somte software for detection and classification of the QRS complexes (R-points) as normal sinus, supraventricular premature or ventricular premature complexes. The automated annotations were reviewed by a trained technician, who made appropriate corrections. Both the NN and the R-to-R (RR) interval time series were analyzed in the present study.

Participants with one or more of the following were excluded: poor signal quality (n = 35), pacemaker (n = 13), atrial fibrillation (AF) at the time of the PSG (n = 22), <2 h of combined sleep periods scored as rapid eye movement (REM), stage 1, 2, 3, or 4 (n = 16), and <75% normal sinus beats between sleep onset and termination (n = 11). Participants with CVEs before the PSG (n = 185) and seven others for whom the last recorded follow-up was prior to the PSG were excluded from the analyses of the associations between HR dynamical metrics and incident CVEs. These participants were included in analyses of CV mortality.

2.2. Clinical Follow-Up and Event Classification

In addition to clinical exams, participants are followed every 9–12 months to inquire about hospital admissions, CV outpatient diagnoses and procedures, and deaths. Discharge diagnosis codes are obtained for all hospitalizations and medical records are obtained when heart failure, myocardial infarction, stroke, or death are reported. For those over age 65 and enrolled in fee-for-service Medicare, claims data are also used to identify diagnosis and procedure codes. Trained personnel abstract any hospital records suggesting possible CVEs, which are then adjudicated by physicians at the coordinating center. Nonfatal endpoints in MESA include congestive heart failure, angina, myocardial infarction, percutaneous coronary intervention, coronary bypass grafting or other revascularization procedure, resuscitated cardiac arrest, peripheral arterial disease, stroke (non-hemorrhagic) and transient ischemia attack (TIA). Cardiovascular deaths, as adjudicated by committee review, included fatalities directly related to stroke or coronary heart disease. For other deaths, the underlying cause are obtained through state or vital statistics departments. The definition and adjudication of these events have been described in detail previously (Bild et al., 2002; Bluemke et al., 2008; Yeboah et al., 2012). The cut-off date for the surveillance period was December 31, 2014.

2.3. Fragmentation Analysis

Fragmentation analysis was performed for 1963 subjects using both NN and RR interval time series. Fragmentation analysis is described in detail in Costa et al. (2017a,b). Briefly, original interbeat interval time series, {s_i}, 1 ≤ i ≤ L (L, time series length) were mapped to a ternary symbolic sequence as follows: “−1” if ΔNN_i < −4 ms, “0” if −4 ≤ ΔNN_i ≤ 4 ms, and “1” if ΔNN_i > 4 ms (Costa et al., 2017b). Note that, since the ECG signals were sampled at 256 Hz, the resolution of the interbeat interval time series is 1/256 ~ 4 ms. Therefore, only NN (or RR) intervals whose difference was >4 ms or <−4 ms were considered different from each other.

Transitions from HR acceleration to HR deceleration (“−1” to “1”) or vice-versa (“1” to “−1”), and from HR acceleration or HR deceleration to no change in HR (“−1” to “0,” “1” to “0,”) or vice-versa (“0” to “−1,” “0” to “1”) were termed “inflection points.” The percentage of inflection points (PIP) (Costa et al., 2017a) constitutes a measure of HRF reflecting its overall degree of prevalence.

To assess the prevalence of dynamical patterns with increasing degrees of fragmentation, the percentages of sequences of 4 consecutive symbols, w_i = (s_i, s_i+1, …, s_i+l−1), 1 ≤ i ≤ L−l+1, termed “words,” with 0, 1, 2, and 3 inflection points were calculated. We refer to these word classes as W₀, W₁, W₂, and W₃, respectively. The full lexicon that comprises 81 different words is given in Costa et al. (2017b). Words derived from the NN (RR) interval time series were termed NN (RR) words. The words in groups W_0 and W_1 are the least fragmented (most “fluent”), those in groups W_2 and W_3 are the most fragmented.

All fragmentation metrics were calculated for the entire sleep period. In the case of PIP, standard sleep stage periods (awake, rapid eye movement [REM], stages N1, N2, and N3), in addition to periods of awake before sleep onset or after sleep termination were also analyzed. Individuals were included if cumulatively they had at least 2000 NN intervals during a given sleep stage. Thus, the number of participants in each of these sub-analyses was a fraction of the number of participants in the analyses of the entire sleep period. Specifically, the number of participants (the number of those with incident CVEs) for the different periods was: awake, 1379 (61); REM, 1451 (54); N1, 1056 (51); N2, 1695 (70); N3, 791 (23); and awake before sleep onset or after sleep termination, 1449 (58).

2.4. Traditional HRV Analysis

The following traditional time domain HRV indices (HRV, 1996) were calculated for 1963 subjects using NN interval time series: (1) the average of all NN intervals (AVNN), (2) mean of the standard deviations (SDs) of NN intervals in all 5-min segments (SDNNIDX), (3) the square root of the mean of the squares of differences between adjacent NN intervals (rMSSD), and (4) the percentage of differences between adjacent NN intervals that are greater than 50 ms (pNN50). The following traditional frequency domain HRV indices were calculated: (1) the total spectral power of all NN intervals between 0.15 and 0.4 Hz (HF) and (2) the ratio of low to high frequency power (LF/HF). Each of these metrics was calculated using a 5-min sliding window (without overlap), with more than 150 beats and more than 75% NN intervals, between sleep onset and sleep termination. Power spectrum estimates were obtained using the Lomb periodogram method, which does not require missing data points (due to removal of ectopic, misdetections, and artifact) in a time series to be interpolated (Moody, 1992, 1993). A total of 170,527 windows were analyzed. For each subject, the values from the different windows were averaged. The source code used for these computations is published on The Research Resource for Complex Signals, (PhysioNet) website (Goldberger et al., 2000; Mietus and Goldberger, 2008).

2.5. Detrended Fluctuation Analysis: Short-Term Scaling Index

Detrended fluctuation analysis (Peng et al., 1995) was developed to quantify the correlation properties of a time series. The methods is based on the assessment of the slope of the regression line of the logarithm of F(n) vs. the logarithm of n, where F(n) is the root-mean-square fluctuation of the integrated and detrended data, computed using windows of length n. For the analysis of heart rate time series, two indices, α₁ and α₂, quantifying short and long-term behavior, respectively, have been proposed. We focused on α₁, defined as the slope of F(n) vs. n, for 4 ≤ n ≤ 11 (Pikkujamsa et al., 1999). Time series for which α~0.5 are uncorrelated (random). Those with α > 0.5 and those with α < 0.5 are long-range correlated and anti-correlated, respectively.

Discontinuities in the NN interval time series due to premature beats, misdetections and artifact were dealt with in the following way. If the gap was < 3 s (typically due to the removal of an ectopic beat), an interpolated beat was inserted. If the gap was wider (typically due to ECG artifact that prevents accurate detection of the peak of the QRS complexes), then the segments that preceded and followed the gap, were “stitched” together. The α₁ exponent was calculated for non-overlapping segments of 1,000 intervals and the results were averaged. The source code used in these computations is published on the PhysioNet website (Goldberger et al., 2000; Mietus et al., 2001).

2.6. Statistical Analysis

Continuous variables are summarized as median, first and third quartiles, unless otherwise indicated. Categorical variables are presented as numbers and percentages.

The associations between independent variables and both incident CVEs and CV mortality were assessed using Cox proportional hazard analysis. Efron's method (Efron, 1977) was used to handle ties. Failure time in the individuals with incident CVEs was the time between the PSG study and the time of event diagnosis. For participants without CVEs, the failure time was the time between the PSG study and the latest follow-up, death, or loss to follow-up. Statistical significance was set at a p-value < 0.05. The independent variables were: the fragmentation indices, PIP, W₀, W₁, W₂, and W₃, derived from both NN and RR interval time series; the HR dynamical indices: AVNN, SDNNIDX, rMSSD, pNN50, HF, and HF/LF; and the short-term fractal index DFA α₁.

Both unadjusted (Model 1) and adjusted models were considered. Adjustments included: (i) traditional CV risk factors: age, sex, systolic blood pressure, total cholesterol, high density lipoprotein (HDL) cholesterol, current smoking status, hypertension medication, diabetes and lipid lowering medication (Model 2); (ii) the Framingham Heart Study 10-year risk index (D'Agostino et al., 2008) (Model 3), and (iii) the MESA CV risk index without coronary calcification (Model 4).

Standardized hazard ratios (per one-SD increase in the independent variable) were calculated with associated 95% confidence intervals (CI). The assumption of proportional hazards was tested using a global test based on Schoenfeld residuals (Grambsch and Therneau, 1994). No violations were noted. The predictive power of the survival models was assessed using Harrell's C statistic. The likelihood ratio test was used to compare the fit of two nested models (null model vs. null model + HR dynamical metric). The three null models were models 2, 3, and 4 described above. The null hypothesis for each of the likelihood ratio tests was that the two nested models fitted the data equally well. Rejection of the null hypothesis implied that the larger model fitted the data better, indicating that a given HR dynamical metric added value to the null model.

Linear regression models with a quadratic term were used to describe the possible nonlinear (U-shaped) relationship between: i) HRF and short-term HRV indices (example of a model with PIP and rMSSD, $\text{PIP}={\beta }_1*\text{ln}(\text{rMSSD})+{\beta }_2*{\left[\text{ln}(\text{rMSSD})\right]}^2+\alpha$, where α is a constant); and ii) short-term HRV indices and the participants' age.

In all analyses, the variables W₀, rMSSD, pNN50, SDNN, or HF were logarithmically transformed since the models with these log-transformed variables fitted the data better than those with untransformed ones. In all other cases, the analyses were performed using untransformed variables.

3. Results

Over a median follow-up period of 1,080 [916–1,259] (median [Q₁–Q₃]) days after the PSG study, 72 out of the 1,771 participants without prevalent CVEs suffered their first adverse event: myocardial infarction (n = 16), resuscitated cardiac arrest (n = 1), angina (n = 14), percutaneous coronary intervention (n = 21), coronary bypass graft (n = 3), other revascularization (n = 6), congestive heart failure (n = 10), peripheral vascular disease (n = 8), transient ischemic attack (n = 5), CV death (n = 14) or stroke (non-hemorrhagic brain infarction, n = 17). From a total of 1,963 participants (1,771 without and 192 with prevalent CVEs), 21 died of CVD.

Characteristics of MESA participants without and with a CVE during follow-up are summarized in Table 1. Individuals who developed CVEs were older and more likely to be male and have diabetes. They tended to have higher seated HR and higher systolic blood pressure. In addition, this risk group tended to have lower sleep efficiency and a higher apnea-hypopnea index. The differences between those who died and did not were qualitatively similar to the differences between those with and without incident CVEs.

TABLE 1

Table 1. Characteristics of MESA participants without and with a CVE during follow-up.

3.1. Relationship of HRF, Traditional HRV and Nonlinear Indices With the Participants' Age

All fragmentation indices, derived from either NN or RR interval time series, were significantly associated with the participants' age. The Pearson correlation coefficients (ρ) for PIP, W₀, W₁, W₂, and W₃ were 0.35, −0.17, −0.31, 0.10, and 0.35, respectively. PIP (Figure 2) and the percentages of fragmented words W₂ and W₃ increased with the participants' age at the estimated rates of 0.28, 0.09, and 0.35% per year of age, respectively. The percentages of fluent words W₀ and W₁ decreased with the participants' age at the rates of −0.06 and −0.39% per year, respectively. Slightly higher rates of change were observed for the indices derived from RR interval time series (not shown).

FIGURE 2

Figure 2. Tukey boxplots of ln(rMSSD) and PIP for participants in successive age groups. PIP, percentage of inflection points; rMSSD, root mean square of the successive differences.

The short-term traditional HRV indices, rMSSD, pNN50 and HF, did not vary linearly with the participants' age. Instead, the association between age and short-term variability depended on the participants' age itself. Figure 2 illustrates, using the representative example of rMSSD, how the amount of short-term variability varied across different age groups. Variability was higher in both the lowest (<55 years) and highest (≥85 years) age groups compared to intermediate ones (U-shape relationship).

The short-term fractal index, α₁ showed a small but significant decrease with cross-sectional age. The Spearman correlation coefficient was −0.18 (p = 0.000).

3.2. Unadjusted Analyses of Risk of Incident CVEs

All fragmentation indices, calculated from both the NN and the RR time series, were significantly associated with incident events (Table 2: Model 1). The association was positive for fragmented words, W₂ and W₃, as well as PIP, and negative for fluent (less fragmented) words, W₀ and W₁. The most discriminatory of the fragmentation indices was the word W₁ derived from the RR interval time series. A one-SD increase in W₁ was associated with a 44% (95% CI: 28–66%) decrease in the rate of CVEs. This variable performed comparably to the Framingham Heart Study and MESA CV risk indices (Figure 3, top panels).

TABLE 2

Table 2. Association of fragmentation, traditional HRV and short-term fractal indices with incident CVEs in unadjusted and adjusted models for standard risk factors.

FIGURE 3

Figure 3. Kaplan-Meier survival curves of analyses of incident CVEs (top panels) and CV mortality (bottom panels), showing the percentage of symbolic words with one inflection point (W₁) derived from RR interval time series (left panels), the Framingham (middle panels) and MESA (right panels) CV risk indices. Q₁–Q₄ indicate first to fourth quartiles. Note that participants in the highest quartile of the Framingham and MESA risk indices had the worst prognosis and those in the highest quartile of word class W₁ (more fluent, less fragmented) had the best prognosis.

None of the traditional time (AVNN, SDNNIDX, rMSSD, pNN50) and frequency domain (HF, LF/HF) HRV variables was significantly associated with incident events. DFA α₁ was also not associated with incident events.

3.3. Analyses of Risk of Incident CVEs Adjusted for Traditional Risk Factors

The models with each of the fragmentation and traditional HRV variables were adjusted for standard CV risk factors: age, sex, systolic blood pressure, total cholesterol, HDL cholesterol, current smoking status, hypertension medication, diabetes and lipid lowering medication. In these analyses (Table 2: Model 2), all fragmentation indices remained significantly associated with incident CVEs. For example, a one-SD increase in PIP was associated with a 43% (15–78%) increase in the rate of incident CVEs. In addition, all fragmentation indices, with the exception of W₃ calculated from the RR interval time series, added significant value to the model with only the risk factors. In contrast, none of the traditional time (AVNN, SDNNIDX, rMSSD, pNN50) and frequency domain (HF, LF/HF) HRV variables was significantly associated with incident events. DFA α₁ was also not associated with incident events.

Age and sex were the only covariates significantly associated with incident CVEs in the adjusted models that included a fragmentation index. Specifically, for a model with PIP, the standardized hazard ratios for each of the covariates were: i) age (SD = 9.1 years), 1.33 (95% CI: 1.03–1.73), p = 0.029; ii) sex, 1.95 (1.15–3.30), p = 0.013; iii) total cholesterol (SD = 36.6 mg/dl), 0.91 (0.69–0.21), p = 0.533; iv) HDL (SD = 16.3 mg/dl), 0.87 (0.65–1.14), p = 0.311; v) diabetes, 1.39 (0.78–2.46), p = 0.260; vi) smoking status, 1.66 (0.74–3.71), p = 0.215; vii) hypertension medication, 1.02 (0.61–1.71), p = 0.938; viii) systolic blood pressure (SD = 20.3 mmHg), 1.12 (0.89–1.42), = 0.342; and ix) lipid lowering medication, 0.97 (0.75–1.25), p = 0.786.

The results did not qualitatively change after further adjusting the analyses for each of the following variables, singly or in combination: race, body mass index, waist circumference, diastolic blood pressure, seated heart rate, use of hypoglycemic agents, total sleep time, sleep efficiency, the apnea-hypopnea index, measures of coronary artery calcification (McClelland et al., 2015) (total Agatston calcium score) and carotid plaque (Gepner et al., 2015, 2017).

For all sleep stages, including awake during sleep, and awake before sleep onset and/or sleep termination, the degree of HRF, measured by PIP, was positively associated with risk of incident CVEs in fully adjusted models. The levels of significance varied slightly per stage. Specifically, the standardized hazard ratios and 95% CI were: REM [1.34 (1.01–1.77)], stage N1 [1.51 (1.13–2.00)], stage N2 [1.36 (1.08–1.70)], and stage N3 [1.55 (0.99–2.43)], awake [1.68 (1.27–2.22)] and awake before sleep onset and/or sleep termination [1.71 (1.34–2.19)].

3.4. Analyses of Risk of Incident CVEs Adjusted for the Framingham and MESA CV Risk Indices

In general, the fragmentation indices remained significantly associated with the risk of CVEs in models adjusted for the Framingham and the MESA CV risk indices (Table 3). Specifically, increased fragmentation, that is, higher PIP, lower percentages of fluent words W₀ and W₁, and higher percentages of fragmented words W₂ and W₃, were significantly associated with increased risk of events. Neither the traditional HRV measures nor DFA α₁ showed significant associations with incident CVEs.

TABLE 3

Table 3. Association of fragmentation, traditional HRV and short-term fractal indices with incident CVEs in models adjusted for the Framingham (model 3) and MESA (model 4) risk indices.

The risk indices in each of these models were also significantly associated with incident CVEs. Specifically, a one-SD increase in the Framingham and in the MESA risk indices, was associated with 80% (95% CI: 43–125%), and 55% (33–81%) increase in the rate of adverse CVEs, respectively. Harrell's C statistic was 0.666 and 0.678 for the Framingham and MESA risk indices, respectively. Overall the best model, with a Harrell's C statistic of 0.703, was the one that combined the word group W₁ derived from RR intervals, with the MESA risk index.

Models incorporating the Framingham risk index and any of the fragmentation measures, except W₂ derived from the NN interval time series, performed significantly better than the Framingham index itself. Similarly, all models that included any of the fragmentation measures in addition to the MESA risk index performed significantly better than the MESA index itself.

The traditional HRV variables and α₁ were not significantly associated with risk of incident CVEs either in unadjusted or adjusted models. Adding one of these indices to a model with a fragmentation index did not improve its performance.

3.5. Analyses of Risk of CV Death: Unadjusted and Adjusted for the Framingham and MESA Risk Indices

Higher PIP and lower percentage of W₁ words were significantly associated with increased risk of CV death in unadjusted analyses as well as in analyses adjusted for the Framingham and the MESA CV risk indices (Table 4, Figure 3). Specifically, a one-SD (~7%) increase in PIP derived from the analysis of NN interval time series was associated with an increase in the rate of CV death of 89% (95% CI: 34–168%) in unadjusted models and of 65% (15–136%) and 67% (19–136%) in models adjusted for the Framingham and the MESA risk indices, respectively. A one-SD (~11%) increase in the percentage of W₁ (“fluent” or least fragmented) words, also derived from the analysis of NN interval time series, was associated with a decrease in the rate of CV death of 59% (37–75%) in unadjusted models and of 52% (21–71%) and 55% (25–72%) in models adjusted for Framingham and the MESA risk indices, respectively. Similar results were obtained from the analyses of the RR interval time series. Further adjusting the models by prevalent CVD did not change the significance of the associations between the fragmentation metrics and risk of CV death.

TABLE 4

Table 4. Association of fragmentation, traditional HRV and short-term fractal indices with CV death in models adjusted for the Framingham and MESA CV risk indices.

Lower percentages of fluent words W₀ were associated with increased risk of CV death in all models. However, these associations were weaker than those with word W₁. The percentages of word W₂ and W₃, the most fragmented, were positively associated with the risk of CV death in all models. However, the significance of the associations depended on the particular model (Table 4). A one-SD increase in the Framingham (~9%) and in the MESA risk (~6%) indices was associated with 137% (95% CI: 48–278%) and 62% (35–96%) increase in the rate of CV death, respectively. Harrell's C statistic for the former and latter models was 0.749 and 0.797, respectively. Both variables, PIP and W₁, calculated from NN and RR interval time series, added significant information to models with either of these risk indices. The results for word groups W₀, W₂, and W₃ depended on the particular model. Overall, the best model, with a Harrell's C statistic of 0.838, was the one that combined the word W₁, derived from RR intervals, with the MESA risk index.

The traditional HRV variables and α₁ were not significantly associated with risk of CV death either in unadjusted or adjusted models. Adding one of these indices to a model with a fragmentation index did not improve its performance.

3.6. Relationship of HRF With Short-Term Traditional HRV Indices

Nonlinear (U-shape) relationships were found between fragmentation indices and traditional HRV measures of short-term variability. Figure 4 shows one representative example, the relationship between PIP and ln(rMSSD). For the first three quartiles of rMSSD values, (i.e., for rMSSD values below the 75th percentile of rMSSD, specifically, ln(rMSSD) < 3.7 ms), the degree of fragmentation and the amount of short-term variability were inversely correlated. In the upper quartile of rMSSD values, however, the degree of fragmentation and the amount of short-term variability were positively associated. Qualitatively similar results were found for pNN50 and HF power.

FIGURE 4

Figure 4. Scatter plot of PIP vs. the natural logarithm of rMSSD. The fitting line is described by the equation: PIP = −26.2*ln(rMSSD)+3.54*[ln(rMSSD)]²+105.3. The 95% CI for the 1st, 2nd and 3rd terms are: (−30.0 − −22.5), (3.03–4.05) and (98.5–112.1), respectively. PIP, percentage of inflection points; rMSSD, root mean square of the successive differences.

4. Discussion

The present investigation was designed to test the association of quantitative measures of HRF, a newly defined property of short-term sino-atrial rhythm dynamics, with adverse CV outcomes in MESA, a large ongoing multicenter study of individuals recruited from the general community. The key findings of this study were that: (1) increased HRF was significantly associated with risk of incident CVEs and CV mortality; (2) measures of fragmentation added value to Framingham and MESA risk prediction indices; and (3) traditional metrics of short-term HRV as well as a nonlinear index (DFA α₁) were not associated with incident CVEs or CV death.

The development of the concept of fragmentation and its quantitative metrics was motivated in part by an apparent paradox in the results of traditional time and frequency domain analyses of a number of studies (de Bruyne et al., 1999; Stein et al., 2005; Huikuri and Stein, 2012; Drawz et al., 2013; Costa et al., 2017a; Raman et al., 2017). While the mechanisms of the relatively slow (i.e., below the respiratory frequency) variations in HR are attributable to complex interactions between the parasympathetic and the sympathetic branches of the autonomic nervous system, faster variations, in the range of 0.15−0.40Hz, are mainly attributed to vagal tone modulation (HRV, 1996). Therefore, short-term (high frequency) measures of HR dynamics, such as rMSSD, pNN50, and HF power, are typically interpreted as surrogate measures of cardiac vagal tone. In contexts where cardiac vagal tone modulation is known to be diminished, for example, with advanced aging and established CVD, these “vagal” measures are expected to be lower. In fact, a monotonic decrease in high frequency variability with increasing age has been reported in multiple cross-sectional studies of ostensibly healthy adults (Pikkujamsa et al., 1999; Bonnemeier et al., 2003; Costa et al., 2017a). Furthermore, the association between extremely low variability and adverse outcomes is well-documented (HRV, 1996; Huikuri and Stein, 2012; O'Neal et al., 2016).

However, in certain cases (de Bruyne et al., 1999; Stein et al., 2008; Almeida-Santos et al., 2016; Raman et al., 2017; Wdowczyk et al., 2018) a paradoxical increase in short-term HR fluctuations was observed in contexts where reduced vagal tone would have been expected based on age and/or advanced heart disease. In the present study, we observed a U-shaped relationship between traditional short-term HRV measures and cross-sectional age (Figure 2). From approximately ages 45–65 years the amount of short-term variability decreased. Subsequently, variability increased despite the well-known decrease in cardiac vagal tone modulation with advancing age. These results provide further evidence that in cohorts of middle-aged to elderly individuals, such as MESA, traditional HRV indices may fail to reflect accurately changes in cardiac vagal tone.

We introduced fragmentation analysis to quantify short-term HR dynamics in such types of cohorts. The term “fragmented heart rate” refers to rhythms in which HR acceleration sign changes at a frequency higher than that attributable to healthy vagal tone modulation of the SAN. These anomalous rhythms include but are not limited to classic sinus alternans and its variants. Due to their anti-correlated structure, i.e., the fact that an increase in HR is likely followed by a decrease and vice-versa, fragmented sinus rhythms are not random/erratic. In fact, they exhibit a much higher degree of predictability/regularity (Figure 1, bottom panels) than random rhythms. The term “fragmented” was also chosen to help convey the putative pathophysiologic concept of regulatory network disintegration and/or breakdown of neuroautonomic coupling between heart rate and respiration.

An intuitive measure of fragmentation is the percentage of changes in HR acceleration sign, that is, PIP, in NN (or RR) time series. Most recently, a symbolic dynamical approach was introduced (Costa et al., 2017b) that quantifies the frequency of occurrence of different patterns of fluctuations, from least fragmented (most fluent) to most fragmented. These fragmentation indices were originally tested in studies of publicly available databases from the Rochester THEW archives. Of note, if the amplitude of the fluctuations is low (e.g., ≲80ms), fragmentation is unlikely to be detected in clinical readings of short (typically 10 s) and long (Holter) ECG recordings. However, as shown in Figure 1, visual detection of HRF is facilitated by inspection of HR time series.

The origins of HRF remain speculative. Possible pathophysiologic mechanisms include increased automaticity in or proximal to the SAN, exit block in the SAN area, modulated sinus/atrial parasystole, abrupt pacemaker shifts in the SAN (Boyett et al., 2000; Kodama et al., 2004), beat-to-beat changes due to perturbations in atrial stretch receptors (Costa et al., 2017a,b) or alterations in membrane and cellular pacemaker clocks (Lakatta et al., 2010). These electrophysiologic perturbations, in turn, may be related to underlying atrial (Sosnowski and Petelenz, 1995; Roberts-Thomson et al., 2007; Goette et al., 2017) or ventricular disease. Systemic and local factors that may also contribute to pathophysiologic dysregulation of SAN dynamics include inflammation, degeneration, fibrosis and calcification (Costa et al., 2017a). Future experimental and mathematical modeling studies will hopefully shed light on the putative links between these and other mechanisms and fragmentation. Possible genomic associations with HRF remain to be explored.

Based on the analyses of the THEW databases (Costa et al., 2017a,b), we hypothesized that increased HRF might be a biomarker of increased risk of incident CVEs and CV death. To explore these hypotheses, we analyzed HR dynamics from a subset of the participants in the MESA. This national study is one of the largest prospective investigations designed to track meticulously the course of CVD in an ethnically diverse population free of overt clinical CVD at study entry. Two types of ECG recordings with detailed follow-up data were available at the time of this analysis: traditional 10-s ECGs and the ECG channel of the PSG studies. We chose to examine the latter given the non-stationary nature of HR dynamics, [which implies that statistical time series analysis tools are most reliable when applied to “long” recordings (HRV, 1996)], and the fact that previous studies (Costa et al., 2017a,b) have shown that the discriminatory power of HRF was generally comparable during awake and sleep periods. In MESA, we confirmed and extended this initial observation. Specifically, we found that the positive association between HRF and incident CVEs (in fully adjusted analyses) could be detected for each of the sleep stages and for awake periods before sleep onset and/or after sleep termination.

A concomitant finding was the absence of associations between the most commonly used traditional short-term HRV measures and incident CVEs and CV death, in unadjusted and adjusted models. These results are not as surprising as they might appear at first glance. First, the U-shaped relationship between traditional short-term HRV measures and the participants' age (Figure 2) was indicative that such measures would be of limited utility in this cohort. Second, as previously mentioned, traditional measures of HRV, in contrast to fragmentation measures, also failed to discriminate patients with CAD from ostensibly healthy subjects in databases provided by the University of Rochester (Costa et al., 2017a,b). Third, HRF, by increasing variability not ascribable to physiologic vagal tone modulation may confound the results of traditional HRV. The nonlinear (U-shaped) relationship between HRF and short-term variability (Figure 4) supports this conjecture. In fact, the subgroups of participants with the lowest and highest amounts of variability, which, in the conventional HRV framework, would be presumed to have the highest and lowest risk of adverse events, respectively, both showed increased HRF. These findings are consistent with the reports of Stein et al. (Stein, 2002; Stein et al., 2008) and others (de Bruyne et al., 1999; Huikuri and Stein, 2012; Drawz et al., 2013).

Of note, fragmentation and traditional HRV indices differ in the following major way. By construction, fragmentation indices do not mathematically depend on mean HR and/or the amplitude of its fluctuations. These salient attributes derive from the fact that accelerations/decelerations are defined as increments/decrements in HR of any magnitude. In contrast, by definition, short-term HRV indices quantify information that is encoded in the amplitude of the fluctuations. As previously mentioned rMSSD, pNN50, and HF power were not associated with risk of incident CVEs and CV mortality. The other widely used HRV metrics, AVNN, SDNNIDX, and LF/HF were also not associated with risk of incident CVEs and CV mortality. Furthermore, none of the traditional indices improved the performance of models that included a fragmentation index.

In this study the short-term DFA index, α₁, decreased with the participants' age. In addition, α₁ was lower in those with adverse events (incident CVEs and CV death) than in those without. However, the associations were not statistically significant in any of the models.

Of potential basic and translational importance is the fact that we used both the NN and RR time series in fragmentation analyses. The NN series were employed using expert edited time series from THEW and MESA to insure that the fragmentation was likely related to beats originating in or near to SAN, therefore not distinguishable from sinus beats, at least from the single lead provided. The RR time series were used to demonstrate that fragmentation analysis, not relying on detailed beat annotation, had comparable (or even superior) discriminatory power to that employing NN time series, substantially facilitating the development of automatable analyses.

Finally, it is worth emphasizing that the Framingham and the MESA indices are composite measures incorporating information related to demographics (age, sex, race), lifestyle (smoking status), vital signs (blood pressure) and blood analytes (lipids, glucose). In contrast, each HRF index is a dynamical measure reflecting the frequency of the changes in HR acceleration sign. How can such a single metric based on a continuous ECG keep “pace” with these other multivariable risk stratification tools? The answer may relate in part to the fact that HRF indices are dynamical measures, not static probes. In contrast, blood pressure, cholesterol, glucose and other common biomarkers are single time point readouts (“snapshot”). Thus, they provide limited information on the dynamics of the underlying control mechanisms.

More generally, HRF metrics belong to the new class of dynamical probes (Cysarz et al., 2000; Goldberger, 2006; Costa et al., 2014; Porta et al., 2015; Makowiec et al., 2015; Hoyer et al., 2017; Zhang et al., 2017) that can report, in real-time, on salient aspects of integrative, multiscale, regulatory systems and of their breakdown with aging and disease. The use of these probes may enhance the clinical utility of traditional risk assessment tools (Tables 3, 4) and of other emerging technologies, such as genomic profiling. In furtherance of the goals of precision medicine, the dynamical property of HRF may also constitute a novel “target” for therapeutic interventions.

5. Conclusion

HRF, a newly defined manifestation of anomalous short-term sino-atrial variability, is associated with increased risk of cardiac adverse events and cardiac mortality in MESA. The measures, readily computable from long, continuous ECGs, added value to the canonical Framingham and MESA CV risk indices. Furthermore, fragmentation measures outperformed conventional HRV measures and the short-term DFA α₁ index. Future studies are needed to confirm the utility of HRF measures in risk stratification, and in prediction of cardiomyopathies and atrial fibrillation.

Author Contributions

The idea for this study was developed by MC and AG. Data acquisition was the responsibility of SR. The analysis of the data was performed by MC. RD and SH directed the statistical analysis. All co-authors participated in the interpretation of the findings. MC and AG wrote the manuscript in collaboration with the other co-authors.

Funding

The MESA study was supported by contracts HHSN268201500003I, N01-HC-95159, N01-HC-95160, N01-HC-95161, N01-HC-95162, N01-HC-95163, N01-HC-95164, N01-HC-95165, N01-HC-95166, N01-HC-95167, N01-HC-95168, N01-HC-95169, and R35HL135818 from the National Heart, Lung, and Blood Institute. The MESA Sleep ancillary study was supported through R01HL098433. The authors gratefully acknowledge support from G. Harold and Leila Y. Mathers Charitable Foundation, the James S. McDonnell Foundation, and the NIH grants R01GM104987, R24HL114473, and R01HL127659.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors thank the other investigators, the staff, and the participants of the MESA study for their valuable contributions. A full list of participating MESA investigators and institutions can be found at http://www.mesa-nhlbi.org.

Abbreviations

AVNN, average value of the NN intervals; BMI, body mass index; BP, blood pressure; CAD, coronary artery disease; CV, cardiovascular; CVD, cardiovascular disease; CVE, cardiovascular event; DFA, detrended fluctuation analysis; ECG, electrocardiogram; HDL, high density lipoprotein; HF, total spectral power of all NN intervals between 0.15 and 0.4 Hz; HR, heart rate; HRF, heart rate fragmentation; HRV, heart rate variability; LF/HF, ratio of low to high frequency power; MESA, Multi-Ethnic Study of Atherosclerosis; NN, normal-to-normal sinus (interval); PIP, percentage of inflection points (changes in heart acceleration sign); pNN50, percentage of differences between adjacent NN intervals that are greater than 50 ms; PSG, polysomnography or polysomnogram; rMSSD, square root of the mean of the squares of differences between adjacent NN intervals; RR, R-to-R (interval); SAN, sino-atrial node; SD, standard deviation; SDNNIDX, mean of the standard deviations of NN intervals in all 5-min segments; W_j, segment, termed “word”, of 4 consecutive differences between adjacent interbeat intervals presenting j changes in heart rate acceleration sign.

References

Almeida-Santos, M. A., Barreto-Filho, J. A., Oliveira, J. L., Reis, F. P., da Cunha Oliveira, C. C., and Sousa, A. C. (2016). Aging, heart rate variability and patterns of autonomic regulation of the heart. Arch. Gerontol. Geriatr. 63, 1–8. doi: 10.1016/j.archger.2015.11.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Bild, D. E., Bluemke, D. A., Burke, G. L., Detrano, R., Diez Roux, A. V., Folsom, A. R., et al. (2002). Multi-Ethnic Study of Atherosclerosis: objectives and design. Am. J. Epidemiol. 156, 871–881. doi: 10.1093/aje/kwf113

PubMed Abstract | CrossRef Full Text | Google Scholar

Binkley, P. F., Eaton, G. M., Nunziata, E., Khot, U., and Cody, R. J. (1995). Heart rate alternans. Ann. Intern. Med. 122, 115–117. doi: 10.7326/0003-4819-122-2-199501150-00007

PubMed Abstract | CrossRef Full Text | Google Scholar

Bluemke, D. A., Kronmal, R. A., Lima, J. A., Liu, K., Olson, J., Burke, G. L., et al. (2008). The relationship of left ventricular mass and geometry to incident cardiovascular events: the MESA (Multi-Ethnic Study of Atherosclerosis) study. J. Am. Coll. Cardiol. 52, 2148–2155. doi: 10.1016/j.jacc.2008.09.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Bonnemeier, H., Richardt, G., Potratz, J., Wiegand, U. K., Brandes, A., Kluge, N., et al. (2003). Circadian profile of cardiac autonomic nervous modulation in healthy subjects: differing effects of aging and gender on heart rate variability. J. Cardiovasc. Electrophysiol. 14, 791–799. doi: 10.1046/j.1540-8167.2003.03078.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Boyett, M. R., Honjo, H., and Kodama, I. (2000). The sinoatrial node, a heterogeneous pacemaker structure. Cardiovasc. Res. 47, 658–687. doi: 10.1016/S0008-6363(00)00135-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Costa, M. D., Davis, R. B., and Goldberger, A. L. (2017a). Heart rate fragmentation: a new approach to the analysis of cardiac interbeat interval dynamics. Front. Physiol. 8:255. doi: 10.3389/fphys.2017.00255

PubMed Abstract | CrossRef Full Text | Google Scholar

Costa, M. D., Davis, R. B., and Goldberger, A. L. (2017b). Heart rate fragmentation: a symbolic dynamical approach. Front. Physiol. 8:827. doi: 10.3389/fphys.2017.00827

PubMed Abstract | CrossRef Full Text | Google Scholar

Costa, M. D., Henriques, T., Munshi, M. N., Segal, A. R., and Goldberger, A. L. (2014). Dynamical glucometry: use of multiscale entropy analysis in diabetes. Chaos 24:033139. doi: 10.1063/1.4894537

PubMed Abstract | CrossRef Full Text | Google Scholar

Cysarz, D., Bettermann, H., and van Leeuwen, P. (2000). Entropies of short binary sequences in heart period dynamics. Am. J. Physiol. Heart Circ. Physiol. 278, H2163–H2172. doi: 10.1152/ajpheart.2000.278.6.H2163

PubMed Abstract | CrossRef Full Text | Google Scholar

D'Agostino, R. B., Vasan, R. S., Pencina, M. J., Wolf, P. A., Cobain, M., Massaro, J. M., et al. (2008). General cardiovascular risk profile for use in primary care: the Framingham Heart Study. Circulation 117, 743–753. doi: 10.1161/CIRCULATIONAHA.107.699579

PubMed Abstract | CrossRef Full Text | Google Scholar

de Bruyne, M. C., Kors, J. A., Hoes, A. W., Klootwijk, P., Dekker, J. M., Hofman, A., et al. (1999). Both decreased and increased heart rate variability on the standard 10-second electrocardiogram predict cardiac mortality in the elderly: the Rotterdam Study. Am. J. Epidemiol. 150, 1282–1288. doi: 10.1093/oxfordjournals.aje.a009959

PubMed Abstract | CrossRef Full Text | Google Scholar

Drawz, P. E., Babineau, D. C., Brecklin, C., He, J., Kallem, R. R., Soliman, E. Z., et al. (2013). Heart rate variability is a predictor of mortality in chronic kidney disease: a report from the CRIC Study. Am. J. Nephrol. 38, 517–528. doi: 10.1159/000357200

PubMed Abstract | CrossRef Full Text | Google Scholar

Efron, B. (1977). The efficiency of Cox's likelihood function for censored data. J. Am. Stat. Assoc. 72, 557–565. doi: 10.1080/01621459.1977.10480613

CrossRef Full Text | Google Scholar

Gepner, A. D., Young, R., Delaney, J. A., Budoff, M. J., Polak, J. F., Blaha, M. J., et al. (2017). Comparison of carotid plaque score and coronary artery calcium score for predicting cardiovascular disease events: the Multi-Ethnic Study of Atherosclerosis. J. Am. Heart Assoc. 6:e005179. doi: 10.1161/JAHA.116.005179

PubMed Abstract | CrossRef Full Text | Google Scholar

Gepner, A. D., Young, R., Delaney, J. A., Tattersall, M. C., Blaha, M. J., Post, W. S., et al. (2015). Comparison of coronary artery calcium presence, carotid plaque presence, and carotid intima-media thickness for cardiovascular disease prediction in the Multi-Ethnic Study of Atherosclerosis. Circ. Cardiovasc. Imaging 8:e002262. doi: 10.1161/CIRCIMAGING.114.002262

PubMed Abstract | CrossRef Full Text | Google Scholar

Goette, A., Kalman, J. M., Aguinaga, L., Akar, J., Cabrera, J. A., Chen, S. A., et al. (2017). EHRA/HRS/APHRS/SOLAECE expert consensus on atrial cardiomyopathies: definition, characterization, and clinical implication. Heart Rhythm 14, e3–e40. doi: 10.1016/j.hrthm.2016.05.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldberger, A. L. (2006). Giles F. Filley lecture. Complex systems. Proc. Am. Thorac. Soc. 3, 467–471. doi: 10.1513/pats.200603-028MS

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldberger, A. L., Amaral, L. A., Glass, L., Hausdorff, J. M., Ivanov, P. C., Mark, R. G., et al. (2000). PhysioBank, PhysioToolkit, and PhysioNet: components of a new research resource for complex physiologic signals. Circulation 101, E215–E220. doi: 10.1161/01.CIR.101.23.e215

PubMed Abstract | CrossRef Full Text | Google Scholar

Grambsch, P. M., and Therneau, T. M. (1994). Proportional hazards tests and diagnostics based on weighted residuals. Biometrika 81, 515–526. doi: 10.1093/biomet/81.3.515

CrossRef Full Text | Google Scholar

Hoyer, D., Zebrowski, J., Cysarz, D., Goncalves, H., Pytlik, A., Amorim-Costa, C., et al. (2017). Monitoring fetal maturation-objectives, techniques and indices of autonomic function. Physiol. Meas. 38, R61–R88. doi: 10.1088/1361-6579/aa5fca

PubMed Abstract | CrossRef Full Text | Google Scholar

HRV (1996). Heart rate variability: standards of measurement, physiological interpretation and clinical use. Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology. Circulation 93, 1043–1065.

Huikuri, H. V., and Stein, P. K. (2012). Clinical application of heart rate variability after acute myocardial infarction. Front. Physiol. 3:41. doi: 10.3389/fphys.2012.00041

PubMed Abstract | CrossRef Full Text | Google Scholar

Kodama, I., Honjo, H., Dobrzynski, H., and Boyett, M. R. (2004). “Cellular mechanisms of sinoatrial activity,” in Cardiac Electrophysiology: From Cell to Bedside, 4th Edn., eds D. P. J. Zipes (Philadelphia, PA: Saunders), 192–202.

Google Scholar

Lakatta, E. G., Maltsev, V. A., and Vinogradova, T. M. (2010). A coupled SYSTEM of intracellular Ca2+ clocks and surface membrane voltage clocks controls the timekeeping mechanism of the heart's pacemaker. Circ. Res. 106, 659–673. doi: 10.1161/CIRCRESAHA.109.206078

PubMed Abstract | CrossRef Full Text | Google Scholar

Lewis, T. (1920). The Mechanism and Graphic Registration of the Heart Beat. New York, NY: PB Hoeber.

Makowiec, D., Wejer, D., Kaczkowska, A., Zarczynska-Buchowiecka, M., and Struzik, Z. R. (2015). Chronographic imprint of age-induced alterations in heart rate dynamical organization. Front. Physiol. 6:201. doi: 10.3389/fphys.2015.00201

PubMed Abstract | CrossRef Full Text | Google Scholar

McClelland, R. L., Jorgensen, N. W., Budoff, M., Blaha, M. J., Post, W. S., Kronmal, R. A., et al. (2015). 10-year coronary heart disease risk prediction using coronary artery calcium and traditional risk factors: derivation in the MESA (Multi-Ethnic Study of Atherosclerosis) with validation in the HNR (Heinz Nixdorf Recall) Study and the DHS (Dallas Heart Study). J. Am. Coll. Cardiol. 66, 1643–1653. doi: 10.1016/j.jacc.2015.08.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Mietus, J. E., and Goldberger, A. L. (2008). Heart Rate Variability Analysis With the HRV Toolkit. Available online at: https://physionet.org/tutorials/hrv-toolkit/

Mietus, J. E., Peng, C.-K., and Moody, G. B. (2001). Detrended Fluctuation Analysis. Available online at: https://www.physionet.org/physiotools/dfa/dfa.c

Moody, G. B. (1992). Lomb. Available online at: https://www.physionet.org/physiotools/lomb

Moody, G. B. (1993). “Spectral analysis of heart rate without resampling,” in Proceedings of Computers in Cardiology Conference (Long Beach, CA: IEEE Computer Society), 715–718.

Google Scholar

O'Neal, W. T., Chen, L. Y., Nazarian, S., and Soliman, E. Z. (2016). Reference ranges for short-term heart rate variability measures in individuals free of cardiovascular disease: the Multi-Ethnic Study of Atherosclerosis (MESA). J. Electrocardiol. 49, 686–690. doi: 10.1016/j.jelectrocard.2016.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, C.-K., Havlin, S., Stanley, H. E., and Goldberger, A. L. (1995). Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. Chaos 5, 82–87. doi: 10.1063/1.166141

PubMed Abstract | CrossRef Full Text | Google Scholar

Pikkujamsa, S. M., Makikallio, T. H., Sourander, L. B., Raiha, I. J., Puukka, P., Skytta, J., et al. (1999). Cardiac interbeat interval dynamics from childhood to senescence. Circulation 100, 393–399. doi: 10.1161/01.CIR.100.4.393

PubMed Abstract | CrossRef Full Text | Google Scholar

Porta, A., Baumert, M., Cysarz, D., and Wessel, N. (2015). Enhancing dynamical signatures of complex systems through symbolic computation. Philos. Trans. A Math. Phys. Eng. Sci. 373:20140099. doi: 10.1098/rsta.2014.0099

PubMed Abstract | CrossRef Full Text | Google Scholar

Raman, D., Kaffashi, F., Lui, L. Y., Sauer, W. H., Redline, S., Stone, P., et al. (2017). Polysomnographic heart rate variability indices and atrial ectopy associated with incident atrial fibrillation risk in older community-dwelling men. JACC. Clin. Electrophysiol. 3, 451–460. doi: 10.1016/j.jacep.2016.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Redline, S., Budhiraja, R., Kapur, V., Marcus, C. L., Mateika, J. H., Mehra, R., et al. (2007). The scoring of respiratory events in sleep: reliability and validity. J. Clin. Sleep. Med. 3, 169–200.

PubMed Abstract | Google Scholar

Redline, S., Sanders, M. H., Lind, B. K., Quan, S. F., Gottlieb, D. J., Bonekat, W. H., et al. (1998). Methods for obtaining and analyzing unattended polysomnography data for a multicenter study. SHHS Research Group. Sleep 21, 759–767.

Google Scholar

Roberts-Thomson, K. C., Sanders, P., and Kalman, J. M. (2007). Sinus node disease: an idiopathic right atrial myopathy. Trends Cardiovasc. Med. 17, 211–214. doi: 10.1016/j.tcm.2007.06.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Sosnowski, M., and Petelenz, T. (1995). Heart rate variability. Is it influenced by disturbed sinoatrial node function? J. Electrocardiol. 28, 245–251. doi: 10.1016/S0022-0736(05)80263-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Stein, P. K. (2002). Heart rate variability is confounded by the presence of erratic sinus rhythm. Comput. Cardiol. 26, 669–672. doi: 10.1109/CIC.2002.1166861

CrossRef Full Text | Google Scholar

Stein, P. K., Domitrovich, P. P., Hui, N., Rautaharju, P., and Gottdiener, J. (2005). Sometimes higher heart rate variability is not better heart rate variability: results of graphical and nonlinear analyses. J. Cardiovasc. Electrophysiol. 16, 954–959. doi: 10.1111/j.1540-8167.2005.40788.x

CrossRef Full Text | Google Scholar

Stein, P. K., Le, Q., and Domitrovich, P. P. (2008). Development of more erratic heart rate patterns is associated with mortality post-myocardial infarction. J. Electrocardiol. 41, 110–115. doi: 10.1016/j.jelectrocard.2007.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Wdowczyk, J., Makowiec, D., Gruchala, M., Wejer, D., and Struzik, Z. R. (2018). Dynamical landscape of heart rhythm in long-term heart transplant recipients: a way to discern erratic rhythms. Front. Physiol. 9:274. doi: 10.3389/fphys.2018.00274

PubMed Abstract | CrossRef Full Text | Google Scholar

Yeboah, J., McClelland, R. L., Polonsky, T. S., Burke, G. L., Sibley, C. T., O'Leary, D., et al. (2012). Comparison of novel risk markers for improvement in cardiovascular risk assessment in intermediate-risk individuals. JAMA 308, 788–795. doi: 10.1001/jama.2012.9624

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X. D., Pechter, D., Yang, L., Ping, X., Yao, Z., Zhang, R., et al. (2017). Decreased complexity of glucose dynamics preceding the onset of diabetes in mice and rats. PLoS ONE 12:e0182810. doi: 10.1371/journal.pone.0182810

PubMed Abstract | CrossRef Full Text | Google Scholar