# Comparison of entropy rate measures for the evaluation of time series complexity: simulations and application to heart rate and respiratory variability

Chiara Bara, Riccardo Pernice, Cristina Angela Catania, Mirvana Hilal, Alberto Porta, Anne Humeau-Heurtier, Luca Faes

<figure><img src="image_1.png" /></figure>

## To cite this version:

Chiara Bara, Riccardo Pernice, Cristina Angela Catania, Mirvana Hilal, Alberto Porta, et al.. Comparison of entropy rate measures for the evaluation of time series complexity: simulations and application to heart rate and respiratory variability. Biocybernetics and Biomedical Engineering, 2024, 44 (2), pp.380-392. 10.1016/j.bbe.2024.04.004. hal-04562987

HAL Id: hal-04562987

https://univ-angers.hal.science/hal-04562987

Submitted on 10 Jun 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comparison of entropy rate measures for the evaluation of time series complexity: simulations and application to heart rate and respiratory variability

Chiara Bara, Riccardo Pernice, Cristina Angela Catania, Mirvana Hilal, Alberto Porta, Anne Humeau-Heurtier, Luca Faes

<figure><img src="image_1.png" /></figure>

## To cite this version:

Chiara Bara, Riccardo Pernice, Cristina Angela Catania, Mirvana Hilal, Alberto Porta, et al.. Comparison of entropy rate measures for the evaluation of time series complexity: simulations and application to heart rate and respiratory variability. Biocybernetics and Biomedical Engineering, 2024, 44 (2), pp.380-392. 10.1016/j.bbe.2024.04.004. hal-04562987

HAL Id: hal-04562987

https://univ-angers.hal.science/hal-04562987v1

Submitted on 10 Jun 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Research Paper

# Comparison of entropy rate measures for the evaluation of time series complexity: Simulations and application to heart rate and respiratory variability

Chiara Baràª, Riccardo Perniceª,*, Cristina Angela Cataniaª, Mirvana Hilalᵇ, Alberto Porta c,†, Anne Humeau-Heurtierᵇ, Luca Faesª

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

### Keywords:

- Complexity

- Conditional entropy

- Cardiorespiratory dynamics

- Linear parametric entropy estimation

- Non-linear model-free entropy estimation

## ABSTRACT

Most real-world systems are characterised by dynamics and correlations emerging at multiple time scales, and are therefore referred to as complex systems. In this work, the complexity of time series produced by complex systems was investigated in the frame of information theory computing the entropy rate via the conditional entropy (CE) measure. A comparative investigation of several CE estimators, based on linear parametric and non-linear model-free representations of the process dynamics, was performed considering simulated linear autoregressive (AR) and mixed non-linear deterministic and linear stochastic dynamics processes, as well as physiological time series reflecting short-term cardiorespiratory dynamics. In simulations, the estimated CE values decreased when reducing the system complexity through an increase in the pole radius of the AR process or with the predominance of the deterministic behaviour in the mixed dynamics. In the application to cardiorespiratory dynamics, a reduction in physiological complexity was observed resulting from a regularization of the time series of heart rate and respiratory volume when decreasing the breathing rate. Our results evidence how simple and fast approaches based on linear parametric or permutation-based model-free estimators allow efficient discrimination of complexity changes in the short-term evolution of complex dynamic systems. However, in the presence of non-linear dynamics, the superiority of the more general but computationally expensive nearest-neighbour method is highlighted. These findings have implications for the assessment of complex dynamics both in clinical settings and in physiological monitoring.

## 1. Introduction

Real-world physical and social systems - such as climatic, physiological, biological, financial or evolutionary systems - are broadly characterised as complex [1–3]. The non-trivial and often controversial nature of this concept is demonstrated by the wide variety of applications covered by the term [4]. Over the years, numerous works have attempted to provide a unique definition of complexity in order to identify valid tools for studying the dynamics of complex systems. The idea of complexity is related to phenomena involving multiple interactions of several independent and intrinsically complicated parts [4,5]. The degree of complexity of these systems is highly dependent on the number of factors influencing the individual parts or their interactions, and consequently on the different directions in which the system can evolve. A property commonly ascribed to complex systems

is non-linearity, manifested when the principles of proportionality and superposition of effects that characterise linear systems are lacking [6]. More generally, complex systems are self-organising, with important emergent dynamics in response to small variations in any individual components [4,7,8], and exhibiting self-similarity across multiple time scales and long-range correlations [9].

The human organism can be seen as a complex system where a large number of control mechanisms operate simultaneously to provide homeostatic balance [6,10]. Indeed, the combined effect of multiple coupling and feedback interaction mechanisms between physiological sub-systems results in non-linear dynamics characterised by high variability and long-term self-similarity [11–13]. These complex physiological dynamics reflect the capability of physiological mechanisms to maintain an overall steady state in response to external stimuli or

<footnote>ªDepartment of Engineering, University of Palermo, Palermo, 90128, Italy</footnote>

<footnote>ᵇUniv Angers, LARIS, SFR MATHSTIC, F-49000, Angers, France</footnote>

<footnote>ᶜDepartment of Cardiothoracic, Vascular Anesthesia and Intensive Care, IRCCS Policlinico San Donato, San Donato Milanese, Milan, 20097, Italy</footnote>

<footnote>ᵈDepartment of Biomedical Sciences for Health, University of Milan, Milan, 20133, Italy</footnote>

<footnote>* Corresponding author.
E-mail address: riccardo.pernice@unipa.it (R. Pernice).</footnote>
Research Paper

# Comparison of entropy rate measures for the evaluation of time series complexity: Simulations and application to heart rate and respiratory variability

Chiara Barà<sup>a</sup>, Riccardo Pernice<sup>a,*</sup>, Cristina Angela Catania<sup>a</sup>, Mirvana Hilal<sup>b</sup>, Alberto Porta<sup>c,d</sup>, Anne Humeau-Heurtier<sup>b</sup>, Luca Faes<sup>a</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

### Keywords:

- Complexity

- Conditional entropy

- Cardiorespiratory dynamics

- Linear parametric entropy estimation

- Non-linear model-free entropy estimation

## ABSTRACT

Most real-world systems are characterised by dynamics and correlations emerging at multiple time scales, and are therefore referred to as complex systems. In this work, the complexity of time series produced by complex systems was investigated in the frame of information theory computing the entropy rate via the conditional entropy (CE) measure. A comparative investigation of several CE estimators, based on linear parametric and non-linear model-free representations of the process dynamics, was performed considering simulated linear autoregressive (AR) and mixed non-linear deterministic and linear stochastic dynamics processes, as well as physiological time series reflecting short-term cardiorespiratory dynamics. In simulations, the estimated CE values decreased when reducing the system complexity through an increase in the pole radius of the AR process or with the predominance of the deterministic behaviour in the mixed dynamics. In the application to cardiorespiratory dynamics, a reduction in physiological complexity was observed resulting from a regularization of the time series of heart rate and respiratory volume when decreasing the breathing rate. Our results evidence how simple and fast approaches based on linear parametric or permutation-based model-free estimators allow efficient discrimination of complexity changes in the short-term evolution of complex dynamic systems. However, in the presence of non-linear dynamics, the superiority of the more general but computationally expensive nearest-neighbour method is highlighted. These findings have implications for the assessment of complex dynamics both in clinical settings and in physiological monitoring.

## 1. Introduction

Real-world physical and social systems - such as climatic, physiological, biological, financial or evolutionary systems - are broadly characterised as complex [1–3]. The non-trivial and often controversial nature of this concept is demonstrated by the wide variety of applications covered by the term [4]. Over the years, numerous works have attempted to provide a unique definition of complexity in order to identify valid tools for studying the dynamics of complex systems. The idea of complexity is related to phenomena involving multiple interactions of several independent and intrinsically complicated parts [4,5]. The degree of complexity of these systems is highly dependent on the number of factors influencing the individual parts or their interactions, and consequently on the different directions in which the system can evolve. A property commonly ascribed to complex systems

is non-linearity, manifested when the principles of proportionality and superposition of effects that characterise linear systems are lacking [6]. More generally, complex systems are self-organising, with important emergent dynamics in response to small variations in any individual components [4,7,8], and exhibiting self-similarity across multiple time scales and long-range correlations [9].

The human organism can be seen as a complex system where a large number of control mechanisms operate simultaneously to provide homeostatic balance [6,10]. Indeed, the combined effect of multiple coupling and feedback interaction mechanisms between physiological sub-systems results in non-linear dynamics characterised by high variability and long-term self-similarity [11–13]. These complex physiological dynamics reflect the capability of physiological mechanisms to maintain an overall steady state in response to external stimuli or

<footnote><sup>a</sup> Department of Engineering, University of Palermo, Palermo, 90128, Italy</footnote>

<footnote><sup>b</sup> Univ Angers, LARIS, SFR MATHSTIC, F-49000, Angers, France</footnote>

<footnote><sup>c</sup> Department of Cardiothoracic, Vascular Anesthesia and Intensive Care, IRCCS Policlinico San Donato, San Donato Milanese, Milan, 20097, Italy</footnote>

<footnote><sup>d</sup> Department of Biomedical Sciences for Health, University of Milan, Milan, 20133, Italy</footnote>

<footnote>* Corresponding author.
E-mail address: riccardo.pernice@unipa.it (R. Pernice).</footnote>
intrinsic physiological alterations [14,15]. Several studies show how performing a specific task results in the predominance and synchronisation of specific physiological mechanisms, with decreasing complexity and increasing regularity of physiological dynamics. For example, it was seen how postural stress induced by tilting or standing leads to lower physiological complexity [16–18], or how the practice of physical activity increases the ability to react to physiological variations resulting in a greater complexity of physiological dynamics even in resting conditions [19]. It has also been observed that age-related physiological changes are associated with a decrease in the complexity of both cardiovascular [20,21] and neural [22] dynamics. More regular dynamics can generally be associated with reduced cooperation between physiological systems, indicating the presence of a pathological state [6,23,24], as shown in many investigations conducted on patients with neurological or psychiatric [25–27], cardiac [28], immune [29] or behavioural [30] disorders.

The studies mentioned above rely on diverse meanings of the term “complexity”, which is in fact elusive and can be quantified in different ways [12]. Given the multi-faceted properties of complex systems, tools such as fractal dimension [31], Lyapunov exponent [32], and Lempel–Ziv complexity [33] have been used to study their dynamical behaviour, as these methods are able to capture self-similarities, emergencies, and numerosity of the systems, respectively. In view of the above, a complex system is rarely completely deterministic or chaotic [5], since it exhibits random dynamics that cannot be trivially represented via mathematical models. As a consequence, the stochastic representation of complex dynamics has emerged as a preferential framework for their analysis. A typical approach is to exploit the statistical description of the data generated by the system to study its dynamical behaviour using information-theoretic measures. In broad terms, the evolution of an isolated system results in a large amount of information and increased entropy, whereas the cooperation of systems as parts of a large network results in a more balanced state and more stable entropy [34]. In this context, the idea of complexity can be associated with that of non-regularity and unpredictability of the system dynamics, and can be measured in terms of entropy rate [35,36].

Thanks to their suitability for the description of short-length data with strong stochastic and noisy components, the measures of entropy rate are of great interest for the practical analysis of physiological time series. Depending on the characteristics of the system under investigation, different methods can be used to estimate the entropy rate of a time series. The linear parametric estimator, widely used for its simplicity, is limited to identifying linear dynamics and is suitable only under certain probability distribution assumptions [37]. For this reason, more versatile but computationally expensive methods have been developed to consider the non-linear dynamics often characterising complex physiological systems. These estimation methods can be related to the coarse-graining of embedding spaces, as for the nearest-neighbours [38] or kernel [39] estimators, or to the symbolization of dynamics, as done by binning [40], permutation [41] or slope [42] estimators. The aim of this work is to assess, through a systematic comparative analysis, the ability of the above-mentioned entropy rate estimators to differentiate complex dynamics in synthetic and real physiological dynamical systems. Specifically, we focus first on the analysis of simulated time series generated by a linear autoregressive system with oscillatory components of different amplitude, and by a modulated convex combination of non-linear deterministic dynamics, consisting of a logistic map, and their linear stochastic counterpart, represented by the Iterative Amplitude Adjusted Fourier Transform surrogate. Then, the investigated estimators are compared applying them on physiological time series describing the beat-to-beat variability of the heart period and the respiration amplitude measured during a controlled breathing task.

## 2. Methods

### 2.1. Entropy as a measure of uncertainty

In the framework of information-theory, the Shannon information content is the central concept for the definition of entropy measures [43]. Considering a d-dimensional random variable $V$, this term quantifies the uncertainty about a specific outcome $v = [v_1, ..., v_d]$ of the variable looking at its probabilistic nature. Specifically, the information content associated to $v$ is quantified as $h(v) = -\log p(v)$, where $p(v) = Pr\{V = v\}$ is the probability that $V$ is equal to $v$ [43]. Looking at all the outcomes of the variable, entropy quantifies information as the average uncertainty about such outcomes, i.e., $H(V) = E[h(v)] = -E[\log p(v)]$, where $E[\cdot]$ is the expectation operator [43]. Thus, $H(V) = 0$ if only one outcome occurs with unit probability, whereas if all outcomes have the same probability of occurrence $H(V)$ takes the maximum value.

The probability $p(v)$ can be either continuous or discrete depending on whether $v$ takes values inside a continuous set of values $D_V$ or a finite-size alphabet $A_V$, respectively [44]. If the variable $V$ is continuous its uncertainty is expressed using the differential entropy computed via integration over the continuous domain $D_V$, whereas if $V$ is discrete the entropy is computed via sums extended to the alphabet $A_V$. Entropy measures are typically expressed in nats for continuous variables and in bits for discrete variables, respectively using the natural or the base 2 logarithm [44].

### 2.2. Conditional entropy as measure of complexity

Let us consider a dynamical system $\mathcal{X}$ and the stationary stochastic process $X$ which describes the states of the system over time. Let us further define $X_n$ as the scalar variable obtained sampling $X$ at the present time $n$, $X_n^k = [X_{n-1}, ..., X_{n-k}]$ as the vector sampling the past of the process over the past $k$ lags, and $X_n^- = \lim_{k \to \infty} X_n^k$ as the infinite-dimensional variable representing its whole past history. The separation between present and past states allows to study the dynamical evolution of the information carried by the system. This concept was introduced with the so-called Kolmogorov–Sinai (KS) entropy [45,46], which quantifies the average rate at which the information is produced by the system:

$$H_{KS} = \lim_{n \to \infty} \frac{1}{n} H(X_n, X_n^-), \quad (1)$$

where $H(X_n, X_n^-)$ is the joint entropy of the states of the process $X$ considered up to time $n$. Thus, if the information carried by the system does not increase with time, the average amount of information gained by the system tends towards zero.

Under the assumption of stationarity of the process $X$, the Shannon entropy quantifying the information carried by the system at the current state is a static quantity, i.e., $H(X) = H(X_n) = -E[\log p(x_n)]$, since it provides the same information content at all times [36]. To have insight on the dynamical evolution of the system, it is possible to quantify the new information carried by the present state $X_n$ of the process which cannot be inferred from its past $X_n^-$, i.e., the Conditional Entropy (CE) [36]:

$$H_X = H(X_n|X_n^-) = H(X_n, X_n^-) - H(X_n^-). \quad (2)$$

The conditional entropy defined for a random process as in Eq. (2) is strongly related to the concept of entropy rate. Indeed, the entropy rate (Eq. (1)) and the CE measure (Eq. (2)) asymptotically converge to the same value if the process is stationary [44]. According to these definitions, the level of complexity of the system $\mathcal{X}$ is directly investigated looking at the predictability of the process $X$ [40]. Specifically, if the average information produced by the system at each point in time is high (i.e., high CE), the system generating the signal is characterized by unpredictable, non-deterministic and possibly non-linear interdependencies, and is therefore considered to be overall complex.
intrinsic physiological alterations [14,15]. Several studies show how performing a specific task results in the predominance and synchronisation of specific physiological mechanisms, with decreasing complexity and increasing regularity of physiological dynamics. For example, it was seen how postural stress induced by tilting or standing leads to lower physiological complexity [16–18], or how the practice of physical activity increases the ability to react to physiological variations resulting in a greater complexity of physiological dynamics even in resting conditions [19]. It has also been observed that age-related physiological changes are associated with a decrease in the complexity of both cardiovascular [20,21] and neural [22] dynamics. More regular dynamics can generally be associated with reduced cooperation between physiological systems, indicating the presence of a pathological state [6,23,24], as shown in many investigations conducted on patients with neurological or psychiatric [25–27], cardiac [28], immune [29] or behavioural [30] disorders.

The studies mentioned above rely on diverse meanings of the term “complexity”, which is in fact elusive and can be quantified in different ways [12]. Given the multi-faceted properties of complex systems, tools such as fractal dimension [31], Lyapunov exponent [32], and Lempel–Ziv complexity [33] have been used to study their dynamical behaviour, as these methods are able to capture self-similarities, emergencies, and numerosity of the systems, respectively. In view of the above, a complex system is rarely completely deterministic or chaotic [5], since it exhibits random dynamics that cannot be trivially represented via mathematical models. As a consequence, the stochastic representation of complex dynamics has emerged as a preferential framework for their analysis. A typical approach is to exploit the statistical description of the data generated by the system to study its dynamical behaviour using information-theoretic measures. In broad terms, the evolution of an isolated system results in a large amount of information and increased entropy, whereas the cooperation of systems as parts of a large network results in a more balanced state and more stable entropy [34]. In this context, the idea of complexity can be associated with that of non-regularity and unpredictability of the system dynamics, and can be measured in terms of entropy rate [35,36].

Thanks to their suitability for the description of short-length data with strong stochastic and noisy components, the measures of entropy rate are of great interest for the practical analysis of physiological time series. Depending on the characteristics of the system under investigation, different methods can be used to estimate the entropy rate of a time series. The linear parametric estimator, widely used for its simplicity, is limited to identifying linear dynamics and is suitable only under certain probability distribution assumptions [37]. For this reason, more versatile but computationally expensive methods have been developed to consider the non-linear dynamics often characterising complex physiological systems. These estimation methods can be related to the coarse-graining of embedding spaces, as for the nearest-neighbours [38] or kernel [39] estimators, or to the symbolization of dynamics, as done by binning [40], permutation [41] or slope [42] estimators. The aim of this work is to assess, through a systematic comparative analysis, the ability of the above-mentioned entropy rate estimators to differentiate complex dynamics in synthetic and real physiological dynamical systems. Specifically, we focus first on the analysis of simulated time series generated by a linear autoregressive system with oscillatory components of different amplitude, and by a modulated convex combination of non-linear deterministic dynamics, consisting of a logistic map, and their linear stochastic counterpart, represented by the Iterative Amplitude Adjusted Fourier Transform surrogate. Then, the investigated estimators are compared applying them on physiological time series describing the beat-to-beat variability of the heart period and the respiration amplitude measured during a controlled breathing task.

## 2. Methods

### 2.1. Entropy as a measure of uncertainty

In the framework of information-theory, the Shannon information content is the central concept for the definition of entropy measures [43]. Considering a d-dimensional random variable $V$, this term quantifies the uncertainty about a specific outcome $v = [v_1, ..., v_d]$ of the variable looking at its probabilistic nature. Specifically, the information content associated to $v$ is quantified as $h(v) = -\log p(v)$, where $p(v) = Pr\{V = v\}$ is the probability that $V$ is equal to $v$ [43]. Looking at all the outcomes of the variable, entropy quantifies information as the average uncertainty about such outcomes, i.e., $H(V) = E[h(v)] = -E[\log p(v)]$, where $E[\cdot]$ is the expectation operator [43]. Thus, $H(V) = 0$ if only one outcome occurs with unit probability, whereas if all outcomes have the same probability of occurrence $H(V)$ takes the maximum value.

The probability $p(v)$ can be either continuous or discrete depending on whether $v$ takes values inside a continuous set of values $D_V$ or a finite-size alphabet $A_V$, respectively [44]. If the variable $V$ is continuous its uncertainty is expressed using the differential entropy computed via integration over the continuous domain $D_V$, whereas if $V$ is discrete the entropy is computed via sums extended to the alphabet $A_V$. Entropy measures are typically expressed in nats for continuous variables and in bits for discrete variables, respectively using the natural or the base 2 logarithm [44].

### 2.2. Conditional entropy as measure of complexity

