Voronoi Decomposition of Cardiovascular Dependency Structures in Different Ambient Conditions: An Entropy Study

Abstract

This paper proposes a method that maps the coupling strength of an arbitrary number of signals D, D ≥ 2, into a single time series. It is motivated by the inability of multiscale entropy to jointly analyze more than two signals. The coupling strength is determined using the copula density defined over a [0 1]^D copula domain. The copula domain is decomposed into the Voronoi regions, with volumes inversely proportional to the dependency level (coupling strength) of the observed joint signals. A stream of dependency levels, ordered in time, creates a new time series that shows the fluctuation of the signals’ coupling strength along the time axis. The composite multiscale entropy (CMSE) is then applied to three signals, systolic blood pressure (SBP), pulse interval (PI), and body temperature (t_B), simultaneously recorded from rats exposed to different ambient temperatures (t_A). The obtained results are consistent with the results from the classical studies, and the method itself offers more levels of freedom than the classical analysis.

1. Introduction

Approximate [1,2] and sample entropies [3], ApEn and SampEn, have been intensively implemented in a range of scientific fields to quantify the unpredictability of time series fluctuations. Contributions that apply ApEn and SampEn are measured by thousands [4], confirming their significance. The cross entropies—XApEn and XSampEn—are designed to measure a level of asynchrony of two parallel time series [3,5,6]. Descriptions of (cross) entropy concepts can be found in numerous articles, but a recent comprehensive review [7] provides an excellent tutorial with the guidelines aimed to help the research society to understand ApEn and SampEn and to apply them correctly [7].

Multiscale entropy (MSE) [8,9], based on SampEn, investigates the changes in complexity caused by a change of the time scale. Composite MSE (CMSE) performs an additional averaging, thus solving the problem of decreased reliability induced by temporal scaling [10,11]. A comprehensive study of fixed and variable thresholds at different scales also presents an excellent review of the MSE improvements [12].

The benefits offered by entropy are explored in cardiovascular data analysis. Entropy was implemented to determine the cardiac variability [13], the complexity changes in cardiovascular disease [14], a level of deterministic chaos of heart rate variability (HRV) [9], HRV complexity in diabetes patients [15], in heart failure [16], in stress, [17,18] or in different aging and gender groups [19,20], while multiscale cross-entropy was applied for health monitoring systems [11].

SampEn- or ApEn-based entropy estimates are designed for one signal, or at most for two signals (cross-entropy), but biomedical studies often require an analysis of three or more simultaneously recorded signals.

We propose a method that maps levels of interaction of two or more time series into a single signal. Levels of interaction are assessed using the copula density [21]. The transformation from the probabilistic copula domain to the beat-to-beat time domain is performed by Voronoi decomposition.

The method is applied to multivariate time series that comprises three simultaneously recorded signals: systolic blood pressure (SBP), pulse interval (PI), and body temperature (t_B) recorded at different ambient temperatures (t_A). It is well known that thermoregulation can affect cardiovascular homeostasis [22]. Analysis of heart rate (HR) and SBP in the spectral domain has shown that changes of ambient temperature modulate vasomotion in the skin blood vessels, reflected in the very-low-frequency range of SBP and reflex changes in HR spectra [23,24]. Thermoregulation is complex and involves autonomic, cardiovascular, respiratory as well as a metabolic adaptation [25,26,27,28]. The key corrector of blood pressure is the baroreceptor reflex (BRR). The disfunction of BRR is the hallmark of cardiovascular diseases with a bad clinical prognosis. Thus, evaluating its functioning is important not only for the diagnosis and prognosis of cardiovascular diseases but also for the evaluation of treatment.

The aims of this study are:

• To propose a method that enables an application of multiscale entropy to an arbitrary number of signals and to analyze the outcome;
• To compare the results of the classical multiscale method and the proposed method when applicable, i.e., in a case of two-dimensional signals;
• To test whether the proposed method recognizes the changes of dependency level (coupling strength, level of interaction) of joint multivariate signals in different biomedical experiments.

The paper is organized as follows: the experimental setting for signal acquisition is explained in Section 2.1, together with surrogate signals and artificially generated control signals. The signal pre-processing that ensures the reliability of the results is explained in Section 2.2. Section 2.3. shows the mathematical tools assembled to create the proposed method: it gives an introduction to the copula theory, it outlines copula advantages and applications, and it discusses the various procedures for density estimation to justify the preference of Voronoi decomposition.

