The AR model in Eq. (1) can be represented in the Z-domain through its Z-transform yielding Y(z) = H(z)U(z), where H ( z ) = [ 1 − ∑ k = 1 p a k z − k ] − 1 = A ̄ ( z ) − 1 is the TF modelling the relationship between the input U(z) and the output Y(z). Applying the residue theorem, the latter can be expressed as (Baselli et al., 1997): H z = z p ∏ k = 1 p z − p k = ∏ k = 1 p H k z , (2) where p _k , k = 1, …, p, are the p poles of the AR process, i.e., the roots of A ̄ ( z ) , while the terms H ( k ) ( z ) = z z − p k ⋅ 1 / z * 1 / z * − p k are pole-specific factors associated each to a given pole p _k , with * the Hermitian transpose. Then, the power spectral density (PSD) of the process can be written in the Z-domain as P ( z ) = H ( z ) σ U 2 H * ( 1 z * ) , and expanded exploiting the Heaviside decomposition with simple fractions relevant to all the poles and weighted by the relevant residuals of P(z), to get (Baselli et al., 1997; Pernice et al., 2021): P z = ∑ q = 1 Q P q z = ∑ q = 1 Q r q p q z − p q − r q p q − 1 z − p q − 1 , (3) where r q = σ U 2 z ∏ h ≠ q ( z − p h ) ⋅ ∏ ( z − 1 − p h ) | z = p q are the residuals of P(z), q = 1, …, Q. Note that the number of components of the PSD is Q ≃ p/2 depending on the number of real poles; specifically, there is one component for each real pole and for each pair of complex conjugate poles. Then, by computing P(z) on the unit circle of the complex plane, P ( f ) = P ( z ) | z = e j 2 π f / f s , where f ∈ [0, f _s /2], with f _s the sampling frequency, it is possible to obtain the spectral profile of the process, P(f), as well as of the qth component, P ^(q)(f). Crucially, each spectral component P ^(q)(f) is described by a specific profile that is shaped by the corresponding TF factor | H ( q ) ( f ) | 2 = H ( q ) ( z ) H * ( q ) ( z ) | z = e j 2 π f / f s , and has an associated frequency (related to the pole frequency, f q = a r g { p q } 2 π ) and power (related to the pole residual, σ q 2 = r q for real poles and σ q 2 = 2 ⋅ R { r q } for complex conjugate poles). Note that the sum of the pole variances σ q 2 , with q = 1, …, Q, equals the total power of the process, which represents its variance σ Y 2 .

Importantly, in comparison to classical approaches based on integrating the PSD profile within the spectral bands of interest to get band-specific time domain powers (Krohova et al., 2020), the proposed method of spectral decomposition allows to focus only on the spectral components with frequencies within those bands, thus avoiding spurious contributions due to broadband oscillations. This peculiarity is exploited in the next subsection to define a pole-specific measure of self-predictability in the frequency domain.

In the previous section we have seen how analysing the spectrum of a physiologic process provides noteworthy information on the frequency-specific location of the oscillations of that process, thus allowing to identify and separate its different spectral components. Here, we show how the characterization of the spectral dynamics of a process can be also carried on by looking at its degree of linear self-predictability, drawing a connection with information theory.

