METHODS

General Procedure

The experimental protocol was approved by the Tohoku University Committee on Animal Care. Seven adult mongrel dogs of both sexes, weighing 15–35 kg, were anesthetized by intravenous thiopental sodium (2.5 mg/kg) and ketamine sodium (5.0 mg/kg) administration and nitrous oxide inhalation under mechanical ventilation. Experimental preparations are shown in Fig. 1. After tracheal tube intubation, the animals were placed on a volume-limited respirator (ARF-850, Acoma, Tokyo, Japan). Electrodes for electrocardiograms were implanted in the left foreleg and both hind legs.

After the left pleural cavity was opened through the fifth intercostal space, the pericardium was partially incised to expose the left atrium. Aortic pressure and left atrial pressure (LAP) were monitored continuously by catheters inserted into the aorta and left atrium through the left femoral artery and the left appendage, respectively. For left artificial heart implantation, the intercostal arteries were dissected to free the descending aorta. A polyvinyl chloride outflow cannula with a T-shaped pipe was inserted into the descending aorta and secured with a ligature. An inflow cannula was inserted into the left atrium through the left atrial appendage. Both cannulas were connected to the authors’ TH-7B pneumatically driven sac-type blood pump via the built-in valve connectors (37, 38). An outflow cannula was inserted into the pulmonary artery, and an inflow cannula was inserted into the right atrium through the right atrial appendage. Both cannulas were connected to the right pump. After the driving of both pumps was initiated, ventricular fibrillation was electrically induced. After ventricular fibrillation was induced, we observed that the pulse-synchronous discharges in the RSNA did not change their periodicity and quantity significantly. With respect to the influence of the fibrillated heart on RSNA, Toda et al. (34) reported that the fibrillated heart had little influence on the RSNA and was presumably not essential for maintaining nervous control of circulation as long as the circulation was maintained by an artificial heart.

The pumps were driven by a pneumatic drive console that we developed (24). The positive and negative pressure and systolic duration of both blood pumps were selected to maintain hemodynamic derivatives within the control range: mean aortic pressure was maintained within 60–140 mmHg, mean LAP within 5–12 mmHg, mean pulmonary artery pressure <30 mmHg, and right atrial pressure <5 mmHg (Fig. 1). Pump output was controlled to maintain a cardiac output similar to that during the control state.

Recording of Nerve Activity

The left flank was opened between the iliac crest and the costovertebral angle, and the left renal artery was then exposed. The left renal sympathetic nerve bundle was separated from the left renal artery and surrounding connective tissues. After the nerve sheath was removed, the nerve was placed on bipolar stainless steel electrodes for recording of RSNA. The renal sympathetic nerve discharges were led into a preamplifier with a gain of 1,000 and an amplifier with a gain of 50 (DPA-21 and DPA-11E, DIA Medical System, Tokyo, Japan). Amplified signals were routed through a band-pass filter with a bandwidth of 30–3,000 Hz and through an amplitude discriminator to a storage oscilloscope. For analysis, the filtered neurogram was integrated by a resistance-capacitance circuit with a time constant of 0.1 s.

To identify the observed signal as the sympathetic nerve discharge, we transiently clamped the descending aorta above the level of the renal arteries, as indicated in Fig. 1. During the clamp RSNA decreased rapidly because BP increased at the aortic arch (Fig.2). However, the reason why BP decreased abruptly in Fig. 2 is that a catheter to measure BP was left below the clamping site. We confirmed that the respiratory rhythm of RSNA was changed significantly in accordance with the changes in the respiratory rhythm of the mechanical ventilator (14). Figure3 shows the recording from the natural heart state to a stationary state of biventricular bypass circulation after electrical induction of ventricular fibrillation. Careful inspection of Fig. 3 suggests a reciprocal relationship between RSNA and BP; that is, RSNA bursts occurred primarily when BP was falling, whereas there was little RSNA when BP was increasing. When the nerve discharge signal was significantly contaminated by environmental electromagnetic noise, the animal was covered with wire netting to avoid production of the noise.

Data Analysis

