0 EURASIP Journal on

Advances in Signal Processing

a SpringerOpen Journal

RESEARCH Open Access

A multivariate time-frequency method to characterize the influence of respiration over heart period and arterial pressure

Michele Orini1-2'3'4*, Raquel Bailon1,2, Pablo Laguna1,2, Luca T Mainardi3 and Riccardo Barbieri4

Abstract

Respiratory activity introduces oscillations both in arterial pressure and heart period, through mechanical and autonomic mechanisms. Respiration, arterial pressure, and heart period are, generally, non-stationary processes and the interactions between them are dynamic. In this study we present a methodology to robustly estimate the time course of cross spectral indices to characterize dynamic interactions between respiratory oscillations of heart period and blood pressure, as well as their interactions with respiratory activity. Time-frequency distributions belonging to Cohen's class are used to estimate time-frequency (TF) representations of coherence, partial coherence and phase difference. The characterization is based on the estimation of the time course of cross spectral indices estimated in specific TF regions around the respiratory frequency. We used this methodology to describe the interactions between respiration, heart period variability (HPV) and systolic arterial pressure variability (SAPV) during tilt table test with both spontaneous and controlled respiratory patterns. The effect of selective autonomic blockade was also studied. Results suggest the presence of common underling mechanisms of regulation between cardiovascular signals, whose interactions are time-varying. SAPV changes followed respiratory flow both in supine and standing positions and even after selective autonomic blockade. During head-up tilt, phase differences between respiration and SAPV increased. Phase differences between respiration and HPV were comparable to those between respiration and SAPV during supine position, and significantly increased during standing. As a result, respiratory oscillations in SAPV preceded respiratory oscillations in HPV during standing. Partial coherence was the most sensitive index to orthostatic stress. Phase difference estimates were consistent among spontaneous and controlled breathing patterns, whereas coherence was higher in spontaneous breathing. Parasympathetic blockade did not affect interactions between respiration and SAPV, reduced the coherence between SAPV and HPV and between respiration and HPV. Our results support the hypothesis that non-autonomic, possibly mechanically mediated, mechanisms also contributes to the respiratory oscillations in HPV. A small contribution of sympathetic activity on HPV-SAPV interactions around the respiratory frequency was also observed.

Keywords: Cross time-frequency analysis, Cardiovascular interactions, Baroreflex, Partial coherence, Phase difference

Introduction

Non-stationary signal processing techniques applied to electrophysiological signals have been successful in the assessment of important features of cardiovascular control physiology.

Correspondence: michele@unizar.es

1 Communications Technology Group (GTC), Aragon Institute of Engineering Research (I3A), University of Zaragoza, M de Luna 1,50018 Zaragoza, Spain

2CIBERde Bioingenieria, Biomaterialesy Nanomedicina (CIBER-BBN), Aragon

Health Sciences Institute, Zaragoza, Spain

Full list of author information is available at the end of the article

Although it is well known that respiratory activity affects cardiovascular regulation [1-4], the mechanisms responsible for the coordination between the autonomic control of circulation and respiration are still unclear and currently matter of debate [5-7]. Respiration affects both heart period length and arterial pressure, which in turn are mutually related. Interconnections between these three physiological parameters can be broadly reduced to four interactions: (i) Blood pressure decreases during inspiration and increases during expiration, following the changes in intrathoracic pressure. (ii) Heart period

Springer

© 2012 Orini et al.; licensee Springer. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

decreases during inspiration and increased during expiration, a phenomenon called respiratory sinus arrhythmia (RSA) [1,2], that may be due to direct central nervous modulation, reflex phenomena and mechanical influence of respiration [5-7]. (iii) A decrease (increase) in blood pressure provokes a decrease (increase) in heart period, via the baroreflex, a negative feedback mechanism that buffers short-term fluctuations in arterial pressure. (iv) A change in heart period causes a change in blood pressure through direct mechanical effects. As a result, the respiratory induced oscillations in the heart period and arterial pressure variability are mutually related.

To assess the functioning of the cardiovascular system in both normal and impaired conditions, it is important to reliably characterize the dynamic interactions between physiological signals involved in its regulation, i.e., heart period, arterial pressure and respiration. Cross spectral analysis has been widely used to study the cardiovascular system [3,8-10]. Cross spectral indices include spectral coherence, partial spectral coherence, phase difference, and time delay. Spectral coherence measures the strength of the coupling between the spectral components of two signals. Partial coherence measures the coherence between two signals after removing the influence that a third signal has over them [11,12]. Phase difference and time delay measure the latencies between two correlated oscillations. A limitation of spectral analysis is that it requires signals to be stationary and does not provide a dynamic characterization of the interactions among them. The purpose of this study is to present a methodology to robustly estimate the time course of cross spectral indices to characterize dynamic interactions between cardiovascular signals in non-stationary conditions. In particular, an assessment of the coupling between respiratory induced oscillations of heart period and blood pressure, as well as their interactions with respiratory activity is proposed. The characterization is based on the estimation of the time course of coherence, partial coherence, phase difference and time delay. The purpose of this study is to present a methodology to robustly estimate the time course of cross spectral indices to characterize dynamic interactions between respiration, heart period and arterial pressure, simultaneously. The characterization of the cardiovascular function is based on the estimation of the time course of coherence, partial coherence, phase difference, and time delay. Recently, distributions belonging to the Cohen's class were used to estimate time-frequency (TF) coherence between cardiovascular signals [13], and time-varying baroreflex sensitivity [14]. In this article, we present a comprehensive framework which also includes TF phase difference and TF partial coherence to simultaneously characterize the dynamic interactions between these signals during tilt

table test in spontaneous and controlled breathing, and after autonomic blockade.

Study populations

In this study, we used data from two experiments. In the first experiment, subjects performed a tilt table test under spontaneous respiration [15]. In the second experiment, subjects performed both a tilt table test and autonomic blockade under controlled respiration [3]. Tilt table test is one of the most well-established test for the assessment of cardiovascular autonomic function [9], while autonomic blockade data offers the unique possibility of studying the influence of sympathetic and parasympathetic activity separately [3]. As explained in the following, in this study we divided these data in three groups.

Tilt table test during spontaneous respiration (TTSR)

Fourteen subjects (9 males, 24-34yr, median 28 yr) without any previous cardiovascular history underwent a head-up tilt table test according to the following protocol: 4min in early supine position (TES), 5min tilted head-up to an angle of 70° (THT) and 4min back to later supine position (21s) [15].

12-lead ECG were recorded with Biopac MP150 system, with a sampling frequency of 1000 Hz. Blood pressure was measured at the finger by means of Finometer device with sampling frequency of 250 Hz. During the procedure, the Finometer was recalibrated at the beginning of THT and TLS. The recalibration took few seconds and introduced artefacts which were detected and corrected by interpolation. Arterial pressure from the finger was not corrected for the hydrostatic gradient change during tilt. Respiratory signal was recorded with a sampling frequency of 125 Hz, by using TSD201 transducer which measure thoracic expansion while breathing. This signal gives a measure correlated with lung volume changes.

Tilt table test during controlled respiration (TTCR)

Data set TTCR [3,4], composed of 14 subjects undergoing a tilt table test, is used as control for the results obtained in TTSR. Fourteen subjects (aged 19-38year, median 21 years) participated in the study according to the following protocol: 13min in supine position and 13 min head-up tilted. During the entire experiment, subjects were instructed to breath at a controlled rate, by following auditory cues. The time between the onset of two consecutive inspirations was randomly chosen from a Poisson distribution with mean 0.2 Hz. These data were originally used in [3], with the purpose of assessing the frequency response of the cardiovascular system. The rationale for randomizing the respiratory cycle was to whiten the power spectral density of the respiratory signal, which is seen as the input of the cardiovascular system.

Lung volume changes were measured with a two belt chest-abdomen inductance plethysmograph (Respitrace Systems). Arterial pressure was measured with a 22-gauge Teflon catheter in the radial artery of the nondominant hand. For each participant, intervals of about 7 min, TES and THT, recorded during supine and head-up tilt conditions, respectively, were selected for the analysis.

Autonomic blockade during controlled respiration (ABCR)

After having been tilted, subjects of TTCR data set were returned to the supine position and given either atropine or propranolol. Atropine (0.03 mg/kg) was used for parasympathetic blockade in seven subjects, while pro-pranolol (0.2 mg/kg) was used to sympathetic blockade in the other seven subjects. Physiological data were recorded in supine position following the same controlled respiratory protocol used in TTCR data set. Intervals of about 7 min, where stationarity was assumed were selected for the analysis of blockade condition. Data from TTCR recorded in supine position (TES) are used as control for autonomic blockade conditions.

Methods

Signal processing

