The control mechanisms of heart rate dynamics in a new heart rate nonlinear time series model

Figure 1. Cardiac sympathetic/parasympathetic nervous system (SNS/PNS) and forces acting on heart rate (HR). (A) The relative strength of PNS against SNS activity (PNS-SNS activity or PNS-SNS) played two roles: a response to cardiac stress (noise or stimulus) and a restoring force to maintain HR around equilibrium. (B) There were three forces acting on heart rhythms: cardiac stress as a disturbance that may cause heart rate to deviate from the equilibrium heart rate, a restoring force to reduce the derivation via PNS-SNS, and a myocardial electrical resistive force against forward action potential propagation that hinders rapid changes in heart rate.

Figure 2. Heart rate (HR) slope and amplitude indicators η1 and η2. (A) HR for chf01, f1y02, and s20641. η1 = 0 (chf01, f1y02) corresponded to a horizontal mean line, a slightly larger absolute slope indicator corresponded to a slightly larger slope of the mean heartbeat line (s20641) (an upward-trending mean line for η1 = 0.0013: 131–300; a downward-trending mean line for η1 = −0.0005: 301–700; and an upward-trending mean line for η1 = 0.0012: 701–900). A slightly larger amplitude indicator corresponded to a larger amplitude (η2 = 0.0451 for chf01, η2 = 0.0630 for f1y02; η2 = 0.1722 for s20641:1–900 points). (B) Scatter plots showing a strong positive correlation between the sample standard deviation from the mean sd and estimated η2 for the YOUNG and ELDER groups (no cardiac diseases), CHF (congestive heart failure), ST (ST segment alteration), SCD (sudden cardiac death), and VT/VF (ventricular tachycardia/ventricular fibrillation) groups.

Figure 3. Prediction performance of the noise-driven HRV responder. The HRV responder within the stable fixed-point range is called the homeostatic HRV responder. The predicted traces by the noise-driven homeostatic HRV responder (B) mimicked real HRV traces for 0209.vt1 (ventricular tachycardia), f1y07 and f1y06 (young healthy), and chf06 and chf01 (congestive heart failure). (A) Moreover, the amplitude of HRV decreased as the amplitude indicator η2 decreased.

Figure 4. Predictions of HR dynamics and long memory. (A) The out-of-sample predicted HRV trace exhibited similar trace to the actual HRV from an elderly subject (f1o04: 1801–3600 points). (B) The out-of-sample predicted ACF trace better mimicked the actual ACF trace based on the heartbeat data from the same subject (f1o04: 1801–3600 points) with a restoration delay of 10. The predicted ACF values were calculated using the data generated by the noise-driven homeostatic HRV responder as a realization of the random variable. The predicted and actual HRV responders in (A,B) had the same parameter values estimated using the heartbeat data from the same subject (f1o04: 1–1800 points). A process has long memory if its autocorrelation function (ACF) decays more slowly than an exponential decay (e.g., at a hyperbolic rate or an oscillation) as the lag increases. These results show the ability of the noise-driven homeostatic HRV responder to predict the HRV dynamics and oscillatory long memory for a long even restoration delay.

Figure 5. Fitting performance of the stimulus-driven homeostatic HRV responder. (A) The real HRV trace plots of 0213.vf1 (ventricular fibrillation), 0067.vt1 (ventricular tachycardia), 0261.vt1 (ventricular tachycardia), 0209.vf1 (ventricular fibrillation), 0040.vt2 (ventricular tachycardia), and 0217.vf1 (ventricular fibrillation). (B) The fitted traces produced by the stimulus-driven homeostatic HRV responder. These fitted traces mimic real HRV, implying that the stimulus-driven homeostatic HRV responder can describe pathological arrhythmias when cardiovascular perturbations are unpredictably pathological stimuli or intense vibration. These arrhythmias reflect the pathological state of cardiovascular and autonomic nervous systems.