The signals of hemodynamic parameters and sympathetic nerve activity were recorded with an ink-jet recorder and on magnetic tape after the stabilization of all hemodynamic derivatives (almost 20–30 min after the preparation). After the control data were recorded, ventricular fibrillation was induced electrically. Time series were recorded in a stationary state of the hemodynamics during the different pacing rates from 60 to 160 beats/min. These data were played back from the magnetic tape.

The occurrence of R waves was detected in each record, and a smoothed instantaneous HR time series was constructed and sampled at 8 Hz using the algorithm proposed by Berger et al. (8). All R-R intervals during control were measured at an accuracy of 1 ms. A time series of R-R intervals that comprised only beats in normal sinus rhythms and in a stationary state was selected for the final analysis.

A time series of BP was constructed as that of arterial blood pressure averaged from beat to beat. Time series of BP and RSNA were splined and sampled at 8 Hz so that values of all the constructed time series occurred simultaneously.

MI Analysis

We calculated MI values according to an algorithm proposed by Fraser and Swinney (15). For a couple of time series,x(t) andy(t), we measured the dependency of the values ofy(t+ T) on the values ofx(t), where T was a time delay and data length was the power of two. We made the assignment (s,q) = [x(t),y(t +T)] to consider a general coupled system (S,Q). MI of this system [I(S,Q)] is defined as the answer to the question, “Given a measurement ofs, how many bits on the average can be predicted about q?”

We denote time series of HR, RSNA, and BP as HR(t), RSNA(t), and BP(t), respectively. The data length was between 2⁹ and 2¹¹ samples, that is, between 64 and 256 s, because the data were sampled at 8 Hz. IfS =Q, the correlation between them should be perfect, and thenI(S,Q) =n, where the data length is 2ⁿ, because the algorithm is developed to the discrete case. This point will be discussed inappendix . The MI value between the same two time series is n. Hence, MI values were normalized by n; in other words, these values were divided by n. We calculated mutual informationI(T) of (S,Q), whereS is a time series of HR(t), RSNA(t), or BP(t) andQ is another time series with a time delay T; for example,S = HR(t) andQ = RSNA(t +T).T was then between −5 and 5 s in a step of 0.125 s. The maximum value ofI(T) between S andQ, whereT was from −5 to 5 s, is denoted byI_max(S,Q). We considered a time delayT_max(S,Q), at which the maximum value ofI(T) of (S,Q) was given, as a physiological delay between these parameters. The definition ofT_max(S,Q) will be discussed later (seediscussion). IfT_max(S,Q) is negative, thenS leadsQ, and ifT_max(S,Q) is positive, vice versa.

Coherency Analysis

To compare MI analyses, we calculated coherency among HR, RSNA, and BP. Ordinary coherence function estimates were computed as the ratio of the magnitude square of the cross spectra divided by the product of the autospectra. Coherence is a measure of the statistical link between two variability series at any given frequency and is expressed as a number between zero and one. The significant value of coherency will be discussed in Simulation 4. The confidence in spectral estimates can be enhanced by dividing data into multiple epochs and ensemble averaging. Although this decreases variance and increases confidence, it reduces resolution because a reciprocal of the length of the data epoch determines the limits of the low-frequency resolution. Auto- and cross-spectra estimates were optimized by ensemble averaging between 9 and 31 data epochs of 128 points (16 s) using Welch’s method (7). To reduce the loss of stability, the original data were sectioned using a 50% overlap. Before a fast Fourier algorithm was used, linear trends were removed from the data, and data were tapered with the use of a Hanning window.

Statistical Procedures

Values of continuous variables are presented as means ± SE. We used nonparametric two-way ANOVA (Friedman’s test) to test for statistical significance between the various stages of the experimental study. A two-tailed P < 0.05 was considered significant.

RESULTS

Simulation Studies

Simulation 1.

To demonstrate an effect of noise added to original simulated data on MI values between the original time series and its noisy series, a couple of time series were generated with noise according to the following formulas (29)