Section 3.1. shows the basic statistical analysis of the experimental SBP, PI, and t_B signals. For the sake of comparison, this section includes the outcomes of classical (X)SampEn and CMSE entropy analysis. Section 3.2. introduces the new signal, created by the proposed method, for a two-dimensional case (SBP and PI mapped into the new D = 2 signal) and a three-dimensional case (SBP, PI, and t_B mapped into the new D = 3 signal). In both cases, the SBP-PI offset (delay) is taken into account ranging from 0 to 5 beats [29]. The wide sense stationarity of the created signals is checked and the correction proposed. The signals’ statistical properties, in terms of skewness and kurtosis, are estimated and discussed. In Section 3.3., the entropy parameters are analyzed and the proper ones that ensure the reliable estimates are selected. Then, the results of experiments performed to justify the consistency with the classical methods (in cases when the comparison is possible) are presented. The results showing that the method recognizes the changes in the level of signal interaction in various experimental environments are presented as well. The results are discussed in Section 4 with respect to the aims of this paper. The same section gives the conclusion and the possibilities for further method applications.

A brief description of well-known entropy concepts—ApEn, SampEn, MSE, and CMSE—is included in the Appendix A.

2. Materials and Methods

2.1. Experimental Setting and Signal Acquisition

All experimental procedures conformed to Directive 2010/63/EU National Animal Welfare Act 2009/6/RS and Rule Book 2010/RS. The protocol was approved by the University of Belgrade Ethics review board (license n°323-07-10519/2013-05/2).

Adult male Wistar outbred rats, weighing 300–350 g, housed under control laboratory conditions (temperature—22 ± 2 °C; relative humidity: 60–70%; lighting: 12:12 h light-dark cycle) with food (0.2% NaCl) and tap water ad libitum were used in experimentation. Vasopressin selective antagonists of V1a or V2 receptors were injected via cannula chronically positioned in the lateral cerebral ventricle of the rat. The concomitant measurement of blood pressure waveforms (BP) and body temperature was performed using TL11M2-C50-PXT (DSI, St. Paul, MN, USA) equipment implanted into the abdominal aorta. The measurements were performed at the neutral ambient temperature (NT), 27 rats at t_A = 22 ± 2 °C, and the increased ambient temperature (HT), 28 rats at t_A = 34 ± 2 °C. The four rats recorded at the low temperature (LT), t_A = 12 ± 2 °C, were included as an illustrative example. There are five subgroups in NT and HT groups: control group, V1a-100 ng, V1a-500 ng, V2-100 ng, and V2-500 ng. The experimental timeline is shown in Figure 1.

The experimental environment includes two types of control signals. The first controls are isodistributional surrogate data [30,31]. Surrogate data are derived from the experimental time series by randomizing the property that needs to be tested, keeping the other signal attributes intact. Thus, isodistributional surrogates randomly permute the signal to destroy the orderliness that is checked by entropy analysis. The signal distribution function remains unchanged. The second controls are artificially generated signals—a series of independent and identically distributed (i.i.d.) samples with Gaussian distribution and with exponential distribution. Gaussian signals possess a unique property in which linear independency implies statistical independency [32], and that it is an asymptotic distribution of the sum of i.i.d. samples (with some constraints) [32], an issue important for the multiscale entropy coarse-graining. Signals with exponential distribution are often implemented when there is a need to test the signals with large variance.

2.2. Signal Pre-Processing

Arterial blood pressure (BP) and body temperature signals were acquired using a sampling frequency of 1000 Hz. Systolic blood pressure (SBP) and pulse interval (PI) time series were derived from the BP waveforms as the local maxima and as the intervals between the successive maximal BP positive changes, respectively. The samples from the body temperature signals were taken simultaneously with SBP to create body temperature beat-to-beat time series t_B. Artifacts were detected semi-automatically using the filter [33] adjusted to the signals recorded from the laboratory rats. A visual examination was then performed to find the residual errors. A very low signal component (trend) was removed by a high-pass filter designed for biomedical time series [34], thus ensuring SBP, PI, and t_B signal stationarity. All the signals were cut to the length of the shortest time series, n = 14,405 samples. The time series X₁ = SBP, X₂ = PI and X₃ = t_B jointly create a single three-dimensional signal (D = 3). Its samples X_1k, X_2k, and X_3k, k = 1,…, N create points in the three-dimensional signal space.