The most popular information-theoretic measure of self-predictability is the information storage (IS), which quantifies, for a random process Y, the amount of information contained in the present state Y _n that can be predicted by the knowledge of its past states, Y n − (Lizier et al., 2012). In the linear signal processing framework, the IS has a straightforward formulation that involves the variance of the process and the variance of the prediction error of its AR representation (Barnett et al., 2009; Faes et al., 2015b): S Y = 1 2 ln σ Y 2 σ U 2 , (4) where σ Y 2 is the variance of Y and σ U 2 is the variance of the residual in Eq. (1), i.e., the partial variance of Y _n given its past Y n p . The quantity defined in Eq. (4) is a measure of LSP in the time domain. To expand it in the frequency domain, we exploit the spectral representation of the AR model in Eq. (1) and the spectral decomposition described in the previous and represented by Eqs. (2, 3). We start noting that the TF H(z) contains spectral information about the predictable dynamics of Y, as it is directly derived from the Z-domain representation of the AR model coefficients a _k , k = 1, …, p, which in turn describe these dynamics in the time domain. Then, we exploit the fact that such frequency-specific information can be particularized to each oscillatory component considering the TF factor |H ^(q)(f)|², so as to retrieve information on the pole-specific self-dynamics of the AR process. This factor is the squared TF associated to a real pole or the squared product of the two TFs associated to a pair of complex conjugated poles. We expect that stronger self-dynamics of Y at the frequency f _q are reflected by higher values of |H ^(q)(f)|², which indeed shows a positive peak at that frequency. Therefore, we define the frequency-specific spectral LSP measure as s Y q f = 1 2 ln σ Y 2 | H q f | 2 σ U 2 . (5)

The spectral in Eq. (5) can be written also as s Y ( q ) ( f ) = 1 2 ln σ Y 2 σ U 2 + 1 2 ln | H ( q ) ( f ) | 2 = S Y + s ̄ Y ( q ) ( f ) , and satisfies the spectral integration property (Geweke, 1982), i.e., it is such that its integral extended over all frequencies returns the time-domain LSP measure: S Y = 1 2 π ∫ − π π s Y q f d f . (6)

Note that Eq. (6) is satisfied because the full-frequency integral of the term s ̄ Y ( q ) ( f ) is null, i.e., ∫ − π π ⁡ ln | H ( q ) ( f ) | 2 d f = 0 (Rozanov 1967; Chicharro, 2011). Therefore, the spectral LSP consists of a frequency-independent part equal to S _Y and a frequency-specific part quantified by s ̄ Y ( q ) ( f ) , which takes negative values at some frequencies, depending on where the informative content is located.

The spectral decomposition of the pole-specific LSP into terms related to the Q oscillations of the AR process, depicted in Eq. (5), allows to locate the self-dynamics of Y in specific spectral bands with given frequency, as well as to compute their strength as the integral of these profiles within the investigated bands. Remarkably, the spectral LSP profiles display peaks as the PSD does, since both are derived from adaptations of the TF of the AR model describing the data. The difference between the two resides in the logarithmic formulation of the LSP in the well-known framework of information theory. Indeed, being a measure of information shared between the present and past states of the investigated processes according to its mathematical definition, it can be quantified in natural units (nats) in the time domain, and in nats/Hz in the spectral domain, thus acquiring a clear meaning in the context of information theory.

This section presents the use of surrogate and bootstrap data analyses to assess the statistical significance of the proposed measures of time domain and spectral LSP.

As far as we know, the task of validating the presence of a significant non-flat component within the spectrum of a process is not straightforward. In the spirit of surrogate data analysis, one should destroy the oscillation under scrutiny without altering the remaining spectral patterns of the process. However, this is undoubtedly challenging and requires in-depth analyses. Therefore, in this work we propose an empirical approach based on bootstrap data analysis (Politis, 2003) to assess the statistical significance of pole-specific measures of self-predictability, while time domain measures are validated through the well-known method of surrogate data analysis (Theiler et al., 1992). Validation is performed at the level of individual realizations of the observed process Y, obtained in the form of the time series y = {y(1), …, y(M)}, where M is the length of the time series. Specifically, a linear model as in Eq. (1) is first identified on the time series y through the vector least square approach; then, estimates of the spectral and the time domain LSP measures, denoted respectively as s y ( q ) ( f ) , with q = 1, …, Q, and S _y , are obtained from the estimated model parameters using Eqs. (5, 6). The spectral profiles s y ( q ) ( f ) are integrated within the spectral band of interest F to get estimates of self-predictability in that band, indicated as s y ( F ) .

This section reports the application of the described methods on arterial compliance time series taken from a larger database recently collected (Švec et al., 2021). The original study, approved by the ethical committee of the Jessenius Faculty of Medicine, Comenius University, included a total of 81 young and healthy Caucasians, aged 18.56 ± 2.88 years.

