## Abstract

A double exogenous autoregressive (XXAR) causal parametric model was used to estimate the baroreflex gain (α_{XXAR}) from spontaneous R-R interval and systolic arterial pressure (SAP) variabilities in conscious dogs. This model takes into account *1*) effects of current and past SAP variations on the R-R interval (i.e., baroreflex-mediated influences), *2*) specific perturbations affecting R-R interval independently of baroreflex circuit (e.g., rhythmic neural inputs modulating R-R interval independently of SAP at frequencies slower than respiration), and *3*) influences of respiration-related sources acting independently of baroreflex pathway (e.g., rhythmic neural inputs modulating R-R interval independently of SAP at respiratory rate, including the effect of stimulation of low-pressure receptors). Under control conditions, α_{XXAR} = 14.7 ± 7.2 ms/mmHg. It decreases after nitroglycerine infusion and coronary artery occlusion, even though the decrease is significant only after nitroglycerine, and it is completely abolished by total arterial baroreceptor denervation. Moreover, α_{XXAR} is comparable to or significantly smaller than (depending on the experimental condition) the baroreflex gains derived from sequence, power spectrum [at low frequency (LF) and high frequency (HF)], and cross-spectrum (at LF and HF) analyses and from less complex causal parametric models, thus demonstrating that simpler estimates may be biased by the contemporaneous presence of regulatory mechanisms other than baroreflex mechanisms.

- cardiovascular variability
- linear parametric modeling
- identification techniques

evaluation of the baroreflex gain is considered an important tool in clinical practice (16). The baroreflex gain is usually estimated by inducing an increase of arterial pressure by means of a vasoconstrictive drug and by evaluating the subsequent lengthening of the heart period (23). Because this procedure is based on a nonphysiological stimulus that may alter the actual value of the baroreflex gain, several techniques have been proposed to evaluate the baroreflex gain based on the spontaneous variability of systolic arterial pressure (SAP) and heart period (R-R interval). Traditional methods are based on *1*) detection of sequences of SAP and R-R interval values characterized by the simultaneous increase or decrease of both variables (10),*2*) calculation of the power spectrum of R-R interval and SAP variability series (19), and *3*) calculation of the magnitude of the SAP-R-R transfer function (22).

The main limitation of these methods is that they are not capable of evaluating the fraction of R-R interval variability driven by SAP changes because they do not impose causality (i.e., a relationship linking the current R-R interval to past SAP values). Therefore, they cannot reveal whether the observed changes in SAP and R-R interval are the result of feedback mechanisms from SAP to R-R interval (i.e., the baroreflex pathway) or of the feedforward relationship from R-R interval to SAP (mainly related to Starling and windkessel effects). As a consequence, these methods can only produce indexes lumping together the properties of both paths, and denervation experiments or pacing are required to disentangle the feedback pathway from the feedforward pathway (2, 17). Moreover, the effects of other variables on the estimate of the baroreflex gain are not explicitly taken into account, although they are acknowledged. For example, Smyth et al. (23) suggested fixing the respiratory phase in assessing the baroreflex gain. Changes of the R-R interval synchronous with respiration can occur independently of variations in SAP and, therefore, of the arterial baroreflex circuit. For example, changes in venous return and central venous pressure related to respiration can produce activation or deactivation of the low-pressure receptors (25) and, therefore, through the Bainbridge reflex, changes in the heart period (28). Moreover, even effects of centrally mediated rhythmic variations in the autonomic drive to the sinus node (14) produce a certain amount of R-R variability not mediated by the baroreceptive circuit. The presence of a fraction of R-R interval variability independent of SAP changes, if not explicitly addressed, may lead to an overestimate of the baroreflex gain.

The aim of this study is to propose a model-based approach to the evaluation of the baroreflex gain capable of imposing causality (i.e., to evaluate the fraction of R-R variability driven by SAP changes) and directly accounting for sources capable of influencing the R-R interval independently of SAP changes (mainly the effects of slow rhythms and respiration directly driving R-R interval).

Three parametric linear model structures are considered: *1*) the R-R-SAP exogenous model (X model), *2*) the R-R-SAP exogenous model with autoregressive (AR) input (XAR model), and*3*) the R-R-SAP exogenous model with exogenous respiration and AR input (XXAR model). In the X model the R-R interval depends on the current and several previous SAP values. In the XAR model the R-R-SAP relationship is described by an X model, but rhythmic sources affecting R-R interval independently of baroreflex pathway are also modeled. In the XXAR model, in addition to XAR features, the nonbaroreflex effects of respiration on R-R interval are described.

