## Abstract

Heart rate (HR) power spectral indexes are limited as measures of the cardiac autonomic nervous systems (CANS) in that they neither offer an effective marker of the β-sympathetic nervous system (SNS) due to its overlap with the parasympathetic nervous system (PNS) in the low-frequency (LF) band nor afford specific measures of the CANS due to input contributions to HR [e.g., arterial blood pressure (ABP) and instantaneous lung volume (ILV)]. We derived new PNS and SNS indexes by multisignal analysis of cardiorespiratory variability. The basic idea was to identify the autonomically mediated transfer functions relating fluctuations in ILV to HR (ILV→HR) and fluctuations in ABP to HR (ABP→HR) so as to eliminate the input contributions to HR and then separate each estimated transfer function in the time domain into PNS and SNS indexes using physiological knowledge. We evaluated these indexes with respect to selective pharmacological autonomic nervous blockade in 14 humans. Our results showed that the PNS index derived from the ABP→HR transfer function was correctly decreased after vagal and double (vagal + β-sympathetic) blockade (*P* < 0.01) and did not change after β-sympathetic blockade, whereas the SNS index derived from the same transfer function was correctly reduced after β-sympathetic blockade in the standing posture and double blockade (*P* < 0.05) and remained the same after vagal blockade. However, this SNS index did not significantly decrease after β-sympathetic blockade in the supine posture. Overall, these predictions were better than those provided by the traditional high-frequency (HF) power, LF-to-HF ratio, and normalized LF power of HR variability.

- heart rate variability
- modeling
- power spectral analysis
- system identification
- transfer function

over two decades ago, Cohen and co-workers demonstrated that power spectral analysis of resting heart rate (HR) variability could provide more specific measures of the cardiac autonomic nervous systems (CANS) than a basic single time statistical analysis (e.g., mean and SD) (1, 2, 26). Using selective pharmacological autonomic nervous blockade, these investigators demonstrated that low-frequency (LF; 0.04–0.15 Hz) fluctuations in HR are jointly mediated by the parasympathetic nervous system (PNS) and β-sympathetic nervous system (SNS), whereas high-frequency (HF; 0.15–0.4 Hz) fluctuations are solely mediated by the PNS. It therefore directly followed that the HF power of HR variability could serve as a simple quantitative index of the PNS. Shortly thereafter, Malliani et al. (18) proposed the ratio of the LF power to HF power of HR variability [LF-to-HF ratio (LF/HF ratio)] as an index of sympathovagal balance as well as normalized LF and HF powers as indexes of the CANS independent of the total power or variance. Since then, these easily obtainable indexes have been extensively studied under a wide variety of pathophysiologic conditions (31), having most notably been shown to be associated with mortality in various patient classes [e.g., postmyocardial infarction (10), cardiac arrest (13), and critically ill (36)] as well as to be altered by various disease processes [e.g., diabetic autonomic neuropathy (14, 23), congestive heart failure (11, 29), and hypertension (15, 16)]. Today, the HF power of HR variability is generally accepted as a useful index of the PNS, whereas the LF/HF ratio and normalized LF power are sometimes considered to be indexes of sympathovagal balance or the SNS (31).

Despite their widespread use and established association with disease, HR power spectral indexes are quite controversial in terms of being actual measures of the CANS (18, 24). Part of this controversy simply stems from the common misconception that these indexes should be indicative of basal or mean autonomic outflow when they are, in fact, intended to reflect the fluctuations in autonomic outflow about the mean value. The remaining part of the controversy surrounding HR power spectral indexes is based on their intrinsic limitations as measures of such fluctuations. One of the major limitations is that HR power spectral indexes, including the LF/HF ratio and normalized LF power, generally cannot provide an effective marker of the SNS (i.e., fluctuations in SNS outflow), because the LF power is mediated by both CANS, as shown in the initial study by Cohen and co-workers and verified in a subsequent study (25). Another major limitation, alluded to above, is that none of the HR power spectral indexes is truly specific to the CANS, since HR is just an output of these systems. For example, it is well known that the ultimate source of the HF power of HR variability is respiration. Thus, changes in the HF power could also reflect changes in the respiratory input. Similarly, changes in the LF power of HR variability could also reflect changes in arterial blood pressure (ABP) fluctuations, which may cause HR fluctuations through the baroreflex (1). So, for example, the HF power may not provide a useful means to monitor the PNS during exercise (18), whereas the LF/HF ratio may not be an effective indicator of sympathovagal balance when local vasomotor activity alters the variability of the ABP input (1, 40).

