Cross-recurrence quantification of cardiovascular signals in newborns is a sensitive marker of health status

Introduction

Currently, to study the properties of complex systems, including in experimental studies, the strategy of studying the signals generated by the system is widely used. For example, researchers often use signals such as electrocardiogram, photoplethysmogram, blood pressure signal, and others to study the dynamics of the cardiovascular system.

Screening diagnostics of health status by analyzing noninvasively recorded signals from compact wearable devices is one of the promising areas for preventing the development of pathological changes and complications of disease severity [1, 2]. This issue is especially important when diagnosing diseases of newborns, in which the use of many instrumental research methods is difficult [3, 4]. A promising direction in such conditions is the analysis of RR-interval signals [5] and photoplethysmograms [6]. These signals can be recorded by small devices with compact ergonomic sensors and can provide valuable diagnostic information about the dynamics of the circuits in autonomic regulation of blood circulation [6-9].

Due to the complexity of the system under study, the set of traditional (linear) research methods needs to be expanded with new methods, including those from the theory of nonlinear dynamics and chaos. These methods were tested in practice; they help obtaining additional information about the properties of the system [7, 8].

One of such methods is cross-recurrence quantification (CRQ) analysis [10-11], which has shown its capabilities in the analysis of signals in living systems: in the diagnosis and classification of cognitive impairment [12], analysis of the relationship between heart rate and systolic blood pressure [13], diagnosis of multiclass tachycardias and ensemble classifiers [14]. However, the use of this method requires its adaptation to work with a specific type of experimental data of specific categories of subjects and a targeted selection of these method parameters.

Hence, the objective of our study was to adapt the method of analyzing signals from newborns and verify the results of its application when comparing very different experimental samples: healthy newborns and preterm neonates. We previously demonstrated the possibility of using CRQ analysis to investigate the degree of connection between the regulatory mechanisms of the cardiovascular system in healthy newborns and preterm neonates. We studied low-frequency (LF) range [0.04-0.14 Hz], reflecting primarily the processes of autonomic regulation of heart rate and arterial vascular tone [7].

Material and Methods

Signal recording

The study was conducted at the Department of Pediatrics and Neonatology of Saratov State Medical University (Saratov, Russia).

The study included two groups of neonates:

i) 10 full-term newborns (of a gestational age at birth of 37-40 weeks) with a physiological course of neonatal adaptation, selected on the basis of the healthy anamnesis of their mothers and normal births. Identification of diseases of any kind constituted grounds for disqualification.

ii) 10 preterm neonates (of a gestational age of 34 weeks gestation or older).

We carried out simultaneous recording of electrocardiographic (ECG) and photoplethysmographic (PPG) signals (blood filling signal recorded by an optical sensor with active backlight and adjustment for reflected light). Signals were recorded during feeding, thereby yielding data for the waking state.

Signals were recorded on a standard polygraph EEGA-21/26 Encephalan131-03 (Medicom MTD LLC, Taganrog, Russia). The bandpass filter was set to 0.04-100 Hz and the sampling rate was set to 250 Hz at 14-bit resolution. ECG signals were recorded in standard lead I (differential signal from the wrists of two hands), while PPG signals were recorded from the earlobe [6] using a reflected infrared light sensor (740 nm) [15]. Experimental signals were recorded at rest in a horizontal position. The duration of each recording was 10 minutes.

Information about the dynamics of the heart rate regulation system was obtained from the ECG signal, previously converted into a signal of RR intervals, which was further reduced to an equidistant form and filtered in the LF range [0.04-0.15] Hz. Information about the dynamics of the vascular tone regulation system was obtained from the PPG signal by filtering in the appropriate LF range [0.04-0.15] Hz.

Signal processing

