|Home | About | Journals | Submit | Contact Us | Français|
This paper presents the concept of perfect matrices of Lagrange differences which are used to analyze relationships between RR and JT intervals during the bicycle ergometry exercise. The concept of the perfect matrix of Lagrange differences, its parameters, the construction of the load function and the corresponding optimization problem, the introduction of internal and external smoothing, embedding of the scalar parameter time series into the phase plane—all these computational techniques allow visualization of complex dynamical processes taking place in the cardiovascular system during the load and the recovery processes. Detailed analysis is performed with one person’s RR and JT records only—but the presented techniques open new possibilities for novel interpretation of the dynamics of the cardiovascular system.
ECG analysis is the basic, the primary, and the most studied noninvasive technique used for the contemporary investigation of the functionality of the cardiovascular system. The relevance of ECG parameters to the clinical practice is unquestionable ECG parameters are effectively used for the identification of various heart rate and conductivity defects, different heart hypertrophies and ischemic processes. Cardiac time intervals are sensitive markers of cardiac dysfunction, even when it goes unrecognized by conventional echocardiography . There exist some opinions that diagnostic features of ECG parameters are completely understood—and that there is no need to seek for any new approaches in ECG analysis.
However, interrelations between EGC parameters is still an active area of research. A typical example is the relationship between RR and QT intervals. This relationship was first described by Bazzet  in a form of a functional correlation and the definition of tolerance intervals for QT. But one hundred years later, discussions about this relationship do not seem to be stopping. A number of new functional relations between these parameters were proposed. For example, functional relations among RR, JT, and QT intervals are evaluated in . Many researchers were proposing new analytical expressions and confidence intervals for inter-parameter relationships—all adapted to their specific cohort [4–6]. Unfortunately, other researchers were not able to prove the validity of these relationships . This is probably due to some of the reasons described below.
For example, it is inappropriate to interpret QT as a characteristic of one physiological process. It is well known that QT is the sum of the depolarization and the repolarization of heart’s myocardium. However, depolarization and repolarization may have completely different variation tendencies. For example, depolarization and repolarization may both become shorter during the load—but sometimes depolarization becomes longer and repolarization becomes shorter (what is a rather common phenomenon during the extreme loads). Thus the sum of the depolarization and the repolarization would not represent a clearly determined physiological process during an exercise with a constantly increasing load—and a simple QT-RR relationship cannot describe the behavior of the human cardiovascular system in general. In other case, the examination of QT and JT intervals under the supervised consumption of dietary acids is presented in . The correlation amongst QT, JT and their alternative counterparts is discussed in . Similar problems exist for other relationships between ECG parameters (for example PQ-RR, etc.).
Relationships between ECG parameters do vary according to various physiological and pathological reasons. Holistic models interpret a human being as a complex system—where nonlinear chaotic processes play an important role in the relationships between different sub-systems and in generating reactions of these sub-systems . Therefore, it is probably illogical to seek a unified deterministic model which could describe the relationships between ECG parameters. It makes sense to observe dynamical processes, dynamical relationships—which could exhibit complex chaotic behavior .
The relationships between ECG parameters do depend on different time scales and other factors. For example, the influence of QRS duration on the JT and QT intervals is assessed in . In addition, it is well known that there exist completely different short-term and long-term adaptations of the cardiovascular system to physical loads . It becomes clear that the width of the observation window may seriously impede the computational results.
The complexity of the relationships, the chaotic nature of the processes, the fractality of time scales—all that yields a necessity to develop such computational techniques which could assess the dynamism of these relationships between RR and JT intervals. As mentioned previously, investigation of relationships between RR and JT intervals continues during the last decades. For example, this relationship is used for the detection of prolonged repolarization in ventricular conduction defects —but the JT interval is constructed as a linear function of the RR interval in this study. Similarly, it is shown in  that heart rate-corrected JT interval is a good estimate of specific repolarization time in a cohort of physically fit university students.
On the contrary, the main objective of this paper is to propose a visualization technique of relationships between RR and JT intervals which could reveal the evolution of complex dynamical processes in the self-organization of the heart system during the load and the recovery processes.
The research met all applicable standards for the ethics of experimentation. Permit to perform biomedical investigation was granted by Kaunas Regional Ethics Committee for Biomedical Investigations, No. BE-2-51, 23.12.2015. ECG bicycle ergometry exercise was used to record cardiac RR and JT intervals. Participants provided written informed consent prior to the experiment.
The assessment of functional ECG parameters was performed by using ECG analysis system “Kaunas-Load” [15–18] developed at the Institute of Cardiology, Lithuanian University of Health Sciences. The second ECG lead is used for the signal processing task. JT interval is evaluated as the time difference between the J point and the end of the T wave. The J point is defined as the moment where depolarization processes expire and repolarization processes begin. The J point is automatically detected by fixing the moment where the speed of ECG electrical processes change just before the repolarization. The end of the T wave is detected at the end of repolarization processes. This moment is identified by computing the intersection between the lead axis and the tangent line to the descending slope of the T wave. Numerous clinical trials were used to assess the measurement accuracy of the JT interval during the development of “Kaunas-Load” system .
The bicycle ergometry system (Fig 1) is used for generating stepwise increasing physical loads—whlist “Kaunas-Load” is used for a synchronous registration of twelve different standard parameters of the ECG. The bicycle ergometry exercise is initiated at 50W load—and the load is increased by 50W every consecutive minute. The patient is asked to maintain a constant 60 revolutions per minute bicycle pedals spinning rate during the whole exercise. The load is increased up to 250W and the exercise is continued until the first clinical indications for the limitation of the load according to AHA (American Hearth Association) are observed. The cohort comprised 10 adult men within the age range of 29 and 36 years. All persons have been regular visitors of a local fitness club in Kaunas, and did not have any registered medical conditions.
Without losing the generality, we select two ECG parameters—RR and JT intervals. Sequences of RR and JT intervals are recorded during the bicycle ergometry experiment and denoted as vectors x = (x1, x2, …, xn) and y = (y1, y2, …, yn) accordingly.
Given a scalar time series, derivatives of variables can be computationally assessed by Lagrange differences at the nodal point of the time series. It is well known that derivatives reconstructed from a scalar time series tend to intensify the noise embedded into that series. However, clever manipulations with derivatives may also help to detect and amplify small changes of the dynamical processes—which may not be visible in the original time series. Keeping this information in mind we use a combination of the scalar values of time series x and y (RR and JT intervals) and cross derivatives between elements of x and y. Our intention is to produce one single scalar attribute which could characterize dynamical relationships between x and y. We use square second order matrices, place different values and different differences as the elements of these matrices, and compute a single parameter representing a certain property of that matrix. The following sections are dedicated to the considerations on the architecture of these matrices and the selection of the parameter representing local relationships described by these matrices.
As mentioned previously, we will consider second order square matrices. We assume that every element of a matrix can be either a single element of the time series x or y—or a difference between elements of time series x and y:
Note that symbols x and y in Eq (1) represent only the time series; indexes representing a particular time moment will be assigned later. Different signs can be assigned to elements of x and y. Also, differences are allowed only between elements of different time series (in order to minimize the straightforward amplification of noise). Second order square matrices comprising such elements are named as matrices of Lagrange differences (zeroth or the first order Lagrange cross differences).
A perfect matrix of Lagrange differences is a second order square matrix whereas elements of that matrix do satisfy the following requirements:
Several examples of not perfect and perfect matrices of Lagrange differences.
A natural question is about the number of different perfect matrices of Lagrange differences. In general we seek to represent six different elements of time series x and y as shown in Table 1.
A graphical visualization of perfect matrices of Lagrange differences would help to classify these matrices and to interpret their structure. Zeroth order differences (±x and ±y) are located on the main diagonal—but indexes of these elements can be different. We draw circles around elements from Table 1 which are selected as zeroth order differences. For example, if ±xn is selected as an element representing a zeroth order difference, then a circle is drawn around the element in the first row and the second column of Table 1.
A first order difference is visualized by drawing an arrow connecting the two elements of the difference. For example, an arrow connecting the element in the first row and the first column and the element in the second row and the first column is used to represent a difference ±(xn−δ − yx−δ) in Table 1.
As noted in Example 1, a matrix is a perfect matrix of Lagrange differences. As previously, symbol δ denotes the time delay; n—the time moment around which the elements in Table 1 are centered around; (1) in the superscript denotes the number of the perfect matrices of Lagrange differences. The graphical representation of this matrix is shown in Fig 2.
Such graphical representation of perfect matrices of Lagrange differences allows an efficient classification of all possible perfect matrices. In fact, there are only 18 types of perfect matrices of Lagrange differences—all of them are illustrated in Fig 3.
One of the main characteristics of any square matrix are its eigenvalues (or the spectrum of the matrix). It is clear that elements on the diagonal and the anti-diagonal can be interchanged and their signs can be switched without affecting the maximal absolute eigenvalue. It means that every graphical representation of a matrix in Fig 2 can be a result of 24 = 16 distinct perfect matrices of Lagrange differences in terms of the maximal absolute eigenvalue. Overall, 18 different representations yield 288 distinct perfect matrices of Lagrange differences.
As mentioned previously, a scalar parameter representing local relationships described by perfect Lagrange matrices must be selected. In general, the selection of such a parameter is an ill-posed inverse problem of parameter identification. Some sort of error function must be defined before different parameters could be assessed and compared in respect of their representativeness.
As mentioned in Section “Methods”, loads of the bicycle ergometry exercise are increased at every consecutive minute. However, the RR intervals do change during the exercise. The number of inter-beat intervals recorder during the first minute (before the exercise was started) and the further minutes are different (Fig 4). The x-axis in Fig 4 represents the index number of the x-time series. Note that the x-axis does not represent time—each inter-beat interval is different. Cumulative summation of the x-variable allows simple reconstruction of the time variable. However, we do not use the time variable in further computational experiments—all computations will be based on the x-time series data.
It is clear that Fig 4 is constructed for a single individual person. The investigated person did manage to reach 250W bicycle ergometry load—and then (almost immediately after reaching the highest load) the ergometry test was stopped. However, the recordings of the person’s ECG parameters continued during the recovery process.
Therefore, first 96 RR intervals fit into the first minute when the load was 0W; 108 RR intervals fit into the second minute (the first minute of the exercise when the load was 50W); 124 RR intervals into the 3-rd minute (100W load); 140 RR intervals into the 4-rd minute (150W load); 158 RR intervals into the 5-th minute (200W load) and only 5 RR intervals into the final part of the exercise when the load was increased to 250W (Fig 4). After the exercise was terminated, the recording of ECG parameters continued for another 5 minutes (Fig 4).
Moreover, we norm the opposite values of the load to the y-axis of Fig 4 and fit it to the interval [0, 1]. Therefore, 0W load is mapped to 1; 50W—to 0.8; 100W—to 0.6; 150W—to 0.4; 200W—to 0.2 and 250W—to 0 (Fig 4). In other words, we construct the target function—and the parameter representing local relationships described by perfect Lagrange matrices should follow this target function as close as possible.
We select the following parameters: the maximal absolute eigenvalue of a perfect matrix of Lagrange differences max |λ|, the minimal absolute eigenvalue min |λ|, the structural coefficient str = max |λ|/min |λ| and the discriminant dsk = (a11 − a22)2 + 4a12a21 (where indexes denote the location of elements). The time delay δ is set to 1.
Let us denote the scaled inverse values of the load (Fig 4) as lk; k = 1, 2, …, 1199. Now, the values of the particular parameter computed for the perfect matrix of Lagrange differences constructed from the RR (x-time series) and JT (y-time series) are denoted as pk; k = 1, 2, …, 1199 (note that δ = 1). The optimization problem is formulated as follows—minimize RMSE (root mean square error) between lk and pk—by selecting the most appropriate parameter:
where the minimization is performed in respect of a particular expression of the parameter; values of pk are normalized into interval [0, 1] before the optimization process is commenced.
Computational optimization of the best parameter is done by using ECG records of 10 persons. Every person did stop the bicycle ergometry exercise at different moments—so the number N in Eq (2) is different for each person (N = 1199 for the first person—his bicycle ergometry load protocol is shown in Fig 4). However, we do not preselect a single perfect matrix of Lagrange differences—we do average RSME for all 18 perfect matrices—the results are presented in Table 2. Note that the person whose load diagram is shown in Fig 4 is represented as Person #1 in Table 2.
It can be seen that min |λ| is the best parameter for the second person—but max |λ| is the best parameter for all other persons. Therefore, we fix max |λ| as the best parameter representing the local relationships described by perfect Lagrange matrices.
Computational experiments are continued with where denotes a perfect matrix of Lagrange differences centered around time moment k. It is well known that moving averaging (MA)  helps to smooth coarse signals. Moreover, MA is frequently considered as a simple yet effective option for the prediction of irregular time series .
However, before dealing with MA, we introduce the procedure of internal smoothing. The radius of internal smoothing Δ defines how many parameter values computed for perfect matrices of Lagrange differences (at different δ) are averaged. For example, inner averaging (without external averaging) reads:
Now, the smoothing of the reconstructed parameter values at different time moments k is performed using the standard MA—such smoothing procedure is denoted as external smoothing. However, we do use only odd widths of the averaging windows (denoted as m) for MA—the averaging is constructed symmetrically around the time moment k. Let us denote the reconstructed parameter value as pk(s, Δ, m), where k is the time moment; Δ is the radius of internal smoothing; m is the radius of external smoothing. Then
Note that the order of internal and external averaging is not important (Eq (4)). However, that would not be the case if weights for different δ would be introduced for internal averaging (all weights set to be equal in this paper). Note that external averaging (without internal averaging) reads:
and the procedure of internal averaging (without external averaging) yields Eq (3).
Let us consider the following time series: x = (x1, …, x7) = (1, 3, 2, 0, 4, 3, 1); y = (y1, …, y7) = (2, 1, 0, 4, 3, 1, 1) and let us compute p4(1, 2, 1).
A schematic diagram illustrating the computational procedures of internal and external smoothing is presented in Fig 5. First of all, the perfect matrix of Lagrange differences is constructed around k = 4: . The eigenvalues of are: λ1 ≈ −0.4495; λ2 ≈ 4.4495. Thus, .
Next, . The eigenvalues of are: λ1 ≈ −0.8284; λ2 ≈ 4.8284. Thus, .
Since Δ = 2, the procedure of internal smoothing yields: p4(1, 2, 0) ≈ 0.5(4.4495 + 4.8284) = 4.6389.
Since m = 1, computations must be repeated for k = 3 and k = 5. ; λ1,2 ≈ 1 ± 2.6458i; . ; λ1,2 = 1; . Now p3(1, 2, 0) ≈ 0.5(2.8284 + 1) = 1.9142.
; λ1,2 ≈ 3.5 ± 2.7839i; . ; λ1 = 3; λ2 = 4; . Now p5(1, 2, 0) ≈ 0.5(4.4721 + 4) = 4.236.
Finally, p4(1, 2, 1) = ⅓(p3(1, 2, 0) + p4(1, 2, 0) + p5(1, 2, 0)) ≈ 3.5964
Now, it is possible to perform the minimization of RMSE in respect of s and Δ (Table 3). Note that this optimization is performed for the first person only (whose load diagram is presented in Fig 4). In order to lessen the extent of this computational experiment we set m = Δ and vary the internal and external smoothing from 1 to 8; the computation of RMSE is performed for all 18 perfect matrices of Lagrange differences (Table 3).
It can be observed that a moderate internal and external smoothing (m = Δ = 3) helps to minimize RMSE for almost all perfect matrices of Lagrange differences. In fact, the average value of RMSE (averaged for s = 1, 2, …, 18) is minimal at m = Δ = 3 (Table 3). Therefore, we will fix m = Δ = 3 for further computations.
Also, it can be observed that the RMSE value has a lower dependency on s than Δ (Table 3). Apparently, effects induced by the smoothing operation have a greater impact to the RMSE than the parameter s (at a fixed smoothing window). Some minimal smoothing helps to remove the noise and to minimize the RMSE (at Δ = 3). However, wider smoothing windows distort the signal  (the signal loses the variability, the effect of a time delay is introduced).
Table 4 shows that the lowest RMSE value is achieved at m = Δ = 3 and s = 3 (the third perfect matrix of Lagrange differences). The variation of pk(3, 3, 3) (normalized into interval [0, 1]) is illustrated in Fig 6.
Fig 6 shows the variation of the optimal parameter pk, reconstructed using the optimal matrix of Lagrange differences and optimal smoothing. Yet, Fig 6 reveals interesting features. It is well known that the “collapse of complexity” happens with the heartbeat time series at different pathologies [23, 24]. The parameter pk reaches zero at the moment when person #1 cannot continue the bicycle ergometry test at the maximal load. But pk is the modulus of the maximal eigenvalue of the perfect matrix of Lagrange differences. The complexity of the system is minimal when the modulus of the maximal eigenvalue of the matrix of differences becomes minimal. That also explains why the load diagram (Fig 4) represents the scaled inverse values of the actual load.
The variation of the parameter pk in Fig 6 also reveals an interesting phenomenon. The system (represented by pk) tries to stabilize around the scaled normalized constant value of the load—that is particularly clear around the load level 0.8 in Fig 6. This constant value of the load can be interpreted as a stable fixed point—which losses its stability due to the fatigue. Similar loss of stability in a quasi-isometric arm-curl exercise reflected by the increase of low-frequency fluctuations is observed in .
It is well known that phase maps of bioelectrical signals my useful for the investigation of dynamical processes . In order to better represent the attractors (and to assess their stability) we also visualize the dynamical processes of Fig 6 in a phase plane. Time-delay embedding  is used to visualize temporary stabilization and the subsequent loss of the stability of attractors in the phase plane (Fig 7); the optimal time delay τ is determined by using the technique presented in  (τ = 4 RR intervals). The normalized loads at 1, 0.8, 0.6, 0.4, 0.2 and 0 are shown as star type markers in Fig 7.
Temporal stabilization of pk around 0.8 fixed-point attractor is especially well expressed in Fig 7. It is clear that person #1 would be able to continue the bicycle ergometry test for a much longer time at the load equal to 0.8. However, the increased load (to 0.6) forces the system to re-organize (Fig 7). Note that these reorganization processes are interpreted and visualized just by using two time series—RR and JT interval sequences.
The reorganization process to the 0.6 fixed-point attractor is rather complex. Initially, the system tries to loop around a non-existing attractor slightly higher than 0.6—then goes through some transients—and settles around 0.6 (Fig 7). Again, it is clear that person #1 could continue the exercise much longer at the load 0.6.
However, the load is again increased to 0.4. Now, the system quickly converges to the 0.4 fixed-point attractor—and the complexity of the transient processes is much lower compared to the previous transitions (Fig 7). Again, person #1 could continue the exercise for a longer time than he is allowed.
But the convergence of the process to 0.2 fixed-point attractor is much more complicated—in fact the sequence of pk never stabilizes around 0.2 (Fig 7). It seems that person #1 would not be able to continue the load at 0.2 much longer than he already did. The transition to the 0 fixed-point attractor never happens. Instead “the collapse of complexity” happens and person #1 terminates the bicycle ergometry test.
It is equally interesting to observe the recovery processes right after the bicycle ergometry test is terminated. The inverse effect of the “collapse of complexity” is observed in Fig 8. The trajectory not only moves away from the point of collapse—the complexity of the trajectory itself is drastically changes over time.
Initially, the trajectory rapidly moves towards a temporarily stable attractor around 0.45 (Fig 8). However, this attractor loses the stability as the recovery processes continue. Finally, the system’s trajectory plots a chaotic attractor around 0.8. Such chaotic behavior of the heart could be considered as the validation of the presented model—since it is well known that heart is a chaotic system [29, 30].
A computational technique for visualization of complex transient processes of the heart parameters (RR and JT) is developed in this paper. It is based on the extensive application of novel algorithms and concepts originating from the broad field of nonlinear systems science and engineering. A natural question is whether such complex computations are required—maybe straightforward visualization of RR and JJ in a phase plane would reveal similar relationships?
The variation of RR—JT intervals during the bicycle ergometry test (for the first person) is visualized in Fig 9. Unfortunately, Fig 9 reveals essentially only a single relationship—RR (and JT) intervals do become shorter at the load.
This paper focuses on the development of a novel visualization technique of relationships between RR and JT intervals. But it is unclear if similar information on the self-organization of the heart system during the load and the recovery could not be retrieved from other cardiac intervals (for example RR and QRS; RR and QT). However, it appears that RR and QRS, RR and QT intervals yield completely different representations (S2 Appendix).
An important aspect of the presented computational technique is that an optimal set of parameters (s, m and Δ) is reconstructed for every individual person based on his RR, JT time series and the load data. The resulting optimal time series p(s, m, Δ) minimizes the individual load target function and is used for visual characterization of the self-organization of the heart system.
From the computational point of view, it would be possible to further minimize the differences between p(s, m, Δ) and the load target function by selecting optimal parameters s, m and Δ for each individual load step (not for the whole load history data). However, such an approach is not used due to two important reasons.
Firstly, the optimal time series p(s, m, Δ) would be a discontinuous function (at the time moments when the load is changed). That would complicate the interpretation of the graphical phase diagrams. Secondly, it would be unclear which parameter set should be used for the recovery. In other words, such splitting of the optimization procedure would not allow to visualize the recovery processes—which are considered to be equally (or even more) important than the load processes .
Another important aspect of this study is that we have identified the best fitting perfect matrix of Lagrange differences for the first person (s = 3). We do speculate that the particular number of the perfect matrix of Lagrange differences may serve as an identifier of dynamical processes taking place in the heart system during load and recovery.
The following computational experiment is performed for the illustration of this hypothesis. We use the same parameter (the maximum absolute value of the eigenvalues of a perfect matrix of Lagrange differences), the internal and the external smoothing (m = Δ = 3)—but instead of fixing the third matrix of Lagrange differences (s = 3) we perform the minimization of RMSE for every person for every possible perfect matrix of Lagrange differences (Table 4).
It is interesting to note that RMSE values for the second person (for all s values) are more than three time higher compared to the RMSE values for the first person (Table 4). But the same effect can be observed in Table 2 and could be probably explained by a guess that the self-organization of the heart system for the second person does function differently.
Nevertheless, it is possible to reconstruct a best fitting perfect matrix of Lagrange differences for every single person (Table 4; S1 Appendix). As mentioned previously, we hypothesize that the sequential number of this matrix is a specific identifier of the dynamical relationships and self-organization taking place in a particular person’s heart system.
However, the proposed technique for visualization of relationships between RR and JT intervals does not propose a single marker which can detect such conditions as atrial fibrillation, supraventricular tachycardia, or sudden cardiac death. Instead, we propose a computational technique which could reveal the complexity of the self-organization of the heart system during the load and the recovery processes. This complexity is represented by the parameter p(s, m, Δ). A physician can observe the “collapse of complexity” at the end of the bicycle stress test, temporary stabilization of transient attractors during the load, rich dynamical behavior of the heart system during the recovery process.
We doubt if it is possible to quantity such complex transient dynamics by a single biomarker. However, development of pattern classification algorithms for automatic analysis of transient orbits generated by p(s, m, Δ) (and relating these patterns to novel markers for early disease diagnosis) remains a definite objective of future research.
Target load functions, phase planes during the load and after the load for persons #2, #3, #4 and #5.
Other cardiac intervals yield different representation of the self-organizational processes.
Target load functions, phase planes during the load and after the load for persons #2 (parts a, b, c), #3 (parts d, e, f), #4 (parts g, h, i) and #5 (j, k, l).
Using RR-QRS intervals: target load function and pk(2, 3, 3) (part a), phase plane during the load (part b), after the load (part c), and straightforward mapping RR-QRS (part d) for person #1. Using RR-QT intervals: target load function and pk(9, 3, 3) (part e), phase plane during the load (part f), after the load (part g), and straightforward mapping RR-QT (part h) for person #1.
The supporting information files include unprocessed data from ECG analysis system “Kaunas-Load” for all subjects who were part of this research.
Financial support from the Lithuanian Science Council under project no. MIP078/15 is acknowledged.
Lithuanian University of Health Sciences is acknowledged for covering the open access publication fee. Financial support from the Lithuanian Science Council under project no. MIP078=15 is acknowledged.
All relevant data are within the paper and its Supporting Information files.