As a result, investigators have recently sought improved indexes of the CANS by employing more sophisticated analyses of beat-to-beat cardiovascular variability. Vetter et al. proposed to selectively quantify the PNS and SNS using a blind source separation analysis of HR and ABP variability (33) as well as HR and QT interval variability (34). Chon and co-workers then introduced a principal dynamic mode analysis of HR variability to separately characterize the two CANS (40, 41). The common idea underlying each of these analyses is to distinguish the SNS from the PNS by invoking some sort of independence or orthogonality assumption, which is generally untenable. Despite this assumption, the results of these analyses were quite promising and sparked our own interest in pursuing improved indexes of the PNS and SNS.

In this article, we derive indexes of the PNS and SNS based on a multisignal analysis of the beat-to-beat variability in potentially noninvasively measured HR, ABP, and instantaneous lung volume (ILV). The basic idea of the analysis, which was briefly proposed by us with co-workers in Ref. 38 (see the discussion), was to identify the transfer functions relating the fluctuations in ILV and ABP to HR so as to effectively eliminate the input contributions to HR and then separate the estimated transfer functions in the time domain into PNS and SNS indexes using physiological insight rather than any independence or orthogonality assumption. We then showed that these new indexes are better than traditional HR power spectral indexes in terms of predicting the known effects of selective pharmacological autonomic nervous blockade on the two CANS in 14 healthy human subjects.

Our multisignal analysis of cardiorespiratory variability aims to derive specific indexes of the PNS and SNS from the dynamic coupling relating fluctuations in ILV to HR (ILV→HR), which represents the autonomically mediated respiratory sinus arrhythmia phenomena (21), and the dynamic coupling relating the fluctuations in ABP to HR (ABP→HR), which characterizes the autonomically mediated HR baroreflex (21). [Note that ILV→HR coupling is actually noncausal, perhaps due to coupling between the respiratory and HR control centers in the brain (21, 30, 32).] The analysis specifically invokes two key linearity assumptions that are justified by previous physiological findings. The first assumption is that ILV→HR and ABP→HR couplings may be well approximated with linear and time-invariant transfer functions for small, resting fluctuations. This assumption is supported by the work of Chon et al., who showed that 70–75% of the variance of resting HR variability could be explained by linear couplings, whereas only an additional 10–15% of the variance could be attributed to second-order nonlinear couplings (12). The second assumption is that the ILV→HR and ABP→HR transfer functions may each be modeled in its time domain impulse response form as the superposition of early and fast PNS component and a delayed and sluggish SNS component, as initially observed by Triedman et al. (32). This assumption is buttressed by known physiology and the comparison shown in Fig. 1 of the impulse responses relating pure vagal and sympathetic stimulation to HR fluctuations that were experimentally determined by Berger et al. (8) in dogs with typical ILV→HR and ABP→HR impulse responses that were identified by analysis of cardiorespiratory fluctuations from humans at rest (21). [Note that the data of Berger et al. actually revealed that the sympathetic impulse response is delayed by ∼2 s with respect to the vagal impulse response (8).] It should be further noted that the comparison shown in Fig. 1 effectively represents an extrapolation of the efferent autonomic nervous limbs in dogs to the afferent, central, and efferent autonomic nervous limbs in humans. The results shown in Fig. 1 therefore suggest that the frequency dependency of the ILV→HR and ABP→HR impulse responses can be largely attributed to efferent autonomic nervous limbs and that the CANS of dogs can well represent the CANS of humans. Based on these two linearity assumptions, the multisignal analysis of cardiorespiratory variability is employed in three steps.