Comparison of the baroreflex gain estimated by the X model with those derived from traditional methods (10, 19, 22) allows us to clarify the role of causality in the estimation of baroreflex gain. Comparison of the baroreflex gains obtained by the proposed parametric models allows us to evaluate the distorting effects of sources driving R-R interval variability independently of SAP changes on the estimate of the baroreflex gain and to propose a model (i.e., the XXAR model) providing a less biased estimate of the baroreflex gain.

The study was carried out on conscious dogs (21) undergoing the following experimental conditions: *1*) control, *2*) nitroglycerine infusion (NT) producing a sympathetic activation by decreasing SAP, *3*) coronary artery occlusion (CAO) increasing sympathetic activity without an important change in mean arterial pressure, and *4*) total arterial baroreceptive denervation (TABD) preventing SAP changes to produce baroreflex-mediated R-R interval variations.

## CAUSAL PARAMETRIC MODELING APPROACH TO ESTIMATE OF BAROREFLEX GAIN

This section describes the three parametric models (X, XAR, and XXAR) capable of extracting the baroreflex gain from spontaneous variability of R-R interval and SAP.

The main characteristics of these models are that *1*) they describe the linear relationship among changes of the variables around a set point (the mean value); *2*) they are dynamic because several past samples are considered; *3*) they are causal because they fix a temporal direction of the influences (R-R interval changes depend on current and past variations of SAP); and*4*) they are defined in the beat-to-beat domain, thus allowing us to reduce the complexity of the model structure (we do not have to reconstruct all the features of the waveform) while maintaining information about the short-term regulatory mechanisms (6,11).

In the discussion of the three models rr_{i}, sap_{i}, and resp_{i} represent the difference between the actual sample RR_{i}, SAP_{i}, and Resp_{i} and their mean values, where RR is R-R interval, Resp is the respiratory signal, and *i* is the progressive cardiac beat number.

#### R-R-SAP X model.

The simplest parametric model describing the rr series and the causal influences of sap on rr can be defined as
Equation 1where *h*
_{rr-sap,k} is the*k*th coefficient setting the influences of sap on rr and*w*
_{rr} is a zero mean white noise with variance λ_{wrr}
^{2}. Therefore, the rr_{i} value depends on the sap (the exogenous signal) at the same cardiac beat (i.e., sap_{i}) and on *p* past values of sap. In this way, causal effects of SAP on R-R interval (i.e., the baroreflex-mediated influences) are explicitly addressed. Fig. 1
*A*represents this relationship in terms of a block diagram. All the rr variability independent of sap is lumped into*w*
_{rr}. The sap series is modeled as an AR process of order *p*
Equation 2where *h*
_{sap-sap,k} is the*k*th coefficient of the AR model and*w*
_{sap} is a zero mean white noise with variance λ_{wsap}
^{2}.

#### R-R-SAP XAR model.

To take into account both direct effects of sap and possible influences (indicated as *u*
_{rr}), including rhythmic influences, affecting rr independently of sap, the model described by*Eq. 1
* can be modified as
Equation 3where sap is described as in *Eq. 2
* and*u*
_{rr} is the AR process
Equation 4where *h*
_{urr-urr,k} is the*k*th coefficient of the AR model. Therefore, rr_{i} depends not only on current and past sap values but also on a source (*u*
_{rr}) describing all the influences unrelated to sap but capable of producing rhythmic changes in rr (Fig. 1
*B*). In this way, oscillations generated by mechanisms involving nonbaroreflex circuits (e.g., cardiogenic reflexes mediated by myocardial ischemia) are explicitly addressed.

#### R-R-SAP XXAR model.

The most complex model is introduced to separate the influences of respiration on rr acting independently of sap. It is defined as
Equation 5where *h*
_{rr-resp,k} is the*k*th coefficient setting the influences of resp on rr, and*u*
_{rr} is described by *Eq. 4
*. The respiratory signal is modeled as an AR process
Equation 6where *w*
_{resp} is zero mean white noise with variance λ_{wresp}
^{2}. Therefore, three inputs are considered to be capable of driving rr changes (Fig. 1
*C*): *1*) current and past sap variations producing R-R interval variability via the baroreflex pathway, *2*) current and past respiratory changes not mediated by sap variations (e.g., oscillations mediated by respiratory centers and by the activation of low-pressure receptors), and*3*) unmeasurable inputs, particularly slow oscillations, independent of sap and resp (e.g., slow rhythms accompanying myocardial ischemia).

#### Evaluation of baroreflex gain.

After the coefficients of the models are identified (see*Identification and validation of X, XAR, and XXAR models* in the
), the block *H*
_{rr-sap}(Fig. 1) is fed by an artificial signal simulating a unitary ramplike rise of SAP and the slope of the relevant increase of R-R interval is taken as an index of the baroreflex gain. The slope of the R-R interval response is calculated with least-squares linear fitting evaluated on the first 15 samples.