Arterial blood pressure (ABP) signal from finger, obtained by the photoplethysmographic volume-clamp method followed by brachial ABP reconstruction (Finometer Pro, FMS Netherlands), and electrocardiogram (ECG, CardioFax ECG-9620, NihonKohden Japan) were recorded during two phases of the experimental protocol: (i) the resting supine position (REST), started 8 min after the beginning of the measurement, and (ii) the upright position reached after passive head-up tilt (TILT), started 3 min after the position change from supine to tilt. Heart period (HP) intervals were extracted as the time distance between consecutive R peaks of the ECG. Hemodynamics parameters including cardiac output (CO) were derived on a beat-to-beat basis exploiting the impedance cardiography (ICG, CardioScreen 2000; Medis, Germany) and exerted a main role in the subsequent determination of the AC time series. The value of AC was quantified through a recently developed method, based on a reliable estimation of the time constant τ, i.e., the rate of the peripheral ABP decay during the diastolic phase, for each heart beat separately, as well as on the exploitation of the relationship between τ, AC and the total peripheral resistance (TPR) based on the two-element Windkessel model (Švec and Javorka, 2021). Since the measurement of hemodynamic parameters using ICG is very sensitive to movement artifacts, skin condition and distribution of fat, in some cases these parameters were then not determined for each heart beat, and then only 39 subjects were selected for further analysis. All the acquired signals were digitized at a sampling rate of 1 kHz. Transient changes in cardiovascular parameters between consecutive phases of the study protocol were excluded from analysis. Then, stationary segments of 300 consecutive beats were extracted from the original recordings in the two phases of the protocol. We refer the reader to (Švec et al., 2021) for further details about data acquisition and time series extraction.

The time series extracted for each subject in the two experimental conditions were regarded as realizations of the AC discrete-time process (in the following, referred to as C), assumed as uniformly sampled with a sampling frequency equal to the inverse of the mean HP. First, classical time domain markers, i.e., the mean and variance of AC ( μ C m L m m H g and σ C 2 m L 2 m m H g 2 ) were computed. Then, the series were pre-processed by removing the mean value. The AR model (1) was fitted on each pre-processed series using vector least-squares identification and setting the model order p according to the Akaike Information Criterion (AIC) (maximum scanned model order equal to 14). Since the use of the AIC sometimes led to duplicate peaks or negative power as a result of spectral decomposition (Pernice et al., 2021), the model order was manually adjusted so as to detect spectral components with positive power. After AR identification, the spectral profiles were computed according to (3). Moreover, the LF and HF components, i.e., P ^(LF)(f) and P ^(HF)(f) respectively, were computed from the poles with central frequency located in the ranges 0.04 − 0.15 Hz and 0.15 − 0.4 Hz, respectively, and the related variance was obtained from the pole residuals ( σ LF 2 and σ HF 2 ). For some subjects, more than one peak was found in these bands; in such cases, the poles with the highest power were selected for further analysis. Finally, the spectral profiles of the LSP measure in (5), computed for the LF and HF oscillations, were integrated in these bands and marked as s C ( L F ) , s C ( H F ) , respectively. The time domain LSP was obtained exploiting (6) and marked as S _C .

To test the statistical significance of the time and frequency domain LSP measures, surrogate and bootstrap data analyses were implemented as described in Sect. 2.3.1 and 2.3.2, respectively.

As regards statistical analysis, the distributions of the computed measures were tested for normality using the Anderson-Darling test. Since the hypothesis of normality was rejected for most of the distributions, and given the small sample size, non-parametric tests were employed. Specifically, the statistical significance of the difference between REST and TILT conditions, as well as between integrated PSD values in LF and HF bands in a given experimental condition, was assessed using the Wilcoxon signed-rank test for paired data. In this work, a significance level α = 0.05 was used to compute confidence intervals of the surrogate and bootstrap distributions as well as to conduct statistical tests.