I[x(t),y(t)] were computed from an original torus series and noisy torus series that were 1,024 samples long and contained 10–90% noise. For each percentage of noise, 10 sets ofI[x(t),y(t)] were computed by generating different seeds of random numbers 10 times.I[x(t),y(t)] at 10, 20, ..., 90% noise were compared by one-way factorial ANOVA. Figure 4 shows that the MI value decreased with an increasing percentage of noise (P < 0.0001). For frequency analysis, these data were sectioned using a 50% overlap, and auto- and cross-spectra estimates were optimized by ensemble averaging 15 data epochs of 128 points. Because the ratiof₁/f₂of the two fundamental frequencies (f₁ = 1.59 Hz andf₂ = 1.43 Hz) is rational, all the peaks are harmonics of the frequencyf =f₁ −f₂ (0.16 Hz) in the power spectrum ofx(t). Because the three main peaks in the spectrum ofx(t) became obscure in the spectrum of noisy torusy(t) containing >80% noise, the average ofI[x(t),y(t)] in whichy(t) contained 70% noise (0.047 ± 0.002) was taken as the threshold value to discriminate correlated from noncorrelated data. In particular, an MI value >0.088, which was the average ofI[x(t),y(t)] wheny(t) contained 50% noise, indicated the strong correlation between them. This result would provide a measure to assess the effect of noise on MI between an original series and the noisy series. With an increasing percentage of noise, the main peaks in the spectrum of the noisy torusy(t) that were consistent with those in the spectrum ofx(t) became obscure and the MI value between them decreased. This suggests that the high MI values corresponded to the fluctuations pertaining to these main peaks. Hence, in cases in which two time series show several consistent peaks in their spectra and have a high MI value, the respective fluctuations pertaining to these main peaks would be strongly correlated with each other between these time series.

In appendix , we describe two simulations. One simulation shows that an MI value of zero means that two time series are neither linearly nor nonlinearly correlated at all. The other simulation shows the reliability of MI analysis for detecting the time delay between two torus series.

Simulation 2.

Suppose thatx(t) andy(t) are both formed from the same purely random processz(t), which has a mean of zero and a variance of${\varsigma }_z^2$, by

Simulation 3.

We examined another case in whichx(t) andy(t) were not linearly correlated at all. The time seriesx(t) was a noise process as shown in Fig. 6, where 1 ≤x(t) ≤ 16,384; if t₁≠ t₂, thenx(t₁) ≠x(t₂). The time seriesy(t) was generated by one-to-one relation from theX-Yplot of Fig. 6. The sampling rate was 1 Hz. The region of theX-Yplot is divided into 16 subregions. In each of the four subregions denoted as A,B, C, and D (Fig. 6),x(t) andy(t) were perfectly linearly correlated. The data length of these series was 16,384 samples. The coherency spectrum was computed in the same manner as in Fig. 5. Coherency ofx(t) andy(t) at any frequency was <0.05. Although this indicates thatx(t) andy(t) were not linearly correlated at all, it provides no information as to whether they were nonlinearly correlated. On the other hand, the MI value between the series was just unity, indicating thatx(t) andy(t) had a perfect correlation. This case demonstrates that the coherency method fails to detect correlation between two time series that are nonlinearly correlated.

Simulation 4.

Using simulated data, we determined that the confidence of coherence depends on the number of data epochs. In the same manner as that shown in Fig. 6,x(t) andy(t) were not linearly correlated at all. The time seriesx(t) was a noise process for which 1 ≤x(t) ≤ 1,024; ift₁ ≠t₂, thenx(t₁) ≠x(t₂). The time seriesy(t) was generated by one-to-one relation from theX-Yplot of Fig. 7. The data length of these series was 1,024 samples. The sampling rate was 1 Hz. For coherency analysis, these data were sectioned without any overlap, and auto- and cross-spectra estimates were optimized by ensemble averaging four data epochs of 256 points. For example, assume that a value ofC(ω) = 0.5 is estimated at a frequency of interest. The normalized random error of the coherence function estimate is then approximated by using the estimateC(ω) in place of the unknown true coherence as follows (Ref. 7, p. 311)

Experimental Data