The CRQ analysis method presumes the presence of two interacting nonlinear dynamic systems. Their interaction can be studied by analyzing the dynamics of these systems along their trajectories in a common phase space. The axes of the general phase portrait show the values of dynamic variables (in our study, x(t) means the LF components of the RR-interval signal, and y(t) means the LF components of the PPG signal) [16] (Figure 1a). In phase space, real biological systems exhibit complex dynamics similar to the dynamics of chaotic systems [17]. The image of the dynamics of chaotic systems in phase space is a set of nonintersecting non-closed trajectories [16] (Figure 1b).

Figure 1. Diagram of interacting dynamic systems X and Y (a) and their common (shared) phase space (b).

Bold dots mark close neighbors in the dotted volume of phase space and show the presence of connections between systems (trajectories remain close to each other for some time). The crosses illustrate an example of divergence of close trajectories at some point.

For complex systems of a biological nature, the connection between systems can be very complex and, in most cases, its mathematical description is difficult or impossible, despite regularly undertaken attempts to construct mathematical models of elements of living systems [18]. However, a universal feature of interacting dynamic systems is the tendency of points of coupled systems located close in phase space to move along similar trajectories for some time (thick points in Figure 1b). In the case of moderate coupling, chaotic systems are characterized by mixing, which over time causes the trajectories to diverge into different regions of phase space (crosses in Figure 1b). Therefore, if in a certain volume of phase space (dotted ellipse in Figure 1b), we identify close neighbors (close points) that are distant in time but located on trajectories close to each other, and if we observe them for some time, then it turns out that they will initially remain close (thick dots Figure 1b), but later they can diverge (crosses in Figure 1b).

CRQ approaches offer a tool for analyzing dynamics in a common nearest neighbor phase space and metrics for quantitatively describing such joint dynamics. The first step of CRQ analysis is to reconstruct the phase portrait of the system under study.

Phase portrait reconstruction: time delay method

As a rule, there are no adequate mathematical models of the physiological systems under study or they are too complex [19, 20, 18]. It is usually not possible to obtain data for all relevant variables. For example, in our studies, only experimental ECG and PPG signals were available. For CRQ analysis, it is necessary to reconstruct the phase portrait of each signal. However, the data collected for this single variable contains information about the dynamics of the entire system. Takens’ theorem [21] and corresponding extensions [22] provide the possibility of reconstructing the topological structure of a trajectory formed by state vectors from data for only one variable. For this purpose, multidimensional vectors are extracted from a univariate time series that captures observations of a single variable at discrete points in time.

The time delay method is habitually used to reconstruct system states from a univariate time series [23]. The dimensionality of the reconstructed vectors matches the number of corresponding variables and varies depending on the system. Figure 2 shows the application of the time delay method to one of the experimental time series of LF components of the PPG.

Figure 2. Time delay method.

A discretization of one of the experimental time series of LF components of the PPG (x-component) is presented. The corresponding univariate time series consists of nine data points (t₁ to t₉). The seven input vectors are reconstructed using embedding dimension m=2 and time delay t=2.

For demonstration, vectors consisting of two components are reconstructed. State reconstruction is based on the following parameters:

1. Embedding dimension (m),
2. Time delay (t).

The embedding dimension describes the dimension of the reconstructed vectors. The time delay parameter specifies the time shift of the vector components within the time series. The delay time can be determined in different ways: in this study, we obtain it by calculating the first zero crossing of the autocorrelation function.

Given a time series of l elements, embedding dimension m, and time delay t, the number of multidimensional vectors N is calculated as defined in Equation 1:

 ,                  (1)

To identify the optimal number of relevant variables for estimating the properties of the phase portrait, various methods, such as the false nearest neighbors (FNN) method, are used.

False nearest neighbor algorithm