#### Hypothesis testing and evaluation of performance of model.

The adequacy of the model in describing the dynamics of the R-R signal and its interactions with SAP and Resp signals is tested by evaluating whether the residual signals (*w*
_{rr} and*w*
_{sap} for X and XAR models and*w*
_{rr}, *w*
_{sap}, and w_{resp} for the XXAR model) are white and not correlated with each other (see *Identification and validation of X, XAR and XXAR models* in the
).

The performance of the model is evaluated by calculating the goodness of fit (ρ) of the R-R interval series. It is defined as
Equation 7where ς_{rr}
^{2} and λ_{wrr}
^{2} represent the variance of R-R interval and of its residue, respectively, thus measuring the percent variance of the R-R interval series explained by the proposed model. ρ_{rr} = 1 means that the R-R interval power is completely described by the model; an insufficient model structure produces ρ_{rr} closer to 0.

## TRADITIONAL APPROACHES TO ESTIMATE OF BAROREFLEX GAIN

This section summarizes the most utilized techniques for estimating the baroreflex gain from spontaneous variability series of R-R interval and SAP. These techniques are based on baroreflex sequence (10), power spectrum (19), and transfer function (22) analysis.

#### Baroreflex sequence analysis.

This method is based on the search for sequences characterized by the contemporaneous increase (positive sequence) or decrease (negative sequence) of R-R interval and SAP. The lag (τ) between R-R interval and SAP samples is chosen as the lag producing the maximum in the normalized R-R-SAP cross-correlation function. Both positive and negative sequences are referred to as baroreflex sequences if they match the following prerequisites: *1*) the total of R-R variations is >5 ms; *2*) the total of SAP variations is >1 mmHg; and *3*) the length of the sequences is four beats (3 variations). For each sequence the slope of the regression line in the plane (SAP_{i}
_{-τ},RR_{i}) is calculated. All the slopes with correlation coefficient >0.85 are averaged, and this average, represented by α_{BS}, has been taken as a measure of the baroreflex gain (10). The reliability of α_{BS} depends on the number of baroreflex sequences detected in the series.

#### Power spectral analysis.

The method relies on the calculation of the power spectrum of the rr and sap series (see *Power Spectrum Estimate* in the
) and on the evaluation of the power inside two bands, the low frequency (LF, from 0.04 to 0.14 Hz) and high frequency (HF, ±0.04 Hz around the respiratory rate) bands (26). Two indexes estimating the baroreflex gain are calculated (19) as
Equation 8where P_{rr(LF)}, P_{sap(LF)}, P_{rr(HF)}, and P_{sap(HF)} represent the power in LF and HF bands detected on R-R interval and SAP series, respectively. These indexes are reliable if the coherence function (*K*
^{2}) sampled at LF and HF is >0.5 (11), thus meaning that the rr and sap signals are significantly correlated at that specific frequency and, therefore, the power calculated in LF and HF bands is not only the effect of independent noises impinging on both signals.

#### Transfer function analysis.

The R-R-SAP transfer function modulus can be calculated from the zero mean sap and rr series as
Equation 9where C_{rr-sap} is the cross spectrum between rr and sap series and P_{sap} is the spectrum of sap series (see*Cross-Spectrum Estimate* in the
). The transfer function modulus can be sampled at LF and HF detected on sap series (22), thus obtaining two indexes of baroreflex gain that are referred as to α_{CS(LF)} and α_{CS(HF)}. This approach also requires us to test whether*K*
^{2} is >0.5 at LF and HF (7, 11).

## METHODS

#### Surgical preparation, experimental protocol, and recorded variables.

We used a subset of experiments carried out on chronically instrumented, conscious dogs to evaluate the effects of several experimental maneuvers on the cardiovascular variability series of R-R interval and SAP (21). The surgical procedure and the experimental protocol were described previously by Rimoldi et al. (21).

Briefly, all the dogs were chronically instrumented to measure electrocardiogram (ECG) and arterial pressure. Under anesthesia, the ECG electrodes were subcutaneously fixed to intercostal muscles (lead II). A catheter with strain-gauge transducers (Statham Instruments, Oxnard, CA) was inserted in the femoral artery and advanced to the abdominal aorta. In addition to ECG and arterial pressure, the respiratory movements were monitored by means of a thoracic belt connected to strain gauges. The respiratory signal was used to extract the breathing rate. In six dogs, a hydraulic occluder was positioned around the left circumflex coronary artery to occlude the artery when inflated with a volume of saline. In four dogs, TABD was performed in two steps. First, after both carotid artery bifurcations were identified, the sinus nerves were located and severed. Second, at the time of the thoracotomy, the adventitia surrounding the aortic arc and its branches was carefully dissected. Completeness of denervation was assessed by the loss of the usual heart rate response to pressure increases mechanically produced by occluding the distal thoracic aorta. Attention was paid to avoid visible damages to the vagi or to the sympathetic nerve branches.

