Stochastic Modeling and Optimal Time-Frequency Estimation of Task-Related HRV

: In this paper, we propose a novel framework for the analysis of task-related heart rate variability (HRV). Respiration and HRV are measured from 92 test participants while performing a chirp-breathing task consisting of breathing at a slowly increasing frequency under metronome guidance. A non-stationary stochastic model, belonging to the class of Locally Stationary Chirp Processes, is used to model the task-related HRV data, and its parameters are estimated with a novel inference method. The corresponding optimal mean-square error (MSE) time-frequency spectrum is derived and evaluated both with the individually estimated model parameters and the common process parameters. The results from the optimal spectrum are compared to the standard spectrogram with different window lengths and the Wigner-Ville spectrum, showing that the MSE optimal spectral estimator may be preferable to the other spectral estimates because of its optimal bias and variance properties. The estimated model parameters are considered as response variables in a regression analysis involving several physiological factors describing the test participants’ state of health, ﬁnding a correlation with gender, age, stress, and ﬁtness. The proposed novel approach consisting of measuring HRV during a chirp-breathing task, a corresponding time-varying stochastic model, inference method, and optimal spectral estimator gives a complete framework for the study of task-related HRV in relation to factors describing both mental and physical health and may highlight otherwise overlooked correlations. This approach may be applied in general for the analysis of non-stationary data and especially in the case of task-related HRV, and it may be useful to search for physiological factors that determine individual differences.


Introduction
Extensive research has been dedicated to understand and describe the complex interplay of factors that affect heart rate variability (HRV), the physiological phenomenon of the variation in the time interval between consecutive heartbeats. The dynamic task-specific vagal control of the heart represented by task-related HRV changes is especially valuable in exploring the flexible adaption of the organism to physical and mental challenges [1,2]. Age has been identified as the most critical variable affecting HRV, with the spectral power of HRV decreasing with increasing age [3][4][5]. Nevertheless, several other variables play a role, such as gender, body-mass index (BMI), anxiety and stress levels [6,7].
State-of-the-art spectral methods have been widely applied in the analysis of HRV, especially the periodogram and the Welch method [8]. These methods are suitable for HRV data measured in resting conditions, and for extracting the traditional high-frequency and low-frequency HRV inference method and the optimal time-frequency estimator, as well as the regression analysis approach. In Section 3, the results from the statistical inference are presented, the different time-frequency estimates are compared, and the significant predictors for each model parameter are discussed. Section 4 is dedicated to comments on the significance of this work and directions for further research.

Data Acquisition and Preprocessing
The study was conducted in correspondence with the Helsinki declaration and has been approved by the central ethical review board at Lund (Dnr 2013/754 and 2010/22). The TPs were volunteers, not using medicines or suffering from any disease affecting the cardiovascular system, and could discontinue their participation at any time. Each TP answered a questionnaire to collect information on their general health. In Table 1, the demographics of the TPs are summarized, showing the median, mean and standard deviation of the variables.
Electrocardiography (ECG) using disposable electrodes was used to record the heart rate, whereas respiration was measured with a respiratory belt including a piezo-electric device. The chirp-breathing task consisted of breathing following a metronome starting at frequency 0.12 Hz and increasing to 0.2 Hz. Both ECG and respiration were recorded at 1 kHz using the ML866 Power Lab data acquisition system and analyzed using its software LabChart8 (ADInstruments Pty Ltd., Dunedin, OTA, NZ) and MATLAB (Math-Works, Inc., Natick, MA, USA). The R-waves were detected with LabChart8 and the time difference between two consecutive heartbeats was used to derive the HRV from the heart rate.
The raw data sequences of both HRV and respiratory data were down-sampled to 4 Hz and the number of samples in each time-series is 840, corresponding to 210 seconds of recordings. An example of the respiratory signal and corresponding HRV is shown in Figure 1, where the frequency increase of the respiratory signal matches with a frequency increase in the HRV. For further processing, the data sequences are adjusted to zero mean.