2.3. Copula Density, Voronoi Regions and Dependency Time Series

A copula is a mathematical concept that provides a multidimensional probability density function, where density reflects the level of signal interaction (dependency, coupling). It is introduced in 1959 [21] as a multivariate distribution function with marginals uniformly distributed on [0 1]^D. If X₁, …, X_D are the source signals with joint distribution function H and univariate marginal distribution functions F₁,…, F_D, then copula C is defined as [21]:

$$H\left(X_1,\dots , X_D\right)=C\left(F_1\left(X_1\right),\dots ,F_D\left(X_D\right)\right)$$

(1)

and vice versa:

$$C\left(U_1,\dots U_D\right)=H\left(F_1^{-1}\left(U_1\right),\dots ,F_D^{-1}\left(U_D\right)\right).$$

(2)

Sklar’s theorem [21] states that any D-dimensional joint distribution H with arbitrary univariate marginals could be decomposed into D independent uniform marginal distributions, bound together by a new joint distribution function C, called copula.

The concept of the copula is based on the classical transformation of a random variable. Any continuous variable $X_i$ with a distribution function $F_{i }(X_i$) and density $f_{i }\left(X_i\right)=\frac{d F_{i }\left(X_i\right)}{d{X}_{i }}, i=1,\dots ,D$ can be transformed using a monotone function $U_{i }={\phi }_{i }\left(X_{i }\right)$. The result is a variable $U_i$ with a probability density function [32] $u_i(U_{i })=\frac{ f_i(X_{i })}{\left|d\left({\phi }_{i }\left(X_{i }\right)\right)/d{X}_{i }\right|}, X_{i }= {\phi }_{i }{}^{-1}\left(U_{i }\right).$ The transformation function ${\phi }_{i }\left(X_i\right)$ that creates a copula is the distribution function $F_i\left(X_{i }\right)$ of the signal $X_i$, i.e., ${\phi }_{i }\left(X_{i }\right)= F_i\left(X_i\right)$. The new variable $U_i$ is then defined in [0, 1], as the following holds: $0\le F_i\left(X_{i }\right)\le 1$. It can be easily shown that the probability density function (PDF) of the new variable $U_{i }$ is uniform:

$$u_i(U_{i })=\frac{ f_i(X_{i })}{\left|d\left(F_i\left(X_{i }\right)\right)/d{X}_{i }\right|}=\frac{ f_i(X_{i })}{ f_i(X_{i })}=1.$$

(3)

The transformation of a random variable using its distribution function is known as probability integral transform, PI-transform, or PIT [35], and it is a core of the copula theory. It should be noted that the distribution function of a continuous variable, by definition, monotonically increases so the denominator in Equation (3) is positive, comprising just a single term.

The copula has been intensively used for the analysis and prediction of financial time series and the prediction of insurance risk [36,37,38,39], in hydrology and climate analysis [40,41,42] and communications [43]. Medical applications include aortic regurgitation study [44] and diagnostic classifiers for neuropsychiatric disorders [45]. A possibility to use a bivariate copula to analyze the cardiovascular dependency structures was introduced in [46] and pharmacologically validated by blocking the feedback paths using Scopolamine, Atenolol, Prazosin, and Hexamethonium. It was shown that Frank’s copula is the most appropriate to quantify the level of dependency of cardiovascular signals.

Copula density $c\left(U\right)=\frac{{\partial }^{N}C\left(U_1,\cdots , U_D\right)}{\partial U_1\cdots \partial U_D}$ is used to visualize the intensity of signal coupling. The regions of increased copula density indicate the regions where the dependency of the signal samples increases. The difference between a classical bivariate probability density function (PDF) and the corresponding copula density of SBP and PI signals is that PDF shows the distribution of amplitude levels, while copula density shows the distribution of coupling strength between these amplitudes, regardless of the absolute amplitude values. An illustration of this difference is shown in Figure 2. SBP and PI signals and their probability integral transformed (PIT) counterparts are separated in time, first by two heartbeats (SBP_k is coupled with PI_{k + 2}, k = 1, 2, …, N−2), and then by ten heartbeats (SBP_k is coupled with PI_{k + 10}, k = 1, 2, …, N−10). The copula density in panel b exhibits a distinct linear positive coupling structure that follows the known physiological relationships [47]. The copula density in panel d shows almost uniform distribution as the time offset between SBP and PI signals is sufficiently large to attenuate their mutual dependency. Contrary to copula density, the joint probability density functions are almost the same in both cases (panels a and c). The temporal separation of SBP and PI signals does not alter the mutual relationship of signal amplitudes, but it significantly alters the intensity of signal coupling.