Recordings were performed after a recovery period of 1–2 wk. All the dogs were acquainted with the laboratory and were trained to lie unrestrained on the recording table. Recordings were carried out:*1*) under control conditions (*n* = 16),*2*) 5–10 min after intravenous NT (32 μg · kg^{−1} · min^{−1}) (*n* = 8), *3*) during brief (2 min) CAO (*n* = 6), and *4*) after TABD (*n* = 4). Experimental protocols were approved by the Committee on Animal Use and Care of the University of Milan.

At control conditions, mean R-R and SAP were 739 ± 128 ms and 126 ± 22 mmHg, respectively. After NT, a moderate hypotension (104 ± 13 mmHg) and a reflex tachycardia (539 ± 139 ms) were observed. During CAO, the decrease of R-R interval was even more marked (499 ± 44 ms), without marked changes in mean arterial pressure. After TABD, the R-R was 621 ± 179 ms and SAP strongly increased (167 ± 9 mmHg).

#### Beat-to-beat variability series extraction and baroreflex gain evaluation.

ECG, arterial pressure, and respiratory signals were analog/digital converted at 300 Hz. The QRS peak was detected by using a threshold on the first derivative of ECG. Jitters in the location of the R peak were minimized by parabolic interpolation. The *i*th R-R measure was obtained as the temporal distance between two consecutive R peaks. The maximum of the arterial pressure inside the *i*th R-R interval was the *i*th SAP. The respiratory signal was sampled in correspondence with the first QRS defining the *i*th R-R interval.

Indexes of the baroreflex gain (α) were derived directly from sequences of ∼300 consecutive R-R intervals and SAP values by means of the proposed models (α_{X}, α_{XAR}, and α_{XXAR} from X, XAR, and XXAR models, respectively). We also calculated α_{BS} (10), α_{PS(LF)} and α_{PS(HF)} (19), and α_{CS(LF)} and α_{CS(HF)} (22).

#### Statistical analysis.

Comparisons between control and other experimental conditions were performed by means of one-way analysis of variance (Tukey test). If the normality test was not passed, Kruskal-Wallis one-way analysis of variance on ranks was used (Dunn test).

The index α_{X} was compared with α_{BS}, α_{PS}, and α_{CS} in the same experimental condition by using one-way repeated-measures analysis of variance (Tukey test). If the hypothesis of normality was not fulfilled, Friedman one-way repeated-measures analysis of variance on ranks was used (Tukey test). This comparison is performed to evaluate the extent to which the estimate of the baroreflex gain changed when causality was imposed and variability sources different from SAP were disregarded. The same statistical analysis was used to compare α_{XXAR}to α_{X} and α_{XAR} in the same experimental condition. This comparison was performed to evaluate the extent to which the estimate of the baroreflex gain varied when sources independent of SAP variations but capable of driving R-R interval changes were explicitly disentangled. A *P* < 0.05 was considered significant.

## RESULTS

#### Estimating baroreflex gain from spontaneous beat-to-beat variability of R-R interval and SAP: an example.

An example of beat-to-beat variability series of R-R interval and SAP under control conditions is depicted in Fig.2, *A* and *B*. The power spectra (Fig. 2, *C* and *D*) point out that the R-R interval and SAP variabilities were dominated by a clear HF rhythm (at 0.23 Hz in this example) synchronous with respiration (not shown). The relevant *K*
^{2} (Fig.3) is close to 1 at HF and even at frequencies higher than HF. In contrast,*K*
^{2} is <0.5 at frequencies lower than HF (from 0 to 0.1 in this example). Both the proposed causal parametric approach estimating α_{X}, α_{XAR}, and α_{XXAR} and traditional methods producing α_{BS}, α_{PS}, and α_{CS} were applied to this sequence of data.

In Fig. 4 the response of the*H*
_{rr-sap} block of the X, XAR, and XXAR models (Fig. 1, *A*, *B*, and *C*, respectively) to the unitary ramplike increase of sap is represented. Their least-squares fittings, the slope of which provided the estimate of the baroreflex gain (α_{X}, α_{XAR}, and α_{XXAR}, respectively), were superposed. The X model provided the largest estimate of the baroreflex gain (α_{X} = 48.9 ms/mmHg), whereas the smallest estimate was furnished by the XXAR model (α_{XXAR} = 31.6 ms/mmHg). α_{XAR}, based on the XAR model, was 41.2 ms/mmHg.