Figure 6. Sensitivity of HR parameters. (A,B) The sensitivity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S(\omega ,\phi )$$\end{document}S(ω,ϕ) of the HR parameter ϕ = α, β, γ, σ, η2 to changes in the mean of unexpected perturbation Δωi (i = 1, ⋯ , 20) based on the HR parameter estimates of heart rate time series f1y02 and s20641. (C,D) The sensitivity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S(\sigma ,\phi )$$\end{document}S(σ,ϕ) of ω, α, β, γ, η1 to changes in the standard deviation of unexpected perturbation Δσi for f1y02 and s20641. Heart rate time series f1y02 (900 samples) followed the NLARI process within the stable fixed-point range, while heart rate time series s20641 (900 samples) followed the NLARI process on the boundary of the stable fixed-point range.

Figure 7. Time scales covered by the HR model. The polylines on the left- and right-hand side reflect changes in HR parameter estimates on average at the second scales m = 1, ⋯ , 10 (s) and the minute scales m = 1, ⋯ , 10 (min), respectively, where α is the myocardial electrical resistance coefficient, β is the PNS-SNS activity coefficient, ω is the the mean of metabolic perturbations, and σ is the standard deviation of metabolic perturbations. The figure below presents HR parameter values of m = 1 for f1y02, f1o04, chf01, s20641, and 0030.vf1 (900 samples).

Figure 8. Comparison of the restoration coefficient with the PSD in actual data. (A) The power spectral density (PSD) of 0030.vf1, f1y02, chf01, f1o04, and s20641 from the VT/VF, YOUNG, CHF, ELDER, and ST groups. (B) HR trajectories and the restoration coefficient (β) for 0030.vf1, f1y02, chf01, f1o04, and s20641. The restoration coefficient estimates were roughly consistent with their PSDs for f1y02, chf01, and f1o04, but significantly inconsistent with the PSD for s20641. Heart rate time series s20641 exhibited very low-frequency fluctuations, which was inconsistent with its very low LF/HF ratio but consistent with its small restoration coefficient.

Figure 9. Comparison of the restoration coefficient with the PSD in simulation data. The realizations of the PSD were obtained by the autoregressive model of order p = 16 at the time scale of m = 3 for the restoration coefficient β from 0.005 to 1.6. The original time series (m = 1) were generated by the NLARI process in Eq. (2), where the time delay κ2 = 1, the resistance coefficient α = 0.7786, the standard deviation σ = 0.0275, and the sample size T = 1500. It is seen that a small restoration coefficient \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta \in [0.01,0.1]$$\end{document}β∈[0.01,0.1] could correspond to a very low LF/HF ratio (s20641 with β = 0.0171); a relatively large restoration coefficient \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta \in (0.1,0.5)$$\end{document}β∈(0.1,0.5) could have a low LF/HF…

Figure 10. Relationship between the frequencies and the restoration coefficient. (A) The trajectories of the time series generated by the NLARI process on different time scales ranging from m = 1 to 50. As the time scale increases, the trajectories of the time series became more smooth. (B) A plot of frequency of the time series against the restoration coefficient β. The time series were generated by the NLARI process in Eq. (2) using the resistance coefficient α = 0.9489, β = 0.02i for i = 1, ⋯ , 50, the standard deviation σ = 0.022, the sample size 3000, and repeated times 3000. It is shown that the frequency of the time series was positively related to the restoration coefficient, but the frequency became insensitive to an increased restoration coefficient for an over enlarged time scale.

Figure 11. HR parameter alterations. Myocardial electrical resistance indicator α had an increase in the ELDER group, the largest increase in the CHF group, an ultraincrease in the SCD and VT/VF groups. PNS-SNS activity indicator β had a decrease in the ELDER group, a mild decrease compared to age in the CHF group, the largest decrease in the ST group, and the largest increase in the VT/VF group. HRV stability indicator γ had a decrease due to decreased PNS-SNS activity in the ELDER group, an overt increase due to ultraincreased myocardial resistance in the CHF group, the largest decrease due to ultradecreased PNS-NS activity in the ST group, and the largest increase due to ultraincreased myocardial resistance/PNS-SNS activity in the VT/VF group. The mean in myocardial metabolic rate ω had the largest increase in the VT/VF group leading to the largest increase HR amplitude indicator η1. The sd in myocardial metabolic rate σ had an ultradecrease in the EDLER group, the largest decrease in the CHF group, and the largest increase in the VT/VF group. HR amplitude indicator η2 had a decrease due to ultradecreased changes in myocardial metabolic rate in the ELDER group, the largest decrease due…