The advantages of copula are numerous. Copula density visualizes the dependency structures of the observed signals, and it quantifies the signal coupling strength (“copula parameter”). It captures both linear and nonlinear relationships between the signals. It can quantify the intensity of signal coupling within the different regions of the copula domain, and, in particular, it can model the tail dependencies of the signals.

Such a visualization, in a case of SBP-PI signals, cannot be achieved by other methods: the Oxford method, the oldest and the referent procedure for the evaluation of the baroreceptor reflex, uses increasing doses of short-acting vasoconstrictors (e.g., phenylephrine) and vasodilators (e.g., nitroprusside) to trigger heart deceleration of acceleration. The SBP and PI relationship is plotted as a fitted sigmoid curve. It is an invasive method, and it does not show spontaneous BRR. The most acknowledged among the non-invasive approaches is the sequence method, with the visualization that shows the scatterplot of the signal points that are elements of BRR sequences (i.e., the scatterplot contains a subset of all signal points). The method quantifies the spontaneous BRR operating range and set point [48], but the visualization is similar to the classical probability density function.

The time offset (delay) between SBP and PI signals is important for the signal coupling, and it depends on species. It was shown [49] that the delay of 0, 1, and 2 beats is the most appropriate for humans, while the delays of 3, 4 and 5 beats are appropriate for rats [29] and mice [47]. In [50], it was shown that, in laboratory rats, the highest level of comonotonic behavior of pulse interval and systolic blood pressure is observed at time lags 0, 3, and 4 beats, while a strong counter-monotonic behavior occurs at time lags of 1 and 2 beats.

Copula density is a probabilistic quantity. To convert it into a time series, to each point in the time domain, an appropriate density (dependency level) $D{L}_k$ should be assigned, thus creating a dependency signal $D{L}_k,k=1,\dots , N$.

A trivial way to estimate a copula density is to create a D-dimensional histogram. The obtained $D{L}_k$ signal would be discrete, as the points within the same histogram bins would get the same value. An increased number of histogram bins would increase the number of discrete signal levels, but the estimation reliability would decrease.

Density estimation based on Markov chains [51] creates a stochastic matrix of “transition probabilities”—scaled distances—between the points, with the steady-state probabilities proportional to the required distribution. The method is computationally inefficient in multidimensional space, except for the short time series.

A D-dimensional sphere (or cube) around a particular signal point defines a local density according to the number of encircled neighbors. The procedure is efficient, but the neighboring spheres overlap inducing the bias, and the result depends on the sphere diameter (i.e., threshold) choice.

The chosen approach expresses the sample density proportionally to the non-overlapping free space surrounding the sample. Such a concept has long been known as the Voronoi region. It can be traced back to the scholars from the 17th and 18th centuries, but it was re-discovered, analyzed, and its applications outlined at the beginning of the 20th century [52].

The concept is simple: Let A be the set of all points in a [0 1]^D copula space. Let $U_k=[U_{1k},\dots , U_{Dk}], k=1, \dots , N$ be a D-dimensional point from a PI-transformed multivariate time series. Then, the Voronoi region $R_k^D$ around the point $U_k$ comprises all the points from A that are closer to the particular point $U_k$ then to any other point $U_j, j=1, \dots , N, j\ne k.$ More formally,

$$R_k^D=\{a\in A | d\left(a,U_k\right) \le d\left(a,U_j\right), \forall j \ne k\}.$$

(4)

A classical Euclidean distance is typically chosen for the distance $d\left(a, U_k\right)$, but any other distance measure can be used as well, resulting in different Voronoi decompositions.

