*Article* **Multivariate and Multiscale Complexity of Long-Range Correlated Cardiovascular and Respiratory Variability Series**

**Aurora Martins 1,2, Riccardo Pernice 3,***∗***, Celestino Amado 1, Ana Paula Rocha 1,2, Maria Eduarda Silva 4,5, Michal Javorka 6,7 and Luca Faes <sup>3</sup>**


Received: 18 February 2020; Accepted: 09 March 2020 ; Published: 11 March 2020

**Abstract:** Assessing the dynamical complexity of biological time series represents an important topic with potential applications ranging from the characterization of physiological states and pathological conditions to the calculation of diagnostic parameters. In particular, cardiovascular time series exhibit a variability produced by different physiological control mechanisms coupled with each other, which take into account several variables and operate across multiple time scales that result in the coexistence of short term dynamics and long-range correlations. The most widely employed technique to evaluate the dynamical complexity of a time series at different time scales, the so-called multiscale entropy (MSE), has been proven to be unsuitable in the presence of short multivariate time series to be analyzed at long time scales. This work aims at overcoming these issues via the introduction of a new method for the assessment of the multiscale complexity of multivariate time series. The method first exploits vector autoregressive fractionally integrated (VARFI) models to yield a linear parametric representation of vector stochastic processes characterized by short- and long-range correlations. Then, it provides an analytical formulation, within the theory of state-space models, of how the VARFI parameters change when the processes are observed across multiple time scales, which is finally exploited to derive MSE measures relevant to the overall multivariate process or to one constituent scalar process. The proposed approach is applied on cardiovascular and respiratory time series to assess the complexity of the heart period, systolic arterial pressure and respiration variability measured in a group of healthy subjects during conditions of postural and mental stress. Our results document that the proposed methodology can detect physiologically meaningful multiscale patterns of complexity documented previously, but can also capture significant variations in complexity which cannot be observed using standard methods that do not take into account long-range correlations.

**Keywords:** multi-scale entropy (MSE); vector autoregressive fractionally integrated (VARFI) models; heart rate variability (HRV); systolic arterial pressure (SAP)

#### **1. Introduction**

Cardiovascular oscillations are influenced by the combined activity of different physiological regulation processes and, as a consequence, exhibit a rich dynamical structure [1]. Such different processes do not usually work at a single time-scale, but instead operate across different temporal scales, that for example reflect thermoregulatory or neural parasympathetic and sympathetic control. For this reason, different methods have been recently developed to evaluate the 'multiscale complexity' of cardiovascular oscillations. These methods allow both to characterize the physiological regulatory systems and to extract diagnostic parameters, and thus can have noteworthy clinical implications: in fact, a decrease of dynamical complexity can be related to an impaired capability of the subsystems composing the organism to interact among each other and it has been already proposed as a marker of pathological conditions [2,3]. Among the proposed methods, the one which is likely most popular is the so-called multiscale entropy (MSE), developed by Costa et al. [4]. This method calculates the conditional entropy (CE) of a single time series (usually the heart period, HP) as a function of the time scale of observation; the change of time scale is achieved through averaging consecutive segments of the time series via a procedure that has been lately recognized to correspond to a two-step process consisting of a low-pass filter followed by downsampling [5]. It is worth noting, however, that the initial formulation of multiscale entropy suffered from drawbacks related both to issues due to the filtering and downsampling steps, and to the unsuitability of CE analysis in conditions of data paucity caused by the availability of short time series and by the needs to explore multivariate time series at coarse time scales. Therefore, in the last years, the definition of MSE has been refined to take into account typical requirements of cardiovascular and cardiorespiratory signal analysis, and specifically: (i) to allow the joint calculation of complexity of multiple variables besides HP, for example systolic arterial pressure (SAP) and respiration (RESP) [6]; (ii) to allow the assessment of the complexity of shorter time series, usually few hundred beats long [7,8]. For short-term physiological time series, complexity has been related to the regularity of the temporal patterns observed in the signals, and thus is usually a measure of the unpredictability of the present sample given a limited number of past samples [9]. However, recent studies have recognized the importance of long-range correlations resulting in slowly varying dynamics also for the analysis of short-term complexity [10], and have started to account for these correlations in multiscale entopy-based analysis [8].

In this work, we introduce a novel method to compute multivariate and multiscale complexity of cardiovascular oscillations. This method fits a multivariate time series using a vector autoregressive fractionally integrated (VARFI) model and then provides the multiscale representation of the VARFI parameters using the theory of state space models, allowing to extract from such parameters multiscale and multivariate measures of complexity. This approach presents several advantages if compared to other previous works in the same field [4–8,10], and in particular: (i) the parametric formulation employed permits to work reliably on short time series; (ii) fractional integration allows to take into account not only short-term dynamics, but also the long-range correlations; (iii) the vector formulation permits to compute the overall complexity of multivariate time series, or the complexity of an individual time series when one or more other series are considered. In this work, such approach is applied to HP, SAP and RESP time series measured in healthy subjects monitored in a resting condition and during two types of physiological stress: postural stress provoked by head-up tilt and mental stress induced by mental arithmetics.

The algorithms for multivariate multiscale analysis of physiological time series presented in this work are collected in the MSE-VARFI MATLAB toolbox, which can be freely downloaded from http://www.lucafaes.net/LMSE-MSE\_VARFI.html.

#### **2. Methods**

#### *2.1. Measures of Complexity in Linear Multivariate Stochastic processes*

Considering a dynamical system X whose activity is defined by the zero-mean stationary multivariate stochastic process **X** = [*X*<sup>1</sup> ··· *XM*], where each *Xj*, *j* = 1, ... , *M*, denotes a scalar stochastic process *Xj* = [··· *Xj*,1 ··· *Xj*,*<sup>n</sup>* ··· ], let us define as **X***<sup>n</sup>* = [*X*1,*<sup>n</sup>* ··· *XM*,*n*] the *M*−dimensional random variable describing the present state of the system, and as **X**− *<sup>n</sup>* = [**X***n*−1**X***n*−2...] the infinite-dimensional vector variable describing its past states. Then, a measure of complexity of the system, typically defined for univariate systems [10] and here extended to the multivariate case, is the entropy rate defined as

$$\mathbf{C}\_{\mathbf{X}} = H(\mathbf{X}\_n | \mathbf{X}\_n^-) = H(\mathbf{X}\_{n+1}^-) - H(\mathbf{X}\_n^-), \tag{1}$$

where *H*(**X***n*|**X**<sup>−</sup> *<sup>n</sup>* ) is the conditional entropy of the present given the past, with *H*(·) denoting the Shannon Entropy. If the observed process **X** is a jointly distributed Gaussian process, it can be described without loss of generality through a vector linear regression model fed by white and uncorrelated innovations **<sup>E</sup>***<sup>n</sup>* = [*E*1,*<sup>n</sup>* ··· *EM*,*n*] such that, for each *<sup>j</sup>* <sup>∈</sup> 1, . . . , *<sup>M</sup>*, *Ej*,*<sup>n</sup>* <sup>=</sup> *Xj*,*<sup>n</sup>* <sup>−</sup> <sup>E</sup>[*Xj*,*n*|**X**<sup>−</sup> *<sup>n</sup>* ] [11]. In such a case, the entropy of the present state of the process and the conditional entropy of the present given the past can be expressed analytically in terms of the covariance of the process, <sup>Σ</sup>**<sup>X</sup>** = E[**X***<sup>T</sup> <sup>n</sup>***X***n*], and the covariance of the innovations, <sup>Σ</sup>**<sup>E</sup>** = E[**E***<sup>T</sup> <sup>n</sup>***E***n*], as [11–13]:

$$H(\mathbf{X}\_n) = \frac{1}{2} \ln \left( (2\pi \mathbf{e})^M |\Sigma\_\mathbf{X}| \right),\tag{2a}$$

$$H(\mathbf{X}\_n|\mathbf{X}\_n^-) = \frac{1}{2} \ln((2\pi\mathbf{e})^M |\Sigma\_\mathbf{E}|),\tag{2b}$$

where |·| stands for matrix determinant. Then, it follows immediately that the complexity of a multivariate linear process with joint Gaussian distribution is given by Equation (2b). In this work we provide an alternative definition of complexity, which includes a normalization of the innovation covariance to the process covariance:

$$\overline{\mathbf{C}}\_{\mathbf{X}} = \frac{1}{2} \ln \left( (2\pi \mathbf{e})^{M} \frac{|\Sigma\_{\mathbf{E}}|}{|\Sigma\_{\mathbf{X}}|} \right) . \tag{3}$$

Such normalization, which is implicitly implemented in studies assessing the complexity of univariate time series where the series is normalized to unit variance before computing the complexity measure, is formalized here in Equation (3) for multivariate series, and will be fundamental in the definition of multiscale complexity where the process covariance intrinsically changes with the time scale.

The measure of global complexity defined above for a multivariate process can be particularized to the characterization of one of its constituent processes. Specifically, the complexity of a scalar process *Xj* ∈ **X**, *j* ∈ 1, . . . *M*, can be defined in an univariate context with respect to its own dynamics only, in a bivariate context with respect to its dynamics and to the dynamics of another scalar process *Xi* ∈ **X**, *i* = *j*, or in a fully multivariate context with respect to the dynamics of the whole observed vector process **X** [14]. The derivations are based on the knowledge that, in Gaussian processes, a linear parametric representation captures all of the entropy differences that define the conditional entropies [11] and that these entropy differences are related to the partial variances of the present of the target given its past and the past of one or more other processes, intended as variances of the prediction errors resulting from linear regression [13,15]. Specifically, let us denote as *Ej*|*j*,*<sup>n</sup>* <sup>=</sup> *Xj*,*<sup>n</sup>* <sup>−</sup> <sup>E</sup>[*Xj*,*n*|*X*<sup>−</sup> *j*,*n*] and *Ej*|*ij*,*<sup>n</sup>* <sup>=</sup> *Xj*,*<sup>n</sup>* <sup>−</sup> <sup>E</sup>[*Xj*,*n*|*X*<sup>−</sup> *<sup>i</sup>*,*n*, *X*<sup>−</sup> *<sup>j</sup>*,*n*] the prediction error of a linear regression of *Xj*,*<sup>n</sup>* performed respectively on *X*− *<sup>j</sup>*,*<sup>n</sup>* and (*X*<sup>−</sup> *<sup>j</sup>*,*n*, *X*<sup>−</sup> *<sup>i</sup>*,*n*), and consider that the prediction error of a linear regression of *Xj*,*<sup>n</sup>* on **X**<sup>−</sup> *<sup>n</sup>* is the *j th* innovation process *Ej*,*<sup>n</sup>* defined above. Then, denoting the variances of the three

prediction errors as <sup>Σ</sup>*Ej*|*<sup>j</sup>* <sup>=</sup> <sup>E</sup>[*E*<sup>2</sup> *<sup>j</sup>*|*j*,*n*], <sup>Σ</sup>*Ej*|*ij* <sup>=</sup> <sup>E</sup>[*E*<sup>2</sup> *<sup>j</sup>*|*ij*,*n*] and <sup>Σ</sup>*Ej* <sup>=</sup> <sup>E</sup>[*E*<sup>2</sup> *<sup>j</sup>*,*n*] = Σ**E**(*j*, *j*), the univariate, bivariate and multivariate normalized complexity measures relevant to the process *Xj* are defined as:

$$\overline{\mathbb{C}}\_{X\_{\overline{j}}|X\_{\overline{j}}} = \frac{1}{2} \ln 2\pi e \frac{\Sigma\_{E\_{\overline{j}\overline{j}}}}{\Sigma\_{X\_{\overline{j}}}},\tag{4a}$$

