# Efficient Computation of Multiscale Entropy over Short Biomedical Time Series Based on Linear State-Space Models

Authors: Luca Faes; Alberto Porta; Michal Javorka; Giandomenico Nollo
Journal: Complexity
 (2017)

---

## Abstract

The most common approach to assess the dynamical complexity of a time series across multiple temporal scales makes use of the multiscale entropy (MSE) and refined MSE (RMSE) measures. In spite of their popularity, MSE and RMSE lack an analytical framework allowing their calculation for known dynamic processes and cannot be reliably computed over short time series. To overcome these limitations, we propose a method to assess RMSE for autoregressive (AR) stochastic processes. The method makes use of linear state-space (SS) models to provide the multiscale parametric representation of an AR process observed at different time scales and exploits the SS parameters to quantify analytically the complexity of the process. The resulting linear MSE (LMSE) measure is first tested in simulations, both theoretically to relate the multiscale complexity of AR processes to their dynamical properties and over short process realizations to assess its computational reliability in comparison with RMSE. Then, it is applied to the time series of heart period, arterial pressure, and respiration measured for healthy subjects monitored in resting conditions and during physiological stress. This application to short-term cardiovascular variability documents that LMSE can describe better than RMSE the activity of physiological mechanisms producing biological oscillations at different temporal scales.


---

## Body

## 1. Introduction

An intrinsic feature of almost all physiological systems, which is clearly visible in the time course of the variables measured from these systems, is their dynamical complexity. It is indeed widely acknowledged that physiological systems such as the brain, the cardiovascular system, and the muscular system produce output signals that exhibit highly complex dynamics, resulting from the combined activity of several mechanisms of physiological regulation which are coupled with each other through structural and functional pathways [1–4]. Since these multiple and simultaneously active mechanisms typically operate across multiple temporal scales, a surge of interest has emerged in the last two decades in computational methods able to connect the complexity of a dynamic oscillation to its temporal scale. The research in this context has been driven by the work of Costa et al. [5], who introduced multiscale entropy (MSE) as a measure of the complexity of a time series evaluated as a function of the time scale at which the series is observed. Since its introduction, MSE has been successfully employed in several fields of science [6], becoming a prevailing method to quantify the complexity of biomedical time series [7, 8] and gaining particular popularity in the analysis of brain [9, 10] and cardiovascular variability [11, 12] signals.In spite of its acknowledged usefulness, MSE has been shown to present some shortcomings which have led many authors to propose improvements and modifications of the original algorithm [6]. The main problems identified in the original MSE formulation are related to both its defining algorithmic steps, the rescaling procedure that changes the temporal scale of the observed series progressively filtering out the shorter scales and keeping the longer ones, and the computation of the entropy rate of the rescaled time series performed by means of the sample entropy (SampEn) metric [13]. Specifically, it has been reported [11, 14] that the rescaling procedure makes use of a suboptimal low-pass filter (i.e., an averaging filter) which cannot prevent aliasing and that the progressive application of SampEn employs a badly designed coarse-graining step which makes MSE artificially decrease as a function of the time scale. These problems were overcome by the introduction of the so-called refined MSE (RMSE) [11], which exploits well-designed procedures for low-pass filtering and coarse graining.Even though RMSE and other refinements and extensions [6] have improved the performance and widened the applicability of this methodology, some limitations still remain which hamper a full exploitation of the concept underlying MSE. The main limitation is in the fact that MSE requires long time series to be reliably computed: as any nonparametric entropy estimator, the SampEn is highly data-demanding, and the issue is exacerbated at increasing time scales where the rescaled signals get progressively shorter. This issue is only partly addressed by versions of MSE that intensify the number of patterns over which SampEn is calculated [15, 16], because the problem ultimately stands in the lower number of cycles of the slowest oscillations that are available for a given data length. Another limitation is in the fact that the length of the temporal scales which can be explored is usually expressed in number of samples and is restricted to integer values, thus limiting the possibility of fine-tuning the elimination of the fast temporal scales. These issues strongly limit the applicability of MSE or RMSE to short biomedical time series, such as those typically considered in short-term cardiovascular variability, where the dynamics of interest are deployed over a few minutes (~300 samples) [17]. Moreover, in this applicative context fine-tuning of the filtering process is fundamental to elicit changes in complexity related to the oscillatory rhythms typically observed in cardiovascular variability, that is, low-frequency (LF, from 0.04 Hz to 0.15 Hz) and high-frequency (HF, from 0.15 Hz to 0.5 Hz) oscillations [18, 19]. A further issue with existing MSE techniques is the lack of theoretical procedures to derive exact values of the multiscale complexity for known dynamical processes which one can study analytically and can fully control; though being theoretical, this issue has practical implications as analytical methods would allow one to assess the estimation bias of existing MSE measures.The present study introduces a new approach for the evaluation of multiscale complexity that is explicitly designed to address the limitations of existing MSE methods described above. The approach builds on recent theoretical advances providing exact techniques for the analytical computation of information-theoretic measures, including entropy rates, for linear autoregressive (AR) stochastic processes [20–23]. The key steps of the derivation of the new measure of multiscale complexity, which we denote as linear MSE (LMSE), are the closed-form representation of filtered and downsampled AR processes as state-space (SS) processes and the analytical quantification of complexity from SS parameters. Following these steps, we devise a procedure which allows an exact computation of the refined version of MSE for linear processes, starting from their AR parameters and from the desired scale factor. Additionally, a technique to set the scale factor at any rational number instead of an integer number is devised, thus allowing the fine-tuning of the filtering step of MSE computation. The proposed approach is tested on simulations of linear stochastic processes, first in a theoretical analysis aimed at relating multiscale complexity to the dynamical properties of the process without the unavoidable errors connected with estimation procedures and then in realizations of these simulated processes generated to assess the computational reliability of the LMSE compared with the traditional RMSE estimator. Moreover, the comparative ability of LMSE and RMSE in assessing the multiscale complexity of short-term cardiovascular variability is assessed on the beat-to-beat time series of the heart period, arterial pressure, and respiration measured in a large group of healthy subjects resting in a relaxed state and during two commonly studied conditions of physiological stress, that is, postural stress and mental stress.


## 2. Approaches to Multiscale Entropy Analysis

### 2.1. Multiscale Entropy

Multiscale entropy (MSE), originally proposed by Costa et al. [5], is a measure assessing the complexity of a process across multiple temporal scales. Its computation relies first on rescaling the observed process (i.e., focusing on a specific range of temporal scales) and then on assessing the dynamical complexity of the rescaled process through the computation of its rate of entropy production.Specifically, let us consider a discrete-time, stationary stochastic processx with zero mean and variance σx2. Let us further define x(n) as the variable obtained sampling the process at the discrete time n and set an integer scale factor τ. The rescaling step is aimed at eliminating the fast temporal scales and is performed applying the following transformation to the original process x, which averages τ consecutive samples to yield the rescaled process xd:(1)xdn=1τ∑i=0τ-1xnτ-i; the rescaled process is denoted here as xd because it is a downsampled version of the original process x; in fact, in (1) one sample of the process xd is obtained by taking one sample every τ values of the averaged version of x, where averaging is performed over τ consecutive samples. The second part of the procedure is implemented by calculating the conditional entropy of the present state of the rescaled process, xd(n), given its past states, xd–(n)=[xdn-1,xdn-2,…]. In MSE, this is accomplished approximating the past history of the rescaled process with the finite-length pattern xdm(n)=[xd(n-1),…,xd(n-m)] and then on estimating the conditional entropy by means of the Sample Entropy (SampEn) metric [13]. SampEn is an estimate of the probability that patterns which are detected as similar in the m-dimensional space remain similar also in the (m+1)-dimensional space when they are incremented with their future values, that is, the probability that xdm+1(n1+1) and xdm+1(n2+1) are similar given that xdm(n1) and xdm(n2) are similar; similarity is assessed through a coarse-graining procedure which applies the Heaviside step function with parameter r on the distance between patterns (i.e., r is a threshold such that xdm(n1) and xdm(n2) are similar if their distance in the m-dimensional space is lower than r, where the distance is assessed using the maximum norm).The MSE measure resulting from the above procedure, which we denote asCx(τ), quantifies the dynamical complexity of the original process x observed at scale τ. The free parameters of the MSE estimator are the length m of the patterns used to approximate the past of the process and the threshold r that sets the similarity between patterns. In the application of MSE, the pattern length is limited to a few samples (typically, m=2), while the threshold distance is a fraction of the standard deviation of the original process (typically, r=0.2·σx).


### 2.2. Refined Multiscale Entropy