Figure 3 shows examples of the Voronoi regions in two and three dimensions. The line segments that separate particular Voronoi cells $R_k^2$ and $R_j^2$ in the left panel of Figure 3 are the sets of the points $a\in A$ that are equidistant to the points $U_k$ and $U_j$, i.e., $d\left(a, U_k\right)=d\left(a,U_j\right).$ The Voronoi vertex $a\in A$ in the same panel is the point equidistant to three (or more) time series points, e.g., $d\left(a, U_k\right)=d\left(a, U_j\right)=d\left(a, U_l\right)$. The right panel (D = 3) also shows Voronoi lines and vertices, but, in the [0 1]³ domain, this is more difficult to visualize. Uncolored Voronoi regions are either unbounded, or the boundaries are outside the [0 1]^D space. These regions are cut to fit the [0 1]^D space.

A series of surface areas in two-dimensional Voronoi regions and a series of volumes in three-dimensional Voronoi regions are a good foundation to quantify the dependency level and to form the time series $D{L}_k,k=1,\dots , N, \mathrm{as}$:

(a)

The surface/volume of $R_k^D$ is inversely proportional to the dependency level of the point $U_k$. An increased density of dependency structures in [0 1]^D space implies a decrease of available space between the points.

(b)

The region $R_k^D$ is shaped like the best distance separation of the point $U_k$, so its surface/volume is unambiguously calculated and unique, without a necessity to include any thresholds.

The drawback of the method is that a change of distance measure changes the shape of regions. We have opted for Euclidian distance as a classical approach for distance measurement, widely used in a wide range of applications.

3. Results

3.1. Source Signal Analysis

The total number of SBP-PI-t_B signal triplets is equal to 59. The basic statistical parameters, shown as a control, are presented in

Table 1. Statistical parameters of the source data (mean ± standard deviation).

Ambient Temperature (°C)	Drug	SBP	(mmHg)	PI	(ms)	tB	(°C)
NT 22 ± 2	Control	112.81	±19.54	179.22	±33.22	38.07	±0.29
	V1a, 100 mg	115.62	±12.17	173.74	±20.69	38.42	±0.10
	V1a, 500 mg	110.28	±15.35	184.77	±28.39	38.05	±0.10
	V2, 100 mg	119.98	±16.53	184.79	±38.30	38.54	±0.38
	V2, 500 mg	108.61	±14.79	176.16	±4.04	38.33	±0.41
HT 34 ± 2	Control	107.26	±4.19	188.63	±8.95	38.27	±0.34
	V1a, 100 mg	107.90	±10.52	197.08	±21.63	38.52	±0.26
	V1a, 500 mg	110.40	±10.07	177.21	±16.34	38.57	±0.57
	V2, 100 mg	113.26	±15.41	193.14	±30.65	38.01	±0.37
	V2, 500 mg	114.28	±6.14	184.23	±12.97	38.33	±0.47
LT 12 ± 2	Control	115.22	±5.23	164.54	±24.31	37.51	±0.43

Note: Results are presented as mean ± standard deviation; SBP: systolic blood pressure; PI: pulse interval; t_B: body temperature; NT: neutral temperature; HT: high temperature; LT: low temperature.

. Results in Table 1 show no significant changes in statistical parameters of SBP, PI, and t_B signals. An earlier study [26] revealed that V1a antagonists increase body temperature. The differences might be the outcome of different measurement procedures: in this study, the temperature is measured using a telemetric probe in the abdominal aorta, while, in [26], the temperature was measured rectally.

Figure 4 shows the results of the classical entropy analysis, performed for the sake of comparison. The left panels show the CMSE of the signals recorded from control animals at different ambient temperatures. The middle panels show the effect of drugs at the neutral temperature. The left panels show the effects of drugs at a high temperature. Each signal is accompanied by ten isodistributional surrogate signals, generated by a random temporal permutation of the signal samples [30,31].

The first three rows in Figure 4 show the classical composite multiscale entropy analysis of a single-dimensional time series, SBP, PI, and t_B, respectively. The last row shows multiscale SBP-PI cross-entropy that can be compared to the multiscale entropy of the new signals.

3.2. Properties of the Dependency Time Series

The created time series are new signals, so their statistical properties need to be checked before entropy analysis.