$$\overline{C}\_{X\_{\bar{j}}|X\_{i},X\_{\bar{j}}} = \frac{1}{2} \ln 2\pi a \frac{\Sigma\_{E\_{\bar{j}|\bar{j}}}}{\Sigma\_{X\_{\bar{j}}}},\tag{4b}$$

$$\overline{C}\_{X\_{\!\!\!/]\!X}} = \frac{1}{2} \ln 2\pi e \frac{\Sigma\_{E\_j}}{\Sigma\_{X\_j}},\tag{4c}$$

where Σ*Xj* = Σ**X**(*j*, *j*) is the *j th* diagonal element of Σ**X**.

#### *2.2. Linear Multivariate Stochastic Processes with Long Range Correlations*

In this section we present the parametric approach to the description of linear multivariate Gaussian stochastic processes exhibiting both short-term dynamics and long-range correlations. The approach is based on representing the *M*-dimensional discrete-time, zero-mean and unit variance stochastic process **X***n* defined in Section 2.1 as a vector autoregressive fractionally integrated (VARFI) process fed by the uncorrelated Gaussian innovations **E***n*. The VARFI process takes the form [16]:

$$\mathbf{A}(L)\text{diag}(\nabla^{\mathbf{d}})\mathbf{X}\_{\mathbf{n}}=\mathbf{E}\_{\mathbf{n}\prime}\tag{5}$$

where *L* is the back-shift operator (*L<sup>i</sup>* **<sup>X</sup>***<sup>n</sup>* <sup>=</sup> **<sup>X</sup>***n*−*i*), **<sup>A</sup>**(*L*) = **<sup>I</sup>***<sup>M</sup>* <sup>−</sup> *<sup>p</sup>* ∑ *i*=1 **A***iL<sup>i</sup>* (**I***<sup>M</sup>* is the identity matrix of size *M*) is a vector autoregressive (VAR) polynomial of order *p* defined by the *M* × *M* coefficient matrices **A**1,..., **A***p*, and

$$\text{diag}(\nabla^{\mathbf{d}}) = \begin{bmatrix} (1-L)^{d\_1} & 0 & \dots & 0\\ 0 & (1-L)^{d\_2} & \dots & 0\\ \vdots & \vdots & & \vdots\\ 0 & 0 & \dots & (1-L)^{d\_M} \end{bmatrix}^\prime$$

where (<sup>1</sup> − *<sup>L</sup>*)*di* , *<sup>i</sup>* = 1, . . . , *<sup>M</sup>*, is the fractional differencing operator defined by [17]:

$$(1 - L)^{d\_i} = \sum\_{k=0}^{\infty} G\_k^{(i)} L^k \quad , \quad G\_k^{(i)} = \frac{\Gamma(k - d\_i)}{\Gamma(-d\_i)\Gamma(k + 1)} \tag{6}$$

with Γ(·) denoting the gamma (generalized factorial) function. The parameter **d** = (*d*1, ... , *dM*) in Equation (5) determines the long-term behavior of each individual process, while the coefficients of the polynomial **A**(*L*) allow description of the short-term dynamics. Note that the process defined in Equation (5) is a particular case of the broader class of VARFIMA(*p*, **d**, *l*) processes [16], which also contains the class of autoregressive processes VAR(p); here we restrict our analysis to the description of the VARFIMA(*p*, **d**, 0) process, which we denote as a VARFI(*p*, **d**) process.

The parameters of the VARFI model (5), namely the coefficients of **A**(*L*) and the variance of the innovations Σ**E**, are typically obtained from process realizations of finite length first estimating the differencing parameters *di* by means of the Whittle semi-parametric local estimator [17] separately for each individual process *Xi*, then defining the filtered data *<sup>X</sup>*(*f*) *<sup>i</sup>*,*<sup>n</sup>* = (<sup>1</sup> − *<sup>L</sup>*)*diXi*,*n*, and finally estimating the VAR parameters from the filtered data **<sup>X</sup>**(*f*) *<sup>n</sup>* using the ordinary least squares method to solve the VAR model **<sup>A</sup>**(*L*)**X**(*f*) *<sup>n</sup>* <sup>=</sup> **<sup>E</sup>***n*, with model order *<sup>p</sup>* assessed through the Bayesian information criterion [18].

#### *2.3. Multiscale Complexity of VARFI Processes*

In this section we report how to compute across multiple temporal scales the complexity measures defined in Section 2.1, under the hypothesis that the analyzed multivariate process is properly described by the VARFI representation provided in Section 2.2. The procedure for multiscale complexity analysis, which extends the approach proposed in Reference [19] to multivariate processes incorporating long-range correlations, is presented here reporting the essential steps and described with more mathematical details in the Appendices. There are three appendices to the main text. Appendix A reports the procedure for computing the coefficients of a finite-order VAR process that approximates an assigned VARFI process. Appendix B recalls the derivations relevant to the identification of multiscale state space (SS) processes defined as filtered and downsampled versions of an assigned VAR process; the parameters of the SS process defined at an assigned time scale are exploited to compute the multivariate complexity measure at that time scale. Appendix C describes how to define restricted SS processes and to rearrange them in order to extract the partial variances needed for the computation of the univariate and bivariate multiscale complexity measures. The derivations described in Appendices B and C are taken from Refs. [7,8,19,20].

Before implementing the rescaling procedure, we approximate the VARFI process (5) with a finite order VAR process by truncating the fractional integration part at a finite lag *q*, that is, we perform the following approximation:

$$\text{diag}(\nabla^{\mathbf{d}}) \approx \mathbf{G}(L) = \sum\_{k=0}^{q} \mathbf{G}\_k L^k,\tag{7a}$$

$$\mathbf{G}\_k = \text{diag}\left[\mathbf{G}\_k^{(1)}, \dots, \mathbf{G}\_k^{(M)}\right]. \tag{7b}$$

This allows us to rewrite the VARFI(*p*, **d**) process as a VAR(*m*) process, where *m* = *p* + *q*:

$$\mathbf{B}(L)\mathbf{X}\_{\mathbb{N}} = \mathbf{E}\_{\mathbb{N}},\tag{8}$$

where the new coefficients are **B**(*L*) = **A**(*L*)**G**(*L*), with **A**(*L*) as in Equation (5) and **G**(*L*) as in Equation (7a). The exact procedure to derive the coefficients of the VAR(*m*) process is explained in Appendix A.

In order to represent a scalar stochastic process at the temporal scale defined by the scale factor *τ*, a two step procedure is typically employed which consists first in filtering the process with a lowpass filter with cutoff frequency *f<sup>τ</sup>* = 1/(2*τ*), and then downsampling the filtered process using a decimation factor *τ* [5,19]. Extending this approach to the multivariate process **X**, we first implement the following linear finite impulse response (FIR):

$$\mathbf{X}\_n^{(r)} = \mathbf{D}(L)\mathbf{X}\_n \quad , \tag{9}$$

where *<sup>r</sup>* denotes the filter order and **<sup>D</sup>**(*L*) = *<sup>r</sup>* ∑ *k*=0 **I***MDkL<sup>k</sup>*, where the coefficients of the polynomial *Dk*, *k* = 1, ... ,*r*, are the same for all scalar processes *Xj* ∈ **X** and are chosen to set up a lowpass FIR configuration with cutoff frequency 1/(2*τ*). The filtering step transforms the VAR(*p+q*) process (8) into a VARMA(*p+q,r*) process with moving average (MA) part determined by the filter coefficients:

$$\mathbf{B}(L)\mathbf{X}\_n^{(r)} = \mathbf{D}(L)\mathbf{B}(L)\mathbf{X}\_{\mathbb{R}} = \mathbf{D}(L)\mathbf{E}\_{\mathbb{R}} \quad . \tag{10}$$

Then, we exploit the connection between VARMA processes and state space (SS) processes [21] to evidence that the VARMA process (10) can be expressed in SS form as:

$$\mathbf{Z}\_{n+1}^{(r)} = \mathbf{B}^{(r)} \mathbf{Z}\_{n}^{(r)} + \mathbf{K}^{(r)} \mathbf{E}\_{n}^{(r)} \tag{11a}$$

$$\mathbf{X}\_{n}^{(r)} = \mathbf{C}^{(r)} \mathbf{Z}\_{n}^{(r)} + \mathbf{E}\_{n}^{(r)},\tag{11b}$$

where **<sup>Z</sup>**(*r*) *<sup>n</sup>* = [**X**(*r*) *<sup>n</sup>*−<sup>1</sup> ··· **<sup>X</sup>**(*r*) *<sup>n</sup>*−*m***E***n*−<sup>1</sup> ··· **<sup>E</sup>***n*−*r*] *<sup>T</sup>* is a (*<sup>m</sup>* + *<sup>r</sup>*)-dimensional state process, **<sup>E</sup>**(*r*) *<sup>n</sup>* = **D**0**E***<sup>n</sup>* is the SS innovation process, and the vectors **K**(*r*) and **C**(*r*) and the matrix **B**(*r*) can be obtained from **B**(*L*) and **D**(*L*) (see Appendix B).

The second step of the rescaling procedure is to downsample the filtered process in order to complete the multiscale representation. This is done exploiting theoretical findings [7,22,23] which allow to describe the filtered SS process after downsampling in the form:

$$\mathbf{Z}\_{n+1}^{(\tau)} = \mathbf{B}^{(\tau)} \mathbf{Z}\_{n}^{(\tau)} + \mathbf{K}^{(\tau)} \mathbf{E}\_{n}^{(\tau)} \tag{12a}$$

$$\mathbf{X}\_{n}^{(\tau)} = \mathbf{C}^{(\tau)} \mathbf{Z}\_{n}^{(\tau)} + \mathbf{E}\_{n}^{(\tau)}.\tag{12b}$$

Equation (12) provides the SS form of the filtered and downsampled version of the original VARFI(*p*, **d**) process, and its parameters (**B**(*τ*) , **C**(*τ*) , **K**(*τ*) , Σ**E**(*τ*)) can be obtained from the SS parameters before downsampling and from the downsampling factor *τ*; moreover, the variance of the downsampled process, Σ**X**(*τ*), can be computed analytically from the parameters of the SS model (12) by solving a discrete-time Lyapunov equation (these steps are shown in the Appendix B). The parameters relevant to the computation of complexity at scale *τ* are the variance of the downsampled process, Σ**X**(*τ*), and the variance of the corresponding innovations, Σ**E**(*τ*). These variances can indeed be combined in a similar way to that of Equation (3) to yield the expression of the complexity of the original process **X***n* when it is observed at scale *τ*:

$$\nabla\_{\mathbf{X}} = \frac{1}{2} \ln \left( (2\pi \mathbf{e})^M \frac{|\boldsymbol{\Sigma}\_{\mathbf{E}^{(\tau)}}|}{|\boldsymbol{\Sigma}\_{\mathbf{X}^{(\tau)}}|} \right) . \tag{13}$$

Note that this measure tends to the theoretical value <sup>1</sup> <sup>2</sup> ln(2*πe*)*<sup>M</sup>* when *<sup>τ</sup>* → <sup>∞</sup>.

Finally, we show how to compute any partial variance appearing in Equation (4) from the parameters of an SS model in the form of (12), so that the three complexity measures relevant to the scalar process *Xj* can be computed at any assigned time scale *τ*. To do this, we exploit the formulations reported in Refs. [22,23], showing that the partial variance <sup>Σ</sup>*E*(*τ*) *j*|*a* , where the subscript *a* denotes any combination of indexes ∈ 1, . . . , *M*, can be derived from the SS representation of the innovations of a submodel obtained removing the variables not indexed by *a* from the observation equation. Specifically, we need to consider the submodel with state Equation (12a) and observation equation:

