Sleep-disordered breathing: statistical characteristics of joint recurrent indicators in ЕЕG activity

Introduction

The search for functional biomarkers for automated detection of early stages of various diseases is one of the complex tasks of interdisciplinary data analysis. However, there are a large number of problems in this area. First of all, we must take into account that physiological signals naturally combine entire classes of signals from different systems of the human body. There are many methods aimed at attempting to separate such signals into relatively pure components [1]. However, these methods are often complex and do not always work correctly. In turn, the complexity of the mentioned methods can lead to false interpretation of the results, as well as to deliberately erroneous recognition of biomarkers of pathological processes by medical personnel [2].

In addition, methods originally developed for the analysis of stationary processes are often used to process physiological signals. The dynamics of living systems differs from model stochastic and chaotic systems in their greater complexity, accompanied by changes in both system control parameters and, apparently, continuous evolution and bifurcations of the internal structure of living systems as such. Finally, when processing signals from living systems, it is necessary to remember the individuality of the subjects. In particular, even the characteristics of invasive recording of brain activity in genetically homogeneous strains of laboratory rats have significant differences, falling into many subtypes [3]. Thus, automatic work with the functional activity of the brain of random patients, whose medical history, besides the diagnosis, may be burdened with a significant amount of comorbid pathology, is significantly complicated by the presence of “invisible” system parameters.

At the moment, most systems for automatic detection of biomarkers process, in fact, images of rather slow processes, such as, for example, well-known and used in medicine algorithms for detecting precursors of the development of benign and malignant neoplasms, skin lesions [4], medical decision support systems for radiography and magnetic -resonance tomography [5], etc. At the same time, the issues of monitoring the main characteristics of the functioning of the human cardiovascular system using the parameters of the basic rhythm of peak R waves on cardiograms have been relatively successfully resolved. The latter are well recorded by signal shape evaluation methods [6]. For example, automatic diagnosis of arrhythmias is now confidently carried out [7], the risk of developing episodes of atrial and ventricular fibrillation is identified [8], etc. In this regard, the interest of researchers in the search for simple informative biomarkers that are resistant to noise and do not require complex medical and computer equipment does not decrease [1, 9].

In this work, one of the simplest and at the same time effective methods of nonlinear dynamics was used - the construction of recurrence diagrams and their further numerical analysis. The first attempts to use recurrent methods of quantitative analysis to solve applied problems appeared in the 1980s in Ekman's publications [10]. Since then, methods based on recurrent analysis have been actively developed in various areas of biomedical signal analysis [11]. Recurrent analysis has a number of advantages over other methods - simplicity of calculations and the ability to work with a small set of signal points [12]. Moreover, other methods often require the assumption of initial stationarity of the experimental series under consideration [13, 14]. It is obvious that real signals of physiological systems, including EEG signals, do not at all satisfy the requirement of stationarity, exhibiting a complex set of features of dynamic systems with linear and nonlinear properties [15-17].

Because this study used polysomnography recordings, a key feature of which is the large experimental batch size, typically involving more than 7 hours of continuous recording at a sampling rate of about 500–1000 Hz, the signals were divided into small time windows, for which joint recurrent indicators and their statistical characteristics were calculated to highlight certain signal features that can be further used as biomarkers.

Material and Methods

Recurrence Analysis

Recurrent analysis allows you to establish relationships and correlations between signals in complex distributed systems. This method has found application in a wide range of tasks for processing complex signals of various natures [18, 19]. The calculation algorithm itself is extremely simple. Let us consider a signal x(t), the values of which are known only at discrete times t_i, where i=1, …, n. Let this signal x(t_i) be equidistant, i. e. t_i+1-t_i=t_i+1-t_i=∆t for any numbers i and j. Then consider constructing a recurrent plane as follows:

 , (1)

where R_i,j – element of the recurrent matrix for the signal x, t_i and t_j – time moments t, ε – an empirically determined threshold value that ensures the necessary accuracy of the method, θ — the Heaviside function, which is defined as:

 , (2)