Mapping the signals into the dependency time series takes into account the delay (offset) between the PI and SBP signals. The time delay (offset) $DEL=0,\cdots , 5$ [beats] applied to each pair of SBP-PI signals resulted in six two-dimensional (2D) time series $\left(SB{P}_k, P{I}_{k+DEL}\right)$, and six three-dimensional (3D) time series $\left(SB{P}_k, P{I}_{k+DEL},t_{\mathrm{B}k}\right), k=1,\dots , N-DEL.$ The total of 354 SBP-PI pairs and 354 SBP-PI-t_B triplets were converted into two-dimensional and three-dimensional Voronoi cell time series. An average percentage of Voronoi cells that had to be cut to fit the [0 1]^D space was 2.68% for two-dimensional, and 16.26% for three-dimensional signals (cf. Figure 3). Additionally, 11 signal points (0.0002%) were too close to the vertices of the [0 1l³ cube to generate the three-dimensional polyhedrons, so they were managed manually.

Examples of Voronoi cell time series are shown in Figure 5.

A wide sense stationarity (WSS) test [53] is then applied, as the stationarity is an obligatory prerequisite for entropy estimation [7,38]. The test checks the stationarity of the first and the second statistical moments. The three-dimensional Voronoi cells time series failed the second-moment test.

Figure 6 shows the negative effects of non-stationarity: a three-dimensional non-stationary Voronoi cell time series is cut into 14 successive segments, each one comprising n = 1000 signal points. Then, mean, variance, and entropy were estimated from each segment and plotted in Figure 6a. Figure 6b shows the same parameters but estimated from the two-dimensional stationary Voronoi cells time series.

The difference between two- and three- dimensional signals is a consequence of coverage. The number of signal points in two-dimensional space is sufficient to ensure good coverage. The same points are sparsely and unevenly scattered in three-dimensional space, so the estimation is unreliable, resulting in different values obtained from the different sections of the same signal.

Panel d in Figure 6 shows the results from artificially generated time series with exponentially distributed samples. It is a usual control example of a signal with large, but time-invariant, variance. It passed the stationarity test and the parameters estimated from it are constant in each segment.

Taking a logarithm is a procedure that ensures the stationarity of the second moment. A negative logarithm corresponds to the inverse of the Voronoi cell volume, and it is proportional to the local sample density. It is always positive as the inverse of any Voronoi cell volume in [0 1]^D domain is greater than 1. Panel c of Figure 6 is a visual confirmation of a successful test outcome.

The dependency level time series, DL, is finally defined as the negative logarithm of the Voronoi cell time series. The number of signal points (N = 14,400) ensures the signal stationarity at least in the wide sense for two-dimensional (D = 2) and three-dimensional (D = 3) signals.

The statistical properties of the new signals—probability density function, skewness, and kurtosis—are shown in Figure 7 and Figure 8.

Empirical probability density functions for different scaling levels are plotted in Figure 7. Signals are normalized and centered so the changes of mean and variance due to the convolution are not visible. There is no significant difference in the estimated probability density functions for the scaling levels greater than five.

Skewness is a third statistical moment that shows a level of signal asymmetry around the mean. The skewness of the two-dimensional dependency signal (SBP and PI interaction) is presented in Figure 8a). It is positive, with a right tail exhibited, indicating the existence of signals with a strong dependency level between systolic blood pressure and a pulse interval. The positive skewness increases with the increasing offset (delay) DEL between SBP in PI. It is in accordance with [29] that located the dominant SBP-PI relationships at offset of 3, 4, and 5 beats.

The skewness of the three-dimensional dependency signal (SBP, PI, and t_B interaction) is presented in Figure 8b). It is close to zero—slightly negative—so the level of the dependency between the three signals is almost symmetric. It may indicate that the inclusion of body temperature into the new signal attenuates the SBP-PI signal coupling. The increase of the SBP-PI offset (delay) DEL results in the increased skewness shifted closer to zero, towards the positive values, again in accordance with [29].

Kurtosis measures the intensity of probability density function “tails”. It is shown in Figure 8b for two-dimensional signals (SBP and PI interaction) and in Figure 8d for three-dimensional signals (SBP, PI, and t_B interaction). The tails of dependency signals are heavy if compared to Gaussian distribution, indicating an increased number of signals with very high and very low dependency levels. It is expected, due to the high variance of three-dimensional dependency signals. The intensity of tails increases with the increased SBP-PI offset DEL, and it also increases with an increase of scale. This is also expected, as scaling convolves the probability density functions of the components that are coarse-grained. The convolution emphasizes the tail parts of the distribution in spite of the normalization, as the convolved samples are not Gaussian.