α_{BS} is based on the detection of the baroreflex sequences present in the two series. In this set of data they are only 4%. All these sequences, when plotted in the plane (SAP, RR), appeared as straight segments (Fig. 5). These segments were covered from the lower left corner to the upper right corner for the positive sequences and in the opposite direction for the negative sequences. The index α_{BS} = 54.2 ms/mmHg was obtained as the average slope of these segments.

α_{PS} is based on power spectral analysis of R-R interval and SAP series and on the calculation of the power associated with spectral components. The HF components are represented in Fig.6. The power associated with the HF components (Fig. 6, *A* and *B*) was used to evaluate an α_{PS(HF)} = 39.2 ms/mmHg, whereas a consistent estimate of α_{PS(LF)} was prevented by the absence of LF components of R-R interval variability. The reliability of the estimate of α_{PS(HF)} was confirmed by the high degree of squared coherence (Fig. 3) at HF (*K*
_{HF}
^{2} = 0.99).

α_{CS} is based on the estimation of the SAP-R-R transfer function magnitude (Fig. 7). It was sampled at HF detected on the SAP series, thus producing α_{CS(HF)} = 41.9 ms/mmHg. α_{CS(LF)} could not be calculated consistently because*K*
_{LF}
^{2} < 0.5 (Fig. 3).*K*
_{HF}
^{2} > 0.5 confirmed the reliability of α_{CS(HF)}.

#### Testing reliability of baroreflex gain estimates.

The reliability of the baroreflex gain estimates based on the causal parametric approach depends on the ability of the model to produce white and noncorrelated residual signals. In the X model,*w*
_{sap} was always white, whereas*w*
_{rr} passed the whiteness test in 14 of 16 animals under control conditions, in 1 of 8 after NT, in 4 of 6 during CAO, and in 2 of 4 in TABD. These results pointed out that the structure of this model is too simple to describe the R-R interval dynamics. However, the two residual signals *w*
_{rr}and *w*
_{sap} in the X model were not correlated, even at zero lag, in any condition. The inability of the X model to describe the R-R interval dynamics is confirmed by the relatively low value of the goodness of fit ρ_{rr} (0.59 ± 0.16 at control, 0.59 ± 0.16 after NT, 0.45 ± 0.14 during CAO, and 0.57 ± 0.19 after TABD). When the structure of the model was rendered more complex (XAR and XXAR models), the tests on the whiteness of the residual signals and on their noncorrelation were fulfilled in all conditions. ρ_{rr} was increased accordingly in the XAR model (0.69 ± 0.13 at control, 0.71 ± 0.17 after NT, 0.75 ± 0.06 during CAO, and 0.76 ± 0.17 after TABD) and revealed the largest values in the XXAR model (0.78 ± 0.12 at control, 0.73 ± 0.16 after NT, 0.79 ± 0.06 during CAO and 0.81 ± 0.12 after TABD).

The reliability of α_{BS} depends of the number of baroreflex sequences detected in the variability series. Under control conditions the mean percentage of baroreflex sequences (over 300 samples) was 5%. In 8 of 16 dogs the percentage was <3%, and in 4 animals no slope was present. An example of R-R interval and SAP beat-to-beat variability in which at least three contemporaneous increases or decreases of R-R interval and SAP were not present (no baroreflex sequence could be detected) is depicted in Fig.8, *A* and *B*. The magnification of two periods of dominant HF rhythm (Fig. 8
*C*) shows that, whereas three consecutive increases can be found in SAP series, R-R interval can increase or decrease only for two beats, thus producing the absence of baroreflex sequences. The percentage of baroreflex sequences was larger after NT (18%; only 1 dog exhibited <3%), during CAO, and after TABD (12% and 24%; no animal exhibited <3%).

The reliability of α_{PS(LF)}, α_{CS(LF)}, α_{PS(HF)}, and α_{CS(HF)} was tested by evaluating the value of the squared coherence at LF and HF (*K*
_{LF}
^{2} and*K*
_{HF}
^{2}). α_{PS(LF)} and α_{CS(LF)} could not be calculated in all animals because LF oscillations were not always found. Moreover, even when LF oscillations were observed, they exhibited a value of*K*
_{LF}
^{2} > 0.5 only in a few cases (7 of 16 at control, 2 of 8 after NT, 5 of 6 during CAO, and 0 of 4 after TABD). In contrast, *K*
_{HF}
^{2} was >0.5 in all animals in all the conditions (0.98 ± 0.03 at control, 0.95 ± 0.05 after NT, 0.91 ± 0.16 after CAO, and 0.98 ± 0.01 after TABD).