The original MSE formulation suffers from two main limitations: the suboptimal procedure for elimination of the fast temporal scales implemented by (1), which tends to introduce spurious oscillations in the rescaled time series, and the fact that the threshold parameter of the coarse graining is kept at a constant value for all time scales, which brings about an artificial reduction of the MSE values with increasing scales. Valencia et al. [11] proposed a refined MSE (RMSE) measure that is aimed at overcoming these limitations. The solution to the first limitation is based on the rationale that the rescaling procedure actually consists of two steps: a filtering step which eliminates the fast temporal scales from the original process and a downsampling step that eliminates the redundancy introduced by the first step. Accordingly, the filtering step was designed to limit as much as possible the aliasing that can occur with the following downsampling step.In fact, the change of scale of a processx is performed first applying a low-pass filter to obtain the filtered process xf(n) and then reducing the sampling rate of the filtered process to obtain the rescaled (downsampled) processxd(n). The two steps yield, respectively, the processes(2a)xfn=∑i=0qbixn-i-∑j=1rcjxfn-j,(2b)xdn=xfnτ,where b(i) and c(j) are the filter coefficients and q and r are the filter orders. As it averages τ consecutive samples, the original MSE formulation [5] implicitly applies a finite impulse response (FIR) low-pass filter [24] of orders q=τ–1 and r=0 (the filter is nonrecursive) and with coefficients b(i)=1/τ for each i=0,1,…,q=τ–1; the cutoff frequency of the filter is constrained to the value fτ=1/(2τ). To improve elimination of the fast temporal scales and satisfy the Shannon theorem in the subsequent downsampling step, the RMSE approach implements an infinite impulse response (IIR) avoiding ripples in the stop band; specifically, a Butterworth filter [24] of order q (with r=q-1) is implemented in which the coefficients b(i) and c(j) are set to realize a low-pass implementation with cutoff frequency fτ [11]. This choice brings also the advantage that since the cutoff frequency fτ can take any real value, in RMSE the scale factor τ=1/(2fτ) is not constrained to integer numbers as in MSE; this allows filtering out the oscillatory components of the original process with a better resolution before computing the complexity of the rescaled process.Besides the type of filter, another crucial difference exists between the original and refined formulations of MSE. Whereas in MSE the parameterr is constant for all scale factors, as it is set at a fraction of the standard deviation of the original process (e.g., r=0.2·σx), in RMSE this parameter is set at a fraction of the standard deviation of the rescaled process (e.g., r=0.2·σxd). This choice is meant to reduce the dependence of the estimated conditional entropy on the decrease of variance due to filtering. As a consequence of this refinement, RMSE does not exhibit a tendency to decline with the time scale as a consequence of the inherent reduction of the variance of the filtered process: for example, for a white noise process MSE decreases with the time scale while RMSE is constant.


### 2.3. Linear Multiscale Entropy

In this section we present a method to assess the multiscale complexity of linear Gaussian stochastic processes. The method is based on the fact that if the variables obtained sampling the considered processx have a joint Gaussian distribution, all the variability that defines the entropy rate of the process is fully captured by a linear autoregressive (AR) model [25]. Accordingly, let us consider the linear AR representation of the process x, defined as(3)xn=∑k=1pakxn-k+en,where p is the order of the process, a(k), k=1,…,p, are linear regression coefficients describing the interactions between the process variables as a function of the lag k, and e(n) is an innovation process with variance σe2. The process x is fully described by the AR model (3) when the innovation processe is formed by uncorrelated Gaussian variables.Given the AR representation, the variance of the process and of the innovations can be used to derive an information-theoretic description of the statistical structure of the process [26]. Specifically, the entropy of the process is related to its variance by [27](4)Hx=12ln⁡2πeσx2, and the entropy rate of the process, that is, the conditional entropy of the present x(n) given the past x–(n)=[xn-1,xn-2,…], is related to the variance of the innovations by [28](5)Hx∣x-=12ln⁡2πeσe2. Then, (4) and (5) can be combined to provide a version of the conditional entropy which quantifies the dynamical complexity of the normalized process as(6)Cx=Hx∣x--Hx+12ln⁡2πe=12ln⁡2πeσe2σx2; note that applying (6) corresponds to computing the conditional entropy of the original process after normalizing it to unit variance.Now we turn to show how to compute analytically the variance of the AR innovations after rescaling the original process, in a way such that this variance can be used as in (6) to assess multiscale complexity. First, we perform a virtual upsampling of the original process x, by setting an integer upsampling factor s and by defining the process(7)xun=∑k=1psaukxun-k+eun,where ps=p·s is the order of the upsampled process and the coefficients au are set to be a zero-padded version of the original AR coefficients; that is, we set au(ls)=a(l) for each l=1,…,p and au(k)=0 for each k=1,…,ps, k≠ls. Note that the innovations for the original and upsampled processes are formally the same, that is, σeu2=σe2. Then, we filter the upsampled process using a linear FIR filter, that is, a filter implemented as in (2a) with c(j)=0 for each j; with some algebraic manipulation we find that, substituting (7) in (2a) applied with r=0 and with xu in place of x, the filtered representation of the upsampled process takes the form(8)xfn=∑k=1psaukxfn-k+∑i=0qbieun-i.Equation (8) shows that the filtering step introduces a moving average (MA) component in the original AR process, transforming it into an ARMA process [21]. Then, we exploit the connection between ARMA processes and state-space (SS) processes [29] to evidence that the ARMA process (8) can be expressed in an SS form as(9)Zfn+1=Af·Zfn+Kfefn,xfn=Cf·Zfn+efn,where Zf(n)=[xf(n-1)⋯xf(n-ps)u(n-1)⋯u(n-q)]T is a (ps+q)-dimensional state process, ef(n) is the scalar SS innovation process defined as ef(n)=b0eu(n), and the vectors Cf and Kf and the matrix Af are defined as(10)Af=au1⋯aups-1aupsb1⋯bq-1bq1⋯000⋯00⋮⋱⋮⋮⋮⋱⋮⋮0⋯100⋯000⋯000⋯000⋯001⋯000⋯0⋮⋮⋱⋮⋮0⋯000⋯10.Kf=10⋮0b00⋮0Cf=au1⋯aups-1aupsb1⋯bq-1bq.Such a connection between ARMA and SS models allows finding, through (10) and the relation σef2=b(0)2σeu2, the parameters of the SS process descriptive of the process xf. Such a process is a filtered version of the upsampled process xu, obtained passing xu through a FIR low-pass filter with cutoff frequency equal to fτ=1/(2τ); therefore, since the sampling frequency of the upsampled process xu is s times higher than that of the original process x, the cutoff frequency of the low-pass filter applied to x for obtaining xf becomes fτs=s/(2τ).The next step of the procedure is to provide a representation for the downsampled processxd, where downsampling is applied to the filtered process xf. To this end, we exploit recent theoretical findings leading to an analytical derivation of the parameters of the SS model which describes the downsampled process [20, 21, 30]. According to these findings, the downsampled process xd(n)=xf(nτ) can be represented in SS form involving a vector state process Y as(11)Yn+1=Ad·Yn+Wnxdn=Cd·Yn+vn,where the parameters of the SS model (11) are the matrix Ad, the vector Cd, the variance of the scalar noise process v(n), σv2, the covariance matrix of the vector state noise process W(n), ΣW, and the cross-covariance between v(n) and W(n), ΣWv. These parameters can be computed at scale τ, in terms of parameters previously obtained, as follows [20, 30]:(12)Ad=Afτ;Cd=Cf;σv2=σef2ΣWv=Afτ-1Kfσef2ΣWτ=AfΣWτ-1AfT+σef2KfKfT,τ≥2ΣWτ=σef2KfKfT,τ=1.Then, the SS model (11) can be transformed in a form similar to that of (9) which evidences the innovations:(13)Zdn+1=Ad·Zdn+Kdedn,xdn=Cd·Zdn+edn,by solving a so-called discrete algebraic Riccati equation [20, 30](14)P=AdPAdT+ΣW-AdPCd+ΣWvCdPCdT+σv2-1CdPAdT+ΣWvTapplied to the parameters of the SS model (11), to complete the derivation of the parameters of (13) as(15)σed2=CdPCdT+σv2Kd=AdPCdT+ΣWvσv2.Finally, the variance of the downsampled process can be computed analytically solving a discrete-time Lyapunov equation [20] as(16)Ω=AdΩAdT+σv2KdKdTσxd2=CdΩCdT+σv2.The derivations above allow computing analytically all the parameters of the state-space process (see (13)) which describes the filtered (see (9)) and downsampled (see (13)) versions of the upsampled process (see (7)), which in turn can be analytically described starting from the original AR process x given by (3). Among these parameters, the relevant ones for our purpose are the variance of the downsampled process xd and that of the relevant innovations ed, which are used to compute the complexity of xd in analogy to (6), thus obtaining our MSE measure:(17)Cxτs=Cxd=12ln⁡2πeσed2σxd2.Note that MSE measure defined in (17), which we denote as linear MSE (LMSE), is relevant to a time scale given by the rational number τs=τ/s, obtained setting integer values for the scale factor τ and the upsampling factor s. As this corresponds to applying to the original process a low-pass filter with cutoff frequency fτs=1/(2τs), such cutoff frequency can be arranged according to the value of τs in order to provide fine-tuning of the filtering process. Moreover we stress that since we filter the original process with a FIR filter that prevents aliasing and we compute at each time the complexity of the normalized process, our state-space MSE measure follows the philosophy of the RMSE method rather than that of the original MSE.


### 2.4. Practical Implementation and Applicability of LMSE and RMSE

