Multivariable Fuzzy Measure Entropy Analysis for Heart Rate Variability and Heart Sound Amplitude Variability

Abstract

Simultaneously analyzing multivariate time series provides an insight into underlying interaction mechanisms of cardiovascular system and has recently become an increasing focus of interest. In this study, we proposed a new multivariate entropy measure, named multivariate fuzzy measure entropy (mvFME), for the analysis of multivariate cardiovascular time series. The performances of mvFME, and its two sub-components: the local multivariate fuzzy entropy (mvFEL) and global multivariate fuzzy entropy (mvFEG), as well as the commonly used multivariate sample entropy (mvSE), were tested on both simulation and cardiovascular multivariate time series. Simulation results on multivariate coupled Gaussian signals showed that the statistical stability of mvFME is better than mvSE, but its computation time is higher than mvSE. Then, mvSE and mvFME were applied to the multivariate cardiovascular signal analysis of R wave peak (RR) interval, and first (S1) and second (S2) heart sound amplitude series from three positions of heart sound signal collections, under two different physiological states: rest state and after stair climbing state. The results showed that, compared with rest state, for univariate time series analysis, after stair climbing state has significantly lower mvSE and mvFME values for both RR interval and S1 amplitude series, whereas not for S2 amplitude series. For bivariate time series analysis, all mvSE and mvFME report significantly lower values for after stair climbing. For trivariate time series analysis, only mvFME has the discrimination ability for the two physiological states, whereas mvSE does not. In summary, the new proposed mvFME method shows better statistical stability and better discrimination ability for multivariate time series analysis than the traditional mvSE method.

1. Introduction

Short-term, beat-to-beat cardiovascular variability reflects the inherent interactions from different components of the cardiovascular system and dynamic interplay between ongoing perturbations to the circulation and compensatory response of neurally mediated regulatory mechanisms [1]. Analysis of cardiovascular variability is a prerequisite for understanding the underlying signal generating mechanisms and detecting the cardiovascular diseases [2,3]. There is also an increasing interest in the application of cardiovascular variability monitoring to improve clinical outcomes [4].

Standard univariate time series analysis, typically as the heart rate variability (HRV), has been applied in a large variety of physiological measurements and clinical evaluations. Univariate time series analysis is only applicable if all the time series are statistically independent or uncorrelated at the very least, which is often not the case. It is well known that the human system is a complex system, and its behavior depends on the interactions of its components [5]. Multivariate analysis approach permits accounting for the dynamical relationships existing between variables [6,7,8]. The developments in sensor technology have also enabled routine recording of multivariate time series. Thus, there are substantial advantages in simultaneously analyzing multivariate time series observed from the cardiovascular system [8,9].

In recent years, entropy-based measures, such as the typical approximate entropy (ApEn) [10] and sample entropy (SampEn) [11], have been widely used in cardiovascular time series analysis. Entropy refers to the degree of regularity or irregularity of a time series and is estimated by counting how many ‘template’ patterns repeat. Repeated patterns imply increased regularity in the time series and lead to low entropy values [11,12,13,14]. SampEn is regarded as a modified version of ApEn to solve the shortcomings, such as bias and relative inconsistency [11]. The existing SampEn method was designed for the analysis of univariate time series, and is not able to reveal the dynamics across the channels measured in experimental and biological systems. When dealing with the multivariate time series, SampEn treats them as a set of individual time series by considering each variable separately. However, this is only applicable if all the data channels are statistically independent or uncorrelated at the very least (which is often not the case) [7,8]. The typical example is that only measuring the z coordinate of the Lorenz system cannot accurately reconstruct the dynamics of the Lorenz system since it does not resolve the x–y symmetry [8,15]. Meanwhile, there are substantial advantages in simultaneously analyzing several variables observed from the same phenomenon, especially if there is a large degree of uncertainty and coupling underlying the system dynamics or data acquisition. Thus, as the generalized form of SampEn, multivariate sample entropy (mvSE) was developed for the analysis of multivariate time series with the multiscale extension [7,8]. The mvSE method has received much attention in the biomedical and mechanical fields [16,17,18].

Although SampEn and its multivariate version of mvSE are powerful and popular algorithms, when applied to short time series, the results may be undefined or unreliable. This is because the vector similarity definitions in these two methods are based on Heaviside function, i.e., binary classification, which makes the boundary very rigid. Only the vectors within the threshold r are treated equally, whereas the vectors outside of this threshold r are ignored. This rigid boundary may induce abrupt changes of entropy values when the threshold r changes slightly, and it may even fail to define the entropy if no vector-matching could be found for the very small threshold r [11,19]. To enhance the statistical stability, fuzzy entropy (FuzzyEn) has been developed for univariate time series analysis by highlighting the notion of entropy with fuzzy theory [19,20]. FuzzyEn employed a fuzzy function to substitute the Heaviside function to make a gradually varied entropy value when threshold r monotonously changes. Its generalized multivariable form, multivariate fuzzy entropy (mvFE), has also been developed for multivariate time series analysis [18,21]. However, no matter whether it is for FuzzyEn or mvFE, the local vector similarity is exceedingly underlined [22]. Thus, both of them might give inaccurate results for the slow change signals since they both neglected the signal global characteristics. In our previous work, we developed a novel fuzzy measure entropy (i.e., FuzzyMEn) that combined both the local and global similarity, and FuzzyMEn has shown better algorithm discrimination ability than FuzzyEn [23,24]. In this study, we generalized the FuzzyMEn method for multivariate time series analysis by analyzing both simulation and cardiovascular multivariate time series. We defined the new multivariate entropy measure as multivariate fuzzy measure entropy (mvFME).