Let us consider a dynamical system $X$ and the stationary stochastic process $X$ which describes the states of the system over time. Let us further define $X_n$ as the scalar variable obtained sampling $X$ at the present time $n$, $X_n^k = [X_{n-1}, ..., X_{n-k}]$ as the vector sampling the past of the process over the past $k$ lags, and $X_n^- = \lim_{k \to \infty} X_n^k$ as the infinite-dimensional variable representing its whole past history. The separation between present and past states allows to study the dynamical evolution of the information carried by the system. This concept was introduced with the so-called Kolmogorov–Sinai (KS) entropy [45,46], which quantifies the average rate at which the information is produced by the system:

$$H_{KS} = \lim_{n \to \infty} \frac{1}{n} H(X_n, X_n^-), \quad (1)$$

where $H(X_n, X_n^-)$ is the joint entropy of the states of the process $X$ considered up to time $n$. Thus, if the information carried by the system does not increase with time, the average amount of information gained by the system tends towards zero.

Under the assumption of stationarity of the process $X$, the Shannon entropy quantifying the information carried by the system at the current state is a static quantity, i.e., $H(X) = H(X_n) = -E[\log p(x_n)]$, since it provides the same information content at all times [36]. To have insight on the dynamical evolution of the system, it is possible to quantify the new information carried by the present state $X_n$ of the process which cannot be inferred from its past $X_n^-$, i.e., the Conditional Entropy (CE) [36]:

$$H_X = H(X_n|X_n^-) = H(X_n, X_n^-) - H(X_n^-). \quad (2)$$

The conditional entropy defined for a random process as in Eq. (2) is strongly related to the concept of entropy rate. Indeed, the entropy rate (Eq. (1)) and the CE measure (Eq. (2)) asymptotically converge to the same value if the process is stationary [44]. According to these definitions, the level of complexity of the system $X$ is directly investigated looking at the predictability of the process $X$ [40]. Specifically, if the average information produced by the system at each point in time is high (i.e., high CE), the system generating the signal is characterized by unpredictable, non-deterministic and possibly non-linear interdependencies, and is therefore considered to be overall complex.
## 2.3. CE estimation approaches

In this section we describe different methods for the practical computation of the entropy rate of a stationary dynamic process X starting from a realization available in the form of the time series $x = \{x_i\}$, $i = 1, \dots, N$. To estimate the CE measure from a single time series of finite length, it is necessary to assume that X is an ergodic Markov process, so as to collect observations over time rather then across realizations and to approximate the infinite-length past history $X_n^-$ with the vector of the m lagged components $X_n^m = [X_{n-1}, \dots, X_{n-m}]$. With these assumptions, the CE in Eq. (2) can be computed as follows:

$$ \hat{H}_X = \hat{H}(X_n, X_n^m) - \hat{H}(X_n^m) = -\mathbb{E}\left[\log \frac{\hat{p}(x_n, X_n^m)}{\hat{p}(X_n^m)}\right], \quad (3) $$

where $\hat{p}(\cdot)$ is an estimate of the probability obtained from the time series x.

Several approaches can be formulated for estimating the probability distributions $p(x_n, x_n^m)$ and $p(x_n^m)$ which, once obtained, can be plugged into Eq. (3) to compute the CE estimates. A first broad categorization can be made between methods assuming a specific probability distribution and using its parameters to derive entropy estimates, i.e., model-based approaches, and methods approximating the probability density directly from the data, i.e., model-free approaches [47,48]. The former approaches are simpler and computationally faster, but lack generality and can fail to describe complex non-linear dynamics. Furthermore, estimators can work either with discrete or continuous variables; since the variables describing physiological dynamics are typically continuous, discrete estimators need to provide a symbolic representation of the values assumed by such variables before attempting to compute entropy measures [40,47,49]. Continuous estimators are more precise but also less data-efficient and more time-consuming, while methods performing discretization can be implemented in faster algorithms but are exposed to accuracy issues related to coarse graining [35]. Both approaches suffer from the curse of dimensionality, i.e., the difficulty of reliably estimating entropies for variables of high dimension (parameter m in Eq. (3)) due to the increasing sparsity of data represented in spaces of growing dimension [35,50]. This issue makes the CE estimates to decrease towards zero even for complex unpredictable time series when the dimension m increases [35,51].

In the following, different CE estimation approaches are taken into account, discussing how the above-mentioned issues can be addressed by a careful choice of the estimation parameters of each method. Later, in Section 5, different applications are discussed to demonstrate the feasibility of successfully exploiting the various estimators in the context of complex physiological dynamics.

The algorithms employed to assess time series complexity have been collected in the univariate Information Decomposition (unID) MAT-LAB toolbox. Freely available on GitHub, the unID toolbox allows to compute conditional entropy and entropy measures describing jointly and separately the dynamics of the past and of the current states of a process, with each of the estimators described in the present work.

### 2.3.1. Linear estimator

The linear estimator is a well-known model-based approach for the estimation of information-theoretic measures based on the assumption that data follow a Gaussian probability distribution [52]. Given a generic d-dimensional Gaussian variable $X_n^d$ with zero-mean, its probability distribution is given by:

$$ p(X_n^d) = \frac{1}{\sqrt{(2\pi)^d |\Sigma_{X_d}|}} e^{-\frac{1}{2} X_n^d \Sigma_{X_d}^{-1} X_n^d^T}, \quad (4) $$

where $\Sigma_{X_d} = \mathbb{E}[X_n^{d^T} X_n^d]$ is the covariance matrix of $X_n^d$. Then, considering the linear parametric representation of a stationary zero-mean Gaussian process X provided by the autoregressive (AR) model of order m, i.e., $X_n = AX_n^m + U_n$ with A a m-dimensional coefficient vector and $U_n$ a zero-mean white Gaussian noise modelling the prediction

error, the entropy of $X_n^m$ and $X_{n+1}^{m+1} = [X_n X_n^m]$ can be computed using Eq. (4) in the definition of differential entropy to get [48]:

$$ H_{\text{lin}}(X_n^m) = \frac{1}{2} \ln(2\pi e |\hat{\Sigma}_{X_m}|), \quad (5) $$

$$ H_{\text{lin}}(X_{n+1}^{m+1}) = \frac{1}{2} \ln(2\pi e |\hat{\Sigma}_{X_{m+1}}|), \quad (6) $$

where $\hat{\Sigma}_{X_m}$ and $\hat{\Sigma}_{X_{m+1}}$ are the estimated covariance matrices of $X_n^m$ and $X_{n+1}^{m+1}$, respectively. Then, the linear estimate of the conditional entropy can be found as [48]:

$$ CE_{\text{lin}} = \frac{1}{2} \ln(2\pi e \hat{\sigma}_U^2), \quad (7) $$

where $\hat{\sigma}_U^2 = \frac{|\hat{\Sigma}_{X^{m+1}}|}{|\Sigma_{X^m}|}$ is the prediction error variance obtained after identifying the regression coefficients of the AR model via the least squares method [53], that coincides with the partial covariance of $X_n^m$ given $X_{n+1}^{m+1}$.

The order m for the AR model can be estimated by the Akaike (AIC) [54] or Bayesian (BIC) [55] Information Criterion to efficiently fit data and avoid overfitting or, conversely, loss of information. Indeed, the parametric representation guarantees high data-efficiency and the possibility to obtain an accurate estimation of entropy measures even from short data sets. Although in many applications the assumption of Gaussianity is satisfied at least in first approximation, when this is not the case the estimator may miss important non-linear structures in the data [56].

### 2.3.2. Nearest-neighbour estimator

The most widely used non-parametric estimator is based on the idea that the probability density around a data point is constant and inversely related to the distance from its nearest samples [38]. Considering the $N-d+1$ observations of the generic d-dimensional continuous random variable $X_n^d$, fixed a number k of neighbour samples, its probability density $p(x_n^d)$ is related to the probability mass of the sphere of diameter $\epsilon_{n,k}^d$ and volume $c_{d,L} \epsilon_{n,k}^d$ centered at the reference point $x_n^d$ and reaching its kth neighbour [40]:

$$ p(x_n^d) = \frac{k}{(N-d) c_{d,L} \epsilon_{n,k}^d}. \quad (8) $$

The estimation of entropy measures by the nearest neighbour metric is not typically performed by directly substituting the derived probability distributions into the equation for entropy, because such distributions are evaluated within embedding spaces of different dimension [38]; actually, several formulations that have been proposed aimed at compensating for the presence of biases. In our work the formalism introduced by Kraskov, Stögbauer, and Grassberger was used [57]. This approach allows to compensate the bias derived from the combination of entropy terms for variables of different dimensions, e.g., $X_n^m$ and $X_{n+1}^{m+1}$, in the estimation of CE, defining the searching distance in the highest dimensional space and then searching for neighbours in the lower dimensional spaces using the same range. The final CE formulation is obtained as [48]:

$$ CE_{knn} = -\psi(k) + \frac{1}{N-m} \sum_{n=1}^{N-m} (\ln \epsilon_{n,k} + \psi(N X_n^m + 1)), \quad (9) $$

where $\psi(\cdot)$ is the digamma function, $\epsilon_{n,k}$ is the maximum-norm distance of $(x_n, x_n^m)$ to its kth neighbour in the $(m+1)$-dimensional space, and $N X_n^m$ the number of points whose distance from $x_n^m$ is lower than $\epsilon_{n,k}/2$ in the m-dimensional space.

The two estimation parameters of the nearest neighbour method are the embedding dimension m, which needs to be set finding a tradeoff between the accuracy of the reconstruction of the history embedding and the curse of dimensionality, and the number of samples k used to evaluate the probability distribution [48], which needs to be set finding a tradeoff between the variance, i.e., the statistical error of the estimates around the mean, and the bias, i.e., the deviation of the estimated mean value from the theoretical one.
## 2.3. CE estimation approaches

In this section we describe different methods for the practical computation of the entropy rate of a stationary dynamic process X starting from a realization available in the form of the time series $x = \{x_i\}$, $i = 1, \dots, N$. To estimate the CE measure from a single time series of finite length, it is necessary to assume that X is an ergodic Markov process, so as to collect observations over time rather then across realizations and to approximate the infinite-length past history $X_n^-$ with the vector of the m lagged components $X_n^m = [X_{n-1}, \dots, X_{n-m}]$. With these assumptions, the CE in Eq. (2) can be computed as follows:

$$ \hat{H}_X = \hat{H}(X_n, X_n^m) - \hat{H}(X_n^m) = -\mathbb{E} \left[ \log \frac{\hat{p}(x_n, x_n^m)}{\hat{p}(x_n^m)} \right], \quad (3) $$

where $\hat{p}(\cdot)$ is an estimate of the probability obtained from the time series x.

Several approaches can be formulated for estimating the probability distributions $p(x_n, x_n^m)$ and $p(x_n^m)$ which, once obtained, can be plugged into Eq. (3) to compute the CE estimates. A first broad categorization can be made between methods assuming a specific probability distribution and using its parameters to derive entropy estimates, i.e., model-based approaches, and methods approximating the probability density directly from the data, i.e., model-free approaches [47,48]. The former approaches are simpler and computationally faster, but lack generality and can fail to describe complex non-linear dynamics. Furthermore, estimators can work either with discrete or continuous variables; since the variables describing physiological dynamics are typically continuous, discrete estimators need to provide a symbolic representation of the values assumed by such variables before attempting to compute entropy measures [40,47,49]. Continuous estimators are more precise but also less data-efficient and more time-consuming, while methods performing discretization can be implemented in faster algorithms but are exposed to accuracy issues related to coarse graining [35]. Both approaches suffer from the curse of dimensionality, i.e., the difficulty of reliably estimating entropies for variables of high dimension (parameter m in Eq. (3)) due to the increasing sparsity of data represented in spaces of growing dimension [35,50]. This issue makes the CE estimates to decrease towards zero even for complex unpredictable time series when the dimension m increases [35,51].

In the following, different CE estimation approaches are taken into account, discussing how the above-mentioned issues can be addressed by a careful choice of the estimation parameters of each method. Later, in Section 5, different applications are discussed to demonstrate the feasibility of successfully exploiting the various estimators in the context of complex physiological dynamics.

The algorithms employed to assess time series complexity have been collected in the univariate Information Decomposition (unID) MATLAB toolbox. Freely available on GitHub, the unID toolbox allows to compute conditional entropy and entropy measures describing jointly and separately the dynamics of the past and of the current states of a process, with each of the estimators described in the present work.

### 2.3.1. Linear estimator

The linear estimator is a well-known model-based approach for the estimation of information-theoretic measures based on the assumption that data follow a Gaussian probability distribution [52]. Given a generic d-dimensional Gaussian variable $X_n^d$ with zero-mean, its probability distribution is given by:

$$ p(X_n^d) = \frac{1}{\sqrt{(2\pi)^d |\Sigma_{X_d}|}} e^{-\frac{1}{2} X_n^d \Sigma_{X_d}^{-1} X_n^d^T}, \quad (4) $$

where $\Sigma_{X_d} = \mathbb{E}[X_n^{d^T} X_n^d]$ is the covariance matrix of $X_n^d$. Then, considering the linear parametric representation of a stationary zero-mean Gaussian process X provided by the autoregressive (AR) model of order m, i.e., $X_n = AX_n^m + U_n$ with A a m-dimensional coefficient vector and $U_n$ a zero-mean white Gaussian noise modelling the prediction

error, the entropy of $X_n^m$ and $X_{n+1}^{m+1} = [X_n X_n^m]$ can be computed using Eq. (4) in the definition of differential entropy to get [48]:

$$ H_{\text{lin}}(X_n^m) = \frac{1}{2} \ln(2\pi e |\hat{\Sigma}_{X_m}|), \quad (5) $$

$$ H_{\text{lin}}(X_{n+1}^{m+1}) = \frac{1}{2} \ln(2\pi e |\hat{\Sigma}_{X_{m+1}}|), \quad (6) $$

where $\hat{\Sigma}_{X_m}$ and $\hat{\Sigma}_{X_{m+1}}$ are the estimated covariance matrices of $X_n^m$ and $X_{n+1}^{m+1}$, respectively. Then, the linear estimate of the conditional entropy can be found as [48]:

$$ CE_{\text{lin}} = \frac{1}{2} \ln(2\pi e \hat{\sigma}_U^2), \quad (7) $$

where $\hat{\sigma}_U^2 = \frac{|\hat{\Sigma}_{X_{m+1}}|}{|\Sigma_{X_m}|}$ is the prediction error variance obtained after identifying the regression coefficients of the AR model via the least squares method [53], that coincides with the partial covariance of $X_n^m$ given $X_{n+1}^{m+1}$.

The order m for the AR model can be estimated by the Akaike (AIC) [54] or Bayesian (BIC) [55] Information Criterion to efficiently fit data and avoid overfitting or, conversely, loss of information. Indeed, the parametric representation guarantees high data-efficiency and the possibility to obtain an accurate estimation of entropy measures even from short data sets. Although in many applications the assumption of Gaussianity is satisfied at least in first approximation, when this is not the case the estimator may miss important non-linear structures in the data [56].

### 2.3.2. Nearest-neighbour estimator

The most widely used non-parametric estimator is based on the idea that the probability density around a data point is constant and inversely related to the distance from its nearest samples [38]. Considering the $N-d+1$ observations of the generic d-dimensional continuous random variable $X_n^d$, fixed a number k of neighbour samples, its probability density $p(x_n^d)$ is related to the probability mass of the sphere of diameter $\epsilon_{n,k}^d$ and volume $c_{d,L}^d$ centered at the reference point $x_n^d$ and reaching its kth neighbour [40]:

$$ p(x_n^d) = \frac{k}{(N-d) c_{d,L}^d \epsilon_{n,k}^d}. \quad (8) $$

The estimation of entropy measures by the nearest neighbour metric is not typically performed by directly substituting the derived probability distributions into the equation for entropy, because such distributions are evaluated within embedding spaces of different dimension [38]; actually, several formulations that have been proposed aimed at compensating for the presence of biases. In our work the formalism introduced by Kraskov, Stögbauer, and Grassberger was used [57]. This approach allows to compensate the bias derived from the combination of entropy terms for variables of different dimensions, e.g., $X_n^m$ and $X_{n+1}^{m+1}$, in the estimation of CE, defining the searching distance in the highest dimensional space and then searching for neighbours in the lower dimensional spaces using the same range. The final CE formulation is obtained as [48]:

$$ CE_{knn} = -\psi(k) + \frac{1}{N-m} \sum_{n=1}^{N-m} (\ln \epsilon_{n,k} + \psi(N X_n^m + 1)), \quad (9) $$

where $\psi(\cdot)$ is the digamma function, $\epsilon_{n,k}$ is the maximum-norm distance of $(x_n, x_n^m)$ to its kth neighbour in the $(m+1)$-dimensional space, and $N X_n^m$ the number of points whose distance from $x_n^m$ is lower than $\epsilon_{n,k}/2$ in the m-dimensional space.

The two estimation parameters of the nearest neighbour method are the embedding dimension m, which needs to be set finding a tradeoff between the accuracy of the reconstruction of the history embedding and the curse of dimensionality, and the number of samples k used to evaluate the probability distribution [48], which needs to be set finding a tradeoff between the variance, i.e., the statistical error of the estimates around the mean, and the bias, i.e., the deviation of the estimated mean value from the theoretical one.
### 2.3.3. Kernel estimator

Another well-known approach for a model-free entropy measures estimation makes use of the kernel method. The probability distribution of the generic variable $X_n^d$ is estimated with this approach giving a different weight to each d-dimensional outcome $x_n^d$ depending on its distance from the reference sample $x_n^d$ [48]:

$$p(x_n^d) = \frac{1}{N-d-1} \sum_{l=1, l \neq i}^{N-d} K(\|x_n^d - x_l^d\|), \quad (10)$$

with $K(\cdot)$ the kernel and $\| \cdot \|$ the norm functions. When adopting the Heaviside kernel function with threshold distance $r$ and the Chebyshev distance (or maximum norm), the kernel-based measure of CE is:

$$CE_{ker} = -\ln \frac{\sum_{n=1}^{N-m} p(x_n^{m+1})}{\sum_{n=1}^{N-m} p(x_n^m)}, \quad (11)$$

where the probabilities $p(x_n^m)$ and $p(x_{n+1}^{m+1})$ are estimated counting the number of observations of $X_n^m$ and $X_{n+1}^{m+1}$ falling within distance $r$ from $x_n^m$ and $x_{n+1}^{m+1}$.

This CE estimator is equivalent to the widely used measure of Sample Entropy [39] and is also related to the previously defined measure of Approximate Entropy [58], where each sample is self-counted in the definition of the probability distributions in Eq. (10) and the sum over all possible patterns is outside the logarithm, unlike in Eq. (11).

The estimation parameters for the kernel method are the embedding dimension $m$ and the width of the Heaviside function $r$. To ensure that the reconstructed space is not too sparsely populated for the evaluation of the probability distribution, which would lead to a bad estimate of the conditional probability, previous work has recommended having at least $10^m$ sample points and large $r$. A common choice for $r$ is to set it as a percentage (between 10 and 30%) of the standard deviation of the analysed time series [39,59].

### 2.3.4. Binning estimator

The most intuitive approach for estimating entropy measures is the so-called binning method, based on the discretization of the continuous random variable representing the process through quantization and on the computation of its probability distribution via the frequentistic approach [40]. Considering a generic stochastic variable $X$ which takes values in the continuous domain $D_X = [x_{min}, x_{max}]$, it is possible to obtain a symbolic representation by associating each sample $x_n$, $n = 1, \dots, N$, with a discrete value $b_n$ belonging to the alphabet $A_{X_h} = \{1, \dots, b\}$. Being $b$ the number of quantization levels (or bins) of amplitude $q = (x_{max} - x_{min})/b$, the discrete variable $X_h$ takes the value $x_{b,n} = i$ when the observation $x_n$ takes a value such that $x_{min} + (i-1)q \le x_n < x_{min}+iq$. Then, the probability of all possible quantized vector variables of length $d$, taking values inside the alphabet $A_{X_h}$ of $b^d$ symbols, is measured as the frequency of occurrence of each possible realization over many observations, and the achieved distribution accounted for the estimation of the entropy term of $X_n^d$. Thus, the binning estimate of conditional entropy, i.e., $CE_{bin}$, is obtained as the difference of the two entropy terms related to the discrete vector variables representing the process dynamics up to time $n$ and $n-1$.