The proposed approach for multiscale complexity analysis is implemented in the LMSE MATLAB® toolbox, which includes the algorithms for computing LMSE and RMSE for the simulated processes and exemplificative realizations of the cardiovascular data studied in this paper. The toolbox is uploaded as supplementary material to this article and is freely available for download fromhttp://www.lucafaes.net/LMSE.html. The codes allow also nonexpert users to implement multiscale complexity analysis without the need to go deep into the mathematical details presented in the previous subsections. Complexity is assessed at a given time scale, after removing the faster temporal scales through low-pass filtering, as the unpredictability of future values of the time series given past observations. Unpredictability is quantified either using a model-free estimation of the conditional entropy (RMSE) or using a linear model (LMSE). As we will see in the next two sections, the analytical expressions underlying the computation of LMSE make this measure much more reliable than RMSE in the presence of short time series and long temporal scales. On the other hand, one should bear in mind that the context of application for LMSE is that of time series which are adequately represented by linear stationary stochastic processes. Linear stochastic processes are nondeterministic dynamic processes for which the relation between samples can be described by linear equations. The range of applicability of this simple generative model for time series data is very broad, including time series collected from many physical, biological, and social systems. Nevertheless, if the observed system is supposed to generate time series with strong nonlinear dynamics or significant departures from the null hypothesis of Gaussianity of the distribution are proven, the LMSE measure may capture only partly the complexity of these time series; in such a case the RMSE measure is a more appropriate choice. Furthermore, the hypothesis of stationarity implies that the statistical properties of the observed process (e.g., mean and variance) do not vary across time. For time series with evident nonstationary behaviors, the potential applicability of the method here proposed to short time series allows its easy adaptation to time-varying formulations.


## 2.1. Multiscale Entropy

Multiscale entropy (MSE), originally proposed by Costa et al. [5], is a measure assessing the complexity of a process across multiple temporal scales. Its computation relies first on rescaling the observed process (i.e., focusing on a specific range of temporal scales) and then on assessing the dynamical complexity of the rescaled process through the computation of its rate of entropy production.Specifically, let us consider a discrete-time, stationary stochastic processx with zero mean and variance σx2. Let us further define x(n) as the variable obtained sampling the process at the discrete time n and set an integer scale factor τ. The rescaling step is aimed at eliminating the fast temporal scales and is performed applying the following transformation to the original process x, which averages τ consecutive samples to yield the rescaled process xd:(1)xdn=1τ∑i=0τ-1xnτ-i; the rescaled process is denoted here as xd because it is a downsampled version of the original process x; in fact, in (1) one sample of the process xd is obtained by taking one sample every τ values of the averaged version of x, where averaging is performed over τ consecutive samples. The second part of the procedure is implemented by calculating the conditional entropy of the present state of the rescaled process, xd(n), given its past states, xd–(n)=[xdn-1,xdn-2,…]. In MSE, this is accomplished approximating the past history of the rescaled process with the finite-length pattern xdm(n)=[xd(n-1),…,xd(n-m)] and then on estimating the conditional entropy by means of the Sample Entropy (SampEn) metric [13]. SampEn is an estimate of the probability that patterns which are detected as similar in the m-dimensional space remain similar also in the (m+1)-dimensional space when they are incremented with their future values, that is, the probability that xdm+1(n1+1) and xdm+1(n2+1) are similar given that xdm(n1) and xdm(n2) are similar; similarity is assessed through a coarse-graining procedure which applies the Heaviside step function with parameter r on the distance between patterns (i.e., r is a threshold such that xdm(n1) and xdm(n2) are similar if their distance in the m-dimensional space is lower than r, where the distance is assessed using the maximum norm).The MSE measure resulting from the above procedure, which we denote asCx(τ), quantifies the dynamical complexity of the original process x observed at scale τ. The free parameters of the MSE estimator are the length m of the patterns used to approximate the past of the process and the threshold r that sets the similarity between patterns. In the application of MSE, the pattern length is limited to a few samples (typically, m=2), while the threshold distance is a fraction of the standard deviation of the original process (typically, r=0.2·σx).


## 2.2. Refined Multiscale Entropy

The original MSE formulation suffers from two main limitations: the suboptimal procedure for elimination of the fast temporal scales implemented by (1), which tends to introduce spurious oscillations in the rescaled time series, and the fact that the threshold parameter of the coarse graining is kept at a constant value for all time scales, which brings about an artificial reduction of the MSE values with increasing scales. Valencia et al. [11] proposed a refined MSE (RMSE) measure that is aimed at overcoming these limitations. The solution to the first limitation is based on the rationale that the rescaling procedure actually consists of two steps: a filtering step which eliminates the fast temporal scales from the original process and a downsampling step that eliminates the redundancy introduced by the first step. Accordingly, the filtering step was designed to limit as much as possible the aliasing that can occur with the following downsampling step.In fact, the change of scale of a processx is performed first applying a low-pass filter to obtain the filtered process xf(n) and then reducing the sampling rate of the filtered process to obtain the rescaled (downsampled) processxd(n). The two steps yield, respectively, the processes(2a)xfn=∑i=0qbixn-i-∑j=1rcjxfn-j,(2b)xdn=xfnτ,where b(i) and c(j) are the filter coefficients and q and r are the filter orders. As it averages τ consecutive samples, the original MSE formulation [5] implicitly applies a finite impulse response (FIR) low-pass filter [24] of orders q=τ–1 and r=0 (the filter is nonrecursive) and with coefficients b(i)=1/τ for each i=0,1,…,q=τ–1; the cutoff frequency of the filter is constrained to the value fτ=1/(2τ). To improve elimination of the fast temporal scales and satisfy the Shannon theorem in the subsequent downsampling step, the RMSE approach implements an infinite impulse response (IIR) avoiding ripples in the stop band; specifically, a Butterworth filter [24] of order q (with r=q-1) is implemented in which the coefficients b(i) and c(j) are set to realize a low-pass implementation with cutoff frequency fτ [11]. This choice brings also the advantage that since the cutoff frequency fτ can take any real value, in RMSE the scale factor τ=1/(2fτ) is not constrained to integer numbers as in MSE; this allows filtering out the oscillatory components of the original process with a better resolution before computing the complexity of the rescaled process.Besides the type of filter, another crucial difference exists between the original and refined formulations of MSE. Whereas in MSE the parameterr is constant for all scale factors, as it is set at a fraction of the standard deviation of the original process (e.g., r=0.2·σx), in RMSE this parameter is set at a fraction of the standard deviation of the rescaled process (e.g., r=0.2·σxd). This choice is meant to reduce the dependence of the estimated conditional entropy on the decrease of variance due to filtering. As a consequence of this refinement, RMSE does not exhibit a tendency to decline with the time scale as a consequence of the inherent reduction of the variance of the filtered process: for example, for a white noise process MSE decreases with the time scale while RMSE is constant.


## 2.3. Linear Multiscale Entropy