For multivariate cardiovascular time series analysis, existing studies included: the evaluation of the differences between HRV and blood pressure variability [3], the differences between systolic and diastolic interval variability [25], the relationship among the R wave peak (RR) interval, systolic and diastolic interval time series variability [26], the coupling between cardiac and respiratory systems through frequency analysis [27], the multi-site pulse oximeter data analysis [28], HRV, systolic interval variability and pulse transit time variability [29], the heart sound morphological variability [30], etc. Heart sounds, also termed as phonocardiogram (PCG), carry information about the mechanical activity of the cardiovascular system. Auscultation of the heart sounds is an essential part of the physical examination and is usually used as one of the first steps in evaluating the cardiovascular system in clinical practice [31,32]. Unlike the wide studies of HRV, there are few variability studies for heart sound signals. Xiao et al. used the regularity of the first (S1) heart sound amplitude to evaluate the practical cardiac contractility variability [33]. They also explored the variability of the S1 and second (S2) heart sound amplitude ratio (S1/S2) [34]. However, the interpretation of the heart sound amplitude changes within different physiological states and the simultaneous monitoring of multivariate cardiovascular time series variability is still an open problem. Thus, in this study, we performed the multivariate analysis for three cardiovascular time series: the RR interval time series from electrocardiogram (ECG) and the S1 and S2 amplitude series from PCG, by using both the traditional mvSE method and newly developed mvFME method.

The rest of the paper is organized as follows. Section 2 gives the definition of mvSE and mvFME to allow the detailed comparison and inspection of these two multivariate entropy measures to be observed. Section 3 discusses the experimental design where the multivariate simulation signals and cardiovascular time series (RR interval, S1 and S2 amplitude series) from 20 healthy subjects were constructed for the multivariate entropy analysis. In Section 4, the results for both simulation and cardiovascular time series are provided. Finally, Section 5 draws the discussions and identifies the future work.

2. Multivariate Entropy Measures

In this section, we first give a brief introduction to mvSE, and then describe the definition of mvFME.

2.1. Multivariate Sample Entropy (mvSE)

Since introduced by Richman and Moorman in 2000 [11], the SampEn method has been widely applied for the univariate short-term physiological time series analysis. Recently, Ahmed and Mandic developed SampEn for the multivariate analysis and produced the mvSE method. The mvSE analysis is performed through the following steps [8,9]:

(1)

For a p-variate time series ${\left\{x_{k,i}\right\}}_{i=1}^N$, $k=1,2,\cdots ,p$, where $N$ is the number of samples in each variate, firstly normalize each time series for $k=1,2,\cdots ,p$, and then form the composite delay vector using a composite delay factor based on the multivariate embedded reconstruction:

$$X_i^m=\left[x_{1,i},x_{1,i+{\tau }_1},\cdots ,x_{1,i+\left(m_1-1\right){\tau }_1},x_{2,i},x_{2,i+{\tau }_2},\cdots ,x_{2,i+\left(m_2-1\right){\tau }_2},\cdots ,x_{p,i},x_{p,i+{\tau }_p},\cdots ,x_{p,i+\left(m_p-1\right){\tau }_p}\right],$$

(1)

where $M=\left[m_1,m_2,\cdots ,m_p\right]$ is the embedding vector, $\tau =\left[{\tau }_1,{\tau }_2,\cdots ,{\tau }_p\right]$ is the time lag vector, and $X_i^m\in R^m$ is the $m$ dimension composite delay vector with $m={\displaystyle {\sum }_{k=1}^p{m}_k}$, and $i=1,2,\cdots ,N-n$ with $n=\max \left(M\right)\times \max \left(\tau \right)$.

(2)

Define the distance $d_{i,j}^m$ between any two composite delay vectors $X_i^m$ and $X_j^m$ as the maximum norm, that is,

$$d_{i,j}^m=d[X_i^m,X_j^m]=\underset{l=1,2,\cdots ,m}{\max }\left(\left|x\left(i+l-1\right)-x\left(j+l-1\right)\right|\right).$$

(2)

(3)

For a given vector $X_i^m$ and a threshold $r$, count the number of instances $P_i$ where $d_{i,j}^m\le r$, $j\ne i$, and then calculate the frequency of occurrence, $B_i^m\left(r\right)={\left(N-n-1\right)}^{-1}{P}_i$, and define a global quantity $B^m\left(r\right)={\left(N-n\right)}^{-1}{\displaystyle \sum }_{i=1}^{N-n}{B}_i^m\left(r\right)$.

(4)

Extend the dimensionality of the multivariate delay vector in Equation (1) from $m$ to $m+1$. This can be performed in $p$ different ways, as the system can evolve to any space with $M=\left[m_1,m_2,\cdots ,m_k+1,\cdots ,m_p\right]$ ($k=1,2,\cdots ,p$). Thus, a total of $p$ vectors $X_i^{m+1}\in R^{m+1}$ are obtained. The k-th vector $X_i^{m+1}$ is:

$$\begin{gathered} X_i^{m+1}=[x_{1,i},x_{1,i+{\tau }_1},\cdots ,x_{1,i+\left(m_1-1\right){\tau }_1},x_{2,i},x_{2,i+{\tau }_2},\cdots ,x_{2,i+\left(m_2-1\right){\tau }_2},\cdots ,x_{k,i},x_{k,i+{\tau }_k}, \\ \cdots ,x_{k,i+\left(m_k-1\right){\tau }_k},x_{k,i+m_k{\tau }_k},\cdots ,x_{p,i},x_{p,i+{\tau }_p},\cdots ,x_{p,i+\left(m_p-1\right){\tau }_p}]. \end{gathered}$$

(3)

(5)

For a given $X_i^{m+1}$, count the number of instances $Q_i$, where $d_{i,j}^{m+1}\le r$, $j\ne i$, and then calculate the frequency of occurrence, $B_i^{m+1}\left(r\right)={\left(p\left(N-n\right)-1\right)}^{-1}{Q}_i$, and define $B^{m+1}\left(r\right)={\left(p\left(N-n\right)\right)}^{-1}{\displaystyle \sum }_{i=1}^{p\left(N-n\right)}{B}_i^{m+1}\left(r\right)$.

