CC BY
23
SPECTRAL INDEPENDENT COMPONENT ANALYSIS OF HEART RATE VARIABILITY
A. V. Martynenko, A.S. Antonova, A.M. Yegorenkov V.N. Karazin Kharkiv National University
SUMMARY
Formal mathematical procedure of developed Spectral Independent Component Analysis (SICA) of heart rate variability (HRV) allows to obtain no more than three components forming the rhythmogramm of healthy person. This fact rationally represents physiological hypotheses about regulation systems taking part in forming HRV phenomenon. Applying SlCA to HRV approves itself on timing intervals form 3 to 15 minutes. Breaking the low limit of this timing interval causes essential worsening of quality of ICA components forming the rhythmogramm. Optimal application of SICA for splitting initial registered HRV signal into components is using SICA with five-minute HRV registering protocol.
KEY WORDS: rythmogramm, heart rate variability, independent component analysis
The phenomenon of heart rate variability (HRV) was noticed by physicians long ago: first mention about HRV in European scientific literature attributed to Stephen Hales in 1733. However, wide development of HRV-based diagnostic methods and their applications in medicine dates from the end of last century, because of development of high-precision digital computer electrocardiography (ECG) and standardization of HRV-technology, carried out by the workgroup of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology
5-minute, 15-minute and 24-hour protocols are most often used for researches of heart rate variability. In HRV analysis the most valuable information can be achieved with the application of dimensional-spectral methods: quick transformation by Furrier, autoregressive, and others. Important spectrum descriptions are: its power, and powers of its separate zones (domains). In 5-minute and 15-minute protocols, 3 inherent zones can be differed: VLF
- zone of very low frequencies (0.0033-0.04 Hz), LF - zone of low frequencies (0.04-0.15 Hz), HF - zone of high frequencies (0.25-0.4 The article is devoted to statistical approach that is based on the method of separating registered composite rhythmogram to its independent components, each formed by corresponding regulating unit of an organism. It’s clear that the application of well known methods (such as factor analysis - principal component analysis) that use a linear decorellation of components, for this purpose is limited and almost incompetent, because of essential non-linearity of the phenomenon being considered - the heart rate variability. The appearance of non-linear statistical analysis method - independent component analysis (ICA) [2] and its successful development in solving of composite medical problems (for example, in encephalogram analysis) show it being very perspective for application in HRV analysis.
Hz) [1]. Thus, clear correlation of the power of every spectral domain with its corresponding regulating unit of an organism was determined. Parasympathic link of the vegetative nervous organism regulation influences fast rate variations of a cardiac rhythm (HF); sympathic link of the vegetative nervous organism regulation influences medium rate variations of a cardiac rhythm (LF); humoral link of the neurohumoral organism regulation influences slow rate variations (VLF). At the same time, traditional approach has certain drawbacks, generally related to formal techniques of selecting regulation domains in the spectrum. Essential progress in solution of the current problem can be attained by application of more complex mathematical methods, which allow separating of initial registered rhythmogram into independent components. Two non-formal approach in selecting spectral-regulating domains in HRV can be noted marked:
- statistical methods of analysis, based on common approaches in time series analysis;
- mathematical modeling, representing physiological features of regulating processes of heart activity.
MATERIALS AND METHODS
To determine ICA application for HRV more strictly, we’ll use statistical “hidden variables” of models [3, 4]. Let’s consider that we have n linear mixes of x\,...,xn independent components (directing influence of each regulating unit).
for each meaning ofj (1)
Timing index t is omitted here. In the ICA model, instead of corresponding timing reckon, both mix Xj(t) and independent component sk supposed to be a random value. Observable values of xj(t) are a sample of this random value (not regarding that the scalar function HR used
in HRV as an initial observable value, due to well known F.Takens theorem [5]). We can suppose without loss of generality, that both mix variables and independent components average value is zero. If that is not so, then the observable variable x, always can be centred by subtraction its sample average of it, that forward us to a model with average value of zero.
Instead of using sums in the previous equation, we can use matrix-vector form of designation. Let x represent random vector which components are xi,...,x„ mixes, and the similar way s - random vector which
components are s1,...,s„. Matrix A is a matrix with elements aij. Bold small letters usually indicate vectors, and bold capital symbols
identify matrixes. All vectors represent column vectors, thus, xT (transponsed x) is a string
vector. With the application of this vector-
matrix nomenclature, the above-mentioned mixing model can be written like this:
x = As (2)
Sometimes we might need the columns of the matrix A, marked as aP The model also can be written like this:
n
x = X as (3)
i=1
Statistical model in equation (3) called the independent component analysis, or the ICA model. The model is an ICA-originative model, that means it describes the way observable data are generated by the mixing process si. Independent components are hidden variables, that means they cannot be observed directly. Also it is accepted that the mixing matrix value is unknown. The only thing we observe is the random vector x, and we have to estimate A and s using it. That must be done in suppositions as common as it is possible.
ICA initial point - the simpliest supposition,
The transformation of the stationary process s(t) in the form of
h (t) = J eat j (l )Fd1 (8)
will be the linear transformation with the spectral characteristic j e L, setting the stationary process n(t) with the spectral measure
G(A) = J |j(l)|2 Fdl . (9)
A
If the initial stationary process s(t) has the spectral density f(X), i.e. there is absolutely continuous spectral measure
that components s, are statistically independent. We also must suppose, that the independent component must have non-Gaussian distribution. Besides, in the initial model we not suppose these distributions to be known (if they are known - the problem can be substantially simplified). For simplicity purpose, we also suppose, that the unknown mixing matrix is a square matrix, but this condition can be widened sometimes. Then, after estimation of matrix A, we can calculate its inverse matrix, marked W, and to obtain an independent component just like this:
s = Wx (4)
We note that for every random stationary process s(t) with a continual covariant function the representation by the way of [6] is allowed:
s(t) = J eatFdl; (5)
where O is a generalized measure on the real axis -oo<A,<oo with values in Hilbert space, O must be such that the ortogonality condition must be carried out for non-overlapping measurable sets:
(F(Ai), F(A 2)) = M{F(Ai)F(A 2)} = 0; (6)
where M is the average of distribution. The orthogonal measure O is specified at o-algebra of the sets, that are measurable relating to nonnegative bounded spectral measure F(A)= M |®(A) |2 ■ and O is concerned with the covariant function B(t) by the correlation:
B(t) =J eatFd1. (7)
F (A) = J f (l )d1, (10)
A
then the corresponding process n(t) has the spectral density
2
g(i)=j(i) f (i). (ii)
The stationary process n(t) allows the spectral representation
h (t) = J eat Ydl, (12)
where the stochastic spectral measure Y concerns to measure ® this way:
Y (A) = J j (l )Fd1 (13)
A
Next, the stationary process s(t) can be obtained with the reverse transformation
s(t) = J eaty (l )Yd1; (14)
where the spectral characteristic of the transformation y(X) is set by formula
y (l) = 1/ j (l) (15)
This way, the linear transformations, being carried out by ICA algorithm, are tottaly equivalent on timing and frequency intervals, i.e. we can transfer ICA process to the spectral area of representation of the required s(t), without loss of commonality. Also, we note important features of the transfer of ICA to the frequency interval, determining advantages of the ICA spectral algorithm (SICA):
1. In case of the ergodic process s(t) spectral moments have the feature of independkelncy at spaces of any measurable sets A1 c and Aj £ L^2’1^ [7]:
M{kJ )(A1 xA 2) = M ^(AjM(^(A 2), (16)
where for the pair of processes s(t) and q(t), concerned with the linear transformation, the correlation exists:
m/^a, xa,)=JJ jik m, M^dUm (17)
A1A2
Next, at the initial timing interval the correlation moments do not have such feature, and the condition of independency is one of the ICA principles, which holds true by completing the minimization process. So, the suggested spectral algorithm is more accurate and has an increased convergence, in comparison with the classic ICA scheme for ergodic process analysis, to which the HRV belongs, too.
2. The covariation at the spectral interval is a power measure of independency of the components, because in the spectral area it always gives us estimation from above. At the timing interval we cannot apply such measure of independency of the ICA components [2].
Hence, in the subsequent text we imply that ICA process shifted from time domain to frequencies domain. The ICA method in frequencies domain we will call Spectral ICA (SICA). Where it is important to separate ICA
and SICA we will take notice.
Intuitively speaking, the key to the estimation of ICA models is the non-Gaussivity. Without the non-Gaussivity, the computation is almost totally impossible. In most applications of the classic statistical theory it is supposed, that random values have Gaussian distribution, so this way it hinders the ICA application. The central limit theorem, the classical result of probability theory, claims that with known initial conditions, the distribution of the sum of independent random values aims the Gaussian distribution. Hence the sum of two independent random values usually has distribution more close to Gaussian distribution than any of both initial random values.Now let’s suppose, that the data vector x is distributed in correspondence with the ICA data model in the equation 4, so that it is a mix of independent components.or simplicity, let’s suppose that all of the independent components have identical distributions. To estimate one of the independent components, we consider the linear combination xi,
y = wTx = £ w x marking i ,
where w is a vector that we must determine. If w would be equal one of the strings of the inversed matrix A, then, in fact, this linear combination equals one of the independent components. Now, here is the question: how we should use the central limit theorem to determine w that way, so it equals one of the inversed A strings? In fact, we cannot determine this w precisely, because we know nothing of matrix A, but we are able to find an estimation function, that gives us a good approach.
To find out, how it leads us to the basic principle of the ICA estimation, let’s replace the
variables, by claiming z = A w . Then we have
y = wTx = wTAs = zTs . Thus, y is the linear combination of si with the data weights zi. The sum of two independent random values is more
T
Gaussian than initial variables, that’s why z s is more Gaussian than any of the si, and it becomes minimally Gaussian, when it actually equals one of the st. In this case, obviously, only one of the zi elements is non-zero (attention: we accept here, that si have identical distributions).
Consequently, we can consider w as a vector, that maximizes the non-Gaussivity of
T
w x. Such vector must certainly correspond (in measures of transformed coordinate system) to vector z, that has only one non-zero
T ___ T
component. It means, w x _ z s equals one of the independent components!
This way, non-Gaussivity maximization of
T
w x gives us one of the independent components. In fact, the optimization relief for non-Gaussivity in the «-dimensional space of w vectors has 2n local maximums, two for every independent component, in correspondence with s, and -Si (confirming that independent components can be precisely estimated only up to a multiplicative sign). To find out a number of independent components, we have to find all of this local maximums. That is not very difficult, because different independent components are non-correlating: we can always look for them in the space only, this way estimations we achieve, are non-correlating with any of the previous. This corresponds to an ortogonalization in the accordingly transformed space.
Mathematical method of Spectral ICA andalgorithm described above were developed based on a standard fast ICA software kit for MATLAB-compatible system [2]. For a dimensional-spectral of the HRV analysis, the quick Furiet transformation were used. We used 24-hour rhythmogramm samples of health volunteers of both male and female gender in the research. We analysed fragments of the rhythmogramm of an arbitrary length and an arbitrary placement into a monitor record. The aim of the research was:
- to determine the adequacy of mathematical method (Spectral ICA) in application to heart rate variability analysis, due to physiological hypotheses about regulating systems, forming HRV phenomenon;
- to determine limits of the ICA application in HRV analysis;
- to find out optimum conditions for using the Spectral ICA in HRV analysis.
RESULTS AND DISCUSSIONS
An example of 8-minute rhythmogramm shown at the figure 1. At figures 2-4, components selected from the initial
rhythmogramm (a), and their spectrums (b) are shown in pairs. Specific of method is a strict satisfaction to the Pirson decorrelation of separated components, this, however, does not provides the decorrelation of their spectrums. Figures 5 and 6 show timing dependencies (heart beat quantities) of correlation coefficients (corr) between the spectrums of separated components of the initial rhythmogramm. Well shown, that the application of ICA in HRV analysis is limited by timing interval from 3 to 15 minutes, with the selectable optimum about 5-7 minutes. Table 1 shows correlation coefficients for the case of separating the initial rhythmogramm into 4 components. It can be found out easily, that adding of the fourth component is not necessary: the offered method determines a high correlation between 3rd and 4th signal, that means impossibility of their separation (fig. 4) (tabl. 1).
Thereby, Spectral ICA applied to HRV analysis used to determine the quantity of the independent variables, in which rhythmogramm can be decomposed by natural appearance. If the quantity of the variables was equal to expected number of regulating systems of an organism, that is, three, then every spectrum of each of these signals appears to be at the only of the regions: low-frequency, medium-frequency or high-frequency one. Also, if there was an attempt to use four or more independent components, then it turned out, that at least two of them appear to be at the same frequency region, and moreover, their spectrums are sufficiently correlated (fig. 4).
Similar results were obtained in all cases, when we use decomposition more then in 3 independent components, - only 3 signals appear to be independent indeed (that is, insignificantly correlative in the spectral area). This sustains the abscence of any other regulating systems of an organism, and also competence of the offered method of analysis.
Table 1
Spectrum domains pairwise correlations
Heart beats Pairwise correlations (corr)
400 0.027 0.038 0.084 0.156 0.285 0.842
450 0.002 0.046 0.053 0.118 0.278 0.9
500 0.004 0.023 0.1 0.147 0.226 0.869
550 0.061 0.075 0.078 0.14 0.167 0.985
600 0.018 0.044 0.121 0.163 0.174 0.998
650 0.036 0.057 0.15 0.136 0.164 0.994
700 0.04 0.043 0.047 0.254 0.266 0.983
750 0.012 0.02 0.075 0.315 0.324 0.997
800 0.061 0.078 0.234 0.234 0.356 0.996
850 0.016 0.037 0.176 0.186 0.224 0.997
900 0.036 0.063 0.079 0.152 0.254 0.979
P,ms2
Pv^.f=650 ms
PTF=1077 ms2
PHF=599 ms
0 100 200 300 400 500 600 700 HB
Fig. 1a. The initial rhythmogramm (S-minute sample)
I
0 0.1 0.2 0.3 0.4 H, Hz
Fig. 1b. The spectrum of the initial rhythmogramm
P, ms2 •
0 100 200 300 400
500 600 ?00eart beat
0 0.1 0.2 0.3 0.4 H, Hz
Fig. 2a. The rhythmogramm, separated in two components Fig. 2b. The spectrum of the components (correlation
coefficient of the components corr=0.1)
Very Low P=909 ms2
P=966 ms2
Hiah P=501 ms2
0 100 200 300 400 500 Heart beat
Fig. 3a. The rhythmogramm, separated in three components
0 0.1 0.2 0.3 0.4 H, Hz
Fig. 3b. The spectrum of the components
(corr=0.32;0.24;-0.04)
VervLow 1
0.1 0.3 0.5 0.1 0.3 H, Hz
0.6
0.4
0.2 -
0.0
7
3
0.3
0.1
0.3
15
0.1
10 -
P
ms
5
0.3
0 -
0.1
0.3
0.1
Fig. 4a. The rhythmogramm, separated in four corHeart beat
Fig. 3b. The spectrum of the components (co^7=0.99;0.31;0.3;0.28;0.26;0.06)
Important questions, that were considered in the work, are the area of application of ICA in HRV analysis and the problem of the optimal rhythmogramm length. To discover this, calculations on a different duration intervals (from 150 to 1350 heart beats) were carried out (total time 15 minutes). Thus, it turned out that when the length of the ICA processed interval is less then 400 points, the low-frequency decomposition component represents
Fig. 5. The correlation between spectrums for two components
separated by SICA(solid), ICA(dashed).
The SICA advantages over ICA have shown at the figures 5, 6, 7 and 8. At the figures 5 and 6 we can see that SICA does not have any data length using limitation unlike ICA. At the figures 7 and 8 presents that ICA cannot separate initial RR (fig.la) for LF and HF signals (the largest corr=0.88) but SICA made separation clearly (the largest corr=0.32). Note, that in the case of ICA we have convergence after 9 iterations, in case of SICA - 6 iterations.
incorrectly due to natural loss of low frequencies. At the same time, increasing the length of the processed interval over than 900 points, in fact do not provided any kind of additional information, but simultaneously the correlation of obtained by ICA signals was growing sharply. That’s why the use of 5-minute protocol for Spectral ICA in HRV analysis is optimal.
Fig. 6. The correlation between spectrums for three components separated by SICA
L.
0.1
0.2 0.3
0.4
H, Hz
Fig. 7. The spectrum of the ICA components
(corr=0.88;0.25;0.21)
Fig. 8. The spectrum of the SICA components (corr=0.32;0.24;-0.041)
REFERENCES
1.
Яблучанский Н.И., Мартыненко А.В., Исаева А.С.. Основы практического применения неинвазивной технологии исследования регуляторных систем человека. - Харьков: «Основа». - 2000. - 88 c.
2. Aapo Hyvarinen and Erkki Oja, Helsinki University of Technology, Laboratory of Computer and Information Science, Independent Component Analysis, P.O. Box 5400, FIN-02015 Espoo, Finland,. - 1999.
3. Jutten C. and Herault J. // Signal Processing. - 1991. - Vol. 24. - P. 1-10,
4. Comon P. // SignalProcessing. - 1994. - Vol. 36. - P. 287-314,
5. Takens F // Lecture Notes in Mathematics. - 1981. - Vol. 898. Springer-Verlag. - P. 366-381.
6. Прохоров Ю.В., Розанов Ю.А. Теория вероятности. - М: Наука. - 1967.- 496 с.
7. Ширяев А.М. // Теор. вероятн. и ее применен. - 1960. - Vol. 3. - С. 293-313.
0.3 -
0.1 -
Р
lF
2
ms
Р
ms
HF
0.3 -
0.1 -
0