QRS complexes in the ECG were detected by using a methodology described in [16]. Instantaneous heart rate was derived by integral pulse frequency modulation (IPFM) model, which also accounts for the presence of ectopic beats [17], and then evenly resampled at 4Hz, using spline interpolation. Instantaneous heart period was then obtained as the reciprocal of instantaneous heart rate. The heart period variability (HPV) signal, xH(t), was obtained by high pass filtering the heart period signal with a cut-off frequency of 0.03 Hz. The systolic arterial pressure series was obtained by taking the maximum of the pressure signal within a short interval following a QRS detection. The time series were subsequently interpolated at the time of occurrence of the systolic peak by splines with a sampling frequency of 4 Hz, and the SAPV signal, xS(t), was obtained by high-pass filtering with a cut-off frequency of 0.03 Hz. The signal from the thoracic belt was decimated to 4Hz to get the respiratory signal, xR(t), which gives a measure correlated with instantaneous lung volume.

Cross time-frequency analysis

In the following, signals {xt(t),xk(t)}e{xH(t),xS(t), — xR(t)}, indicate the complex analytical signal representation of HPV, SAPV and inverse respiratory signal, respectively. Owing to the inversion of the respiratory signal, and under physiological conditions, all the members of the triplet {xH(t),xS(t), —xR(t)} are expected to increase and decrease together.

The TF cross spectrum, Sik(t,f), is estimated by using a TF distribution (TFD) belonging to the Cohen's class [18]:

Sk(t,f) = &-d(t, v)A,k(r, v)e!27l(tv—ft)dvdx

Ak(T, v) = ^ 00 x, (t + 2) x* (t — 2) e—j2nvtdt,

(1) (2)

where Alk (t, v) is the narrow-band symmetric ambiguity function [19] of signals x(t) and xk(t), that in (1) is windowed by an elliptical exponential kernel, defined in the ambiguity function domain as [13]:

<Ai-d(t, v) = exp < — n

The kernel function 0d_D (t, v) can be equivalently defined in the TF domain as:

(t,f) = jj ^d-D(t, v)ej2n(tv—Tf) dTdv

The choice of the kernel function is discussed in Section 'Time-frequency filtering. Time-frequency coherence is estimated as [13]:

Yk(t,f) =

\Sk(t,f )|

A(t,f )Skk(t,f)

; Yk(t,f) e [0,1]

Time-frequency coherence quantifies the strength of the local coupling between two non-stationary signals, being Yk(t,f) = 1 in the TF regions where the signals are perfectly coupled and Yk(t,f) = 0 in the regions where signals are uncorrelated. Note that in [13] it was shown that the local averaging performed in (1) can be used to estimate non-stationary spectra, whose 4 definition includes expectation over different realizations of a given process. This implies that (5) can be used to estimate coherence function from only one pair of signals. Partial coherence is used to assess the coupling of two signals xi(t) and xk(t), after having removed the influence of a third signal xz(t) [20]. Time-frequency partial coherence can be defined as:

Ym(t,f)

Mt,f )\

y/Sm(t,f )Snz(t,f )

=_\Ski(t,f )Szz(t,f )—Skz(t,f )Szl(t,f )\_

V(Skk(t,f )Szz(t,f) — \Skz(t,f W)(S,i(t,f )Szz(t,f ) — \Sz(tf )|2)

In the TF regions where xz(t) is uncorrelated with xi(t) and xk(t), i.e., wherever Szi(t,f) = 0 and Szk(t,f) = 0, TF partial coherence is equal to TF coherence, Ywz(t,f) = Yk(t,f). Furthermore, partial coherence vanished, Yik/z(t,f) = 0, wherever xi(t) = axz(t) and

xk(t) = bxz(t), with {a, b} e R. For any triplet of signals, the interpretation of TF partial coherence is as follows: if around a TF point (t0,f0) yik/z(t,f) << yik(t,f), then it follows that the TF structure of signal xz (t) matches with that of xi(t) and xk(t) around (t0,f0). Moreover, if in a given TF region yikJz(t,f) < yzlk(t,f) < ykzn(t,f), it follows that xz(t) better represents the communality shared by the three signals [21].

Partial spectra in (6) can be obtained as:

Sik/z (t,f) = Sik(t,f) —

Szz(t,f )Szk(t,f ) Szz(t,f )

= ( |slk(t,f) | - lS'z(t,f)||Sk(t,f^(tf)

Su/z(t,f) = S(t,f) -

Szz(t,f )

Sz(t,f )Sz(t,f )

Szz(t,f )

= (1 - y^(t,f))S11(t,f)

Expression (7) shows that Sik (t,f) and Sik/z (t,f) are complex functions characterized by same phase and different magnitude. Time-frequency phase difference (TFPD) spectrum is defined as:

0lk(t,f) = arg [Sik(t,f)]

= arctan

3 [Sik(t,f)] M [Sik(t,f)]

flk(t,f) e[-n,n]

Within this framework, in a given TF region, a change in x^t) precedes (leads) a correlated change in xk (t) wherever dik(t,f) e[0,n], while a change inxt(t) lags behind a correlated change inxk(t) wherever dik(t,f) e[ -n, 0].

Note that according to the open-loop assumption of cross spectral analysis, 9ik(t,f) = -0ki(t,f). Moreover, phase spectra flRS(t,f), 9m(t,f) and 0SH(t,f) are related by flsH (t,f) = 9 RH(t,f) - flRS(t,f ). Finally, it is worth mentioning that, given two signals {xi(t), xk(t)},ifxz(t) = -xk(t) ^ 0,z(t,f) = 9,k(t,f) ± n.

Time-frequency filtering

The kernel function determines the degree of TF filtering, and, consequently, the TF resolution of the spectra and the interference terms (ITs) reduction. In this study, time resolution is quantified by At, the full width at half maximum of 0t-f(t, 0), while frequency resolution is quantified by Af, the full width at half maximum of </>t-f(0,f). These quantities measure the degree of spreading of a line in the TF domain: At and Af are equal to the full width at half maximum of the TFD of a Dirac impulse, evaluated along t for a given frequency, and of a sinusoid,

evaluated along f for a given time instant, whose ideal TF representations would be straight lines [13]. The kernel function used in this study gives a TF resolution of {At, Af} = {10.9 s, 0.039 Hz}.

The choice of the kernel is especially important in coherence analysis, because to obtain reliable coherence estimates, i.e., yik(t,f) e[0,1], the filtering provided by the kernel should be able to completely remove the ITs that characterize the Wigner-Ville distribution [13,22]. A necessary, but not sufficient, condition to obtain reliable coherence estimates is the positiveness of the auto-spectra [13]. This condition is not sufficient, since residual (oscillating) ITs present in positive TFDs may cause coherence estimates to be higher than 1. Reducing ITs and obtaining reliable coherence estimates are two aspects of the same problem. Thus, our strategy consists in finding a kernel function able to provide yik(t,f) e[0,1], which in turn implies ITs canceling.

In this study, we used a kernel function (3) that is a particular case of the multiform, tiltable exponential kernel proposed in [23]. Among all the possibilities given by this kernel, we used a function whose isocontours are, in the ambiguity function domain, ellipses with major and minor axes along t and v [24]. The choice of an elliptical shape is motivated by its good concentration around the origin of the ambiguity function domain (where auto-terms are located). Parameters v0 and t0 are used to change the length of the ellipse axes aligned along v (the degree of time filtering) and t (the degree of frequency filtering), respectively. The parameter X sets the roll-off of the filter as well as the size of the tails of the kernel. Several simulation studies demonstrated the effectiveness of this kernel in reducing ITs [13,14,25,26]. In particular, in [13], it was shown that this kernel offers the possibility of obtaining TF coherence estimates bounded between 0 and 1, as well as spectra characterized by a better TF resolution than multitaper spectrogram and continuous wavelet transform. The simulation study carried out in Section 'Simulation study' confirms these previous results.

The choice of the parameters was made as follows: First, the desired TF resolution {At,Af}, corresponding to the minimum amount of TF filtering, is decided based on a-priori information about the signals and the experimental settings. The set of parameters {t0, v0, X} that provides the desired TF resolution is used as starting point. If using this set of parameters yik(t,f) e[0,1], the degree of time (or frequency) filtering is maintained constant, while the frequency (or time) filtering is increased until reaching meaningful estimates over the entire TF domain. If at the end of the process, the frequency (or time) resolution is not satisfactory, the time (or frequency) resolution is decreased, i.e., the corresponding v0 (or t0) is increased, and the process iterates. This process allows to adjust the TF filtering to the specific needs of analysis.

In this study, we used K, v„, X} = {0.05,0.046,0.300}, with 2048 frequency points (corresponding to [At, Af} = [10.9 s, 0.039 Hz}), because among all the explored combination of parameters which gave yik(t,f) e[0,1], this was associated to the minimum degree of TF filtering. An example of this scheme was given in [14].