$$\mathbf{X}\_{a,n}^{(\tau)} = \mathbf{C}\_{a}^{(\tau)} \mathbf{Z}\_{n}^{(\tau)} + \mathbf{E}\_{a,n}^{(\tau)},\tag{14}$$

where the additional subscript *a* denotes the selection of the rows with indices *a* in a vector or a matrix. This restricted SS model can be rearranged to extract the partial variance <sup>Σ</sup>*E*(*τ*) *j*|*a* from the covariance matrix of its innovations (see Appendix C), so that with this procedure the univariate, bivariate and multivariate normalized complexity measures can be computed inserting the partial variances <sup>Σ</sup>*E*(*τ*) *j*|*j* , in Equation (4a) and Equation (4b), and using Σ(*τ*)

<sup>Σ</sup>*E*(*τ*) *j*|*ij Ej* = <sup>Σ</sup>**E**(*τ*)(*j*, *<sup>j</sup>*) in Equation (4c), together with Σ(*τ*) *Xj* = <sup>Σ</sup>**X**(*τ*)(*j*, *<sup>j</sup>*) from Equation (A9b). These individual measures tend to the theoretical value 1 <sup>2</sup> ln(2*πe*) when *τ* → ∞.

#### **3. Application to Cardiovascular Variability Processes**

In this section the proposed method is illustrated using cardiovascular and respiratory series: the heart period (HP), systolic arterial pressure (SAP), and respiration (RESP). Several studies report the interaction between the dynamics of these three time series [24–28], which motivates their use in a multivariate context. The variation of heart period, usually referred as heart rate variability (HRV), reflecting cardiovascular complexity and representing the capability of the organism to react to environmental and psychological stimuli, is the most studied variable and main target variable in cardiovascular and cardiorespiratory spontaneous variability [29–31]. For this reason, the conditional measures of a single scalar process will focus on the HP time series, as it has been shown that SAP and RESP have an effect (direct or indirectly) on this process [24–28].

#### *3.1. Experimental Protocol*

The HP, SAP and RESP time series were measured in a group of 62 healthy subjects (19.5 ± 3.3 years old, 37 females) monitored in the resting supine position (SU1), in the upright position (UP) reached through passive head-up tilt, in the recovery supine position (SU2) and during mental stress induced by mental arithmetics (MA) [32]. During the measurements, the subjects were free-breathing. For each subject and condition, the multivariate process **X** is defined as **X** = [*X*HP, *X*SAP, *X*RESP]. The acquired signals were the surface electrocardiogram (ECG), the finger arterial blood pressure recorded noninvasively by the photoplethysmographic method, and the respiration signal recorded through respiratory inductive plethysmography. For each subject and experimental condition, the values of HP, SAP and RESP were measured on a beat-to-beat basis respectively as the sequences of the temporal distances between consecutive R peaks of the ECG, the maximum values of the arterial pressure waveform taken within the consecutively detected heart periods, and the values of the respiratory signal sampled at the onset of the consecutively detected heart periods. The study was approved by Ethical Committee of the Jessenius Faculty of Medicine, Comenius University and all participants signed a written informed consent. A detailed description of the experimental protocol and signal measurement is reported in Ref. [32].

The analysis was performed on segments of at least 400 consecutive points, free of artifacts and deemed as weak-sense stationary through visual inspection, extracted from the time series for each subject and condition. Missing values and outliers were corrected through linear interpolation and, for HP and when possible, erroneous/missing intervals were substituted by pulse intervals measured as the difference in time between two consecutive SAP measurements (Δ*t*SAP (*n*) = *t*SAP(*n* + 1) − *t*SAP(*n*)). The three time series were normalized to zero mean and unit variance before multiscale analysis.

#### *3.2. Data Analysis*

Two different approaches were followed to compute multiscale complexity: (i) the "eVAR" approach, based on pure VAR model identification, that is, performing the whole procedure described in Section 2 after forcing **d** = [0, 0, 0] in Equation (5); (ii) the "eVARFI" approach, based on complete VARFI model identification, that is, following the whole procedure described in Section 2 with *di*, *i* = 1, ... , 3, estimated individually from the original time series and considered in the computations. Pursuing these approaches we compare, respectively, (i) the traditional complexity analysis where long-range correlations are neither removed nor modeled, and (ii) the complexity analysis performed by modeling the long-range correlations and considering them together with the short-term dynamics. Such a comparison is meant to infer the role of long-range correlations in contributing to the multiscale complexity and to its variation between conditions. The VARFI model fitting the time series **X** was identified first estimating the fractional differencing parameter *di*, *i* = 1, ... , 3, individually for each time series using the Whittle estimator, then filtering the time series with the fractional integration polynomial truncated at a lag *q* = 50, and finally estimating the parameters of the polynomial relevant

to the short-term dynamics through least squares VAR identification. The value of *q* has to be selected to approximate the VARFI process, which is theoretically of infinite order, with a finite order VAR process. According to previous studies [8,33] *q* = 50 is an appropriate value for truncating the VARFI process. By increasing *q*, we can obtain a more precise approximation of the fractional integration part but with a higher computational cost, while a reduced value (and thus an excessive truncation) causes an underestimation of the complexity and the smoothing of the nonmonotonic trends with the timescale [8]. The VAR model order *p* was selected as the minimum of the Bayesian information criterion (BIC) figure of merit [34]. Then, multiscale complexity measures were computed implementing a FIR lowpass filter of order *r* = 48, for time scales *τ* in the range (1,. . . ,30), which corresponds to lowpass cutoff normalized frequencies *f<sup>τ</sup>* = (0.5, ... , 0.01667). The order of the FIR filter determines its selectivity around the cutoff frequency. In this study, an order *r* = 48 was selected according to previous settings [8]. The changes in complexity related to the oscillatory rhythms are typically evaluated in cardiovascular variability in low-frequency (LF, from 0.04 Hz to 0.15 Hz) and high-frequency (HF, from 0.15 Hz to 0.4 Hz) spectral bands [29]. To account for the different mean heart rate for each subject and condition, multiscale analysis was performed determining the cutoff frequency of the rescaling filter based on the Nyquist frequency in Hz of the HP time series, which for the time scale *τ* is given by *f<sup>τ</sup>* Hz = 1/(2*τ*HP), where HP stands for the mean of the HP process. In this work, considering the physiological LF and HF frequency ranges, four cutoff frequencies *f* Hz were chosen to perform the analysis: 0.4 Hz, 0.15 Hz, 0.1 Hz and 0.04 Hz. To obtain the complexity values at such frequencies, linear interpolation was performed on the profiles *f<sup>τ</sup>* Hz .

The differencing parameters *di* were estimated individually for each time series in the interval [−0.5, 1] since the VARFI model is stationary for −0.5 < *di* < 0.5, while it is nonstationary but mean reverting for 0.5 ≤ *di* < 1. As such, the subjects with an estimation of *di* ≈ 1 in at least one condition were removed given that the series is possibly non mean reverting; three subjects were removed, so that a total of 59 subjects was considered for the statistical analysis.

#### *3.3. Statistical Analysis*

Significant changes in complexity across the pairs of experimental conditions SU1 vs. UP and SU2 vs. MA were assessed via a linear mixed-effects model, that is, a linear regression model that incorporates both fixed and random effects [35]. In our case, the fixed-effects (or factors) were condition and scale, while the random-effect was the subject-dependent intercept that allows for the random variation between subjects. Additionally, the interaction between the factors is also considered. The significance of both the effects (fixed and random) and the interaction between fixed effects was assessed by the significance of the corresponding estimated coefficients at the level of *p* < 0.05. Residuals were checked for whiteness.

To evaluate the changes of interest, estimated marginal means (EMM) [36] were computed for each difference, or contrast, UP-SU1 and MA-SU2, at each frequency of interest (0.4 Hz, 0.15 Hz, 0.1 Hz and 0.04 Hz). A Z-test was applied to check the significance of these differences with a significance level of *p* < 0.05 and Tukey correction for multiple comparisons. The packages lme4 [37] and emmeans [38] of the R software [39] were used to build the model and to compute EMM, respectively.

#### **4. Results**

This section presents the results of multiscale analysis performed for the multivariate complexity *C***<sup>X</sup>** as defined in Equation (3), as well as for the complexity computed for the scalar process *X*HP in four different ways: the univariate complexity of HP with respect to its own dynamics, *CX*HP|*X*HP (Equation (4a)), the bivariate complexity of HP with respect to its own dynamics and to the dynamics of SAP, *CX*HP|*X*SAP,*X*HP , or RESP, *CX*HP|*X*RESP,*X*HP (Equation (4b)), and the multivariate complexity of HP with respect to the dynamics of the whole trivariate process, *CX*HP|**<sup>X</sup>** (Equation (4c)).

Figure 1 presents the median and quartiles across subjects of the multivariate complexity *C***<sup>X</sup>** computed for eVAR (first row) and eVARFI (second row) as a function of the time scale *τ* = 1, . . . , 30,

for SU1 vs. UP (left column) and SU2 vs. MA (right column). The measure always tends to its theoretical value <sup>1</sup> <sup>2</sup> ln(2*πe*)<sup>3</sup> when computed for eVAR at high values of *<sup>τ</sup>*. In fact, the complexity value is in general lower for eVARFI than for eVAR, and presents more variability across subjects. From a visual inspection of the multiscale patterns one can infer that for eVAR there is a decrease in complexity moving from SU1 to UP, evident at short time scales (*τ* < 10), and that the complexity seems not to change while moving from SU2 to MA. For eVARFI, the decrease from SU1 to UP is visible across the whole range of time scales and there seems to be a decrease from SU2 to MA for scales larger than 4.

**Figure 1.** Distribution across subjects of the multivariate complexity measure *C***<sup>X</sup>** as a function of the time scale *τ* for eVAR (first row) and eVARFI (second row), for SU1 vs. UP (left column) and SU2 vs. MA (right column).

Figure 2 presents the same as the previous figure but for the univariate (Figure 2a), bivariate (Figure 2b with SAP and Figure 2c with RESP) and multivariate (Figure 2d) complexity measures relevant to the HP time series. Again, eVARFI estimation presents more variability than eVAR, visible across all measures, and does not reach the theoretical value <sup>1</sup> <sup>2</sup> ln(2*πe*) at long time scales. Although the four measures are similar with each other, slight changes occur when the complexity of HP is computed accounting for the dynamics of the other time series. When compared to the univariate measure *CX*HP|*X*HP , the median complexity value decreases for the bivariate measures at lower time scales, being lower with RESP (*CX*HP|*X*RESP,*X*HP ) than with SAP (*CX*HP|*X*SAP,*X*HP ); the multivariate measure *CX*HP|**<sup>X</sup>** presents even lower median complexity values for lower time scales. For both eVAR and eVARFI, the complexity of HP decreases at short time scales for all four measures. The differences are more subtle during MA, where only a slight decrease in the median is noticeable at very short time scale for eVARFI.

Figure 3 reports the distribution of the multivariate complexity measure *C***<sup>X</sup>** (blue dots) (values and boxplot) computed at the four predetermined cutoff frequencies of the multiscale filter (0.4 Hz, 0.15 Hz, 0.1 Hz and 0.04 Hz) for each experimental condition. Statistically significant changes (*p* < 0.05) in complexity across the pairs SU1 vs UP or SU2 vs MA, as assessed by the estimated marginal means based on the linear mixed-effects model, are marked with ∗. Comparing SU and UP, the multivariate

