SYSTEM FOR THE ANALYSIS OF THE DAILY HEART RHYTHM AUTONOMIC NERVOUS SYSTEM BALANCE

Abstract

An analysis system for analysis of the balance of circadian heart rhythm autonomic nervous system, the analysis system having tools which are designed to perform the following steps after artifact and arrhythmias removal from RR interval series: identification of the chaotic part of the RR interval series and analysis of the chaotic part; separation of the chaotic part from RR interval series and acquisition of the clean circadian RR series; acquisition of the normalized RR interval series using—interpolation, resampling and normalization; analysis of the heart rate bimodal distribution and modes of the normalized RR interval series; identification of the heart rate circadian period in the normalized RR data series; analysis of the heart rhythm variable part SNS and PNS regulation indicators in the normalized RR data series; generation of the report on obtained results in graph and tabular form.

Claims

1. An analysis system for analysis of the balance of circadian heart rhythm autonomic nervous system comprising a RR interval collection block (20), tools (110) for artifact and arrhythmias removal from RR interval series (65), an identifier for circadian period (360) and RR interval series (180) data analyzer finding sympathetic and parasympathetic nervous system balance, characterized in that the analysis system has tools which are designed to perform the following steps after artifact and arrhythmias removal from (110) RR interval series (65): identification of the chaotic part of the RR interval series (120) and analysis of the chaotic part (190); separation of the chaotic part from RR interval (140) series and acquisition of the clean circadian RR series (150); acquisition of the normalized RR interval series (180) using—interpolation, resampling and normalization (160); analysis (200) of the heart rate bimodal distribution and modes of the normalized RR interval series (180); identification of the heart rate circadian period (300) in the normalized RR data series (180); analysis (400) of the heart rhythm variable part SNS and PNS regulation indicators in the normalized RR data series (180); generation of the report (70) on obtained results in graph (92) and tabular (94) form.

2. The analysis system according to claim 1, characterized in that tools (110) for artifact and arrhythmias removal from RR series (65) performing chaotic series separation from RR series (140) and obtaining clean circadian RR interval series (150) use wavelet transform method.

3. The analysis system according to claim 1, characterized in that to obtain normalized data series (180) interpolation of circadian RR interval series (150) is performed using cubic Hermite interpolation.