3.3. Entropy Analysis of the Dependency Time Series

Entropy is a parametric method, with parameters determined to ensure its reliable estimation. The statistical analysis from Section 3.2, however, is insufficient to provide the guidelines for entropy parameter selection for the new dependency time series.

It has already been pointed out [18,54,55,56,57,58,59] that the threshold (filter) r (cf. Equation (A3)) is one of the major causes of inconsistency in entropy estimation and that its choice is related to the series length N. Thus, the threshold and length profiles of the dependency time series are plotted in Figure 9. Although the multiscale entropy is defined on the SampEn basis, the figure also includes ApEn as the worst-case example.

In a multiscale entropy approach, the time series length step-wise decreases with the increased scaling level. Maximal series length of our signals is equal to n = 14,400. If the scaling level is equal to 15, then the minimal series length would be equal to n = 960. Panels a and c of Figure 9 show the threshold profile of ApEn and SampEn for the maximal and for the minimal lengths. Stable results are achieved for threshold r = 0.3 [18].

The length profile is plotted in panels b and d for the typical threshold value r = 0.15 and the chosen threshold value r = 0.3. It can be seen that the results are not consistent for lengths below n = 900, so the choice of 15 scaling levels is justified.

Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14 show the main result of the composite multiscale entropy study of dependency level signals, with the scaling level set to 15, and the threshold level set to r = 0.3. These results are discussed in Section 4.

Figure 10a presents CSME estimated from the dependency level signals of control rats at different ambient temperatures, for two-dimensional signals (SBP and PI interaction, D = 2). Each signal is accompanied by 10 isodistributional surrogates (the estimated entropy is averaged). These results can be compared to the results of the classical entropy analysis shown in Figure 6j. The figure also presents CSME estimated from the artificially generated two-dimensional signals with Gaussian distribution. Figure 10b shows the same entropies but estimated from the three-dimensional signals.

The aim of Figure 11 is to show whether the CMSE estimated from the signals of control rats at different ambient temperatures can distinguish different temporal offsets (DEL) between the SBP and PI. Such an analysis cannot be performed by classical entropy study.

Differences between experimental groups were analyzed by a Mann–Whitney U-test. Statistical significance was considered at p < 0.05.

Figure 12 presents the entropy estimates after the administration of vasopressin antagonists at neutral ambient temperature, for two different SBP-PI temporal offsets (delays), DEL = 0, and DEL = 3. Figure 13 presents the same entropy estimates but at the high ambient temperature. Both Figure 12 and Figure 13 are accompanied by isodistributional surrogate data controls. The purpose of these two figures is to investigate whether CMSE can distinguish the V1a and V2 antagonist administration if compared to the control case (without the drugs). The two-dimensional cases can be compared to the classical entropy study shown in Figure 4.

The purpose of Figure 14 is to check whether the proposed entropy can distinguish the signals after administering the different doses of V1a and V2 antagonists. The same experiments are repeated separating the signals according to the SBP-PI offset (delay), and these results are presented as Supplementary data.

4. Discussion

The first aim of our study is to create a single time series that reflects the level of interaction between the arbitrary number of simultaneously recorded signals. It is accomplished by mathematical tools: copula density captures the level of signal interaction, and the Voronoi decomposition maps the interaction levels into the temporal signal.

Figure 10a shows that the CMS entropy of two-dimensional SBP-PI dependency signal decreases at high ambient temperature. The classical CMSE analysis provided the same result (Figure 4, panel j). The high ambient temperature in this experiment exceeds the boundary set to 29.5 °C [60], inducing heat dissipation in rats such as vasodilatation, evaporation, sweating, panting and affecting the blood vessel circulatory strain [60]. The existence of the dominant component reduces the influence of the other mechanism, the signals and their mutual interactions become less complex, and the entropy decreases. On the other hand, neutral and low ambient temperatures in our experiments are within the normal boundaries, and entropy estimates overlap, both in the proposed and classical entropy estimates—Figure 10a and Figure 4j—respectively.