complexity decreases significantly both for eVAR and eVARFI at 0.4 Hz, while no significant differences are observed at lower cutoff frequencies. Moving from SU2 to MA, the measure increases significantly for eVAR at all frequencies except 0.04 Hz, but decreases for eVARFI at frequencies 0.1 Hz and 0.04 Hz.

**Figure 2.** Distribution across subjects of the univariate (**a**), bivariate (**b**) with SAP and (**c**) with RESP, and multivariate (**d**) complexity measures as a function of the time scale *τ* for eVAR (first row) and eVARFI (second row), for SU1 vs. UP (left column) and SU2 vs. MA (right column).

**Figure 3.** Distribution of the multivariate complexity measure *C***X**, depicted as boxplot (mean and confidence intervals, yellow filled box; standard deviation, black vertical line) and original values (dots) for selected frequencies (*f<sup>τ</sup>* Hz = 0.4 Hz; 0.15 Hz; 0.1 Hz; 0.04 Hz), computed for the four experimental conditions using eVAR (first row) and eVARFI (second row) identification methods. Statistically significant differences between pairs of conditions are marked with an asterisk.

Figure 4 presents the same information as the previous figure but considering the HP only as target process and computing the corresponding univariate (Figure 4a), bivariate (Figure 4b with SAP and Figure 4c with RESP) and multivariate (Figure 4d) complexity measures. We find that, moving from SU to UP, the univariate measure *CX*HP|*X*HP evaluated at a time scale corresponding to the cutoff frequency of 0.4 Hz decreases significantly for both eVAR and eVARFI; this decreased complexity of HP in the UP condition is observed identically when the dynamics of SAP, of RESP, and of both SAP and RESP are included in the conditional entropy measure. Moreover, the inclusion of one or both the other time series makes such a decrease statistically significant also with a cutoff of 0.15 Hz when eVAR analysis is performed. Moving from SU2 to MA, we observe that only the eVARFI analysis detects statistically significant differences. Specifically, using eVARFI the univariate complexity of HP decreases during MA at time scales compatible with the cutoff frequency of 0.4 Hz, and also with cutoff 0.15 Hz when SAP dynamics (but not RESP dynamics) are considered. For both protocols and using both estimation methods, the analysis at longer time scales (cutoff 0.1 Hz and 0.04 Hz) does not lead to any significant variation in any of the HP complexity measures.

**Figure 4.** Distribution of the univariate (**a**), bivariate (**b**) with systolic arterial pressure (SAP) and (**c**) with respiration (RESP), and multivariate (**d**) complexity measures depicted as boxplot (mean and confidence intervals, yellow filled box; standard deviation, black vertical line) and original values (dots) for selected frequencies (*f<sup>τ</sup>* Hz = 0.4 Hz; 0.15 Hz; 0.1 Hz; 0.04 Hz), computed for the four experimental conditions using eVAR (first row) and eVARFI (second row) identification methods. Statistically significant differences between pairs of conditions are marked with an asterisk.

Figure 5 depicts the distribution across subjects of the differencing parameter **d** = [*d*HP, *d*SAP, *d*RESP] computed using the Whittle semiparametric estimator for the three time series and for each condition. The differencing parameter is positive for HP and SAP, while it is lower and centered around zero for RESP. The statistical significance of the variations of this parameter moving from SU1 to UP, or from SU2 to MA, was tested with the same linear mixed-effects model used for the complexity measures. Results indicate that there are statistically significant differences in the values of *d* only for SAP, documenting that conditions of stress evoked by UP or MA determine an increase of the differencing parameter estimated for this time series.

**Figure 5.** Distribution of the long-rang parameter *di* for each of the time series considered (heart period (HP), SAP and RESP) and for the four conditions.

#### **5. Discussion**

Building on recent developments regarding linear multiscale time series analysis [7,8,19,20], the present study introduces a novel approach to assess the dynamical complexity of multivariate processes accounting for the presence of short term dynamics and long-range correlations. The adopted linear parametric framework retains the advantage of previous formulations [7,19,20] related to the computational reliability of complexity estimated even at coarse time scales and over relatively short time series (few hundred points), and is here integrated with the incorporation of long-range dynamics (fundamental to a proper evaluation of complexity at coarse time scales [8]) and extended for the first time towards a fully multivariate implementation.

The usefulness of the proposed multivariate and multiscale measure of complexity was demonstrated in practice by the analysis of cardiovascular and cardiorespiratory interactions. In particular, the multivariate measure of complexity allowed us to evaluate the overall impact that the physiological mechanisms related to postural and mental stress have on the joint dynamics of HP, SAP and RESP. Simultaneously, the multiscale assessment of complexity led to elicit the contribution of mechanisms operating across multiple time scales (evaluation at fine time scales, with cutoff frequency of 0.4 Hz encompassing VLF, LF and HF oscillations) and that of mechanisms confined to slower oscillations (evaluation at coarser time scales, with cutoff at 0.15 Hz excluding HF oscillations and cutoffs at 0.1 Hz and 0.04 Hz excluding progressively also LF oscillations). In addition the evaluation of the complexity of HP, either considering its own dynamics only or the dynamics of SAP and/or RESP, was exploited to relate the overall complexity trends reflected in the multivariate measure to variations specifically related to physiologically well-known mechanisms of cardiovascular and cardiorespiratory interactions. To aid interpretation of the results and the discussion of the related physiological mechanisms, we report in Table 1 the statistically significant increase or decrease in complexity observed without (eVAR method) or with (eVARFI method) modeling long-range correlations during head-up tilt (comparison SU1 vs. UP) or during mental arithmetics (comparison SU2 vs. MA) at the time scales encompassing the whole spectral content (0.4 Hz), filtering the HF oscillations (0.15 Hz), and filtering in part (0.1 Hz) or completely (0.04 Hz) the LF oscillations. Complexity is linearly dependent on tilt inclination: a steeper inclination of tilt table produces higher degrees of sympathetic activation, with corresponding lower values of entropy-based complexity measures [40]. Studying the trends of the multivariate complexity measure, its decrease with tilt (Figure 3) documents a simplification of the overall dynamics of HP, SAP and RESP. This result can be mainly driven by the less complex dynamics of heart rate variability in the upright position, that we document calculating the complexity of HP (Figure 4a). Indeed it is well known that the dynamics of HP tend to become less complex as a consequence

of the rise of LF oscillations and the weakening of HF oscillations related to the shift of the sympathovagal balance towards sympathetic activation and vagal deactivation during head-up tilt [15,40]. The fact that the decrease of multivariate complexity is not statistically significant for cutoff 0.15 Hz and lower suggests that it is related mainly to the weakening of HF oscillations; indeed, HF oscillations of both HP and SAP are known to be blunted with tilt [15,41,42]. In addition, the finding that the decrease in complexity is observed in the same way for VAR and VARFI analysis suggests that the impact of long-range correlations does not change substantially from rest to tilt. Overall, these results suggest that cardiovascular and respiratory oscillations in the LF and VLF band considered together do not alter significantly their complexity in the transition from the supine to the upright body position.


**Table 1.** Significant differences (*p-value* < 0.05) between pair of conditions for each measure and frequency. The arrows indicate if the measure increases or decreases from rest to stress.

On the other hand, the increase of the global complexity with MA observed for eVAR modeling across multiple time scales (cutoffs 0.4, 0.15, 0.1) is in agreement with previous findings relevant to the HP time series [7]. The fact that such increase is not observed for eVARFI indicates that long range correlations have a different impact on the cardiovascular and respiratory dynamics during MA compared with the relaxed condition SU2. In particular, eVARFI estimation reports values of the multivariate complexity which are unchanged at low time scales (cutoff 0.4 Hz), tend to decrease at cutoff 0.15 Hz, and decrease significantly at cutoffs 0.1 Hz and 0.04 Hz. We ascribe this lower complexity to a larger contribution of long-range correlations acting in LF and especially VLF bands, which is supported by the fact that the differencing parameter increases significantly, for HP and especially for SAP, during MA. This confirms the regularizing role of long range correlations on physiological dynamics [8,10].

The complexity of HP assessed at lower time scales (cutoff 0.4 Hz) always decreases in the UP position, for all measures (univariate, both bivariate, and multivariate), as shown in Figure 4. This documents the well-known simplification of heart rate variability induced by head-up tilt, which is known to evoke sympathetic activation and vagal withdrawal making the cardiac dynamics more regular [31,40,43,44]. The fact that this result is observed identically for the eVAR and eVARFI approaches indicates that long-range correlations do not impact significantly the evaluation of complexity performed at short time scales. On the other hand, focusing on intermediate time scales for which HF oscillations are cut off (*f<sup>τ</sup>* = 0.15 Hz), we find that the decrease in complexity is lost when considering the HP dynamics individually, but is maintained when SAP and RESP, taken individually or together, are used to reduce the complexity of HP; moreover, this holds for eVAR but not for eVARFI. Taken together, these results suggest that slow oscillations of SAP and/or RESP reduce the complexity of HP in the transition from rest to tilt, and such regularizing action is related to long-range correlations. This finding is compatible with the known increase of cardiovascular interactions during tilt, documented in previous studies using a number of causality measures including information-theoretic ones [13,24,26,27], and also with the existence of slowly varying

respiratory patterns that enhance their effect on the variability of HP during postural stress [24]. Here, the presence of cardiorespiratory and cardiovascular interactions is documented indirectly observing the lower values during tilt of the bivariate and multivariate complexity measures compared with the univariate one, through an approach similar to that proposed in Ref. [14].

The complexity of HP tends to decrease also with the mental stress induced by MA. Contrary to the postural stress induced by HUT, the decrease is observed only using the VARFI approach, and not using the VAR. The absence of significant changes in complexity using methods that do not model long-range correlations is in agreement with previous findings [7,14]. Our results document that the decrease in complexity is due to the different impact of long-range correlations, detected only modeling them through the VARFI approach, and again are confirmed by the higher values of the differencing parameter estimated for HP and for SAP (but indeed not for RESP) during MA compared with the relaxed condition (Figure 5). A differencing parameter *d* = 0 indicates that there are no long-range correlations. The rather small values of *d* observed for the RESP time series likely reflect the fact that, even if not controlled, the respiratory activity is confined in a narrow band of the frequency spectrum and thus it does not exhibit the slow trends typical of long-range correlated dynamics. Moreover, the changes of the differencing parameter *d* observed in the upright position for the systolic pressure (and only to a smaller extend for the heart period) document that the sympathetic activation produced by the tilt table inclination may modulate the impact that long range correlations have on the cardiovascular variability. We note also that, when SAP is considered, the significant changes extend to the second cutoff frequency (0.15 Hz), indicating the increasing contribution of SAP long-memory dynamics in reducing the complexity of HP. Therefore, we state the importance of accounting for long-range correlations in the assessment of the changes in the complexity of heart rate variability induced by mental stress. This may have relevance for the practical applications focused on the detection of mental workload or stress [45–48]).

Future extensions of the present work can be targeted first at investigating methodological aspects, such as the possibility of describing interactions between processes at the level of long-range correlations (possibly with the modeling of cointegration [49]) or the extension of the framework to the direct evaluation of Granger-causal interactions [20] in a multivariate context where long-range correlations are modeled [8]. Regarding applicative contexts, the analysis of multivariate multiscale dynamics is of particular interest in econometrics [50] and in the neurosciences [51,52], where dynamics spanning several temporal scales are commonly observed and multichannel data acquisition is ubiquitous.