In the first step, ILV→HR and ABP→HR impulse responses are estimated by applying parametric system identification to resting beat-to-beat fluctuations in HR, ABP, and ILV according to the block diagram shown in Fig. 2. In this way, the known correlation between ILV and ABP inputs is accounted for by simultaneously considering the contributions of both inputs to HR (i.e., dual-input identification), and the feedback baroreflex effects of fluctuations in ABP on HR may be disentangled from the feedforward mechanical effects of fluctuations in HR on ABP through the imposition of causality (i.e., open-loop identification of closed-loop systems) (35). The block diagram also includes a perturbing noise source (*N*_{HR}), which is also estimated and represents the residual HR variability not accounted for by the two inputs. This step is mathematically achieved with the following double-exogenous model with autoregressive (AR) input (27) originally proposed by Baselli et al. (4, 5): (1) where *n* is discrete time. The sets of unknown finite impulse response (FIR) parameters {*g*[*i*], *h*[*i*]}, respectively, define ILV→HR and ABP→HR impulse responses, whereas the unmeasured residual white noise error (*W*_{HR}) and the set of unknown AR parameters {*a*[*i*]} specify *N*_{HR}. The terms *m*, *M*, and *N* represent the maximum expected durations of the FIRs, whereas *P* limits the number of AR parameters. Note that the causality of the ABP→HR impulse response is enforced here by virtue of setting the *h*[*i*] parameters to 0 for negative values of *i*, whereas the noncausality of the ILV→HR impulse response is accounted for by setting *m* to a negative value. The parameter sets are estimated from ∼6-min intervals of zero-mean fluctuations in HR, ABP, and ILV sampled at 1.5 Hz based on least-squares minimization of *W*_{HR} using a weighted-principal component regression method (37) (see the appendix). Prior to the application of this method, HR and ABP are normalized by their mean values and ILV is normalized by its SD, as described in Ref. 38. Thus, the identified impulse responses will have units of s^{−1} and represent the relationship between relative fluctuations in the input about its mean value and relative fluctuations in the output about its mean value.

In the second step, the estimated ILV→HR and ABP→HR impulse responses are each separated into PNS and SNS components as generally shown in Fig. 1. More specifically, the typical ILV→HR impulse response here indicates that, in response to an impulse increase in ILV at *time 0*, HR would first increase due to early (noncausal) and fast parasympathetic withdrawal and then decrease due to delayed and sluggish β-sympathetic withdrawal. Thus, the PNS and SNS components are of opposite direction in amplitude and delayed in time with respect to each other. As a result, the PNS component of the ILV→HR impulse response is precisely determined as the initial portion of the impulse response extending up to the first 0 crossing (as this marker specifies the amplitude direction change) following the time of the peak amplitude that occurs within the first 5 s (as the PNS is fast), whereas the SNS component is specified as the remaining latter portion of the impulse response (as the SNS is slow and delayed). The typical ABP→HR impulse response shown in Fig. 1 indicates that, in response to an impulse increase in ABP at *time 0*, HR would decrease first due to early (causal) and fast parasympathetic activation and then due to delayed and sluggish β-sympathetic withdrawal. Thus, the PNS and SNS components here are also delayed in time with respect to each other but are of the same direction in amplitude. As a result, the separation procedure is not as straightforward for the ABP→HR impulse response. In particular, the initial downstroke of this impulse response is regarded as the first part of the PNS component. The time interval of the initial downstroke is precisely defined up to the minimum amplitude that occurs within the first 5 s (as the PNS is fast). The remaining part of the PNS component (i.e., the return of this stable component to 0), which may be obscured by the subsequent SNS component, is assumed to be symmetric to the visible first part, as suggested by the data of Berger et al. shown in Fig. 1. The SNS component is then established by subtracting the total PNS component from the entire ABP→HR impulse response (as the SNS is slow and delayed). Note that the separation procedure for the ILV→HR impulse response assumes that the PNS and SNS components do not overlap in time, whereas this separation procedure makes no such assumption.

