## Abstract

We applied cardiovascular system identification (CSI) to characterize closed-loop cardiovascular regulation in patients with diabetic autonomic neuropathy (DAN). The CSI method quantitatively analyzes beat-to-beat fluctuations in noninvasively measured heart rate, arterial blood pressure (ABP), and instantaneous lung volume (ILV) to characterize four physiological coupling mechanisms, two of which are autonomically mediated (the heart rate baroreflex and the coupling of respiration, measured in terms of ILV, to heart rate) and two of which are mechanically mediated (the coupling of ventricular contraction to the generation of the ABP wavelet and the coupling of respiration to ABP). We studied 37 control and 60 diabetic subjects who were classified as having minimal, moderate, or severe DAN on the basis of standard autonomic tests. The autonomically mediated couplings progressively decreased with increasing severity of DAN, whereas the mechanically mediated couplings were essentially unchanged. CSI identified differences between the minimal DAN and control groups, which were indistinguishable based on the standard autonomic tests. CSI may provide a powerful tool for assessing DAN.

- cardiovascular modeling
- diabetes mellitus
- hemodynamics
- mathematical modeling
- nervous system

autonomic neuropathy is a well-recognized complication of diabetes mellitus (9). The clinical manifestations of diabetic autonomic neuropathy (DAN) include postural hypotension, gastrointestinal symptoms, hypoglycemic unawareness, and sweating disturbances (7, 9). These clinical manifestations are slowly progressive, usually irreversible (7), and are associated with considerable mortality (8). It is important to quantify the degree of DAN to obtain a physiological measure of the progression of the neuropathy as a guide for clinical management of the diabetes. Consequently, autonomic tests based on cardiovascular reflexes to various physiological perturbations (e.g., heart rate response to deep breathing, postural change, and the Valsalva maneuver and blood pressure response to sustained hand grip and postural change) are commonly employed. However, these tests do not provide a comprehensive, quantitative assessment of autonomic function, are insensitive to early sympathetic neuropathy (7), and generally require the active cooperation of the patient, which may make the test results difficult to reproduce.

In this study, we applied cardiovascular system identification (CSI) for the quantification of DAN. This method, which is described in detail in Ref. 18, mathematically analyzes the beat-to-beat fluctuations in noninvasively measured heart rate signals [derived from the surface electrocardiogram (ECG)], arterial blood pressure (ABP), and instantaneous lung volume (ILV) to characterize quantitatively the physiological mechanisms responsible for the couplings between these variables. Figure1 illustrates the model upon which the method is based. The model includes four physiological couplings [circulatory mechanics, heart rate tachogram (HR) baroreflex, ILV→HR, and ILV→ABP] whose transfer properties are determined from analysis of the fluctuations in the physiological signals. The transfer properties are represented below in terms of impulse response functions. The impulse response represents the response of the coupling mechanisms to an arbitrarily narrow input pulse of unit area at *time 0*.

Circulatory mechanics depicts the coupling of ventricular contraction, represented by a pulsatile heart rate (PHR) signal (see below), to the generation of pulsatile ABP. Circulatory mechanics reflects the mechanical properties of the heart, great vessels, and peripheral circulation and may also include autonomic reflex control of the vascular mechanical properties. The circulatory mechanics impulse response represents the ABP wavelet that results after each cardiac contraction. HR baroreflex represents the autonomically mediated coupling between ABP and HR, which reflects baroreflex control of heart rate. ILV→HR represents the autonomically mediated coupling between respiration, measured in terms of the time variation of ILV, and heart rate. ILV→HR is responsible for mediating the respiratory sinus arrhythmia. ILV→ABP represents the mechanically mediated coupling between respiration and ABP, reflecting the effect of changes in intrathoracic pressure on the filling of intrathoracic structures and venous return.