Estimation of synchronization indices

The time course of coherence, partial coherence, and phase difference is extracted by averaging the corresponding TF representation in specific TF regions.

The region where yik(t,f) is significant, i.e., that where the two signals are sharing approximately the same instantaneous frequencies, is defined as:

^<Y) = (t,f) e (R+ x R+) If = f() ±

^ = {^(Y) n o R(t,f );

Yik(t) = — Yk(t,f)df;

Af Jq(y)

Yik/z (t) = -1 f Ykii(t,f )df Af Jq(y)

Index dik(t) is estimated (in radians) by averaging the TFPD spectrum in Q^:

0,k(t) =

I w)df]/\Ldf

The time delay associated to dik(t) is estimated (in seconds) by index Dk(t), defined as:

0,k(t) 2nf()

^ = yt,f) e (R+ x R+) | yk(t,f) > Ym(i,f) j; (10)

where yTH(t,f) is a threshold function, estimated by using surrogate data [13], and which depends on the TF resolution of the spectra. Briefly, yTH (t,f) is obtained by estimating the TF coherence between several realizations of uncorrelated white Gaussian noises, and taking at each TF point the 95th percentile of TF coherence estimates [13].

Given that we are interested in assessing the influence of respiration on HPV and SAPV, the TF region where the time course of spectral coherence is estimated is centered around respiratory rate, fR(t), and is defined as:

where Af is a term related to the frequency resolution. Respiratory rate, fR(t), is estimated from the TF spectrum of the respiratory signal, as the frequency corresponding to the maximum of the instantaneous spectral peak.

The TF region n(k) where the time course of phase difference index is estimated is composed of those part of ) in which coherence is statistically significant:

In this expression, R(t,f) is a rectangle of sides 2sx A2f Hz and o denotes the opening (processing technique which involves erosion and dilation). The opening excludes from [n(Y) n nlk} the portions of TF domain that are smaller than R(t,f), thus adding robustness to the final estimates.

The time course of the band coherence, as well as the time course of partial coherence, is then obtained by averaging y^tf) and yikiz(t,f) in ):

A method to reduce the uncertainty of phase difference in non-stationary signals

The cross spectra Sk(t,f) are complex functions, and as such, their phase is dik(t,f) = arg[Sk(t,f)exp(j2nn)], with n e Z. The periodicity of 2n introduces an uncertainty over the actual value of 9k(t,f), which may prevent one from drawing conclusions about the temporal sequence of events described by xt(t) and xk(t). For instance, it is not possible to determine whether xi(t) precedes or lags behind xk(t), since values for dik(t) + n 2n and Vik(t) + n TR(t), where TR(t) = 1/fR(t), are positive for n > 1 and negative for n < 1. Although in cardiovascular applications the range of values for n is usually reduced to n = {-1,0,1} by considering physiological information [10], the uncertainty still remains. In the contest of non-stationary signals, we propose to use TF coherence estimates as control parameters to reduce this uncertainty. The idea is that, the lowest the time delay between non-stationary spectral components, the highest the TF coherence. This is due to the fact that TF coherence is a measure of local correlation. Thus, to determine the actual time delay among [Vlk(t) - Tr(1), V<k(t), V<k(t) + Tr©}, that is associated to phase difference 9ik(t), one can use the following procedure: (i) Generate a set of pairs of delayed signals [xi(t),xk(t + rt0)}, with r e Z and being t0 a small time delay. (ii) Estimate yik(t; r), dik(t; r) and Dik(t; r) between each pair of signals for each r, as well as their temporal median yt\r), et\r) and D,(km)(r). (iii) Find rm, as the sample for which yt}(r) is maximal. (iv) Among the possible time delays [Dik(t) - Tr(0, Vik(t), Vik(t) + Tr(^}, the closest to rmt0 is the actual one.

Statistical analysis

Statistical analysis of each data set is performed as follows. The time course of each general index from one subject s is denoted as I(t, s), with I(t, s) e [yik(t, s), Yikiz(t, s), 6k(t,s), Dik(t,s)}, with [xi(t),Xk(t)} e [Xh(t),Xs(t), —xR(t)}, and s e[1,...,N}, being N the number of subjects.

The median time course of I(t, s), I(m)(t), estimated across subjects, as well as the interquartile range, is used to describe the pattern of response of the population during a given condition.

The temporal median values of indices I(t, s), denoted as I(m)(s), are estimated during epochs where stationar-ity is assumed, and are used to: (i) Assess inter-conditions differences, i.e., whether I(m)(s) estimated during a given condition are statistically different from I(m)(s) estimated during another condition. (ii) Assess the inter-indices differences, i.e., whether Ik"' (s) estimated during a given condition are statistically different from Ipm)(s) estimated during the same condition, with (i, k) = (p, q).

Pairwise comparisons between the same indices evaluated in different epochs are performed by using the Wilcoxon signed rank test, while comparisons between different indices or between the same indices but in different data sets are performed by using the Wilcoxon ranksum test. Statistical significance is assumed for P < 0.05.

Results

Simulation study

A simulation study was carried out to verify that TFDs can be used to correctly estimate the time delay between two non-stationary spectral components of signals xi(t) and xk(t), even when 0k(t) > 2n. The SAPV signal from subjects of TTSR data were used to create pairs of signals characterized by a known time delay T:

x1(t) = xS(t, s) + ut(t) x2(t) = xs(t — T, s) + U2(t)

where T = -Jm = 3TRm), xS(t,s) is the complex analytical signal representation of SAPV of subject s, fR(m) is the median value of the respiratory rate, and ux(t) and u2(t) are two zero-mean complex white Gaussian noises associated to a SNR equal to 10 dB. The procedure described in the previous section was applied, and coherence, phase difference and time delay between {x(),x2(t + rt0)}, with r e {—30, —29,..., 30} and t0 = 0.5 s, were estimated.

Results are shown in Figure 1. In panel (a), the results from only two signals {x(), x2(t)}, obtained from the same xS(t), are shown. The highest coherence estimates (median value Y1(2m) = 0.997) were obtained for a time shift rt0 = 7.5 s, that is very close to the actual time delay T = 7.28 s. As expected, for rt0 = 7.5 s, phase differences and time delay were almost zero.

For n e {—1,0,1}, median time delay estimates V™ + nTRm) were equal to {—2.49,2.37,7.22} s, respectively. Given that coherence estimates assessed for the closest rt0 to Vm + nTRm), {—2.50,2.50,7.00} s, were equal to {0.884,0.936,0.997}, phase difference and time delay were chosen for n = 1, i.e., 012(t) = 012(t) + 2n and V12(t) = V12(t) + TRm) = 7.22 s. In Figure 1b, global results obtained considering all 14 pairs of signals are shown. Time delays T are plotted over their estimates V-™1+nTRm),

with n e {—1,0,1}. Owing to the assessment of the coherence estimates, the correct time delay values, those for n = 1, were estimated.

These results confirm those from previous studies, and demonstrate that the methodology described in Sections 'Cross time-frequency analysis' offers the possibility of estimating TF coherence and time-varying phase difference accurately [14,25,26].

Dynamic interactions during tilt table test Results from TTSR data set

An illustrative example of subject-based cross TF analysis of signals {xi(t),xk(t)} e {xH(t),xS(t), —xR(t)} by TFD is shown in Figure 2. Time-frequency coherence functions Yk(t,f), partial coherence Ywz(t,f) and phase difference spectra 9k(t,f) are shown in panels a-c, d-f, and g-i, respectively. The time course of the synchronization indices are shown below each TF representation. These indices changed during time, showing that correlations between cardiovascular oscillations are time-varying. Although the respiratory rate fR(t) changed fast during the entire test, the signals were highly correlated inside Y), which in a-c and d-f is delimited by white and black contour, respectively. The time-course of coherence indicates that the strength of the coupling, although high, was intermittent (see Yik(t) around 250 and 400 s). As expected, in (t,f) e partial coherence Ysh/r(t,f) was lower than YSH(t,f). Importantly, the influence of respiration was removed even when the respiratory rate was inside the traditional LF range, fR(t) < 0.15 Hz. The same held true for yrs/h(t,f), coherence between —xR(t) and xS(t) after having removed xH(t), and YRH/S(t,f), coherence between —xR(t) and xH(t) after having removed xS(t). Region ^*k), which in panels g-i is delimited by a black contour, only includes portions of ), where coherence estimates were statistically significant. In this particular example, differences between ) and ^ are visible during the transition from supine to head-up tilt positions and around 400 s. Interestingly, in this subject, the position changes caused fast changes in partial coherence, phase difference and time delay indices.