Thus, the recurrent matrix, constructed according to expressions (1) and (2), is formed from elements of two types – «0» and «1». The matrix element is equal to «1» if the value of the signal x(t_i) at time t_i falls into the ε-neighborhood of the signal value x(t_j) at time t_j. At the same time, the matrix element is equal to «0» if the values of the signal x at times t_i и t_j are far from each other. These recurrent matrices (1) are often shown graphically in the form of recurrence diagrams, in which colored dots correspond to one values and white dots correspond to zero values of the matrix. Thus, the recurrent properties of the time series x(t_i) are represented in the form of geometric structures and allow us to visualize the dynamics of the series in the form of a simple graphical convolution.

Recurrence analysis includes methods for studying the location of points on the constructed surface of a recurrence diagram [20], which have been used in recent years to process stochastic time series of various natures [10, 21]. Further, with the development of machine learning methods, convolutional neural networks began to be used to directly recognize geometric structures appearing on recurrence diagrams [22, 23].

Note that in the case of single-frequency periodic dynamics in the recurrence diagram one can observe the resemblance of a grating, the period of which will correspond to the period of oscillation of the system [24]. In the case of multi-frequency periodic dynamics, superposition of gratings with different periods is observed. The greater the number of repetitions of a particular value, the more corresponding elements of the recurrent matrix are equal to «1». From this it is easy to conclude that the higher the oscillation frequency, the more points we get on the recurrence diagram. This fact makes it easy to identify the most frequently occurring values in the signal [11]. Therefore, recurrent analysis, although it does not belong to the group of frequency methods, allows you to automatically take into account the frequency of signal oscillations. To estimate the number of repetitions in the signal as a whole, the following recurrent indicator is used:

 ,       (3)

Such an indicator can be calculated for each analyzed signal x over the entire recorded length or over the required time fragment.

To compare two signals, a similar calculation of joint recurrence diagrams and joint recurrence indicators can be used. Formula (1) changes quite slightly:

,    (4)

where JR_i,j – element of the joint recurrent matrix for signals x and у, t_i and t_j - times t, ε and θ — have the same meaning as in the formula (1). Formula (4), thus, gives the values 1 for the elements of the joint recurrent matrix only if at moments t_i and t_j both signals x and у are in their ε-neighborhoods. Then, by analogy with formula (3), we can calculate the joint recurrent indicator:

 ,        (5)

This indicator is very useful, as it shows how often these signals demonstrate similar dynamics - returns of signal variables at the same points in time.

Statistical Metrics

In the course of this work, a large number of standard statistical metrics were calculated for the calculated joint recurrent indicators. This section provides the corresponding calculation formulas for all of them.

Average value

The formula for calculating the average is simple:

  ,     (6)

Harmonic mean

The harmonic mean is one of the ways in which one can understand the “average” value of some set of numbers.

,       (7)

Geometric mean

The geometric mean of several positive real numbers is a number that can be used to replace each of these numbers so that their product does not change.

,       (8)

Cubic mean

The cubic average is a characteristic of volumetric features. This is a special case of the power mean and therefore obeys the inequality about means. In particular, for any numbers it is not less than the arithmetic mean.

,      (9)

Median

The median or middle value of a set of numbers is the number that is in the middle of this set, if ordered in ascending order, that is, a number such that half of the elements of the set are not less than it, and the other half are not greater. Another equivalent definition: the median of a set of numbers is the number whose sum of distances (or, more strictly, moduli) from all numbers from the set is minimal.

Dispersion

Dispersion of a random variable is a measure of the dispersion of the values of a random variable relative to its mathematical expectation. The formula for calculating a biased estimate of the variance of a random variable from a sequence of realizations of this random variable has the form:

,    (10)

Standard deviation

The standard deviation is the most common indicator of the dispersion of the values of a random variable relative to its mathematical expectation (an analogue of the arithmetic mean with an infinite number of outcomes). Usually means the square root of the variance of a random variable, but sometimes it can mean one or another version of estimating this value.

,      (11)

The coefficient of variation

In probability theory and statistics, the coefficient of variation, also known as relative standard deviation, is a standard measure of the dispersion of a probability or frequency distribution.

,          (12)

Asymmetry indicator

The asymmetry indicator is a value in probability theory that characterizes the asymmetry of the distribution of a given random variable.

,          (13)

Kurtosis indicator

The kurtosis indicator in probability theory is a measure of the sharpness of the peak of the distribution of a random variable.

,        (14)

Polysomnography data obtained

The subjects were individuals with nocturnal respiratory dysfunction in the form of obstructive sleep apnea syndrome (N=30, age 48,0±19,1, median 43 years, male to female ratio = 18/12). Sleep duration was 6-9 hours, с 21.30-23.30 until the patient’s usual time of awakening.

