Multivariate autoregressive modeling of autonomic cardiovascular control in neonatal lamb

COMPUTERS

AND

BIOMEDICAL

RESEARCH

21, 512-530 (1988)

Multivariate Autoregressive Modeling of Autonomic Cardiovascular Control in Neonatal Lamb SEPPOKALLI,*

JUHANI ~RONLuND,t ILKKA VALIMAKi,t

HEIMO IHALAINEN,~ ANJA SIIMES,? AND KARI ANTILA~

*Medical Engineering Laboratory, Technical Reseurch Centre c;f’Finland; *Cardiorespirutor Research Unit, Uniuersity of Turku; and SMeasurement Engineering Laboratory, Tampere University of Technology. Tampere, Finland Received June 22, 1987

The neonatal cardiovascular control system is a complicated interactive system which is under vigorous development at birth. From the measurement point of view the cardiovascular control is a closed-loop system. However, it can be examined on a beat-by-beat basis by analyzing circulatory-controlled variables with advanced signal analysis techniques. This paper proposes to use a multivariate autoregressive modeling technique in the analysis of several simultaneous physiological signals in order to examine interactions and inherent properties in the system. With the proposed multivariate autoregressive modeling technique, a signal is modeled as a linear combination of its own past and the past values of the other simultaneous signals plus a predictive error term of the model. The interactions in the system after the model identification are analyzed in frequency domain utilizing power spectrum estimates of the signals and signal contributions. The applicability of the proposed method was examined by a three-variable model between heart rate. blood pressure and respiration in the study of autonomic cardiovascular control in a chronic neonatal lamb model, in which the cardiovascular status was changed by using a P-adrenergic autonomic nervous blockade. The study showed that the multivariate autoregressive modeling technique is a feasible technique in studying complicated interactions within the cardiovascular control system. 0 1988 Academic Press. Inc.

1. INTRODUCTION The cardiovascular control system is a complicated system that comprises numerous parallel subsystems. The purpose of the control system is to guarantee adequate tissue perfusion in each part of the body in all conditions of life. In principle, the cardiovascular control comprises a negative feedback loop with a time delay. Therefore, spontaneous oscillations and fluctuations occur in blood pressure and heart rate. These are mediated by neural and humoral control mechanisms. From the systems analysis point of view the cardiovascular control system is a closed black-box system. However, the analysis of the oscillations and interrelationships between some cardiovascular variables, such as blood pressure 512 OOlO-4809/88 $3.00 Copyright B 1988 by Academic Press. Inc. All rights of reproduction in any form reserved.

AR MODELING

OF CARDIOVASCULAR

CONTROL

513

and heart rate, may be used to provide some criteria to interpret the interplay between the complex control mechanisms. Recent advances in computer systems and signal analysis techniques have greatly enhanced the study of the cardiovascular control (3, 8, 34). In recent years growing interest has been aroused to quantify the fast beat-to-beat variability in heart rate and blood pressure (10, 15, 17,25). Variability in heart rate has been studied with a heart rate variability signal corresponding to the instantaneous heart rate or heart interval. One major advantage of heart rate variability signal is that it can be easily derived noninvasively from a single lead ECG. Heart rate variability has been usually assessed by simple indices (24, 39). Conventional power spectrum analysis methods have been applied to study various components in heart rate variability (2, 4, 32). Data modeling has obtained an important role in signal processing and spectral analysis (23). In parametric modeling, a signal is represented by a set of parameters which are then used for subsequent analysis and description (9,28). Parametric autoregressive identification techniques have been applied to the analysis of heart rate and blood pressure signals (7, 26). Multivariate autoregressive modeling has proved to be a promising tool for quantitative analysis and description of cardiovascular control system in human adults (20,21). The multivariate autoregressive model describes a system where all signals are related to each other. The dynamics of the system can then be examined on the basis of the model. Several scientists have shown that there are specific regions of interest in the frequency spectrum of heart rate and blood pressure (2, 16, 17, 36). One is around 0.25 Hz, at the frequency of respiration, another caused by the baroregulatory oscillation is at about 0.1 Hz. A third one in the lower-frequency region from 0.04 to 0.0s Hz results from the thermoregulatory activity and from local adjustment of resistance in individual vascular beds. The closed-loop feedback system coupling blood pressure, heart rate, and respiration is depicted in Fig. 1. The transfer function Hlz describes the heart rate’s influence on blood pressure and &i, vice versa. The transfer function H3i describes respiration’s influence on blood pressure and Hj2, influence on heart rate. All the variables have inherent sources (ni, n2, n3). Short-term control of the cardiovascular system can be described as a closedloop system, where blood pressure is continuously monitored by the baroreceptors. These send afferent information to the central nervous system, which processes it and accordingly adjusts the length of the heart interval, i.e., the heart rate, as well as the degree of peripheral vascular resistance. Heart rate has respectively an influence on blood pressure through cardiac output and peripheral vascular resistance. Respiration acts primarily as an external disturbance inducing heart rate and blood pressure changes. The classical assumption is that the influence of respiration on blood pressure is mainly mediated by direct mechanical effects as the Bainbridge reflex (6). Thus breathing would act directly on the intrathoracic vessels and stroke volume and at the same time on the baroreceptors. Another assumption is that variability arises from the cen-