The results of the time domain analysis are reported in Table 1, revealing that both the mean μ _C and the variance σ C 2 of the AC time series decreased significantly with head-up tilt (p < 0.001). This is in accordance with previous findings (Hasegawa and Rodbard, 1979; Huijben et al., 2012; Švec et al., 2021) and suggests that, when higher sympathetic activity is assumed, i.e., during the orthostatic challenge, the well-known changes of heart rate and total peripheral resistance occur rapidly through baroreflex mechanisms (Cooper and Hainsworth, 2001; Sugawara et al., 2012; Nardone et al., 2018), and are accompanied by a simultaneous rise in arterial stiffness.

Figure 3 shows the boxplot distributions of the spectral power of AC in the REST (Figure 3A) and TILT (Figure 3B) conditions computed within the LF ( σ L F 2 , green circles) and HF ( σ H F 2 , purple circles) bands, and depicted in a way such that subject-specific information relevant to the frequency location of the LF and HF spectral peaks is also provided (each circle has coordinates (f _q , σ q 2 ), where q represents the LF or HF band). While the tendency of the LF power is towards an increase moving from REST to TILT (p = 0.068), the HF power significantly decreases (p = 0.002). Furthermore, the assessment of the significance of the difference between power values integrated in LF and HF bands in a given condition revealed that the latter ones are predominant (p < 0.001) during the supine rest (Figure 3A, σ L F 2 vs. σ H F 2 ). This finding may reflect the fact that HF oscillations of AC can be heavily affected by several respiration-related mechanisms, including (i) the direct mechanical effect of intrathoracic pressure oscillations on the arterial wall stretch, and (ii) the effect of HF oscillations in AC modulators such as heart rate and ABP, with former bringing information about the mechanisms of respiratory sinus arrhythmia (RSA) (Elstad and Walløe, 2015; Švec and Javorka, 2021; Wiszt and Javorka, 2022). The predominance of HF withdraws with tilt due to a slight increase of LF power (p = 0.068) and a significant decrease of HF power (p = 0.002), as depicted in Figure 3B. An increase of magnitude of LF oscillations can reflect the sympathetically-driven vasomotion as a result of its baroreflex-mediated activation associated with orthostasis (Cooper and Hainsworth, 2001; Nardone et al., 2018; Czippelova et al., 2019). Conversely, a decrease of magnitude of HF oscillations could be attributed to the parasympathetic inhibition during orthostasis reflected by decreased RSA magnitude (Berntson et al., 1994; Javorka et al., 2018). Quantifying the effects of potential drivers of AC oscillations, such as changes of heart rate and TPR, could improve our understanding of the observed changes in AC variability (Czippelova et al., 2019; Švec et al., 2021).

In Figure 4, the spectral representation of the LSP is shown in terms of boxplot distributions of the integrated measure over all frequencies (S _C , Figure 4A), as well as in the LF ( s C ( L F ) , Figure 4B, left) and HF ( s C ( H F ) , Figure 4B, right) bands, computed in the REST (left boxplots, cyan circles) and TILT (right boxplots, magenta circles) conditions. The significant increase of the time domain LSP moving from REST to TILT (p = 0.002) is confirmed only in the HF band of the spectrum (p = 0.007). This suggests that the overall increase of regularity of the process cannot be generalized to the whole frequency content of arterial compliance, but is rather confined to the HF band of the spectrum and may have different origins.

One of them is related to the mathematical nature of the LSP measure s C ( F ) ( f ) (F represents the LF or HF band), whose spectral profile is given by the sum of the frequency-independent term S _C , and the zero-mean term s ̄ C ( F ) ( f ) showing a peak in band F. Potential tilt-induced significant increases of S _C are thus frequency-independent and distributed uniformly along the whole frequency axis. Then, the spectral LSP s C ( F ) ( f ) is influenced by this increase even in the case when there is no significant change in fluctuations of s ̄ C ( F ) ( f ) ; this influence has major weight in the larger spectral band due to the higher number of integrated frequency bins (i.e., the HF band), and may be thus responsible for the observed change of self-predictability in this band.