Figure 8 shows the representative tachograms and power spectra of HR, RSNA, and BP during the natural heart state in a dog. The coherency spectra and phase-angle spectra between each two variables are also shown in Fig. 8. The frequency of the maximum peak power in these power spectra is almost the same (0.12 Hz). The coherency values between HR and BP and between RSNA and BP were >0.5 at 0.12 Hz, indicating that they are significantly linearly correlated with each other, whereas coherency between HR and RSNA was not significant. In the time-domain analysis, the correlation coefficient between HR and BP, that between RSNA and BP, and that between HR and RSNA were −0.17, 0.25, and −0.19, respectively. Although these correlation coefficients were small, the hypothesis of zero correlation for each pair was rejected (Ref. 7, p. 101). On the other hand, the maximum MI values between these variables were I_max(HR, RSNA) = 0.25,I_max(HR, BP) = 0.17, andI_max(RSNA, BP) = 0.13. These values indicate that these variables have mutual dependence on one another and, especially, that HR and RSNA have a very strong correlation. From the cross-spectrum phase-angle estimates at 0.12 Hz, RSNA led HR by 2.9 ± 0.8 s, BP led HR by 2.1 ± 0.3 s, and RSNA led BP by 0.9 ± 0.6 s. These standard deviations of time delays at the frequency were estimated with the use of a formula to compute those of the phase angles (Ref. 7, p. 301). From MI analysis, RSNA led HR by 2.1 s, BP led HR by 2.0 s, and RSNA led BP by 0.5 s. The time delay between each pair of these variables as estimated by MI analysis was within its range as estimated by phase-angle analysis, and the relationship between them with respect to leading was the same.

Figures 9 and10 show tachograms, power spectra of RSNA and BP, and coherency spectra between RSNA and BP during the artificial heart states with pacing rates of 60 and 160 beats/min, respectively, in the same dog as in Fig. 8. With a pacing rate of 60 beats/min, the low-frequency peak in the RSNA spectrum was consistent with that in the BP spectrum. The considerable components appeared in a broader frequency band of spectra of RSNA and BP for a pacing rate of 160 beats/min compared with those for a pacing rate of 60 beats/min. The coherency value between RSNA and BP at 0.12 Hz was higher with a pacing rate of 60 beats/min than with a pacing rate of 160 beats/min. The correlation coefficient between RSNA and BP was higher with a pacing rate of 60 beats/min than with a pacing rate of 160 beats/min (−0.58 and 0.32, respectively). These findings indicate that the linear correlation between RSNA and BP decreased at a higher pacing rate. Because HR was constant during the artificial heart state, only I(RSNA, BP) was calculated.I_max(RSNA, BP) at pacing rates of 60 beats/min and 160 beats/min were 0.16 and 0.08, respectively. This indicates that the correlation between RSNA and BP decreased with an increasing pacing rate. RSNA and BP, however, were significantly correlated even with a pacing rate of 160 beats/min. At a pacing rate of 60 beats/min, RSNA led BP by 3.6 ± 0.2 s according to the cross-spectrum phase-angle estimates at 0.12 Hz, whereas RSNA led BP by 0.6 s according to MI analysis. At a pacing rate of 160 beats/min, the coherencies between RSNA and BP at some frequencies >0.25 Hz were significant and the phase angles at the respective frequencies were roughly zero radian. These results were compatible with those from MI analysis in that they were strongly correlated and the time delays between them were <0.25 s.

During the control state,I_max(HR, RSNA) > I_max(HR, BP) > I_max(RSNA, BP) as shown by ANOVA (P < 0.05) (Fig. 11). This indicates that the correlation between HR and RSNA is stronger than that between HR and BP and that between RSNA and BP.I_max(RSNA, BP) was significantly larger at a pacing rate of 60 beats/min than at the other pacing rates (80–160 beats/min) or in control (P < 0.05) (Fig.12).

The time delay of HR and RSNA [T_max(HR, RSNA)] was significantly different fromT_max(HR, BP) andT_max(RSNA, BP) during the control state (P < 0.05) (Fig. 13). These results show that RSNA leads BP, BP leads HR, and RSNA leads HR during the control state. The values ofT_max(RSNA, BP) were significantly different at the pacing rates (60–160 beats/min) or in control (P < 0.05) (Fig. 14).

DISCUSSION