#### Causal modeling approach vs. traditional methods.

The results are summarized in Table 1. At control, α_{BS}, α_{PS(HF)}, and α_{CS(HF)} were very high (>40 ms/mmHg). α_{PS(LF)} and α_{CS(LF)} were smaller and similar to α_{X} and α_{XAR}. α_{XXAR}showed the smallest value, significantly smaller than α_{X}and α_{XAR}.

All indexes diminished after NT, with α_{PS(LF)}, α_{CS(LF)}, α_{X}, α_{XAR}, and α_{XXAR} comparable and smaller than α_{BS}, α_{PS(HF)}, and α_{CS(HF)}. During CAO, the indexes α_{BS}, α_{PS(HF)}, α_{CS(HF)}, and α_{X} decreased significantly whereas α_{PS(LF)}, α_{CS(LF)}, and indexes based on more complex models (α_{XAR} and α_{XXAR}) exhibited a nonsignificant fall. It is worth noting that, during CAO, α_{XAR} and α_{XXAR} were comparable but significantly smaller than α_{PS(LF)} and α_{CS(LF)}. After TABD, all the indexes showed a drastic reduction. However, only the indexes based on a causal parametric approach (α_{X}, α_{XAR}, and α_{XXAR}) were close to zero. After TABD, α_{PS(LF)} and α_{CS(LF)} could not be calculated because of the lack of coherence in all animals.

## DISCUSSION

All the methods considered here for estimating the baroreflex gain from spontaneous R-R interval and SAP variabilities are able to detect the decrease of the baroreflex sensitivity during baroreceptive unloading provoked by NT, during sympathetic activation induced by CAO, and after TABD, thus confirming the capability of these parameters to measure the baroreflex gain (Table 1). However, Table 1 points out that the values of the estimated baroreflex gains can be very different. These differences are the result of the different capabilities of the various methods to separate the baroreflex pathway from mechanisms capable of producing R-R interval and SAP variability independently of baroreflex circuit (Table 2). Three main mechanisms are responsible for R-R interval and SAP variabilities independent of baroreflex pathway: *1*) Starling and windkessel effects producing SAP variations driven by R-R interval changes (i.e., the feedforward mechanisms), *2*) neural modulations at HF capable of driving the R-R interval independently of SAP changes [e.g., neural influences projecting the activity of central respiratory oscillators (14) and neural reflexes involving low-pressure receptors activated by respiration-related changes in blood volume and central venous pressure (13,28)], and *3*) neural modulations affecting the sinus node at LF independently of SAP variations [e.g., neural influences projecting the activity of central slow oscillators (20) and cardiogenic reflexes initiated by ischemia]. Also, mechanical influences of respiration on sinus node produce R-R interval changes independent of SAP variations, but they are usually irrelevant [they become apparent in the denervated heart (9)].

#### Traditional methods for estimation of baroreflex gain.

The indexes α_{PS} and α_{CS} (19,22) are calculated without separating all these mechanisms from the baroreflex pathway (Table 2). Indeed, they are estimated without imposing any specific causal direction in the interactions between R-R interval and SAP variabilities (i.e., the feedforward mechanical path from R-R interval to SAP is merged with the feedback baroreflex pathway from SAP to R-R interval). In contrast, the index α_{BS}(10) is actually calculated by separating the baroreflex sequences from the nonbaroreflex sequences, but no memory is allowed in feedback and feedforward relationships, thus preventing the full separation of feedback from feedforward influences (Table 2). As a consequence, the traditional indexes α_{BS}, α_{PS}, and α_{CS} mix the gains of both paths and are closer to the baroreflex gain only if the R-R interval is mainly driven by SAP variations. This condition is more likely to be fulfilled at LF (2, 11). In contrast, R-R interval variability leads SAP oscillations at HF (2, 27), thus rendering α_{PS(HF)} and α_{CS(HF)} sensible to the mechanical pathway (i.e., the Starling and windkessel effects) and introducing a bias when these indexes are used to evaluate the baroreflex gain. According to this consideration, α_{PS(HF)}and α_{CS(HF)} can be different from α_{PS(LF)}and α_{CS(LF)}, respectively, and can be significantly larger than 0 even after TABD (Table 1). Because the index α_{BS} is not conceptually different from α_{PS}and α_{CS}, it could be more similar to α_{PS(LF)} and α_{CS(LF)} or to α_{PS(HF)} and α_{CS(HF)}, depending on the predominant dynamics in the R-R interval and SAP variabilities. In dogs, as a result of a dominant HF dynamic in every experimental condition, α_{BS} is similar to α_{PS(HF)} and α_{CS(HF)} (Table 1). These observations suggest the use of α_{PS(LF)} and α_{CS(LF)} instead of α_{BS}, α_{PS(HF)}, and α_{CS(HF)} in dogs and in other experimental preparations in which respiratory influences might directly drive the sinus node. Unfortunately, the presence of a small amount of R-R power in the LF band in dogs (21) makes it difficult to obtain a robust estimation of α_{PS}(LF) and α_{CS}(LF) (*K*
_{LF}
^{2} is usually <0.5), thus rendering the estimate of baroreflex gain unreliable even in the LF band and necessitating new tools based on a causal modeling approach.