In the third step, scalar indexes of the PNS and SNS are determined from the corresponding components of the ILV→HR and ABP→HR impulse responses. Specifically, the two norms (i.e., square root of the energy) of the PNS and SNS components are calculated to arrive at a pair of PNS and SNS indexes denoted as P_{1} and S_{1} as determined from the ILV→HR impulse response and P_{2} and S_{2} as determined from the ABP→HR impulse response. Note that these indexes will be dimensionless due to the aforesaid signal normalization.

## MATERIALS AND METHODS

#### Experimental data.

We evaluated P_{1} and P_{2} indexes of the PNS and S_{1} and S_{2} indexes of the SNS with respect to a previously collected human cardiorespiratory dataset. The materials and methods utilized in the collection of this dataset are described in detail elsewhere (30). The experimental protocol was approved by the MIT Committee on the Use of Humans as Experimental subjects, and written informed consent was obtained from the subjects.

Briefly, 14 healthy male subjects (age: 19–38 yr) participated. Surface ECG, ABP (radial artery catheter), and ILV (chest-abdomen inductance plethysmography) were continuously recorded at a sampling rate of 360 Hz. Throughout the recordings, the subjects initiated each respiratory cycle on cue to a sequence of audible tones spaced at random intervals of 1–15 s with a mean of 5 s (7). The purpose of this random-interval breathing protocol was to produce a broadband ILV signal so as to enhance the reliability of subsequent system identification analyses, especially pertaining to the estimation of couplings involving ILV as the input. Approximately 13-min intervals of the random-interval breathing cardiorespiratory measurements were recorded for each of six experimental conditions. First, control measurements were recorded from the subjects in supine and standing postures. Then, seven of the subjects received atropine at 0.03 mg/kg iv (vagal blockade), and measurements were recorded from the subjects in the two postures. Finally, these subjects received propranolol at 0.2 mg/kg iv (β-sympathetic blockade) to achieve total cardiac autonomic nervous blockade (double blockade), and measurements were again recorded from the subjects in both postures. The remaining seven subjects were studied similarly but received the drugs in reverse order. A 5-min period for hemodynamic equilibration was allowed between each experimental condition.

#### Data analysis

As described in Refs. 6 and 30, ABP and ILV were digitally filtered and decimated to 1.5 Hz, and synchronized 1.5 Hz HR tachograms were constructed from surface ECGs. PNS and SNS indexes were then derived as described above from two nonoverlapping intervals of each ∼13-min record in the cardiorespiratory dataset. The resulting pairs of each index from each data record were then averaged. For comparison, the HF power, LF/HF ratio, and normalized LF power [LF/(LF+HF)] of HR tachograms were analogously computed for each data record using AR power spectral analysis (31) with postfiltering (6). [HF and LF bands were defined here to range from 0.04 to 0.15 Hz and from 0.15 to 0.40 Hz, respectively (31).] For each of the seven investigated indexes, outliers defined to be >2 SDs from the mean value were removed, wherein means and SDs were robustly estimated with the Gaussian-based method of Armoundas et al. (3). Finally, paired *t*-tests were performed to determine if the indexes were significantly altered from the control condition to each drug condition.

## RESULTS

Figures 3 and 4 show group average results (mean ± SE) for the estimated ILV→HR and ABP→HR impulse responses along with the corresponding P_{1} and P_{2} indexes of the PNS and S_{1} and S_{2} indexes of the SNS before (control) and after vagal, β-sympathetic, and double blockade in the standing posture. Generally speaking, it can be seen visually that the early and fast PNS components of the two impulse responses were blunted much more following vagal and double blockade than β-sympathetic blockade, whereas the delayed and slower SNS components were more strongly attenuated after β-sympathetic and double blockade than vagal blockade.