Indices extracted from population TTSR are given in Figure 3, where black lines and shadowed area represent median trends and interquartile ranges across the 14 subjects, respectively. For all indices, the interquartile range is low, indicating that in all subjects cardiovascular signals are related by a common mechanism. In median, spectral coherence, shown in Figure 3a-c, was very high during the entire test and decreased after the position changes. The interquartile range of YRS(t, s) was lower than that of YSH(t, s) and YRH(t, s), indicating that the coupling between respiration and SAPV was less affected by inter subject variability. Partial spectral coherence, Ywz(t,s) shown in Figure 3d-f, was always lower than spectral coherence

(a) One subject

(b) Global results

-2 0 2

Time shift rt.

Time delay T

Figure 1 Results of the simulation study. (a) Results from only two signals. Actual time delay was T = 7.28 s. Highest coherence estimates y12(f; r) were obtained for rt0 = 7.5 s. For rt0 = 7.5 s, phase difference dr2(t;r) and time delay D12(t; r) indices were almost zero. (b) Global results obtained considering all the 14 pairs of signals. Estimates D^ + nTRm) are plotted over actual time delay T, for n e {-1,0,1}. Owing to coherence analysis, the correct time delay values, those for n = 1, were selected. Red cross indicates time delay estimate for the signals whose results are shown in panel (a).

Yik(t,s). Importantly, during the test, ylk/z(i,s) changed regardless of ylk(i, s), thus showing that partial coherence reveals new information that is independent from coherence estimates. Orthostatic stress due to position changes caused phase difference and time delay, shown in Figure 3g-l, to change. The most pronounced change was that between {—xR(t), xH (t)}. Notably, changes observed in 0ik(t) were not due to changes infR(t), since Dik(t) changed as well.

To further portray the comparison between the estimated indices, the median trends X{m)(t) are shown in Figure 4, while the distributions of the temporal median, I(m)(s), are given in Figure 5. Note that in TTSR, median values X{m)(s) were estimated in the intervals which are denoted by bold gray lines in Figure 4. Numeric values are reported in Table 1. Spectral coherence between the three pairs of signals, shown in Figures 4a and 5a, followed very similar patterns. After upward and downward movement, from TES to THT and from THT to TLS, coherence abruptly decreased. Previous values were restored in about 1-2 min. During head-up tilt, yRH^s) was statistically higher than during TES and TLS (P < 5 • 10-4). During head-up tilt, the coupling between —xR(t) and xS(t) was also higher than the coupling between —xR(t) and xH(t),

yR?(s) > YR?>(s) (P < 0.023).

Partial coherence estimates, shown in Figures 4b and 5b, were always lower than the corresponding coherence estimates Yikm)(s) < y(S(s)(P < 2 • 10-5). This is not surprising since according to the physiological model, respiratory activity affects both HPV and SAPV. During head-up tilt, yR?H(s) was higher than during both

TES and TLS (P < 0.043), YR(Hm/)S(s) was lower than during TLS (P < 0.011), and ySh/R(s) was lower than during TES (P < 0.025). The increase in yRH/H(s) during head-up tilt indicates that, owing to the position change, the respiratory component in xH(t), was no longer able to explain the coupling between the respiratory components in xR(t) and in xS(t). The fact that ySHHuO was higher during TES indicates that in supine position, the respiratory induced oscillations in SAPV and HPV were not perfectly fitted by —xR(t), and that after removing respiration, the residual coupling represents mutual SAPV-HPV interactions.

The comparison between Yik/z(t, s) for different {i, k, z} shows that during Tes, yRH?s(s) < yHR(s) (P < 0.011). During Tht, YRm)s(s) < YRmH(s) (P < 3 • 10-4), and ySHR(s) < YRH/H(s) (P < 0.026). During TLS, no differences were found between partial coherence estimates. The fact that the coupling between {xS(t),xH(t)} after removing —xR(t), YsH/R(t, s), was not lower than both YRS/H(t, s) and Yrh/s(t, s) indicates that RSA and the respiratory induced oscillation in SAPV were not merely reflecting —xR(t), and that closed loop interactions between SAPV and HPV had an important role in determining the strength of the coupling between them around the respiratory frequency.

In Figures 4c and 5c, a comparison between phase difference estimates is shown. Head-up tilt caused 0{RHl(s) to statistically change from both TES and TLS (P < 0.002). Similarly, gRf (s) changes from Tes to Tht (P < 0.035). During supine position, no differences were found between 0Rf(s) and 0RRH)(s), while during head-up tilt the differences between 0i(km)(s) and ^(s), for (i, k) = (p, q), were

(a) 7sh(î, /)

(b) 7rh {t, f)

(c)7rs (£,/)

Time [si

Time [si

Time [si

Figure 2 Cross TF representations between SAPV (S), respiration (R) and HPV (H) from a subject of TTSR data set. (a)-(c) Time-frequency coherence yik(f, f); &y1 is delimited by white contour. The time course of yik(f) is reported below. (d)-(f) Partial TF coherence yik/z(t, f); 1 is delimited by black contour. The time course of yuk(t) is reported below.(g)-(i) Phase difference spectra 6ik(t, f) with (t, f) e fij®1. The time course of 9ik(t) and ^k(t) are reported below.

significant. The same changes observed in 0ik(t) are also observed in Dik(t).

Comparison between TTSR and TTCR data set

The distribution of the temporal median X{m)(s) extracted from TTCR data set is shown in Figure 5e-h, while numerical results are reported in Table 1. A segment of two representative respiratory signals, xR(t), from TTSR and TTCR data sets, as well as their TF spectra are shown in Figure 6a,b. Respiratory signals from TTSR are non-stationary and relatively narrow band, while xR(t) from TTCR are characterized by a wider instantaneous spectra.

The distribution of the median respiratory frequency,/Rm), reported in Figure 6c,d, shows that respiratory frequency had a higher inter-subject variability in TTSR, (respiratory protocol in TTCR was the same for all the participants), while intra-subject variability was higher in TTCR.

Although in the two data set the respiratory pattern differs substantially, results obtained with the same procedure are consistent. Pairwise comparisons between all indices I(m) (s) extracted from TTSR and TTCR show that differences were never significant, except for coherence estimates obtained during head-up tilt, that in TTCR were lower than in TTSR (P < 0.001).

f^ 0.8

(a) 7sn(i) [SAPV-RRV]

(b) 7rH(i) [RESP-RRV]

(c) 7rs (i) [RESP-SAPV]

f—r**\

(d) 7sH/R(i)

(e) 7rh/s (i)

(f) 7Rs/H(i)

a . U. 0.5

■1.5

r^ 0.5

n -0.5

( g) Osu(t)

ULUifcjwi IW* lyvW

(j) ®sh (i)

(h) eKU(t)

(k) VRH(t)

(1) zwt)

400 600

Time is

200 400 600 800

Time [si

Figure 3 Results from population TTSR (see Section 'Tilt table test during spontaneous respiration (TTSR)'). Black line and shadowed area denote the median trend and the interquartile range of estimates across subjects. From left to right: interactions between {xH(t), xS(f)| (SAPV-HPV), {—xR(t),xH(t)} (RESP-HPV), {—xR(t),xS(t)} (RESP-SAPV). From top to bottom: time course of: (a)-(c) coherence yik(f), (d)-(f) partial coherence yikJz(t), (g)-(i) phase difference 9ik(t) and (j)-(l) time delay ^k(t). Vertical lines denote early supine TES, head-up tilt THT and later supine TLS conditions. Dashed lines in (a)-(c) represent the mean value of the threshold function used to assess the significance of coherence estimates.

Changes during autonomic blockade

The distribution of median indices X{m)(s) from the autonomic blockade data set are shown in Figure 7, while numerical values are reported in Table 2. The comparisons between control and blockade conditions included only those subjects in which a given index was estimated for at least 3min. As expected, parasympathetic blockade due to atropine injection had a much stronger effect than propranolol sympathetic blockade, especially on the interactions that involve HPV. Atropine had virtually no effect on the interaction between {— xR(t),xS(t)}, while it reduced the strength of coupling between

{xS(t),xH(t)} (P < 0.047). Although the coupling between {—xR(t),xH(t)} decreased, this decrease was not statistically significant (P < 0.078). The decrease in partial coherence yHO from 0.623 ± 0.297 to 0.365 ± 0.172 (P < 0.047) indicates that the amount of coupling between {xS(t),xH(t)} that still remained after atropine injection was mainly due to respiration. After atropine injection, Ki^S/H(s) increased toward values of yR^s). This is likely due to the fact that after parasympathetic blockade, respiratory HPV oscillations decreased; consequently, xH(t) was no longer able to estimate the communality shared by {—xR(t),xS(t)} around respiratory frequency.

(a) Coherence

(c) Phase difference

Tes I I Tht I I TLS

1 l/^T^ ¥ U I I JzL -w ............ 2*

(b) Partial coherence

(d) Time delay

Time [s]

Figure 4 Results from TTSR data set. Median trend estimated across 14 subjects of indices. (a) Coherence yfW; (b) Phase difference 0(m)(f); (c) Partial coherence y^(f); (d) Time delay ^km)(f). Interactions between {xS(f),xH(f)}, {-xR(f),xH(f)}, |-xR(f),xS(f)} are reported in black, red, and blue, respectively. Solid gray lines in panel (a) mark the intervals in which temporal median X*m) (s) are estimated.