The number of bins $b$ used to discretize the values of the original time series determines the level of coarse graining, so that when $b$ is higher the description of the dynamics of the process is more detailed; this parameter is fixed at 2 when the dynamics are binarized and is typically set to 6 when the CE needs to be estimated on time series of 300 points [40]. Nevertheless, the choice of both the number of bins and the embedding dimension influences the size of the symbolic alphabet, in a way such that it is appropriate to keep $b^{m+1}$ smaller than the time series length $N$ [35,49]. To increase the effectiveness of this estimator when working with a limited number of samples, some approaches have been introduced, such as that of the Corrected Conditional Entropy (CCE) [51], where the embedding dimension $m$ is not a priori defined.

### 2.3.5. Permutation estimator

Among the model-free approaches for entropy estimation, the permutation-based estimator is one of the simplest, most efficient, and robust to noise. This method does not take into account the information carried in the samples of the process in terms of amplitude, but only in their temporal ranking compared to neighbouring points [41]. Unlike the binning method, the discretization performed via permutations works directly on vector variables. Specifically, given a realization $x_n^d = [x_{n,1}, \dots, x_{n,d}]$ of the generic d-dimensional continuous variable $X_n^d$, the rank-ordering procedure allows to obtain the symbolic representation $r_n^d = [r_{n,1}, \dots, r_{n,d}] \in A_{R_X}$ where $r_{n,i} \in \{1, \dots, d\}$ is the rank order of $x_{n,i}$ inside the sequence $x_n^d$ rearranged in ascending order, i.e., $r_{n,i} = 1$ if $x_{n,i} = \min(x_n^d)$ and $r_{n,i} = d$ if $x_{n,i} = \max(x_n^d)$; for two equal components of $x_n^d$, the smallest rank is assigned to the component appearing last. Also in this case, the probability of the $d!$ symbolic vectors belonging to the alphabet $A_{R_X}$ is used to estimate the entropy measure of the variable $X_n^d$. Similarly to the binning approach, the entropy estimates of the variables $X_{n+1}^{m+1}$ and $X_n^m$ are subtracted to quantify the permutation-based CE measure, i.e., $CE_{perm}$. This estimation approach differs from that commonly used in the literature to obtain a permutation-based evaluation of the CE measure [60].

The only estimation parameter used by this approach is the embedding dimension $m$. Even though larger values of $m$ lead to a more detailed description of the dynamics of the process, it is recommended to limit the size of the alphabet used for estimating conditional entropy, i.e., $(m+1)!$, below the length $N$ of the analysed series [49,61].

### 2.3.6. Slope estimator

According to this approach, the discrete vector representing the slope dynamics of the generic d-dimensional variable $X_n^d$ is defined via the $(d-1)$-dimensional variable $S_n^d = [s_{n,1}, \dots, s_{n,d-1}] \in A_{S_X}$ where $s_{n,i} \in \{-2, -1, 0, 1, 2\}$ is the symbol assigned to the difference between two consecutive samples, i.e., $x_{n,i+1} - x_{n,i}$. Specifically, being $\gamma$ and $\delta$ the two parameters indicating respectively the vertical increment and the vicinity to the zero-difference condition when comparing $x_{n,i+1}$ and $x_{n,i}$, the symbols are assigned such that $s_{n,i} = 2$ if $x_{n,i+1} - x_{n,i} > \gamma$, $s_{n,i} = 1$ if $\gamma \ge x_{n,i+1} - x_{n,i} > \delta$, $s_{n,i} = 0$ if $|x_{n,i+1} - x_{n,i}| \le \delta$, $s_{n,i} = -1$ if $-\delta > x_{n,i+1} - x_{n,i} \ge -\gamma$, and $s_{n,i} = -2$ if $x_{n,i+1} - x_{n,i} < -\gamma$. In this case, the slope-based estimate of the entropy of $X_n^d$ is computed as the Shannon entropy of the probability of the $5^{d-1}$ possible realizations of $S_n^d$. For the computation of CE, this approach is applied to the continuous variables $X_{n+1}^{m+1}$ and $X_n^m$ to obtain the corresponding discrete variables $S_{n+1}^{m+1}$ and $S_n^m$, which are then subtracted according to Eq. (3) to estimate the conditional entropy term $CE_{slope}$.

Using this method, apart from the embedding dimension parameter which should be set in order to have $5^m$ lower than $N$, the gradient threshold parameters $\gamma$ and $\delta$ have to be set, even if previous work demonstrates how this approach is fairly stable to variations of these parameters [42].

## 3. Simulation studies

In this section we investigate the behaviour of the different approaches for the estimation of the entropy rate on two simulated dynamic systems. Linear and chaotic dynamics are examined to assess the ability of these methods to capture variations in the complexity of the simulated dynamics as both system and estimation parameters change.

Specifically, the performances of CE estimators are investigated as the complexity and the length of the time series representing system realizations vary. The analyses were carried out first at varying the embedding dimension (or the model order for the linear approach) in the range $m = [2, 3, 4, 6]$ while fixing the other estimation parameters ($k = 10$ for $CE_{knn}$ [62], $r = 0.3$ for $CE_{ker}$ [59], $b = 4$ for $CE_{bin}$ [49], and $\delta = 0.001$ and $\gamma = 1$ for $CE_{slope}$ [42]), and then keeping constant $m = 3$ while varying the number of neighbours in the range $k =$
### 2.3.3. Kernel estimator

Another well-known approach for a model-free entropy measures estimation makes use of the kernel method. The probability distribution of the generic variable $X_n^d$ is estimated with this approach giving a different weight to each d-dimensional outcome $x_n^d$ depending on its distance from the reference sample $x_n^d$ [48]:

$$p(x_n^d) = \frac{1}{N-d-1} \sum_{l=1, l \neq i}^{N-d} K(\|x_n^d - x_l^d\|), \quad (10)$$

with $K(\cdot)$ the kernel and $\| \cdot \|$ the norm functions. When adopting the Heaviside kernel function with threshold distance $r$ and the Chebyshev distance (or maximum norm), the kernel-based measure of CE is:

$$CE_{ker} = -\ln \frac{\sum_{n=1}^{N-m} p(x_n^{m+1})}{\sum_{n=1}^{N-m} p(x_n^m)}, \quad (11)$$

where the probabilities $p(x_n^m)$ and $p(x_{n+1}^{m+1})$ are estimated counting the number of observations of $X_n^m$ and $X_{n+1}^{m+1}$ falling within distance $r$ from $x_n^m$ and $x_{n+1}^{m+1}$.

This CE estimator is equivalent to the widely used measure of Sample Entropy [39] and is also related to the previously defined measure of Approximate Entropy [58], where each sample is self-counted in the definition of the probability distributions in Eq. (10) and the sum over all possible patterns is outside the logarithm, unlike in Eq. (11).

The estimation parameters for the kernel method are the embedding dimension $m$ and the width of the Heaviside function $r$. To ensure that the reconstructed space is not too sparsely populated for the evaluation of the probability distribution, which would lead to a bad estimate of the conditional probability, previous work has recommended having at least $10^m$ sample points and large $r$. A common choice for $r$ is to set it as a percentage (between 10 and 30%) of the standard deviation of the analysed time series [39,59].

### 2.3.4. Binning estimator

The most intuitive approach for estimating entropy measures is the so-called binning method, based on the discretization of the continuous random variable representing the process through quantization and on the computation of its probability distribution via the frequentistic approach [40]. Considering a generic stochastic variable $X$ which takes values in the continuous domain $D_X = [x_{min}, x_{max}]$, it is possible to obtain a symbolic representation by associating each sample $x_n$, $n = 1, \dots, N$, with a discrete value $b_n$ belonging to the alphabet $A_{X_b} = \{1, \dots, b\}$. Being $b$ the number of quantization levels (or bins) of amplitude $q = (x_{max} - x_{min})/b$, the discrete variable $X_b$ takes the value $x_{b,n} = i$ when the observation $x_n$ takes a value such that $x_{min} + (i-1)q \le x_n < x_{min} + iq$. Then, the probability of all possible quantized vector variables of length $d$, taking values inside the alphabet $A_{X_b}$ of $b^d$ symbols, is measured as the frequency of occurrence of each possible realization over many observations, and the achieved distribution accounted for the estimation of the entropy term of $X_n^d$. Thus, the binning estimate of conditional entropy, i.e., $CE_{bin}$, is obtained as the difference of the two entropy terms related to the discrete vector variables representing the process dynamics up to time $n$ and $n-1$.

The number of bins $b$ used to discretize the values of the original time series determines the level of coarse graining, so that when $b$ is higher the description of the dynamics of the process is more detailed; this parameter is fixed at 2 when the dynamics are binarized and is typically set to 6 when the CE needs to be estimated on time series of 300 points [40]. Nevertheless, the choice of both the number of bins and the embedding dimension influences the size of the symbolic alphabet, in a way such that it is appropriate to keep $b^{m+1}$ smaller than the time series length $N$ [35,49]. To increase the effectiveness of this estimator when working with a limited number of samples, some approaches have been introduced, such as that of the Corrected Conditional Entropy (CCE) [51], where the embedding dimension $m$ is not a priori defined.

### 2.3.5. Permutation estimator

Among the model-free approaches for entropy estimation, the permutation-based estimator is one of the simplest, most efficient, and robust to noise. This method does not take into account the information carried in the samples of the process in terms of amplitude, but only in their temporal ranking compared to neighbouring points [41]. Unlike the binning method, the discretization performed via permutations works directly on vector variables. Specifically, given a realization $x_n^d = [x_{n,1}, \dots, x_{n,d}]$ of the generic d-dimensional continuous variable $X_n^d$, the rank-ordering procedure allows to obtain the symbolic representation $r_n^d = [r_{n,1}, \dots, r_{n,d}] \in A_{R_X}$ where $r_{n,i} \in \{1, \dots, d\}$ is the rank order of $x_{n,i}$ inside the sequence $x_n^d$ rearranged in ascending order, i.e., $r_{n,i} = 1$ if $x_{n,i} = \min(x_n^d)$ and $r_{n,i} = d$ if $x_{n,i} = \max(x_n^d)$; for two equal components of $x_n^d$, the smallest rank is assigned to the component appearing last. Also in this case, the probability of the $d!$ symbolic vectors belonging to the alphabet $A_{R_X}$ is used to estimate the entropy measure of the variable $X_n^d$. Similarly to the binning approach, the entropy estimates of the variables $X_{n+1}^{m+1}$ and $X_n^m$ are subtracted to quantify the permutation-based CE measure, i.e., $CE_{perm}$. This estimation approach differs from that commonly used in the literature to obtain a permutation-based evaluation of the CE measure [60].

The only estimation parameter used by this approach is the embedding dimension $m$. Even though larger values of $m$ lead to a more detailed description of the dynamics of the process, it is recommended to limit the size of the alphabet used for estimating conditional entropy, i.e., $(m+1)!$, below the length $N$ of the analysed series [49,61].

### 2.3.6. Slope estimator

According to this approach, the discrete vector representing the slope dynamics of the generic d-dimensional variable $X_n^d$ is defined via the $(d-1)$-dimensional variable $S_n^d = [s_{n,1}, \dots, s_{n,d-1}] \in A_{S_X}$ where $s_{n,i} \in \{-2, -1, 0, 1, 2\}$ is the symbol assigned to the difference between two consecutive samples, i.e., $x_{n,i+1} - x_{n,i}$. Specifically, being $\gamma$ and $\delta$ the two parameters indicating respectively the vertical increment and the vicinity to the zero-difference condition when comparing $x_{n,i+1}$ and $x_{n,i}$, the symbols are assigned such that $s_{n,i} = 2$ if $x_{n,i+1} - x_{n,i} > \gamma$, $s_{n,i} = 1$ if $\gamma \ge x_{n,i+1} - x_{n,i} > \delta$, $s_{n,i} = 0$ if $|x_{n,i+1} - x_{n,i}| \le \delta$, $s_{n,i} = -1$ if $-\delta > x_{n,i+1} - x_{n,i} \ge -\gamma$, and $s_{n,i} = -2$ if $x_{n,i+1} - x_{n,i} < -\gamma$. In this case, the slope-based estimate of the entropy of $X_n^d$ is computed as the Shannon entropy of the probability of the $5^{d-1}$ possible realizations of $S_n^d$. For the computation of CE, this approach is applied to the continuous variables $X_{n+1}^{m+1}$ and $X_n^m$ to obtain the corresponding discrete variables $S_{n+1}^{m+1}$ and $S_n^m$, which are then subtracted according to Eq. (3) to estimate the conditional entropy term $CE_{slope}$.

Using this method, apart from the embedding dimension parameter which should be set in order to have $5^m$ lower than $N$, the gradient threshold parameters $\gamma$ and $\delta$ have to be set, even if previous work demonstrates how this approach is fairly stable to variations of these parameters [42].

## 3. Simulation studies

In this section we investigate the behaviour of the different approaches for the estimation of the entropy rate on two simulated dynamic systems. Linear and chaotic dynamics are examined to assess the ability of these methods to capture variations in the complexity of the simulated dynamics as both system and estimation parameters change.

Specifically, the performances of CE estimators are investigated as the complexity and the length of the time series representing system realizations vary. The analyses were carried out first at varying the embedding dimension (or the model order for the linear approach) in the range $m = [2, 3, 4, 6]$ while fixing the other estimation parameters ($k = 10$ for $CE_{knn}$ [62], $r = 0.3$ for $CE_{ker}$ [59], $b = 4$ for $CE_{bin}$ [49], and $\delta = 0.001$ and $\gamma = 1$ for $CE_{slope}$ [42]), and then keeping constant $m = 3$ while varying the number of neighbours in the range $k =$
[5, 10, 15, 30] for $CE_{knn}$, the Heaviside kernel threshold in the range $r = [0.1, 0.2, 0.3, 0.5]$ for $CE_{ker}$, the number of bins in the range $b = [2, 4, 6, 8]$ for $CE_{bin}$, and the gradient thresholds in the ranges $\gamma = [0.0001, 0.001, 0.01]$ and $\delta = [0.8, 1, 1.2]$ for $CE_{slope}$.

### 3.1. Stationary AR process

An autoregressive (AR) stationary stochastic process $X$ was considered as defined by the linear regression model [48]:

$$X_n = \sum_{k=1}^{m} a_k X_{n-k} + U_n, \quad (12)$$

where $a_k$, $k = 1, \dots, m$, are the linear regression coefficients, $m$ is the order of the AR model, and $U_n$ is a white Gaussian noise with zero-mean and unit variance. In this work, we simulated a third-order AR process featuring an oscillatory dynamic with frequency $f = 0.25$ Hz and amplitude proportional to the modulus $\rho$ of the complex-conjugate poles describing the process in the complex-plane. The amplitude of the stochastic oscillation determines the regularity of the AR process, thus also the complexity of the system dynamics: as the value of the parameter $\rho$ increases, the regularity of the process also increases, shifting from complete unpredictability for $\rho = 0$ to high predictability for $\rho$ close to 1. The true theoretical values of CE were obtained from the known model parameters [63]. However, it is relevant to highlight that theoretical values cannot be used as reference for estimators using symbolic approaches due to the necessity of a corrective term related to the discrete representation of the space [35,64,65].

The analysis was performed on one-hundred realizations of the process obtained for values of $\rho$ ranging from 0 to 0.9 with a step of 0.1, and time series lengths of 300 points. Fig. 1 shows the expected decrease in the CE as $\rho$ increases. This decrease is observed using all estimators, with some exceptions depending on the estimation-specific parameters. We find that both $CE_{lin}$ and $CE_{knn}$ trends agree with the theoretical value of conditional entropy for any $m$ and, for the model-free estimator, also for any $k$, with a general slight increase of the bias with both parameters. The setting of the estimation parameters has a more evident impact on the reliability of the CE estimates obtained using the methods based on discretization. For example, $CE_{ker}$ could not be computed with a high embedding dimension $m$ or similarity threshold $r$ due to the increasing sparsity of the explored spaces, observed particularly for higher system complexity. Moreover, the $CE_{bin}$, $CE_{perm}$, and $CE_{slope}$ approaches fail to reflect the decrease in complexity with increasing $\rho$ when the alphabet dimension becomes larger than the length of the analysed time series, i.e., when $m > 4$ and $b = 4$ or $m = 3$ and $b > 6$ for the binning method (alphabet size of 1024 and 1296, respectively), and $m = 6$ for the permutation method (alphabet size of 5040). Regarding $CE_{slope}$, although the alphabet size is still larger than $N = 300$ for $m = 4$ (i.e., 625), the expected trend is still observed; in addition, the parameters $\delta$ and $\gamma$ have less influence on the trend of the measure. Furthermore, for $m = 2$ the permutation- and slope-based methods do not capture system dynamics variations as the simulation involves a third-order process.

The ability of the estimators to capture variations in magnitude of the system complexity was then studied quantifying the percentage variations obtained for the CE measure when varying the parameter $\rho$ between pairs of consecutive values taken in the set $\rho = [0, 0.3, 0.5, 0.7, 0.8, 0.9]$, as well as computing the statistical significance of such variations through a paired Student's t-test. The results reported in Fig. 2 show how variations of the parameter $\rho$ associated with higher process complexity, i.e., 0.5 vs. 0.7, 0.7 vs. 0.8 and 0.8 vs. 0.9, lead to changes in CE that are detected as significant by all estimators. Since these are differential entropy metrics, the values obtained with the linear and the nearest neighbour estimators behave more similarly to the true theoretical values, with also lower $p$-values than the other approaches. Otherwise, for lower complexity variations, i.e., 0 vs. 0.3 and 0.3 vs. 0.5, not all estimators are able to reveal significance.

Finally, the influence of the time series length on the reliability of the complexity estimates was investigated computing the CE for increasing values of $N$ in the range [50, 100, 200, 300, 500, 700, 1000], after fixing the amplitude of the stochastic oscillation to $\rho = 0.8$. From the results reported in Fig. 3, for each estimator, increasing the time series length produced the expected effect that the estimates of the conditional entropy stabilize and present lower variance. Moreover, we observe that the rate of convergence to stable values at increasing $N$ is faster when the embedding dimension is the one imposed for the simulated system ($m = 3$) or lower, and decreases at increasing $m$. In particular, the estimation of $CE_{lin}$ and $CE_{knn}$ tends towards the theoretical value of conditional entropy. The convergence values of $CE_{ker}$ and $CE_{bin}$ are also strongly dependent on the parameters $r$ and $b$, respectively.

### 3.2. Mixed non-linear deterministic and linear stochastic dynamics process

A combination of linear stochastic and non-linear deterministic dynamics was used to test the capability of the CE estimators to capture complexity changes related to the randomness of the system. Introduced by May to describe the evolutionary dynamics of population or phenomena [66], the logistic map was used to simulate a chaotic deterministic system. According to this model, the temporal evolution of a logistic process $Y$ is described by the non-linear finite-difference equation [66]:

$$Y_n = c Y_{n-1}(1 - Y_{n-1}), \quad (13)$$

where $c$ is the growth rate of the chaotic dynamics. Specifically, the internal dynamics of the system become increasingly non-trivial as the parameter $c$ rises up, turning into chaotic in the range $c = [3.57, 4)$. After generating a chaotic dynamic imposing $c = 3.9$, a process $Z$ with the same spectrum and linear correlation was generated by using the Iterative Amplitude Adjusted Fourier Transform (IAAFT) technique [67], and then the two were combined as follows:

$$X_n = \lambda Y_n + (1 - \lambda) Z_n, \quad (14)$$

The parameter $\lambda$ regulates the convex combination of the two opposite types of dynamics: as the value of $\lambda$ increases, the system shifts from purely linear stochastic dynamics (obtained with $\lambda = 0$) to non-linear deterministic (obtained with $\lambda = 1$).

One-hundred 300-sample realizations of the process described by Eq. (14) were generated, by varying the parameter $\lambda$ from 0 to 1 in steps of 0.1. Fig. 4 shows that the estimated CE values follow an overall tendency to decrease as the value of $\lambda$ increases, after being almost constant up to $\lambda = 0.5$. Thus, high levels of stochasticity result in highly complex dynamics, while the complexity decreases as the deterministic component of the system becomes more dominant. Interesting observations arise from the trend of $CE_{lin}$, which remains constant at varying $\lambda$. This trend reflects the fact that, since the linear temporal correlation is preserved by the IAFFT surrogate process, the complexity related to the linear dynamics remains unchanged. Observations similar to the first simulation can be made regarding the influence of the estimation parameters on the CE trends for all the methods. Parameter settings leading to work in high dimensional spaces are, also in this case, unable to detect the expected CE changes. Nonetheless, in this case, the lowest embedding dimension is also able to track the expected CE trends using discretization-based approaches, as the model of the simulated process involves one-lag dynamics.