(6)

Finally, mvSE is defined by

$$\mathrm{mvSE}\left(M,\tau ,r,N\right)=-\ln \left[\frac{B^{m+1}\left(r\right)}{B^m\left(r\right)}\right].$$

(4)

2.2. Multivariate Fuzzy Measure Entropy (mvFME)

mvFME is generalised from the FuzzyMEn method [23,24]. The calculation process of mvFME was summarized as follows:

(1)

For a $p$-variate time series ${\left\{x_{k,i}\right\}}_{i=1}^N$, $k=1,2,\cdots ,p$, where $N$ is the number of samples in each variate, firstly normalize each time series for $k=1,2,\cdots ,p$, and then form the local composite delay vector and global composite delay vector using a composite delay factor based on the multivariate embedded reconstruction:

$${\overline{X}}_i^m=\left[{\overline{x}}_{1,i},{\overline{x}}_{1,i+{\tau }_1},\cdots ,{\overline{x}}_{1,i+\left(m_1-1\right){\tau }_1},{\overline{x}}_{2,i},{\overline{x}}_{2,i+{\tau }_2},\cdots ,{\overline{x}}_{2,i+\left(m_2-1\right){\tau }_2},\cdots ,{\overline{x}}_{p,i},{\overline{x}}_{p,i+{\tau }_p},\cdots ,{\overline{x}}_{p,i+\left(m_p-1\right){\tau }_p}\right],$$

(5)

$${\stackrel{=}{X}}_i^m=\left[{\stackrel{=}{x}}_{1,i},{\stackrel{=}{x}}_{1,i+{\tau }_1},\cdots ,{\stackrel{=}{x}}_{1,i+\left(m_1-1\right){\tau }_1},{\stackrel{=}{x}}_{2,i},{\stackrel{=}{x}}_{2,i+{\tau }_2},\cdots ,{\stackrel{=}{x}}_{2,i+\left(m_2-1\right){\tau }_2},\cdots ,{\stackrel{=}{x}}_{p,i},{\stackrel{=}{x}}_{p,i+{\tau }_p},\cdots ,{\stackrel{=}{x}}_{p,i+\left(m_p-1\right){\tau }_p}\right],$$

(6)

where the local composite delay vector ${\overline{X}}_i^m$ represents $m$ consecutive values ${\overline{x}}_{k,i}$; however, removing the local baseline of ${\overline{x}}_{k,i}=x_{k,i}-\frac{1}{m_k}{{\displaystyle \sum }}_{t=1}^{m_k}{x}_{k,i+\left(t-1\right){\tau }_k}$, the global composite delay vector ${\stackrel{=}{X}}_i^m$ represents $m$ consecutive values ${\stackrel{=}{x}}_{k,i}$ but removing global mean value of ${\stackrel{=}{x}}_{k,i}=x_{k,i}-\frac{1}{N}{{\displaystyle \sum }}_{t=1}^N{x}_{k,t}$, $M=\left[m_1,m_2,\cdots ,m_p\right]$ is the embedding vector, $\tau =\left[{\tau }_1,{\tau }_2,\cdots ,{\tau }_p\right]$ is the time lag vector, $m={\displaystyle {\sum }_{k=1}^p{m}_k}$, $i=1,2,\cdots ,N-n$ and $n=\max \left(M\right)\times \max \left(\tau \right)$.

(2)

Define the distance $d{L}_{i,j}^m$ between any two local composite delay vectors ${\overline{X}}_i^m$ and ${\overline{X}}_j^m$, and the distance $d{G}_{i,j}^m$ between any two global composite delay vectors ${\stackrel{=}{X}}_i^m$ and ${\stackrel{=}{X}}_j^m$, as the maximum norm, that is,

$$d{L}_{i,j}^m=d\left[{\overline{X}}_i^m,{\overline{X}}_j^m\right]=\underset{l=1,2,\cdots ,m}{\max }\left(\left|\overline{x}\left(i+l-1\right)-\overline{x}\left(j+l-1\right)\right|\right),$$

(7)

$$d{G}_{i,j}^m=d\left[{\stackrel{=}{X}}_i^m,{\stackrel{=}{X}}_j^m\right]=\underset{l=1,2,\cdots ,m}{\max }\left(\left|\stackrel{=}{x}\left(i+l-1\right)-\stackrel{=}{x}\left(j+l-1\right)\right|\right),$$

(8)

(3)

Given the parameters’ local vector similarity weight $n_L$, local tolerance threshold $r_L$, global vector similarity weight $n_G$ and global tolerance threshold $r_G$, calculate the similarity degree $D{L}_{i,j}^m\left(n_L,r_L\right)$ between the local composite delay vectors ${\overline{X}}_i^m$ and ${\overline{X}}_j^m$ by the fuzzy function $\mu L(d{L}_{i,j}^m,n_L,r_L)$, as well as calculate the similarity degree $D{G}_{i,j}^m\left(n_L,r_L\right)$ between the global composite delay vectors ${\stackrel{=}{X}}_i^m$ and ${\stackrel{=}{X}}_j^m$ by the fuzzy function $\mu G(d{G}_{i,j}^m,n_G,r_G)$:

$$D{L}_{i,j}^m\left(n_L,r_L\right)=\mu L\left(d{L}_{i,j}^m,n_L,r_L\right)=\exp \left(-\frac{{\left(d{L}_{i,j}^m\right)}^{n_L}}{r_L}\right),$$

(9)

$$D{G}_{i,j}^m\left(n_G,r_G\right)=\mu G\left(d{G}_{i,j}^m,n_G,r_G\right)=\exp \left(-\frac{{\left(d{G}_{i,j}^m\right)}^{n_G}}{r_G}\right).$$

(10)