From a physiological point of view, the degree of complexity of arterial compliance could be a result of the combined effects of external influences modulating its dynamic activity and operating over different temporal scales, such as direct mechanical or neural influences arising from central oscillators (respiratory and vasomotor oscillators), feedback loops (e.g., baroreceptive closed-loop circuit), and complex physiologic mechanisms adjusting TPR, ABP and heart rate.

At first sight, the unaltered regularity of AC in the LF band can be attributed to a hidden tilt-induced sympathetic activation, due to the high degree of co-ordination and synchronicity of several simultaneous control mechanisms regulating AC in both the resting state and tilt conditions (e.g., heart rate, blood pressure and TPR). However, bootstrap data analysis yielded opposite results, since we found that the significance of AC self-predictability in the LF band increased with tilt (from 8/39, i.e., 20.51%, in the supine rest to 21/39, i.e., 53.85%, in tilt), as depicted in Figure 4B (left). This finding is of great importance and confirms the slight activation of sympathetic vasomotor control observed for the power spectral density of compliance (Cooper and Hainsworth, 2001; Nardone et al., 2018; Czippelova et al., 2019), besides possibly reflecting an increased amplitude of LF oscillations in external modulators such as ABP or TPR (Cooke et al., 1999; Elstad et al., 2011). Noteworthy, the augmented number of significant spectral measures in this band can be explained by considering the subject-specific frequency profiles of s C ( L F ) ( f ) , which are likely to show more prominent peaks in LF during tilt, in accordance with wider fluctuations of s ̄ C ( L F ) ( f ) and in spite of an overall frequency-independent increase of S _C (as shown in Figure 4A), i.e., the threshold for assessing significance in bootstrap data analysis (Sect. 2.3.2).

The tilt-induced physiologic responses resulting in decreased respiratory rates and increased tidal volumes (Brown et al., 1993; Porta et al., 2000; Javorka et al., 2018), as well as in a slight diminished complexity of the respiration signal (Valente et al., 2017), may be responsible for the increase of regularity of AC in the HF band. Indeed, the increased mechanical influences on arterial vessels due to an augmented tidal volume are likely to produce an augmented coupling between arterial stiffness and respiration, reflected by an increase of AC self-predictability in the HF band. Moreover, an increased HF-related regularity in AC could be attributed also to the effect of increased magnitude of ABP variability in this band (Cooke et al., 1999), probably resulting from the tilt-induced suppression of buffering effect of RSA on ABP variability at the respiratory frequency (Cooke et al., 1999; Elstad et al., 2001). Noteworthy, the latter findings highlight one important limitation of the LSP measure, which is its formulation in absence of a multivariate context taking into account potential oscillatory external drivers of AC variability, such as ABP and respiration. A spectral measure of autonomous dynamics defined as Granger Autonomy has been proposed recently, which is an extension of the LSP to the bivariate case and takes into account potential confounding mechanisms deriving from external sources (Sparacino et al., 2023). It is worth noting that the significance of HF self-predictability decreases from 30/39 (76.92%) in the supine rest to 27/39 (69.23%) during tilt, as depicted in Figure 4B (right). One more time, this result could be interpreted by looking at the spectral profiles of s C ( H F ) ( f ) : while the increase of self-predictability in HF may be associated to the frequency-independent increase of the term S _C , the tilt-induced diminished significance of the spectral LSP in the same band could be the result of less prominent peaks due to dampened fluctuations of s ̄ C ( H F ) ( f ) . Again, if combined with the increase of significance of LF regularity, this result confirms a parasympathetic withdrawal related to heart rate control and suggests the importance of LF fluctuations when the process has to cope with the physiological perturbations due to the orthostatic challenge.

The application of the proposed approach to arterial compliance data has demonstrated the significance of computing frequency-specific self-predictability measures in the case of physiological variables rich of oscillatory components with different frequencies and shape, suggesting that the overall changes of self-predictability in the time domain may be confined to specific bands of the spectrum. Moreover, investigating the spectral self dynamics of physiological processes may have a great impact in understanding their role in multivariate contexts.