Fig. 5 depicts the percentage variations in the CE estimates obtained by changing the parameter $\lambda$ between pairs of values in the range [0, 0.1, 0.3, 0.5, 0.7, 0.9, 1]. A statistically significant decrease in the estimated CE is reported for all the estimation approaches (except for the linear parametric one) for values of $\lambda$ higher than 0.5. Even with higher $p$-values, slight complexity variations are deemed statistically significant for values of $\lambda$ lower than 0.5 for some estimators, i.e., the kernel, binning and slope based. The nearest-neighbour approach seems
[5, 10, 15, 30] for CE$_{knn}$, the Heaviside kernel threshold in the range r = [0.1, 0.2, 0.3, 0.5] for CE$_{ker}$, the number of bins in the range b = [2, 4, 6, 8] for CE$_{bin}$, and the gradient thresholds in the ranges γ = [0.0001, 0.001, 0.01] and δ = [0.8, 1, 1.2] for CE$_{slope}$.

### 3.1. Stationary AR process

An autoregressive (AR) stationary stochastic process X was considered as defined by the linear regression model [48]:

$$X_n = \sum_{k=1}^{m} a_k X_{n-k} + U_n, \qquad (12)$$

where $a_k$, $k = 1, \dots, m$, are the linear regression coefficients, $m$ is the order of the AR model, and $U_n$ is a white Gaussian noise with zero-mean and unit variance. In this work, we simulated a third-order AR process featuring an oscillatory dynamic with frequency $f = 0.25$ Hz and amplitude proportional to the modulus $\rho$ of the complex-conjugate poles describing the process in the complex-plane. The amplitude of the stochastic oscillation determines the regularity of the AR process, thus also the complexity of the system dynamics: as the value of the parameter $\rho$ increases, the regularity of the process also increases, shifting from complete unpredictability for $\rho = 0$ to high predictability for $\rho$ close to 1. The true theoretical values of CE were obtained from the known model parameters [63]. However, it is relevant to highlight that theoretical values cannot be used as reference for estimators using symbolic approaches due to the necessity of a corrective term related to the discrete representation of the space [35,64,65].

The analysis was performed on one-hundred realizations of the process obtained for values of $\rho$ ranging from 0 to 0.9 with a step of 0.1, and time series lengths of 300 points. Fig. 1 shows the expected decrease in the CE as $\rho$ increases. This decrease is observed using all estimators, with some exceptions depending on the estimation-specific parameters. We find that both CE$_{lin}$ and CE$_{knn}$ trends agree with the theoretical value of conditional entropy for any m and, for the model-free estimator, also for any k, with a general slight increase of the bias with both parameters. The setting of the estimation parameters has a more evident impact on the reliability of the CE estimates obtained using the methods based on discretization. For example, CE$_{ker}$ could not be computed with a high embedding dimension m or similarity threshold r due to the increasing sparsity of the explored spaces, observed particularly for higher system complexity. Moreover, the CE$_{bin}$, CE$_{perm}$, and CE$_{slope}$ approaches fail to reflect the decrease in complexity with increasing $\rho$ when the alphabet dimension becomes larger than the length of the analysed time series, i.e., when $m > 4$ and $b = 4$ or $m = 3$ and $b > 6$ for the binning method (alphabet size of 1024 and 1296, respectively), and $m = 6$ for the permutation method (alphabet size of 5040). Regarding CE$_{slope}$, although the alphabet size is still larger than $N = 300$ for $m = 4$ (i.e., 625), the expected trend is still observed; in addition, the parameters $\delta$ and $\gamma$ have less influence on the trend of the measure. Furthermore, for $m = 2$ the permutation- and slope-based methods do not capture system dynamics variations as the simulation involves a third-order process.

The ability of the estimators to capture variations in magnitude of the system complexity was then studied quantifying the percentage variations obtained for the CE measure when varying the parameter $\rho$ between pairs of consecutive values taken in the set $\rho = [0, 0.3, 0.5, 0.7, 0.8, 0.9]$, as well as computing the statistical significance of such variations through a paired Student's t-test. The results reported in Fig. 2 show how variations of the parameter $\rho$ associated with higher process complexity, i.e., 0.5 vs. 0.7, 0.7 vs. 0.8 and 0.8 vs. 0.9, lead to changes in CE that are detected as significant by all estimators. Since these are differential entropy metrics, the values obtained with the linear and the nearest neighbour estimators behave more similarly to the true theoretical values, with also lower p-values than the other approaches. Otherwise, for lower complexity variations, i.e., 0 vs. 0.3 and 0.3 vs. 0.5, not all estimators are able to reveal significance.

Finally, the influence of the time series length on the reliability of the complexity estimates was investigated computing the CE for increasing values of N in the range [50, 100, 200, 300, 500, 700, 1000], after fixing the amplitude of the stochastic oscillation to $\rho = 0.8$. From the results reported in Fig. 3, for each estimator, increasing the time series length produced the expected effect that the estimates of the conditional entropy stabilize and present lower variance. Moreover, we observe that the rate of convergence to stable values at increasing N is faster when the embedding dimension is the one imposed for the simulated system ($m = 3$) or lower, and decreases at increasing $m$. In particular, the estimation of CE$_{lin}$ and CE$_{knn}$ tends towards the theoretical value of conditional entropy. The convergence values of CE$_{ker}$ and CE$_{bin}$ are also strongly dependent on the parameters r and b, respectively.

### 3.2. Mixed non-linear deterministic and linear stochastic dynamics process

A combination of linear stochastic and non-linear deterministic dynamics was used to test the capability of the CE estimators to capture complexity changes related to the randomness of the system. Introduced by May to describe the evolutionary dynamics of population or phenomena [66], the logistic map was used to simulate a chaotic deterministic system. According to this model, the temporal evolution of a logistic process Y is described by the non-linear finite-difference equation [66]:

$$Y_n = c Y_{n-1}(1 - Y_{n-1}), \qquad (13)$$

where c is the growth rate of the chaotic dynamics. Specifically, the internal dynamics of the system become increasingly non-trivial as the parameter c rises up, turning into chaotic in the range $c = [3.57, 4)$. After generating a chaotic dynamic imposing $c = 3.9$, a process Z with the same spectrum and linear correlation was generated by using the Iterative Amplitude Adjusted Fourier Transform (IAAFT) technique [67], and then the two were combined as follows:

$$X_n = \lambda Y_n + (1 - \lambda) Z_n. \qquad (14)$$

The parameter $\lambda$ regulates the convex combination of the two opposite types of dynamics: as the value of $\lambda$ increases, the system shifts from purely linear stochastic dynamics (obtained with $\lambda = 0$) to non-linear deterministic (obtained with $\lambda = 1$).

One-hundred 300-sample realizations of the process described by Eq. (14) were generated, by varying the parameter $\lambda$ from 0 to 1 in steps of 0.1. Fig. 4 shows that the estimated CE values follow an overall tendency to decrease as the value of $\lambda$ increases, after being almost constant up to $\lambda = 0.5$. Thus, high levels of stochasticity result in highly complex dynamics, while the complexity decreases as the deterministic component of the system becomes more dominant. Interesting observations arise from the trend of CE$_{lin}$, which remains constant at varying $\lambda$. This trend reflects the fact that, since the linear temporal correlation is preserved by the IAFFT surrogate process, the complexity related to the linear dynamics remains unchanged. Observations similar to the first simulation can be made regarding the influence of the estimation parameters on the CE trends for all the methods. Parameter settings leading to work in high dimensional spaces are, also in this case, unable to detect the expected CE changes. Nonetheless, in this case, the lowest embedding dimension is also able to track the expected CE trends using discretization-based approaches, as the model of the simulated process involves one-lag dynamics.