Due to coherence reduction, d™(s) and 0^(s) were estimated in only four subjects. In these subjects, parasympathetic blockade caused the latency between respiratory oscillations in {xS(t), xH(t)} to change of about 1.2 s. These changes were not significant, probably due to the small number of subjects involved in the comparison. Interestingly, within this framework, this change may imply a reverse of the causality in HPV-SAPV interactions, that after atropine injection would be characterized by respiratory oscillations in HPV preceding respiratory oscillations in SAPV.

Propranolol caused ySH?* (s) to slightly decrease (P < 0.078). This may indicate that sympathetic modulation is also involved in the HPV-SAPV interactions. Phase difference 6»]RIS1)(s) also slightly changed (P < 0.078), i.e., after sympathetic blockade phase differences between -xR(t) andxS(t) decreased.

Discussion

The main purpose of this study was to describe a new framework for characterizing the dynamic interactions between non-stationary signals and to demonstrate its usefulness in the analysis of cardiovascular control. The most important contributions of the study are: For the first time, non-parametric time-frequency analysis is used to explore the interactions between the three most relevant cardiovascular signals (heart period, arterial pressure and respiration) simultaneously. Non-parametric

time-frequency partial coherence is used, for the first time, in the assessment of the cardiovascular regulation. A new method for reducing the limitation of phase periodicity in the estimation of the latencies between two non-stationary spectral components is proposed. Additionally, in this article we had the unique opportunity to analyze the changes in the interactions between cardiovascular signals after autonomic blockade.

Cross time-frequency analysis

The methodology presented in this article is based on TFD and provides a robust and fast tracking of coherence, partial coherence, phase difference and time delay, by means of a procedure composed of the following steps:

(i) TF representations of spectra, coherence, partial coherence and phase differences are estimated through (1)-(9);

(ii) A statistical test is performed to localize TF regions, nik, where the strength of the local coupling between signals {xt(t),xk(t)} is significant [13]; (iii) A time-varying spectral band, n1», whose bandwidth is related to the frequency resolution of the TF representation, is defined around the respiratory rate. (iv) Time-varying spectral bands, njf, which include the portions of n(Y' in which coherence estimates are significant are localized. A minimum bandwidth and duration, given by R(t,f), is imposed to n® by the opening. (v) The time course of coherence and partial coherence is estimated by averaging the corresponding TF representations in n(Y), while phase

^r 0.4

Ï 1 ~ 0 a j» r-1

(a) Coherence [TTSR]

X......X...........x x.....X X X......>S<

(b) Partial Coherence [TTSR]

<....." N 1/ r i le *

. . , , , ,

(c) Phase difference [TTSR]

(d) Time delay [TTSR]

X.........X............x

x^S< x-X X

r 11 ri i rTi / # i t f r # t rri f § I f t f

-ÎES -l HT J LS J ES -i HT J LS -t ES -1 HT -t LS

0.9 0.8

0.8 0.6 0.4

2 1 0 -1 -2

2 1 0 -1 -2

(e) Coherence [TTCR]

X^h' X 7rh*

X 7rs^

(f) Partial Coherence [TTCR]

S <..........) C

* TSH^R X7rh/S

X 7rs/h

(g) Phase difference [TTCR]

X............x

X...........x

Xfl(m) S H

Xfl(m) "rh

Xfl(m)

(h) Time delay [TTCR]

X...........x

X X ^.........x

X2>£>

l(m) RS

X x>(m)

/ # i / i j f i ! / # i / # ] f i !

J ES -i HT -ÎES J HT J ES -i HT

Figure 5 (a)-(d): Results from the TTSR data-set. (e)-(h) Results from the TTCR data-set. Circles and bars represent the median and the interquartile range of the temporal median of indices I(t) e {yik (t), wz(t), 9]k (t), Dik(t)}. In TTSR, on the left, temporal median are estimated in early supine (TES), head-up tilt (THT) and later supine (TLS) positions. In TTCR, on the right, temporal median are estimated in supine (TES) and head-up tilt (THT) positions. Dashed lines in (a) and (e) represent the mean value of the threshold function used to assess the significance of coherence estimates.

difference changes are estimated by averaging the phase difference spectrum in fi®. (vi) Time delay is estimated to reduce the influence of instantaneous frequency over phase difference estimates.

In this study, the respiratory rate is used to localize fi00 and to convert phase difference into time delay estimates. However, in the case in which the respiratory signal was not available or the interest focused on interactions around another spectral component, the respiratory frequency should be replaced by the central frequency of the instantaneous spectral peak of interest in \Sk(t,f) |.

One of the most important aspects of this procedure is the localization of specific TF regions to estimate any given index. This offers the possibility of solving the problem of estimating the time course of indices that are reliable only around specific time-varying oscillations. Moreover, the statistical assessment of coherence level in

the phase difference estimation add further robustness to the estimates.

This TFD has been used because it offers the possibility of independently controlling the time and frequency filtering, and it is characterized by fine resolution. Its capability of accurately estimating the TF structure of physiological and synthetic real-like cardiovascular signals has been demonstrated in different studies [27-31]. In a recent study, TFD has been shown to provide more accurate coherence estimates than multitaper spectrogram and continuous wavelet transform [13]. The estimation of the phase differences in the joint TF domain was used in few studies, where it was performed by wavelet [32,33], Rihaczek [34], and reduced interference distributions [35]. Among them, no one focuses on the characterization of cardiovascular dynamics. The simulation study performed in this article shows that TFD gives accurate

Table 1 Results of tilt table test

TTSR data set TTCR data set

Index Tes Tht Tls Tes Tht

YsH (5) [nu] 0.965/0.062 0.974/0.033 0.957/0.067 0.950/0.072+ 0.906/0.078

YrH (s) [nu] 0.962/0.037 0.965/0.049 0.950/0.038 0.940/0.028+ 0.883/0.101

YrH (s) [nu] 0.974/0.044+ 0.986/0.021 0.954/0.061+ 0.952/0.032 0.940/0.052

YsHr (s) [nu] 0.650/0.166+ 0.531/0.238 0.595/0.266 0.612/0.317 0.647/0.217

ÄM [nu] 0.508/0.176 0.407/0.181 0.607/0.256+ 0.550/0.170+ 0.493/0.166

YiH (s) [nu] 0.534/0.268+ 0.710/0.280 0.577/0.197+ 0.714/0.121 0.719/0.161

[rad] 0.340/0.595 0.835/0.850 0.427/0.710 0.511/0.396+ 0.835/0.231

eRH'(s) [rad] —0.527/0.920+ 0.170/1.089 —0.621/0.705+ —0.274/0.559 —0.154/0.995

eRm'M [rad] — 1.110/1.465+ —0.447/1.291 —0.369/1.339 —0.867/0.498+ —0.546/0.490

DsH(s) [s] 0.308/1.015 0.577/0.797 0.270/0.637 0.485/0.413 0.798/0.371

D'HH(s) [s] —0.345/1.203+ 0.133/0.914 —0.419/0.852+ —0.281/0.764 —0.153/0.966

D|H)(s) [s] —0.751/2.191 + —0.326/1.740 —0.276/1.058 —0.862/0.817+ —0.509/0.547

at: values are statistically different from head-up tilt condition P < 0.05.

bFor each index I(t, s), values are given as the median and the interquartile range (med/iq) of I(m)(s), where I(m)(s) is the temporal median of I(t, s).

estimates of phase difference and time delay, and confirm previous results [14,26].

To the extent of our knowledge, this is the first time that TF partial coherence is used. Partial coherence includes information from three signals sharing a similar oscillation [20]. It was shown that in a given triplet, the lowest

partial coherence estimate is obtained when the signal which contains the shared oscillation with the better signal to noise ratio is extracted [21]. In other words, the better the extracted signal estimates (explains) the other two, the lower the partial coherence. In this study, partial coherence was the most sensitive index to orthostatic

(a) Example of xR(t) from TTCR (c) Distribution from TTCR

Time [s]