#### Role of causality in estimate of baroreflex gain.

The proposed modeling approach fixes the temporal direction of the influences of SAP on R-R interval (i.e., causality), thus exploring the baroreflex pathway directly. Indeed, the R-R interval depends on SAP inside the same cardiac beat (fast effect) and on past SAP values (dynamic effect). The fast effect has been ascribed to the action of vagal circuits and the dynamic, slow effect to the sympathetic branch (12). This feature gives the model the possibility of separating the baroreflex pathway (the feedback path) from the mechanical relationship (mainly windkessel and Starling effects, the feedforward path) (Table 2). In other words, only that part of the R-R variability causally related to SAP changes is exploited to evaluate the baroreflex gain. As a result, α_{X}, α_{XAR}, and α_{XXAR} are smaller than α_{BS}, α_{PS(HF)}, and α_{CS(HF)} in all experimental conditions, comparable to α_{PS(LF)} and α_{CS(LF)} under control conditions and after NT but significantly smaller during CAO and close to 0 after TABD.

#### Role of rhythmic influences independent of sap changes in estimate of baroreflex gain

The simplest model (i.e., the X model) is able to evaluate the amount of R-R interval variability driven by SAP changes (i.e., the baroreflex-mediated R-R variability) but cannot model rhythmic inputs capable of driving R-R interval independently of SAP variations (Table2). Therefore, α_{X} is calculated without separating the baroreflex pathway from mechanisms [e.g., central slow oscillators (20), central respiratory oscillators (14), and cardiopulmonary reflexes (13, 28)] producing LF and HF oscillations on R-R interval not mediated by the baroreflex circuit (Table 2). In contrast, the XAR and XXAR models can disentangle the baroreflex pathway from mechanisms driving R-R interval independently of SAP (Table 2). The XXAR model can be considered a refinement of the XAR model (Table 2). Indeed, among the rhythmic sources driving R-R interval independently of SAP the XXAR model distinguishes those driven by respiration and, therefore, in the HF band, from those independent of respiration (i.e., in the LF band).

Under control conditions, α_{XXAR} is significantly smaller than α_{X} and α_{XAR}, as a result of the effect of taking into account the respiration-related changes of R-R interval independent of SAP variations, which are particularly prominent in dogs (2, 4). During NT, α_{X}, α_{XAR}, and α_{XXAR} are similar, as a result of a more important role of the baroreflex circuit in regulating R-R interval than that of central respiratory oscillators and low-pressure areas. During CAO, α_{X} is significantly smaller than α_{XAR} and α_{XXAR} because of the technical inability of the X model to fit the data. The use of XAR and XXAR models allows us to solve the identification problem correctly and to find out that α_{XAR} and α_{XXAR} are similar. During CAO, the use of an XAR or XXAR model is necessary to take into account LF oscillations directly impinging on the sinus node as previously reported (4). In this experimental condition, XAR and XXAR models perform similarly because of the reduced importance of HF direct effects on the sinus node. It is worth noting that α_{PS(LF)} and α_{CS(LF)} remain significantly larger than α_{XAR} and α_{XXAR} as a result of the influence of LF fluctuations on R-R interval independent of SAP, thus pointing out that even in the LF band traditional indexes may be biased and making obvious the improvement of the causal model-based approach.

In conclusion, our data suggest that α_{XXAR} provides a less biased value of baroreflex gain than any other indexes derived from sequence, power spectrum, and transfer function analyses or from simpler causal models. α_{XXAR} is calculated by exploiting the R-R variability driven by SAP changes after disentangling different variability sources capable of producing changes in the R-R interval variability independent of the baroreflex circuit. If R-R variability was completely driven by SAP changes, traditional indexes would be equal to α_{XXAR}. However, the former indexes may be larger than the latter because of the bias of direct effects of respiration on R-R variability, of slow fluctuations of R-R variability independent of the baroreflex circuit, and of the mechanical feedforward pathway. The estimation of α_{XXAR} requires the additional recording of a respiratory signal but, in contrast to the sequence analysis, it uses all the information content of R-R interval series and, unlike spectral and transfer function analyses, it does not require the setting of specific frequency bands.