In this section we present a method to assess the multiscale complexity of linear Gaussian stochastic processes. The method is based on the fact that if the variables obtained sampling the considered processx have a joint Gaussian distribution, all the variability that defines the entropy rate of the process is fully captured by a linear autoregressive (AR) model [25]. Accordingly, let us consider the linear AR representation of the process x, defined as(3)xn=∑k=1pakxn-k+en,where p is the order of the process, a(k), k=1,…,p, are linear regression coefficients describing the interactions between the process variables as a function of the lag k, and e(n) is an innovation process with variance σe2. The process x is fully described by the AR model (3) when the innovation processe is formed by uncorrelated Gaussian variables.Given the AR representation, the variance of the process and of the innovations can be used to derive an information-theoretic description of the statistical structure of the process [26]. Specifically, the entropy of the process is related to its variance by [27](4)Hx=12ln⁡2πeσx2, and the entropy rate of the process, that is, the conditional entropy of the present x(n) given the past x–(n)=[xn-1,xn-2,…], is related to the variance of the innovations by [28](5)Hx∣x-=12ln⁡2πeσe2. Then, (4) and (5) can be combined to provide a version of the conditional entropy which quantifies the dynamical complexity of the normalized process as(6)Cx=Hx∣x--Hx+12ln⁡2πe=12ln⁡2πeσe2σx2; note that applying (6) corresponds to computing the conditional entropy of the original process after normalizing it to unit variance.Now we turn to show how to compute analytically the variance of the AR innovations after rescaling the original process, in a way such that this variance can be used as in (6) to assess multiscale complexity. First, we perform a virtual upsampling of the original process x, by setting an integer upsampling factor s and by defining the process(7)xun=∑k=1psaukxun-k+eun,where ps=p·s is the order of the upsampled process and the coefficients au are set to be a zero-padded version of the original AR coefficients; that is, we set au(ls)=a(l) for each l=1,…,p and au(k)=0 for each k=1,…,ps, k≠ls. Note that the innovations for the original and upsampled processes are formally the same, that is, σeu2=σe2. Then, we filter the upsampled process using a linear FIR filter, that is, a filter implemented as in (2a) with c(j)=0 for each j; with some algebraic manipulation we find that, substituting (7) in (2a) applied with r=0 and with xu in place of x, the filtered representation of the upsampled process takes the form(8)xfn=∑k=1psaukxfn-k+∑i=0qbieun-i.Equation (8) shows that the filtering step introduces a moving average (MA) component in the original AR process, transforming it into an ARMA process [21]. Then, we exploit the connection between ARMA processes and state-space (SS) processes [29] to evidence that the ARMA process (8) can be expressed in an SS form as(9)Zfn+1=Af·Zfn+Kfefn,xfn=Cf·Zfn+efn,where Zf(n)=[xf(n-1)⋯xf(n-ps)u(n-1)⋯u(n-q)]T is a (ps+q)-dimensional state process, ef(n) is the scalar SS innovation process defined as ef(n)=b0eu(n), and the vectors Cf and Kf and the matrix Af are defined as(10)Af=au1⋯aups-1aupsb1⋯bq-1bq1⋯000⋯00⋮⋱⋮⋮⋮⋱⋮⋮0⋯100⋯000⋯000⋯000⋯001⋯000⋯0⋮⋮⋱⋮⋮0⋯000⋯10.Kf=10⋮0b00⋮0Cf=au1⋯aups-1aupsb1⋯bq-1bq.Such a connection between ARMA and SS models allows finding, through (10) and the relation σef2=b(0)2σeu2, the parameters of the SS process descriptive of the process xf. Such a process is a filtered version of the upsampled process xu, obtained passing xu through a FIR low-pass filter with cutoff frequency equal to fτ=1/(2τ); therefore, since the sampling frequency of the upsampled process xu is s times higher than that of the original process x, the cutoff frequency of the low-pass filter applied to x for obtaining xf becomes fτs=s/(2τ).The next step of the procedure is to provide a representation for the downsampled processxd, where downsampling is applied to the filtered process xf. To this end, we exploit recent theoretical findings leading to an analytical derivation of the parameters of the SS model which describes the downsampled process [20, 21, 30]. According to these findings, the downsampled process xd(n)=xf(nτ) can be represented in SS form involving a vector state process Y as(11)Yn+1=Ad·Yn+Wnxdn=Cd·Yn+vn,where the parameters of the SS model (11) are the matrix Ad, the vector Cd, the variance of the scalar noise process v(n), σv2, the covariance matrix of the vector state noise process W(n), ΣW, and the cross-covariance between v(n) and W(n), ΣWv. These parameters can be computed at scale τ, in terms of parameters previously obtained, as follows [20, 30]:(12)Ad=Afτ;Cd=Cf;σv2=σef2ΣWv=Afτ-1Kfσef2ΣWτ=AfΣWτ-1AfT+σef2KfKfT,τ≥2ΣWτ=σef2KfKfT,τ=1.Then, the SS model (11) can be transformed in a form similar to that of (9) which evidences the innovations:(13)Zdn+1=Ad·Zdn+Kdedn,xdn=Cd·Zdn+edn,by solving a so-called discrete algebraic Riccati equation [20, 30](14)P=AdPAdT+ΣW-AdPCd+ΣWvCdPCdT+σv2-1CdPAdT+ΣWvTapplied to the parameters of the SS model (11), to complete the derivation of the parameters of (13) as(15)σed2=CdPCdT+σv2Kd=AdPCdT+ΣWvσv2.Finally, the variance of the downsampled process can be computed analytically solving a discrete-time Lyapunov equation [20] as(16)Ω=AdΩAdT+σv2KdKdTσxd2=CdΩCdT+σv2.The derivations above allow computing analytically all the parameters of the state-space process (see (13)) which describes the filtered (see (9)) and downsampled (see (13)) versions of the upsampled process (see (7)), which in turn can be analytically described starting from the original AR process x given by (3). Among these parameters, the relevant ones for our purpose are the variance of the downsampled process xd and that of the relevant innovations ed, which are used to compute the complexity of xd in analogy to (6), thus obtaining our MSE measure:(17)Cxτs=Cxd=12ln⁡2πeσed2σxd2.Note that MSE measure defined in (17), which we denote as linear MSE (LMSE), is relevant to a time scale given by the rational number τs=τ/s, obtained setting integer values for the scale factor τ and the upsampling factor s. As this corresponds to applying to the original process a low-pass filter with cutoff frequency fτs=1/(2τs), such cutoff frequency can be arranged according to the value of τs in order to provide fine-tuning of the filtering process. Moreover we stress that since we filter the original process with a FIR filter that prevents aliasing and we compute at each time the complexity of the normalized process, our state-space MSE measure follows the philosophy of the RMSE method rather than that of the original MSE.


## 2.4. Practical Implementation and Applicability of LMSE and RMSE

The proposed approach for multiscale complexity analysis is implemented in the LMSE MATLAB® toolbox, which includes the algorithms for computing LMSE and RMSE for the simulated processes and exemplificative realizations of the cardiovascular data studied in this paper. The toolbox is uploaded as supplementary material to this article and is freely available for download fromhttp://www.lucafaes.net/LMSE.html. The codes allow also nonexpert users to implement multiscale complexity analysis without the need to go deep into the mathematical details presented in the previous subsections. Complexity is assessed at a given time scale, after removing the faster temporal scales through low-pass filtering, as the unpredictability of future values of the time series given past observations. Unpredictability is quantified either using a model-free estimation of the conditional entropy (RMSE) or using a linear model (LMSE). As we will see in the next two sections, the analytical expressions underlying the computation of LMSE make this measure much more reliable than RMSE in the presence of short time series and long temporal scales. On the other hand, one should bear in mind that the context of application for LMSE is that of time series which are adequately represented by linear stationary stochastic processes. Linear stochastic processes are nondeterministic dynamic processes for which the relation between samples can be described by linear equations. The range of applicability of this simple generative model for time series data is very broad, including time series collected from many physical, biological, and social systems. Nevertheless, if the observed system is supposed to generate time series with strong nonlinear dynamics or significant departures from the null hypothesis of Gaussianity of the distribution are proven, the LMSE measure may capture only partly the complexity of these time series; in such a case the RMSE measure is a more appropriate choice. Furthermore, the hypothesis of stationarity implies that the statistical properties of the observed process (e.g., mean and variance) do not vary across time. For time series with evident nonstationary behaviors, the potential applicability of the method here proposed to short time series allows its easy adaptation to time-varying formulations.


## 3. Simulation Study

To investigate the theoretical profiles of the dynamical complexity of a stochastic process as a function of the parameters determining its dynamics, as well as to assess the computational reliability of the proposed estimator of the refined MSE, in this section we consider a set of linear processes simulated with varying oscillatory components, for which we compute exact values of MSE and compare LMSE and RMSE estimates obtained from short process realizations.


### 3.1. Simulation Design

Simulations are designed considering AR processes described byl pairs of complex-conjugate poles with assigned modulus ρi and phase φi=2πfi, i=1,…,l (the order of the process is p=2l). Type 1 simulations are designed with one pair of complex-conjugate poles (l=1, p=2) according to two configurations: (a) fixing the pole phase to φ1=π/5 (frequency f1=0.1Hz) and varying the modulus in the range ρ1={0,0.5,0.8,0.9}; (b) fixing the pole modulus to ρ1=0.8 and varying the frequency in the range f1={0.1,0.2,0.3,0.4}. Type 2 simulations are designed with two pairs of complex-conjugate poles (l=2, p=4) according to two configurations both with a fixed low-frequency pole obtained setting ρ1=0.8 and f1=0.1; (c) fixing the phase of the second pole to φ2=π/2 (frequency f2=0.25Hz) and varying the modulus in the range ρ2={0,0.3,0.6,0.8}; (d) fixing the modulus of the second pole to ρ2=0.8 and varying the frequency in the range f2=0.15,0.25,0.3. Type 1 simulations (configurations (a, b)) are set to study how the complexity of a process with a single stochastic oscillatory component varies with the amplitude and the frequency of such component. Type 2 simulations (configurations (c, d)) are set to study how the complexity of a process with two stochastic oscillatory components varies with the mismatch in the amplitude and in the frequency of these two components.Given the configurations described above, first we determine the theoretical values of the MSE using the procedure described in Section2.3. To this end, the exact values of the AR coefficients are determined as the coefficients of the polynomial with roots given by the poles set in the complex plane as zi=ρi·(cos⁡φi±jsin⁡φi); these coefficients, together with the innovation variance σe2=1, are set as parameters of the AR model of (3). From these parameters, MSE values are computed for specific values of the upsampling factor s and the scale factor τ, that is, (s,τ)={(1,1),(8,9),(4,5),(7,10),(3,5),(5,9),(1,2),(8,18),(2,5), (7,20), (3,10),(1,4),(1,5),(3,20),(1,10),(1,20)}, chosen in order to yield an approximately uniform range of values for the cutoff frequency of the rescaling low-pass filter, that is, fτs=s/(2τ)={0.5, 0.444, 0.4, 0.35, 0.3, 0.278, 0.25, 0.222, 0.2, 0.175, 0.15, 0.125, 0.1, 0.075, 0.05, 0.025}; the filter order is set to q=48.After theoretical analysis, practical estimation of MSE was performed choosing two representative cases of the parameter setting for Type 1 simulation(ρ1=0.8, f1=0.1) and Type 2 simulation (ρ1=0.8, f1=0.1; ρ2=0.8, f2=0.2), generating for each case 100 realizations of the simulation by feeding the model of (3) with realizations of 300 white noise samples taken from the Gaussian distribution with zero mean and unit variance and computing LMSE and RMSE estimates. For each realization, LMSE was obtained through the procedure described in Section 2.3, identifying an AR model through the standard least squares method, setting the model order according to the Bayesian Information Criterion [31], and using a FIR low-pass filter of order q=48 [24]; RMSE was obtained through the procedure described in Section 2.2 and using the parameters of [11], that is, filtering the data with a Butterworth low-pass filter of order q=6 [24] and calculating SampEn with embedding and tolerance parameters m=2 and r=0.2·σ^xd (σ^xd is the variance of the downsampled signal).