Polysomnographic recording included electrocardiogram (ECG), respiratory function, oculography (OCG), electromyogram (EMG) and two electroencephalogram (EEG) signals recorded during night sleep. The ECG signal was recorded in standard lead I according to Einthoven. Respiration signals were recorded using a flow-through oronasal temperature sensor and a snoring sensor. EMG signals were recorded on the patient's chin, right forearm and left shin. OСG signals included recordings of horizontal and vertical eye movements.

EEG signals were recorded in 2 standard leads according to the 10-20 scheme. EEG signals were bandpass filtered 0.1-40 Hz and sampled at 500 Hz, ∆t=0,002 seconds. Registration of each EEG channel can be considered as a separate one-dimensional signal x(t_i) for subsequent recurrent analysis.

Based on these signals, joint recurrent indicators were constructed between EEG channels. The initial marking of the dependences of the joint recurrent indicator on time was carried out to identify special sleep events. If the joint recurrent indicator for a time of more than 150 seconds exceeded its average value by an amount greater than the variance (JRI_i>J̅R̅I̅+σ_JRR), then this event was marked, that is, the start and end times of the event were recorded. It was noted as a positive sleep abnormality. Events were similarly marked when the indicator for a time of more than 150 seconds was lower by an amount greater than the dispersion of its average value (JRI_i>J̅R̅I̅-σ_JRR). Such sleep events were noted as negative anomalies.

It is worth noting that the signal may contain individual extremes that go beyond the established boundaries, but they were not marked as events, since their length does not exceed 150 seconds. Subsequently, statistical characteristics were calculated separately for positive and negative anomalies in order to divide patient groups into clusters and identify possible biomarkers from these clusters.

Results

The work calculated a large number of statistical characteristics of joint recurrent indices (JRIs) corresponding to special sleep events. Characteristics for positive and negative sleep anomalies were considered separately. In this case, the patients were initially divided into two groups. The first group included patients with an apnea-hypopnea index less than 25. In the second group, the apnea-hypopnea index exceeds 25. Special sleep events for these groups were calculated separately to compare their statistical patterns and identify, if possible, simple linear classifiers for clustering according to these patterns.

Clustering was carried out based on the support vector machine and the k-means algorithm for each pair of statistical characteristics, constructed using joint recurrent indicators. The essence of the complete method: the centers of mass of the distribution for both groups are located in the two-dimensional space of statistical characteristics. Points located at a distance of more than three dispersion values from the center of mass are removed from further consideration to avoid the influence of statistical outliers on the results. Then a straight line is constructed through the centers of mass and the central point between the centers of mass lying on this straight line is calculated. A perpendicular is built through this point, which will separate the resulting clusters. It is this line that will be considered a linear classifier.

To assess the accuracy of this method, it is proposed to calculate a certain specially introduced coefficient µ. Its calculation involves measuring the distance to the linear classifier for each point of each group. The normal is lowered onto the classifier from the point and using the Euclidean measure, the distance from the point to the line is calculated. Since different characteristics have different distribution widths, the resulting distance is normalized to the distance from the center of mass to the origin of coordinates. The coefficient is the sum of all distances, however, if the point is on the side opposite the center of mass relative to the linear classifier, then the distance is taken with a minus sign. Thus, the greater the value of the coefficient µ, the better the resulting linear classifier divides groups into clusters. Figure 1 shows examples of successful and failed clustering.

Figure 1. Examples of successful and failed clustering.

а – an example of successful separation by harmonic mean and median for negative sleep anomalies. б – an example of erroneous division based on kurtosis and median for positive sleep anomalies. Red and green dots show the obtained statistical characteristics in two-dimensional space, and the black line shows the linear classifier.

For the example shown in Figure 1,а the coefficient µ takes the value 186.9, while for the example in Figure 1,б µ=14,29. In this case, the coefficient µ can also take negative values if the division into clusters was very unsuccessful.

Tables 1 and 2 show all coefficient values for negative and positive sleep anomalies. From these tables, we can identify those pairs of statistical characteristics for each type of sleep event that best separate these events and are the most likely candidates for the role of biomarkers.

Table 1. Values of the coefficient µ for all pairs of statistical characteristics for negative sleep anomalies