The model also includes a fifth predefined coupling, sinoatrial (SA) node, which represents the transformation of HR, which is closely related to the net autonomic input to the SA node, to the PHR, which is a train of impulses corresponding to the times of onset of ventricular contractions. The SA node in this model is represented as an integrate-and-fire device and is not identified from the experimental data, since its dynamics are predefined.

Along with the five coupling mechanisms, the model incorporates two perturbing noise sources, N_{HR} and N_{ABP}. N_{HR} represents the perturbations to HR not attributable to the HR baroreflex and ILV→HR couplings, such as autonomically mediated perturbations driven by cerebral activity. N_{ABP} represents perturbations to ABP not attributable to the circulatory mechanics and ILV→ABP couplings, such as fluctuations in peripheral vascular resistance as tissue beds adjust local vascular resistance to match local blood flow to demand. The power spectra of the noise sources are determined by analysis of the physiological signals.

The CSI technique utilizes a linear, time-invariant systems analysis for examining the couplings between fluctuations in each of the physiological signals about their mean values. Although previous studies have shown that these couplings have nonlinear components, the contributions of the nonlinear components are small when considering small fluctuations in these physiological variables about their mean values (6). In particular, the model in Fig. 1 (except for the SA node, which is predefined) is represented by the following two autoregressive exogenous input (ARX) difference equations
where *t* is discrete time, and W_{HR} and W_{ABP} are residual terms. Using continuous records of the measured ECG, ABP, and ILV signals, the six sets of unknown parameters (*a _{i}
*,

*b*,

_{i}*c*,

_{i}*d*,

_{i}*e*, and

_{i}*f*) in these equations are determined by the least-squares minimization of the residual terms in conjunction with an ARX parameter reduction algorithm (20). The algorithm involves choosing an initial, maximal ARX model order that is believed to contain the true ARX parametrization and efficiently reducing this maximal model down to the true parametrization via Rissanen’s minimum description length criterion. We define the initial, maximal model order with the set of variables {

_{i}*m*,

*n*,

*p*′,

*p*} and {

*q*,

*r*,

*s*′,

*s*} in the above equations. We chose these parameters to be {

*m*,

*n*,

*p*′,

*p*} = {10, 15,−10, 15} and {

*q*,

*r*,

*s*′,

*s*} = {10,15,0,35}. The parametrization chosen by the ARX parameter reduction algorithm, along with the residual terms, fully define the transfer properties of the physiological coupling mechanisms and the power spectra of the perturbing noise sources.

The ARX formulation allows the imposition of causality on the closed-loop relationship between ABP and heart rate. Therefore, unlike nonparametric techniques, identification of the distinct feedforward (circulatory mechanics) and feedback (HR baroreflex) coupling mechanisms may be achieved. For reliable identification to be attained, the sampling rates of the ABP, ILV, PHR, and HR signals may be selected to match the bandwidth of the specific coupling mechanisms involved. Precise mathematical definitions of the signals, the details of the identification procedure, and the computation of N_{HR} and N_{ABP} are given in Ref. 18.

In Ref. 18, we demonstrated the ability of the CSI method to characterize the effects of pharmacological blockade and postural changes. Now, we apply this method to physiological data collected noninvasively from diabetic and control subjects to assess the alterations in cardiovascular control mechanisms associated with DAN.

## METHODS

### Subjects

Sixty diabetic patients [age: 47 ± 17 (SD) yr] referred to the Autonomic Evaluation Unit of the New England Deaconess Hospital and 37 control volunteer subjects [age: 42 ± 14 (SD) yr] participated in this study. Experimental protocols were approved by the Institutional Review Boards, and informed, written consent was obtained from all subjects.

### Grading of Subjects by Standard Autonomic Tests