### 3.2. Theoretical Analysis

Results of the theoretical analysis are reported in Figure1. As a general result, the complexity of the simulated AR processes tends to increase at increasing the time scale τs=τ/s (i.e., at decreasing the cutoff frequency of the low-pass filter applied to the original process, fτs=1/(2τs)=s/(2τ)). This result is expected, because filtering removes the oscillatory components that bring regular dynamics to the signal, in a way such that when all stochastic oscillations are removed the signal is left with no regular dynamics and the maximum complexity level is achieved (i.e., Cx=0.5ln⁡2πe, corresponding to uncorrelated white noise).Figure 1
Theoretical profiles of the multiscale entropy (MSE) computed applying the proposed approach to the true parameters of simulated AR processes. Plots depict the exact values of MSE (Cx(τs)) computed as a function of the cutoff frequency of the low-pass filter used to eliminate the fast temporal scales of the process (fτs) for different values of the modulus ρ1 (a) and of the frequency f1 (b) of the pole set in Type 1 simulations, and for different values of the modulus ρ2 (c) and frequency f2 (d) of the second pole set in Type 2 simulations. The multiscale complexity of AR processes increases with decreasing the amplitude (a) or moving the frequency of a single stochastic oscillation closer to half the Nyquist frequency (b), as well as adding a second stochastic oscillation with increasing amplitude (c) or with increasing frequency mismatch compared to the first one (d).
(a)
Theoretical MSE
(b)
Theoretical MSE
(c)
Theoretical MSE
(d)
Theoretical MSEFigures1(a) and 1(b) report the exact values of MSE computed for simulated stochastic processes featuring a single oscillatory component. Generating oscillations with fixed frequency and varying amplitude (Figure 1(a)), the MSE of the process increases at decreasing the pole modulus, documenting the higher complexity of stochastic oscillations associated with poles with higher distance from the unit circle in the complex plane [32]. Looking at the multiscale behavior, MSE reaches its maximum value when the regular component, oscillating at f1=0.1 Hz in this example, is completely removed by the filtering procedure; the slope of the rise in complexity decreases with the pole modulus, becoming zero for ρ1=0 when the process is a Gaussian white noise without regular dynamics. The behavior of MSE for oscillations with fixed amplitude (ρ1=0.8) and varying frequency is more complicated (Figure 1(b)). At scale one (fτs=0.5), the MSE is the same if the frequency of the oscillatory dynamics, f1, has the same distance from half the Nyquist frequency of 0.25 Hz and decreases with such a distance; the same symmetric behavior, with maximum complexity at half the Nyquist frequency, was observed in [33]. Then, the multiscale behavior of the complexity measure is related to the frequency of the stochastic oscillation in a way such that faster oscillations are removed more easily than slower oscillations, and thus MSE reaches its higher values for higher values of fτs (see the trends of f1=0.4 versus f1=0.1 and f1=0.3 versus f1=0.2, in Figure 1(b)).Figures1(c) and 1(d) show the exact values of MSE computed for simulated stochastic processes featuring two oscillatory components. Simulating the rise of a high-frequency component in the presence of a stable low-frequency stochastic oscillation (Figure 1(c), where ρ1=0.8, f1=0.1, and ρ2 increases with fixed f2=0.25) determines an increase of MSE that is revealed across multiple temporal scales. Simulating the progressive separation of stochastic oscillations within a process (Figure 1(d), where ρ1=0.8, f1=0.1, and f2 increases with fixed ρ2=0.8) determines again an increase of MSE across multiple scales. These results suggest that the simultaneous presence of multiple oscillatory mechanisms tends to produce more complex dynamics than in the case of single mechanisms, with a complexity degree that increases with the strength of the stochastic oscillations and with the mismatch of their frequency. These findings are supported by the results of recent studies [32, 33] and are observed over a range of time scales that comprises the characteristic periods of all oscillations.


### 3.3. Estimation Performance

Figure2 reports the results of the practical estimation of complexity over short realizations of the simulations, performed using the refined approach and the linear approach to MSE computation. For both the parameter settings chosen to simulate an individual stochastic oscillation (Figures 2(a) and 2(b); ρ1=0.8, f1=0.1) and a pair of stochastic oscillations (Figures 2(c) and 2(d); ρ1=0.8, f1=0.1; ρ2=0.8, f2=0.2), it is evident that the LMSE estimator proposed in this study outperforms the RMSE estimator. In fact, LMSE estimates are substantially unbiased, since the median estimated profile of multiscale complexity overlaps with the true profile, and exhibit a low variability around their median value, particularly at high time scales where they converge to the expected value 0.5ln⁡2πe (Figures 2(a) and 2(c)). On the contrary, RMSE estimates are strongly biased at all time scales and also exhibit a substantial variance that increases dramatically with the time scale (Figures 2(b) and 2(d)); note that RMSE could not be computed at the higher time scales (fτ=0.05, fτ=0.025) due to the limited number of data points. The poor performance of RMSE confirms the known difficulty of yielding accurate model-free complexity estimates using the SampEn estimator applied to short time series [33, 34].Figure 2
Estimation of multiscale entropy (MSE) over finite-length realizations of simulated AR processes. Plots depict the exact values (red lines) and the distributions (median and 10th–90th percentiles) of the MSE estimates (Cx(τs)) computed as a function of the cutoff frequency of the low-pass filter used to eliminate the fast temporal scales of the process (fτs) for representative parameter settings of Type 1 simulation (a, b) and of Type 2 simulation (c, d). Estimates are obtained using the linear MSE (LMSE) method proposed in this study (a, c) and using the refined MSE (RMSE) method proposed in [11] (b, d). In simulated AR processes, RMSE estimates exhibit high bias and a variance increasing with the time scale, while LMSE estimates display high computational reliability at all time scales.
(a)
LMSE
(b)
RMSE
(c)
LMSE
(d)
RMSE


## 3.1. Simulation Design

Simulations are designed considering AR processes described byl pairs of complex-conjugate poles with assigned modulus ρi and phase φi=2πfi, i=1,…,l (the order of the process is p=2l). Type 1 simulations are designed with one pair of complex-conjugate poles (l=1, p=2) according to two configurations: (a) fixing the pole phase to φ1=π/5 (frequency f1=0.1Hz) and varying the modulus in the range ρ1={0,0.5,0.8,0.9}; (b) fixing the pole modulus to ρ1=0.8 and varying the frequency in the range f1={0.1,0.2,0.3,0.4}. Type 2 simulations are designed with two pairs of complex-conjugate poles (l=2, p=4) according to two configurations both with a fixed low-frequency pole obtained setting ρ1=0.8 and f1=0.1; (c) fixing the phase of the second pole to φ2=π/2 (frequency f2=0.25Hz) and varying the modulus in the range ρ2={0,0.3,0.6,0.8}; (d) fixing the modulus of the second pole to ρ2=0.8 and varying the frequency in the range f2=0.15,0.25,0.3. Type 1 simulations (configurations (a, b)) are set to study how the complexity of a process with a single stochastic oscillatory component varies with the amplitude and the frequency of such component. Type 2 simulations (configurations (c, d)) are set to study how the complexity of a process with two stochastic oscillatory components varies with the mismatch in the amplitude and in the frequency of these two components.Given the configurations described above, first we determine the theoretical values of the MSE using the procedure described in Section2.3. To this end, the exact values of the AR coefficients are determined as the coefficients of the polynomial with roots given by the poles set in the complex plane as zi=ρi·(cos⁡φi±jsin⁡φi); these coefficients, together with the innovation variance σe2=1, are set as parameters of the AR model of (3). From these parameters, MSE values are computed for specific values of the upsampling factor s and the scale factor τ, that is, (s,τ)={(1,1),(8,9),(4,5),(7,10),(3,5),(5,9),(1,2),(8,18),(2,5), (7,20), (3,10),(1,4),(1,5),(3,20),(1,10),(1,20)}, chosen in order to yield an approximately uniform range of values for the cutoff frequency of the rescaling low-pass filter, that is, fτs=s/(2τ)={0.5, 0.444, 0.4, 0.35, 0.3, 0.278, 0.25, 0.222, 0.2, 0.175, 0.15, 0.125, 0.1, 0.075, 0.05, 0.025}; the filter order is set to q=48.After theoretical analysis, practical estimation of MSE was performed choosing two representative cases of the parameter setting for Type 1 simulation(ρ1=0.8, f1=0.1) and Type 2 simulation (ρ1=0.8, f1=0.1; ρ2=0.8, f2=0.2), generating for each case 100 realizations of the simulation by feeding the model of (3) with realizations of 300 white noise samples taken from the Gaussian distribution with zero mean and unit variance and computing LMSE and RMSE estimates. For each realization, LMSE was obtained through the procedure described in Section 2.3, identifying an AR model through the standard least squares method, setting the model order according to the Bayesian Information Criterion [31], and using a FIR low-pass filter of order q=48 [24]; RMSE was obtained through the procedure described in Section 2.2 and using the parameters of [11], that is, filtering the data with a Butterworth low-pass filter of order q=6 [24] and calculating SampEn with embedding and tolerance parameters m=2 and r=0.2·σ^xd (σ^xd is the variance of the downsampled signal).