Fig. 5 depicts the percentage variations in the CE estimates obtained by changing the parameter $\lambda$ between pairs of values in the range [0, 0.1, 0.3, 0.5, 0.7, 0.9, 1]. A statistically significant decrease in the estimated CE is reported for all the estimation approaches (except for the linear parametric one) for values of $\lambda$ higher than 0.5. Even with higher p-values, slight complexity variations are deemed statistically significant for values of $\lambda$ lower than 0.5 for some estimators, i.e., the kernel, binning and slope based. The nearest-neighbour approach seems
<figure><img src="image_4.png" /><figcaption>Fig. 1. Performance of the six conditional entropy estimators as a function of the estimation parameters (left column: embedding dimension *m*; right column: estimator-specific parameter) for a stationary AR process featuring a stochastic oscillation, studied at varying the amplitude *ρ* of the stochastic oscillation. Errorbar plots (mean ± std) depict the distributions of CE computed over 100 realizations of length *N* = 300 using the (a) linear (CE<sub>lin</sub>), (b) nearest neighbour (CE<sub>knn</sub>), (c) kernel (CE<sub>ker</sub>), (d) binning (CE<sub>bin</sub>), (e) permutation (CE<sub>perm</sub>), and (f) slope (CE<sub>slope</sub>) approaches. For the linear and nearest neighbour estimates, a black line representing the theoretical values of the conditional entropy is also reported.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. Performance of the six conditional entropy measures in detecting the magnitude of complexity variations in the dynamics of a stationary AR process with fixed frequency *f* = 0.25 Hz. Errorbar plots depict the percentage variation (mean + std.dev. across 100 realizations of length *N* = 300 points) of the CE computed by increasing the amplitude of the AR stochastic oscillation taking consecutive values in the range *ρ* = [0, 0.3, 0.5, 0.7, 0.8, 0.9]; specific indications about the pair of *ρ* values associated with the compared distributions are reported on the x-axis. The estimation parameters are set to *m* = 3 for all methods, and *k* = 10 for CE<sub>knn</sub>, *r* = 0.3 for CE<sub>ker</sub>, *b* = 4 for CE<sub>bin</sub>, and *δ* = 0.001 and *γ* = 1 for CE<sub>slope</sub>. The horizontal black lines represent the true theoretical percentage variation of the conditional entropy measure. For each estimation approach, asterisk denotes a statistically significant difference between the distributions of the CE obtained for the pair of values selected for *ρ* (paired Student's *t*-test, *p* < 0.05).</figcaption></figure>

to be the most sensitive in revealing complexity changes, as percentage
differences are evaluated using very low initial absolute CE values.

In addition, also for these simulated dynamics, the influence of the
time series length was investigated in the complexity index reliability
once a value of λ = 0.5 was fixed, reporting similar results to those

achieved for the first simulation, as shown in Fig. 6. Indeed, it is
observed that the estimated values of CE tend to converge faster using
smaller values of the embedding parameter *m*, and report a strong
dependence on *r* with regard to the kernel approach and on *b* for the
binning one.

3.3. *Computational costs*

Further analyses were conducted to compare the computational
costs of the investigated estimation approaches. Specifically, compu-
tational complexity and times have been evaluated according to the
implementation of the algorithms of the unID toolbox.

Table 1 reports the averaged single-run computational times taken to estimate CE and the other information dynamic measures by each approach, expressed as the average over one-hundred realizations of the AR process (setting parameters: *f* = 0.25 Hz and *ρ* = 0.8) and of the mixed non-linear deterministic and linear stochastic dynamics process (setting parameter: *c* = 3.9 and *λ* = 0.5), at varying time series length *N* = [300, 500, 1000]. For all methods, the embedding parameter was fixed to *m* = 3, setting *k* = 10 for *CE*<sub>knn</sub>, *r* = 0.3 for *CE*<sub>ker</sub>, *b* = 4 for *CE*<sub>bin</sub>, and *δ* = 0.001 and *γ* = 1 for *CE*<sub>slope</sub>. Computational times were calculated using computer equipped with an Intel Core i5-6200U CPU (2.30 GHz), 8 GB RAM, Windows 10 Home, MATLAB R2020a. The percentage variations of the computational time increasing the time series length from 300 to 1000 samples are also reported in the table. The computational times measured for the linear approach are sharply lower compared to nearest-neighbour and kernel methods, with differences of two or even three orders of magnitude. Discretization-based methods require a computational time higher by one order of magnitude than the linear approach, with the binning method being the most demanding one. As expected, for all estimation approaches the computational time increases with the time series length, most strongly for the nearest-neighbour and the kernel estimators.

These results are partially related to the computational complexity
of the estimation algorithms, which was investigated by using the O
<figure><img src="image_4.png" /><figcaption>Fig. 1. Performance of the six conditional entropy estimators as a function of the estimation parameters (left column: embedding dimension *m*; right column: estimator-specific parameter) for a stationary AR process featuring a stochastic oscillation, studied at varying the amplitude *ρ* of the stochastic oscillation. Errorbar plots (mean ± std) depict the distributions of CE computed over 100 realizations of length *N* = 300 using the (a) linear (CE<sub>lin</sub>), (b) nearest neighbour (CE<sub>knn</sub>), (c) kernel (CE<sub>ker</sub>), (d) binning (CE<sub>bin</sub>), (e) permutation (CE<sub>perm</sub>), and (f) slope (CE<sub>slope</sub>) approaches. For the linear and nearest neighbour estimates, a black line representing the theoretical values of the conditional entropy is also reported.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. Performance of the six conditional entropy measures in detecting the magnitude of complexity variations in the dynamics of a stationary AR process with fixed frequency *f* = 0.25 Hz. Errorbar plots depict the percentage variation (mean + std.dev. across 100 realizations of length *N* = 300 points) of the CE computed by increasing the amplitude of the AR stochastic oscillation taking consecutive values in the range *ρ* = [0, 0.3, 0.5, 0.7, 0.8, 0.9]; specific indications about the pair of *ρ* values associated with the compared distributions are reported on the x-axis. The estimation parameters are set to *m* = 3 for all methods, and *k* = 10 for CE<sub>knn</sub>, *r* = 0.3 for CE<sub>ker</sub>, *b* = 4 for CE<sub>bin</sub>, and *δ* = 0.001 and *γ* = 1 for CE<sub>slope</sub>. The horizontal black lines represent the true theoretical percentage variation of the conditional entropy measure. For each estimation approach, asterisk denotes a statistically significant difference between the distributions of the CE obtained for the pair of values selected for *ρ* (paired Student's *t*-test, *p* < 0.05).</figcaption></figure>

to be the most sensitive in revealing complexity changes, as percentage
differences are evaluated using very low initial absolute CE values.

In addition, also for these simulated dynamics, the influence of the
time series length was investigated in the complexity index reliability
once a value of λ = 0.5 was fixed, reporting similar results to those

achieved for the first simulation, as shown in Fig. 6. Indeed, it is
observed that the estimated values of CE tend to converge faster using
smaller values of the embedding parameter *m*, and report a strong
dependence on *r* with regard to the kernel approach and on *b* for the
binning one.

3.3. *Computational costs*

Further analyses were conducted to compare the computational
costs of the investigated estimation approaches. Specifically, compu-
tational complexity and times have been evaluated according to the
implementation of the algorithms of the unID toolbox.

Table 1 reports the averaged single-run computational times taken to estimate CE and the other information dynamic measures by each approach, expressed as the average over one-hundred realizations of the AR process (setting parameters: *f* = 0.25 Hz and *ρ* = 0.8) and of the mixed non-linear deterministic and linear stochastic dynamics process (setting parameter: *c* = 3.9 and *λ* = 0.5), at varying time series length *N* = [300, 500, 1000]. For all methods, the embedding parameter was fixed to *m* = 3, setting *k* = 10 for *CE*<sub>knn</sub>, *r* = 0.3 for *CE*<sub>ker</sub>, *b* = 4 for *CE*<sub>bin</sub>, and *δ* = 0.001 and *γ* = 1 for *CE*<sub>slope</sub>. Computational times were calculated using computer equipped with an Intel Core i5-6200U CPU (2.30 GHz), 8 GB RAM, Windows 10 Home, MATLAB R2020a. The percentage variations of the computational time increasing the time series length from 300 to 1000 samples are also reported in the table. The computational times measured for the linear approach are sharply lower compared to nearest-neighbour and kernel methods, with differences of two or even three orders of magnitude. Discretization-based methods require a computational time higher by one order of magnitude than the linear approach, with the binning method being the most demanding one. As expected, for all estimation approaches the computational time increases with the time series length, most strongly for the nearest-neighbour and the kernel estimators.

These results are partially related to the computational complexity
of the estimation algorithms, which was investigated by using the O
<figure><img src="image_6.png" /><figcaption>Fig. 3. Performance of the six conditional entropy estimators as a function of the estimation parameters (left column: embedding dimension *m*; right column: estimator-specific parameter) for time series of length *N* representing a stationary AR process featured by a stochastic oscillation of frequency *f* = 0.25 Hz and amplitude *ρ* = 0.8. Errorbar plots (mean ± std) depict the distributions of CE computed over 100 realizations using the (a) linear (CE<sub>lin</sub>), (b) nearest neighbour (CE<sub>knn</sub>), (c) kernel (CE<sub>ker</sub>), (d) binning (CE<sub>bin</sub>), (e) permutation (CE<sub>perm</sub>), and (f) slope (CE<sub>slope</sub>) approaches. For the linear and nearest neighbour estimates, a black line representing the theoretical values of the conditional entropy is also reported.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Performance of the six conditional entropy estimators as a function of the estimation parameters (left column: embedding dimension *m*; right column: estimator-specific parameter) for a mixed non-linear deterministic and linear stochastic process studied at varying the *λ* parameter. Errorbar plots depict the distributions of CE computed over 100 realizations of length *N* = 300 using the (a) linear (CE<sub>lin</sub>), (b) nearest neighbour (CE<sub>knn</sub>), (c) kernel (CE<sub>ker</sub>), (d) binning (CE<sub>bin</sub>), (e) permutation (CE<sub>perm</sub>), and (f) slope (CE<sub>slope</sub>) approaches.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 3. Performance of the six conditional entropy estimators as a function of the estimation parameters (left column: embedding dimension *m*; right column: estimator-specific parameter) for time series of length *N* representing a stationary AR process featured by a stochastic oscillation of frequency *f* = 0.25 Hz and amplitude *ρ* = 0.8. Errorbar plots (mean ± std) depict the distributions of CE computed over 100 realizations using the (a) linear (CE<sub>lin</sub>), (b) nearest neighbour (CE<sub>knn</sub>), (c) kernel (CE<sub>ker</sub>), (d) binning (CE<sub>bin</sub>), (e) permutation (CE<sub>perm</sub>), and (f) slope (CE<sub>slope</sub>) approaches. For the linear and nearest neighbour estimates, a black line representing the theoretical values of the conditional entropy is also reported.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Performance of the six conditional entropy estimators as a function of the estimation parameters (left column: embedding dimension *m*; right column: estimator-specific parameter) for a mixed non-linear deterministic and linear stochastic process studied at varying the *λ* parameter. Errorbar plots depict the distributions of CE computed over 100 realizations of length *N* = 300 using the (a) linear (CE<sub>lin</sub>), (b) nearest neighbour (CE<sub>knn</sub>), (c) kernel (CE<sub>ker</sub>), (d) binning (CE<sub>bin</sub>), (e) permutation (CE<sub>perm</sub>), and (f) slope (CE<sub>slope</sub>) approaches.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 5. Performance of the six conditional entropy measures in detecting the magnitude of complexity variations in the dynamics of a mixed non-linear deterministic and linear stochastic process. Errorbar plots depict the percentage variation of the CE computed by increasing the parameter regulating the convex combination of the two dynamics taking consecutive values in the range $\lambda = [0,0.1,0.3,0.5,0.7,0.9,1]$; specific indications about the pair of $\lambda$ values associated with the compared distributions are reported on the x-axis. The estimation parameters are set to $m = 3$ for all methods, and $k = 10$ for $CE_{knn}$, $r = 0.3$ for $CE_{ker}$, $b = 4$ for $CE_{bin}$, and $\delta = 0.001$ and $\gamma = 1$ for $CE_{slope}$. For each estimation approach, asterisk denotes a statistically significant difference between the distributions of the CE obtained for the pair of values selected for $\lambda$ (paired Student's t-test, $p < 0.05$).</figcaption></figure>

<caption>Table 1</caption>

Single-run computational times (expressed in milliseconds) for the estimation of CE measure and the other information dynamic measures, averaged on one-hundred realizations of both autoregressive and mixed non-linear deterministic and linear stochastic dynamics processes. The estimation parameters are set to $m = 3$ for all methods, and $k = 10$ for $CE_{knn}$, $r = 0.3$ for $CE_{ker}$, $b = 4$ for $CE_{bin}$, and $\delta = 0.001$ and $\gamma = 1$ for $CE_{slope}$. In the last row, the percentage variations of computational time when increasing the time series length from $N = 300$ to $N = 1000$ are reported.

<table><thead><tr><th rowspan="3">Time series<br>length N</th><th colspan="6">Estimation approaches</th></tr><tr><th colspan="6">Computational times [ms]</th></tr><tr><th>Lin</th><th>Knn</th><th>Ker</th><th>Bin</th><th>Perm</th><th>Slope</th></tr></thead><tbody><tr><td>300</td><td>0.19</td><td>12.37</td><td>9.18</td><td>0.72</td><td>0.53</td><td>0.71</td></tr><tr><td>500</td><td>0.27</td><td>33.85</td><td>21.71</td><td>1</td><td>0.69</td><td>0.73</td></tr><tr><td>1000</td><td>0.31</td><td>121.33</td><td>80.53</td><td>2.43</td><td>1.74</td><td>1.48</td></tr><tr><td colspan="7">From N = 300 to N = 1000</td></tr><tr><td colspan="7">Time increment [%]</td></tr><tr><td colspan="7">59 881 777 238 225 109</td></tr></tbody></table>

Big notation. Specifically, the linear parametric and the discretization-based estimation approaches are of $O(N)$ and $O(N \log(N))$ complexity, respectively. Otherwise, the methods based on computing distances are more demanding, with computational complexities of $O(N^2 \log(N))$ and $O(N^2)$, for the nearest-neighbour and the kernel approaches, respectively. Considering that reliable estimates of entropy measures are obtained for large values of the time series length $N$ and that this is quite larger than the embedding size $m$, i.e., $N \gg m$, the only parameter affecting the computational complexity of the employed algorithms is indeed $N$. On the other hand, the influence of estimator-specific parameters, i.e., $k$ for nearest-neighbour, $r$ for kernel and $b$ for binning, can be neglected in the asymptotic complexity behaviours, even if they in practice cause variations in computational time.

## 4. Application to physiological time series

In this section, conditional entropy estimation approaches are investigated in the context of short-term physiological variability analysis, where heart period and respiratory time series are studied to assess variations in the complexity of cardiorespiratory dynamics during spontaneous and controlled breathing.

### 4.1. Experimental protocol and time series extraction

The investigated physiological signals include electrocardiographic (ECG) and respiratory flow recordings (sampling frequency: 300 Hz) acquired on 19 young healthy subjects (8 males; age range 27–35 years) in resting supine position during different breathing patterns. The analysed experimental conditions were spontaneous breathing (SB) and paced breathing at 10, 15, and 20 breaths/min (C10, C15, and C20, respectively) [18,68,69]. All procedures were approved by the ethical committee of the “L. Sacco” Hospital (Milan, Italy) and of the Department of Technologies for Health, University of Milan (Italy) [68].

The analysis was carried out on the sequence of heart periods (RR intervals) extracted as the time intervals between two consecutive ECG R peaks and on the breathing (RESP) time series obtained by sampling the respiratory flow signal at the onset of every heart period. For each subject and condition, stationary windows of 256 samples were selected synchronously for the two time series. Before estimating CE measures, all time series were preprocessed applying a high-pass AR filter (zero phase and cut-off frequency 0.0156 Hz) and by removing and interpolating samples differing more than three standard deviations from the mean value. Finally, the series were normalized to zero mean and unit variance. Further details on data acquisition and pre-processing can be found in [68].

### 4.2. Data analysis

The CE estimation approaches described in Section 2.3 were applied setting the estimation parameters on the basis of the results of the simulation studies and considering the most widely employed values in literature [40,42,49]. Specifically, the Bayesian Information Criterion (BIC) was employed to define the optimal order of the AR model used to approximate the past dynamics of the investigated process (maximum order fixed to 10), while an embedding dimension $m = 2$ was set to apply nearest-neighbour, kernel and binning approaches. The other estimation parameters set to implement these methods were $k = 10$, $r = 0.3$, and $b = 6$, respectively. The embedding dimension was set to $m = 4$ for the permutation-based method, and to $m = 3$ for the slope-based one in combination with $\delta = 0.001$ and $\gamma = 1$.

### 4.3. Statistical analysis

After testing the hypothesis of gaussianity via the Kolmogorov-Smirnov test, parametric statistical tests were used for each estimator to identify differences in the CE distributions among breathing states. Specifically, the ANOVA analysis of variance followed by a post-hoc paired Student's t-test between the CE values obtained in the spontaneous breathing condition and in each of the controlled one were employed; the Bonferroni correction for multiple comparison ($n = 3$) was applied. For all the statistical tests the significance level was set to 5%.

Moreover, the statistical significance of all the estimated CE value was assessed for each subject, condition and estimator using 100 random shuffle surrogates. Specifically, each measure was deemed as statistically significant if its value computed on the original series was lower than the 5th percentile of the values computed on surrogates.

### 4.4. Results

Fig. 7 reports the CE distributions computed on the heart period time series in the four breathing phases using all the estimation approaches. All methods show a decrease of the entropy rate of RR during C10 and C15 compared to SB, as well as similar values during C20. This trend is always statistically significant for the slowest paced breathing rate, i.e., C10, being instead not noticeable for the C15 breathing rate using kernel and binning estimation approaches (Fig. 7(c,d)). Only the kernel method detects a significant increase of CE during C20, even if
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Performance of the six conditional entropy measures in detecting the magnitude of complexity variations in the dynamics of a mixed non-linear deterministic and linear stochastic process. Errorbar plots depict the percentage variation of the CE computed by increasing the parameter regulating the convex combination of the two dynamics taking consecutive values in the range $\lambda = [0,0.1,0.3,0.5,0.7,0.9,1]$; specific indications about the pair of $\lambda$ values associated with the compared distributions are reported on the x-axis. The estimation parameters are set to $m = 3$ for all methods, and $k = 10$ for $CE_{knn}$, $r = 0.3$ for $CE_{ker}$, $b = 4$ for $CE_{bin}$, and $\delta = 0.001$ and $\gamma = 1$ for $CE_{slope}$. For each estimation approach, asterisk denotes a statistically significant difference between the distributions of the CE obtained for the pair of values selected for $\lambda$ (paired Student's t-test, $p < 0.05$).</figcaption></figure>

<caption>**Table 1**</caption>

Single-run computational times (expressed in milliseconds) for the estimation of CE measure and the other information dynamic measures, averaged on one-hundred realizations of both autoregressive and mixed non-linear deterministic and linear stochastic dynamics processes. The estimation parameters are set to $m = 3$ for all methods, and $k = 10$ for $CE_{knn}$, $r = 0.3$ for $CE_{ker}$, $b = 4$ for $CE_{bin}$, and $\delta = 0.001$ and $\gamma = 1$ for $CE_{slope}$. In the last row, the percentage variations of computational time when increasing the time series length from $N = 300$ to $N = 1000$ are reported.

<table><thead><tr><th rowspan="3">Time series<br>length N</th><th colspan="6">Estimation approaches</th></tr><tr><th colspan="6">Computational times [ms]</th></tr><tr><th>Lin</th><th>Knn</th><th>Ker</th><th>Bin</th><th>Perm</th><th>Slope</th></tr></thead><tbody><tr><td>300</td><td>0.19</td><td>12.37</td><td>9.18</td><td>0.72</td><td>0.53</td><td>0.71</td></tr><tr><td>500</td><td>0.27</td><td>33.85</td><td>21.71</td><td>1</td><td>0.69</td><td>0.73</td></tr><tr><td>1000</td><td>0.31</td><td>121.33</td><td>80.53</td><td>2.43</td><td>1.74</td><td>1.48</td></tr></tbody><tbody><tr><td>From N = 300<br>to N = 1000</td><td colspan="6">Time increment [%]</td></tr><tr><td></td><td>59</td><td>881</td><td>777</td><td>238</td><td>225</td><td>109</td></tr></tbody></table>

Big notation. Specifically, the linear parametric and the discretization-based estimation approaches are of $O(N)$ and $O(N \log(N))$ complexity, respectively. Otherwise, the methods based on computing distances are more demanding, with computational complexities of $O(N^2 \log(N))$ and $O(N^2)$, for the nearest-neighbour and the kernel approaches, respectively. Considering that reliable estimates of entropy measures are obtained for large values of the time series length $N$ and that this is quite larger than the embedding size $m$, i.e., $N \gg m$, the only parameter affecting the computational complexity of the employed algorithms is indeed $N$. On the other hand, the influence of estimator-specific parameters, i.e., $k$ for nearest-neighbour, $r$ for kernel and $b$ for binning, can be neglected in the asymptotic complexity behaviours, even if they in practice cause variations in computational time.

## 4. Application to physiological time series

In this section, conditional entropy estimation approaches are investigated in the context of short-term physiological variability analysis, where heart period and respiratory time series are studied to assess variations in the complexity of cardiorespiratory dynamics during spontaneous and controlled breathing.

### 4.1. Experimental protocol and time series extraction

The investigated physiological signals include electrocardiographic (ECG) and respiratory flow recordings (sampling frequency: 300 Hz) acquired on 19 young healthy subjects (8 males; age range 27–35 years) in resting supine position during different breathing patterns. The analysed experimental conditions were spontaneous breathing (SB) and paced breathing at 10, 15, and 20 breaths/min (C10, C15, and C20, respectively) [18,68,69]. All procedures were approved by the ethical committee of the “L. Sacco” Hospital (Milan, Italy) and of the Department of Technologies for Health, University of Milan (Italy) [68].

The analysis was carried out on the sequence of heart periods (RR intervals) extracted as the time intervals between two consecutive ECG R peaks and on the breathing (RESP) time series obtained by sampling the respiratory flow signal at the onset of every heart period. For each subject and condition, stationary windows of 256 samples were selected synchronously for the two time series. Before estimating CE measures, all time series were preprocessed applying a high-pass AR filter (zero phase and cut-off frequency 0.0156 Hz) and by removing and interpolating samples differing more than three standard deviations from the mean value. Finally, the series were normalized to zero mean and unit variance. Further details on data acquisition and pre-processing can be found in [68].

### 4.2. Data analysis

The CE estimation approaches described in Section 2.3 were applied setting the estimation parameters on the basis of the results of the simulation studies and considering the most widely employed values in literature [40,42,49]. Specifically, the Bayesian Information Criterion (BIC) was employed to define the optimal order of the AR model used to approximate the past dynamics of the investigated process (maximum order fixed to 10), while an embedding dimension $m = 2$ was set to apply nearest-neighbour, kernel and binning approaches. The other estimation parameters set to implement these methods were $k = 10$, $r = 0.3$, and $b = 6$, respectively. The embedding dimension was set to $m = 4$ for the permutation-based method, and to $m = 3$ for the slope-based one in combination with $\delta = 0.001$ and $\gamma = 1$.

### 4.3. Statistical analysis

After testing the hypothesis of gaussianity via the Kolmogorov-Smirnov test, parametric statistical tests were used for each estimator to identify differences in the CE distributions among breathing states. Specifically, the ANOVA analysis of variance followed by a post-hoc paired Student's t-test between the CE values obtained in the spontaneous breathing condition and in each of the controlled one were employed; the Bonferroni correction for multiple comparison ($n = 3$) was applied. For all the statistical tests the significance level was set to 5%.

Moreover, the statistical significance of all the estimated CE value was assessed for each subject, condition and estimator using 100 random shuffle surrogates. Specifically, each measure was deemed as statistically significant if its value computed on the original series was lower than the 5th percentile of the values computed on surrogates.

### 4.4. Results

Fig. 7 reports the CE distributions computed on the heart period time series in the four breathing phases using all the estimation approaches. All methods show a decrease of the entropy rate of RR during C10 and C15 compared to SB, as well as similar values during C20. This trend is always statistically significant for the slowest paced breathing rate, i.e., C10, being instead not noticeable for the C15 breathing rate using kernel and binning estimation approaches (Fig. 7(c,d)). Only the kernel method detects a significant increase of CE during C20, even if
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Performance of the six conditional entropy measures as a function of the estimation parameters (left column: embedding dimension *m*; right column: estimator-specific parameter) for time series of length *N* representing a mixed non-linear deterministic and linear stochastic process with parameter *λ* = 0.5. Errorbar plots (mean ± std) depict the distributions of CE computed over 100 realizations using the (a) linear (CE<sub>lin</sub>), (b) nearest neighbour (CE<sub>knn</sub>), (c) kernel (CE<sub>ker</sub>), (d) binning (CE<sub>bin</sub>), (e) permutation (CE<sub>perm</sub>), and (f) slope (CE<sub>slope</sub>) approaches.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Boxplots distributions and individual values of the CE computed for the RR time series during spontaneous (SB) and controlled breathing at 10, 15, and 20 breaths/min (C10, C15, and C20) with (a) linear, (b) nearest-neighbours, (c) kernel, (d) binning, (e) permutation, and (f) slope approaches. In each panel, empty circles correspond to values deemed as statistically non-significant using surrogate data analysis; black open circles refer to the sample mean of each distribution. The symbol # indicates statistically significant differences across the four protocol phases (ANOVA test, *p* < 0.05), while the symbol * denotes statistically significant differences between SB and C10/C15/C20 phases (paired Student's t-test with Bonferroni correction for multiple comparison, *p* < 0.05/3).</figcaption></figure>

the measure is not significant according to surrogate analysis for about 26% of the cohort protocol (Fig. 7(c)).

Fig. 8 reports the distributions of the CE measure computed on the respiratory time series during the four breathing phases using all estimation approaches. Except for the slope estimation method (Fig. 8(f)), trends similar to those reported for RR time series are also reported for RESP, with a statistically significant decrease in complexity shown during C10 with all estimators and during C15 with the linear and the permutation ones. Like for the heart period time series, the permutation

method shows the largest relative variations in magnitude of the CE measure during C10 compared with SB.

## 5. Discussion

In this work, a comparative study of several methods for estimating complexity of time series was carried out on both simulated and cardiorespiratory dynamics. Measures of entropy rate computing the CE of time series have been proven to be an useful and reliable tool for the analysis of complex system dynamics in real data in several fields [48],
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Performance of the six conditional entropy measures as a function of the estimation parameters (left column: embedding dimension *m*; right column: estimator-specific parameter) for time series of length *N* representing a mixed non-linear deterministic and linear stochastic process with parameter *λ* = 0.5. Errorbar plots (mean ± std) depict the distributions of CE computed over 100 realizations using the (a) linear (CE<sub>lin</sub>), (b) nearest neighbour (CE<sub>knn</sub>), (c) kernel (CE<sub>ker</sub>), (d) binning (CE<sub>bin</sub>), (e) permutation (CE<sub>perm</sub>), and (f) slope (CE<sub>slope</sub>) approaches.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Boxplots distributions and individual values of the CE computed for the RR time series during spontaneous (SB) and controlled breathing at 10, 15, and 20 breaths/min (C10, C15, and C20) with (a) linear, (b) nearest-neighbours, (c) kernel, (d) binning, (e) permutation, and (f) slope approaches. In each panel, empty circles correspond to values deemed as statistically non-significant using surrogate data analysis; black open circles refer to the sample mean of each distribution. The symbol # indicates statistically significant differences across the four protocol phases (ANOVA test, *p* < 0.05), while the symbol * denotes statistically significant differences between SB and C10/C15/C20 phases (paired Student's *t*-test with Bonferroni correction for multiple comparison, *p* < 0.05/3).</figcaption></figure>

the measure is not significant according to surrogate analysis for about 26% of the cohort protocol (Fig. 7(c)).

**Fig. 8** reports the distributions of the CE measure computed on the respiratory time series during the four breathing phases using all estimation approaches. Except for the slope estimation method (Fig. 8(f)), trends similar to those reported for RR time series are also reported for RESP, with a statistically significant decrease in complexity shown during C10 with all estimators and during C15 with the linear and the permutation ones. Like for the heart period time series, the permutation

method shows the largest relative variations in magnitude of the CE measure during C10 compared with SB.

## 5. Discussion

In this work, a comparative study of several methods for estimating complexity of time series was carried out on both simulated and cardiorespiratory dynamics. Measures of entropy rate computing the CE of time series have been proven to be an useful and reliable tool for the analysis of complex system dynamics in real data in several fields [48],
<figure><img src="image_11.png" /><figcaption>Fig. 8. Boxplots distributions and individual values of the CE computed for the RESP time series during spontaneous (SB) and controlled breathing at 10, 15, and 20 breaths/min (C10, C15, and C20) with (a) linear, (b) nearest-neighbours, (c) kernel, (d) binning, (e) permutation, and (f) slope methods. In each panel, empty circles correspond to values deemed as statistically non-significant using surrogate data analysis; black open circles refer to the sample mean of each distribution. The symbol # indicates statistically significant differences across the four protocol phases (ANOVA test, *p* < 0.05), while the symbol * denotes statistically significant differences between SB and C10/C15/C20 phases (paired Student's *t*-test with Bonferroni correction for multiple comparison, *p* < 0.05/3).</figcaption></figure>

e.g., finance [70,71], climate [72,73], or biomedicine [12,21,35]. Here, we have applied methods for estimating CE based on continuous and discrete representations of the system dynamics, as well as parametric and model-free estimation of probability distributions.

In both simulated settings, involving linear stochastic and chaotic deterministic dynamics, all estimators show a general CE decrease as the complexity of the system, understood in terms of the unpredictability of the dynamics, is reduced. Furthermore, the application to physiological data shows how more regular dynamics associated with the strengthening of specific regulatory mechanisms lead to a decrease in the complexity of the system. Specifically, the reduced complexity of heart period dynamics during controlled breathing at a rate slower than the spontaneous breathing rate can be related to the respiratory sinus arrhythmia (RSA) mechanism, responsible for the modulation of the cardiac dynamics by breathing activity [74], as this mechanism is enhanced during forced ventilation at lower respiratory rates, resulting in a more regular heart rhythm [75]. On the other hand, the cardio-ventilatory coupling [76] and the respiratory stroke volume synchronization [77] mechanisms, which induce changes in respiratory dynamics similar to those of the heart, can be responsible for the comparable trends obtained when observing the complexity of the breathing patterns. These results are well established in the literature and have been reported in terms of several complexity, e.g., CE [18, 49], information storage [63], non-linear prediction [78,79], and interactions, e.g., cross-entropy [63], transfer entropy [49,80], mutual information rate [49], indices for the investigation of physiological network dynamics.

Nevertheless, our findings highlight that selecting the most appropriate method and estimation parameters to provide an accurate measure according to system characteristics is a non trivial task. The linear parametric estimator is widely used to study the complexity of dynamical systems, showing its efficiency in different physiological frameworks [16,47,62,81]. Several works emphasize that this approach is mostly suited for Gaussian processes, being not able to take into account non-linear dynamics [48,56]. This was also confirmed by our results of the second simulation, in which the linear estimator was unable to detect the complexity reduction associated only with non-linear dynamics.

In the literature, various model-free approaches have been developed to move from a linear and simplified representation of the process

to a quantification of probability distributions directly from the data. Among these approaches, the one based on the nearest neighbour metric exhibits high performance for non-Gaussian processes [12,21, 48,82,83]. These findings are confirmed by our results, which show that the CE estimates obtained using the *knn* approach on simulated AR linear dynamics follow the theoretical value of the measure and reveal statistically significant variations in the complexity of the system that are also detected by the linear estimator, allowing as well to strongly discriminate complexity changes also on mixed linear and chaotic dynamics. Despite the bias related to the dimensionality of the space in the search for neighbour samples, this metric allows to compute consistent CE values at varying of the embedding size *m* and the number of neighbours *k*. Moreover, the bias related to the parameter *N* is negligible compared to the model-free estimation approaches based on discretization and symbolization of the process dynamics. The high dependence of the latter approaches on the time series length also results in a strong dependence on the estimation parameters [35].

Although being very popular for calculating approximate entropy and sample entropy to assess changes in physiological complexity under different physiological and pathological conditions [13,19,22,27,29, 84-87], the kernel method shows a very high estimation variability. For the two simulated dynamics, our results highlight that when the embedding dimension *m* increases and the threshold parameter *r* decreases, and thus the searching space becomes increasingly sparse, the kernel estimator may become unfeasible for estimating the CE, especially for shorter time series. This is likely the reason behind the large variability of the estimates, which causes the kernel estimator to fail in detecting small complexity variations.

Regarding model-free estimation methods based on discretization, the most commonly used is the binning approach [20,30,51,88], mainly thanks to its simplicity. Nonetheless, our results evidence that this approach does not allow the proper identification of small complexity changes. The reason behind this inaccuracy may be in the excessive coarse graining performed by binning, which cannot be avoided when dealing with short time series due to the curse of dimensionality [35]. As observed in previous works [49,61], the permutation-based approach uses a smaller number of patterns for the same embedding size compared to the binning one, making it more efficient in assessing the dynamics of a process. Moreover, this method correctly discriminates
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Boxplots distributions and individual values of the CE computed for the RESP time series during spontaneous (SB) and controlled breathing at 10, 15, and 20 breaths/min (C10, C15, and C20) with (a) linear, (b) nearest-neighbours, (c) kernel, (d) binning, (e) permutation, and (f) slope methods. In each panel, empty circles correspond to values deemed as statistically non-significant using surrogate data analysis; black open circles refer to the sample mean of each distribution. The symbol # indicates statistically significant differences across the four protocol phases (ANOVA test, *p* < 0.05), while the symbol * denotes statistically significant differences between SB and C10/C15/C20 phases (paired Student's *t*-test with Bonferroni correction for multiple comparison, *p* < 0.05/3).</figcaption></figure>

e.g., finance [70,71], climate [72,73], or biomedicine [12,21,35]. Here, we have applied methods for estimating CE based on continuous and discrete representations of the system dynamics, as well as parametric and model-free estimation of probability distributions.

In both simulated settings, involving linear stochastic and chaotic deterministic dynamics, all estimators show a general CE decrease as the complexity of the system, understood in terms of the unpredictability of the dynamics, is reduced. Furthermore, the application to physiological data shows how more regular dynamics associated with the strengthening of specific regulatory mechanisms lead to a decrease in the complexity of the system. Specifically, the reduced complexity of heart period dynamics during controlled breathing at a rate slower than the spontaneous breathing rate can be related to the respiratory sinus arrhythmia (RSA) mechanism, responsible for the modulation of the cardiac dynamics by breathing activity [74], as this mechanism is enhanced during forced ventilation at lower respiratory rates, resulting in a more regular heart rhythm [75]. On the other hand, the cardio-ventilatory coupling [76] and the respiratory stroke volume synchronization [77] mechanisms, which induce changes in respiratory dynamics similar to those of the heart, can be responsible for the comparable trends obtained when observing the complexity of the breathing patterns. These results are well established in the literature and have been reported in terms of several complexity, e.g., CE [18, 49], information storage [63], non-linear prediction [78,79], and interactions, e.g., cross-entropy [63], transfer entropy [49,80], mutual information rate [49], indices for the investigation of physiological network dynamics.

Nevertheless, our findings highlight that selecting the most appropriate method and estimation parameters to provide an accurate measure according to system characteristics is a non trivial task. The linear parametric estimator is widely used to study the complexity of dynamical systems, showing its efficiency in different physiological frameworks [16,47,62,81]. Several works emphasize that this approach is mostly suited for Gaussian processes, being not able to take into account non-linear dynamics [48,56]. This was also confirmed by our results of the second simulation, in which the linear estimator was unable to detect the complexity reduction associated only with non-linear dynamics.

In the literature, various model-free approaches have been developed to move from a linear and simplified representation of the process

to a quantification of probability distributions directly from the data. Among these approaches, the one based on the nearest neighbour metric exhibits high performance for non-Gaussian processes [12,21, 48,82,83]. These findings are confirmed by our results, which show that the CE estimates obtained using the *knn* approach on simulated AR linear dynamics follow the theoretical value of the measure and reveal statistically significant variations in the complexity of the system that are also detected by the linear estimator, allowing as well to strongly discriminate complexity changes also on mixed linear and chaotic dynamics. Despite the bias related to the dimensionality of the space in the search for neighbour samples, this metric allows to compute consistent CE values at varying of the embedding size *m* and the number of neighbours *k*. Moreover, the bias related to the parameter *N* is negligible compared to the model-free estimation approaches based on discretization and symbolization of the process dynamics. The high dependence of the latter approaches on the time series length also results in a strong dependence on the estimation parameters [35].

Although being very popular for calculating approximate entropy and sample entropy to assess changes in physiological complexity under different physiological and pathological conditions [13,19,22,27,29, 84-87], the kernel method shows a very high estimation variability. For the two simulated dynamics, our results highlight that when the embedding dimension *m* increases and the threshold parameter *r* decreases, and thus the searching space becomes increasingly sparse, the kernel estimator may become unfeasible for estimating the CE, especially for shorter time series. This is likely the reason behind the large variability of the estimates, which causes the kernel estimator to fail in detecting small complexity variations.

Regarding model-free estimation methods based on discretization, the most commonly used is the binning approach [20,30,51,88], mainly thanks to its simplicity. Nonetheless, our results evidence that this approach does not allow the proper identification of small complexity changes. The reason behind this inaccuracy may be in the excessive coarse graining performed by binning, which cannot be avoided when dealing with short time series due to the curse of dimensionality [35]. As observed in previous works [49,61], the permutation-based approach uses a smaller number of patterns for the same embedding size compared to the binning one, making it more efficient in assessing the dynamics of a process. Moreover, this method correctly discriminates
relevant variations in the complex behaviour of stochastic system dynamics. However, given its limitation in neglecting the amplitude of the dynamics, it has a greater difficulty in detecting small variations in non-linear dynamics [89]; this aspect was verified in our simulations. The recently developed slope-based method overcomes this limitation by taking into account the amplitude deviation between successive time series samples [42]. Our results evidence that this approach, working on the shorter slope vectors, in turn allows to detect smaller complexity variations when compared to the other methods. Nevertheless, it is worth emphasising that this approach investigates predictability by looking not at the time series samples, but at the differences between them, and how this may result in smoothing variations.

The remarks concerning the simulations are consistent with the results obtained analysing the complexity of the cardiorespiratory dynamics. Strong mechanisms such as those coordinating the cardiac and respiratory dynamics during C10 are mostly discriminated by all estimators, while more weakly synchronized dynamics in C15 are only detected by the linear estimator and the permutation-based method for the respiratory dynamics, and also by the nearest-neighbour and the slope estimators for the cardiac ones. As reported in literature, the predominance of linear dynamics exhibited by the mechanism of respiratory sinus arrhythmia during slow breathing conditions [90] is also evidenced in our results, given the superiority of the linear estimator and the permutation-based approach in discriminating physiological complexity variations. The inability of the slope CE estimator in detecting changes in respiratory dynamics across the different breathing conditions can be related to the above reported methodological difference of this approach compared to the others. It is worth underlying that our analyses have only taken into account physiological states related to parasympathetic activity, highlighting how the measure of conditional entropy is efficient in detecting strong variations in these dynamics. More evidences are reported in the literature on the ability to distinguish, through complexity measures, dynamics related to the sympathetic activity [16,91].

The analysis of the computational complexity for the different CE estimation methods shows that the linear and discretization-based approaches are less demanding compared to the nearest-neighbour and kernel estimation methods. In addition, the overall computational times follow the complexity trend, although these are only indicative of our specific algorithms implementation and strongly depend on the hardware and software computer equipment. Nevertheless, as previously discussed, entropy estimators based on computing distances often provide a more reliable complexity assessment even in presence of non-linear dynamics. For this reason, several works have aimed at implementing faster and more computationally efficient algorithms for both the nearest-neighbour [92–94] and kernel [95–98] approaches. Another important remark about computational times reported in this study is on the use of MATLAB-based algorithms. Being an interpreted environment, MATLAB has higher execution times if compared to compiled languages, e.g., C and C++. For this reason, a possible alternative to increase MATLAB codes performances is to switch to C/C++ algorithms and functions, i.e., MEX files.

The simplicity of the linear and permutation-based approaches and their speed in the context of short-term analysis, i.e., using time series of 300 samples (~5 min), make such measures attractive for their implementation within the firmware of wearable biomedical devices. For example, CE can provide a feasible assessment of cardiac dynamics variability starting from electrocardiographic or photoplethysmographic signals acquired during everyday activities and outside the hospital environment, e.g., home health care services aimed at long-term monitoring or follow-up of frail patients [99–101]. In this regard, both power cost and stationarity of physiological signals become important issues which have pushed the research towards the acquisition and the analysis of even shorter time-series. Several recent works have thus focused on the assessment of physiological parameters using ultra-short term analysis (i.e., <300 samples), highlighting the feasibility of

obtaining reliable $CE_{lin}$ estimates on 120-samples heart rate variability series (~2 min) [81,102]. This would not be achievable through the discretization approach as using shorter series and the same estimation parameters would lead to unreliable estimates [35]. Despite its lower computational costs, the linear approach can however compute the CE value once the signal acquisition has been completed. On the other hand, the dynamics discretization of the permutation-based approach can be carried out in parallel with the signal acquisition, leaving the evaluation of the probability distribution at the end, thus improving the real-time performance and reducing the overall computational time.

## 6. Conclusion

The analyses carried out in this work have demonstrated that the linear parametric approach remains the most valid and quick estimation method for assessing complexity, being able to discriminate complexity variations in real systems such as those characterising cardiorespiratory dynamics in healthy subjects. Among the model-free approaches, the nearest neighbour metric allows to adequately study linear and non-linear dynamics. Together with the stability of the estimated measures when varying the estimation parameters, these features make such method one of the most widely adopted despite its higher computational cost compared to parametric and other model-free approaches. Among the latter, the permutation-based approach is simple and, setting suitable embedding parameters, can be used effectively to assess the complexity of system dynamics.

Although a broad overview and description of several valid and useful approaches to complexity assessment have been provided, this work does not take into account some recently introduced and investigated estimation methods that try to overcome some of the previously mentioned drawbacks. Further analysis should therefore include the conditional entropy implementation of other estimation approaches that limit the dependence on the estimation parameters, e.g., Rank [103] or Bubble Entropy [104], or on the intrinsic signal features, e.g., Range [105] or Diversity Entropy [106], or improve the consistency and validity of the measure, e.g., Fuzzy [107], Distribution [108], Dispersion [109] and Phase Entropy [110]. Future studies are also needed to explore the complementarity of the investigated estimation methods in different physiological and pathological scenarios characterized by different dynamical behaviours, such as non-linearity [28,111–113] or non-stationarity [83,114]. Moreover, in order to appropriately assess the complexity of physiological dynamics, it is envisaged to investigate multivariate approaches for evaluating the measure of conditional entropy, which allow to take into account not only the internal regulatory mechanisms of individual physiological systems, but also their interconnection [12,115].

## Funding

The present work was supported by the Italian MUR PRIN 2022 "High-Order Dynamical Networks in Computational Neuroscience and Physiology: an Information-Theoretic Approach" (HONEST), project code 2022YMHNPY, CUP: B53D23003020006, and by the Italian Complementary National Plan PNC-I.1 "Research initiatives for innovative technologies and pathways in the health and welfare sector", D.D. 931 of 06/06/2022, "DARE - DigitAl lifelong pRevEntion" initiative, code PNC0000002, CUP B53C22006460001. R.P. was partially supported by European Social Fund (ESF)—Complementary Operational Programme (POC) 2014/2020 of the Sicily Region.

## CRediT authorship contribution statement

Chiara Barà: Formal analysis, Investigation, Software, Visualization, Writing – original draft. Riccardo Pernice: Software, Supervision, Visualization, Writing – review & editing. Cristina Angela Catania: Formal analysis, Validation, Writing – review & editing. Mirvana Hilal:
relevant variations in the complex behaviour of stochastic system dynamics. However, given its limitation in neglecting the amplitude of the dynamics, it has a greater difficulty in detecting small variations in non-linear dynamics [89]; this aspect was verified in our simulations. The recently developed slope-based method overcomes this limitation by taking into account the amplitude deviation between successive time series samples [42]. Our results evidence that this approach, working on the shorter slope vectors, in turn allows to detect smaller complexity variations when compared to the other methods. Nevertheless, it is worth emphasising that this approach investigates predictability by looking not at the time series samples, but at the differences between them, and how this may result in smoothing variations.

The remarks concerning the simulations are consistent with the results obtained analysing the complexity of the cardiorespiratory dynamics. Strong mechanisms such as those coordinating the cardiac and respiratory dynamics during C10 are mostly discriminated by all estimators, while more weakly synchronized dynamics in C15 are only detected by the linear estimator and the permutation-based method for the respiratory dynamics, and also by the nearest-neighbour and the slope estimators for the cardiac ones. As reported in literature, the predominance of linear dynamics exhibited by the mechanism of respiratory sinus arrhythmia during slow breathing conditions [90] is also evidenced in our results, given the superiority of the linear estimator and the permutation-based approach in discriminating physiological complexity variations. The inability of the slope CE estimator in detecting changes in respiratory dynamics across the different breathing conditions can be related to the above reported methodological difference of this approach compared to the others. It is worth underlying that our analyses have only taken into account physiological states related to parasympathetic activity, highlighting how the measure of conditional entropy is efficient in detecting strong variations in these dynamics. More evidences are reported in the literature on the ability to distinguish, through complexity measures, dynamics related to the sympathetic activity [16,91].

The analysis of the computational complexity for the different CE estimation methods shows that the linear and discretization-based approaches are less demanding compared to the nearest-neighbour and kernel estimation methods. In addition, the overall computational times follow the complexity trend, although these are only indicative of our specific algorithms implementation and strongly depend on the hardware and software computer equipment. Nevertheless, as previously discussed, entropy estimators based on computing distances often provide a more reliable complexity assessment even in presence of non-linear dynamics. For this reason, several works have aimed at implementing faster and more computationally efficient algorithms for both the nearest-neighbour [92–94] and kernel [95–98] approaches. Another important remark about computational times reported in this study is on the use of MATLAB-based algorithms. Being an interpreted environment, MATLAB has higher execution times if compared to compiled languages, e.g., C and C++. For this reason, a possible alternative to increase MATLAB codes performances is to switch to C/C++ algorithms and functions, i.e., MEX files.

The simplicity of the linear and permutation-based approaches and their speed in the context of short-term analysis, i.e., using time series of 300 samples (~5 min), make such measures attractive for their implementation within the firmware of wearable biomedical devices. For example, CE can provide a feasible assessment of cardiac dynamics variability starting from electrocardiographic or photoplethysmographic signals acquired during everyday activities and outside the hospital environment, e.g., home health care services aimed at long-term monitoring or follow-up of frail patients [99–101]. In this regard, both power cost and stationarity of physiological signals become important issues which have pushed the research towards the acquisition and the analysis of even shorter time-series. Several recent works have thus focused on the assessment of physiological parameters using ultra-short term analysis (i.e., <300 samples), highlighting the feasibility of

obtaining reliable $CE_{lin}$ estimates on 120-samples heart rate variability series (~2 min) [81,102]. This would not be achievable through the discretization approach as using shorter series and the same estimation parameters would lead to unreliable estimates [35]. Despite its lower computational costs, the linear approach can however compute the CE value once the signal acquisition has been completed. On the other hand, the dynamics discretization of the permutation-based approach can be carried out in parallel with the signal acquisition, leaving the evaluation of the probability distribution at the end, thus improving the real-time performance and reducing the overall computational time.

## 6. Conclusion

The analyses carried out in this work have demonstrated that the linear parametric approach remains the most valid and quick estimation method for assessing complexity, being able to discriminate complexity variations in real systems such as those characterising cardiorespiratory dynamics in healthy subjects. Among the model-free approaches, the nearest neighbour metric allows to adequately study linear and non-linear dynamics. Together with the stability of the estimated measures when varying the estimation parameters, these features make such method one of the most widely adopted despite its higher computational cost compared to parametric and other model-free approaches. Among the latter, the permutation-based approach is simple and, setting suitable embedding parameters, can be used effectively to assess the complexity of system dynamics.

Although a broad overview and description of several valid and useful approaches to complexity assessment have been provided, this work does not take into account some recently introduced and investigated estimation methods that try to overcome some of the previously mentioned drawbacks. Further analysis should therefore include the conditional entropy implementation of other estimation approaches that limit the dependence on the estimation parameters, e.g., Rank [103] or Bubble Entropy [104], or on the intrinsic signal features, e.g., Range [105] or Diversity Entropy [106], or improve the consistency and validity of the measure, e.g., Fuzzy [107], Distribution [108], Dispersion [109] and Phase Entropy [110]. Future studies are also needed to explore the complementarity of the investigated estimation methods in different physiological and pathological scenarios characterized by different dynamical behaviours, such as non-linearity [28,111–113] or non-stationarity [83,114]. Moreover, in order to appropriately assess the complexity of physiological dynamics, it is envisaged to investigate multivariate approaches for evaluating the measure of conditional entropy, which allow to take into account not only the internal regulatory mechanisms of individual physiological systems, but also their interconnection [12,115].

## Funding

The present work was supported by the Italian MUR PRIN 2022 "High-Order Dynamical Networks in Computational Neuroscience and Physiology: an Information-Theoretic Approach" (HONEST), project code 2022YMHNPY, CUP: B53D23003020006, and by the Italian Complementary National Plan PNC-I.1 "Research initiatives for innovative technologies and pathways in the health and welfare sector", D.D. 931 of 06/06/2022, "DARE - DigitAl lifelong pRevEntion" initiative, code PNC0000002, CUP B53C22006460001. R.P. was partially supported by European Social Fund (ESF)—Complementary Operational Programme (POC) 2014/2020 of the Sicily Region.

## CRediT authorship contribution statement

Chiara Barà: Formal analysis, Investigation, Software, Visualization, Writing – original draft. Riccardo Pernice: Software, Supervision, Visualization, Writing – review & editing. Cristina Angela Catania: Formal analysis, Validation, Writing – review & editing. Mirvana Hilal:
Validation, Writing – review & editing. **Alberto Porta:** Data curation, Writing – review & editing. **Anne Humeau-Heurtier:** Methodology, Supervision, Writing – review & editing. **Luca Faes:** Conceptualization, Methodology, Resources, Software, Supervision, Validation, Writing – review & editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References

[1] Ivanov PC. The new field of network physiology: building the human physiome. Front Netw Physiol 2021;1:711778.

[2] Boccaletti S, Latora V, Moreno Y, Chavez M, Hwang D-U. Complex networks: Structure and dynamics. Phys Rep 2006;424(4-5):175-308.

[3] Siegenfeld AF, Bar-Yam Y. An introduction to complex systems science and its applications. Complexity 2020;20(2):1-16.

[4] Ladyman J, Lambert J, Wiesner K. What is a complex system? Eur J Philos Sci 2013;3:33-67.

[5] Foote R. Mathematics and complex systems. science 2007;318(5849):410-2.

[6] Goldberger AL. Giles f. filley lecture. complex systems. Proc Am Thorac Soc 2006;3(6):467-71.

[7] Mihailović DT, Mimić G, Nikolić-Djorić E, Arsenić I. Novel measures based on the kolmogorov complexity for use in complex system behavior studies and time series analysis. Open Phys 2015;13(1).

[8] Arshinov V, Fuchs C. Causality, emergence, self-organisation. NIA-Priroda Moscow; 2003.

[9] Cannon J. The fractal geometry of nature, by benoit b. mandelbrot. Amer Math Monthly 1984;91(9):594-8.

[10] Kaplan D, Glass L. Understanding nonlinear dynamics. Springer Science & Business Media; 1997.

[11] Glass L. Synchronization and rhythmic processes in physiology. Nature 2001;410(6825):277-84.

[12] Valente M, Javorka M, Porta A, Bari V, Krohova J, Czippelova B, Turianikova Z, Nollo G, Faes L. Univariate and multivariate conditional entropy measures for the characterization of short-term cardiovascular complexity under physiological stress. Physiol Meas 2018;39(1):014002.

[13] Goldberger AL, Peng C-K, Lipsitz LA. What is physiologic complexity and how does it change with aging and disease? Neurobiol Aging 2002;23(1):23-6.

[14] Shaffer F, McCraty R, Zerr CL. A healthy heart is not a metronome: an integrative review of the heart's anatomy and heart rate variability. Front Psychol 2014;5:1040.

[15] Lehrer P, Eddie D. Dynamic processes in regulation and some implications for biofeedback and biobehavioral interventions. Appl Psychophysiol biofeedback 2013;38:143-55.

[16] Javorka M, Krohova J, Czippelova B, Turianikova Z, Lazarova Z, Wiszt R, Faes L. Towards understanding the complexity of cardiovascular oscillations: Insights from information theory. Comput biol Med 2018;98:48-57.

[17] Javorka M, Turianikova Z, Tonhajzerova I, Javorka K, Baumert M. The effect of orthostasis on recurrence quantification analysis of heart rate and blood pressure dynamics. Physiol Meas 2008;30(1):29.

[18] Porta A, Bari V, Ranuzzi G, De Maria B, Baselli G. Assessing multiscale complexity of short heart rate variability series through a model-based linear approach. Chaos 2017;27(9).

[19] Heffernan KS, Fahs CA, Shinsako KK, Jae SY, Fernhall B. Heart rate recovery and heart rate complexity following resistance exercise training and detraining in young men. Am J Physiol-Heart Circ Physiol 2007;293(5):H3180-6.

[20] Takahashi AC, Porta A, Melo RC, Quitério RJ, da Silva E, Borghi-Silva A, Tobaldini E, Montano N, Catai AM. Aging reduces complexity of heart rate variability assessed by conditional entropy and symbolic analysis. Intern Emerg Med 2012;7:229-35.

[21] Porta A, Faes L, Bari V, Marchi A, Bassani T, Nollo G, Perseguini NM, Milan J, Minatel V, Borghi-Silva A, et al. Effect of age on complexity and causality of the cardiovascular control: comparison between model-based and model-free approaches. PLoS One 2014;9(2):e89463.

[22] Jia Y, Gu H, Luo Q. Sample entropy reveals an age-related reduction in the complexity of dynamic brain. Sci Rep 2017;7(1):7990.

[23] Pincus SM. Greater signal regularity may indicate increased system isolation. Math Biosci 1994;122(2):161-81.

[24] Romero-Ortuño R, Martínez-Velilla N, Sutton R, Ungar A, Fedorowski A, Galvin R, Theou O, Davies A, Reilly RB, Claassen J, et al. Network physiology in aging and frailty: the grand challenge of physiological reserve in older adults. 2021.

[25] Méndez MA, Zuluaga P, Hornero R, Gómez C, Escudero J, Rodríguez-Palancas A, Ortiz T, Fernández A. Complexity analysis of spontaneous brain activity: effects of depression and antidepressant treatment. J Psychopharmacol 2012;26(5):636-43.

[26] Fernández A, Quintero J, Hornero R, Zuluaga P, Navas M, Gómez C, Escudero J, García-Campos N, Biederman J, Ortiz T. Complexity analysis of spontaneous brain activity in attention-deficit/hyperactivity disorder: diagnostic implications. Biol Psychiatry 2009;65(7):571-7.

[27] Hornero R, Abásolo D, Jimeno N, Sánchez CI, Poza J, Aboy M. Variability, regularity, and complexity of time series generated by schizophrenic patients and control subjects. IEEE Trans Biomed Eng 2006;53(2):210-8.

[28] Chen C, Jin Y, Lo IL, Zhao H, Sun B, Zhao Q, Zheng J, Zhang XD. Complexity change in cardiovascular disease. Int J Biol Sci 2017;13(10):1320.

[29] Trunkvalterova Z, Javorka M, Tonhajzerova I, Javorkova J, Lazarova Z, Javorka K, Baumert M. Reduced short-term complexity of heart rate and blood pressure dynamics in patients with diabetes mellitus type 1: multiscale entropy analysis. Physiol Meas 2008;29(7):817.

[30] Tobaldini E, Nobili L, Strada S, Casali KR, Braghiroli A, Montano N. Heart rate variability in normal and pathological sleep. Front Physiol 2013;4:294.

[31] Theiler J. Estimating fractal dimension. J Opt Soc Amer A 1990;7(6):1055-73.

[32] Wolf A, Swift JB, Swinney HL, Vastano JA. Determining lyapunov exponents from a time series. Physica D 1985;16(3):285-317.

[33] Lempel A, Ziv J. On the complexity of finite sequences. IEEE Trans Inform Theory 1976;22(1):75-81.

[34] Baranger M, Institute NECS. Chaos, complexity, and entropy: A physics talk for non-pharmacists. New England Complex Systems Institute; 2001, URL https://books.google.sk/books?id=duzMgEACAAJ.

[35] Faes L, Porta A. Conditional entropy-based evaluation of information dynamics in physiological systems. Direct Inf Meas Neurosci 2014;61-86.

[36] Lizier JT. The local information dynamics of distributed computation in complex systems. Springer Science & Business Media; 2012.

[37] Faes L, Erla S, Nollo G. Measuring connectivity in linear multivariate processes: definitions, interpretation, and practical analysis. Comput Math Methods Med 2012;2012.

[38] Kozachenko LF, Leonenko NN. Sample estimate of the entropy of a random vector. Probl Peredachi Inform 1987;23(2):9-16.

[39] Richman JS, Moorman JR. Physiological time-series analysis using approximate entropy and sample entropy. Am J Physiol-Heart Circ Physiol 2000.

[40] Azami H, Faes L, Escudero J, Humeau-Heurtier A, Silva LE. Entropy analysis of univariate biomedical signals: Review and comparison of methods. Front Entropy across Discipl: Panorama Entropy: Theory Comput Appl 2023;233-86.

[41] Bandt C, Pompe B. Permutation entropy: a natural complexity measure for time series. Phys Rev Lett 2002;88(17):174102.

[42] Cuesta-Frau D. Slope entropy: A new time series complexity estimator based on both symbolic patterns and amplitude information. Entropy 2019;21(12):1167.

[43] Shannon CE. A mathematical theory of communication. Bell Syst Tech J 1948;27(3):379-423.

[44] Cover TM. Elements of information theory. John Wiley & Sons; 1999.

[45] Kolmogorov AN. A new metric invariant of transitive dynamical systems and automorphisms of lebesgue spaces. Trudy Mat Inst Imeni VA Steklova 1985;169:94-8.

[46] Sinai YG. On the notion of entropy of a dynamical system. In: Doklady of Russian academy of sciences. Vol. 124, 1959, p. 768-71.

[47] Porta A, De Maria B, Bari V, Marchi A, Faes L. Are nonlinear model-free conditional entropy approaches for the assessment of cardiac control complexity superior to the linear model-based one? IEEE Trans Biomed Eng 2016;64(6):1287-96.

[48] Xiong W, Faes L, Ivanov PC. Entropy measures, entropy estimators, and their performance in quantifying complex dynamics: Effects of artifacts, nonstationarity, and long-range correlations. Phys Rev E 2017;95(6):062114.

[49] Barà C, Sparacino L, Pernice R, Antonacci Y, Porta A, Kugiumtzis D, Faes L. Comparison of discretization strategies for the model-free information-theoretic assessment of short-term physiological interactions. Chaos 2023;33(3):033127.

[50] Runge J, Heitzig J, Petoukhov V, Kurths J. Escaping the curse of dimensionality in estimating multivariate transfer entropy. Phys Rev Lett 2012;108(25):258701.

[51] Porta A, Baselli G, Liberati D, Montano N, Cogliati C, Gneчи-Ruscone T, Malliani A, Cerutti S. Measuring regularity by means of a corrected conditional entropy in sympathetic outflow. Biol Cybern 1998;78:71-8.

[52] Barrett AB, Barnett L, Seth AK. Multivariate granger causality and generalized variance. Phys Rev E 2010;81(4):041907.

[53] Mardia KV, Marshall RJ. Maximum likelihood estimation of models for residual covariance in spatial regression. Biometrika 1984;71(1):135-46.

[54] Akaike H. A new look at the statistical model identification. IEEE Trans Autom Control 1974;19(6):716-23.

[55] Schwarz G. Estimating the dimension of a model. Ann Stat 1978;461-4.

[56] Barnett L, Barrett AB, Seth AK. Granger causality and transfer entropy are equivalent for gaussian variables. Phys Rev Lett 2009;103(23):238701.

[57] Kraskov A, Stögbauer H, Grassberger P. Estimating mutual information. Phys Rev E 2004;69(6):066138.
Validation, Writing – review & editing. **Alberto Porta:** Data curation, Writing – review & editing. **Anne Humeau-Heurtier:** Methodology, Supervision, Writing – review & editing. **Luca Faes:** Conceptualization, Methodology, Resources, Software, Supervision, Validation, Writing – review & editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References

[1] Ivanov PC. The new field of network physiology: building the human physiome. Front Netw Physiol 2021;1:711778.

[2] Boccaletti S, Latora V, Moreno Y, Chavez M, Hwang D-U. Complex networks: Structure and dynamics. Phys Rep 2006;424(4-5):175-308.

[3] Siegenfeld AF, Bar-Yam Y. An introduction to complex systems science and its applications. Complexity 2020;20(2):1-16.

[4] Ladyman J, Lambert J, Wiesner K. What is a complex system? Eur J Philos Sci 2013;3:33-67.

[5] Foote R. Mathematics and complex systems. science 2007;318(5849):410-2.

[6] Goldberger AL. Giles f. filley lecture. complex systems. Proc Am Thorac Soc 2006;3(6):467-71.

[7] Mihailović DT, Mimić G, Nikolić-Djorić E, Arsenić I. Novel measures based on the kolmogorov complexity for use in complex system behavior studies and time series analysis. Open Phys 2015;13(1).

[8] Arshinov V, Fuchs C. Causality, emergence, self-organisation. NIA-Priroda Moscow; 2003.

[9] Cannon J. The fractal geometry of nature, by benoit b. mandelbrot. Amer Math Monthly 1984;91(9):594-8.

[10] Kaplan D, Glass L. Understanding nonlinear dynamics. Springer Science & Business Media; 1997.

[11] Glass L. Synchronization and rhythmic processes in physiology. Nature 2001;410(6825):277-84.

[12] Valente M, Javorka M, Porta A, Bari V, Krohova J, Czippelova B, Turianikova Z, Nollo G, Faes L. Univariate and multivariate conditional entropy measures for the characterization of short-term cardiovascular complexity under physiological stress. Physiol Meas 2018;39(1):014002.

[13] Goldberger AL, Peng C-K, Lipsitz LA. What is physiologic complexity and how does it change with aging and disease? Neurobiol Aging 2002;23(1):23-6.

[14] Shaffer F, McCraty R, Zerr CL. A healthy heart is not a metronome: an integrative review of the heart's anatomy and heart rate variability. Front Psychol 2014;5:1040.

[15] Lehrer P, Eddie D. Dynamic processes in regulation and some implications for biofeedback and biobehavioral interventions. Appl Psychophysiol biofeedback 2013;38:143-55.

[16] Javorka M, Krohova J, Czippelova B, Turianikova Z, Lazarova Z, Wiszt R, Faes L. Towards understanding the complexity of cardiovascular oscillations: Insights from information theory. Comput biol Med 2018;98:48-57.

[17] Javorka M, Turianikova Z, Tonhajzerova I, Javorka K, Baumert M. The effect of orthostasis on recurrence quantification analysis of heart rate and blood pressure dynamics. Physiol Meas 2008;30(1):29.

[18] Porta A, Bari V, Ranuzzi G, De Maria B, Baselli G. Assessing multiscale complexity of short heart rate variability series through a model-based linear approach. Chaos 2017;27(9).

[19] Heffernan KS, Fahs CA, Shinsako KK, Jae SY, Fernhall B. Heart rate recovery and heart rate complexity following resistance exercise training and detraining in young men. Am J Physiol-Heart Circ Physiol 2007;293(5):H3180-6.

[20] Takahashi AC, Porta A, Melo RC, Quitério RJ, da Silva E, Borghi-Silva A, Tobaldini E, Montano N, Catai AM. Aging reduces complexity of heart rate variability assessed by conditional entropy and symbolic analysis. Intern Emerg Med 2012;7:229-35.

[21] Porta A, Faes L, Bari V, Marchi A, Bassani T, Nollo G, Perseguini NM, Milan J, Minatel V, Borghi-Silva A, et al. Effect of age on complexity and causality of the cardiovascular control: comparison between model-based and model-free approaches. PLoS One 2014;9(2):e89463.

[22] Jia Y, Gu H, Luo Q. Sample entropy reveals an age-related reduction in the complexity of dynamic brain. Sci Rep 2017;7(1):7990.

[23] Pincus SM. Greater signal regularity may indicate increased system isolation. Math Biosci 1994;122(2):161-81.

[24] Romero-Ortuño R, Martínez-Velilla N, Sutton R, Ungar A, Fedorowski A, Galvin R, Theou O, Davies A, Reilly RB, Claassen J, et al. Network physiology in aging and frailty: the grand challenge of physiological reserve in older adults. 2021.

[25] Méndez MA, Zuluaga P, Hornero R, Gómez C, Escudero J, Rodríguez-Palancas A, Ortiz T, Fernández A. Complexity analysis of spontaneous brain activity: effects of depression and antidepressant treatment. J Psychopharmacol 2012;26(5):636-43.

[26] Fernández A, Quintero J, Hornero R, Zuluaga P, Navas M, Gómez C, Escudero J, García-Campos N, Biederman J, Ortiz T. Complexity analysis of spontaneous brain activity in attention-deficit/hyperactivity disorder: diagnostic implications. Biol Psychiatry 2009;65(7):571-7.

[27] Hornero R, Abásolo D, Jimeno N, Sánchez CI, Poza J, Aboy M. Variability, regularity, and complexity of time series generated by schizophrenic patients and control subjects. IEEE Trans Biomed Eng 2006;53(2):210-8.

[28] Chen C, Jin Y, Lo IL, Zhao H, Sun B, Zhao Q, Zheng J, Zhang XD. Complexity change in cardiovascular disease. Int J Biol Sci 2017;13(10):1320.

[29] Trunkvalterova Z, Javorka M, Tonhajzerova I, Javorkova J, Lazarova Z, Javorka K, Baumert M. Reduced short-term complexity of heart rate and blood pressure dynamics in patients with diabetes mellitus type 1: multiscale entropy analysis. Physiol Meas 2008;29(7):817.

[30] Tobaldini E, Nobili L, Strada S, Casali KR, Braghiroli A, Montano N. Heart rate variability in normal and pathological sleep. Front Physiol 2013;4:294.

[31] Theiler J. Estimating fractal dimension. J Opt Soc Amer A 1990;7(6):1055-73.

[32] Wolf A, Swift JB, Swinney HL, Vastano JA. Determining lyapunov exponents from a time series. Physica D 1985;16(3):285-317.

[33] Lempel A, Ziv J. On the complexity of finite sequences. IEEE Trans Inform Theory 1976;22(1):75-81.

[34] Baranger M, Institute NECS. Chaos, complexity, and entropy: A physics talk for non-pharmacists. New England Complex Systems Institute; 2001, URL https://books.google.sk/books?id=duzMgEACAAJ.

[35] Faes L, Porta A. Conditional entropy-based evaluation of information dynamics in physiological systems. Direct Inf Meas Neurosci 2014;61-86.

[36] Lizier JT. The local information dynamics of distributed computation in complex systems. Springer Science & Business Media; 2012.

[37] Faes L, Erla S, Nollo G. Measuring connectivity in linear multivariate processes: definitions, interpretation, and practical analysis. Comput Math Methods Med 2012;2012.

[38] Kozachenko LF, Leonenko NN. Sample estimate of the entropy of a random vector. Probl Peredachi Inform 1987;23(2):9-16.

[39] Richman JS, Moorman JR. Physiological time-series analysis using approximate entropy and sample entropy. Am J Physiol-Heart Circ Physiol 2000.

[40] Azami H, Faes L, Escudero J, Humeau-Heurtier A, Silva LE. Entropy analysis of univariate biomedical signals: Review and comparison of methods. Front Entropy across Discipl: Panorama Entropy: Theory Comput Appl 2023;233-86.

[41] Bandt C, Pompe B. Permutation entropy: a natural complexity measure for time series. Phys Rev Lett 2002;88(17):174102.

[42] Cuesta-Frau D. Slope entropy: A new time series complexity estimator based on both symbolic patterns and amplitude information. Entropy 2019;21(12):1167.

[43] Shannon CE. A mathematical theory of communication. Bell Syst Tech J 1948;27(3):379-423.

[44] Cover TM. Elements of information theory. John Wiley & Sons; 1999.

[45] Kolmogorov AN. A new metric invariant of transitive dynamical systems and automorphisms of lebesgue spaces. Trudy Mat Inst Imeni VA Steklova 1985;169:94-8.

[46] Sinai YG. On the notion of entropy of a dynamical system. In: Doklady of Russian academy of sciences. Vol. 124, 1959, p. 768-71.

[47] Porta A, De Maria B, Bari V, Marchi A, Faes L. Are nonlinear model-free conditional entropy approaches for the assessment of cardiac control complexity superior to the linear model-based one? IEEE Trans Biomed Eng 2016;64(6):1287-96.

[48] Xiong W, Faes L, Ivanov PC. Entropy measures, entropy estimators, and their performance in quantifying complex dynamics: Effects of artifacts, nonstationarity, and long-range correlations. Phys Rev E 2017;95(6):062114.

[49] Barà C, Sparacino L, Pernice R, Antonacci Y, Porta A, Kugiumtzis D, Faes L. Comparison of discretization strategies for the model-free information-theoretic assessment of short-term physiological interactions. Chaos 2023;33(3):033127.

[50] Runge J, Heitzig J, Petoukhov V, Kurths J. Escaping the curse of dimensionality in estimating multivariate transfer entropy. Phys Rev Lett 2012;108(25):258701.

[51] Porta A, Baselli G, Liberati D, Montano N, Cogliati C, Gneчи-Ruscone T, Malliani A, Cerutti S. Measuring regularity by means of a corrected conditional entropy in sympathetic outflow. Biol Cybern 1998;78:71-8.

[52] Barrett AB, Barnett L, Seth AK. Multivariate granger causality and generalized variance. Phys Rev E 2010;81(4):041907.

[53] Mardia KV, Marshall RJ. Maximum likelihood estimation of models for residual covariance in spatial regression. Biometrika 1984;71(1):135-46.

[54] Akaike H. A new look at the statistical model identification. IEEE Trans Autom Control 1974;19(6):716-23.

[55] Schwarz G. Estimating the dimension of a model. Ann Stat 1978;461-4.

[56] Barnett L, Barrett AB, Seth AK. Granger causality and transfer entropy are equivalent for gaussian variables. Phys Rev Lett 2009;103(23):238701.

[57] Kraskov A, Stögbauer H, Grassberger P. Estimating mutual information. Phys Rev E 2004;69(6):066138.
[58] Pincus SM. Approximate entropy as a measure of system complexity. Proc Natl Acad Sci 1991;88(6):2297-301.

[59] Delgado-Bonal A, Marshak A. Approximate entropy and sample entropy: A comprehensive tutorial. Entropy 2019;21(6):541.

[60] Unakafov AM, Keller K. Conditional entropy of ordinal patterns. Physica D 2014;269:94-102.

[61] Kugiumtzis D. Partial transfer entropy on rank vectors. Eur Phys J Spec Top 2013;222(2):401-20.

[62] Pernice R, Javorka M, Krohova J, Czippelova B, Turianikova Z, Busacca A, Faes L, Member I. Comparison of short-term heart rate variability indexes evaluated through electrocardiographic and continuous blood pressure monitoring. Med Biol Eng Comput 2019;57:1247-63.

[63] Faes L, Porta A, Nollo G. Information decomposition in bivariate systems: theory and application to cardiorespiratory dynamics. Entropy 2015;17(1):277-303.

[64] Cover TM, Thomas JA. Differential entropy. Elem Inf Theory 1991;224-38.

[65] Santamaría-Bonfil G, Fernández N, Gershenson C. Measuring the complexity of continuous distributions. Entropy 2016;18(3):72.

[66] May RM. Simple mathematical models with very complicated dynamics. Nature 1976;261:459-67.

[67] Schreiber T, Schmitz A. Improved surrogate data for nonlinearity tests. Phys Rev Lett 1996;77(4):635.

[68] Porta A, Bassani T, Bari V, Pinna GD, Maestri R, Guzzetti S. Accounting for respiration is necessary to reliably infer granger causality from cardiovascular variability series. IEEE Trans Biomed Eng 2011;59(3):832-41.

[69] Cairo B, Bari V, Gelpi F, De Maria B, Porta A. Assessing cardiorespiratory interactions via lagged joint symbolic dynamics during spontaneous and controlled breathing. Front Newt Physiol 2003;3.

[70] Darbellay GA, Wuertz D. The entropy as a tool for analysing statistical dependences in financial time series. Phys A 2000;287(3-4):429-39.

[71] Zhou R, Cai R, Tong G. Applications of entropy in finance: A review. Entropy 2013;15(11):4909-31.

[72] Leung L-Y, North GR. Information theory and climate prediction. J Clim 1990;3(1):5-14.

[73] Hlinka J, Hartman D, Vejmelka M, Runge J, Marwan N, Kurths J, Paluš M. Reliability of inference of directed climate networks using conditional mutual information. Entropy 2013;15(6):2023-45.

[74] Berntson GG, Cacioppo JT, Quigley KS. Respiratory sinus arrhythmia: Automatic origins, physiological mechanisms, and psychophysiological implications. Psychophysiology 1993;30(2):183-96.

[75] Saul JP, Berger RD, Chen M, Cohen RJ. Transfer function analysis of autonomic regulation, ii. respiratory sinus arrhythmia. Am J Physiol-Heart Circ Physiol 1989;256(1):H153-61.

[76] Tzeng Y, Larsen P, Galletly D. Cardioventilatory coupling in resting human subjects. Exp Physiol 2003;88(6):775-82.

[77] Elstad M, O'Callaghan EL, Smith AJ, Ben-Tal A, Ramchandra R. Cardiorespiratory interactions in humans and animals: rhythms for life. Am J Physiol-Heart Circ Physiol 2018;315(1):H6-17.

[78] Porta A, Guzzetti S, Furlan R, Gneccchi-Ruscone T, Montano N, Malliani A. Complexity and nonlinearity in short-term heart period variability: comparison of methods based on local nonlinear prediction. IEEE Trans Biomed Eng 2006;54(1):94-106.

[79] Porta A, Bari V, Gelpi F, Cairo B, De Maria B, Tonon D, Rossato G, Faes L. On the different abilities of cross-sample entropy and k-nearest-neighbor cross-unpredictability in assessing dynamic cardiorespiratory and cerebrovascular interactions. Entropy 2023;25(4):599.

[80] Nuzzi D, Stramaglia S, Javorka M, Marinazzo D, Porta A, Faes L. Extending the spectral decomposition of granger causality to include instantaneous influences: application to the control mechanisms of heart rate variability. Phil Trans R Soc A 2021;379(2212):20200263.

[81] Volpes G, Barà C, Busacca A, Stivala S, Javorka M, Faes L, Pernice R. Feasibility of ultra-short-term analysis of heart rate and systolic arterial pressure variability at rest and during stress via time-domain and entropy-based measures. Sensors 2022;22(23):9149.

[82] Lombardi D, Pant S. Nonparametric k-nearest-neighbor entropy estimator. Phys Rev E 2016;93(1):013310.

[83] Trujillo LT. K-th nearest neighbor (knn) entropy estimates of complexity and integration from ongoing and stimulus-evoked electroencephalographic (EEg) recordings of the human brain. Entropy 2019;21(1):61.

[84] Kaplan DT, Furman MI, Pincus SM, Ryan SM, Lipsitz LA, Goldberger AL. Aging and the complexity of cardiovascular dynamics. Biophys J 1991;59(4):945-9.

[85] Lewis M, Short A. Sample entropy of electrocardiographic rr and qt time-series data during rest and exercise. Physiol Meas 2007;28(6):731.

[86] Lake DE, Richman JS, Griffin MP, Moorman JR. Sample entropy analysis of neonatal heart rate variability. Am J Physiol-Regul Integr Comp Physiol 2002;283(3):R789-97.

[87] Porta A, Gneccchi-Ruscone T, Tobaldini E, Guzzetti S, Furlan R, Montano N. Progressive decrease of heart period variability entropy-based complexity during graded head-up tilt. J Appl Physiol 2007;103(4):1143-9.

[88] Porta A, Faes L, Masé M, D'addio G, Pinna G, Maestri R, Montano N, Furlan R, Guzzetti S, Nollo G, et al. An integrated approach based on uniform quantization for the evaluation of complexity of short-term heart period variability: application to 24h holter recordings in healthy and heart failure humans. Chaos 2007;17(1).

[89] Cuesta-Frau D, Miró-Martínez P, Oltra-Crespo S, Molina-Picó A, Dakappa PH, Mahabala C, Vargas B, González P. Classification of fever patterns using a single extracted entropy feature: A feasibility study based on sample entropy. Math Biosci Eng 2020;17:235.

[90] Morales J, Borzée P, Testelmans D, Buyse B, Van Huffel S, Varon C. Linear and non-linear quantification of the respiratory sinus arrhythmia using support vector machines. Front Physiol 2021;12:623781.

[91] Porta A, Guzzetti S, Montano N, Furlan R, Pagani M, Malliani A, Cerutti S. Entropy, entropy rate, and pattern classification as tools to typify complexity in short heart period variability series. IEEE Trans Biomed Eng 2001;48(11):1282-91.

[92] Arya S, Mount DM, Netanyahu NS, Silverman R, Wu AY. An optimal algorithm for approximate nearest neighbor searching fixed dimensions. J ACM 1998;45(6):891-923.

[93] Merkwirth C, Parlitz U, Luterborn W. Fast nearest-neighbor searching for nonlinear signal processing. Phys Rev E 2000;62(2):2089.

[94] Samet H. K-nearest neighbor finding using maxnearestdist. IEEE Trans Pattern Anal Mach Intell 2007;30(2):243-52.

[95] Manis G. Fast computation of approximate entropy. Comput Methods Programs Biomed 2008;91(1):48-54.

[96] Pan Y-H, Wang Y-H, Liang S-F, Lee K-T. Fast computation of sample entropy and approximate entropy in biomedicine. Comput Methods Programs Biomed 2011;104(3):382-96.

[97] Manis G, Aktaruzzaman M, Sassi R. Low computational cost for sample entropy. Entropy 2018;20(1):61.

[98] Liu W, Jiang Y, Xu Y. A super fast algorithm for estimating sample entropy. Entropy 2022;24(4):524.

[99] Dias D, Paulo Silva Cunha J. Wearable health devices—vital sign monitoring, systems and technologies. Sensors 2018;18(8):2414.

[100] Kakria P, Tripathi N, Kitipawang P. A real-time health monitoring system for remote cardiac patients using smartphone and wearable sensors. Int J Telegmed Appl 2015;2015:8.

[101] Georgiou K, Larentzakis AV, Khamis NN, Alsuhaibani GI, Alaska YA, Gialafos EJ. Can wearable devices accurately measure heart rate variability? a systematic review. Folia Med 2018;60(1):7-20.

[102] Pernice R, Javorka M, Krohova J, Czippelova B, Turianikova Z, Busacca A, Faes L. A validity and reliability study of conditional entropy measures of pulse rate variability. In: 2019 41st annual international conference of the IEEE engineering in medicine and biology society. EMBC, IEEE; 2019, p. 5568-71.

[103] Citi L, Guffanti G, Mainardi L. Rank-based multi-scale entropy analysis of heart rate variability. In: Computing in cardiology 2014. IEEE; 2014, p. 597-600.

[104] Manis G, Aktaruzzaman M, Sassi R. Bubble entropy: An entropy almost free of parameters. IEEE Trans Biomed Eng 2017;64(11):2711-8.

[105] Omidvarnia A, Mesbah M, Pedersen M, Jackson G. Range entropy: A bridge between signal complexity and self-similarity. Entropy 2018;20(12):962.

[106] Wang X, Si S, Li Y. Multiscale diversity entropy: A novel dynamical measure for fault diagnosis of rotating machinery. IEEE Trans Ind Inf 2020;17(8):5419-29.

[107] Chen W, Zhuang J, Yu W, Wang Z. Measuring complexity using fuzzyen, apen, and sampa. Med Eng Phys 2009;31(1):61-8.

[108] Li P, Liu C, Li K, Zheng D, Liu C, Hou Y. Assessing the complexity of short-term heartbeat interval series by distribution entropy. Med Biol Eng Comput 2015;53:77-87.

[109] Rostaghi M, Azami H. Dispersion entropy: A measure for time-series analysis. IEEE Signal Process Lett 2016;23(5):610-4.

[110] Rohila A, Sharma A. Phase entropy: A new complexity measure for heart rate variability. Physiol Meas 2019;40(10):105006.

[111] Platiša MM, Radovanović NN, Pernice R, Barà C, Pavlović SU, Faes L. Information-theoretic analysis of cardio-respiratory interactions in heart failure patients: Effects of arrhythmias and cardiac resynchronization therapy. Entropy 2023;25(7):1072.

[112] Nollo G, Faes L, Antolini R, Porta A. Assessing causality in normal and impaired short-term cardiovascular regulation via nonlinear prediction methods. Phil Trans R Soc A 2009;367(1892):1423-40.

[113] Cerutti S, Corino VD, Mainardi L, Lombardi F, Aktaruzzaman M, Sassi R. Non-linear regularity of arterial blood pressure variability in patient with atrial fibrillation in tilt-test procedure. Europea 2014;16(suppl_4):iv141-7.

[114] Chen Y, Yang H. Multiscale recurrence analysis of long-term nonlinear and nonstationary time series. Chaos Solitons Fractals 2012;45(7):978-87.

[115] Faes L, Porta A, Nollo G, Javorka M. Information decomposition in multi-variate systems: definitions, implementation and application to cardiovascular networks. Entropy 2016;19(1):5.
[58] Pincus SM. Approximate entropy as a measure of system complexity. Proc Natl Acad Sci 1991;88(6):2297-301.

[59] Delgado-Bonal A, Marshak A. Approximate entropy and sample entropy: A comprehensive tutorial. Entropy 2019;21(6):541.

[60] Unakafov AM, Keller K. Conditional entropy of ordinal patterns. Physica D 2014;269:94-102.

[61] Kugiumtzis D. Partial transfer entropy on rank vectors. Eur Phys J Spec Top 2013;222(2):401-20.

[62] Pernice R, Javorka M, Krohova J, Czippelova B, Turianikova Z, Busacca A, Faes L, Member I. Comparison of short-term heart rate variability indexes evaluated through electrocardiographic and continuous blood pressure monitoring. Med Biol Eng Comput 2019;57:1247-63.

[63] Faes L, Porta A, Nollo G. Information decomposition in bivariate systems: theory and application to cardiorespiratory dynamics. Entropy 2015;17(1):277-303.

[64] Cover TM, Thomas JA. Differential entropy. Elem Inf Theory 1991;224-38.

[65] Santamaría-Bonfil G, Fernández N, Gershenson C. Measuring the complexity of continuous distributions. Entropy 2016;18(3):72.

[66] May RM. Simple mathematical models with very complicated dynamics. Nature 1976;261:459-67.

[67] Schreiber T, Schmitz A. Improved surrogate data for nonlinearity tests. Phys Rev Lett 1996;77(4):635.

[68] Porta A, Bassani T, Bari V, Pinna GD, Maestri R, Guzzetti S. Accounting for respiration is necessary to reliably infer granger causality from cardiovascular variability series. IEEE Trans Biomed Eng 2011;59(3):832-41.

[69] Cairo B, Bari V, Gelpi F, De Maria B, Porta A. Assessing cardiorespiratory interactions via lagged joint symbolic dynamics during spontaneous and controlled breathing. Front Newt Physiol 2003;3.

[70] Darbellay GA, Wuertz D. The entropy as a tool for analysing statistical dependences in financial time series. Phys A 2000;287(3-4):429-39.

[71] Zhou R, Cai R, Tong G. Applications of entropy in finance: A review. Entropy 2013;15(11):4909-31.

[72] Leung L-Y, North GR. Information theory and climate prediction. J Clim 1990;3(1):5-14.

[73] Hlinka J, Hartman D, Vejmelka M, Runge J, Marwan N, Kurths J, Paluš M. Reliability of inference of directed climate networks using conditional mutual information. Entropy 2013;15(6):2023-45.

[74] Berntson GG, Cacioppo JT, Quigley KS. Respiratory sinus arrhythmia: Automatic origins, physiological mechanisms, and psychophysiological implications. Psychophysiology 1993;30(2):183-96.

[75] Saul JP, Berger RD, Chen M, Cohen RJ. Transfer function analysis of autonomic regulation, ii. respiratory sinus arrhythmia. Am J Physiol-Heart Circ Physiol 1989;256(1):H153-61.

[76] Tzeng Y, Larsen P, Galletly D. Cardioventilatory coupling in resting human subjects. Exp Physiol 2003;88(6):775-82.

[77] Elstad M, O'Callaghan EL, Smith AJ, Ben-Tal A, Ramchandra R. Cardiorespiratory interactions in humans and animals: rhythms for life. Am J Physiol-Heart Circ Physiol 2018;315(1):H6-17.

[78] Porta A, Guzzetti S, Furlan R, Gneccchi-Ruscone T, Montano N, Malliani A. Complexity and nonlinearity in short-term heart period variability: comparison of methods based on local nonlinear prediction. IEEE Trans Biomed Eng 2006;54(1):94-106.

[79] Porta A, Bari V, Gelpi F, Cairo B, De Maria B, Tonon D, Rossato G, Faes L. On the different abilities of cross-sample entropy and k-nearest-neighbor cross-unpredictability in assessing dynamic cardiorespiratory and cerebrovascular interactions. Entropy 2023;25(4):599.

[80] Nuzzi D, Stramaglia S, Javorka M, Marinazzo D, Porta A, Faes L. Extending the spectral decomposition of granger causality to include instantaneous influences: application to the control mechanisms of heart rate variability. Phil Trans R Soc A 2021;379(2212):20200263.

[81] Volpes G, Barà C, Busacca A, Stivala S, Javorka M, Faes L, Pernice R. Feasibility of ultra-short-term analysis of heart rate and systolic arterial pressure variability at rest and during stress via time-domain and entropy-based measures. Sensors 2022;22(23):9149.

[82] Lombardi D, Pant S. Nonparametric k-nearest-neighbor entropy estimator. Phys Rev E 2016;93(1):013310.

[83] Trujillo LT. K-th nearest neighbor (knn) entropy estimates of complexity and integration from ongoing and stimulus-evoked electroencephalographic (EEg) recordings of the human brain. Entropy 2019;21(1):61.

[84] Kaplan DT, Furman MI, Pincus SM, Ryan SM, Lipsitz LA, Goldberger AL. Aging and the complexity of cardiovascular dynamics. Biophys J 1991;59(4):945-9.

[85] Lewis M, Short A. Sample entropy of electrocardiographic rr and qt time-series data during rest and exercise. Physiol Meas 2007;28(6):731.

[86] Lake DE, Richman JS, Griffin MP, Moorman JR. Sample entropy analysis of neonatal heart rate variability. Am J Physiol-Regul Integr Comp Physiol 2002;283(3):R789-97.

[87] Porta A, Gneccchi-Ruscone T, Tobaldini E, Guzzetti S, Furlan R, Montano N. Progressive decrease of heart period variability entropy-based complexity during graded head-up tilt. J Appl Physiol 2007;103(4):1143-9.

[88] Porta A, Faes L, Masé M, D'addio G, Pinna G, Maestri R, Montano N, Furlan R, Guzzetti S, Nollo G, et al. An integrated approach based on uniform quantization for the evaluation of complexity of short-term heart period variability: application to 24h holter recordings in healthy and heart failure humans. Chaos 2007;17(1).

[89] Cuesta-Frau D, Miró-Martínez P, Oltra-Crespo S, Molina-Picó A, Dakappa PH, Mahabala C, Vargas B, González P. Classification of fever patterns using a single extracted entropy feature: A feasibility study based on sample entropy. Math Biosci Eng 2020;17:235.

[90] Morales J, Borzée P, Testelmans D, Buyse B, Van Huffel S, Varon C. Linear and non-linear quantification of the respiratory sinus arrhythmia using support vector machines. Front Physiol 2021;12:623781.

[91] Porta A, Guzzetti S, Montano N, Furlan R, Pagani M, Malliani A, Cerutti S. Entropy, entropy rate, and pattern classification as tools to typify complexity in short heart period variability series. IEEE Trans Biomed Eng 2001;48(11):1282-91.

[92] Arya S, Mount DM, Netanyahu NS, Silverman R, Wu AY. An optimal algorithm for approximate nearest neighbor searching fixed dimensions. J ACM 1998;45(6):891-923.

[93] Merkwirth C, Parlitz U, Luterborn W. Fast nearest-neighbor searching for nonlinear signal processing. Phys Rev E 2000;62(2):2089.

[94] Samet H. K-nearest neighbor finding using maxnearestdist. IEEE Trans Pattern Anal Mach Intell 2007;30(2):243-52.

[95] Manis G. Fast computation of approximate entropy. Comput Methods Programs Biomed 2008;91(1):48-54.

[96] Pan Y-H, Wang Y-H, Liang S-F, Lee K-T. Fast computation of sample entropy and approximate entropy in biomedicine. Comput Methods Programs Biomed 2011;104(3):382-96.

[97] Manis G, Aktaruzzaman M, Sassi R. Low computational cost for sample entropy. Entropy 2018;20(1):61.

[98] Liu W, Jiang Y, Xu Y. A super fast algorithm for estimating sample entropy. Entropy 2022;24(4):524.

[99] Dias D, Paulo Silva Cunha J. Wearable health devices—vital sign monitoring, systems and technologies. Sensors 2018;18(8):2414.

[100] Kakria P, Tripathi N, Kitipawang P. A real-time health monitoring system for remote cardiac patients using smartphone and wearable sensors. Int. J Telegmed Appl 2015;2015:8.

[101] Georgiou K, Larentzakis AV, Khamis NN, Alsuhaibani GI, Alaska YA, Gialafos EJ. Can wearable devices accurately measure heart rate variability? a systematic review. Folia Med 2018;60(1):7-20.

[102] Pernice R, Javorka M, Krohova J, Czippelova B, Turianikova Z, Busacca A, Faes L. A validity and reliability study of conditional entropy measures of pulse rate variability. In: 2019 41st annual international conference of the IEEE engineering in medicine and biology society. EMBC, IEEE; 2019, p. 5568-71.

[103] Citi L, Guffanti G, Mainardi L. Rank-based multi-scale entropy analysis of heart rate variability. In: Computing in cardiology 2014. IEEE; 2014, p. 597-600.

[104] Manis G, Aktaruzzaman M, Sassi R. Bubble entropy: An entropy almost free of parameters. IEEE Trans Biomed Eng 2017;64(11):2711-8.

[105] Omidvarnia A, Mesbah M, Pedersen M, Jackson G. Range entropy: A bridge between signal complexity and self-similarity. Entropy 2018;20(12):962.

[106] Wang X, Si S, Li Y. Multiscale diversity entropy: A novel dynamical measure for fault diagnosis of rotating machinery. IEEE Trans Ind Inf 2020;17(8):5419-29.

[107] Chen W, Zhuang J, Yu W, Wang Z. Measuring complexity using fuzzyen, apen, and sampa. Med Eng Phys 2009;31(1):61-8.

[108] Li P, Liu C, Li K, Zheng D, Liu C, Hou Y. Assessing the complexity of short-term heartbeat interval series by distribution entropy. Med Biol Eng Comput 2015;53:77-87.

[109] Rostaghi M, Azami H. Dispersion entropy: A measure for time-series analysis. IEEE Signal Process Lett 2016;23(5):610-4.

[110] Rohila A, Sharma A. Phase entropy: A new complexity measure for heart rate variability. Physiol Meas 2019;40(10):105006.

[111] Platiša MM, Radovanović NN, Pernice R, Barà C, Pavlović SU, Faes L. Information-theoretic analysis of cardio-respiratory interactions in heart failure patients: Effects of arrhythmias and cardiac resynchronization therapy. Entropy 2023;25(7):1072.

[112] Nollo G, Faes L, Antolini R, Porta A. Assessing causality in normal and impaired short-term cardiovascular regulation via nonlinear prediction methods. Phil Trans R Soc A 2009;367(1892):1423-40.

[113] Cerutti S, Corino VD, Mainardi L, Lombardi F, Aktaruzzaman M, Sassi R. Non-linear regularity of arterial blood pressure variability in patient with atrial fibrillation in tilt-test procedure. Europea 2014;16(suppl_4):iv141-7.

[114] Chen Y, Yang H. Multiscale recurrence analysis of long-term nonlinear and nonstationary time series. Chaos Solitons Fractals 2012;45(7):978-87.

[115] Faes L, Porta A, Nollo G, Javorka M. Information decomposition in multi-variate systems: definitions, implementation and application to cardiovascular networks. Entropy 2016;19(1):5.