To divide the diabetic patients into groups with differing degrees of autonomic neuropathy, all 97 subjects underwent a battery of standard autonomic tests. This battery included the following tests: heart rate response to respiration (expiratory − inspiratory [E − I] ratio and maximum − minimum), heart rate response to the Valsalva maneuver (Valsalva ratio), diastolic blood pressure increase to isometric exercise, and systolic blood pressure (SBP) fall in response to standing. Procedures for performing these tests are specified in Ref. 13. Two of the tests were selected as measures of parasympathetic and sympathetic function. The logarithm of the high-frequency power in the heart rate power spectrum measured in the supine position [log supine high-frequency power (LSHFP)] was used as a measure of parasympathetic function (1), and the SBP fall (ΔSBP) in response to passive tilting was used as a measure of sympathetic function (10).

*Supine high-frequency power determination*. ECG signals were recorded from each subject in the supine position as described below. The heart rate power spectrum was computed as the Fourier transform of the product of the autocorrelation function of HR with a Gaussian window (3). Supine high-frequency power (SHFP) was computed as the total power in the frequency band between 0.15 and 0.50 Hz.

*ΔSBP determination*. SBP was measured noninvasively with a Critikon Dinamap Vital Signs Monitor 1846SX/P (Johnson and Johnson, Critikon, Tampa, FL) connected to an IBM PC AT via an RS232 serial interface board. SBP was first measured for each subject in the supine posture before passive tilting to 60° (relative to the supine posture) by an electrically driven tilt table. After passive tilting, SBP was measured each minute for a 5-min period. The difference between the lowest SBP measured during this period and the SBP measured in the supine posture was recorded as the ΔSBP.

*Grading of subjects*. SHFP and ΔSBP were measured in all 97 subjects. The SHFP and ΔSBP measures were then age corrected as follows. First, each measure for the control subjects was regressed with age using the following exponential equations
where the*A _{i}
*and

*B*values are unknown parameters, and W

_{i}_{SHFP}and W

_{ΔSBP}are the age-independent contributions to SHFP and ΔSBP in each subject, respectively. These equations were developed to fit empirically the observed dependence of the measures on age in the control subjects. The equations may be interpreted to indicate that each measure has an age-independent component and an age-dependent component. The age-dependent component was found to vary exponentially with age. The variation of SHFP about its mean value for a given age was found to be proportional to the mean value of SHFP. The variation of ΔSBP about its mean value for a given age was found to be independent of its mean value. In order that the regression equation yield nonnegative values of SHFP and that SHFP approach an asymptotic value with increasing age, the following constraints were imposed:

*A*

_{1},

*A*

_{1}+

*A*

_{2}, and

*A*

_{3}≥ 0. These unknown parameters were obtained by minimizing the mean-squared value of the W’s subject to the above constraints using the Hooke-Jeeves Pattern Moves method. Next, each measure for each control and diabetic subject was adjusted to the mean age of all the subjects in the study (45 yr) as follows where SHFP′ and ΔSBP′ are the age-corrected measures. LSHFP′ was computed as the natural logarithm of SHFP′.

A parasympathetic score for each subject was calculated as the difference between LSHFP′ and the mean of LSHFP′ for the control group, divided by the SD of LSHFP′ for the control group. A sympathetic score for each subject was computed similarly using ΔSBP′. Hence, negative parasympathetic and sympathetic scores, respectively, reflect reduced parasympathetic and sympathetic nervous function. Because it has been found that in diabetic patients there is a degradation in both sympathetic and parasympathetic function (9), we computed a combined autonomic score for each subject referred to as the A score. The A score was determined by summing the parasympathetic and sympathetic scores for each subject and then computing the difference between this sum and the mean of this sum for the control group, divided by the SD of this sum for the control group.

Diabetic subjects were divided into three groups based on their A scores. Specifically, diabetic subjects with A scores greater than −2 were placed in the minimal autonomic neuropathy group, those with scores between −4 and −2 were placed in the moderate autonomic neuropathy group, and those with scores less than −4 were placed in the severe autonomic neuropathy group (see Table1). These cutoff values were predetermined independently of the patient data. Diabetic subjects in the minimal group, by standard measures, have autonomic function indistinguishable from the control group. The moderate and severe groups represent increasing degrees of autonomic neuropathy.

### CSI