Locally Stationary Chirp Processes
A zero mean real valued stochastic process , where E denotes the expected value and t 1 , t 2 time points in a time interval [T 0 , T f ] ⊆ R. A LSCP covariance is defined as where q is a non-negative function, also called "power schedule", r a stationary covariance function and f I a chirp covariance function defined by with m, d ∈ R. This means that the covariance C(t 1 , t 2 ) is modulated by the chirp f I (t 1 , t 2 ), whose frequency varies linearly with a delayed time axis t 1 +t 2 2 − d. If f I (t 1 , t 2 ) ≡ 1, then Equation (1) reduces to a LSP covariance [19].
To unequivocally identify an LSCP, the required functions and their parameters must be specified, and the resulting function C must be positive semidefinite. The proposed LSCP model for the task related HRV is identified by where t = t 1 +t 2 2 , τ = t 1 − t 2 and with a, b, c > 0, t ∈ [0, 210] s. The instantaneous power of a process X(t), defined as motivates the choice of the function q, see Figure 2. The value of the chirp parameters m, d and the parameters a, b, c must be estimated from the data.

Inference Method
Since a closed-form expression of the process distribution is not known, a maximum likelihood approach for estimating the model parameters is not feasible. In a previous contribution [23], we addressed the problem of fitting a real-valued LSP to sampled data by proposing the inference method HAnkel-Toeplitz Separation (HATS). The method hereby described is based on the same principle of dividing the inference problem into lower-dimensional sub-problems. The estimation procedure is summarized in Algorithm 1 and described in detail below.

Algorithm 1: Inference for LSCP parameters
Input : Zero mean time-series data, both respiratory signal and corresponding HRV. Output : Estimated chirp parametersm,d and estimated LSP parametersâ,b,ĉ.