An administration of V2-500ng significantly decreases the multiscale entropy of a two-dimensional SBP-PI dependency signal, both at neutral (Figure 12a,b) and at high ambient temperature (Figure 13a,b). Additionally, V1a-500 significantly increases the entropy of a two-dimensional SBP-PI dependency signal, signals at high temperatures (Figure 13a,b). These analyses were performed at the particular offsets (delays) between SBP and PI signals. The results correspond to the spectral analysis of the signals after V1a and V2 antagonist administration: it was shown [27] that V2-500 ng administration increases the low-frequency signal component of the SBP signal; the signal becomes smoother, the number of repetitive patterns increases, and the signal becomes more predictive so the entropy decreases. On the other hand, V1a-500 ng increases the high-frequency signal component [27], the signal, and its interactions become more turbulent and the entropy increases.

While spectral analysis separates high and low signal components, a classical cross-entropy observes the signal as a whole. The template matching procedure (cf. Equation (A3)) averages all the temporal SBP and PI positions, so different offsets (delays) of SBP-PI signals cannot be distinguished. Panels k and l of Figure 4 show that classical SBP-PI cross-entropy does slightly decrease after V2-500 ng administration, but without the statistical significance, especially at the neutral temperature (Figure 4k).

It is consistent with the second aim of our contribution: the proposed analysis corresponds to the classical two-dimensional entropy results and with the results of spectral analysis. In cases when the consistency is not statistically significant, the possible explanation is that the classical entropy observes the complete signal, while the proposed method includes signals separated according to the SBP-PI signal offset.

When the body temperature time series is included to create three-dimensional SBP-PI-t_B dependency signals, the entropy decreases and differences caused by ambient temperature are attenuated (Figure 10b). It should be noted (Figure 4g,h,j) that signal t_B is a low-entropy signal. However, the significant entropy decrease induced by V2-500 ng is preserved in three-dimensional signals, but only at the high temperatures and DEL = 0 (Figure 13c).

The proposed method can make a distinction between the different SBP-PI offsets (delays), as shown in Figure 11. Regardless of the ambient temperature, the entropy of the two-dimensional dependency signals at DEL = 5 is significantly lower than at DEL = 0. It corresponds to [29], as the SBP-PI dependency in rats is the greatest for delays of 3, 4, and 5 beats. When three-dimensional signals are observed, the significant decrease of entropy for DEL = 5 occurs at high ambient temperatures only.

The proposed entropy can distinguish the signals after administering the different doses of V1a and V2 antagonists. It is shown in Figure 14. The distinction is statistically significant for the two-dimensional SBP-PI dependency signal, regardless of antagonist and ambient temperature, and for three-dimensional signals after V2 administration. V1a administration induces a statistically significant difference in three-dimensional SBP-PI-t_B dependency signal at the neutral temperature, and only for a couple of scaling levels.

Figure 14 presents the entropy estimates averaged over all the signals. The Supplementary data comprise Figures S1 and S2, where the dosages are separated according to the SBP-PI offset (delay).

The entropy of the artificial two- and three-dimensional Gaussian control signals is shown in Figure 10, but not repeated in the subsequent figures as it is always the same. The entropy of the surrogate data converges towards the Gaussian, but never reaches it: although the surrogate signals can be regarded as streams of i.i.d. random variables, their distribution remains equal to the distribution of the original dependency signals.

5. Conclusions

The signal framework created in this contribution provides a possibility for an easy analysis of signal dependency structures mapped into a single time series. The estimated entropies of two-dimensional dependency signals correspond to the classical cross-entropy and spectral analysis. The method can recognize entropy changes at different temperature levels, at different SBP-PI offsets, at different administered V1a and V2 dosages. The recognition of these differences is not random: the surrogate data analysis destroys the temporal coupling of the observed dependency signals, yielding in all cases entropy estimates that are almost identical and that approach (but do not reach) the entropy of Gaussian time series.

This method can be applied to any multivariate signals. It is necessary to conduct a deep analysis to find the minimal signal length that provides reliable results for the arbitrary number of signals and to check the possibility of applying other analytical tools besides the entropy. A comparative study of the various mapping procedure (Voronoi cells, eigenvalues of the transition matrix, classical multidimensional histograms) should be performed. Further analysis on coupling the cardiovascular data with a body and ambient temperature could reveal more adverse effects, an issue that could be extremely important regarding the climatic changes.