**Author Contributions:** Conceptualization, A.P.R. and L.F.; methodology, L.F., R.P., A.P.R., A.M. and M.E.S.; software, L.F., R.P., A.P.R., A.M. and C.A.; validation, R.P., A.P.R., A.M. and M.E.S.; formal analysis, A.P.R., A.M. and M.E.S.; investigation, M.J.; resources, L.F. and M.J.; data curation, A.M. and M.J.; writing–original draft preparation, A.M.; writing—review and editing, L.F, R.P., A.P.R., C.A., A.M. and M.E.S.; visualization, A.P.R., A.M., M.E.S. and R.P.; supervision, L.F.; project administration, L.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partially supported by UIDB/00144/2020 (CMUP) and UIDB/04106/2020 (CIDMA), which are funded by FCT with national (MEC) and European structural funds through the program FEDER, under PT2020, and project STRIDE-NORTE-01-0145-FEDER-000033 funded by ERDF-NORTE 2020), PRIN 2017 PRJ-0167, "Stochastic forecasting in complex systems", and PON R&I 2014-2020 AIM project no. AIM1851228-2, University of Palermo, Italy.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

Since the computation of complexity works for finite-order VAR processes, while the VARFI process is equivalent to a VAR of infinite order, it is necessary to devise an approximation procedure. To do this, we truncate the fractional integration part of the VARFI process (5) at a finite lag *q*:

$$\begin{aligned} \text{diag}(\nabla^{\mathbf{d}}) \approx \mathbf{G}(L) &= \text{diag}\left[\sum\_{k=0}^{q} \mathbf{G}\_{k}^{(1)} L^{k} \quad \dots \quad \sum\_{k=0}^{q} \mathbf{G}\_{k}^{(M)} L^{k}\right] = \\ &\sum\_{k=0}^{q} \text{diag}\left[\mathbf{G}\_{k}^{(1)}, \dots, \mathbf{G}\_{k}^{(M)}\right] L^{k} = \sum\_{k=0}^{q} \mathbf{G}\_{k} L^{k}, \end{aligned} \tag{A1}$$

so that the VARFI process of order *p* and with differencing parameter **d** can be rewritten as a VAR process of order *m* = *p* + *q* with parameters given by:

$$\mathbf{B}(L) = \mathbf{A}(L)\mathbf{G}(L) = \left(\mathbf{I}\_M - \sum\_{i=1}^p \mathbf{A}\_i L^i\right) \left(\sum\_{k=0}^q \mathbf{G}\_k L^k\right). \tag{A2}$$

The coefficients of the VAR polynomial **<sup>B</sup>**(*L*) = **<sup>I</sup>***<sup>M</sup>* <sup>−</sup> *<sup>p</sup>*+*<sup>q</sup>* ∑ *k*=0 **B***kL<sup>k</sup>* result from the multiplication of the two polynomials in (A2), which in the case *q* ≥ *p* yields:

$$\begin{aligned} \mathbf{B}\_0 &= \mathbf{I}\_M \\ \mathbf{B}\_k &= \begin{cases} -\mathbf{G}\_k + \sum\_{i=1}^k \mathbf{A}\_i \mathbf{G}\_{k-i\prime} & k = 1, \dots, p \\ -\mathbf{G}\_k + \sum\_{i=1}^p \mathbf{A}\_i \mathbf{G}\_{k-i\prime} & k = p+1, \dots, q \\ \sum\_{i=0}^p \mathbf{A}\_{i+k-q} \mathbf{G}\_{q-i\prime} & k = q+1, \dots, q+p \end{cases} &\text{(A3)} \end{aligned} \tag{A3}$$

Then, the Equation (A3) are exploited to identify the finite-order VAR process, finding its coefficients **B***<sup>k</sup>* from the differencing parameters *di*, which inserted in (6) yield the parameters **G***k*, and from the VAR parameters **A***k*.

#### **Appendix B**

The procedure that represents how the structure of a VAR process is altered when the process is represented at an assigned scale *τ* is based on the two steps of filtering the VAR process and of downsampling the filtered process [19,20]. The first step transforms the VAR process into a VARMA process (Equation (10)), which is equivalent of a SS process (eq. 11), for which the parameters **K**(*r*) , **C**(*r*) and **B**(*r*) are given by [19]:

$$\begin{aligned} \mathbf{C}^{(r)} &= \begin{bmatrix} \mathbf{B}\_1 & \cdots & \mathbf{B}\_m & \mathbf{D}\_1 & \cdots & \mathbf{D}\_r \end{bmatrix} \\ \mathbf{K}^{(r)} &= \begin{bmatrix} \mathbf{I}\_M & \mathbf{0}\_{M \times M(m-1)} & \mathbf{D}\_0^{-T} & \mathbf{0}\_{M \times M(r-1)} \end{bmatrix} \\ \mathbf{B}^{(r)} &= \begin{bmatrix} \mathbf{C}^{(r)} \\ \mathbf{I}\_{M(m-1)} & \mathbf{0}\_{M(m-1) \times M(r+1)} \\ \mathbf{0}\_{M \times M(m+r)} & \mathbf{I}\_{M(r-1)} & \mathbf{0}\_{M(r-1) \times M} \end{bmatrix} \end{aligned} \tag{A4}$$

The second step exploits theoretical findings showing that the downsampled version of an SS process has itself an SS representation [22,23]. Here, downsampling the SS process (11) with a factor *τ* yields the process **<sup>X</sup>**(*τ*) *<sup>n</sup>* <sup>=</sup> **<sup>X</sup>**(*r*) *<sup>n</sup><sup>τ</sup>* , which has the SS representation:

$$\mathbf{Y}\_{n+1} = \mathbf{B}^{(\tau)} \mathbf{Y}\_n + \mathbf{W}\_n \tag{A5a}$$

$$\mathbf{X}\_{\mathfrak{n}}^{(\tau)} = \mathbf{C}^{(\tau)} \mathbf{Y}\_{\mathfrak{n}} + \mathbf{V}\_{\mathfrak{n}\_{\tau}} \tag{A5b}$$

where **V***<sup>n</sup>* and **W***<sup>n</sup>* are different white noise processes with variances Σ**<sup>W</sup>** and Σ**<sup>V</sup>** and covariance Σ**VW**, respectively serving as innovations for the downsampled process **<sup>X</sup>**(*τ*) *<sup>n</sup>* and for the state process **<sup>Y</sup>***n*. Thus, the process (A5) is an SS(**B**(*τ*) , **C**(*τ*) , Σ**W**, Σ**V**, Σ**VW**) process whose parameters can be obtained as [22,23]:

$$\begin{aligned} \mathbf{B}^{(r)} &= \left(\mathbf{B}^{(r)}\right)^{\mathsf{T}} \\ \mathbf{C}^{(r)} &= \mathbf{C}^{(r)} \\ \boldsymbol{\Sigma}\_{\mathsf{V}} &= \boldsymbol{\Sigma}\_{\mathsf{E}^{(r)}} \\ \boldsymbol{\Sigma}\_{\mathsf{V}\mathsf{W}} &= \left(\mathbf{B}^{(r)}\right)^{\mathsf{T}-1} \mathbf{K}^{(r)} \boldsymbol{\Sigma}\_{\mathsf{E}^{(r)}} \\ \boldsymbol{\Sigma}\_{\mathsf{W}}(1) &= \boldsymbol{\Sigma}\_{\mathsf{E}^{(r)}} \left(\mathbf{K}^{(r)}\right)^{\mathsf{T}} \mathbf{K}^{(r)}, \quad \tau = 1 \\ \boldsymbol{\Sigma}\_{\mathsf{W}}(\tau) &= \mathbf{B}^{(r)} \boldsymbol{\Sigma}\_{\mathsf{W}}(\tau - 1) \left(\mathbf{B}^{(r)}\right)^{\mathsf{T}} \\ &+ \boldsymbol{\Sigma}\_{\mathsf{E}^{(r)}} \left(\mathbf{K}^{(r)}\right)^{\mathsf{T}} \mathbf{K}^{(r)}, \quad \tau \ge = 2 \end{aligned} \tag{A6}$$

Then, the SS(**B**(*τ*) , **C**(*τ*) , Σ**W**, Σ**V**, Σ**VW**) process can be converted in the form of Equation (12), which evidences the innovations and has parameters (**B**(*τ*) , **C**(*τ*) , **K**(*τ*) , Σ**E**(*τ*)). To move to this representation from (A5) it is necessary to solve the discrete algebraic Ricatti equation [22,23]:

.

$$\begin{aligned} \mathbf{P} &= \mathbf{B}^{(\tau)} \mathbf{P} (\mathbf{B}^{(\tau)})^T + \boldsymbol{\Sigma}\_{\mathbf{W}} - (\mathbf{B}^{(\tau)} \mathbf{P} \mathbf{C}^{(\tau)} + \boldsymbol{\Sigma}\_{\mathbf{V} \mathbf{W}}) \cdot \\ &\cdot (\mathbf{C}^{(\tau)} \mathbf{P} (\mathbf{C}^{(\tau)})^T + \boldsymbol{\Sigma} \mathbf{v})^{-1} (\mathbf{C}^{(\tau)} \mathbf{P} (\mathbf{B}^{(\tau)})^T + \\ &+ (\boldsymbol{\Sigma}\_{\mathbf{V} \mathbf{W}})^T), \end{aligned} \tag{A7}$$

which leads to the derivation of the two unknown parameters of the downsampled process:

$$
\Sigma\_{\mathbf{E}^{(\tau)}} = \mathbf{C}^{(\tau)} \mathbf{P} (\mathbf{C}^{(\tau)})^T + \Sigma\_\mathbf{V} \tag{A8a}
$$

$$\mathbf{K}^{(\tau)} = \frac{\mathbf{B}^{(\tau)} \mathbf{P} (\mathbf{C}^{(\tau)})^T + \Sigma\_{\mathbf{VW}}}{\Sigma \mathbf{v}}.\tag{A8b}$$

Finally, the covariance of the downsampled process can be computed from the process parameters by solving the following discrete-time Lyapunov equation:

$$\boldsymbol{\Omega} = \mathbf{B}^{(\tau)} \boldsymbol{\Omega} (\mathbf{B}^{(\tau)})^T + \boldsymbol{\Sigma}\_{\mathbf{E}^{(\tau)}} (\mathbf{K}^{(\tau)})^T \mathbf{K}^{(\tau)} \tag{A9a}$$

$$
\Sigma\_{\mathbf{X}^{(\tau)}} = \mathbf{C}^{(\tau)} \mathbf{O} (\mathbf{C}^{(\tau)})^T + \Sigma\_{\mathbf{E}^{(\tau)}}.\tag{A9b}
$$

The two covariance matrices Σ**X**(*τ*) and Σ**X**(*τ*) are used in Equation (14) to derive the multivariate complexity of the process **X** observed at scale *τ*.

#### **Appendix C**

In order to derive the partial variances used in Equation (4a) and Equation (4b) for the computation of the complexity of scalar processes, SS submodels need to be formed in a way such that the observation equation (Equation (14)) contains only some of the scalar processes. It is important to note that these submodels are not in innovations form, but are rather SS models as in (A5) with parameters (**B**, **C**(*a*) , **K**Σ**XK***<sup>T</sup>*, Σ**X**(*a*, *a*), **K**Σ**X**(:, *a*)) [19]. Such SS models can be converted to the SS form as in (12), with innovation covariance <sup>Σ</sup>**E**(*τ*) *<sup>a</sup>* , solving the discrete algebraic Ricatti Equations (A7) and (A8), so that the partial variance <sup>Σ</sup>*E*(*τ*) *j*|*a* is derived as the diagonal element of <sup>Σ</sup>**E**(*τ*) *<sup>a</sup>* corresponding to the position of the target *Xj*,*n*.

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Application of Time-Scale Decomposition of Entropy for Eye Movement Analysis**