Define the function $\varnothing L^m\left(n_L,r_L\right)$ and $\varnothing G^m\left(n_G,r_G\right)$ as:

$$\varnothing L^m\left(n_L,r_L\right)=\frac{1}{N-n}{{\displaystyle \sum }}_{i=1}^{N-n}(\frac{1}{N-n-1}{{\displaystyle \sum }}_{j=1,j\ne i}^{N-n}D{L}_{i,j}^m\left(n_L,r_L\right)),$$

(11)

$$\varnothing G^m\left(n_G,r_G\right)=\frac{1}{N-n}{{\displaystyle \sum }}_{i=1}^{N-n}(\frac{1}{N-n-1}{{\displaystyle \sum }}_{j=1,j\ne i}^{N-n}D{G}_{i,j}^m\left(n_G,r_G\right)).$$

(12)

(4)

Extend the dimensionality of the multivariate delay vectors in Equations (5) and (6) from $m$ to $m+1$. This can be performed in $p$ different ways, as the system can evolve to any space with $M=\left[m_1,m_2,\cdots ,m_k+1,\cdots ,m_p\right]$ ($k=1,2,\cdots ,p$). Thus, a total of $p$ vectors ${\overline{X}}_i^{m+1}\in R^{m+1}$ and a total of $p$ vectors ${\stackrel{=}{X}}_i^{m+1}\in R^{m+1}$ are obtained. The kth vectors ${\overline{X}}_i^{m+1}$ and ${\stackrel{=}{X}}_i^{m+1}$ are respectively:

$$\begin{gathered} {\overline{X}}_i^{m+1}=[{\overline{x}}_{1,i},{\overline{x}}_{1,i+{\tau }_1},\cdots ,{\overline{x}}_{1,i+\left(m_1-1\right){\tau }_1},{\overline{x}}_{2,i},{\overline{x}}_{2,i+{\tau }_2},\cdots ,{\overline{x}}_{2,i+\left(m_2-1\right){\tau }_2},\cdots ,{\overline{x}}_{k,i},{\overline{x}}_{k,i+{\tau }_k}, \\ \cdots ,{\overline{x}}_{k,i+\left(m_k-1\right){\tau }_k},{\overline{x}}_{k,i+m_k{\tau }_k},\cdots ,{\overline{x}}_{p,i},{\overline{x}}_{p,i+{\tau }_p},\cdots ,{\overline{x}}_{p,i+\left(m_p-1\right){\tau }_p}] \end{gathered}$$

(13)

$$\begin{gathered} {\stackrel{=}{X}}_i^{m+1}=[{\stackrel{=}{x}}_{1,i},{\stackrel{=}{x}}_{1,i+{\tau }_1},\cdots ,{\stackrel{=}{x}}_{1,i+\left(m_1-1\right){\tau }_1},{\stackrel{=}{x}}_{2,i},{\stackrel{=}{x}}_{2,i+{\tau }_2},\cdots ,{\stackrel{=}{x}}_{2,i+\left(m_2-1\right){\tau }_2},\cdots ,{\stackrel{=}{x}}_{k,i},{\stackrel{=}{x}}_{k,i+{\tau }_k}, \\ \cdots ,{\stackrel{=}{x}}_{k,i+\left(m_k-1\right){\tau }_k},{\stackrel{=}{x}}_{k,i+m_k{\tau }_k},\cdots ,{\stackrel{=}{x}}_{p,i},{\stackrel{=}{x}}_{p,i+{\tau }_p},\cdots ,{\stackrel{=}{x}}_{p,i+\left(m_p-1\right){\tau }_p}] \end{gathered}$$

(14)

where ${\overline{x}}_{k,i}$ and ${\stackrel{=}{x}}_{k,i}$ have the same meanings as Equations (5) and (6).

(5)

Similarly, define the function $\varnothing L^{m+1}\left(n_L,r_L\right)$ for the local composite delay vectors ${\overline{X}}_i^{m+1}$ and ${\overline{X}}_j^{m+1}$ and the function $\varnothing G^{m+1}\left(n_G,r_G\right)$ for the global composite delay vectors ${\stackrel{=}{X}}_i^{m+1}$ and ${\stackrel{=}{X}}_j^{m+1}$ as:

$$\varnothing L^{m+1}\left(n_L,r_L\right)=\frac{1}{p\left(N-n\right)}{{\displaystyle \sum }}_{i=1}^{p\left(N-n\right)}(\frac{1}{p\left(N-n\right)-1}{{\displaystyle \sum }}_{j=1,j\ne i}^{p\left(N-n\right)}D{L}_{i,j}^{m+1}\left(n_L,r_L\right)),$$

(15)

$$\varnothing G^{m+1}\left(n_G,r_G\right)=\frac{1}{p\left(N-n\right)}{{\displaystyle \sum }}_{i=1}^{p\left(N-n\right)}(\frac{1}{p\left(N-n\right)-1}{{\displaystyle \sum }}_{j=1,j\ne i}^{p\left(N-n\right)}D{G}_{i,j}^{m+1}\left(n_G,r_G\right)).$$

(16)

(6)

Then, the local multivariate fuzzy entropy (mvFEL) and global multivariate fuzzy entropy (mvFEG) are defined by:

$$\mathrm{mvFEL}\left(M,\tau ,n_L,r_L,N\right)=-\ln \left[\varnothing L^{m+1}\left(n_L,r_L\right)/\varnothing L^m\left(n_L,r_L\right)\right],$$

(17)

$$\mathrm{mvFEG}\left(M,\tau ,n_G,r_G,N\right)=-\ln \left[\frac{\varnothing G^{m+1}\left(n_G,r_G\right)}{\varnothing G^m\left(n_G,r_G\right)}\right].$$

(18)

(7)

Finally, mvFME is defined by