*Data acquisition*. The ECG, ABP, and ILV signals were measured noninvasively and were recorded continuously for each subject on a Teac C-71 FM tape recorder (Teac America, Montebello, CA). The ECG signal was measured with a Hewlett-Packard EKG Monitor 78203A (Andover, MA); the ABP signal was measured with a 2300 Finapres Continuous Blood Pressure Monitor (Ohmeda, Englewood, CO); and the ILV signal was measured with a two-belt chest-abdomen inductance plethysmograph (Respitrace System, Ambulatory Monitoring, Ardsley, NY). These signals were recorded first while each subject was supine and then after passive tilting to 60° (relative to the supine posture) by an electrically driven tilt table. Each subject was provided with a period of hemodynamic equilibration before each of the recordings. After each equilibration period, the respiratory signal was calibrated by having the subject alternately fill and empty an 800-ml calibrated Spirobag. Next, ∼8 min of data were recorded. During both recordings, subjects were instructed to breathe on cue according to a randomly spaced sequence of auditory tones. The interbreath times ranged from 1 to 15 s, with a mean of 5 s. This random breathing technique provides broadband respiratory activity while preserving normal ventilation (4). The surface ECG, ABP, and ILV signals for each subject were low-pass filtered (cutoff frequency = 180 Hz) and digitized with a sampling rate of 360 Hz.

*Data analysis*. The 8-min segments of the three digitized signals were analyzed to determine the transfer relations and noise sources in the closed-loop model in the Fig. 1 (see Ref. 18 for data analysis.) The transfer relations, which are presented in their time domain form of impulse response functions, were characterized with the following three parameters: peak amplitude, area, and characteristic time (Fig. 2 and Ref. 18). The perturbing noise sources are presented in terms of their power spectra. The power spectra were characterized by three parameters, which were determined by computing the power in the following three frequency bands: 0–0.50 Hz (total power), 0–0.15 Hz (low-frequency power), and 0.15–0.50 Hz (high-frequency power; see Ref. 18).

*Age correction of parameters*. Parameters characterizing the impulse response functions and perturbing noise sources were age corrected using a method similar to that used to age correct the parasympathetic and sympathetic nervous system measures (see *Grading of subjects*). Each parameter for the control subjects was regressed with age using one of the following two equations
where the*C _{i}
*and

*D*values are unknown parameters and W

_{i}_{parameter}is the age-independent contribution to the parameter. The first equation was used for the area and characteristic time parameters, whereas the second equation was used for the remaining parameters. Constraints were imposed to ensure that the regression equations yield nonnegative values for the characteristic time and peak amplitude parameters, as well as for the spectral powers. The exponents (

*C*

_{3},

*D*

_{3}) also were constrained to be nonnegative except in the case of the peak and area parameters for circulatory mechanics. The age adjustments were made as discussed above.

### Statistical Analysis

Comparisons of the age-adjusted parameters of the impulse response functions and noise power spectra were performed across the subject groups using one-way ANOVA. For those parameters for which the ANOVA test revealed a significant difference (*P* < 0.05), multiple pairwise comparisons between groups were performed using the Tukey honest significant difference for unequal sample sizes test. To make the parameters more normally distributed [as evaluated by normal probability plots of residuals (11)], all parameters, except the area parameter which took on both positive and negative values, were log transformed before statistical analysis.

## RESULTS

In Fig. 1, we present the averaged impulse response functions and power spectra of the perturbing noise sources for data obtained in the tilted posture for each of the subject groups. These averages do not reflect age adjustment, but there were no significant age differences among the groups. In Tables2-4, we present the results of the statistical analyses of the differences across subject groups in the parameters characterizing the impulse response functions HR baroreflex and ILV→HR and the power spectrum of the perturbing noise source N_{HR}. The results indicate that the magnitudes of these impulse response functions and this power spectrum incrementally diminish as one progresses from the control group to the minimal, moderate, and severe autonomic neuropathy groups. The parameters characterizing the circulatory mechanics and ILV→ABP impulse response functions and the power spectrum of N_{ABP} were not statistically different across the groups (except for the characteristic time parameter for circulatory mechanics, which had a positive ANOVA test but no significant pairwise comparison between groups).