514

KALLI

ET AL.

RESPIRATION X3

X2 HEART RATE

FIG. 1. Three-variable respiration

X, BLOOD PRESSURE

closed-loop feedback system coupling blood pressure, heart rate, and

tral nervous system affecting both the heart rate and the respiratory rate at the same time (27). The autonomic cardiovascular control system is not fully developed in newborns. It has been found in neonatal lambs (35) and newborn infants (15) that some maturation occurs during the first weeks of life. Also, the baroreceptor sensitivity has been found to be lower in a newborn lamb than in an adult sheep (14). The same phenomenon has been found also in humans. Monitoring of the heart rate variability has been found to be useful for the evaluation of surveillance of the newborn at risk of asphyxia or in respiratory distress (5, 24). Only a few studies have concerned the role of the autonomic nervous system in neonates (36, 39, 40). This paper introduces a new parametric multivariate autoregressive (MAR) method for the analysis of the cardiovascular control in neonates. The proposed MAR method is applied to the analysis of interactions between blood pressure, heart rate, and respiration in experimental studies performed in chronically instrumented neonatal lambs using /?-adrenergic receptor blockade. 2. ANALYSISOFA 2.1. Multivariate

CLOSED-LOOPSYSTEM WITH MULTIVARIATE AUTOREGRESSIVEMODELING

Autoregressive Model

The multivariate autoregressive (AR) model describes each signal as a linear combination of its own past values and the past values of the other signals plus a prediction (modeling) error term. The model describes a system where all the signals explain each other. The multivariate autoregressive (AR) stochastic process x(k) = [x](k), x2(k), . . . , I,]* of m variables is described by the equation

AR MODELING

OF CARDIOVASCULAR

515

CONTROL

HI

x(k) = 2 a(i)x(k - i) + e(k), i=I where

are the AR coefficient matrices, and e(k) is a white noise vector process e = [Q(& 4V, . . . , e,(~?)]~. According to [l] the process x(k) depends only on the past values of the variables. The model order M denotes the number of the past values that are used to generate a new value for the process. The signal x(k) can be considered the output of a multivariate filter whose input is the white process e(k). The correlation function of this driving white source e(k) is (33) E{e(i + k)eWT} = 6(k)&

[31

where I Cl1 z

=

CT]&7 . . .

Ulm \

j: : . “”

\ Ui?ll

722

urn2

: . .

* . .

and

S(k) =

1, when k = 0; 0,

when k f 0.

141

~mm I

The diagonal elements (+ii of the variance-covariance matrix 2 are the variances of the white sources ei and the nondiagonal elements oij, i # j are the covariances between sources ei and ej. When the sampling interval is Af = Il.&, where fS is the sampling frequency, the spectral matrix of the white process e(k) is S,(f)

[51

= AtI;.

Using the delay operator z-r and the Z transform, the convolution sum in Eq. [I] becomes a multiplication of the matrix of the polynoms of z-r X(z) = A(z)X(z)

[61

+ E(z),

where /A,,(z)

A12tz)

. . .

A,,(z)\

. &I(Z)

Amz(z)

. . .

&m(z)

i

[71

516

KALLI

ET AL

The Z representations of the vector processes x(k) and e(k) are X(z) = CT=-,z-‘x(i) and E(z) = Cy=“=,z-‘e(i), respectively. Now X(z) = C(z)Ek),

L81

C(z) = 11 - A(z)]-',

]91

where C(z) is the transfer function from E(z) to X(z). The element C,-(z) of the transfer function matrix C(z) describes the effect of the source signal ej on the signal xi. The spectrum of the autoregressive process can now be derived from Eqs. [S] and [g]: S,(f)

= C(ei2~rlr,)S,(f)C(ei2~fifS* = At[I - A(eiz~.f~r’/,)]-‘~[I

- A(ef2r.f/f,)]--*.

LlOl

In the above equation * denotes the complex conjugate transpose of a matrix and -*, the inverse of the conjugate transpose matrix. The spectrum can be described in the units of frequency fby substituting z = ei2rf’f\ on the unit circle of the z plane. Equation [I] describes a one-step prediction model where the value of the vector signal x(k) at a certain time instant can be estimated by using M previous values of the same signal. The estimation error is the white signal e(k). If the signal x(k) is a true autoregressive process, then Eq. [l] can be used as the one step ahead predictor. The problem of fitting an autoregressive model to the data has been studied already in (12, 22, 38), and several identification algorithms are available today. One of the most used methods is the Levinson algorithm, which is based on multivariate Yule Walker equations utilizing covariances of the autoregressive processes. The model order selection can be based on several algorithms optimizing the size of the prediction error and the model order (I, 33). 2.2. Analysis

of Transfer

Paths in a Closed-Loop

System

A multivariate closed-loop system, as the one in Fig. 1, can be represented with a transfer function model r

x(k) = c h(i)x(k - i) + n(k),

1111

i=O

where 0 h(i) =

h(j)

W

0

,h,i(i)

h,;(i)

h1

.

.

.

hdi)

.

.

.

hdj)