Figure 5 summarizes group average results (mean ± SE) for the two PNS indexes of P_{1} and P_{2} derived from the estimated ILV→HR and ABP→HR impulse responses and the conventional HF power of HR variability. As expected, the HF power generally performed well in predicting the known effects of the pharmacological autonomic nervous blockade agents on the PNS. In particular, this index was virtually abolished following vagal and double blockade with *P* < 0.05. Furthermore, the index did not change after β-sympathetic blockade in the supine posture. However, counter to expectation, the PNS index did decrease following β-sympathetic blockade in the standing posture with *P* < 0.05. Overall, the P_{1} index was better in predicting the known drug effects on the PNS. That is, this index reduced considerably after vagal and double blockade with *P* < 0.10 and did not change following β-sympathetic blockade in either posture. The P_{2} index provided even better predictive capacity in that it decreased after vagal and double blockade with *P* < 0.01 while remaining unchanged after β-sympathetic blockade. However, the extent of the reductions was not as marked as the other indexes.

Figure 6 summarizes group average results (mean ± SE) for the two SNS indexes of S_{1} and S_{2} derived from the two estimated impulse responses and the conventional LF/HF ratio and normalized LF power of HR variability. As expected, the LF/HF ratio did not perform as well in predicting the known effects of the autonomic nervous blockade drugs on the SNS as the HF power did in predicting the drug effects on the PNS. Whereas the LF/HF ratio did decrease following β-sympathetic blockade with *P* < 0.05 and did not change after vagal blockade, it did not diminish after double blockade. In fact, this index actually increased following double blockade in the standing posture with *P* < 0.10. The LF/HF ratio appeared to be a somewhat better index of sympathovagal balance, as it doubled following vagal blockade, albeit not to any level of statistical significance. Overall, the normalized LF power did not provide an improved ability to predict the drug effects on the SNS. This index did decrease after double blockade with *P* < 0.10 and after β-sympathetic blockade in the standing posture with *P* < 0.001 but did not diminish after β-sympathetic blockade in the supine posture. Furthermore, the extent of the reductions was always very small. In addition, in contrast to expectation, the index increased after vagal blockade in the supine posture but not the standing posture with *P* < 0.05. The S_{1} index likewise did not offer enhanced predictive capacity. This index did substantially decrease after double blockade with *P* < 0.005 but did not change after β-sympathetic blockade in both postures. Moreover, the SNS index actually decreased after vagal blockade in the supine posture with *P* < 0.10. On the other hand, the S_{2} index afforded overall improved ability in predicting the drug effects on the SNS. In particular, this index did not change after vagal blockade and decreased after double blockade in both postures and β-sympathetic blockade in the standing posture with *P* < 0.05. Its only blemish was that the decrease after β-sympathetic blockade in the supine posture was not statistically significant.

## DISCUSSION

Conventional HR power spectral indexes are indisputably of significant physiological and prognostic value but suffer from limitations as measures of the CANS. Most notably, these indexes neither offer an effective marker of the SNS due to its joint mediation with the PNS in the LF band nor afford specific measures of the CANS due to the input contributions to HR (e.g., ILV and ABP). In this study, we aimed to overcome these limitations by deriving indexes of the PNS and SNS via multisignal analysis of cardiorespiratory variability. More specifically, first, we applied parametric system identification to resting beat-to-beat fluctuations in HR, ABP, and ILV so as to estimate the autonomically mediated ILV→HR and ABP→HR transfer functions (see Fig. 2). Then, we separated each estimated transfer function in the time domain into an early and fast PNS component and delayed and sluggish SNS component (see Fig. 1). Finally, we computed scalar indexes of the PNS and SNS from the respective time domain or impulse response components. In this way, the overlap of the SNS and PNS in the frequency domain may be circumvented, and the input contributions to HR may be eliminated. Analogous to the HR power spectral indexes, the derived indexes represent the responses of the CANS to unity changes in their inputs about the mean value and therefore do not reflect basal autonomic tone. Whereas there have been many studies of ILV→HR and ABP→HR couplings (see, e.g., Ref. 39), this study goes one step further by attempting to separately quantify the two CANS from the couplings. We note that this basic idea was initially introduced by us with co-workers in the appendix of Ref. 38. In that study, we specifically proposed the method described herein to separate the ILV→HR impulse response into PNS and SNS indexes (P_{1} and S_{1}). In this study, we introduced a method to decompose the less straightforward ABP→HR impulse response into PNS and SNS indexes (P_{2} and S_{2}). In addition, and significantly, we evaluated, for the very first time, both pairs of indexes, specifically in the context of selective pharmacological autonomic nervous blockade in 14 human subjects. Our results showed that, at the expense of requiring additional measurements, the newly derived indexes, particularly P_{2} and S_{2}, were better than the HF power, LF/HF ratio, and normalized LF power of HR variability in correctly predicting the known drug effects on the PNS and SNS (see Figs. 5 and 6).