## DISCUSSION

In this study, we applied a parametric system identification technique (18) to determine the physiological couplings and perturbing noise sources of a closed-loop model of short-term cardiovascular control. We specifically applied this method to study diabetic patients and control subjects to assess quantitatively and noninvasively the alterations in cardiovascular control mechanisms that result from DAN. The CSI method that we implemented characterizes short-term cardiovascular regulation in terms of the impulse response function of four physiological coupling mechanisms and the power spectra of two perturbing noise sources.

The CSI results are consistent with accepted physiological findings as demonstrated by the averaged impulse responses and the power spectra of the perturbing noise sources of the control group in Fig. 1. The circulatory mechanics impulse response is of the form of a blood pressure wavelet. This impulse response is delayed by ∼0.25 s, which reflects the interval between the R wave of the ECG and the onset of the finger artery ABP pulse. The average pulse pressure of the control group is ∼40 mmHg. The HR baroreflex impulse response decreases in response to a transient increase in ABP and then returns to baseline, which is the expected negative feedback result. The ILV→HR impulse response shows that HR transiently increases in response to inspiration. Furthermore, this increase begins before the impulse inspiration at *time 0*, indicating that HR rises in anticipation of corresponding increases in ILV. This result supports the hypothesis that both respiratory and heart control is neurally coupled within the central nervous system, resulting in HR changes preceding changes in ILV. The ILV→ABP impulse response shows that ABP first abruptly decreases in response to a transient increase in ILV. After ∼2 s, ABP rapidly increases. The immediate decrease in ABP may be attributed to the decrease in intrathoracic pressure with inspiration being transmitted to arterial structures and thus leading to a rapid transient decrease in ABP. Within the next couple of beats, ABP recovers and, in fact, is increased due to increased diastolic filling associated with the decrease in intrathoracic pressure. The N_{HR}power spectrum increases with decreasing frequency. The low-frequency peak in the N_{ABP} power spectrum may correspond to the frequencies of Mayer wave activity (15).

We found that in comparing control subjects with diabetic subjects with minimal, moderate, and severe autonomic neuropathy, as determined by standard autonomic testing methods, the impulse response functions of the autonomically mediated coupling mechanisms (HR baroreflex and ILV→HR) and the power spectrum of the autonomically mediated perturbing noise source (N_{HR}) all progressively decreased in magnitude with increasing autonomic neuropathy across these groups. Importantly, these impulse response functions and the power spectrum showed statistically significant differences even between the control group and the minimal autonomic neuropathy group. These two groups were indistinguishable on the basis of the standard autonomic tests used for grading the extent of neuropathy in the diabetic subjects, suggesting that the CSI method is more sensitive than these standard autonomic tests. Furthermore, comparison of the control group and the minimal autonomic neuropathygroup on the basis of the entire battery of conventional autonomic tests (including the 2 tests used for grading the diabetic patients) by means of unpaired *t*-tests adjusted for multiple comparisons also revealed no statistically significant differences between these two groups (this comparison was not adjusted for age, but there were no significant age differences between the 2 groups).

The impulse response functions of the primarily mechanically mediated coupling mechanisms (circulatory mechanics and ILV→ABP) and the power spectrum of the remaining perturbing noise source (N_{ABP}) essentially showed no differences across the groups. The only marginally significant difference was in circulatory mechanics, which is known to be subject to autonomic modulation.