#### **Katarzyna Harezlak \* and Pawel Kasprowski**

Silesian University of Technology, Akademicka 16, 44-100 Gliwice, Poland; pawel.kasprowski@polsl.pl

**\*** Correspondence: katarzyna.harezlak@polsl.pl; Tel.: +48-32-237-1339

Received: 23 December 2019; Accepted: 30 January 2020; Published: 1 February 2020

**Abstract:** The methods for nonlinear time series analysis were used in the presented research to reveal eye movement signal characteristics. Three measures were used: approximate entropy, fuzzy entropy, and the Largest Lyapunov Exponent, for which the multilevel maps (MMs), being their time-scale decomposition, were defined. To check whether the estimated characteristics might be useful in eye movement events detection, these structures were applied in the classification process conducted with the usage of the kNN method. The elements of three MMs were used to define feature vectors for this process. They consisted of differently combined MM segments, belonging either to one or several selected levels, as well as included values either of one or all the analysed measures. Such a classification produced an improvement in the accuracy for saccadic latency and saccade, when compared with the previously conducted studies using eye movement dynamics.

**Keywords:** eye movement events detection; nonlinear analysis time series analysis; approximate entropy; fuzzy entropy; multilevel entropy map; time-scale decomposition

#### **1. Introduction**

Biological signals representing the electrical, chemical, and mechanical activities that occur during biological events attracted the attention of many researchers. These interests are aimed at discovering patterns which may prove useful in understanding the underlying physiological mechanisms or systems. The range of biological signals explored in various studies include, inter alia: electroencephalogram (EEG), electrocardiogram (ECG), surface electromyogram (sEMG), galvanic skin response (GSR), and arterial blood pressure (ABP). Obtaining these characteristics is of great importance in medicine, as this it may enable the differentiation of typical and atypical behaviors, supporting in this way diagnosing and treatment. However, analysis of some of them may also be applied, for example, in cognitive or psychological studies.

Acquiring bio-signal recordings requires various biomedical instruments to be applied. For example, eye movement signals, which are of interest in this article, can be measured by means of the electrical potential between electrodes placed at points close to the eye (electro-oculography (EOG)) or with the usage of less intrusive video-oculography (VOD), using video cameras and image processing algorithms for eye movement tracking.

Biological signals are naturally analogous; however, during measurement process conducted with a specific sampling rate, they are converted to a discrete-time form and constitute a biological time series. The distinct difficulty in the analysis of bio-signals is their nonlinear nature, therefore in order to extract useful features and components of the recorded signal, it is important to use appropriate processing

methods. Frequently, methods of nonlinear time series analysis, which can be deployed for many types of bio-signals, are chosen for this purpose. Some of these are described in the following section.

#### *1.1. Methods Used in Biological Signal Analysis*

The approaches used for quantifying biological signal characteristics use fractal and dynamic feature analysis. In the former group of traits, fractal power spectra and the fractal dimension are investigated, usually, through detrended fluctuation analysis (DFA) [1]. This approach enables the assessment of the existence of self-similarity and long-range correlations. The main goal of these kinds of studies is revealing the long-term memory of an explored system and predicting the system's future states. This method was successfully applied to such biological processes as ECG, EEG, EMG [2–5].

Nonlinear system dynamics may also be represented by various entropy measures [6], and the Largest Lyapunov Exponent (LLE) [7]. The entropy-based algorithms ascertain the degree of disorder and uncertainty in the underlying system described by an observed variable. Two primary purposes of these studies can be mentioned. The first of them is the discovery of the dynamic characteristics of biological signals of healthy people to use them as reference data in revealing anomalous cases. The second one is to determine events occurring within these signals.

For example, Chen et al. in [8] provided a review of entropy measures used in cardiovascular diseases. Furthermore, *approximate entropy* (ApEn) was used in [9] to obtain a better understanding of abnormal dynamics in the brain in the case of Alzheimer's disease (AD). Yents et al. [10] used spatio-temporal gait data from young and elderly subjects to investigate the performance of ApEn and *sample entropy* (SampEn). Cao et al. [11] presented the *fuzzy entropy* (FuzEn) metrics comparison and pointed the FuzEn out as a better tool than ApEn and SampEn, for distinguishing EEG signals of people with AD from those obtained for non-AD. Additionally, the comparison of various fuzzy entropy measures was conducted in the work [12]. Some synthetic datasets, together with EEG and ECG clinical datasets, as well as the gait maturation database, were used for comparison purposes.