Figure 6 Respiratory protocols in TTSR and TTCR. Representative examples of respiratory signal xR(f) as well as its TF spectrum from (a) TTCR; (b) TTSR. Signal segments were recorded in supine position. (c), (d): circles and bars represent medians and the interquartile ranges of median respiratory rate (fRmed)) and interquartile range of respiratory rate (fRIQ)), during supine and head-up positions for both TTSR and TTCR data-sets.

(a) Coherence [Atropine]

(e) Phase difference [Atropine]

(g) Time delay [Atropine]

(c) Partial Coherence [Atropine]

(b) Coherence [Propranolol]

(d) Partial Coherence [Propranolol]

____^ . . .

(f) Phase difference [Propranolol]

=* ^ iM<

(h) Time delay [Propranolol]

* ^ t ^

Contr. Atrop. Contr. Atrop. Contr. Atrop.

Contr. Porpr. Contr. Porpr. Contr. Porpr.

Tsh/r Trh/s * 7rs/h

■XrO^ -XrO^

\.s -t-|(m)

"TV^sh

Figure 7 Results from autonomic blockade experiment. Circles represent the median value of a given index from a single subject, X(m)(s). Crosses represent the median of I(m)(s). On the left: comparison between indices from control condition and after atropine injection. On the right: comparison between indices from control condition and after propranolol injection. Dashed line in (a), (b) represent the mean value of the threshold function used to assess the significance of coherence estimates.

stress, in the sense that in ym(t) changes due to head-up tilt were more pronounced than in other indices, at least in spontaneous breathing.

Dependence of phase difference on signal representations

Phase difference estimates depend on how the information carried by a signal is organized on the temporal axis. Thus, the phase of a cross spectrum is sensitive to any transformation that may involve the analyzed signals, such as change of sign, inversion, derivative etc. Continuous cardiovascular signals are often obtained by interpolation of discrete values unevenly distributed in time. The interpolation process directly affects phase difference estimates. Figure 8 shows the distribution of ^m'fa) and SRJCs) in TTSR data set obtained by using three different representations of HPV. The HPV signal is estimated by interpolating on the first beat (INT1), on the second beat

(INT2) or by using the IPFM model as in this study. It is shown that by using INT1 instead of INT2, phase differences are reduced of about 1 rad, which supposing a respiratory rate of 0.2 Hz corresponds to about 800 ms, consistent with an average heart beat. The IPFM model provides estimates in between INT1 and INT2. Importantly, these differences are sufficient to change the sign of dik(t) and consequently the interpretation of the temporal sequence of events in cardiovascular interactions. A reduction of the dependence of phase difference estimates on interpolation may be obtained by modeling the heart period as a point process [36], thus providing a continuous assessment of HPV and cardiovascular interactions [37,38]. The estimation of SAPV is also critical. For instance, when arterial pressure is estimated in the periphery, its temporal structure also includes the pulse transit time [30]. Taking together, these considerations

Table 2 Results of autonomic blockade study

suggest that caution should be used in interpreting phase difference estimates, especially when they are used to infer causality. However, these issues do not affect the relative changes of the indices over time, which describe how the synchronization of physiological oscillations responds to changes in external conditions.

Limitations

The main limitation of multivariate analysis based on cross-spectral identification is the assumption of an implicit open loop model for the system. This implies that, although causal information may be inferred, to a

ABCR data set

Propranolol

0.880/0.094 0.921/0.116 0.938/0.044 0.513/0.138 0.508/0.123 0.659/0.095 0.419/0.254 -0.064/0.546 -0.609/0.662 0.327/0.272 -0.038/0.526 -0.461/0.833

certain degree, by phase difference estimation [9], our approach is not able to describe feedback and feedforward pathways of closed loops separately, neither to discriminate between direct and indirect interactions [39]. Causal coherence [40-42] or partial directed coherence [43,44], based on multivariate autoregressive modeling may be used for this purpose. However, multivariate autoregressive analysis is usually applied in stationary conditions and strongly depends on the goodness of fit of the model.

Another limitation is that spectral analysis, contrary to parametric modeling [45] or information based algorithms [46], accounts for linear interactions only.

2 1 ¥ . sh , 3 0 s* -1 —9 6$ (s) [SAPV-HPV] 0&h(«) [RESP-HPV]

< > ■a . eh ,

< < > cd « 9ï" 1 > • INT1