Other techniques such as heart rate power spectral analysis (2, 13, 14,19), pharmacological measurement of baroreflex gain (16, 17), and nonparametric analysis of the cardiorespiratory coupling (5, 12) have been previously performed to assess autonomic neuropathy in patients with diabetes mellitus. These techniques have shown decreased heart rate variability, baroreflex gain, and cardiorespiratory coupling in diabetics. These techniques, however, may not be as effective as parametric system identification techniques for the direct characterization of physiological control mechanisms. Power spectral analysis of heart rate provides an examination of the fluctuations in heart rate, which is one output of the cardiovascular control system. The heart rate power spectrum does not provide direct information about how heart rate changes in response to varying inputs to the control system. It would seem preferable to determine directly the couplings of the control and feedback loops of the cardiovascular control system. Pharmacological measurement of baroreflex gain provides a measure of the coupling between ABP and heart rate. However, it requires intravenous administration of a vasoactive agent that disrupts normal operating conditions. Nonparametric analysis of the cardiorespiratory coupling also provides a technique for characterizing physiological couplings. However, this technique cannot impose causality, which is essential for consistent least-squares estimation when the data are obtained in a closed loop (18). Parametric system identification techniques can impose causality conditions and, consequently, these techniques can properly identify the feedforward (circulatory mechanics) and feedback (HR baroreflex) couplings between heart rate and ABP, even when the data are obtained in a closed loop.

Although the objective of this study was not to examine in detail the effects of posture (see Ref. 18), we do note that the peak amplitude tends to be diminished and characteristic time prolonged in the autonomically mediated impulse response functions (HR baroreflex and ILV→HR) obtained in the tilted versus supine positions. The reduced amplitude is expected due to parasympathetic withdrawal, and prolonged characteristic time is expected due to the slower reaction time of the sympathetic nervous system, which plays a greater role in mediating these couplings in the tilted position.

In this study, a new method for the assessment of the autonomic neuropathy associated with diabetes mellitus was introduced. The results indicate an incremental diminution in the magnitude of the impulse response functions characterizing the two autonomically mediated couplings (HR baroreflex and ILV→HR) and the power spectrum of the autonomically mediated perturbing noise source (N_{HR}) as one progresses from the control to the minimal, moderate, and severe autonomic neuropathy groups. The difference between the control and minimal autonomic neuropathy group was not identified by the standard autonomic tests. In contrast, the impulse response functions characterizing the primarily mechanically mediated couplings (circulatory mechanics and ILV→ABP) and the power spectrum of the remaining perturbing noise source (N_{ABP}) were essentially unchanged across subject groups.

The present study confirms that system identification can provide a powerful, sensitive, and quantitative means of assessing the extent of DAN in individual patients. The fact that the autonomically mediated impulse response functions and noise source power spectrum were selectively affected in patients with diabetes, which affects neurally mediated cardiovascular control mechanisms, helps confirm that the couplings and perturbing noise sources characterized by system identification truly reflect the underlying physiological mechanisms.

System identification promises to provide an important new tool for the characterization of alterations of cardiovascular regulation caused by a variety disease processes and environmental factors (such as microgravity) and provides a new clinical methodology for guiding therapy. This method is very general and could be used to create a more detailed description of the cardiovascular regulatory state when more signals are made available. For example, measurements of second-to-second values of stroke volume and central pressures and volumes could enable the characterization of physiological mechanisms involved in the control of peripheral vascular resistance and cardiac contractility. The approach of system identification, involving the analysis of fluctuations in multiple biological signals, is a general one that may have applications to the study of regulatory processes from the level of molecules to organ systems to intact organisms to ecosystems.

## Acknowledgments

We thank Christopher Broadbridge for expert technical assistance.

## Footnotes

Address for reprint requests: R. J. Cohen, Massachusetts Institute of Technology, 77 Massachusetts Ave., Cambridge, MA 02139.

This work was supported by a grant from the National Space Biomedical Research Institute, National Aeronautics and Space Administration Grant NAGW-3927, and National Heart, Lung, and Blood Institute Grant 5-R01-HL39291. J. M. Mathias is grateful for support from a SERC/NATO Postdoctoral Fellowship.

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. §1734 solely to indicate this fact.

- Copyright © 1999 the American Physiological Society