.

.

. . .

0

9 I

for i = 0, 1, . . . , 2,

[12]

AR MODELING

OF CARDIOVASCULAR

517

CONTROL

and n(4

=

[n~(fV,

n2W,

. . . , ~&)lT.

[I31

The direct transfer function from signal Xj to xi is described by the impulse response function hij(k). Inherent components ni of each signal xi can be considered the input signal sources of the system. Usually the signals ni are not white. If the white noise processes ei(k) are mutually independent then the 2 matrix (Eqs. [4] and [5]) is diagonal. The transfer function model with the inherent signal sources can be expressed in terms of an autoregressive model (29-31). The inherent signals ni form a vector signal n(k), which is expressed n(k) = 2 f(i)n(k - i) + e(k), i=l

where

i: : . *I .fllG)

f(i) =

..a

0

*

: * *

0 .

0

. . . . . fmm(i)

0

0

.

f-22(4

0

The inherent component ni of inherent components nj, j # i. univariate autoregressive model Defining aD = diag(a) and aN =

,

fori=

. . ,M.

[I51

1,2,.

. . ,M

[I61

. . ,M;

1171

for i = 0;

0, aN(i) + x:If

aD(i)h(i -j),

i x.j”=1 aD(i)h(i - j), With the Z representation

1,2,.

each signal xi is independent from the other Each inherent component ni is modeled by a with coefficients &(l), 1 = 1, 2, . . . , M. a - aD, the matrices h(i) and f(i) become

f(i) = aD(i),

h(i) =

fori=

fori=

1,2,.

for i L M.

Eqs. [16] and [I71 become F(z) = AD(z),

[181

and H(z) = [I - AD(z)]-‘AN(z),

1191

where AD is the diagonal and AN the nondiagonal part of the matrix A. The element Hij(z) of the transfer function matrix H(z) is the transfer function from the signal Xj to the signal Xi. It characterizes the transfer path from the signal xj to the signal xi. The total effect of the inherent source nj of the signal xj to the signal xi is called the power contribution of j to i. Assuming that the inherent sources are mutually independent, the power contribution of the source n.i to the variable xi

518

KALLI

ET AL

is equal to the power contribution of the white noise source ej to the variable xi. Thus the power contributions of the inherent sources can be calculated using the transfer functions C,(z) from Eq. 191. The power contribution of the source n,i to the variable xi is @j(f)

= /C;j(ei?~flf~)I’S,jj(f)

[201

= AtiC,;( eM,)12q,.

[211

The spectrum of the signal xi can now be represented as a sum of the spectral effects from all the sources in the system sii(fJ

= 2 Qij(,f). ,=I

1221

The power contribution Qii(f) of the inherent source to the signal itself is the spectrum of the source. It is called the autocontribution of the signal xi. The contribution of the other sources are called crosscontributions. Scaling Q,(f) by dividing it with the autospectrum of the signal -Xi, the power contribution ratio is obtained

The sum of the contribution signal Xi is equal to 1.

ratios of all the inherent

3. EXPERIMENTAL 3.1. Materials

system sources to the

PROCEDURE

and Signal Acquisition

To illustrate the applicability and the potentials of the proposed MAR method, it was applied to the study of the cardiovascular control in an 8-dayold healthy Finnish breed lamb. Preparation of the lamb had been done on the third postnatal day (for details see (35)). ECG and transthoracic electrical impedance respirogram were then recorded from implanted silver electrodes at the midaxillary line of the thoracic wall. Arterial blood pressure was recorded continuously using a polyvinyl catheter inserted into the carotid artery. In order to study the autonomic nervous control, p-adrenergic blockade was induced with an intravenous propranolol dose of 1.O mg/kg. Two-minute segments of the recording were chosen in order to obtain sufficiently stationary signals. During the control period before the blockade eleven segments were recorded and during the blockade after the drug administration, five successive segments were examined. 3.2. Signal Analysis

The blood pressure and respiratory signals were digitized at a sampling frequency of 16.7 Hz using a minicomputer system. From the analog ECG signal

AR MODELING

OF CARDIOVASCULAR

CONTROL

519

R waves were detected with a hardware preprocessor and transformed to heart interval time series in a minicomputer (4). This heart interval time series was then transformed to a time series of instantaneous heart rate, the sampling rate of which was further equalized with a linear interpolation to 16.7 Hz. Heart rate time series was synchronized with the other time series. All the signals, blood pressure, heart rate, and respiration were then low-pass filtered and samples with 3 Hz sampling frequency (13, 18). All the subsequent signal analysis procedures were performed on these time series. The above-described MAR modeling technique was utilized for the analysis of the heart rate, respiration, and blood pressure time series. The length of the MAR model was selected to be 20 in all the analyzed cases. The comparison between the control and the blockade recordings and between individual signal segments was then performed by the analysis of the power spectrum density (PSD) estimate and the signal source contribution ratios. For power spectrum density and contribution analysis the investigation was performed separately on five spectral regions: 0.0-0.02 Hz (Fl), 0.02-0.08 Hz (F2), 0.08-0.25 Hz (F3), 0.25-0.40 Hz (F4), and 0.40-0.75 Hz (FS). The interrelationships between the signals were analyzed by signal source contribution analysis in these frequency bands. For the comparison between the groups, the PSD estimates were band-integrated over the above-mentioned frequency bands. Comparisons between recorded signal segments were done with a two-tailed t test. 4. EXPERIMENTALRESULTS The mean heart rate (&SD) in the control segments (n = 11) was 226 + 4.6 bpm and in the P-blockade segments (n = 5) 182 + 2.1 bpm, respectively. The decrease in the heart rate level during the blockade was statistically significant (P < 0.001). The p-adrenergic blockade caused also a marked decrease in the beat-to-beat heart rate variability (Fig. 2). HR(bpm) 260. ’

FIG. 2. Sample plots of heart rate recordings during the control condition and during the /3adrenergic blockade.

520

KALLI

ET AL.

Sample plots of respiration and blood pressure from the two conditions are depicted in Fig. 3. The rhythm and the amplitude of respiration during the blockade was found to be more regular than in the control condition. Visual examination revealed no clear differences in the blood pressure recordings between the blockade and control conditions. 4.1. Power Spectrum

Estimate

The analysis of power spectrum density estimates revealed a marked decrease in the heart rate variability after the blockade at all frequencies (Fig. 4a) The slightly more regular respiratory rhythm during the blockade could be seen as a more prominent spectral peak at the respiratory frequency (Fig. 4b). The power spectrum estimates of blood pressure showed a power decrease in all but the respiratory frequency region after the blockade (Fig. 4~). Table 1 represents the average power spectrum density estimates derived with integration over the five different frequency bands (FI-FS). The decrease in the heart rate power after the blockade in all but the lowest frequency band was statistically significant. There was a 35 to llO-fold decrease in the heart rate variability at the frequency bands F2-FS (0.02-0.75 Hz). The decrease was

TABLE AVERAGE

POWER

SPECTRUM

DENSITY

ESTIMATES

FIVE

FREQUENCY

F1

F2

I IN THE Two

ANALYZED

CONDITIONS

F4

F3 _----.--.-

HR (bpm) Control Blockade

RESP (arbitrary Control Blockade

BP (arbitrary Control Blockade

Note. 0.40-0.75

AT

REGIONS

Mean SD Mean SD P

36 49 1.2 0.9 0.135

5Y 57 I.7 2.0 0.042

56 35 0.9 1.0 0.004

8.8 5.9 0.1 0. I 0.006

units) Mean SD Mean SD P

558 648 229 362 0.313

1615 1676 508 540 0.175

2972 2050 1394 1516 0.146

3278 2175 2773 2086 0.673

units) Mean SD Mean SD P

725 629 430 476 0.628

1575 1780 430 554 0.186

397 251 232 245 0.238

76 84 20 9.1 0.167

FS .- __-- ~~-~~~~~ 22 5.1 0.1 0.1 0.001 51.050 33.671 79,270 12,048 0.092 192 60 104 20 0.007

Fl denotes 0.0-0.02 Hz; F2, 0.02-0.08 Hz; F3, 0.08-0.25 Hz; F4, 0.25-0.40 Hz; and F5, Hz, respectively. P denotes the statistical significance derived with the two-tailed t test.

arb.unit. 500. '

-500.

! 0.0

20.

40.

60.

I 100.

80.

Time/k arbunit 500.

b

300. 100. -100. -300.

-500. !’ 0.0

20.

40.

60.

80.

100. Time/k

arb.unit 0.0 -25. -50. -75. -100. -125. -150. 0.0

20.

40.

60.

80.

100. Time/Set

arb.unit O.Oi

,

d 1

-25. -50. -75. -100. -125.

Time/%

FIG. 3. Sample recordings of respiratory and blood pressure signals. (a) Respiration during the control condition. (b) Respiration during the blockade. (c) Blood pressure during the control condition. (d) Blood pressure during the blockade. 521

522

KALLI

ET AL

Frequency

(Hz)

&wit

m

0.0

0.25

0.38

0.50

i

0.63

0.75

Frequency

(Hz)

ob. tit XlOK

i cl

IO K

Frequency

(Hz)

FIG. 4. Power spectrum estimates during the control and blockade condition. (a) Spectrum estimates for heart rate. (b) Spectrum estimates for respiration. (c) Spectrum estimates for blood pressure.

greatest at the highest frequency bands. The respiratory power decreased in all but the highest frequency band. The differences between the control and the blockade conditions were not statistically significant considering the respiratory power spectrum density and its distribution over the analyzed five fre-

AR MODELING

OF CARDIOVASCULAR

CONTROL

523

80. 60.

0.0

0.0 0.10

0.40

0.20 0.30

0.60 0.50

O.‘O

4 .8( “..

FIG. 5. Relative (O-100%) signal source contribution to blood pressure during a control condition. The lower open area corresponds to the contribution caused by blood pressure, the hatched area corresponds to respiration, and the upper open area corresponds to heart rate.

quency regions. The power in the blood pressure variability decreased approximately 2-fold in all frequency regions as a response to the p-blockade. The decrease at the highest frequency band, F5, was significant (P < 0.01). 4.2. Contribution

Ratios

A sample plot of contribution ratios to blood pressure in the control condition is depicted in Fig. 5. Generally it was found that the relative crosscontribution decreased and the relative autocontribution increased after the /3-adrenergic blockade. Table 2 represents the average signal source contribution ratios and Fig. 6 depicts the results in a stacked bar graph form. The respiratory source contribution increased to both heart rate and blood pressure after the blockade. The contribution of respiration to heart rate increased from 22 to 40% of total variability in the frequency region F5 after the blockade (P < 0.05). Also the heart rate autocontribution in the frequency band F3 increased from 68 to 78% of the variability (P < 0.05). Also in the F2 region, the heart rate autocontribution increased profoundly from 51 to 78%. The respiratory autocontribution increased during the blockade in all but the lowest frequency region. The increase was marked in the F5 region, from 66 to 77% of the total variability. The blood pressure autocontribution decreased in the F5 region from 45 to 30% of the total variability in the blockade group. The crosscontribution of respiration to blood pressure, on the other hand, increased from 44 to 62% of the total variability (P < 0.05). Also the crosscontribution of heart rate to blood pressure decreased in all but the Fl region after the blockade. The decrease was marked in the F3 region from 25 to 14% of the total variability. 5.

DISCUSSION

Simple stochastic point process analysis is of limited value in featuring of the autonomic cardiac control although it has been successfully applied in clinical

;5 P

P

P

to BP (%) Mean SD Mean SD

to RESP Mean SD Mean SD P

to HR (%) Mean SD Mean SD

SIGNAL

(%I

Norm>. Fl denotes 0.0-0.02 two-tailed I test (* denotes

Blockade

Contribution Control

Blockade

Contribution Control

Blockade

Contribution Control

AVERAGE

BP

14 14 20 20 0.56

62 25 59 25 0.87

I3 II 15 II 0.67

RESP __.__

Fl

CONTRIBUTION

- Hz: F?, 0.02-0.08 P < 0.05).

56 27 48 29 0.58

21 14 24 20 0.78

49 31 28 24 0.22

SOURCE

62 18 60 II 0.83

24 12 22 12 0.80

HR

26 I8 I6 10 0.2s

I? 8.4 24 13 0.04*

Hz: F4. 0.2iXl.40

__--

15 8.1 12 9.5 0.57

51 24 78 19 0.05”

__

CONTROL

61 13 66 16 0.58

12 7.6 9 9.4 0.50

RESP

F2

FOR THE

37 25 I3 I? 0.07

BP

Hz: Fi. 0,08-O.?

30 21 32 I9 0.92

17 16 17 16 0.97

38 35 57 35 0.65

HR

RATIOS

BP

0.4lI-0.77

57 I1 64 II 0.25

18 12 13 12 0.55

F3

25 8. I I4 II 0.05*

I3 7.1 I3 II 0.92

68 8.1 78 8.2 0.03*

HR

SEGMENTS

H/. w\pwrivcl!

I8 I5 22 12 0.71

69 IO 74 21 0.56

13 12 6 3.5 0.20

RESP

BLOCKADE

2

I9 8.7 16 8.1 0.58

-.

I-17. 1”

AND

TABLE THE

I’dcru~tc\

62 19 78 10 0.09

23 14 13 14 0.26

13 8.8 11 8.5 0.78

BP

OVER

I9 I3 10 7.7 0.19

9 6.6 10 7.3 0.88

74 9.7 77 4.9 0.63

HR

4s I? 30 I9 0.08

24 9.5 I5 II 0.1 I

21 9.7 12 7.3 0.06

BP

BANDS

\Ignitii:incc

FREQUENCY

the* \I:itl\llc,1/

I9 IO I2 7.1 0.16

68 I5 77 21 0.64

13 7.6 12 7.1 0.81

RESP

F4

FIVE

dcrivsd

44 I? 62 18 0.03”

66 IO 77 11 0.05*

22 9.5 40 20 0.03*

RESP

FS

-

II 5.b 8 2.6 0.26

IO 5.2 8 3.0 0.31

57 II 48 I7 0.29

HR

with the

(Fl-F.5)

AR MODELING

a

SIQNAL

OF CARDIOVASCULAR SOURCE

525

CONTROL

CONTRIBUTION SIGNAL 0

SOURCE HI? RESP

m

Fl FZF3F4F5

BP

SIaNAL

b

Fl FZF3F4F5 CONTRIBUTION

Fl

F2F3F4F5

TO

RESP

SOURCE

HR

CONTRIBUTION SIGNAL

90

0

m

g

SOURCE HI? RESP

60 8\ -

BP

BP

70 60 50

3 E

40

8

30 20 10 0 Fl

F2F3F4F5

BP

F1 F2F3F4F5 CONTRIBUTION

RESP

Fl F2F3F4F5 TO

HR

FIG. 6. Average signal source contribution ratios to blood pressure (BP), respiration (RESP), and heart rate (HR) over five spectral regions (Fl-FS). Heart rate signal source is denoted with open bars, respiration with hatched bars, and blood pressure with solid bars, respectively. (a) Contribution during the control condition. (b) Contribution during the blockade.

monitoring (5, 37). Spectral analysis of heart rate variability has made it possible to quantify the impact of external disturbance factors (respiration, vasomotor thermoregulation) or blood pressure on heart rate (25, 26). This study deals with autoregressive modeling of autonomic /3-adrenergic cardiovascular control on heart rate, respiration, and blood pressure in a neonatal lamb. The recordings were analyzed with a three-variable autoregressive model. The method provides power spectrum estimates and source contribu-

526

KALLI

ET AL.

tion ratios, which form the basis for the characterization of the system under study. In general, the use of the proposed autoregressive modeling technique in the analysis of a closed-loop control system is done in the following step-like fashion: (1) Select the data. (2) Fit the autoregressive model to the data (the Levinson algorithm). (3) Compute the power spectrum estimates. (4) Compute the signal power contribution ratios. (5) Characterize the system with the computed results, i.e., the power spectrum estimates describe the variability and the signal power contribution ratios describe the interactions. The experimental results showed pronounced differences between the control condition and the P-adrenergic blockade. The most profound difference was found in heart rate. Its mean level and variability obtained by power spectrum density analysis decreased significantly after the blockade. The power of respiration at the respiratory frequency itself, the F5 region, seemed to increase after the blockade. Probably the respiratory rate was sligthly more regular after the blockade. The contribution analysis revealed that there was a large (over 50%) inherent variation in all the variables, heart rate, respiration, and blood pressure, in all the investigated frequency regions. The autocontribution increased during the blockade, especially in heart rate and blood pressure. Accordingly, the crosscontribution of heart rate to blood pressure and the one from blood pressure to heart rate was smaller after the blockade. This was significant in the F2 and F3 regions, i.e., in the vasomotor thermoregulatory and baroregulatory regions, respectively. The results seem to be related to the fact that the used p-adrenergic receptor blocking agent, propranolol, was of nonspecific type affecting the P-receptors both in the heart and in the peripheral vasomotor musculature. Featuring of a complicated biological system is more reliable if a large number of measured variables can be used for the study. On the other hand, the properties of the system are individual and subject-to-subject variability is considerable. Therefore, we decided to examine the same animal during the blockade and under control conditions. In general, cardiovascular signals, such as blood pressure and heart rate, have high mutual correlation (8, II) because arterial blood pressure is instantaneously regulated by heart rate. stroke volume, and peripheral resistance. In practice, the signal sources of the closedloop model of Fig. 1 correlate always to some extent. Therefore the covariance matrix of the corresponding MAR model is not diagonal. The proposed contribution analysis, however, presumes that the sources are uncorrelated and takes account only of the diagonal part of the covariance matrix. This causes a certain bias to the derived contributions. Often it is not possible to record respiration as we did in our study. A typical situation is the ambulatory blood pressure monitoring (19), where we have to

AR MODELING a

SIQNAL

OF CARDIOVASCULAR SOURCE

CONTROL

527

CONTRIBUTION

100 90 60 i

M -

70..

g

50..

60..

*O t 10 0 f Fl

F2

F3

F4 F5 CONTRIBUTION

Fl

F2 TO

SIQNAL

F4

F5

HR

BP b

F3

CONTRIBUTION

SOURCE

SIONAL SOURCE 90

E::

60

K % 5

70

E

40

0

30

60 50

20 10 0 Fl

F2

F3

BP

F4 F5 CONTRIBUTION

Fl

F2 TO

F3

F4

F5

HR

FIG. 7. Average signal source contribution derived with a two-variable autoregressive model between heart rate (HR) and blood pressure (BP). Contribution of heart rate signal source is denoted with open bars and that of blood pressure with solid bars, respectively. (a) Contribution during the control condition. (b) Contribution during the blockade.

analyze the properties of the cardiovascular control system on the basis of recordings of blood pressure and heart rate only. Then respiration is a clear external noise source that affects both of the recorded signals. The signals have therefore a high intercorrelation at respiratory rate. To investigate the behavior of the MAR modeling method with correlated signal sources, we also applied it to the analysis of two signal variables, heart rate and blood pressure. On the basis of the obtained average signal source contribution ratios (Fig. 7) it could be found that the autocontribution increases and the crosscontribution decreases during the blockade as in the three-vari-

528

KALLIETAL

able MAR model above. The respiratory impact of the crosscontribution both to heart rate and to blood pressure in the three-variable model seems now to be included in the inherent signal sources in the two-variable model. The results show that the MAR method provides reasonable results also with correlated sources. This confirms the applicability of the proposed MAR model also in the analysis of interactions between heart rate and blood pressure as done in (20, 21). 6.

CONCLUSIONS

This paper introduces a new parametric analysis method to the study of autonomic cardiovascular control. The method utilizes a multivariate autoregressive model to characterize a closed-loop system coupling the measured variables. The method was applied to the study of cardiovascular control in a neonatal lamb. The method proved to be a feasible way for studying complicated interactions between heart rate, blood pressure, and respiration. REFERENCES 1. AKAIKE. H. A Bayesian extension of the minimum AIC procedure of autoregressive model fitting. Biometrika 66, 237 (1979). 2. AKSELROD, S.. CORDON, D., UBEL, F. A., SHANNON, D. C., BARGER, A. C.. AND COHEN. K. J. Power spectrum analysis of heart rate fluctuation: A quantitative probe of beat-to-beat cardiovascular control. Science 213, 220 (1981). 3. AKSELROD, S., CORDON, D., MADWED, J., SNIDMAN, N.. SHANNON, D. C.. AND COHEN, R. J. Hemodynamic regulation: Investigation by spectral analysis. Amer. J. Physiol. 249, H867 (1985). 4. ANTILA, K. J. Quantitative characterization of heart rate during exercise. Stand. J. C/in. Lab. Inuesf. Suppl. 39, 1 (1979). 5. ANTTILA, H., V.&L~M~~KI. I., GADZINOWSKI, J.. KERO, P., AND ANTILA, K. Power spectrumof heart rate in neonates with respiratory distress. In Progress Reports on Electronics in Medicine and Biology” (K. Copeland, Ed.), p. 319. The Institution of Electronic and Radio Engineers. London, 1986. 6. BAINBRIDGE, F. A. The influence of venous filling upon the rate of the heart. J. Physiol. (London) 50, 65 (1915). 7. BASELLI. G., BOLIS. D., CERUTTI. S.. AND FRESCHI, C. Autoregressive modeling and power spectral estimate of R-R interval times series in arrhythmic patients. Cornput. Biomed. Res. 18, 531 (1985). 8. BASELLI, G.. CERUTTI, S., CIVARDI, S.. LIBERATI. D.. LOMBARDI, F.. MALLIANI, A., AND PAGANI. M. Spectral and cross-spectral analysis of heart rate and arterial blood pressure variability signals. Comput. Biomed. Res. 19, 520 (1986). 9. BENDAT, J. S.. AND PIERSOL, A. G. “Random data: Analysis, measurements, procedures.” Wiley-Interscience, New York, 1971. 10. DE BOER, R. W., KAREMAKER, J. M., AND STRACKEE, J. Beat-to-beat variability of heart interval and blood pressure. Automedica 4, 217 (1983). Il. DE BOER, R. W., KAREMAKER, J. M., AND STRACKEE, J. Relationships between short-term blood pressure fluctuations and heart-rate variability in resting subjects I: A spectral analysis approach. Med. Biol. Eng. Comput. 23, 352 (1985). 12. CHATFIELD, C. “The analysis of time series: An introduction.” Chapman & Hall, London, 1984.

AR MODELING

OF CARDIOVASCULAR

CONTROL

529

13. CROCHIERE, R. E., AND RABINER, L. R. Optimum FIR digital filter implementation for decimation, interpolation, and narrow band filtering. IEEE Trans. Acoust. Speech Signal Process. ASSP-23(5), 444 (1975). 14. DAWES, G. S, JOHNSTON, M., AND WALKER, D. W. Relationship of arterial pressure and heart rate in fetal, new-born and adult sheep. .I. Physiol. 309, 405 (1980). 15. GIDDENS, D. P., AND KITNEY, R. I. Neonatal heart rate variability and its relation to respiration. J. Theor. Biol. 113, 759 (1985). 16. HYNDMAN, B. W. The role of rhythms in homeostasis. Kybernetik 15, 227 (1974). 17. HYNDMAN, B. W., KITNEY, R. I., AND SAYERS, B. McA. Spontaneous rhythms in physiological control systems. Nature (London) 233, 339 (1971). 18. JOKIPII, M., KALLI, S., AND SARANUMMI, N. Software interface to ECG databases for arrhythmia detector evaluation. In “Computers in Cardiology,” p. 387. IEEE Computer Society, Long Beach, CA, 1985. 19. KALLI, S., JKRVILEHTO, P., NARIAINEN, K., SARANUMMI, N., SUORANTA, R., TURIANMAA, V., AND UUSITALO, A. Computerized analysis of long-term invasive blood pressure signals. In: “Proceedings, Fifth International Symposium on Ambulatory Monitoring” (C. Dal Palu and A. C. Pessina, Eds.), p. 93. Cleup Editore, Padua, 1985. 20. KALLI, S., SUORANTA, R., AND JOKIPII, M. Applying a multivariate autoregressive model to describe interactions between blood pressure and heart rate. In “Measurement in Clinical Medicine,” p. 77. The Instititute of Measurement and Control, London. 1986. 21. KALLI, S., TURJANMAA, V., SUORANTA, R., AND UUSITALO, A. Assessment of blood pressure and heart rate variability using multivariate autoregressive modelling. Biomed. Meas. Inf. Control 1, 46 (1987-1988). 22. KAILATH, T. A View of three decades of linear filtering theory. IEEE Trans. Znf. Theory 20, 146 (1974). 23. KAY, S., AND MARPLE, S. Spectrum analysis--A modern perspective. Proc. IEEE 69, 1380 (1981). 24. KERO, P. Heart-rate variation in infants with respiratory distress syndrome. Acta Pediatr. Stand. Suppl. 250, 1 (1974). 25. KITNEY, R. I. A nonlinear model for studying oscillations in the blood pressure control system. J. Biomed. Eng. 1, 89 (1979). 26. KITNEY, R. I., FULTON, T., MCDONALD, A.. AND LINKENS, D. A. Transient interactions between blood pressure, respiration and heart rate in man. J. Biomed. Eng. 7, 217 (1985). 27. LANGHORST, P., SCHULTZ, B., LAMBERTZ, M., SCHULTZ, G., AND CAMERER, H. Dynamic characteristics of the “unspecific brain stem system.” In “Central Interaction between Respiratory and Cardiovascular Control Systems” (H. P. Koepchen, S. M. Milton, and A. Trzebski, Eds.), p. 30. Springer, Berlin, 1980. 28. MAKHOUL, J. Linear prediction: A tutorial review. Proc. IEEE 63, 561 (1978). 29. OGUMA, R. A method to estimate mechanical state inside fuel rod based on noise analysis. J. Nucl. Sci. Tech. 17, 811 (1980). 30. OGUMA, R. Method to evaluate signal transmission paths in dynamic systems and its application to reactor noise analysis. J. Nucl. Sci. Tech. 18, 835 (1981). 3J. OGUMA, R. Extended partial and multiple coherence analyses and their application to reactor noise investigation. J. Nucl. Sci. Tech. 19, 543 (1982). 32. PAGANI, M., LOMBARDI, F., GUZZETTI, S., RIMOLDI, O., FURLAN, R., PIZZINELLI, P., SANDRONE, G., MALFATTO, G.. DELL’ORTO. S., PICCALUGA, E., TURIEL, M., BASELLI, G., CERUTTI, S., AND MALLIANI, A. Power spectral analysis of heart rate and arterial pressure variabilities as a marker of sympatho-vagal interaction in man and conscious dog. Circ. Res. 59, 178 (1986). 33. PRIESTLEY, M. “Spectral Analysis and Time Series,” Vols. 1 and 2. Academic Press, London, 1981. 34. SAYERS, B., McA. Analysis of heart rate variability. Ergononics 16, 17 (1973). 35. SIIMES, A., VKLIMKKI, I., SARAJAS, H., SAKONE, K., AND OJA, R. Heart rate variation in relation to age and sleep state in neonatal Iambs. Acta Physiol. Stand. Suppl. 537, 7 (1984).

530

KALLl

ETAL.

36.

SIIMES, A., VALIMAKI, I., OJA, R.. AND ANTILA. K. Detection of components of autonomic cardiac control by time series analysis of heart rate in lambs: technical report. In “Fetal Physiological Measurements” (P. Rolfe, Ed.). p, 2.59. Butterworths. London, 1986. 37. TARLO. P.. VALIMKKI, I.. AND RAUTAHARJU. P. Quantitative computer analysis of cardiac and respiratory activity in newborn infants. J. Appl. Physiol. 31, 70 (1971). 38. UPADHYAYA. B.. KITAMURA, M., AND KERLJN. T. Multivariate signal analysis algorithms for process monitoring and parameter estimation in nuclear reactors. Ann. NIIL./. i?nergy 7. 243

(1980). 39. ZUGAIB,

M.. FORSYTHE. A. B., NUWAYHID. B.. LIEB, S. M.. ‘I‘ABSH. K.. ERKKOL.A. K.. USHIODA. E., BRINKMAN. C. R., AND ASSALI. N. S. Mechanisms of beat-to-beat variability in heart rate of the neonatal lamb. I. Influence of the autonomic nervous system. Amer. J. Obster. Gynecol. 138, 444 (1980). 40. ZWIENER, U.. WITTE. H.. ROTHER, M.. WAGNER. H.. EISELT. M., AND GLASER. S. Sleep state related dynamics of neonatal heart rate fluctuations and respiration rhythms within several frequency ranges. Auromedica 5, 77 (1984).