The FNN method is based on Takens’ embedding theorem mentioned earlier. It follows that with an appropriate choice of the time step and the dimension of the embedding space, the original and reconstructed attractors must be topologically equivalent. Since the trajectories of the original attractor do not intersect themselves, the trajectories in the reconstructed attractor should not intersect themselves as well. The condition for the absence of self-intersections is that all neighboring points of an attractor reconstructed in a space of dimension m will also be neighboring in a space of dimension m+1. The FNN method allows determining the smallest value of the dimension m of the embedding space, so that when moving to the dimension (m+1), the number of false neighbors, attractor points close to each other in Rm and far apart in Rm+1, will be zero or relatively small. The value of dimension m obtained in this way determines the smallest dimension of the embedding space in which the attractor can be reconstructed without self-intersections. 

The false nearest neighbor algorithm consists of the following steps:

1. Let m=1. For each point x(i) of the time series, we find the nearest neighbor x(j) in m-dimensional space.

2. Then, we calculate the distance:

 ,                     (2)

3. In the next step, we find the distance between these points.

4. If Ri>Rt, where Rt is the appropriate threshold, then the point x(j) is a false near neighbor of the point x(i). As a result, the number of such false nearest neighbors P is calculated for each point x(i).

5. P/N is calculated, and the algorithm is repeated for m=m+1.

6. The algorithm continues until P/N becomes close to zero. In this study, the criterion for stopping the algorithm was the following condition: P/N<10%.

Estimation of the optimal dimension for the signals of the LF components in RR and PPG, as well as the reconstruction of their phase portraits, allows using CRQ.

Cross-recurrence plot

After reconstructing the attractors of RR intervals and PPG in the same phase space, we checked the proximity of each point of the first trajectory of RR intervals (i = 1, ... , N) with each point of the second trajectory of RR intervals (i = 1,. ..,N) with each point of the second trajectory PPG signal (j = = 1, ... , N), as a result of which an N×N array was formed known as cross-recurrence plot (CRP) (3). Examples of CRP are shown below in Figure 3.

 ,         (3)

Figure 3. Graphic representation of diagonal lines on cross-recurrence plot.

where Q is the Heaviside step function, ε is the threshold distance and ||…|| is the norm. ε was set to 10% and the norm was a fixed number of neighbors. Using this criterion, the neighborhood is a ball with a radius that changes for each dynamic state in order to encompass a fixed number of states. Then the radius ϵ changes for each dynamic state, since the neighborhoods of different states do not have to be the same. Using this criterion, the parameter ϵ has a predetermined recurrence density (in our case, ϵ was set to ensure that 5% of all phase space vectors are in the neighborhood of each state). The type of neighborhood that should be used depends on the application. For example, FNN selection is useful for non-stationary data and for bivariate recurrence studies via CRQ analysis.

The calculated cross-recurrence plot values are used to evaluate various characteristics. In our study, we used two main characteristics related to assessing the degree of connectivity between the systems under study.

The mean diagonal line L is defined as:

 ,          (4)

while the length of the longest diagonal line found in CRP is:

 ,      (5)

where  is the total number of diagonal lines.

Below, we present a variant of the CRP procedure with a graphical explanation of the appearance of diagonal lines on it (Figure 3).

Results

Experimental signals of RR intervals and PPG for two samples (healthy newborns and preterm neonates as shown in Figure 4 (a-b, g-h)) were filtered in the LF range of [0.04-0.15] Hz. This procedure allowed us leaving only those components that were traditionally associated primarily with the activity of the sympathetic system for regulating heart rate and vascular tone, respectively. In Figure 4, the column marks the frequencies that remained in the LF components of the RR intervals and PPG. Examples of signals after filtering for two samples of newborns are shown in Figure 4 (e-f, k-l).

Figure 4. Examples of initial signals in RR intervals and PPG are shown for two samples: (a-b, healthy newborns; g-h, preterm neonates), spectra of initial signals (c-d, healthy newborns; i-j, preterm neonates), filtered signals in the LF range (e-f, healthy newborns; k-l, preterm neonates).

The required dimension of the phase space for the LF components of the RR interval and PPG signals was estimated. The results are presented in Figure 5.