In this study we proposed and tested the MI technique to assess correlation between two time series. In the case in which the MI value between a pair of variables is >0.088, which is the average ofI[x(t),y(t)] wherey(t) contains 50% noise, it is indicated that these variables are correlated even if they are contaminated by noise. On the other hand, some researchers took a coherency value of 0.5 as the threshold value above which there is a significant correlation between oscillations in different variables (28, 31). Although the confidence of coherence can be enhanced by dividing data into multiple epochs and ensemble averaging, it might be practically difficult to achieve the normalized random error, which equals 0.1. For example, when a value of coherence equals 0.5, >100 epochs are necessary to achieve ∼95% confidence intervals (7). This reduces resolution because the length of the data epoch determines the limits of the low-frequency resolution. However, the coherency analysis might have the advantage of gauging the linear correlation between two time series at any frequency. Hence, to quantify the relationship between two time seriesx(t) andy(t), coherency analysis as well as MI analysis should be used. Although the MI value does not quantify the relationship between the time series at a frequency of interest, it can quantify the nonlinear as well as linear dependence between them. In particular, the findings inSimulation 1 indicate that, in the case of two time series showing several consistent peaks in their spectra and a high MI value, the respective fluctuations pertaining to these main peaks would be strongly correlated with each other between these time series.

During the control state, the frequency of the maximum power in the HR spectrum was the same as that of the maximum power in the RSNA spectrum. Moreover, the MI value between them was 0.25, indicating that they had a strong correlation. On the other hand, the peak power at 0.12 Hz in the RSNA spectrum was in the frequency range from 0.04 to 0.15 Hz, which has been defined as the low-frequency (LF) component in the power spectra of heart rate variability in humans and dogs (2, 27). The high-frequency component, from 0.15 to 0.4 Hz, has been considered to reflect the parasympathetic activity modulated by respiration. Recently, there have been some debates about whether the LF component reflects the fluctuation of sympathetic activity (28) or not (1, 6,10). However, the frequency of the peak power was 0.12 Hz in both the HR and the RSNA spectra andI_max(HR, RSNA) was high (0.273 ± 0.022). According to the analysis inSimulation 1, these findings suggest that the LF component in the HR spectrum is modulated at least by the sympathetic activity.

Gebber and colleagues (16, 33) reported that a rhythmic fluctuation in the LF range persisted after denervation of all cardiovascular nerves. This suggests that the rhythmic discharge of the peripheral sympathetic nerve might result from a sympathetic rhythm of central origin that is normally entrained to the cardiac cycle by the pressor receptor input (35). To examine this hypothesis, we thought that the artificial heart model might be useful. If the cardiovascular regulatory system causes the rhythmic fluctuations of hemodynamic parameters through the natural heart, sympathetic nerve activity should not cause rhythmic fluctuations with the total artificial heart. The rhythmic fluctuations (0.12 Hz) of RSNA and BP were observed in a pacing rate of 60 beats/min (1 Hz). It is considered that these fluctuations are spontaneous because the rate of the fluctuations was different from the pacing rate. With this finding, the high values ofI_max(RSNA, BP) suggest that these cyclic fluctuations are correlated strongly and are not induced by driving the artificial heart.

Using coherency analysis, Brown et al. (10) have reported that sympathetic interactions are more reliably reflected in the BP power spectrum than in the HR spectrum in the conscious rat (10). Indeed, Fig. 8 shows that the coherence between RSNA and BP was stronger than that between RSNA and HR in the LF band. We consider that these coherencies might reflect only the linear correlations between the series. In the present study the MI values among HR, RSNA, and BP during control indicate that the correlation between HR and RSNA is stronger than that between HR and BP and that between RSNA and BP. This might be because the sinus node is affected by efferent sympathetic nerve directly, although it responds to both sympathetic and vagal inputs.