## 3.2. Theoretical Analysis

Results of the theoretical analysis are reported in Figure1. As a general result, the complexity of the simulated AR processes tends to increase at increasing the time scale τs=τ/s (i.e., at decreasing the cutoff frequency of the low-pass filter applied to the original process, fτs=1/(2τs)=s/(2τ)). This result is expected, because filtering removes the oscillatory components that bring regular dynamics to the signal, in a way such that when all stochastic oscillations are removed the signal is left with no regular dynamics and the maximum complexity level is achieved (i.e., Cx=0.5ln⁡2πe, corresponding to uncorrelated white noise).Figure 1
Theoretical profiles of the multiscale entropy (MSE) computed applying the proposed approach to the true parameters of simulated AR processes. Plots depict the exact values of MSE (Cx(τs)) computed as a function of the cutoff frequency of the low-pass filter used to eliminate the fast temporal scales of the process (fτs) for different values of the modulus ρ1 (a) and of the frequency f1 (b) of the pole set in Type 1 simulations, and for different values of the modulus ρ2 (c) and frequency f2 (d) of the second pole set in Type 2 simulations. The multiscale complexity of AR processes increases with decreasing the amplitude (a) or moving the frequency of a single stochastic oscillation closer to half the Nyquist frequency (b), as well as adding a second stochastic oscillation with increasing amplitude (c) or with increasing frequency mismatch compared to the first one (d).
(a)
Theoretical MSE
(b)
Theoretical MSE
(c)
Theoretical MSE
(d)
Theoretical MSEFigures1(a) and 1(b) report the exact values of MSE computed for simulated stochastic processes featuring a single oscillatory component. Generating oscillations with fixed frequency and varying amplitude (Figure 1(a)), the MSE of the process increases at decreasing the pole modulus, documenting the higher complexity of stochastic oscillations associated with poles with higher distance from the unit circle in the complex plane [32]. Looking at the multiscale behavior, MSE reaches its maximum value when the regular component, oscillating at f1=0.1 Hz in this example, is completely removed by the filtering procedure; the slope of the rise in complexity decreases with the pole modulus, becoming zero for ρ1=0 when the process is a Gaussian white noise without regular dynamics. The behavior of MSE for oscillations with fixed amplitude (ρ1=0.8) and varying frequency is more complicated (Figure 1(b)). At scale one (fτs=0.5), the MSE is the same if the frequency of the oscillatory dynamics, f1, has the same distance from half the Nyquist frequency of 0.25 Hz and decreases with such a distance; the same symmetric behavior, with maximum complexity at half the Nyquist frequency, was observed in [33]. Then, the multiscale behavior of the complexity measure is related to the frequency of the stochastic oscillation in a way such that faster oscillations are removed more easily than slower oscillations, and thus MSE reaches its higher values for higher values of fτs (see the trends of f1=0.4 versus f1=0.1 and f1=0.3 versus f1=0.2, in Figure 1(b)).Figures1(c) and 1(d) show the exact values of MSE computed for simulated stochastic processes featuring two oscillatory components. Simulating the rise of a high-frequency component in the presence of a stable low-frequency stochastic oscillation (Figure 1(c), where ρ1=0.8, f1=0.1, and ρ2 increases with fixed f2=0.25) determines an increase of MSE that is revealed across multiple temporal scales. Simulating the progressive separation of stochastic oscillations within a process (Figure 1(d), where ρ1=0.8, f1=0.1, and f2 increases with fixed ρ2=0.8) determines again an increase of MSE across multiple scales. These results suggest that the simultaneous presence of multiple oscillatory mechanisms tends to produce more complex dynamics than in the case of single mechanisms, with a complexity degree that increases with the strength of the stochastic oscillations and with the mismatch of their frequency. These findings are supported by the results of recent studies [32, 33] and are observed over a range of time scales that comprises the characteristic periods of all oscillations.


## 3.3. Estimation Performance

Figure2 reports the results of the practical estimation of complexity over short realizations of the simulations, performed using the refined approach and the linear approach to MSE computation. For both the parameter settings chosen to simulate an individual stochastic oscillation (Figures 2(a) and 2(b); ρ1=0.8, f1=0.1) and a pair of stochastic oscillations (Figures 2(c) and 2(d); ρ1=0.8, f1=0.1; ρ2=0.8, f2=0.2), it is evident that the LMSE estimator proposed in this study outperforms the RMSE estimator. In fact, LMSE estimates are substantially unbiased, since the median estimated profile of multiscale complexity overlaps with the true profile, and exhibit a low variability around their median value, particularly at high time scales where they converge to the expected value 0.5ln⁡2πe (Figures 2(a) and 2(c)). On the contrary, RMSE estimates are strongly biased at all time scales and also exhibit a substantial variance that increases dramatically with the time scale (Figures 2(b) and 2(d)); note that RMSE could not be computed at the higher time scales (fτ=0.05, fτ=0.025) due to the limited number of data points. The poor performance of RMSE confirms the known difficulty of yielding accurate model-free complexity estimates using the SampEn estimator applied to short time series [33, 34].Figure 2
Estimation of multiscale entropy (MSE) over finite-length realizations of simulated AR processes. Plots depict the exact values (red lines) and the distributions (median and 10th–90th percentiles) of the MSE estimates (Cx(τs)) computed as a function of the cutoff frequency of the low-pass filter used to eliminate the fast temporal scales of the process (fτs) for representative parameter settings of Type 1 simulation (a, b) and of Type 2 simulation (c, d). Estimates are obtained using the linear MSE (LMSE) method proposed in this study (a, c) and using the refined MSE (RMSE) method proposed in [11] (b, d). In simulated AR processes, RMSE estimates exhibit high bias and a variance increasing with the time scale, while LMSE estimates display high computational reliability at all time scales.
(a)
LMSE
(b)
RMSE
(c)
LMSE
(d)
RMSE


## 4. Application to Short-Term Cardiovascular Variability Series

To illustrate the application of the proposed approach for the computation of MSE over short biomedical time series, this section reports the analysis of cardiovascular and respiratory variability series. Specifically, we compare the abilities of LMSE and RMSE in detecting the multiscale complexity of heart period (HP), systolic arterial pressure (SAP), and respiration (RESP) time series measured from a large group of healthy subjects in a resting state condition as well as during two types of physiological stress, that is, postural stress and mental stress [35].


### 4.1. Experimental Protocol and Data Analysis

The analyzed time series belong to a database collected to assess the dynamics of HP, SAP, and RESP during two types of physiological stress commonly studied in cardiovascular variability, that is, postural stress induced by head-up tilt (HUT) and mental stress induced by mental arithmetics (MA); we refer to [35, 36] for a detailed description of the population and the experimental setup. Briefly, a group of 61 young healthy subjects (37 females, 17.5±2.4 years old) were monitored in a relaxed state in the resting supine position (SU), in the upright position during HUT, and in the supine position during MA, measuring the surface electrocardiogram (ECG), the finger arterial blood pressure collected noninvasively by the photoplethysmographic method, and the respiration signal measured through respiratory inductive plethysmography. From these signals, the beat-to-beat time series of HP, SAP, and RESP were measured, respectively, as the sequences of the temporal distances between consecutive R peaks of the ECG, the maximum values of the arterial pressure waveform measured inside the consecutively detected heart periods, and the values of the respiratory signal sampled at the onset of the consecutively detected heart periods.The analysis was performed on segments of 300 consecutive points, free of artifacts and satisfying stationarity requirements, extracted from the three time series for each subject and condition. The preprocessing steps consisted in removing the linear trend from each sequence and in reducing the series to zero mean. Then, for each individual time series, a linear AR model was identified using the standard least squares method and using the Bayesian Information Criterion to set the model order within the range{1,…,12} [37]. From the estimated AR parameters, LMSE was computed by means of the procedure described in Section 2.3, using a low-pass FIR filter of order q=48. The computation of RMSE was performed as described in Section 2.2, using a sixth-order Butterworth low-pass filter before resampling and computing SampEn with parameters m=2 and r equal to 20% of the variance of the resampled signal. As for simulations, in the computation of both LMSE and RMSE, the parameters s and τ determining the time scale were set to obtain the following values for the cutoff frequency of the resampling filter: fτs=1/(2τs)=s/(2τ)={0.5, 0.444, 0.4, 0.35, 0.3, 0.278, 0.25, 0.222, 0.2, 0.175, 0.15, 0.125, 0.1, 0.075, 0.05, 0.025}.Statistically significant differences among the MSE profiles obtained in the three conditions (i.e., SU, HUT, and MA) were first assessed by means of the multivariate ANOVA. Then, if the null hypothesis that the means of MSE computed across time scales for each condition are the same multivariate vector was rejected, the univariate ANOVA was applied to the three distributions of MSE obtained during SU, HUT, and MA at any assigned time scale. Furthermore, if at a given time scale the null hypothesis that the means of MSE computed in the three conditions are the same number was rejected, a post hoc pairwise test (i.e., the Studentt-test for paired data) was performed to assess the statistical significance of the differences between rest and stress conditions (i.e., HUT versus. REST, or MA versus. REST). A p value < 0.05 was always assumed as statistically significant; both in the univariate ANOVA and in the pairwise tests, a Bonferroni-Holm correction for multiple comparisons was employed.


### 4.2. Results and Discussion

The results of multiscale complexity analysis of HP, SAP, and RESP are depicted, respectively, in Figures3, 4, and 5; results are presented as median and interquartile range of the distributions of LMSE and RMSE computed for the 61 subjects during the three considered experimental conditions (SU, HUT, and MA). As a general result, we find that LMSE estimates exhibit lower intersubject variability than RMSE estimates, especially when increasing the scale factor (i.e., when reducing the cutoff frequency of the low-pass filter). This result holds for all time series and experimental conditions, likely reflecting the problems associated with the computation of RMSE for short time series and at long time scales. On the contrary, the variability of LMSE estimates decreases when increasing the scale factor, allowing eliciting statistically significant differences that in RMSE are masked by the larger intersubject variability.Figure 3
Estimation of multiscale entropy (MSE) for the time series of the heart period. Plots depict the distributions (median and 25th–75th percentiles) of the MSE estimates(Cx(τs)) computed as a function of the cutoff frequency of the low-pass filter used to eliminate the fast temporal scales of the process (fτs) in the resting supine position (SU, diamonds), during postural stress induced by head-up tilt (HUT, squares), and during mental stress induced by mental arithmetics (MA, circles). While MA does not induce changes, HUT evokes a significant reduction of the complexity of heart period variability, which is observed using LMSE across a wide range of time scales including both low- and high-frequency oscillations. #Statistically significant difference, HUT versus SU or MA versus SU.
(a)
Heart period, LMSE
(b)
Heart period, RMSEFigure 4
Estimation of multiscale entropy (MSE) for the time series of the systolic arterial pressure. Plots depict the distributions (median and 25th–75th percentiles) of the MSE estimates(Cx(τs)) computed as a function of the cutoff frequency of the low-pass filter used to eliminate the fast temporal scales of the process (fτs) in the resting supine position (SU, diamonds), during postural stress induced by head-up tilt (HUT, squares), and during mental stress induced by mental arithmetics (MA, circles). The multiscale complexity of systolic pressure dynamics increases consistently across multiple time scales during MA and increases for time scales associated with low-frequency oscillations during HUT. #Statistically significant difference, HUT versus SU or MA versus SU.
(a)
Systolic pressure, LMSE
(b)
Systolic pressure, RMSEFigure 5
Estimation of multiscale entropy (MSE) for the time series of respiration. Plots depict the distributions (median and 25th–75th percentiles) of the MSE estimates(Cx(τs)) computed as a function of the cutoff frequency of the low-pass filter used to eliminate the fast temporal scales of the process (fτs) in the resting supine position (SU, diamonds), during postural stress induced by head-up tilt (HUT, squares), and during mental stress induced by mental arithmetics (MA, circles). The multiscale complexity of respiratory dynamics increases consistently across multiple time scales during MA and decreases for time scales associated with high-frequency oscillations during HUT. #Statistically significant difference, HUT versus SU or MA versus SU.
(a)
Respiration, LMSE
(b)
Respiration, RMSEThe analysis of LMSE and RMSE computed for the HP time series, reported in Figure3, documents that the complexity of HP is significantly lower during HUT than during REST, while no significant differences are observed between REST and MA. The lower complexity during HUT is detected up to fτs=0.225 Hz using RMSE (Figure 3(b)) and up to fτs=0.1 Hz using LMSE (Figure 3(a)). At scale 1 (fτs=0.5Hz), these results confirm a number of previous investigations documenting the decreased complexity of heart rate variability during postural stress and its unaltered complexity during mental stress [36, 38–40]. Here we extend these findings showing that during MA the complexity of HP is left unchanged at any time scale and during HUT it is kept at lower levels until both the high-frequency (HF, >0.15 Hz) and the low-frequency (LF, ~0.1 Hz) oscillations are removed by filtering. This latter result, which is documented only using the LSME estimator, is in agreement with a recent study showing that the reduction in complexity induced by HUT is a result of the presence of more regular HP oscillations in the LF band rather than in the HF band and is thus associated more to an enriched sympathetic modulation than to a vagal withdrawal [32].The MSE analysis performed for the SAP time series, reported in Figure4, evidences a clear increase in the MSE computed during MA. This result, which is again more evident using LMSE than RMSE, has been previously observed at scale 1 [36] and is possibly related to complex patterns of autonomic activation following cognitive load [41, 42]; our results, showing the persistence of higher MSE values up to long time scales (fτs=0.15Hz), point to an involvement of LF oscillations in the generation of higher complexity during mental stress. The multiscale pattern of complexity induced by postural stress is more complicated, indicating lower MSE during HUT than during REST for short time scales (fτs=0.5 Hz), no significant alterations for intermediate time scales (up to fτs=0.275 Hz), and higher MSE during HUT than during REST for longer time scales (fτs between 0.275 Hz and 0.15 Hz). The emergence of a significantly higher complexity of SAP during HUT when HF fluctuations are filtered out is a novel finding, suggesting that the postural stress decreases the regularity of LF arterial pressure oscillations, while HF oscillations display similar or even increased regularity.The multiscale complexity analysis of respiration variability, reported in Figure5, documents that MSE is significantly decreased during HUT and significantly increased during MA. These variations are consistently observed across multiple scales, going up to fτs=0.3 Hz for the lower LMSE and RMSE during HUT and up to fτs=0.15 Hz for the higher LMSE during MA. The higher regularity of the respiratory dynamics during postural stress, which seems to be confined to the HF band where respiration usually occurs, may be related to an increased tidal volume during this condition. On the other hand, the lower regularity of respiration during mental stress might be explained by the appearance of long pauses or sighs in the respiration pattern while performing mental arithmetics [43], making it more erratic and thus complex. The existence of more erratic patterns, which likely span a broad band of the frequency spectrum, may also explain the fact that we observe higher complexity of respiration up to long time scales, going beyond the common frequency bands at which respiratory activity is observed.


## 4.1. Experimental Protocol and Data Analysis

The analyzed time series belong to a database collected to assess the dynamics of HP, SAP, and RESP during two types of physiological stress commonly studied in cardiovascular variability, that is, postural stress induced by head-up tilt (HUT) and mental stress induced by mental arithmetics (MA); we refer to [35, 36] for a detailed description of the population and the experimental setup. Briefly, a group of 61 young healthy subjects (37 females, 17.5±2.4 years old) were monitored in a relaxed state in the resting supine position (SU), in the upright position during HUT, and in the supine position during MA, measuring the surface electrocardiogram (ECG), the finger arterial blood pressure collected noninvasively by the photoplethysmographic method, and the respiration signal measured through respiratory inductive plethysmography. From these signals, the beat-to-beat time series of HP, SAP, and RESP were measured, respectively, as the sequences of the temporal distances between consecutive R peaks of the ECG, the maximum values of the arterial pressure waveform measured inside the consecutively detected heart periods, and the values of the respiratory signal sampled at the onset of the consecutively detected heart periods.The analysis was performed on segments of 300 consecutive points, free of artifacts and satisfying stationarity requirements, extracted from the three time series for each subject and condition. The preprocessing steps consisted in removing the linear trend from each sequence and in reducing the series to zero mean. Then, for each individual time series, a linear AR model was identified using the standard least squares method and using the Bayesian Information Criterion to set the model order within the range{1,…,12} [37]. From the estimated AR parameters, LMSE was computed by means of the procedure described in Section 2.3, using a low-pass FIR filter of order q=48. The computation of RMSE was performed as described in Section 2.2, using a sixth-order Butterworth low-pass filter before resampling and computing SampEn with parameters m=2 and r equal to 20% of the variance of the resampled signal. As for simulations, in the computation of both LMSE and RMSE, the parameters s and τ determining the time scale were set to obtain the following values for the cutoff frequency of the resampling filter: fτs=1/(2τs)=s/(2τ)={0.5, 0.444, 0.4, 0.35, 0.3, 0.278, 0.25, 0.222, 0.2, 0.175, 0.15, 0.125, 0.1, 0.075, 0.05, 0.025}.Statistically significant differences among the MSE profiles obtained in the three conditions (i.e., SU, HUT, and MA) were first assessed by means of the multivariate ANOVA. Then, if the null hypothesis that the means of MSE computed across time scales for each condition are the same multivariate vector was rejected, the univariate ANOVA was applied to the three distributions of MSE obtained during SU, HUT, and MA at any assigned time scale. Furthermore, if at a given time scale the null hypothesis that the means of MSE computed in the three conditions are the same number was rejected, a post hoc pairwise test (i.e., the Studentt-test for paired data) was performed to assess the statistical significance of the differences between rest and stress conditions (i.e., HUT versus. REST, or MA versus. REST). A p value < 0.05 was always assumed as statistically significant; both in the univariate ANOVA and in the pairwise tests, a Bonferroni-Holm correction for multiple comparisons was employed.


## 4.2. Results and Discussion

The results of multiscale complexity analysis of HP, SAP, and RESP are depicted, respectively, in Figures3, 4, and 5; results are presented as median and interquartile range of the distributions of LMSE and RMSE computed for the 61 subjects during the three considered experimental conditions (SU, HUT, and MA). As a general result, we find that LMSE estimates exhibit lower intersubject variability than RMSE estimates, especially when increasing the scale factor (i.e., when reducing the cutoff frequency of the low-pass filter). This result holds for all time series and experimental conditions, likely reflecting the problems associated with the computation of RMSE for short time series and at long time scales. On the contrary, the variability of LMSE estimates decreases when increasing the scale factor, allowing eliciting statistically significant differences that in RMSE are masked by the larger intersubject variability.Figure 3
Estimation of multiscale entropy (MSE) for the time series of the heart period. Plots depict the distributions (median and 25th–75th percentiles) of the MSE estimates(Cx(τs)) computed as a function of the cutoff frequency of the low-pass filter used to eliminate the fast temporal scales of the process (fτs) in the resting supine position (SU, diamonds), during postural stress induced by head-up tilt (HUT, squares), and during mental stress induced by mental arithmetics (MA, circles). While MA does not induce changes, HUT evokes a significant reduction of the complexity of heart period variability, which is observed using LMSE across a wide range of time scales including both low- and high-frequency oscillations. #Statistically significant difference, HUT versus SU or MA versus SU.
(a)
Heart period, LMSE
(b)
Heart period, RMSEFigure 4
Estimation of multiscale entropy (MSE) for the time series of the systolic arterial pressure. Plots depict the distributions (median and 25th–75th percentiles) of the MSE estimates(Cx(τs)) computed as a function of the cutoff frequency of the low-pass filter used to eliminate the fast temporal scales of the process (fτs) in the resting supine position (SU, diamonds), during postural stress induced by head-up tilt (HUT, squares), and during mental stress induced by mental arithmetics (MA, circles). The multiscale complexity of systolic pressure dynamics increases consistently across multiple time scales during MA and increases for time scales associated with low-frequency oscillations during HUT. #Statistically significant difference, HUT versus SU or MA versus SU.
(a)
Systolic pressure, LMSE
(b)
Systolic pressure, RMSEFigure 5
Estimation of multiscale entropy (MSE) for the time series of respiration. Plots depict the distributions (median and 25th–75th percentiles) of the MSE estimates(Cx(τs)) computed as a function of the cutoff frequency of the low-pass filter used to eliminate the fast temporal scales of the process (fτs) in the resting supine position (SU, diamonds), during postural stress induced by head-up tilt (HUT, squares), and during mental stress induced by mental arithmetics (MA, circles). The multiscale complexity of respiratory dynamics increases consistently across multiple time scales during MA and decreases for time scales associated with high-frequency oscillations during HUT. #Statistically significant difference, HUT versus SU or MA versus SU.
(a)
Respiration, LMSE
(b)
Respiration, RMSEThe analysis of LMSE and RMSE computed for the HP time series, reported in Figure3, documents that the complexity of HP is significantly lower during HUT than during REST, while no significant differences are observed between REST and MA. The lower complexity during HUT is detected up to fτs=0.225 Hz using RMSE (Figure 3(b)) and up to fτs=0.1 Hz using LMSE (Figure 3(a)). At scale 1 (fτs=0.5Hz), these results confirm a number of previous investigations documenting the decreased complexity of heart rate variability during postural stress and its unaltered complexity during mental stress [36, 38–40]. Here we extend these findings showing that during MA the complexity of HP is left unchanged at any time scale and during HUT it is kept at lower levels until both the high-frequency (HF, >0.15 Hz) and the low-frequency (LF, ~0.1 Hz) oscillations are removed by filtering. This latter result, which is documented only using the LSME estimator, is in agreement with a recent study showing that the reduction in complexity induced by HUT is a result of the presence of more regular HP oscillations in the LF band rather than in the HF band and is thus associated more to an enriched sympathetic modulation than to a vagal withdrawal [32].The MSE analysis performed for the SAP time series, reported in Figure4, evidences a clear increase in the MSE computed during MA. This result, which is again more evident using LMSE than RMSE, has been previously observed at scale 1 [36] and is possibly related to complex patterns of autonomic activation following cognitive load [41, 42]; our results, showing the persistence of higher MSE values up to long time scales (fτs=0.15Hz), point to an involvement of LF oscillations in the generation of higher complexity during mental stress. The multiscale pattern of complexity induced by postural stress is more complicated, indicating lower MSE during HUT than during REST for short time scales (fτs=0.5 Hz), no significant alterations for intermediate time scales (up to fτs=0.275 Hz), and higher MSE during HUT than during REST for longer time scales (fτs between 0.275 Hz and 0.15 Hz). The emergence of a significantly higher complexity of SAP during HUT when HF fluctuations are filtered out is a novel finding, suggesting that the postural stress decreases the regularity of LF arterial pressure oscillations, while HF oscillations display similar or even increased regularity.The multiscale complexity analysis of respiration variability, reported in Figure5, documents that MSE is significantly decreased during HUT and significantly increased during MA. These variations are consistently observed across multiple scales, going up to fτs=0.3 Hz for the lower LMSE and RMSE during HUT and up to fτs=0.15 Hz for the higher LMSE during MA. The higher regularity of the respiratory dynamics during postural stress, which seems to be confined to the HF band where respiration usually occurs, may be related to an increased tidal volume during this condition. On the other hand, the lower regularity of respiration during mental stress might be explained by the appearance of long pauses or sighs in the respiration pattern while performing mental arithmetics [43], making it more erratic and thus complex. The existence of more erratic patterns, which likely span a broad band of the frequency spectrum, may also explain the fact that we observe higher complexity of respiration up to long time scales, going beyond the common frequency bands at which respiratory activity is observed.


## 5. Conclusions and Future Developments

The present study introduces for the first time a multiscale entropy measure that is based on theoretical rather than empirical grounds and can thus be analytically computed from the parametric representation of an observed stochastic process. As a matter of fact, the proposed LMSE method is highly data-efficient because it stems from simple linear parametric modeling and is thus much more reliable than MSE or its modifications [6], including RMSE, in assessing the complexity of short time series at long time scales. Compared with another recently proposed method to assess multiscale complexity from short data sequences [32], LMSE shares the philosophy of being derived for linear AR processes, but differs in the fact that it is designed closely following the two algorithmic steps of MSE computation, that is, low-pass filtering and downsampling. The high computational reliability of LMSE comes from the fact that filtering and downsampling are not actually implemented on the measured time series, but result from the analytical quantification of the impact that they have on the state-space parameters of the observed AR process.Our approach can be fruitfully exploited, as done in the present study, to relate the multiscale complexity of a stochastic process to the parameters that establish its dynamical features, or to estimate patterns of multiscale complexity from short process realizations. In this work, we have formalized the dependence of the multiscale complexity of an AR process on the amplitude and frequency of its stochastic oscillatory components and have assessed multiscale patterns of short-term cardiovascular complexity which cannot be fully retrieved using standard MSE methods. Our results emphasize the role of the sympathetic control in driving the increased regularity of low-frequency heart rate oscillations and the increased complexity of low-frequency arterial pressure oscillations during postural stress. LMSE analysis stresses also the importance of dynamics occurring within the low-frequency band in determining the increased complexity of arterial pressure, and even that of respiration, during mental stress.The main strength of the proposed approach, that is, the linear parametric formulation, constitutes also one of its major limitations. In fact, the computation of LMSE holds exactly only if the observed process has a Gaussian distribution; in such a case, the linear AR description fully captures all of the variability in the process that determines the measured entropy rates, and model-free formulations as the one implemented by SampEn have no additional utility [25]. On the contrary, departures from linearity leading to non-Gaussian distributions may generate dynamics which are captured only partially by LMSE, and thus important features of multiscale complexity may be missed in this case. The suitability of LMSE for the applicative context of this study is supported by the fact that linear methods are ubiquitously exploited in short-term cardiovascular variability studies (e.g., performing parametric spectral analysis [17, 44] or parametric coupling and causality analysis [37, 45]). Nevertheless, future studies should deepen the comparison of LMSE with RMSE or other nonparametric MSE techniques in order to clarify the importance of accounting for nonlinear dynamics in the multiscale analysis of biomedical time series. On the other hand, keeping the linear analysis framework, a desirable extension of the present study would be to integrate the standard AR representation with fractionally integrated (FI) innovation modeling [46]; this would allow the computation of multiscale complexity for ARFI processes, thus incorporating in the approach the ability to describe long-range correlations phenomena which are typically assessed in long-term multiscale analyses [47].



---

*Source: 1768264-2017-12-07.xml*