$$\mathrm{mvFME}\left(M,\tau ,n_L,r_L,n_G,r_G,N\right)=\mathrm{mvFEL}\left(M,\tau ,n_L,r_L,N\right)+\mathrm{mvFEG}\left(M,\tau ,n_G,r_G,N\right).$$

(19)

In this study, we used the following parameter setting as suggested in our previous studies [22,35]: the local similarity weight n_L = 2 and global vector similarity weight n_G = 2, and the local tolerance threshold r_L was set equal to the global threshold r_G, i.e., r_L = r_G = r. Thus, Equation (19) becomes:

$$\mathrm{mvFME}\left(M,\tau ,r,N\right)=\mathrm{mvFEL}\left(M,\tau ,r,N\right)+\mathrm{mvFEG}\left(M,\tau ,r,N\right).$$

(20)

For both mvSE and mvFME, the entropy results were only based on the four parameters: embedding dimension vector $M$, time lag vector $\tau$, tolerance threshold $r$ and time series length $N$. In addition, if we set $p=1$, the multivariate entropy method becomes the univariate analysis. Thus, mvSE is the same as the traditional SampEn method and mvFME is the same as the traditional FuzzyMEn method.

3. Experiment Design

To test the practical applications for these two multivariate entropy measures mvSE and mvFME, we compared their performances on both simulation and real cardiovascular signals. Since the mvFME is the sum of two independent sub-components: mvFEL (considering only the local composite vectors) and mvFEG (considering only the global composite vectors), we also included the performances of these two measures: mvFEL and mvFEG.

3.1. Simulation Signals

Coupled Gaussian noises were used as the simulation multivariate time series. The coupled Gaussian noise model described in [36] was used. Firstly, four independent Gaussian noises $n_1$, $n_2$, $n_3$ and $n_4$ were generated. Then, the Gaussian noise $n_1$ was used as the public sector and the three other Gaussian noises $n_2$, $n_3$ and $n_4$ were mixed with the Gaussian noise $n_1$ as private sectors into three coupled Gaussian noises $x$, $\mathrm{y}$ and $z$ by