In theory, the HF power of HR variability, which is mediated solely by the PNS, should be a useful vagal index, unless there is a significant change in respiratory effort. The experimental results of the present study confirm this theory. In particular, HF power was able to correctly predict the effects of vagal and double blockade in the supine and standing postures and β-sympathetic blockade in the supine posture wherein no significant change in the HF power of ILV fluctuations occurred. However, this PNS index actually decreased after β-sympathetic blockade in the standing posture, likely due to the significant reduction (*P* < 0.05) in the HF power of ILV fluctuations. The change in HF power of ILV fluctuations during only this particular condition may be related to postural factors, bronchoconstriction resulting from β-sympathetic blockade, and/or the random-interval breathing protocol. Regardless of the mechanism, this result suggests that the HF power will be an ineffective marker of the PNS in situations where breathing is altered (e.g., exercise).

As discussed above, a basic idea of the P_{1} and P_{2} indexes is to eliminate the contributions of respiration so as to result in more specific measures of the PNS. Our results showed that both of these indexes were, in fact, able to correctly predict the effects of all studied autonomic nervous blockade conditions. The P_{2} index was a better predictor than the P_{1} index in terms of statistical significance but did not decrease as much after vagal and double blockade.

In theory, without confounding input alterations, the LF/HF ratio of HR variability should generally be a good measure of the SNS when only the SNS changes or the SNS and PNS change in opposite directions and a poor marker when only the PNS changes or both CANS change in the same direction. Indeed, our results showed that this index was able to correctly predict the effects of β-sympathetic blockade but not double blockade on the SNS. The index even showed a tendency to increase after double blockade in the standing posture. In addition, the LF/HF ratio increased by a factor of two after vagal blockade, although not to any level of statistical significance due to the large variance resulting from a division by zero-like effect. Overall, the LF/HF ratio appeared to be a somewhat better marker of sympathovagal balance in the present study, which is consistent with its original intention (22). However, we note that even if it were a perfect index of sympathovagal balance and the HF power were a perfect index of the PNS, they could not, in general, be utilized in tandem to quantify the SNS. For example, if the LF/HF ratio increased while the HF power decreased, this could indicate that the SNS was enhanced, remained the same, or was diminished but not to the extent of the PNS.

The normalized LF power of HR variability differs from the LF/HF ratio in that it includes LF power in the denominator [i.e., LF/(LF+HF) (31)]. In theory, this extra term should improve the ability to predict the effects of double blockade on the SNS (as the term will not go exactly to 0 in practice) and should markedly blunt changes following vagal and β-sympathetic blockade. Consistent with this theory, the normalized LF power was superior to the LF/HF ratio in predicting the effects of double blockade on the SNS and inferior in predicting the effects of β-sympathetic blockade. Moreover, the normalized LF power did not change as much as the LF/HF ratio after vagal blockade. However, the increase in the former index following vagal blockade in the supine posture actually showed statistical significance because of the reduced variance due to the presence of the extra denominator term. Overall, the normalized LF power appeared to be a better index of the SNS than sympathovagal balance in our study.

Like the HF power of HR variability, the predictive capacity of both the LF/HF ratio and normalized LF power can, in theory, be adversely affected by alterations in the inputs to HR. In our study, the LF and HF powers of the ABP and ILV input fluctuations did significantly decrease during some of the experimental conditions. However, the contributions of these input changes to the LF/HF ratio and normalized LF power were generally more difficult to interpret.