We calculated the time delay of each pair of signals derived from the HR, RSNA, and BP signals by regarding a time delay at which the maximum MI value between the signals was given as a physiological delay. Such a time delay T does not indicate a time delay between two time series at a certain frequency of fluctuations originating from respective physiological inputs, for instance, sympathetic or vagal input, but instead means that one of these time series is most strongly correlated with the time series that is delayed for the time delay T from the other time series. The results on detection of the time delay between two torus series showed the superiority of MI analysis to coherency analysis. However, we must be very careful in interpreting the delay relations between these parameters because, during closed-loop operation of the control state, the relationship between HR and BP reflects both the feedback from BP to HR through the baroreflex and the moderating effect of HR on this response by changing BP through changes in cardiac output. Moreover, efferent sympathetic and vagal activities directed to the sinus node are characterized by discharge largely synchronous with each cardiac cycle that can be modulated by central (vasomotor and respiratory centers) and peripheral (oscillation in arterial pressure and respiratory movements) oscillators (23). These oscillators generate rhythmic fluctuation in efferent neural discharge and construct regulatory feedback loops. Although our results are difficult to account for precisely, they suggest that RSNA leads BP, BP leads HR, and RSNA leads HR during the control state in anesthetized dogs. Our results are compatible with the findings reported by Pagani et al. (27); however, there have been only a few reports available regarding the physiological delay between these parameters (31).

Perhaps the most difficult finding to explain is the small lags between RSNA and BP during the artificial heart state with pacing rates >60 beats/min. The demonstration in other studies (10, 12) that sympathetic activity can follow respiration has suggested that sympathetic nerve activity (SNA) can respond to rapid disturbances. Persson et al. (30) have reported that the phase lag of SNA and BP was roughly 0.2 s. Similarly, Fig. 2 shows that SNA can respond to the rapid increase of BP. On the other hand, there were considerable variabilities of RSNA and BP at a broader frequency band with a pacing rate of 160 beats/min than with a pacing rate of 60 beats/min (Figs. 9 and 10). With a pacing rate of 160 beats/min, the coherencies between RSNA and BP at some frequencies >0.25 Hz were significant, and the phase angles at the respective frequencies were roughly zero radian. These results were compatible with those from MI analysis in that they were strongly correlated and the time delays between them were <0.25 s. Although we should take the sampling frequency (8 Hz) into consideration to interpret “roughly zero” lags, these observations raise the possibility of a rapid response of RSNA to BP with an increasing pacing rate.

During the artificial heart state,I_max(RSNA, BP) was significantly larger with a pacing rate of 60 beats/min than with other pacing rates. This finding indicates that the fluctuation of sympathetic nerve is strongly correlated with that of BP with a pacing rate of 60 beats/min.I_max(RSNA, BP) decreased with increasing pacing rate. This might suggest that a response of the baroreflex system to BP depends on pacing rates of the artificial heart and that the response can hardly keep up with the increasing pacing rate. Thus we believe that the neurological effect of the total artificial heart must be considered when its driving conditions are selected.

Limitations

We should take effects of anesthesia on the cardiovascular control system into consideration. Most of the gaseous anesthetic agents, including the nitrous oxide used in the present study, depress sympathetic activities in different regions of the sympathetic nervous system (13, 20). However, Keyl et al. (18) demonstrated, by means of spectral analysis of heart rate variability, that general anesthesia had no disadvantageous effects compared with local anesthesia on the perioperative cardiac autonomic tone during ophthalmic surgical procedures in otherwise healthy patients. Even if the effects of anesthesia on cardiac autonomic tone must be considered, we think that the correlation among HR, RSNA, and BP could be assessed in this study because the frequencies of peak power were consistent in the spectra of these parameters (Figs. 8 and 9). This is supported by Brown et al. (10), who found that anesthesia with pentobarbital lowered spectral power for HR, SNA, and BP but essentially did not change the coherence between SNA and BP.

In conclusion, MI analysis has many features that are advantageous in assessing the relationship between various cardiovascular variables.1) The method consists of probabilistic processes that do not require linearity of data and thus is applicable to the assessment of data generated not only by a linear system but also by a nonlinear one.2) MI analysis can sensitively quantify the relationship between two time series and assess a time delay of the interaction between them.3) The MI algorithm is simple and quick to compute. We applied MI analysis to quantify the correlation among HR, RSNA, and BP in dogs during the natural heart state and during the total artificial heart state. We conclude, from both theoretical and practical perspectives, that MI analysis has a potentially wide range of applications in studies of cardiovascular variation.