## Appendix

In this appendix we summarize the methods used to estimate the coefficients of the X, XAR, and XXAR models, the power spectrum and the cross spectrum.

#### Identification and validation of X, XAR, and XXAR models.

The coefficients of the X model (i.e., *Eq. 1
*) and of the separate AR models describing the dynamics of sap and resp signals (i.e., *Eqs. 2
* and *
6
*) are identified via a least-squares procedure (8, 15). The coefficients of the XAR model (i.e., *Eqs.* 3 and *
4
*) and of the XXAR model (i.e., *Eqs. 5
* and *
4
*) are identified using a generalized least-squares approach (8, 24). This iterative procedure is stopped when the current iteration does not produce a significant percent decrease in the variance of the residue*w*
_{rr} with respect to the previous one (the threshold is 0.001). The solution of both the traditional and generalized least-squares problems is performed by means of the Cholesky decomposition method (15).

The model adequately describes the dynamics of the signal if the residual signals (*w*
_{rr} and*w*
_{sap} for X and XAR models and*w*
_{rr}, *w*
_{sap}, and*w*
_{resp} for the XXAR model) are white (all the information is captured by the model parameters). The whiteness of the residual signals is verified by the Anderson test (15). This test checks that the normalized autocorrelation functions ρ_{wrr}(τ), ρ_{wsap}(τ), and ρ_{wresp}(τ) (the autocorrelation functions divided by λ_{wrr}
^{2}, λ_{wsap}
^{2}, and λ_{wresp}
^{2}, respectively) are equal to 0 for τ≠0 and τ≤40. If this is true for 95% of the values of τ, the test is fulfilled with 5% confidence, [e.g., with τ_{max} = 40, ρ(τ) ≠ 0 is allowed for <3 values of τ]. The structure of the model is adequate to describe interactions among the signals if the residues are not correlated with each other (all relationships among the series are explained by the model structure). The same test used to check the whiteness of the residual signals is carried out to verify their noncorrelation; it checks that the normalized cross-correlation functions ρ_{wrr}
_{-wsap}(τ) and ρ_{wrr}
_{-wresp}(τ) (the cross-correlation functions divided by the product of the standard deviation λ_{wrr}λ_{wsap}and λ_{wrr}λ_{wresp}, respectively) are 0 for τ ≤ 40 (even at τ = 0). Obviously, in the case of X and XAR models, it is sufficient to test the whiteness of *w*
_{rr} and*w*
_{sap} and their noncorrelation. The choice of a high model order can help to fulfill these hypotheses. However, this solution should be avoided because the fitting is performed even on the noise superposed on the data. To favor small model orders, the model order is selected according to the minimum of the Akaike figure of merit for multivariate processes (1). The best model order is searched in the range from 6 to 16. After choosing the best model order, whiteness and noncorrelation on the residual signals are tested and the goodness of fit ρ_{rr} (i.e., *Eq. 7
*) is calculated.

#### Power spectrum estimate.

The power spectrum is calculated by using a parametric approach (18) instead of a nonparametric approach (3). The rr and sap series are described as AR processes considering the current value of the series as a linear combination of its *p*past values (e.g., see *Eq. 2
* for the AR model of sap). The coefficients of this linear combination are estimated via Levinson-Durbin recursion (15), and the number of the coefficients is chosen according to the Akaike criterion (1). The Anderson test is used to test whether the residue is white (15). After the power spectral densities of these two series are calculated, the power spectral decomposition procedure (29) is carried out to evaluate the contribution of each oscillation (described by a real pole or a pair of complex and conjugated poles) to the total power (i.e., the variance) of the process, thus allowing us to estimate the power in LF and HF bands.

#### Cross-spectrum estimate.

The calculation of the cross spectrum between rr and sap series requires a bivariate approach instead of the monovariate approach required by spectral analysis. We choose a parametric approach based on a bivariate AR model (5) to estimate the autospectrum of sap series and the cross spectrum between sap and rr instead of a nonparametric method (11). The model order is fixed at 10, and the coefficients of the bivariate AR model are identified via least-squares methods (8, 15).

## Footnotes

Address for reprint requests and other correspondence: A. Porta, Universitá degli Studi di Milano, Dipartimento di Scienze Precliniche, LITA di Vialba, Via G.B. Grassi 74 20157 Milano, Italy (E-mail: alberto.porta{at}unimi.it).

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 © 2000 the American Physiological Society