Table 2. Values of the coefficient µ for all pairs of statistical characteristics for positive sleep anomalies

Tables 1 and 2 show that the coefficient values are symmetrical relative to the main diagonal, which is left empty for obvious reasons. The filling symmetry is logical, since in this case the clusters will be identical up to the change of variables. It is worth noting, however, that in the work only values below the main diagonal were calculated, and values above the main diagonal were filled in mirror image. This remark is important, since Tables 3 and 4 show the coefficients for the linear classifier obtained for pairs of statistical characteristics located below the main diagonal.

For simplicity, to reduce the number of tables due to the symmetry of the data in Tables 3 and 4, the coefficients а of the linear classifier are located above the main diagonal, and the coefficients b below it. Thus, using Table 1, you can determine the most suitable statistical characteristics for dividing patients into groups of patients for negative sleep anomalies, while using Table 3 you can restore the type of linear classifier that was used to separate the data.

Table 3. Values of the linear classifier’s coefficients а and b (in the form y=ax+b) for all pairs of statistical characteristics for negative sleep anomalies

Table 4. Values of the linear classifier’s coefficients а and b (in the form y=ax+b) for all pairs of statistical characteristics for positive sleep anomalies

From the analysis of Table 1 it is clear that the most suitable characteristics for a biomarker are the average values (arithmetic, harmonic), as well as the median and, to a lesser extent, variance. Between them, as a rule, the coefficient µ takes values exceeding 100, which is a good result. In Table 2 it can be observed that for positive sleep anomalies, instead of dispersion, the asymmetry indicator allows clustering with average values. However, the coefficient µ achieves its greatest value precisely when constructing clusters based on dispersion and asymmetry index.

Discussion

The results list statistical characteristics that can be used as biomarkers to distinguish between patients with and without sleep apnea. However, these results are not enough to form full-fledged biomarkers. Firstly, the results can be improved by using nonlinear functions instead of a linear classifier. Some pairs that did not give a good separation into clusters according to the results of Tables 1 and 3 may give a good result when using non-linear classifiers, which will shift the priority for creating stable biomarkers.

Secondly, before forming a final opinion on the effectiveness of the division, it is necessary to consider the division into clusters according to the given classifiers in the multidimensional space of statistical characteristics. Thus, for negative anomalies it makes sense to construct a five-dimensional space from the arithmetic mean, harmonic mean, geometric mean, median and variance. If the division into clusters remains in five-dimensional space, then using these characteristics it will be possible to build a system for recognizing apnea syndrome.

It also makes sense to cluster according to three or four statistical characteristics, building a multidimensional space. As with a nonlinear classifier, in this case the result may change. Thus, this article is the first important step in identifying biomarkers for the early diagnosis of apnea, however, a lot of additional research is still required to create an effective system for recognizing apnea syndrome in the early stages.

Conclusion

In this article, joint recurrent indices were calculated for polysomnographic recordings of patients with apnea and statistical metrics of special sleep events. Metrics were calculated separately for positive sleep anomalies (when the joint recurrent indicator for a long time exceeds the average value of the indicator by an amount greater than its variance), and separately for negative anomalies (when the recurrent indicator is more than a variance less than its average value). Each calculation included mean, geometric mean, harmonic mean, dispersion, standard deviation, coefficient of variation, skewness index, kurtosis index. Each pair of statistical metrics was used to find a linear classifier by which patients with different apnea-hypnea index can be distinguished. To assess the quality of separation using this linear classifier, the coefficient µ was calculated based on the distance from the points of each cluster to the linear classifier. In addition to the calculated coefficients µ, the used parameters a and b for the linear classifier are given (with the classifier equation y=ax+b). A comparative analysis of the coefficient µ showed that the classifier associated with the median values works best. The asymmetry index and harmonic mean are also often effective for dividing patients into groups. These metrics are the most promising for searching for biomarkers of early diagnosis of apnea using polysomnography data.

Conflict of interest

The authors declare that they have no competing interests.

Funding

Study has been supported by of the Government Procurement of the Russian Federation Ministry of Healthcare within the state assignment “Development of algorithms for recognizing markers of breathing disorders during sleep in patients with various forms of cardiovascular pathology” No 122013100209-5 (2022-2024), performed in National Medical Research Center for Therapy and Preventive Medicine.