As discussed above, the basic idea of the S_{1} and S_{2} indexes is to circumvent the frequency domain overlap of the SNS and PNS, while eliminating input contributions, by separating autonomically mediated transfer functions in the time domain using physiological knowledge and thereby arrive at more effective markers of the SNS. Indeed, the S_{1} index decreased to a larger degree after double blockade and with a greater level of statistical significance than the normalized LF power. However, this SNS index did not decrease following β-sympathetic blockade and actually decreased after vagal blockade in the supine posture. Thus, the S_{1} index did not provide an overall improved measure of the SNS in our study. On the other hand, the S_{2} index did enhance the ability to predict the known drug effects on the SNS. This index correctly decreased following double blockade in both postures and β-sympathetic blockade in the standing posture and did not change after vagal blockade in the two postures. It was only unable to identify a statistically significant decrease after β-sympathetic blockade in the supine posture, perhaps because sympathetic nervous outflow is very low in this posture (26). Note that even the LF/HF ratio only slightly decreased during this experimental condition. We hypothesize that the improved performance of the S_{2} index over the S_{1} index was due to differences in the frequency content of ABP and ILV. In particular, the LF ILV fluctuations may have been insufficient to accurately identify the sluggish SNS component of the ILV→HR impulse response. That is, even though the ILV fluctuations here were induced by a random-interval breathing protocol, the probability density of the respiratory intervals decreased exponentially with interval length (7). In contrast, the naturally occurring LF ABP fluctuations appear to have been enough for more reliable estimation of the SNS component of the ABP→HR impulse response.

The major aim of this study was to highlight the enhanced ability of our multisignal analysis of cardiorespiratory variability over HR power spectral analysis in correctly predicting large changes in the CANS induced by selective pharmacological autonomic nervous blockade. However, our results also revealed the capacity of these indexes to predict postural changes. In particular, all three PNS indexes correctly decreased upon standing in the control condition (see Fig. 5), whereas only the LF/HF ratio and normalized LF power correctly increased upon standing in the control condition (see Fig. 6). This result further suggests that the S_{1} and S_{2} indexes may not be as sensitive to smaller SNS changes. This insensitivity could be partly due to any nonautonomic nervous contributions to the ILV→HR and ABP→HR impulse responses such as mechanical modulation of the sinoatrial node by respiration (9) and variations in baroreceptor and sinoatrial node responsiveness (28). Future attempts to eliminate nonautonomic nervous contributions from the impulse responses [e.g., as proposed by Porta et al. (28)] may increase the sensitivity of the indexes.

Random-interval breathing was employed in this study to enhance the reliability of the estimated ILV→HR and ABP→HR impulse responses. This protocol is not utilized in standard HR variability analysis and could have negatively impacted the performance of the HR power spectral indexes. However, as described above, the results of these indexes reported herein are consistent with expectation based on previous spontaneous or controlled breathing studies. We therefore believe that these results are generally representative of what would have been obtained with conventional breathing protocols. However, we do speculate that some of the detailed aspects of the results may have been affected. For example, conventional breathing may have reduced the variance of the LF/HF ratio during vagal blockade and therefore result in a statistically significant increase in this index following vagal blockade. This result would have made the LF/HF ratio a much better index of sympathovagal balance but an even poorer index of the SNS and thereby better demonstrate the superiority of the S_{2} index as a marker of the SNS.

Our multisignal analysis of cardiorespiratory variability has certain advantages and disadvantages with respect to other recently proposed analyses mentioned above aiming to likewise derive improved indexes of the PNS and SNS. The blind source separation analyses of HR and ABP variability as well as HR and QT interval variability by Vetter et al. (33, 34) are based on the generally invalid assumption that the two CANS operate independently of each other. Our multisignal analysis is theoretically advantageous in that it makes no assumption about the relationship of the PNS and SNS. However, the blind source separation analyses have the practical advantage of requiring fewer measurements. The principal dynamic mode analysis of HR variability by Chon and co-workers (40, 41) is based on a nonlinear expansion. Remarkably, these investigators were able to show that the first two principal dynamic modes of the estimated first- and second-order Volterra kernels (calculated by invoking orthogonality) just happened to correspond to the PNS and SNS in their random-interval breathing datasets. Our multisignal analysis is theoretically advantageous in that it is physiologically based. Furthermore, it potentially has the practical advantage of not being reliant on random-interval breathing, since it is not a nonlinear analysis, which generally has more stringent demands on the frequency content of the inputs (19), employs the parsimonious weighted-principal component regression method (see the appendix), and is most effective with respect to the ABP→HR impulse response (i.e., the P_{2} and S_{2} indexes). On the other hand, the principal dynamic mode analysis has the practical advantage of requiring only a surface ECG measurement and may be theoretically advantageous in terms of extracting nonlinear information from HR variability. Future comparisons of our multisignal analysis with the blind source separation and principal dynamic mode analyses during spontaneous breathing and random-interval breathing are certainly warranted.