$$\{\begin{array}{c}x=c\times n_1+\left(1-c\right)\times n_2\\ y=c\times n_1+\left(1-c\right)\times n_{3,}\\ z=c\times n_1+\left(1-c\right)\times n_4\end{array}$$

(21)

where $c$ is the coupling degree parameter. $c=0$ means the three multivariate time series $x$, $\mathrm{y}$ and $z$ are totally independent and $c=1$ means that they are the same Gaussian noise $n_1$.

3.2. Cardiovascular Signals

To test the practical applications of the multivariate entropy measures, we compared their performances on the multivariate cardiovascular time series: RR interval, S1 and S2 heart sound amplitude series. Twenty healthy young male subjects (24.2 ± 1.9 years) were enrolled in this study. All subjects gave their written informed consent to participate in the study, and confirmed that they had not participated in any other “clinical trial” within the previous three months. The study obtained a full approval from the Clinical Ethics Committee of the Second Hospital of Dalian Medical University (201203267), and all clinical investigations were conducted according to the principles of expressed in the Declaration of Helsinki.

Table 1. Basic clinical characteristics of all 20 subjects.

Variable	Value
Age (year)	24.2 ± 1.9
Height (cm)	174 ± 4
Weight (kg)	64 ± 7
Heart rate (beats/min)	69 ± 9
Systolic blood pressure (mmHg)	121 ± 9
Diastolic blood pressure (mmHg)	65 ± 7

Note: data are expressed as mean ± standard deviation (SD).

depicts the details for the involved subjects.

Data collection was performed in a quiet room. For each subject, three PCG signals were simultaneously recorded with the standard limb lead-I ECG using microphone sensors (MLT201, AD instrument, Bella Vista, Australia) at three common auscultation positions: aortic area, mitral area and tricuspid area [37,38]. The signal sample rate was 2000 Hz. ECG and PCG signals were firstly recorded for about 5 min at rest state. Then, the subject was asked to do the activity of stair climbing for 120 stairs in an office building from the seventh floor to the twelfth floor. The stair climbing usually lasts about 100–120 s. After stair climbing, the subject was asked immediately to lay on his back in an examination bed and start the signal recording for another 5 min. The multivariate entropy values were compared between these two physiological states, i.e., rest state and after stair climbing state.

Before signal recording, each subject had a rest for 10 min to permit the cardiovascular stabilization. Heart rate (HR) and blood pressure (BP) values of each subject, including systolic blood pressure (SBP) and diastolic blood pressure (DBP), were measured by an automatic electronic sphygmomanometer (HEM-7200, OMRON, Osaka, Japan) at the left brachial artery three times, i.e., before the signal recording, during the gap between the rest state and after stair climbing state measurements, and after the signal recording. HR and BP values from the three measurements were averaged to obtain the final results for each subject (see Table 1).

Signal processing was performed offline. ECG signal was firstly filtered using a 0.05–40 Hz band-pass filter. The Pan and Tompkins method [39] was used to locate the R peaks, and thus the RR interval time series were constructed. RR intervals with ectopic beats were detected and excluded using the combination method [40]. PCG signals were firstly filtered using a 20–200 Hz band-pass filter. Then, Springer’s hidden semi-Markov model (HSMM) segmentation method [41] was used to segment each PCG signal to generate the time durations for S1 and S2 heart sounds. Figure 1 shows waveform examples of the recorded ECG and PCG signals, and the corresponding R peak location, S1 and S2 heart sounds segmentation results. S1 and S2 heart sound amplitude series were constructed by calculating the PCG signal amplitudes in each S1 and S2 heart sound state. Figure 2 shows examples of RR intervals, and S1 and S2 amplitude series from the ECG and aortic PCG signals, from both rest and after stair climbing states. For each time series, only the first 200 beats were shown.

3.3. Statistical Analysis

Mean ± standard deviation (SD) of the four multivariate entropy measures were obtained across all 20 subjects. The results of mvSE, mvFEL, mvFEG and mvFME were firstly tested as normal distribution by the Kolmogorov–Smirnov test. If the entropy results met the normal distribution, the paired t-test was used to test the statistical difference between the rest and after climb states. If not, non-parametric test was used. All statistical analyses were performed using the SPSS software (Version 20, IBM, New York, NY, USA). Statistical significance was set a priori at p < 0.05.

4. Results

4.1. Results on Simulation Signals

First, we tested the change of multivariate entropy measures mvSE, mvFEL, mvFEG and mvFME values with the increase of coupling degree parameter $c$ by using different numbers of coupled Gaussian noises, i.e., univariate analysis ($p=1$) using time series $x$ only, bivariate analysis ($p=2$) using time series $x$ and $y$, and trivariate analysis ($p=3$) using time series $x$, $y$ and $z$. For $p$ from 1 to 3, the settings for $M$ are$M=\left[2\right]$, $M=\left[2,2\right]$ and $M=\left[2,2,2\right]$, for $\tau$ are $\tau =\left[1\right]$, $\tau =\left[1,1\right]$ and $\tau =\left[1,1,1\right],$ respectively. The setting for $r$ is $r=0.15$ for all situations. Simulation time series length is set as N = 300. To eliminate the influence of random factor, for each coupling degree $c$ under each multivariate analysis, 100 realizations were generated and the mean values were used as the final results.

Figure 3 shows the mean ± SD results of mvSE, mvFEL, mvFEG and mvFME from 100 repeats when coupling degree c increases. For univariate series analysis, the coupling degree parameter c does not influence the entropy results of the time series x since x is the superposition of two independent Gaussian noises and is still a pure Gaussian noise. Thus, the mvSE, mvFEL, mvFEG and mvFME results do not change with the increase of coupling degree c. In this case, mvSE, mvFEL, mvFEG and mvFME all measure the inherent complexity of the Gaussian noise. It is clear that the SD values in mvFEL, mvFEG and mvFME are much lower than those in mvSE, indicating the statistical stability of the fuzzy function-based entropy method. For bivariate and trivariate analyses, all mvSE, mvFEL, mvFEG and mvFME results show increase trends with the increase of c. However, the increases in mvFEL, mvFEG and mvFME are strictly monotonous and the trend curves are smooth, whereas mvSE shows ups and downs with the local increase of c. Again, the SD values in mvFEL, mvFEG and mvFME are much lower than those in mvSE, confirming the statistical stability of the new proposed mvFME method. It is also worth noting that for analyzing the bivariate and trivariate Gaussian time series, with the increase of coupling degree c, the monotonous increase trend in mvFME is much more obvious than those in mvFEL and mvFEG, confirming the benefit of combining both the local (in mvFEL) and global (in mvFEG) vector similarity to make the new mvFME method and again verifying the better algorithm discrimination ability of mvFME than the single versions of mvFEL and mvFEG.

Table 2. Comparison of the computation time for the four multivariate entropy measures (multivariate Sample Entropy (mvSE), local multivariate fuzzy entropy (mvFEL), global multivariate fuzzy entropy (mvFEG) and multivariate fuzzy measure entropy (mvFME)) on simulation Gaussian time series.

Gaussian Time Series	Time (s)
	mvSE	mvFEL	mvFEG	mvFME
Univariate analysis	6.58	16.67	16.61	33.28
Bivariate analysis	19.42	51.20	51.03	102.23
Trivariate analysis	44.01	108.84	108.47	217.31

shows the computation time for the four multivariate entropy measures (mvSE, mvFEL, mvFEG and mvFME). The results from the 100 realizations of Gaussian time series for all varied coupling degree $c$ situations. mvSE reports the smallest computation time for each of the three multivariate types (univariate, bivariate and trivariate). As expected, mvFEL and mvFEG have similar computation times, and the computation time of mvFME method is the sum of mvFEL and mvFEG methods. The calculation was performed using MATLAB software (Version R2009a) on Windows XP platform (CPU: Intel Core i5, 2.66 GHz).

4.2. Results on Cardiovascular Signals

We set $p=1$ to measure the univariate results of mvSE, mvFEL, mvFEG and mvFME for the RR interval, S1 amplitude and S2 amplitude series, respectively, set $p=2$ to measure the bivariate mvSE, mvFEL, mvFEG and mvFME for any two combinations of these three time series, and set $p=3$ to measure the trivariate mvSE, mvFEL, mvFEG and mvFME for three series. The settings of $M$, $\tau$ and $r$ are the same as Section 3.1. Time series length is the actual beat number from the 5-min measurement.

All mvSE, mvFEL, mvFEG and mvFME results, from both rest and after stair climbing states, had normal distribution from the Kolmogorov–Smirnov test. Figure 4 gives their statistical plots by analyzing the univariate, bivariate and trivariate cardiovascular time series, respectively.

Table 3. Statistical results of mvSE and mvFME between the rest and after stair climbing states by analyzing the univariate, bivariate and trivariate cardiovascular time series from electrocardiogram (ECG) and three phonocardiogram (PCG) signals, respectively.

Signals	Time Series	mvSE			mvFME
		Rest	Climb	p-Value	Rest	Climb	p-Value
ECG + aortic PCG	RR	2.15 ± 0.38	1.11 ± 0.30	6 × 10−9	2.17 ± 0.34	0.93 ± 0.37	3 × 10−10
	S1	2.26 ± 0.26	1.53 ± 0.37	6 × 10−7	2.60 ± 0.18	1.65 ± 0.50	3 × 10−7
	S2	2.33 ± 0.32	2.16 ± 0.29	0.1	2.59 ± 0.23	2.32 ± 0.31	3 × 10−3
	RR & S1	1.52 ± 0.20	1.13 ± 0.13	2 × 10−6	1.81 ± 0.13	1.16 ± 0.19	9 × 10−11
	RR & S2	1.46 ± 0.19	1.22 ± 0.19	9 × 10−4	1.78 ± 0.19	1.32 ± 0.11	1 × 10−8
	S1 & S2	1.60 ± 0.25	1.16 ± 0.22	6 × 10−6	2.00 ± 0.21	1.38 ± 0.30	2 × 10−9
	RR & S1 & S2	1.04 ± 0.42	1.06 ± 0.09	0.8	1.50 ± 0.21	1.09 ± 0.16	1 × 10−8
ECG + mitral PCG	RR	2.15 ± 0.38	1.11 ± 0.30	6 × 10−9	2.17 ± 0.34	0.93 ± 0.37	3 × 10−10
	S1	2.16 ± 0.42	1.66 ± 0.62	5 × 10−3	2.43 ± 0.25	1.81 ± 0.59	2 × 10−4
	S2	2.27 ± 0.38	2.01 ± 0.57	0.1	2.50 ± 0.34	2.27 ± 0.48	0.1
	RR & S1	1.52 ± 0.14	1.14 ± 0.16	2 × 10−8	1.80 ± 0.20	1.18 ± 0.20	3 × 10−11
	RR & S2	1.44 ± 0.25	1.18 ± 0.15	5 × 10−4	1.76 ± 0.23	1.32 ± 0.14	3 × 10−6
	S1 & S2	1.52 ± 0.26	1.20 ± 0.43	9 × 10−3	1.96 ± 0.27	1.46 ± 0.44	3 × 10−4
	RR & S1 & S2	0.98 ± 0.27	1.06 ± 0.17	0.3	1.50 ± 0.25	1.12 ± 0.19	2 × 10−5
ECG + tricuspid PCG	RR	2.15 ± 0.38	1.11 ± 0.30	6 × 10−9	2.17 ± 0.34	0.93 ± 0.37	3 × 10−10
	S1	2.19 ± 0.37	1.75 ± 0.48	7 × 10−3	2.47 ± 0.37	1.92 ± 0.49	2 × 10−3
	S2	2.25 ± 0.29	2.05 ± 0.42	0.1	2.47 ± 0.32	2.37 ± 0.39	0.3
	RR & S1	1.47 ± 0.23	1.17 ± 0.16	2 × 10−6	1.81 ± 0.22	1.22 ± 0.16	2 × 10−9
	RR & S2	1.47 ± 0.18	1.16 ± 0.13	1 × 10−6	1.72 ± 0.21	1.33 ± 0.12	3 × 10−8
	S1 & S2	1.54 ± 0.23	1.27 ± 0.34	9.9 × 10−3	1.93 ± 0.31	1.58 ± 0.33	2 × 10−3
	RR & S1 & S2	1.09 ± 0.34	1.10 ± 0.14	0.97	1.47 ± 0.24	1.18 ± 0.15	1 × 10−4

Note: data are expressed as mean ± standard deviation (SD). Rest: rest state; Climb: after stair climbing state; RR: RR interval time series; S1: S1 heart sound amplitude series; S2: S2 heart sound amplitude series.

shows the numerical values for mvSE and mvFME methods. For cardiovascular time series from ECG and aortic PCG signals, paired t-test results showed that all mvSE and mvFME values in the after stair climbing state were significantly lower than those in the rest state (all p < 0.01, except mvSE of univariate S2 amplitude series with p = 0.1 and mvSE of trivariate RR interval and S1 amplitude and S2 amplitude series with p = 0.8). The results from ECG and mitral PCG signals, and from ECG and tricuspid PCG signals, show similar results compared with the results from ECG and aortic PCG signals. The difference is only that the decreases in mvFME in after stair climbing state were not statistically significant for univariate S2 amplitude series (p = 0.1 and p = 0.3, respectively). It is worth noting that the differences between the rest and after stair climbing states in trivariate RR interval and S1 amplitude and S2 amplitude series from each of the three different position PCG signals were identified as statistically significant by mvFME, whereas not by mvSE. The mean mvSE and mvFME values decreased from using univariate to multivariate time series for both rest and after stair climbing states. In addition, the mvFEL and mvFEG reported similar trends with the mvFME method.

5. Discussions

This study proposed a new multivariate entropy measure, named multivariate fuzzy measure entropy (mvFME), for the analysis of cardiovascular multivariate time series. As an application, the new proposed mvFME, as well as mvSE, has been applied for the multivariate cardiovascular time series to detect their discrimination ability between two physiological states. Short-term, beat-to-beat cardiovascular variability reflects the dynamic interplay between ongoing perturbations to the circulation and the compensatory response of neurally mediated regulatory mechanisms. While univariate time series analysis may be employed to quantify the variability itself, the bivariate and trivariate time series analysis permits the dynamic characterization of the cardiovascular regulatory mechanisms. Compared with univariate time series, bivariate and trivariate time series analysis may be even more illuminating, as it can provide a quantitative characterization of the cardiovascular regulatory mechanisms responsible for coupling the beat-to-beat variability between signals rather than merely the variability that is elicited [22].

For the two multivariate entropy measures mvSE and mvFME, they are both based on the phase space reconstruction theory and are defined as the negative natural logarithm of the conditional probability that $p$-variate time series of length N, having similar patterns for $m$ points composite delay vector $X_i^m$ within a boundary r, will also repeat for m + 1 points composite delay vector $X_i^{m+1}$. The smaller values of mvSE and mvFME in univariate time series correspond to more regular characteristics in the time series structure, while larger values indicate more complex characteristics in their time series structure. As the examples show in Figure 2, the RR interval time series shows an obvious upward tendency after the stair climbing. Thus, the complexity characteristic in the time series reduces, and the time series structure becomes more regular. The same phenomenon exists in the S1 amplitude series, which shows an obvious downward tendency after the stair climbing. Thus, their univariate mvSE and mvFME values significantly decrease (see Table 3). However, the S2 amplitude series is not largely influenced by the stair climbing, and its univariate mvSE and mvFME values do not decrease significantly for all three PCG signal observations, except the mvFME result from ECG and aortic PCG signals. The change trend in RR interval time series is expected since the heart rate of the subject should decrease during the measurement process after the stair climbing, inducing the obvious upward tendency in the RR interval time series. The decrease of the S1 heart sound amplitude in the S1 amplitude series is also expected during the measurement process after the stair climbing. With the cardiovascular stabilization, the cardiac contractility decrease. Precious studies have shown that the relationship between the S1 amplitude and the cardiac contractility is close [42,43]. The decrease of the cardiac contractility after the stair climbing induces the decrease of the S1 heart sound amplitude since the previous study showed that S1 heart sound amplitude change links to the change of cardiac contractility and cardiac reserve [33,34,44]. Unlike S1 amplitude, the S2 amplitude is regarded to be linked with the peripheral vascular resistance and arterial blood pressure rather than the cardiac contractility and cardiac reserve [34,45]. Thus, there is no obvious amplitude change in the S2 amplitude series after the stair climbing measurement.

For bivariate and trivariate time series, composite delay vector $X_i^m$ consists of two or three signal episodes from different time series. For the simulation Gaussian signals, as shown in Equation (21), if $c=1$, the composite delay vector $X_i^m$ will consist of the same univariate time series as:

$$X_i^m=\left[x_{1,i},x_{1,i+{\tau }_1},\cdots ,x_{1,i+\left(m_1-1\right){\tau }_1},x_{1,i},x_{1,i+{\tau }_2},\cdots ,x_{1,i+\left(m_2-1\right){\tau }_2},\cdots ,x_{1,i},x_{1,i+{\tau }_p},\cdots ,x_{1,i+\left(m_p-1\right){\tau }_p}\right]$$

(22)

In this situation, the counted number of composite delay vector $X_j^m$ that meet $d_{i,j}^m\le r$ increases. Thus, the values of the defined $B_i^m\left(r\right)$ and $B^m\left(r\right)$ increase. When the dimensionality of the multivariate delay vector is from $m$ to $m+1$, $B^{m+1}\left(r\right)$ also increases. However, compared with the situation at dimensionality $m$, the increase at dimensionality $m+1$ was divided by the much larger factor of $p\left(N-n\right)-1$, thus the ration of $B^{m+1}\left(r\right)/B^m\left(r\right)$decreases. Therefore, the mvSE defined in Equation (4) increases. The similar fact also exists in the calculation process of mvFME method in spite of the fact that the true distance in the composite delay vectors, rather than the similarity vector number, was used. Thus, the larger values in mvSE and mvFME indicate the more similar or coupled components in the bivariate and trivariate time series since both mvSE and mvFME increase with the increase of coupling degree parameter $c$ as shown in Figure 3 for the simulation Gaussian signals. For the practical cardiovascular time series analysis, whether using mvSE or using mvFME, after stair climbing state output has significantly lower multivariate entropy values than the rest state (see Table 3), suggesting that the coupled relationships between each two of the RR intervals, S1 amplitude and S2 amplitude series, and also among these three time series, declines after the stair climbing. The reason is because the obvious upward or downward tendency exists in the RR interval or S1 amplitude series after the stair climbing, resulting in weak coupled relationships for the bivariate and trivariate time series. From all of the three PCG signals results, we can find that the most statistical differences between the rest and after stair climbing states are reported by the bivariate time series of RR interval and S1 amplitude, further confirming that the decline of the coupled relationship in the multivariate time series induces the decrease of multivariate entropy values.

Unlike mvSE, where the decision rule for composite delay vector similarity is based on Heaviside function, mvFME uses fuzzy function (herein, the exponential function) to redefine the decision rule for composite delay vector similarity. The difference between Heaviside and fuzzy functions is shown in Equations (23) and (24). The rigid membership degree determination in Heaviside function could induce the poor statistical stability in mvSE, which means that the entropy value may have a sudden change when the threshold value $r$ changes slightly. However, fuzzy function exhibits the gentle boundary effect. This phenomenon has been reported in recent research studies [23,24,46,47,48,49]. In this study, as shown in Figure 3, the smaller SD values of the new proposed mvFME method for the 100 Gaussian signal repeats confirm its better statistical stability than the mvSE method. In addition, mvSE could not distinguish the two physiological states (rest and after stair climbing) for the trivariate analysis of RR interval and S1 amplitude and S2 amplitude, whereas the new proposed mvFME method can, again showing the improvement of mvFME and providing a potential solution to understand the multivariate cardiovascular time series between different physiological states:

$$\mathrm{For} \mathrm{Heaviside} \mathrm{function}: \mathrm{Similarity}(d_{i,j}^m,r)=\{\begin{array}{c}1, if d_{i,j}^m\le r \\ 0, if d_{i,j}^m>r\end{array},$$

(23)

$$\mathrm{For} \mathrm{fuzzy} \mathrm{function}: \mathrm{Similarity}(d_{i,j}^m,r)=\exp (-\frac{{(d_{i,j}^m)}^{n_x}}{r}),$$

(24)

where $d_{i,j}^m$ means the distance of two composite delay vectors $X_i^m$ and $X_j^m$, with $r$ as the threshold value and $n_x$ is the vector similarity weight.

In summary, this study demonstrated that compared with mvSE, mvFME has better statistical stability and better discrimination ability for different physiological states, from the analysis of both multivariate simulation and cardiovascular time series. In future experiments, we will include more challenging examples for discriminating the different physiological/pathological states, rather than only the discrimination between rest and stair climbing states in the current study. Moreover, we expect that the newly proposed mvFME will be useful in the practical clinical applications for not only multivariate cardiovascular signals but also other multivariate or multi-channel physiological signals analysis.