, -9 ( > 1 ■ IPFM ♦ INT2

hi rr? ÍT1 rri /71 -íes -íht -íls Jes -i ht -íls Figure 8 The influence of three different representations of HPVover phase difference estimates in the TTSR data set. HPV is estimated by interpolating on the first (INT1) and on the second (INT2) QRS complex, or by using the IPFM model, as in the rest of the study. On the left: phase difference between |xS(f),xH(f)|. On the right: phase difference between {— xR(f),xH(f)|. Markers and bars represent medians and interquartile ranges of median phase differences 0(m)(s).

Index Control Atropine Control

YsH (5) [nu] 0.942/0.039 0.814/0.1331 0.958/0.072

YrH (s) [nu] 0.940/0.028 0.871/0.154 0.939/0.039

yH(S) [nu] 0.961/0.031 0.961/0.018 0.944/0.022

YsHr (s) [nu] 0.623/0.297 0.365/0.172 + 0.565/0.299

yRH/'sM [nu] 0.555/0.190 0.509/0.293 0.530/0.156

yRH/H (s) [nu] 0.730/0.221 0.847/0.148 0.711/0.090

[rad] 0.646/0.385 -0.749/1.313 0.497/0.316

C'(s) [rad] -0.283/0.370 -1.495/1.285 -0.232/0.633

eRSH'M [rad] -0.836/0.240 -0.923/0.408 -0.951/0.511

DsH(s) [s] 0.652/0.391 -0.499/0.887 0.379/0.404

D'mH(s) [s] -0.281/0.379 -0.934/1.428 -0.142/0.795

D'm)(s) [s] -0.842/0.316 -0.987/0.501 -0.882/0.882

at: values are statistically different from control condition P < 0.05.

bFor each index X(t, s), values are given as the median and the interquartile range (med/iq) of I(m)(s), where I(m)(s) is the temporal median of X(t, s).

Although cardiovascular variability is known to be affected by non linear phenomena, such as saturation or hysteresis, a linear model characterizes a system by means of indices, such as power, coherence, phase difference, etc. that have a straightforward physical meaning, which in turn can be more useful for the interpretation of the ongoing processes [47].

Finally, cross spectral analysis is not able to assess hemodynamic changes during inspiration and expiration separately, neither to capture morphological information from the analyzed oscillations, thus disregarding potentially valuable information [48,49].

Physiological results

Major findings were: Intra-subject variability in indices characterizing the dynamic interactions between cardiovascular parameters shows that coupling between physiological rhythms is time-varying and intermittent. In spite of this intra-subject variability, these indices were consistent among subjects, thus suggesting the presence of common underling mechanisms of regulation between HPV, SAPV, and respiration. Systolic arterial pressure changes followed respiratory flow (see Section 'Interactions between respiration and SAPV') both in supine and standing positions and even after selective autonomic blockade. During head-up tilt, latencies between inspiratory/expiratory flow and SAPV reduction/increase augmented (see Section 'Interactions between respiration and SAPV'). Latencies between respiration and HPV were comparable to those between respiration and SAPV during supine position, and significantly increased during standing. As a result, respiratory oscillations in SAPV preceded respiratory oscillations in HPV during standing. In spontaneous breathing, partial coherence was the most sensitive index to orthostatic stress. Phase difference estimates were consistent among spontaneous and controlled breathing patterns, whereas coherence was higher in spontaneous breathing, especially during tilt.

Although in line with other studies [3,4,48,50], caution should be used in drawing conclusions from the results shown in this study, due to the small number of subjects involved in the experiments (14 in TTSR and TTCR, and 7 in ABCR data sets).

Interactions between respiration and SAPV

During inspiration and expiration, intrathoracic pressure changes cause fluctuations in the venous return and cardiac output, which in turn cause respiratory oscillations in arterial pressure [3]. Our results show that xS(t) preceded -xR(t) both in supine and head-up tilt position, both in spontaneous and controlled breathing and even after parasympathetic or sympathetic blockade. This temporal relationship was also documented in other studies

[3,48,50], and has been explained by suggesting a close mechanical dependence between arterial pressure and the derivative of respiratory volume, i.e., the respiratory flow [3]. Considering the phase shift due to the derivative, —d(xR(t))/dt, if —xR(t) and xS(t) were synchronous dRS(t) = —n/2. Since we found that dRS(t) > —n/2, our results suggest that SAPV always follows —d(xR(t))/dt with a short time delay. This delay increased in standing position [3,50], and the pattern reported in Figure 3c,d show that this change was fast. The local coupling between respirations and SAPV increases after 1-2 min from tilting during spontaneous breathing but decreased during controlled breathing. Interestingly, during head-up tilt, respiratory SAPV better explained the coupling between all the three variables. Autonomic blockade did not affect the interactions between SAPV and respiration significantly, even if it seems that propranolol slightly increased the latencies between them.

Interactions between respiration and HPV

Respiratory sinus arrhythmia, the heart rate variability in synchrony with respiration [2], has both clinical and physiological relevance [1]. Several hypothesis have been made to explain the origin of the respiratory oscillations in HPV. They include a central mechanism [5,11], the baroreflex [7,8,51], the mechanical stretch of the sinus node [3], which is enhanced during exercise [52], and a mixture of them [6,48,50]. However, most of the mechanisms underlying RSA are still unclear, and there is no general consensus about its origin [5]. Importantly, the estimation of the time delay between HPV and SAPV around the respiratory frequency is a central point of this debate [5]. Although the discussion of the origin of RSA is beyond the scope of this study, it is important to stress that the presented methodology, which provides robust and reliable phase difference estimates, may be useful to gain some insight into this debate.

Our results show that xH(t) preceded —xR(t), probably due to a similar mechanism as that illustrated in the previous section, i.e., respiratory oscillations in HPV followed changes in the respiratory flow. The fact that changes in HPV anticipated changes in the instantaneous lung volume was also documented in [3,48,50]. Change from supine to standing position increased the time delay between respiratory flow and HPV. This increase may be related to an increase of baroreflex contribution to the respiratory oscillations in HPV [3].

Partial coherence analysis shows that during standing, HPV was not able to explain the coupling between respiration and SAPV, likely due to a decrease in the amplitude of respiratory HPV. Parasympathetic blockade reduced the strength of the coupling between HPV and respiration. However, this decrease was not statistically significant and did not reduce coherence estimates below the

significance threshold. The amount of HPV not abolished by atropine injection was characterized by a phase difference with respect to —xR(t) close to —n/2, which indicates very short time delay with respect to respiratory flow. This supports the hypothesis that non-autonomic, possibly mechanically mediated, mechanisms also contributes to RSA [3].

Interactions between SAPVand HPV

With intact autonomic modulation the removal of respiration influence reduced, but did not abolish, the coupling between respiratory oscillation in SAPV and in HPV. This indicates that some correlation that cannot be perfectly captured by respiration is present in the closed loop between HPV and SAPV around the respiratory frequency. Parasympathetic blockade dampened the local coupling below significance level. Since during parasympathetic blockade vagal baroreflex is abolished, while mechanical influences of respiration on SAPV are intact, the residual interactions between HPV and SAPV are likely due to correlation between residual mechanical oscillations in HPV and respiratory induced oscillations in SAPV. In absence of vagal control, phase differences were consistent to the direction of the feedforward path. During parasympathetic blockade, the amount of residual coupling almost disappeared after removing the influence of respiration. Interestingly, propranolol slightly reduced the strength of the coupling, suggesting a small contribution of sympathetic activity on HPV-SAPV interactions around the respiratory frequency.

Conclusions

In this paper, a new framework for characterizing the dynamic interactions between non-stationary signals was presented and its usefulness in the analysis of cardiovascular control was demonstrated. The most important contributions of this methodology are: For the first time, non-parametric time-frequency analysis is used to explore the interactions between the three most relevant cardiovascular signals (heart period, arterial pressure and respiration) simultaneously. Non-parametric time-frequency partial coherence is used, for the first time, in the assessment of the cardiovascular regulation. A new method for reducing the limitation of phase periodicity in the estimation of the latencies between two non-stationary spectral components is proposed. Additionally, in this article we had the unique opportunity to analyze the changes in the interactions between cardiovascular signals after autonomic blockade.

Competing interests

The authors declare that they have no competing interests. Acknowledgements

Thanks to Dr. JP Saul and Dr. A Minchole for providing the authors with TTSR, TTCB, and ABCRdata sets. This study was supported in part by NIH Grants

R01-HL084502, Ministerio de Ciencia e Innovación, Spain, under Projects TEC2010-21703-C03-02 and TRA2009-0127, Diputación General de Aragón, Spain, through Grupos Consolidados GTC ref:T30, Instituto de Salud Carlos III, Spain, through CIBER CB06/01/0062, and ARAID and Ibercaja under project "Programa de apoyo a la I+D+i".

Author details

1 Communications Technology Group (GTC), Aragon Institute of Engineering Research (I3A), University of Zaragoza, M de Luna 1, 50018 Zaragoza, Spain. 2CIBERde Bioingeniería, Biomaterialesy Nanomedicina (CIBER-BBN), Aragon Health Sciences Institute, Zaragoza, Spain.3Department of Bioengineering, Politecnico di Milano, Via Golgi 39,20131 Milano, Italy. 2Neuroscience Statistics Research Laboratory, Massachusetts General Hospital, Harvard Medical School, Boston, MA 02114, USA.

Received: 13 December 2011 Accepted: 10 September 2012

Published: 4 October 2012

References

1. P Grossman, EWTaylor, Toward understanding respiratory sinus arrhythmia: relations to cardiac vagal tone, evolution and biobehavioral functions. Biol. Psychol. 74(2), 263-285 (2007)

2. F Yasuma, J Hayano, Respiratory sinus arrhythmia*: why does the heartbeat synchronize with respiratory rhythm? Chest. 125(2), 683-690 (2004)

3. JP Saul, RD Berger, P Albrecht, SP Stein, MH Chen, RJ Cohen,Transfer function analysis of the circulation: unique insights into cardiovascular regulation. Am. J. Physiol. Heart Circ. Physiol. 261 (4), H1231-H1245 (1991)

4. R Barbieri, A Bianchi, J Triedman, L Mainardi, S Cerutti, J Saul, Model dependency of multivariate autoregressive spectral analysis. IEEE Eng. Med. Biol. Mag. 16(5), 74-85 (1997)

5. DL Eckberg, Pointcounterpoint: respiratory sinus arrhythmia is due to a central mechanism vs. respiratory sinus arrhythmia is due to the baroreflex mechanism. J. Appl. Physiol. 106(5), 1740-1742 (2009). discussion 1744

6. C Julien, MJ Parkes, SYCTzeng, PYW Sin, PN Ainslie, P van de Borne, JO Fortrat, MA Custaud, C Gharib, A Porta, F Vallais, G Baselli, M Pagani, D Lucini, RL Hughson, JA Taylor, CO Tan, DM Baekey, TE Dick, JFR Paton, B Taha, Comments on point:counterpoint: respiratory sinus arrhythmia is due to a central mechanism vs. respiratory sinus arrhythmia is due to the baroreflex mechanism. J. Appl. Physiol. 106(5), 1745-1749 (2009)

7. JM Karemaker, Last word on point:counterpoint: respiratory sinus arrhythmia is due to a central mechanism vs. respiratory sinus arrhythmia is due to the baroreflex mechanism. J. Appl. Physiol. 106(5),

1750 (2009)

8. RW deBoer, JM Karemaker, J Strackee, Hemodynamic fluctuations and baroreflex sensitivity in humans: a beat-to-beat model. Am. J. Physiol. 253(3 Pt 2), H680-H689 (1987)

9. WH Cooke, JB Hoag, AA Crossman, TA Kuusela, KUOTahvanainen, DL Eckberg, Human responses to upright tilt: a window on central autonomic integration. J. Physiol. 517(2), 617-628 (1999)

10. A Porta, AM Catai, ACM Takahashi, V Magagnin, T Bassani, E Tobaldini, P van de Borne, N Montano, Causal relationships between heart period and systolic arterial pressure during graded head-up tilt. Am. J. Physiol. Regul. Integr. Comp. Physiol. 300(2), R378-R386 (2011)

11. LJ Badra, WH Cooke, JB Hoag, AA Crossman, TA Kuusela, KUO Tahvanainen, DL Eckberg, Respiratory modulation of human autonomic rhythms. Am. J. Physiol. Heart Circ. Physiol. 280(6), H2674-H2688 (2001)

12. PD Larsen, CD Lewis, GL Gebber, S Zhong, Partial spectral analysis of cardiac-related sympathetic nerve discharge. J. Neurophysiol. 84(3), 1168-1179(2000)

13. M Orini, R Bailon, L Mainardi, P Laguna, P Flandrin, Characterization of the dynamic interactions between cardiovascular signals by time-frequency coherence. IEEE Trans. Biomed. Eng. 59(3), 663-673 (2012)

14. M Orini, P Laguna, LT Mainardi, R Bailon, Assessment of the dynamic interactions between heart rate and arterial pressure by the cross time-frequency analysis. Physiol. Meas. 33(3), 315-331 (2012)

15. A Minchole, E Pueyo, JF Rodríguez, E Zacur, M Doblare, P Laguna, Quantification of restitution dispersion from the dynamic changes of the T-wave peak to end, measured at the surface ECG. IEEE Trans. Biomed. Eng. 58(5), 1172-1182 (2011)

16. JP Martinez, R Almeida, S Olmos, AP Rocha, P Laguna, A wavelet-based ECG delineator: Evaluation on standard databases. IEEE Trans. Biomed. Eng. 51(4), 570-581 (2004)

17. J Mateo, P Laguna, Analysis of heart rate variability in the presence of ectopic beats using the heart timing signal. IEEE Trans. Biomed. Eng. 50, 334-343 (2003)

18. F Hlawatsch, GF Boudreaux-Bartels, Linear and quadratic time-frequency signal representations. IEEE Signal Process. Mag. 9(2), 21-67 (1992)

19. P Flandrin, Time-Frequency Signal Analysis and Processing, chap. Ambiguity Function. (Philadelphia, Elsevier, 2003)

20. J Bendat, A Piersol, Random Data: Analysis and Measurement Procedures. (New York, John Wiley & Sons, 2010). Wiley Series in Probability and Statistics

21. LA Baccala, K Sameshima, Comments on 'Is partial coherence a viable technique for identifying generators of neural oscillations?': Why the term 'Gersch Causality' is inappropriate: common neural structure inference pitfalls. Biol. Cybern. 95(2), 135-141 (2006)

22. G Matz, F Hlawatsch, Time-frequency coherence analysis of nonstationary random processes. in Proc. Tenth IEEE Workshop Statistical Signal and Array Processing. (Pennsylvania, PA, USA, 2000) pp. 554-558

23. A Costa, G Boudreau-Bartels, Design of time-frequency representations using a multiform, tiltable exponential kernel. IEEE Trans. Signal Process. 43(10), 2283-2301 (1995)

24. F Hlawatsch, P Flandrin, The Wigner Distribution - Theory and Applications in Signal Processing. (Philadelphia, Elsevier, 1997) pp. 59-113. Elsevier 1997 chap. The interference structure of the Wigner distribution and related time-frequency signal representations

25. M Orini, R Bailon, L Mainardi, A Minchole, P Laguna, Continuous quantification of spectral Coherence using Quadratic Time-Frequency distributions, error analysis and application. in Internat. Conf. Computers in Cardiology. (Park City, UT, USA, 2009) pp. 681-684

26. M Orini, R Bailon, LT Mainardi, P Laguna,Time-frequency phase differences and phase locking to characterize dynamic interactions between cardiovascular signals. in Conf. Proc. IEEE Eng. Med. Biol. Soc. (vol 2011, Boston, MA USA, 2011) pp. 4689-4692

27. S Pola, A Macerata, M Emdin, C Marchesi, Estimation of the power spectral density in nonstationary cardiovascular time series: assessing the role of the time-frequency representations (TFR). IEEE Trans. Biomed. Eng. 43, 46-59 (1996)

28. M Orini, R Bailon, R Enk, S Koelsch, LT Mainardi, P Laguna, A method for continuously assessing the autonomic response to music-induced emotions through HRV analysis. Med. Biol. Eng. Comput. 48(5), 423-433 (2010)

29. R Bailon, L Mainardi, M Orini, SP Land Laguna, Analysis of heart rate variability during exercise stress testing using respiratory information. Biomed. Signal Process Control. 5,299-310 (2010)

30. E Gil, M Orini, R Bailon, JM Vergara, L Mainardi, P Laguna, Photoplethysmography pulse rate variability as a surrogate measurement of heart rate variability during non-stationary conditions. Physiol. Meas. 31(9), 1271 (2010)

31. M Orini, R Bailon, L Mainardi, P Laguna, Synthesis of HRV signals characterized by predetermined time-frequency structure by means of time-varying ARMA models. Biomed. Signal Process Control. 7, 141-150(2012)

32. JP Lachaux, E Rodriguez, J Martinerie, FJ Varela, Measuring phase synchrony in brain signals. Hum Brain Mapp. 8(4), 194-208 (1999)

33. MLV Quyen, J Foucher, J Lachaux, E Rodriguez, A Lutz, J Martinerie, FJ Varela, Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony. J. Neurosci. Methods. 111(2), 83-98(2001)

34. S Aviyente, EM Bernat, WS Evans, SR Sponheim, A phase synchrony measure for quantifying dynamic functional integration in the brain. Hum Brain Mapp. 32(1), 80-93 (2010)

35. YJ Shin, D Gobert, SH Sung, EJ Powers, JB Park, Application of cross time-frequency analysis to postural sway behavior: the effects of aging and visual systems. IEEE Trans. Biomed. Eng. 52(5),

859-868 (2005)

36. R Barbieri, EC Matten, AA Alabi, EN Brown, A point-process model of human heartbeat intervals: new definitions of heart rate and heart rate variability. Am. J. Physiol. Heart Circ. Physiol. 288,

H424-H435 (2005)

37. Z Chen, EN Brown, R Barbieri, Assessment of autonomic control and respiratory sinus arrhythmia using point process models of human heart beat dynamics. IEEE Trans. Biomed. Eng. 56(7), 1791-1802 (2009)

38. Z Chen, P Purdon, G Harrell, E Pierce, J Walsh, E Brown, R Barbieri, Dynamic assessment of baroreflex control of heart rate during induction of propofol anesthesia using a point process method. Annals Biomed. Eng. 39, 260-276(2011)

39. R Barbieri, G Parati, J Saul, Closed- versus open-loop assessment of heart rate baroreflex. IEEE Eng. Med. Biol. Mag. 20(2), 33-42 (2001)

40. H Zhao, WA Cupples, KH Ju, KH Chon, Time-varying causal coherence function and its application to renal blood pressure and blood flow data. IEEE Trans. Biomed. Eng. 54(12), 2142-2150 (2007)

41. A Porta, R Furlan, O Rimoldi, M Pagani, A Malliani, P van de Borne, Quantifying the strength of the linear causal coupling in closed loop interacting cardiovascular variability signals. Biol. Cybern. 86(3), 241-251 (2002)

42. G Nollo, L Faes, A Porta, R Antolini, F Ravelli, Exploring directionality in spontaneous heart period and systolic pressure variability interactions in humans: implications in the evaluation of baroreflex gain. Am. J. Physiol. Heart Circ. Physiol. 288(4), H1777-H1785 (2005)

43. LA Baccala, K Sameshima, Partial directed coherence: a new concept in neural structure determination. Biol. Cybern. 84,463-474 (2001)

44. L Faes, A Porta, G Nollo, Testing frequency-domain causality in multivariate time series. IEEE Trans. Biomed. Eng. 57(8), 1897-1906 (2010)

45. Z Chen, EN Brown, R Barbieri, Characterizing nonlinear heartbeat dynamics within a point process framework. IEEE Trans. Biomed. Eng. 57(6), 1335-1347(2010)

46. L Faes, G Nollo, A Porta, Information domain approach to the investigation of cardio-vascular, cardio-pulmonary and vasculo-pulmonary causal couplings. Frontiers Physiol. 2(80), 1-13 (2011)

47. X Xiao, TJ Mullen, R Mukkamala, System identification: a multi-signal approach for probing neural cardiovascular regulation. Physiol. Meas. 26(3), R41-71 (2005)

48. PYW Sin, DC Galletly, YCTzeng, Influence of breathing frequency on the pattern of respiratory sinus arrhythmia and blood pressure: old questions revisited. Am. J. Physiol. Heart Circ. Physiol. 298(5), H1588-H1599 (2010)

49. YCTzeng, PYW Sin, DC Galletly, Human sinus arrhythmia: inconsistencies of a teleological hypothesis. Am. J. Physiol. Heart Circ. Physiol. 296, H65-H70 (2009)

50. M Elstad, KToska, KH Chon, EA Raeder, RJ Cohen, Respiratory sinus arrhythmia: opposite effects on systolic and mean arterial pressure in supine humans. J. Physiol. 536(Pt 1), 251-259 (2001)

51. JM Karemaker, Counterpoint: respiratory sinus arrhythmia is due to the baroreflex mechanism. J. Appl. Physiol. 106(5), 1742-1743 (2009). discussion 1744

52. G Blain, O Meste, S Bermon, Influences of breathing patterns on respiratory sinus arrhythmia in humans during exercise. Am. J. Physiol. Heart Circ. Physiol. 288(2), H887-95 (2005)

doi:10.1186/1687-6180-2012-214

Cite this article as: Orini etal.: A multivariate time-frequency method to characterize the influence of respiration over heart period and arterial pressure. EURASIP Journal on Advances in Signal Processing 2012 2012:214.