EEG, ECG, and gait characteristics were also explored with the usage of the LLE. The works conducted in [13,14] indicated that this measure can be used to identify sleep stages, whereas in [15] the LLE was applied for the analysis of nonlinear changes of neural dynamics in the processing of stressful anxiety-related memories. Furthermore, automated diagnosis of electrocardiographic changes was successfully explored in [16] by using the Lyapunov Exponents for the detection of four types of ECG beats. Additionally, some works were devoted to investigating the application of the LLE in human locomotion for gait stability assessment. For example, Josinski et al., in their research [17] focused on finding LLE-based, easy-to-measure, objective biomarkers that could classify PD (Parkinson's disease) patients in early (preclinical) stages of the disease.

Nonlinear analysis was also conducted in regard to eye movement signal. However, its dynamics in this scope was not so widely explored as in other biological signals. Nonetheless, in several studies, it was discovered, based on the LLE [18–20] and entropy measures: ApEn and SampEn ([21,22]), that various phases of eye movement signal reveal different dynamic characteristics.

#### *1.2. Dynamic Eye Movement Signal Characteristics*

Eye movements are triggered by brain activity, which is a response to visual stimuli or an intention to gain knowledge of the surrounding world. The biological signal evoked this way consists of two main components: fixations, occurring when eyes are focused on part of a scene, and saccades taking place when eyes move to another scene area. During both events, eyes are in constant motion, even during fixations, generating micro-movements such as tremors, microsaccades, and drifts [23,24]. Additionally, in the reaction to the observed stimulus movement, the brain needs some time to commit a saccade. This

occurrence called *saccadic latency* [25] is shown, together with the preceding fixation and the following saccade, in Figure 1.

**Figure 1.** An example of eye movement events. At the beginning the eye is focused on the stimulus. When the stimulus moves to another location, the eye moves as well, yet after a period of saccadic latency [22].

Fixation and saccade differ in their characteristics. The first element is described by slow motions, small spatial dispersion, and relatively long duration (200–300 ms up to several sec.), while during the second one eyes move very quickly, reaching up to even 500◦/s, with a short duration (30–80 ms) and a higher amplitude of movement [26]. These differences are of prominent importance when event detection algorithms are considered. The most common approaches for extracting fixations and saccades from registered signal use dispersion and velocity thresholds [27]. One of the major drawbacks of these solutions is the dependency of the obtained results on user-defined thresholds and the lack of their commonly accepted values. The alternative approach presented in [28] assumes the use of machine learning and 14-features vector constructed based on the 100–200-ms surroundings of each sample. This vector consists of features describing the data in terms of dispersion, velocity, sampling frequency, and precision. The results of these studies showed that the machine learning technique is a promising solution. It was taken into consideration when a new approach was being developed using the findings of the previously described nonlinear dynamic systems analysis. For eye movement events detection, in [22], a multilevel map of ApEn was defined, and entropy values calculated at each map level were used for the feature vectors construction. Determining particular eye movement segments was conducted using the kNN classifier for various values of the k parameter. The classification outcomes were promising as well, especially in the scope of saccadic period detection, confirming the usefulness of ApEn measure in such an application. However, because one measure is only a single index describing the behavior of a dynamic time series, further investigations are still required. The currently presented research is the answer to this need.

#### *1.3. Contribution*

These studies introduce an extended application of the multilevel map by adding two additional measures: FuzEn and the LLE. According to the authors' knowledge, the fuzzy entropy has not been applied in exploring eye movement dynamics so far. The idea of the application of the multilevel map, built based on FuzEn and the LLE, for eye movement time series analysis is also a novel approach, as well as combining three multilevel maps in order to define feature vectors for machine learning eye movement event detection.

#### **2. Materials and Methods**

#### *2.1. The Dataset Description*

The presented research was conducted with the usage of the same dataset, which was applied in the studies described in [22]. It was collected during the"jumping point" experiment, using 29 dark points (*Npp*), distributed over a white 1280 × 1024 (370 mm × 295 mm) screen. The points' layout was designed in such a way to ensure both covering a screen area evenly and to obtain varying lengths of saccadic movements. Although the point was displayed in each location for 3 sec, only the first *NEMr* = 1024 ms were taken into account during further analysis. Eye movement signals were acquired employing the head-mounted JAZZ-novo eye tracker [29], with a sampling rate equal to 1000 Hz. It uses Direct Infra Red Oculography (IROG), which is embedded in the Cyclop ODS sensor measuring the resultant rotations of the left and the right eye. During registration, eyes are illuminated with a low intensity infrared (IR) light. The difference between the amounts of IR reflected back from the eye surfaces, carries information pertaining to the eye position changes.

24 participants (*Np*) aged between 22 and 24, with normal vision, took part in the experiment, which consisted of two sessions, separated by three weeks (*Nspp*). As a result,*NEMs* participants' sessions were gathered, *Nps* = *Np* ∗ *Nspp*, each of which comprised subsequent eye positions for all points' locations. Such a dataset was divided into *NEMs* eye movement series (*NEMs* = *Np* ∗ *Nspp* ∗ *Npp*), which are later referred to as *EM series*. Each of them included one participant's eye positions registered between the moment of the appearance of the stimulus and the time of its position change. Due to some registration and processing problems four participant session were removed from further analysis. Thus, the *Nps* value decreased by 4, and *NEMs* by 4 × *Npp*. Finally, *NEMs*—the number of eye movement time series—was 1276 (44 *participants*\_*sessions* × 29 *points*)

Subsequently, by applying the standard procedure of the two-point signal differentiation, the first derivative of horizontal coordinates for each EM series was obtained. Thus, each series corresponds to the velocity of eye movement in a period between two appearances of the stimulus. This quantity was used to make research independent to the point positions.

#### *2.2. The Method*

The first step of the nonlinear time series analysis was the evaluation of the three chosen measures: ApEn, FuzEn and the LLE, for each EM series. These calculations were conducted taking into account the structure of the multilevel map (MM) introduced in [22]. It represents the concept of a time-scale decomposition of the particular measure, whose values are calculated for various segments of EM series following the schema presented in Figure 2.

**Figure 2.** Segments in the multilevel time-scale structure.

A signal at the (*l* + 1)-th level of the map is divided into two shorter segments at the *l*-th level. The smallest size of segments was used at the 1st level of the map, and for this research, it was set to *Nss*=64 ms. As mentioned, the time series length *NEMr* was set to 1024 ms, which resulted in *l* = 5 map levels. This means that at the 2nd level the length of a segment equaled 128 ms, at the 3rd: 256 ms, at the 4th: 512 ms and finally at the 5th: 1024 ms. For each map cell, all three measures were estimated independently; thus, three MMs were obtained for each EM series (Figure 4).

#### *2.3. Approximate Entropy*

The approximate entropy method was proposed in [30] as a measure of regularity for quantifying complexity within a time series. For a given integer parameter *m* and a positive real one *r*, for a time series *<sup>u</sup>*(1), ... , *<sup>u</sup>*(*N*) in **<sup>R</sup>**, the sequence of vectors *<sup>x</sup>*(1), ... , *<sup>x</sup>*(*<sup>N</sup>* − *<sup>m</sup>* + <sup>1</sup>) in **<sup>R</sup>***<sup>m</sup>* can be defined, where *x*(*i*)=[*u*(*i*),..., *u*(*i* + *m* − 1)].

For 1 ≤ *i* ≤ *N* − *m* + 1, using probability that any two vectors of length *m* are similar within a tolerance given by *r*:

$$\mathbb{C}\_{i}^{m}(r) = \frac{\text{number of } j \le N - m + 1 : d[\mathbf{x}(i), \mathbf{x}(j)] \le r}{N - m + 1} \tag{1}$$

where the distance *d* is defined as follows:

$$d[\mathbf{x}(i), \mathbf{x}(j)] = \max\_{k \in \{1, m\}} \left( |u(i+k-1) - u(j+k-1)| \right) \tag{2}$$

and the average logarithmic probability over all m-patterns:

$$\Phi^m(r) = \left(N - m + 1\right)^{-1} \sum\_{i=1}^{N-m+1} \log \mathbb{C}\_i^m(r) \tag{3}$$

approximate entropy can be obtained:

$$\text{ApEn}(m, r) = \lim\_{N \to \infty} [\Phi^m(r) - \Phi^{m+1}(r)]. \tag{4}$$

The asymptotic formula in Equation (4) cannot be used directly, and the estimator of approximate entropy (for sufficiently large values of *N*) is defined as:

$$\text{ApEn}(m, r, N) = \Phi^m(r) - \Phi^{m+1}(r). \tag{5}$$

There is usually significant variation in *ApEn*(*m*,*r*, *N*) over the range of *m* and *r*, when a particular system is taken into consideration [31]. Thus, based on the analysis conducted in [22], in this research *m* was set to 2 and *r* as 0.2 × *SD* (the standard deviation of the entire time series) in order to define the first MM.

#### *2.4. Fuzzy Entropy*

The second MM was prepared based on the fuzzy entropy values calculated according to the method introduced in [32], which consists of several steps.

For a given time series *u*(1), ... , *u*(*N*) and integer parameter *m*, the sequence of vectors *X<sup>m</sup>* <sup>1</sup> ,..., *<sup>x</sup><sup>m</sup>* (*N*−*m*+1) is defined as:

$$X\_i^{\mathfrak{m}} = u(i), \dots, u(i+m-1) - u0(i) \tag{6}$$

where *X<sup>m</sup> <sup>i</sup>* represents *m* consecutive *u* values, commencing with the *i* − *th* point and generalized by removing a baseline:

$$
u0(i) = \frac{1}{m} \sum\_{j=0}^{m-1} u(i+j)\tag{7}$$

Subsequently, the distance between two vectors is defined as:

$$d\_{ij}^{m} = d[X\_i^{m}, X\_j^{m}] = \max\_{k \in \{0, m-1\}} |(u(i+k) - u0(i)) - (u(j+k) - u0(j))|\tag{8}$$

and given *n* and *r* the similarity degree *D<sup>m</sup> ij* of *<sup>X</sup><sup>m</sup> <sup>j</sup>* to *<sup>X</sup><sup>m</sup> <sup>i</sup>* through a fuzzy function *<sup>μ</sup>*(*d<sup>m</sup> ij* , *n*,*r*) is calculated:

$$D\_{ij}^{\rm m} = \mu(d\_{ij}^{\rm m}, n, r) \tag{9}$$

where:

$$\mu(d\_{ij}^{\rm m}, n, r) = \exp\left(-\frac{\left(d\_{ij}^{\rm m}\right)^{\rm n}}{r}\right) \tag{10}$$

Finally, the fuzzy entropy for finite sets is estimated as:

$$FuzEn(m, n, r, N) = \ln \phi^m(n, r) - \ln \phi^{m+1}(n, r) \tag{11}$$

where

$$\phi^{\rm m} = \frac{1}{N - m} \sum\_{i=1}^{N-m} \left( \frac{1}{N - m - 1} \sum\_{j=1, j \neq i}^{N-m} D\_{ij}^{\rm m} \right) \tag{12}$$

*N* corresponds to the size of the time series, which depends on the MM level under consideration. *m* is the length of sequences to be compared and was set to the same value as for evaluating approximate entropy. The parameters *n* and *r*, determining the width and the gradient of the boundary of the exponential function, respectively, were chosen based on experimental data, following suggestions provided by the method's authors. Fuzzy entropy was evaluated for 20 time series with different *r* and *n* values. Subsequently, standard deviation (SD) of the results was calculated. The parameter values for which a slow decrease in SD was observed were chosen for the purpose of these studies. They were 2 for *n* and 0.075 × *SD* (of the particular EM series) for *r*.

#### *2.5. The Largest Lyapunov Exponent*

The Largest Lyapunov Exponent is a measure which estimates the amount of chaos in dynamic systems. Chaos is observed when neighbouring paths followed by the system, starting from very close initial conditions, rapidly—exponentially fast—move to different states. The reconstruction of the system states is feasible based on measurements of a single observed property when Takens' embedding theorem is applied [33]. A time delay embedded series—with time lag denoted by *τ* and embedding dimension by *m*—can be defined as a transformation of the original time series *u*, such that:

$$\mathbf{x}\begin{pmatrix} i \end{pmatrix} = \left[ \boldsymbol{u}\begin{pmatrix} i \end{pmatrix}, \boldsymbol{u}\begin{pmatrix} i+\tau \end{pmatrix}, \cdots, \boldsymbol{u}\begin{pmatrix} i+(m-1)\ \tau \end{pmatrix} \right], \mathbf{i} = \{1, 2\cdot \cdots \cdot, M\} \tag{13}$$

where *M* = *N* − (*m* − 1) × *τ*.

A pair [*x*(*i*), *x*(*j*)] of nearest neighbours, starting close to one another in a chaotic system, diverges approximately at a rate given by the Largest Lyapunov Exponent *λ* [34]:

$$d\_j(i) \approx d\_{j0} e^{\lambda(i\Delta t)}\tag{14}$$

where *dj*(*i*) is the Euclidean distance after *i* time steps, Δ*t* is the sampling rate of the time series and *dj*<sup>0</sup> is the initial pair separation. Solving the Equation (14) by taking the logarithm of both sides, the Largest Lyapunov Exponent can be calculated as follows:

$$
\lambda \approx \frac{1}{i\Delta t} \ln(\frac{d\_j(i)}{d\_{j0}}) \tag{15}
$$

Equation (15) provides the way for evaluating *λ* for two specific neighbouring points over a specific interval of time. Thus, to approximate the Lyapunov Exponent for a whole dynamic system, it has to be averaged for all *j* = 1, 2, . . . , *M*, where *M* = *N* − (*m* − 1) × *τ*. Negative Lyapunov Exponent values indicate convergence, while positive demonstrate divergence and chaos. The parameters *m* and *τ* are calculated with the usage of the False Nearest Neighbours (FNN) [35] and Mutual Information Factor [36], respectively. The above-described methods were used in the current research to define the third MM.

#### *2.6. Detection of Eye Movement Events*

Following the studies presented in [22], the elements of three MMs were used to define feature vectors for the classification process. Various vector types were defined. They consisted of differently combined MM segments, belonging either to one or several selected levels, as well as included values either of one or all the analysed measures (Figure 3). However, before the creation of the vector, data from all MMs levels were rescaled using min-max normalization:

$$\mathbf{x}\_{\text{new}} = \frac{\mathbf{x} - \mathbf{x}\_{\text{min}}}{\mathbf{x}\_{\text{max}} - \mathbf{x}\_{\text{min}}} \tag{16}$$

Finally, the following sets were prepared (see Figure 3 for examples):

	- setX\_ApEn, setX\_FuzEn, setX\_LLE, including one feature
	- setX\_ApEn\_FuzEn\_LLE, including three features
	- setX\_Y\_ApEn, setX\_Y\_FuzEn, setX\_Y\_LLE, including two features
	- setX\_Y\_ApEn\_FuzEn\_LLE, including six features
	- setX\_Y\_Z\_ApEn, setX\_Y\_Z\_FuzEn, setX\_Y\_Z\_LLE, including three features
	- setX\_Y\_Z\_ApEn\_FuzEn\_LLE, including nine features

**Figure 3.** Sample one-, two-, and three-dimensional feature vectors embedded in the MM structure [22].

The number of classes (*Ncl*)—the number of segments at the lowest level in the given feature vector —differed in particular feature vectors, depending on the segment used ([22]):

$$N\_{cl} \left( \text{set} X\_{-} \text{Y\\_Z} \right) = \frac{N\_{EMr}}{N\_{ss} \times 2^{l\chi - 1}} \tag{17}$$

where *lX* is the map level and *X* = min(*X*,*Y*, *Z*). For example, if the set64\_128\_256 is considered: *X* = 64 and *lX* = 1. This means that for *lX* equal to:


These sets were separately used for feeding the *kNN* classifier, run with several values for the *k* parameter: *k* = {3, 7, 15, 31, 63, 127, 255}. Each feature vector set was divided into training and test sets, with the usage of the *leave-one-out cross-validation* approach: in each classifier run, *Npp* maps—calculated for one participant and for one session—were always left for defining a test set.

#### **3. Results**

The multilevel maps were estimated for each EM series and each analysed measure separately. Subsequently, they were averaged over all *NEMs* series, and finally, three resulting maps were created, as shown in Figure 4.


**Figure 4.** The multilevel maps obtained for the three measures: ApEN, FuzEn, the LLE, with values averaged for all EM series. The table at the bottom presents the maps' structure.

The maps cells were coloured according to their contents. Green indicates the lowest values, while red the highest ones. A repeating pattern can be noticed in each MM, in the 3rd and the 4th segments at the first level and the 2nd segment at the second one. These cells contain the lowest values calculated for each measure and relate to 128 ms of eye movement signal commencing with 128 ms. This dependency is easily noticeable in Figure 5, where the results obtained for the first two MMs levels are juxtaposed in the charts. However, to make this comparison feasible, the results had to be rescaled with the min-max normalization (Equation (16)).

**Figure 5.** The charts presenting values from the two levels of the ApEn, FuzEn and LLE maps.The first level at the top and the second at the bottom.

Statistical significance of the outcomes was checked by means of ANOVA and Tukey's Honest Significance Difference tests. The two result sets—for ApEn and the LLE measures—verified by the Shapiro–Wilk test, turned out to be normally distributed, yet for the third one—FuzEn—the null hypothesis was rejected. In this case, the Kruskal-Wallis test was additionally run. All differences in the results obtained for the 3rd and the 4th times scopes at the first MM level, compared to the remaining set of segments turned out to be significant. However, there was one exception for the LLE—the difference between the 3rd and 5th segments—and a few for FuzEn—the differences between the 3rd, 4th, and 5th time scopes—which occurred not to be significant. While analysing the results for the second level of the maps, it was revealed that all measures provided significant differences for the 2nd segment when collated with all other time scopes. On the contrary, the outcomes from the second part of the EM series—starting from the 6th (ApEn) and the 8th (FuzEn and the LLE) segments at the first MM level, as well as from the 5th segment at the second level—disclosed the lack of significant differences.

In the second stage of the tests, the usefulness of the three MMs in recognizing eye movement segments was verified with the usage of the kNN method. During the first classification runs, feature vectors consisting of one, two, or three elements calculated only for one measure, were used. Because the attention was focused on the highest possible signal resolution, only features based on segments 64 and 128, and their combinations with values from other levels were taken into account. As can be seen in Figure 6, the performance of the segment recognition is rather low, independently of *k* − *value* used,

especially when time scopes later than 256 ms are considered. Thus, in further analysis, only the first four segments were addressed.

**Figure 6.** The classification accuracy based on the first level of the ApEn, FuzEn and LLE multilevel maps, and for different the k–parameter values: k = 15 (top), k = 63 (middle), k = 255 (bottom).

The exploration of the presented charts revealed that a higher *k* parameter did not always provide better results. For this reason, the attention was focused on *k* − *value* equal to 63, for which the classification accuracy was similar to that, disclosed for k=255, and which ensures better performance of the execution.

The effect of the classification with the usage of feature vectors defined based on one measure, yet using several levels of the MMs is presented in Table 1. There are nine sets taken into account: three consisting of one feature: *set64\_[ApEN, FuzEN, LLE]*, three with two elements: *set64\_128\_[ApEN, FuzEN, LLE]*, and three including three components: *set64\_128\_256\_[ApEN, FuzEN, LLE].*


**Table 1.** The percentage of the correctly classified samples for various feature vectors defined with the usage of one measure. The smallest segment size = 64 ms.

The table content shows that the best results were obtained for the LLE measure in the first segment of eye movement signal, yet less than that got for k=255 (86% see Figure 6). However, its classification accuracy decreases significantly when other time scopes are considered. The opposite situation takes place for the two remaining groups of sets. As can be observed, both entropy measures at the beginning of the EM series revealed a low accuracy for vectors defined with only one element. The performance increases for later eye movement scopes when features are merged into two- or three-elements vectors (*set64\_128\_*, *set64\_128\_256\_*).

The individual measure outcomes can be further collated with those presented in Table 2, obtained for the mixed feature vectors.

**Table 2.** The percentage of the correctly classified samples for the feature vectors defined with the usage of three measures: ApEN\_FuzEn\_LLE. The smallest segment size = 64 ms.


Similarly, the best performance, slightly worse than previously, was achieved for the first 64 ms of the signals, classified with samples coming from the first level of the MMs. Subsequently, the accuracy decreases, especially for the *set64\_ApEN\_FuzEn\_LLE*. Thus, *set64\_128\_256\_ApEN\_FuzEn\_LLE* seems to be the proper choice for the classification purpose, despite lower than the best accuracy produced for the first segment. It provides the best or close to the best efficiency for all analysed segments.

In the next two Tables 3 and 4, outcomes for the subsequent segment size (128 ms) were shown. They, in terms of the result pattern, are similar to those discussed above. The best detection was delivered by the LLE measure in the first 128 ms. In the further time scope, the efficiency of this measure dropped, and ApEn and FuzEn disclosed better performance; the best for vectors with three features *(set64\_128\_256\_ApEn*, *set64\_128\_256\_FuzEn*). While comparing the classification accuracy to the previous results (for 64ms segment), it turned out to be generally higher. However, it should not be surprising, because of the less by half the number of classes in a feature set, which facilitates the sample recognition. Additionally, once again, the *set128\_256\_512\_ApEn\_FuzEn\_LLE* despite worse efficiency in the first time scope, provides, on average, the best results.


**Table 3.** The results–the percentage of the correctly classified samples for various feature vectors defined with the usage of one measure. The smallest segment size = 128 ms.

**Table 4.** The results–the percentage of the correctly classified samples for the feature vectors defined with the usage of three measures: ApEn\_FuzEn\_LLE. The smallest segment size = 128 ms.


The last step in the data exploration was the juxtaposition (Table 5) of the best results obtained in [22], using only ApEn (the third column), with the corresponding to them, achieved in these studies with the usage of three measures (the fourth column). The first column in the table represents the feature vector type, while in the second one, the eye movement segments were included. The last column in this table shows the best classification accuracy received in these studies for a feature set displayed in the first column and for the segment from the second one. All values presented in this column were estimated for the LLE measure in the first EM series scope.

In almost all cases, *Current accuracy* exposes better or similar performance; however, in one case (*set*64\_128), the worse accuracy was achieved.

**Table 5.** The best accuracy obtained in [22] and in the current research for the same feature vector types and the same EM series segment.


#### **4. Discussion**

Eye movement analysis has been conducted over many years, but research in this domain intensified over the last ten years. It resulted from the fact that acquiring knowledge concerning visual patterns has proved useful in many areas, because they provide essential data regarding people's behaviour, intentions, and interests. Proper reasoning in these fields requires an appropriate interpretation of collected recordings, which corresponds to correctly recognised eye movement events. Elements that constitute visual patterns

are fixations and saccades, and extracting them adequately from the eye movement signal is crucial. For many years, approaches based on spatial recordings dispersion and velocity have been used. However, as mentioned, they rely on user-defined thresholds, which is their disadvantage as most of the biological patterns vary between individuals (Figure 7). Therefore, some steps were taken to search for other algorithms, ensuring more adjusted event detection.

**Figure 7.** Example signal courses for eye movement velocity. For two more, refer to [22].

The studies presented here base their research on the dynamical aspects of eye movement signal. The main purpose of these investigations was to find characteristics of the particular signal parts, short enough to enable searching for the smallest eye movement events, saccades.

#### *4.1. Multilevel Maps Analysis*

The results presented in Figure 4 showed that measures calculated for one EM series, including three eye movement events—saccadic latency, saccade and fixation—have different values in various time scopes of the EM series. The previous studies [18,19,21,22] already provided some evidence in this regard for ApEn, SampEn, and the LLE; however, only during independent analyses. One of the advantages of the current studies is a comparison of these measures and introducing the additional signal description through the FuzEn usage.

When analysing the three multilevel maps from Figure 4, it may be noticed that two segments at the first level (128–192 ms and 192–256 ms) and one segment at the second level (128–256 ms) include the lowest values among all those calculated for the particular map. This means that all three measures detected, within the same period, an activity different than in the remaining EM series parts. This period is characterized by lower entropy values and a lack of chaotic behaviour. Collating MMs segments with signal recordings, for example, in Figure 7 shows that this period corresponds to the saccade. Further exploration of this figure unveils another distinction visible in the LLE map, in its first two segments at the first two levels. The positive values in this period show that neighboring points of EM series, which start close to one another, diverge with time approximately at a rate given by the LLE, as shown in Figure 8.

**Figure 8.** Example of divergent paths.

Such behaviour characterises chaotic dynamics, which, as it was revealed in [18,19], is related to saccadic latency. Fixations provide more homogeneous results than the saccadic latency and saccade, which is visible in the second part of the time series (512–1024 ms). The segments representing this component contain comparable values in the case of each of the studied metrics.

The above analysis demonstrate that given the presented MMs, it is feasible to indicate time series scopes within which three eye movement events take place.

#### *4.2. Classification Performance*

Promising visual inspection of the MMs found also a partial confirmation in the classification results, although overall accuracy was rather low. This outcome is a consequence of the fact that at least half of the samples—corresponding to a fixation—have very similar characteristics. It is especially visible in the charts presented in Figure 6. However, this is not the only reason, as for saccadic latency and fixation the similar entropy values were estimated (see Figure 4). Therefore, entropy metrics, at the lowest level of the MMs, do not introduce the distinction among these two events. The situation changes at the second level, where the fixation parts differ more noticeably from the first signal scopes.

Analysis of the results presented in Tables 1 and 3 shows that if saccadic latency is being searched for, the best approach is the usage of the LLE measure, for which accuracy at levels 76% and 84%, respectively, were achieved. Taking this event duration into account (approximately 120-200 ms), the segment of size equal to 128 ms appears to be sufficient. The advantage of the LLE over ApEn and FuzEn results from the fact that, while all three metrics quantifies irregularity in a time series, the first one is also able to describe how a system evolves in a particular period of its evolution. For saccadic latency segments, it indicates divergence, whereas during saccade, convergence. In the case of fixation, behaviour changes from convergent to chaotic and conversely. Therefore, this quantity is less efficient in fixation scopes, in which both entropy types yielded interchangeably better results. Yet, changing the measure during the classification process, depending on the segment under consideration seems not to be a proper solution. Thus, it is more reasonable to apply a feature vector consisting of all measures. Confirmation for such an approach provides the analysis of the third columns in Tables 2 and 4. The accuracy achieved for the combination of ApEn, FuzEn and the LLE was better than for any quantity except for the LLE. The decrease in its efficiency was approximately 3% for the longer segment (128 ms) and 10% for the shorter one (64 ms).

Generally, more satisfying classification outcomes were obtained for feature vectors with samples based on 128 ms-segments (Tables 3 and 4). As was previously noted, one of the reasons is the lower

number of classes in the training and testing sets. However, it may also stem from the fact that for longer segments, entropy values for different eye movement events are more distinguishable.

Finally, the comparison between findings presented in [22] and provided here indicates that the usage of more than one measure may advance the events detection efficiency (Table 5). It especially refers to the shorter segments at the beginning of the eye movement signal. Thus, by applying the presented approach, better classification has become possible for the higher resolution of signal segments, even when using fewer than ten features. The high efficacy of the saccadic latency detection facilitates recognizing the following saccade, while finding saccade period opens the possibilities for more detailed exploration in this scope.

However, the results of this method still need further improvement, which can be achieved by, for example, extending the feature vector and usage of other classifiers. Additionally, decreasing the number of fixation segments is worth considering. In this research, the studied fixation period was approximately 700 ms and could be shortened by 200–300 ms. Furthermore, within the explored eye movement data, there was a lack of movement called *smooth pursuit*, when the eyes follow a moving object. The dynamics of such a motion should also be explored.

The suggestions mentioned above are plans for future work. Moreover, because the proposed method was verified only for one dataset gathered with the usage of a particular eye tracker type, other signal recordings are intended to be analysed.

#### **5. Conclusions**

The methods for nonlinear time series analysis were used in the presented research in order to reveal eye movement signal characteristics. The primary motivation for these investigations was to find patterns, which can characterise the main components of eye movements: fixation, saccade and saccadic latency. Their proper recognition in registered eye movement signals supports studies in many areas using eye tracking technology. Fixations uncover the gazed areas of observed scenes, while saccades reveal changes in gaze positions. Sometimes, an eye movement is caused by the motion of an object belonging to the scene. Such a situation introduces an additional period in the signal called saccadic latency: the time needed for the brain to react to the motion of the observed stimulus. Discovering this period is the indication that the subsequent signal segment will represent a saccade in the direction of a new stimulus position. Recognizing this element is, in turn, equivalent to revealing the subsequent fixation. Although there are some commonly used methods for eye movement events detection, new solutions, which will provide more accurate event detection, are still being sought.

In the presented approach, three measures were used: approximate entropy, fuzzy entropy and the Largest Lyapunov Exponent, for which multilevel maps, being their time-scale decomposition, were defined. To check if the estimated characteristics may be useful in distinguishing eye movement components, these structures were applied in the classification process. It produced an improvement in the accuracy, for saccadic latency and saccade, when compared to previously conducted studies using eye movement dynamics. Nonetheless, further investigations towards more precise results are planned in the future.

**Author Contributions:** Conceptualization, K.H.; Data curation, P.K.; Formal analysis, K.H.; Methodology, K.H.; Software, K.H. and P.K.; Validation, K.H. and P.K.; Writing – original draft, K.H.; Writing – review & editing, K.H. and P.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Statutory Research funds of Faculty of Automatic Control, Electronics and Computer Science, Grant Number BK/2020.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**



c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