Chirp estimation
1 Estimate the instantaneous frequency of the chirp from the the respiratory signal, by taking the frequency of maximum power in the chirp range from a TF representation of the signal; 2 Findm,d and the corresponding covariance matrixF I by fitting the linear chirp to the instantaneous frequency of the respiratory data; LSP estimation 3 Compute the instantaneous powerP from the HRV data; 4 Findâ,b by fitting the model function q toP; 5 Findĉ by fitting R •F I toR F =Ĉ SCM Q where R is the stationary covariance matrix of the model function r, C SCM is the estimated sample covariance matrix of the HRV data andQ is the corresponding Hankel matrix based on the model q with parametersâ,b; Return :m,d,â,b,ĉ.
The first step is estimating the instantaneous frequency of the chirp from the respiratory signal, by taking the frequency of maximum power in the chirp range (frequencies above 0.12 Hz) of the spectrogram with Hanning window of 128 samples (32 s). The parametersm andd are estimated through a robust least square fitting of the linear chirp mt − md to the estimated instantaneous frequency of the respiratory signal. A corresponding chirp covariance matrixF I is constructed according to the model Equation (2), where the covariance matrix values areF I (t 1 , t 2 ) = f I (t 1 , t 2 ). The estimation of the chirp from the respiratory signal of each TP allows accounting for the individual differences in the interpretation of the task of breathing accordingly to the metronome. An example of linear chirp fitted to the instantaneous frequency extracted from the spectrogram is presented in Figure 3.
At this point, a variation of the inference methods HATS can be used to compute the LSP parameters. First the function q in Equation (4) is fitted to the instantaneous power estimateP(t) = x 2 (t), where x(t) is the zero-mean adjusted HRV at time t, as in Equation (6), obtaining the estimated parametersâ andb, from which the Hankel matrixQ, with valuesQ (t 1 , t 2 ) =q ((t 1 + t 2 )/2), can be calculated. The next step is computing the ordinary sample covariance matrix of the HRV dataĈ SCM . Thereafter, the Hadamard division ofĈ SCM byQ provides the estimateR F =Ĉ SCM Q , where R F corresponds to a non-stationary covariance matrix with values R F (t 1 , . The parameterĉ is obtained by fitting the model stationary covariance matrix R obtained from Equation (5), Hadamard multiplied with the chirp matrixF I , toR F . Note that the division by the chirp matrix must be avoided for numerical issues as the chirp matrixF contains values close to 0, differently fromQ, which is a strictly positive matrix.

Mean Square Error Optimal Time-Frequency Kernel
The MSE optimal model based TF smoothing kernel can be derived and evaluated once the model parameters of the LSCP are estimated. In the following, F f denotes the Fourier transform of the function f .
The Wigner-Ville spectrum (WVS) of an LSP is defined as [ 19,20], and the corresponding ambiguity spectrum as Any TF representation member of the Cohen's class can be expressed as where Φ is an ambiguity kernel [11]. In [20], the MSE optimal kernel for circularly symmetric LSPs with Gaussian marginals has been derived as with The separable functions in Equation (11) are calculated according to the introduced model Equations (4) and (5) as The modification of the LSP kernel for LSCP is based on the fact that if a signal g 0 (t) is multiplied as g(t) = g 0 (t) · e i·( m 2 )t 2 then the WVS is coordinate transformed as W(t, ω) → W(t, ω − m · t). From [20], we accordingly find the optimal kernel to become The multitaper approach allows efficient estimation as a weighted sum of windowed spectrograms with weights w k , and windows h k (t), k = 1 . . . K, given by the eigenvalues α and eigenvectors q respectively solution to where the Hermitian rotated time-lag kernel is defined as with as in [20,21].

Regression Analysis
We performed a regression analysis [24] to search for the physiological factors that determine individual differences in the LSCP model parameters, exploring the predictive power of several physiological factors of interest: Gender (male, female), Age (years), Weight (kg), BMI (kg/m 2 ), State and Trait scores from the Spielberg State-Trait Anxiety Inventory (STAI) [25] and score in the Shirom-Melamed Burnout Questionnaire (SMBQ) [26][27][28]. Previous studies have validated the Swedish version of the STAI [29][30][31] and the SMBQ [32,33]. It is interesting considering both STAI and SMBQ measures since burnout consists of a combination of physical fatigue, emotional exhaustion, and cognitive weariness, not interchangeable with depression and anxiety [33,34].
Simple regression (i.e., regression models with a single explanatory variable) was tested first to isolate the effect of every factor, and then multivariate models were evaluated based on the statistical significance of the predictors (significance level of 0.1), the adjusted coefficient of determination R 2 adj , and the Akaike Information Criterion (AIC). Regression diagnostics included residual analysis, F-test for testing inclusion of variables, detection and treatment of outliers and influential observations.

Inference on the Model Parameters
The optimization steps in Algorithm 1 for estimating the parameters a, b, c were constrained so that a ∈ (0, 5], b ∈ (0, 1], c ∈ (0, 10]. In Table 2, we report the mean, median, standard deviation and coefficient of variation c v of the individually estimated LSCP parameters for the 92 TPs (LSCP i ). These values can be compared with the parameters estimated from a common LSCP (LSCP c ) with 92 realizations, giving one set of parameters for all TPs. These parameter estimates are similar to the medians and means of the LSCP i parameters, with the exception of the parameter b related to the rate of the instantaneous power decrease during the chirp-breathing task, whose mean is almost 3 times larger than the median. The parameters with larger variations are those that determine the strongest individual differences among TPs.

Comparison between Time-Frequency Estimates
The MSE optimal eigenvectors and eigenvalues, corresponding to multitapers and weights of the MSE optimal LSCP c spectral estimator, are shown in Figure 4. These will in general differ from the MSE optimal eigenvectors and eigenvalues based on the individually estimated parameters, used for computing the MSE optimal LSCP i spectral estimator, an example of which can be seen in Figure 5. In particular, the number of eigenvalues significantly different from zero varies depending on the parameters.
In Figure 6, we present examples of TF estimates for three TPs, obtained with the spectrogram with 64 samples Hanning window (16 s) S64, the spectrogram with 256 samples Hanning window (64 s) S256, the Wigner-Ville spectrum WVS, the MSE optimal LSCP c spectrum, and the MSE optimal LSCP i spectrum. The number of tapers K used for computing the LSCP i spectrum is dependent on the individual parameters. In the examples shown in Figure 6, K is equal to 10, 19, 3 for Ex. 1, Ex. 2, and Ex. 3 respectively.  The longer window spectrogram S256 offers a lower temporal resolution, compared to the shorter window spectrogram S64, while the opposite holds for the frequency resolution. Additionally to the resolution trade-off, the spectrogram suffers from a large variance, which is a significant problem when the estimate is based on a single realization of the underlying process, see the panels in the first two rows of Figure 6. The typical cross-terms of the WVS, large oscillating terms located in between the actual signal components, are visible in the panels in the middle row of Figure 6. By construction, the MSE optimal spectral estimators minimize bias and variance of the spectral estimates for LSCPs. LSCP i is preferred for capturing features that vary between the TPs, see panels in the last row of Figure 6. On the other hand, the LSCP c parameters are more robust than the individually estimated ones. This must be taken into account when the spectral estimates are used for extracting further features, e.g., frequency bands content from the time-frequency marginals or other spectral measures. It is worth noticing that for non-stationary signals of this kind, spectral measures such as the power in the low-and high-frequency bands and their ratio have limited interest, since the spectral content of the chirp signal is found above 0.12 Hz, which is above the usual low-frequency band. Figure 6. Examples of the TF estimates obtained for three TPs with the spectral estimators considered: S64 (first row); S256 (second row); WVS (third row); the MSE optimal LSCP c spectrum (fourth row); the MSE optimal LSCP i spectrum (fifth row). All the spectra for one TP are in the same column.

Regression Analysis
In the following, we report and discuss the regression analysis results for the LSCP parameters.

Parameter m
The chirp parameters m and d determine the linear chirp describing the instantaneous frequency increase of the respiratory signal. These values differ among the TPs as a result of the individual interpretation of the chirp-breathing task.
The significant predictors of m in simple regression models are Age, State, Trait, and SMBQ. When considering multiple regression models including the predictor Age, only SMBQ is still significant, whereas the variables State and Trait cease to be significant. On the other hand, the variable Gender becomes significant in models including Age. The best multiple regression model according to both AIC and R 2 adj values includes Gender, Age, and SMBQ and achieves R 2 adj = 0.63 (Table 3). According to this model, an increase in age corresponds to an increase in the value of m. An opposite effect on m is given by being male respect to female and having a higher SMBQ.

Parameter d
The significant covariates in simple regression models for the parameter d are the same as for m, a fact that relates to the correlation between the two chirp parameters. None of the variables is still significant in multivariate models including the variable Age, therefore the selected model for d includes only the variable Age, with R 2 adj = 0.59 (Table 4). For both m and d, even more explanatory than the numerical variable Age is the categorical variable Age Groups, that alone achieves R 2 adj = 0.79 and R 2 adj = 0.73 for m and d respectively. According to these models, being over the age 40 is what matters in affecting the value of the chirp parameters. However, only 10 TPs are older than 40 years, limiting the reliability of the last result.
To understand the meaning of the chirp parameters, a person older than 40 years would typically have the chirp parameter m approximately the double of the average values (i.e., m ≈ 0.0006) and the parameter d, representing the chirp delay, about half of the average value (i.e., d ≈ −200). These values correspond to an angular coefficient of the estimated linear chirp twice as large of the typical one, as a result of a steeper instantaneous frequency increase of the respiratory signal and consequently of the HRV.

Parameter a
The parameter a corresponds to the power at time zero, with a larger value of a corresponding to higher power, and it represents the amplitude multiplier that scales the exponential function. We consider the logarithmic transformation of a in the regression models to avoid positively skewed residuals.
Only the covariate Age is a significant predictor when considering models with a single explanatory variable. The regression model with only Age as predictor achieves R 2 adj = 0.24. The same R 2 adj is obtained with the categorical variable Age, with all Age groups significantly different from the baseline. The strong correlation of the scaling parameter a with the variable Age is expected, since it is known that the HRV power decreases with age.
In multivariate models, both Age and Stress play a role, and the best predictive model includes them as categorical variables with R 2 adj = 0.27 (Table 5). According to this model, the value of log(a) is affected similarly by being in the age group 30-40 years old, compared to the baseline 20-30 years old, and having an SMBQ value compatible with pre-ED.
The fact that both Age and Stress are significant predictors of the value of log(a) with a similar effect (negative slope and same scale) suggests an analogy between the effect of aging and higher stress levels on the HRV instantaneous power. However, being in the age group over 40 years old has a 3.5 times stronger effect than having an SMBQ value compatible with pre-ED.

Parameter b
The value of b is related to the rate of the power decrease consequent to the increase in breathing frequency. The estimates of the parameter b show considerable large variability among the TPs, as indicated in Table 2. A higher value of b corresponds to a faster power decrease and does not point to a fault in the estimation; however, this large variability prevented the regression models from achieving high values of R 2 adj . The only significant covariate in simple regression models is SMBQ. Step-wise model selection based on AIC led to a model including only SMBQ (Table 6), which however has a very low value of R 2 adj = 0.02, reflecting the fact that much of the variability of the parameter b remains unexplained. After outliers and influential observations treatment, leading to the removal of 4 TPs, also the variables Gender, Age and Stress become significant in simple regression. The best multivariate model after removal of the 4 outliers has the covariates Gender, Age and Stress, with only the category pre-ED significantly different from the baseline ( Table 7). This model has a higher R 2 adj = 0.11.

Parameter c
The parameter c relates to the local stationarity of the underlying stochastic process, with smaller values of c indicating a long-lasting autocorrelation. Conversely, a larger value of the parameter c corresponds to a smaller standard deviation of the Gaussian bell in Equation (5), meaning faster decaying autocorrelation of the process.
In simple regression, all the variables Gender, Age, Weight, BMI, Fit, State, Trait, SMBQ and Stress are significant predictors of the parameter c. The step-wise selected multivariate model includes Gender, Age, and Fit, with only the category Obese significantly different from the baseline category Normal Weight (Table 8). This model has a coefficient of determination R 2 adj = 0.25. The most significant and clear predictor is Age, which has an effect of decreasing the value of c of 0.063 for each year of age. Fitness levels and being Male respect to Female have stronger effect than Age but they are not equally significant. It should be noted that Gender and BMI are correlated, with men having higher BMI than women, and recent studies suggest that the traditional BMI value might not be a good indicator of the fitness of a person.

Conclusions
In this paper, we propose a time-varying stochastic model, based on the definition of LSCPs, for task-related HRV measured during a novel chirp-breathing task. The proposed inference method allows the evaluation of the MSE TF spectrum, derived from the proposed stochastic model, using either individual model parameters or the parameters estimated from a common process. For comparison, other TF spectral estimates are computed, namely the spectrogram with different window lengths and the WVS without any smoothing kernel. The MSE optimal spectral estimator evaluated with the estimated parameters minimizes bias and variance of the spectral estimates for LSCPs, and therefore its use is preferable for extracting further features from the spectral estimates.
The LSCP model parameters are used as response variables in a regression analysis using several physiological covariates as predictors. The regression analysis shows correlation of the LSCP parameters with gender, age, levels of stress and fitness. The presented results are consistent with those of more extensive studies examining the relationship between HRV and physiological variables. Since each model parameter relates to a different aspect of the underlying LSCP, this approach may be useful to search for physiological factors that determine individual differences in the HRV.
The proposed complete framework for the study of task-related HRV in relation to factors describing both mental and physical health has general validity for the analysis of non-stationary data, and especially in the case of task-related HRV.
Further research can be addressed to identify new spectral measures to be extracted from the improved TF estimates since the traditional estimation of the spectral power divided in the low-and the high-frequency band is of limited interest for time-varying signals. Acknowledgments: A preliminary version of this work was presented as part of the Ph.D. thesis "Statistical inference and time-frequency estimation for non-stationary signal classification" by Rachele Anderson, defended at Centre for Mathematical Sciences, Lund University, Sweden, on the 4th of October 2019. The content here published has been used according to Lund University policy.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations (in alphabetic order) are used in this manuscript: