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.
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 process with zero mean and variance . Let us further define as the variable obtained sampling the process at the discrete time 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 , which averages consecutive samples to yield the rescaled process : the rescaled process is denoted here as because it is a downsampled version of the original process ; in fact, in (1) one sample of the process is obtained by taking one sample every values of the averaged version of , 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, , given its past states, . In MSE, this is accomplished approximating the past history of the rescaled process with the finite-length pattern 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 -dimensional space remain similar also in the -dimensional space when they are incremented with their future values, that is, the probability that and are similar given that and are similar; similarity is assessed through a coarse-graining procedure which applies the Heaviside step function with parameter on the distance between patterns (i.e., is a threshold such that and are similar if their distance in the -dimensional space is lower than , where the distance is assessed using the maximum norm).
The MSE measure resulting from the above procedure, which we denote as , quantifies the dynamical complexity of the original process observed at scale . The free parameters of the MSE estimator are the length of the patterns used to approximate the past of the process and the threshold that sets the similarity between patterns. In the application of MSE, the pattern length is limited to a few samples (typically, ), while the threshold distance is a fraction of the standard deviation of the original process (typically, ).
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 process x is performed first applying a low-pass filter to obtain the filtered process and then reducing the sampling rate of the filtered process to obtain the rescaled (downsampled) process xd(n). The two steps yield, respectively, the processeswhere and are the filter coefficients and and 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 and (the filter is nonrecursive) and with coefficients for each ; the cutoff frequency of the filter is constrained to the value . 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 (with ) is implemented in which the coefficients and are set to realize a low-pass implementation with cutoff frequency [11]. This choice brings also the advantage that since the cutoff frequency can take any real value, in RMSE the scale factor 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 parameter is constant for all scale factors, as it is set at a fraction of the standard deviation of the original process (e.g., ), in RMSE this parameter is set at a fraction of the standard deviation of the rescaled process (e.g., ). 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 process 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 , defined aswhere is the order of the process, , , are linear regression coefficients describing the interactions between the process variables as a function of the lag , and is an innovation process with variance . The process is fully described by the AR model (3) when the innovation process e 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] and the entropy rate of the process, that is, the conditional entropy of the present given the past , is related to the variance of the innovations by [28] 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 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 , by setting an integer upsampling factor and by defining the processwhere is the order of the upsampled process and the coefficients are set to be a zero-padded version of the original AR coefficients; that is, we set for each and for each , . Note that the innovations for the original and upsampled processes are formally the same, that is, . Then, we filter the upsampled process using a linear FIR filter, that is, a filter implemented as in (2a) with for each ; with some algebraic manipulation we find that, substituting (7) in (2a) applied with and with in place of , the filtered representation of the upsampled process takes the formEquation (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 aswhere is a -dimensional state process, is the scalar SS innovation process defined as , and the vectors and and the matrix are defined asSuch a connection between ARMA and SS models allows finding, through (10) and the relation , the parameters of the SS process descriptive of the process . Such a process is a filtered version of the upsampled process , obtained passing through a FIR low-pass filter with cutoff frequency equal to ; therefore, since the sampling frequency of the upsampled process is times higher than that of the original process , the cutoff frequency of the low-pass filter applied to for obtaining becomes .
The next step of the procedure is to provide a representation for the downsampled process , where downsampling is applied to the filtered process . 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 can be represented in SS form involving a vector state process aswhere the parameters of the SS model (11) are the matrix , the vector , the variance of the scalar noise process , , the covariance matrix of the vector state noise process , , and the cross-covariance between and , . These parameters can be computed at scale , in terms of parameters previously obtained, as follows [20, 30]:Then, the SS model (11) can be transformed in a form similar to that of (9) which evidences the innovations:by solving a so-called discrete algebraic Riccati equation [20, 30]applied to the parameters of the SS model (11), to complete the derivation of the parameters of (13) asFinally, the variance of the downsampled process can be computed analytically solving a discrete-time Lyapunov equation [20] as
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 given by (3). Among these parameters, the relevant ones for our purpose are the variance of the downsampled process and that of the relevant innovations , which are used to compute the complexity of in analogy to (6), thus obtaining our MSE measure: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 , obtained setting integer values for the scale factor τ and the upsampling factor . As this corresponds to applying to the original process a low-pass filter with cutoff frequency , 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 from http://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 by pairs of complex-conjugate poles with assigned modulus and phase , (the order of the process is ). Type 1 simulations are designed with one pair of complex-conjugate poles , according to two configurations: (a) fixing the pole phase to (frequency ) and varying the modulus in the range ; (b) fixing the pole modulus to and varying the frequency in the range . Type 2 simulations are designed with two pairs of complex-conjugate poles , according to two configurations both with a fixed low-frequency pole obtained setting and ; (c) fixing the phase of the second pole to (frequency ) and varying the modulus in the range ; (d) fixing the modulus of the second pole to and varying the frequency in the range . 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 Section 2.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 ); these coefficients, together with the innovation variance , are set as parameters of the AR model of (3). From these parameters, MSE values are computed for specific values of the upsampling factor and the scale factor , that is, ,, , , chosen in order to yield an approximately uniform range of values for the cutoff frequency of the rescaling low-pass filter, that is, , 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, ; the filter order is set to .
After theoretical analysis, practical estimation of MSE was performed choosing two representative cases of the parameter setting for Type 1 simulation , and Type 2 simulation , ; , , 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 [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 [24] and calculating SampEn with embedding and tolerance parameters and ( is the variance of the downsampled signal).
3.2. Theoretical Analysis
Results of the theoretical analysis are reported in Figure 1. As a general result, the complexity of the simulated AR processes tends to increase at increasing the time scale (i.e., at decreasing the cutoff frequency of the low-pass filter applied to the original process, . 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., , corresponding to uncorrelated white noise).
(a) Theoretical MSE
(b) Theoretical MSE
(c) Theoretical MSE
(d) Theoretical MSE
Figures 1(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 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 when the process is a Gaussian white noise without regular dynamics. The behavior of MSE for oscillations with fixed amplitude and varying frequency is more complicated (Figure 1(b)). At scale one , the MSE is the same if the frequency of the oscillatory dynamics, , 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 (see the trends of versus and versus , in Figure 1(b)).
Figures 1(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 , , and increases with fixed ) 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 , , and increases with fixed ) 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
Figure 2 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); , ) and a pair of stochastic oscillations (Figures 2(c) and 2(d); , ; , ), 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 (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 , 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].
(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, 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 [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 . 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 and equal to 20% of the variance of the resampled signal. As for simulations, in the computation of both LMSE and RMSE, the parameters and determining the time scale were set to obtain the following values for the cutoff frequency of the resampling filter: , 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, .
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 Student -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 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 Figures 3, 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.
(a) Heart period, LMSE
(b) Heart period, RMSE
(a) Systolic pressure, LMSE
(b) Systolic pressure, RMSE
(a) Respiration, LMSE
(b) Respiration, RMSE
The analysis of LMSE and RMSE computed for the HP time series, reported in Figure 3, 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 Hz using RMSE (Figure 3(b)) and up to Hz using LMSE (Figure 3(a)). At scale 1 , 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 Figure 4, 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 , 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 ( Hz), no significant alterations for intermediate time scales (up to Hz), and higher MSE during HUT than during REST for longer time scales ( 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 Figure 5, documents that MSE is significantly decreased during HUT and significantly increased during MA. These variations are consistently observed across multiple scales, going up to Hz for the lower LMSE and RMSE during HUT and up to 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].
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
Luca Faes and Giandomenico Nollo were partly supported by funding from the Healthcare Research and Implementation Program (IRCS) of the Autonomous Province of Trento (PAT), Italy. Michal Javorka was supported by Grants APVV-0235-12, VEGA 1/0117/17, and VEGA 1/0202/16 and project “Biomedical Center Martin,” ITMS code 26220220187, the project cofinanced from EU sources.
Supplementary Materials
The LMSE MATLAB toolbox is available as supplementary material to this article. The package contains functions for the computation of multiscale entropy based on the linear state-space approach and on the refined approach proposed in [11], as well as scripts that allow performing multiscale complexity analysis for the simulated and real data considered in the present study. (Supplementary Materials)