Potential applications of our multisignal analysis of cardio-respiratory variability include any research or clinical setting in which continuous ABP, HR, and ILV or at least continuous ABP and surface ECGs are available [as ILV may potentially be estimated from surface ECGs (20)]. An example of such a clinical setting is the intensive care unit wherein HR power spectral indexes, as mentioned above, have been shown to be predictive of patient outcome (36). Further demonstrations of the multisignal analysis in providing improved markers of the PNS and SNS as well as enhanced association with mortality and disease are needed, particularly under spontaneous breathing conditions, to ultimately realize these applications.

## APPENDIX

The weighted-principal component regression method utilized herein to estimate the parameter sets specifying the ILV→HR and ABP→HR impulse responses in *Eq. 1* is described in detail in Ref. 37. We briefly review the major steps of this method below.

First, the AR parameters {*a*[*i*]} in *Eq. 1* are initialized to 0, and the terms *m*, *M*, and *N* are set to values that are believed to encompass the actual durations of the ILV→HR and ABP→HR impulse responses (i.e., FIRs). (In this study, we set *m*, *M*, and *N* to −7, 42, and 50, respectively, thereby assuming that each FIR is no greater than 50 samples in duration and accounting for the noncausality of the ILV→HR impulse response.) Then, the upper model in *Eq. 1* is formed into a matrix equation by stacking each individual equation, corresponding to each discrete time *n*, one on top of the other as follows: **y** = **Xh** + **w**, where **y**, **h**, and **w** are vectors, respectively, comprising the samples of HR, the two FIRs to be estimated ({*g*[*i*]} and {*h*[*i*]}), and *W*_{HR}, whereas **X** is a matrix consisting of samples of delayed versions of ILV and ABP. Next, matrix **X** is postmultiplied by diagonal matrix **W**, whose diagonal elements decay exponentially so as to assume that the current output sample is more dependent on recent input samples than distant input samples. Then, the principal components of matrix **XW** (“weighted-principal components”) are added into a regression structure, one at a time, in order of descending principal components and regressed on to vector **y** until the minimum description length criterion is minimized (17) so as to arrive at an estimate of vector **x** (i.e., the two FIRs). This calculation is repeated for various exponential weighting decay rates, with the one minimizing vector **w** in the least squares sense ultimately selected. Next, *N*_{HR} is calculated (i.e., **y** − **Xh**), and the AR parameters are estimated by standard linear least-squares minimization of *W*_{HR} in the lower model in *Eq. 1* (17). (For simplicity, in this study, *P* was set to 20. This value was justified by observing that *W*_{HR} was indeed generally white via standard autocorrelation analysis.) Finally, HR, ABP, and ILV are prefiltered using the estimated AR parameters, and the above steps are repeated until the estimated FIRs converge. Importantly, this method effectively represents the FIRs, asymptotically, with damped sinusoidal basis functions that reflect the dominant frequency content of the inputs. Thus, only a small number of basis functions (or weighted-principal components) should be needed, thereby decreasing the number of parameters to be estimated and thus the estimation error. (In this study, 20–25 basis functions were nominally used, which represents a >75% reduction in the number of estimated parameters.)

## GRANTS

This work was supported by National Institutes of Health Grants EB-004444 and HL-080568.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2008 by the American Physiological Society