4. The analysis system according to claim 1, characterized in that the following steps are performed by the heart rhythm bimodal distribution and analysis of the modes detected in the normalized RR data series (180) block (200): calculation of the input histogram sequence (210) out of normalized RR data series (180) and histogram graphic presentation (histGraf, 270); clustering (220) of the normalized RR data series (180) into two RR clusters (modes) of which the RR cluster 1 (221) (first mode) involves low heart rhythm frequencies and RR cluster 2 (222) (second mode) involves high heart rhythm frequencies; calculation of the statistical characteristics of the first mode series (230) and second mode series (240); calculation of the probability density functions for both modes and their histograms (250); summation of the probability density histograms (260) and their graphic presentation (histGraf, 270); determination of the reliability of the compliance of the distributions between summary histogram (260) and input histogram (210) using cji-square criterion and calculation of the summary histogram for both modes (290); determination of the divergence coefficient in the summary histogram (290) of both modes; determination of the skewness coefficient of the heart rate and skewness coefficient of the total time in the input histogram (290); tabular presentation of the statistical calculations (LaikPask_Lent, 296); graphic presentation of circadian time distribution of the RR interval within both modes (histLaikGraf, 297); graphic presentation of the normalized time course of the circadian tonic sympathovagal balance ((laikGraf, 298); tabular presentation of the values of results of summary calculations of the values of RR interval analysis (laikLent2, 299).

5. The analysis system according to claim 4, characterized in that heart rhythm bimodal distribution and analysis of the modes detected in normalized RR data series (180) block (200) is clustering (220) normalized RR data series (180) into two groups using k-means clustering method.

6. The analysis system according to claim 1, characterized in that block for heart rhythm circadian period identification within RR data series (180) performs: isolation of the heart rhythm circadian period performing continual empirical mode decomposition (Huang transform) involving detection of interim modes (IMF) (320) and their gradual elimination from every residual RR interval series; obtaining of the heart rhythm circadian period series (cirkRR, 360) and analysis of the circadian period structure (370).

7. The analysis system according to claim 6, characterized in that block (300) for heart rhythm circadian period identification within RR data series (180) performs analysis of the circadian period structure using Cosinor calculation method and includes statistical evaluation of parameters obtained presenting them in the tabular (cirkLent, 380) and graph form (cirkGraf, 390).

8. The analysis system according to claim 1, characterized in that block for the analysis of the heart rhythm variable part SNS and PNS regulation criteria (400) performs the following steps: separation of the RR variable part which is achieved by subtraction of the heart rhythm circadian period (cirkRR, 365) from normalized RR data sequence (normRR, 180); filtering of the variable heart rhythm series (420) when three RR interval time series of different frequencies are separated; calculation of the dispersion and mean standard deviations and of three RR intervals time series (430); evaluation of the heart rhythm variability (450); graphic (laikSpekGraf, 435) and tabular (laikSpekLent, 455) presentation of the evaluation of the heart rhythm variability.

9. The analysis system according to claim 8, characterized in that variable heart rhythm series filter (420) employs infinite impulse response three-filter restriction method determining very low frequency (VLF), low frequency (LF) and high frequency (HF) RR interval time series.

10. The analysis system according to claim 8, characterized in that block for the analysis of the heart rhythm variable part SNS and PNS regulation criteria (400) performs the following evaluating of heart rhythm variability (450): for the very low frequency filter (VLF, from 0 Hz to 0.04 Hz) determining passband edge frequency equal to 0.035 Hz and stopband edge frequency equal to 0.04 Hz; for low frequency band filter (LF 0.04-0.15 Hz) determining stopband edge frequency at low frequency side equal to 0.035 Hz, and at high frequency side equal to 0.0155 Hz, passband edge frequency at low frequency side equal to 0.04 Hz, and at high frequency side equal to 0.15 Hz; for the high frequency band filter (HF, 0.15-0.40 Hz) determining stopband edge frequency at low frequency side equal to 0.145 Hz, and at high frequency side equal to 0.45 Hz, passband edge frequency at low frequency side equal to 0.15 Hz, and at high frequency side equal to 0.40 Hz.

11. The analysis system according to claim 10, characterized in that block for the analysis of the heart rhythm variable part SNS and PNS regulation criteria (400) evaluation of the heart rhythm variability (450) is performed by calculation of the dispersion ratio for every time frame within low frequency and high frequency bands presenting phase balance of SNS and PNS changes.

Description

SHORT DESCRIPTION OF THE FIGURES

[0057] Best options for the implementation of the invention are detailed below together with references to the attached diagrams involving:

[0058] FIG. 1: System for the analysis of the circadian heart rhythm autonomic nervous system balance.

[0059] FIG. 2: Flow chart of the analysis system block for heart rhythm artifacts and arrhythmias removal, chaotic part detection and series interpolation.

[0060] FIG. 3: Flow chart of the block for the heart rhythm bimodal distribution and analysis of the detected modes.

[0061] FIG. 4: Flow chart of the block for the identification of the circadian period of the heart rhythm.

[0062] FIG. 5: Block for the analysis of the SNS and PNS regulation criteria of the heart rhythm variable part.

[0063] FIG. 6: Content of the data presented by the reporting block.

[0064] FIG. 7: RR interval series reporting graph.

[0065] FIG. 8: Histogram of the block according to FIG. 2 chaotic series.

[0066] FIG. 9: Table containing data obtained according to FIG. 2 and FIG. 3 blocks.

[0067] FIG. 10: Image of the histogram of the sequence according to FIG. 3.

[0068] FIG. 11: Histogram of the RR interval bimodal distribution according to FIG. 3 block.

[0069] FIG. 12: Statistical estimation of the RR interval bimodal time distribution every 2 hours.

[0070] FIG. 13: Time course of the circadian tonic sympatho-vagal balance according to FIG. 3.

[0071] FIG. 14: Statistical evaluation of the structure of circadian period of the circadian heart rhythm.

[0072] FIG. 15: Diagram of the circadian period sequence of the circadian heart rhythm.

[0073] FIG. 16: Classical infinite impulse response three-filter system.

[0074] FIG. 17: Circadian time course of the indicators of heart rhythm variability.

DETAILED DESCRIPTION OF THE INVENTION (DESCRIPTION OF THE PREFERRED EMBODIMENT)

[0075] Subject's 10 physiological signals enter heart rhythm RR interval collection block 20 (FIG. 1) where digital values of signals are placed with at least at 512 Hz discretization frequency. Block 20 is able to collect only signal values or employ software that could be able to produce durations of the heart action periods. Signal transmission from the block 20 to the data formatter 30 is performed after completion of daily data collection adding starting time and completion time (timer). Signals in the format of the selected data from the data formatter 30 though local network sending equipment 40 are transmitted to various computer network receiving equipment 50 which enables their storage in the computer 60 memory 65. Subsequently computer processes data received and transmit them to the report block 70 and printer.

[0076] The order of the processing of data downloaded to the memory 65 is detailed below. RR interval series file from the memory 65 is transferred directly to the report block 70 in the form of RR interval report graph RR-Graf, 66 (FIG. 6, FIG. 7) and enters block for heart rhythm purification. chaotic part identification and series interpolation 100 (FIG. 1, FIG. 2) where artifacts and arrhythmias are removed from RR series (100), RR series chaotic parts are identified (120) and chaotic series chaosRR, 130 is identified and separated from RR series (140); interpolation, resample and normalization (160) and analysis of chaotic parts (190) are also performed resulting in chaotic series histogram chaos_HistGraf, 195 (FIG. 8) and completion of time distribution table laik_Lent 1, 199 (FIG. 9).

[0077] Artifact and arrhythmias removal (110) from circadian RR interval series (65) is executed using wavelet transform method (Stephane Mallat. A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way. Academic press of Elsevier, 805p., 2009, ISBN-13: 978-0123743701) which separates slow fluctuations from RR interval series and store them in the memory (65). Besides slow fluctuations series mean (LxRR) of frame and standard deviation (LsRR) are calculated in every RR interval series with moving time frame (L≧300 s) while moving through single memory address. RR interval values overriding LxRR+3*LsRR boundaries are adjusted at every step by replacing them with new LxRR values. After completion of movement slow fluctuation RR interval series obtained after wavelet transform is added to the new adjusted RR interval series and clean circadian RR interval series is obtained (den RR, 150) which is stored in the memory 65.

[0078] Obtained RR interval series (denRR, 150) using obtained real time series (timer) is interpolated by cubic Hermite interpolation method (described Burden, Richard L.; Faires, J. Douglas (2004). Numerical Analysis. Belmont: Brooks/Cole. P. 872. 2011) with the period equal or shorter than 0.5 second. Mean (xRR) and standard deviation (sRR) of the obtained new RR data sequence is calculated and stored in the memory 65 for future calculations. This series is normalized (average is subtracted and result is divided by two) and final RR data series (norm RR) ready for analysis is obtained.

[0079] Block for the heart rhythm bimodal distribution and analysis of the detected modes 200 (FIG. 3) executes calculation of the initial histogram series (210) from normalized RR data series (norm RR, 180) with 60-80 graduation marks, their number is determined by selected histogram interval delta with values between 12 and 15 ms. Results are presented in graph form (histGraf, 270) (FIG. 10). Using k-means clustering method (Spaeth, Helmuth. Cluster-Analyse-Algorithmen zur Objektklassifizierung and Datenreduktion, Verfahren der Datenverarbeitung, Muenchen: R.Oldenbourg, P. 71, 1975) RR series clustering is carried out when normRR series (180) is divided into two clusters (modes): first mode is RR cluster 1 (221) with low heart rate frequency, and second mode is RR cluster 2 (222) with high heart rate frequency. For this mode (230) No 1 mean xRR1, standard deviation sdRR1 and number of RRs are calculated. Using the same interval delta out of xRR1 and sdRR1 probability density function PDF1 is calculated which is normalized with respect to the whole normRR series length (No 0) (240) using normalization coefficient NO 1/No 0 (240). For the second mode these calculations are repeated, however, using their mean xRR2, standard deviation sdRR2 and normalization coefficient NO 2/No1 (250). Both density functions are analyzed in the form of histogram; they are summarized (260) and presented as the graph together with histogram graph (histGraf, 270) (FIG. 10). Using chi-square criteria distribution consistency between aggregated and initial (210) histogram is determined and calculation of the aggregate histogram of both modes is performed (290). Maximum values of peaks equal to xRR1 and xRR2 are found at aggregated histogram as well as their equivalence in the graduation marks of the RR value histogram. Difference between xRR1 and xRR2 is equal to the coefficient of divergence between PNS and SNS (FIG. 10) which shows level of the dominance of both parts of ANS in the heart rhythm tonic regulation. Level point of intersection of probability densities of xRR1 and xRR2 values is determined in the aggregate histogram as well as its equivalent to the graduation mark of the RR histogram (FIG. 10), which divides initial normRR sequence histogram histGraf (270) into two parts where numbers of RR intervals nrRR1 and nrRR2 are counted as well as their summarized durations sekRR1 and sekRR2. Ratio of nrRR1 and nrRR2 is calculated which is defined as heart rate coefficient of asymmetry, and calculated ratio of sekRR1 and sekRR2 represent total duration asymmetry coefficient. Point of intersection of both distributions (level point) shows value of the RR interval dividing level of dominance of both parts of ANS (SMS and PNS) or value at which impact of both parts of ANS is equally balanced.

[0080] Summary duration times with desirable time frames (e.g., L≧1 hour) and RR interval histogram distribution per day presented in the histogram histLaikGraf, 297 (FIG. 11) showing ANS heart rhythm asymmetry 24-hour time course are calculated from the RR intervals of every mode. Statistical calculations of RR intervals are performed at every time frame, and results are presented in the table LaikPaisk_Lent, 296 (FIG. 12).

[0081] After figure of one was assigned to any value of RR of the first mode and figure of two was assigned to any value of RR of the second mode they are summarized in the new value frame (I≧300 s), result is multiplied by two, divided by the length of frame and after figure of three is subtracted from the result normalized sympatho-vagal balance time course laikGraf, 298 (FIG. 13) is obtained. When value equals to plus one (+1) is obtained, SNS is totally dominating in the heart rhythm regulation. Values equal to minus one (−1) show that PNS is totally dominating (FIG. 13). Additionally RR interval values in every time frame are also presented in the general table laikLent2, 299 (FIG. 9) as summary calculation.

[0082] In the block for the identification of the circadian period of the heart rhythm (FIG. 4) heart rhythm circadian periods are extracted from normalized RR data series (normRR, 180) using continual empirical mode decomposition (Huang transform) (Huang N. E., Shen, Z., Long, S. et al., The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-Stationary Time Series Analysis. Proc. Royal Soc. London Series A, 1998, vol. 454, pp. 903-995; Ming-Chya Wu and Chin-Kun Hu. Application of Empirical Mode Decomposition to Cardiorespiratory Synchronization. S. K. Dana et al. (eds.), Complex Dynamics in Physiological Systems: From Heart to Brain, Understanding Complex Systems, Springer Science+Business Media B. V. 2009). Intrinsic mode functions (IMF) are identified in the RR data series (normRR, 180) and they are gradually removed from every occurred residual interval series. When less than three peaks and less than four valleys remain in the residual series operation is terminated. Residual series is the heart rhythm circadian period series (cirkRR). By applying empirical mode decomposition (EMD) (Huang transform) method to RR data series (180) its two components are extracted: RR values tonic circadian HR period and residual fluctuations of RR values that are shorter than circadian period.

[0083] The following stages and steps are performed applying empirical mode decomposition: [0084] 1. By applying first derivative local extreme points (305) are identified in the whole RR data (normRR, 180) time line x(t). [0085] 2. All maxima are separately connected (310) by cubic splines using cubic Hermite polynomial interpolation (Burden, Richard L.; Faires, J. Douglas (2004). Numerical Analysis. Belmont: Brooks/Cole. 872 p. 2011) resulting in the formation of upper envelope curve, u(t).

[0086] 3. All local minima are connected (315) applying Hermite interpolation procedure for local minima resulting in the formation of lower envelope curve I(t). [0087] 4. Data mean for both envelope curves is calculated using operation m(t)=[u(t)+I(t)/2] and interim intrinsic mode function IMF(t), h(t) (320) is determined by calculating the difference between data x(t) and obtained envelope values of mean m(t), h(t)=x(t)−m(t). IMF (1) is subtracted from the series normRR, 180 for the first time. [0088] 5. INF(i) series, h(t) dispersion and ratio with previously calculated IMF(i−1) dispersion Var(i−1) are calculated in order to achieve reliable accuracy of calculation with the probability of 0.01 (Fisher's dispersion distribution coefficient). [0089] 6. Steps 1-5 are repeated with interim intrinsic mode function IMF(i), h(t) for as long as its magnitude of dispersion compared (330) with previously obtained IMF(i−1) match the criteria (331) (more than 0.01). [0090] 7. When restriction (332) indicated in the 6.sup.th step is matched after previous steps, this stage and extraction of the interim IMF(i) is completed; this interim IMF(i) taken as a last IMF(i) component, c(t) found at this stage. [0091] 8. Next stage repeats steps 1-7 for remaining series r(t): r(t)=x(t)−c(t) (340), where r(t) is considered as an new time series, x(t). [0092] 9. Number of peaks and valleys in the new time series x(t) is determined.

[0093] Stages repeat steps 1-9 every time checking number of peaks and valleys (355) as long as restriction (355) is matches, and these steps are completed when residual time series contains less than three peaks and less than four valleys (357), and thus heart rhythm circadian period series (cirkRR, 360) is obtained which may have close to sinusoid wave shape. To obtain heart rhythm circadian period sinusoid-shaped curve multicomponent Cosinor calculation method is applied (Nelson W, Tong Y L, Lee J K, Halberg F. Methods for cosinor-rhythmometry. Chronobiologia. 1979, 6(4), 305-323; Cornelissen G. Cosinor-based rhythmometry. Theor Biol Med Model. 2014, 11; 11(1):16.; Bingham C, Arbogast B, Cornelissen Guillaume G, Lee J K, Halberg F: Inferential statistical methods for estimating and comparing cosinor parameters. Chronobiologia 1982, 9:397-439). Cosinor method allows obtaining heart rhythm circadian period curve and performs its structure analysis (370) extracting additional elements of smaller period. Cosinor method allows calculating structure indicators of the circadian heart rhythm circadian period: mean amplitudes of the circadian period and its components (MESOR), waves amplitudes, acrophases, their duration and dispersion. Statistical evaluation of obtained parameters is performed using Fisher's reliability criterion (F-test), and results are presented in the table cirkLent, 380 (FIG. 6) and also presented in graph form cirkGraf, 390 (FIG. 15).

[0094] In the block for the analysis of the SNS and PNS regulation criteria of the heart rhythm variable part (FIG. 5) first RR variable part separation is performed (405), which is obtained by subtraction of heart rhythm circadian period and its components series (cirkRR, 365) (FIG. 3) from normalized RR data series (normRR, 180) (FIG. 2). New interim RR data series (flucRR, 410) obtained represents RR interval frequencies that are higher than circadian fluctuations.

[0095] Using classic infinite impulse response three-filter restriction method (Butterworth, Kaiser et al.) first of all filtering of the variable heart rhythm sequence (420) is performed determining three RR interval time series: very low frequency (VLF), low frequency (LF) and high frequency (HF).

[0096] Further by applying three-filter restriction system (FIG. 16) analysis of the variability of fluctuating RR data series is performed: [0097] 1. For the very low frequency filter (VLF, from 0 Hz to 0.04 Hz) determining passband edge frequency equal to 0.035 Hz and stopband edge frequency equal to 0.04 Hz. [0098] 2. For low frequency band filter (LF 0.04-0.15 Hz) determining stopband edge frequency at low frequency side equal to 0.035 Hz, and at high frequency side equal to 0.0155 Hz, passband edge frequency at low frequency side equal to 0.04 Hz, and at high frequency side equal to 0.15 Hz. [0099] 3. For the high frequency band filter (HF, 0.15-0.40 Hz) determining stopband edge frequency at low frequency side equal to 0.145 Hz, and at high frequency side equal to 0.45 Hz, passband edge frequency at low frequency side equal to 0.15 Hz, and at high frequency side equal to 0.40 Hz. [0100] 4. Equal conditions for all filters determining stopband percent ripple equal or larger than 60 dB and passband percent ripple equal and lower than 0.01 dB. [0101] 5. Specified filter characteristics allow obtaining constant result regardless investigator's choice of available classic infinitive response filters for the filtering of circadian heart rhythm except filters for short time series because all these filters possess different distortions of the results initial sections of series. [0102] 6. In order to avoid these distortions every defined time series is filtered twice, first direct filtering is performed, and second filtering is performed by turning around memory addressing of the reversely obtained series and by turning around filtering results once more. This helps to avoid violation of the structure of initial time in the RR interval series. [0103] 7. After filtering every series is multiplied by median standard deviation of the initial RR sequence.

[0104] Dispersion and mean standard deviation values that selectively may be expressed by dispersion values (ms.sup.2) or amplitude values (ms) are calculated in the selected time frames (I≧100 s) of the three series obtained after filtering. To evaluate these dispersions and mean standard deviations their distribution in the desirable time frame within circadian time frame is calculated. Obtained results represent SNS and PNS involvement level in the regulation of the phase changes of heart rhythm throughout 24 hours.

[0105] To evaluate variability of heart rhythm (450) ratio of dispersions found in low frequency (LF) and high frequency (HF) bands is calculated showing changes in SNS and PNS phase balance within 24 hours in every time frame (No Authors Listed. Heart Rate Variability: Standards of Measurement, Physiological Interpretation, and Clinical Use, Circulation, 93, (1996), pp. 1043-1065; Heart Rate Variability: Standards of Measurement, Physiological Interpretation, and Clinical Use, European Heart Journal, 17, Prepared by the Task Force of The European Society of Cardiology and The North American Society of Pacing and Electrophysiology; published by the American Heart Association, Inc.; European Society of Cardiology, (1996), pp. 354-381).

[0106] Time course of the heat rhythm variability is further presented in the graph (laikSpekGraf, 435) (FIG. 17) and table (laikSpek_Lent, 445) (FIG. 12).