Figure 5. Embedding dimension estimation for 10 RR intervals and 10 PPG signals (thin black lines correspond to the estimate of each individual time series, while thick red line corresponds to the mean value in each embedding dimension): a and c are RR intervals for healthy newborns and preterm neonates, respectively; b and d are PPG for healthy newborns and preterm neonates, respectively.

It can be seen that for all groups of newborns, dimension 4 is sufficient to restore phase spaces, which can be explained by the relatively poor dynamics of signals in the LF range.

We calculated CRP for signal pairs of RR intervals and PPG for two samples: healthy newborns and preterm neonates. Examples of the resulting CRP are shown in Figure 6. It shows diagonal lines corresponding to the close movement of the phase trajectories of the two analyzed systems in a shared phase space.

Figure 6. The result of calculating a cross-recurrence plot for two paired signals of RR intervals and PPG.

In this study, we calculated quantitative metrics characterizing the tendency of phase trajectories to remain close in order to quantitatively describe the collective dynamics of the systems under consideration. The results of the quantitative analysis are presented in Table 1.

Table 1. Quantitative CRP estimates for two samples: healthy newborns and preterm neonates

Discussion

Our previous publication [24] presented the results of a study on phase synchronization of LF processes in RR intervals and PPG in healthy newborns. We demonstrated the presence of functional interaction between elements of the autonomic control of heart rate and vascular tone. However, the degree of synchronization between these systems was significantly lower than in adults. This could reflect the immaturity of the elements of autonomic regulation immediately after birth. Also, in another publication [25], we showed the results of a study on the dependence of S (total fraction of phase synchronization) on gestational age. The results of the work did not reveal a pronounced dependence of S on the gestational age. S values in newborns were significantly lower than in healthy adult volunteers, indicating a significantly weaker connection between the studied processes [26].

In our present study, we examined connectivity using the method of CRQ analysis, which in the well-known publications of other researchers was confirmed to be a sensitive method for diagnosing weak connections in the study of non-stationary noisy signals of a physiological nature [27]. The results of a quantitative assessment of the connectivity between the studied systems expressed via L and Lmax yielded higher values for the group of healthy newborns. The values of these indicators can be interpreted as indices characterizing the strength of effective communication between systems [27].

Thus, our research demonstrated that the relationship between the studied processes was very weak in preterm neonates and much stronger in healthy newborns, which brings the latter group closer to healthy adult volunteers. These findings may indicate that in healthy adults, there is significant coordination of oscillations in the circuits of autonomic control of heart rate and arterial tone, which is less pronounced in healthy newborns and minimal in preterm neonates due to the immaturity of the autonomic control components. This, on the one hand, is a potential biomarker of the state of the regulatory circuits in newborns with pathology. On the other hand, insufficient coordination of autonomic regulation processes in preterm neonates may be an independent cause of aggravated pathological processes in their body, which may require special attention from doctors.

Conclusion

In our study, we tested a nonlinear method for analyzing the connectivity of systems in order to assess the degree of synchronization of the autonomic control processes in the cardiovascular system of healthy newborns and preterm neonates. In contrast to previously published data on the absence of a pronounced trend in the overall fraction of phase synchronization S in preterm neonates vs. healthy newborns, quantitative assessment of the mean diagonal lines on CRP (L) showed a slight but distinct increase in the group of healthy newborns vs. preterm neonates. Hence, the method of CRQ analysis exhibited greater sensitivity compared with the previously used method for estimating the overall fraction of phase synchronization, which could be caused by its greater robustness to noise.

Funding

The presented research was funded by the Russian Science Foundation under scientific project No. 22-75-00086.

Ethical approval

All procedures performed in the studies involving human participants were in accordance with the ethical standards of the institutional and national research committee and with 1964 Declaration of Helsinki and its later amendments or comparable ethical standards. The design of this study was approved by the Ethics Committee of Saratov State Medical University (Saratov, Russia) in 2021. All parents were informed about the nature and objectives of the study and provided their consent before participation.

Conflict of interest

The authors declare no conflicts of interest.