Journal of Applied Physiology
Vol. 83, No. 3, pp. 975-993, September 1997
CONTROL OF BREATHING, CIRCULATION, AND TEMPERATURE

Nonlinear systems identification: autocorrelation vs. autoskewness

Michel Sammon and Frederick Curley

Division of Pulmonary and Critical Care Medicine, University of Massachusetts Medical Center, Worcester, Massachusetts 01655

ABSTRACT

INTRODUCTION

METHODS

LINEAR VS. NONLINEAR AUTOCORRELATION

AUTOREGRESSIVE SYSTEMS WITH CORRELATED, NON-GAUSSIAN INPUTS

METASTABLE SYSTEMS WITH THRESHOLD NONLINEARITIES

PERIODIC AND CHAOTIC OSCILLATIONS

DISCUSSION

APPENDIX

ACKNOWLEDGEMENTS

FOOTNOTES

REFERENCES

ABSTRACT

Sammon, Michel, and Frederick Curley. Nonlinear systems identification: autocorrelation vs. autoskewness. J. Appl. Physiol. 83(3): 975-993, 1997.Autocorrelation function (C¹) or autoregressive model parameters are often estimated for temporal analysis of physiological measurements. However, statistical approximations truncated at linear terms are unlikely to be of sufficient accuracy for patients whose homeostatic control systems cannot be presumed to be stable local to a single equilibrium. Thus a quadratic variant of C¹ [autoskewness function (C²)] is introduced to detect nonlinearities in an output signal as a function of time delays. By use of simulations of nonlinear autoregressive models, C² is shown to identify only those nonlinearities that "break" the symmetry of a system, altering the mean and skewness of its outputs. Case studies of patients with cardiopulmonary dysfunction demonstrate a range of ventilatory patterns seen in the clinical environment; whereas testing of C¹ reveals their breath-by-breath minute ventilation to be significantly autocorrelated, the C² test concludes that the correlation is nonlinear and asymmetrically distributed. Higher-order functionals [e.g., autokurtosis (C³)] are necessary for global analysis of metastable systems that continuously "switch" between multiple equilibrium states and unstable systems exhibiting nonequilibrium dynamics.

structural stability; metastability; symmetry breaking ; respiratory failure; congestive heart failure; nonlinear dynamics

INTRODUCTION

RECENT REVIEWS by Bruce and Daubenspeck (2, 3) on time-series analysis of respiratory measurements emphasized the issue of "temporal scaling of respiratory pattern variability." This observation is indicative of the complexity and state dependence of the respiratory control system: metabolism, distribution of brain stem neurotransmitters, blood gas tensions, ventilatory dead space, sensations of dyspnea, cardiac output, pulmonary mechanics, and strength of respiratory muscles can vary with different time constants (4, 5, 7, 9, 10, 13-15, 17, 23, 24, 26-28, 30-32, 35-37). Under conditions of respiratory distress, highly complex ventilatory dynamics may unfold that would not be observable in the unstressed system. Thus comparative analyses of the ventilatory patterns of healthy adults and patients with cardiopulmonary disease are critical toward understanding mechanisms of respiratory failure.

However, such analyses are complicated by the fact that various assumptions of classical linear statistics may not be valid when the control system is under duress and not operating locally to a single equilibrium state (a ventilatory "set point") (19, 22, 29). For example, many investigators have found that variables such as depth or rate of breathing exhibit a positive breath-to-breath autocorrelation (1, 13, 23, 26, 36), a result that has been interpreted as evidence of "memory" within the central pattern generator (1) or "noise" within the peripheral feedback control loop (23). Regardless of physiological interpretations, nonzero autocorrelation implies that a fraction of the signal variance has a nonrandom component (see APPENDIX A). Some underlying assumptions of linear, Gaussian systems analysis are as follows: 1) mean output of the controller and the system equilibrium are equivalent, 2) first- and second-order moments (e.g., mean and variance) are independent, 3) autocorrelation and frequency responses are symmetrically distributed about the mean, and 4) system response to perturbation is independent of the amplitude and direction of the disturbance (i.e., memory is independent of noise). By definition, however, autocorrelation in nonlinear dynamical systems is more positive in some directions than in others; when the distribution is asymmetric, the mean output becomes skewed away from the system equilibrium toward those domains where autocorrelation is more positive. Furthermore, when the deviation from linearity is continuous, changes in the signal-to-noise ratio (SNR) alter not only the mean and variance but also the skewness and correlated fraction of the system output.

Therefore, statistical analyses of the output from a nonlinear control system must take into consideration that such parameters (mean, variance, autocorrelation, and skewness) can be interdependent quantities, particularly as the system becomes less stable. Here, dynamic interactions between the first- and second-order moments are referred to as autoskewness and quantified using a quadratic variant of the autocorrelation function (C¹), which will be referred to as the autoskewness function (C²). Reviews of nonlinear systems analysis using functional series expansion can be found elsewhere (16, 19, 20, 29); our analysis adopts a pulse amplitude modulation that restricts estimates of higher-order statistical cumulants (29) or Wiener kernels (39) to a single time-delay axis. This greatly simplifies visualization and statistical testing of nonlinearities in an output signal. Probabilistic (8, 25, 29) and geometric (6, 11, 21, 34) analyses of various nonlinear systems are relied on to outline the distinctions of linear vs. nonlinear autocorrelation and local vs. global stability. Inasmuch as deviations from linearity can be specified for these systems, they will be used to generate test signals to demonstrate the utility of C² analysis in quantifying the differential responses of a system to weak vs. strong perturbations, thus identifying "symmetry-breaking" nonlinearities (21, 31). These analytic techniques will be shown to have diagnostic value for distinguishing the types of nonlinearities typically seen in the ventilatory patterns of patients with cardiorespiratory dysfunction. Higher-order functionals [e.g., autokurtosis (C³)] are necessary for analysis of systems with structural instabilities (6, 34), e.g., where dynamics vary bimodally relative to a threshold nonlinearity or intermittently "switch" between multiple equilibria.

METHODS

Ventilatory patterns for the seven subjects of Figs. 7, 9, and 12 were acquired by inductance plethysmography (Respitrace Plus, NonInvasive Monitoring, W. Palm Beach, FL), with tidal volume calibrated via spirometry. These cases were selected for the range of variability and system memory in their respiratory patterns to demonstrate issues related to the statistical methods. Ventilatory records were collected during an afternoon nap from patients A and B (57- and 74-yr-old men) with mild upper airway obstructions and audible snoring. Five patients were monitored with Respitrace on days of successful "T-piece" weaning from mechanical ventilation in the University of Massachusetts Medical Center Intensive Care Unit (ICU): patients C, D, and E (67-, 74-, and 82-yr-old men) had congestive heart failure (CHF) with left ventricular ejection fractions <50% and marked Cheyne-Stokes respiration (CSR); patients F and G (a 73-yr-old man and a 57-yr-old woman) were recovering from pneumonia. Approval was obtained from our institutional review board for human subjects for off-line analysis of these measurements, which were recorded as part of standard clinical protocols. The algorithms listed in APPENDIXES A AND B are written in the pseudo-C language of the software Matlab (MathWorks, Natick, MA). Inasmuch as the "autofunctionals" are estimated by standardizing each record according to its own probability density function, the analytic results are independent of the absolute volume calibration of the Respitrace signal. Those signals were also not subjected to any filtering before analysis.

Statistical Method and Definition of Terms

t  1

t  2

t  D

t  1

t  D

t  d

t  d

(1)

(2)

(3)

(4)

t  0

t  

t  

LINEAR VS. NONLINEAR AUTOCORRELATION

A statistical method for detecting signal nonlinearities is effective only if it quantifies these phenomena independently of the system equations (presumed unknown) that actually generated the signal (20). Consequently, two generalized first-order, nonlinear autoregressive systems are evaluated, for which the deviations from linearity can be specified using a single control parameter (b). These maps will first be used to familiarize readers with some statistical properties of nonlinear dynamical systems and then as signal generators to evaluate the C² test for linearity. Accurate interpretations of C^1,2 necessitate a comparison of systems that have some properties in common but also differ in other aspects.

Nonlinear Autoregressive Systems (Maps)

t  1

t  1

Parameter Space (a and b)

For parameter values outside its stable domain, system 1 grows unbounded to ± (domain A of Fig. 1). However, for |a| > 1, system 2 has regions where periodic and chaotic oscillations occur (domain B); elsewhere, X_t grows unbounded or converges to an "apnea" equilibrium located at the origin (a > 1, b < a); this equilibrium is important for global analysis of system 2.

Our present focus is on systems that are stable locally to a single equilibrium; within a neighborhood of (a,b) = (0,0), the qualitative dynamics depend on where the control parameters lie (I-IV). These quadrants are not unique to systems 1 and 2; to a degree of approximation, local stability of any nonlinear control system can be categorized similarly. With regard to system identification, estimates within the C^1,2 plane will correspond directly with the (a,b) plane. In the parameter space the null hypothesis of a memoryless system corresponds with a single point (a,b) = (0,0), whereas the null hypothesis of linearity is a restriction to the line (a,0). For local analysis, the statistical question is: How far must the systems be perturbed from the point/line in order for the C^1,2 tests to reject the null hypothesis? The flip-side question [Is the system sufficiently close to (a,b) = (0,0) that C^3,4,... can be disregarded?] is deferred until sections on global analysis.

Perturbation Analysis

In Fig. 2, left, the systems are configured such that the slope of the map at is positive (a = 0.5); in response to small perturbations in either direction, the systems exhibit monotonic decays back to equilibrium (overdamped response). Because our primary interest is in analysis of ventilatory patterns, systems with negative nonlinear coefficients (b < 0) will be referred to as hypoventilatory, indicating that the slope of the map is more positive below (+++ in Fig. 2); as a consequence, these systems have impulse responses that are (on average) below that expected from a linear system. In contrast, the simulations in Fig. 2, right, are configured such that the slope is negative at (a = 0.5); the systems respond to small perturbations with underdamped oscillations above and below . With positive curvature (b > 0), the autocorrelation becomes more positive above the equilibrium (+++ in Fig. 2), and the responses to larger perturbations become skewed in that direction (hyperventilatory). Deviation from the linear response is strictly negative for overdamped hypoventilatory systems, whereas underdamped hyperventilatory systems tend to oscillate about the linear perturbation response.

Stochastic Analysis

Because the parameter space for system 1 is "skewed" (Fig. 1), deviation between its mean and equilibrium depends on linear and nonlinear coefficients. The orthogonality of a and b for system 2 results in the difference depending only on b. However, the mean and equilibrium are equivalent only in the linear case for both systems (b = 0) On average, outputs of nonlinear control systems tend toward domains where autocorrelation is more positive: below for b < 0 (hypoventilatory) and above for b > 0 (hyperventilatory).

The coefficient of determination relates the percent increase in variance between input and output signals due to the memory of the system: r² = 1  ²/²_x. For stable, first-order linear systems, r² = a². For a fixed value of a, deviations from linearity generally result in an increase in the correlated fraction of the signal variance (r² > a²), although slight reductions can occur for systems with nonzero second derivatives (r² < a², see arrow for system 2 in Fig. 3).

In response to a symmetrically distributed input, the output of a linear system is also symmetrical about its mean ( = 0). However, skewness of the output of a nonlinear system varies in proportion to its deviation from linearity and in the direction where autocorrelation is more positive: hypoventilatory nonlinearities perturb skewness toward more negative values, whereas hyperventilatory nonlinearities make  more positive.

In Fig. 3, deviations between the equilibrium and mean output due to nonlinear autocorrelation depend on the SNR of the system; i.e., increasing stochastic variability amplifies the deviation for systems 1 and 2. Changes in skewness and coefficient of determination also vary with respect to the SNR, but only in systems with nonzero second derivatives (system 2). Note that an asymmetric P_X (|| > 0) does not imply that the dynamics of X are nonlinearly autocorrelated any more than nonzero variance implies that X is linearly autocorrelated. These "lag-zero" statistics have no diagnostic value with respect to the nature of the system memory but are essential for scaling the C^1,2 functions so that M^1,2 can be quantified (APPENDIX A).

Sensitivity of C^1,2 to SNR

Relation Between C^1,2 and (a,b) Planes

System Nonlinearities at Multiple Time Delays

With regard to system identification, a minimum of five parameter estimates must be specified to (empirically) characterize a nonlinear oscillation: the slope and curvature at the smallest time lag, C^1,2₁, a feedback delay (D), and slope and curvature of the vector field at this delay, C^1,2_D, i.e., two points on the C^1,2 plane, with D representing time to traverse between them. Depending on curvature of the nonlinear feedback, the minimum or maximum of C² generally occurs at the "true" delay (arrows in Fig. 6), whereas the minimum of C¹ occurs at a higher time lag.

Uniform vs. Nonuniform Sampling

Fortunately, the debate is often unnecessary, because analytic results using the two estimation methods are generally not contradictory, as in Fig. 7, where 30-min records of breath-by-breath minute ventilation (I) and breath duration (TT) are shown for four patients who exhibited oscillations in their ventilatory pattern. The discrete and continuous estimates of C^1,2 for I are shown for each patient. When the coefficient of variation of TT is <1/4 (Fig. 7A), the two estimates of C^1,2 are nearly identical, because the breath-by-breath intervals are relatively constant over time. When TT is more variable (coefficient of variation > 1/4), greater differences between estimates are seen; C² dynamics tend to dampen out more quickly for discrete estimates.

Which algorithm should be used for breath-by-breath data? Although more computative, the continuous version is preferred, particularly for patients in whom TT is highly variable. If comparisons are to be made between patients with different respiratory rates, it is important that C^1,2 be estimated on equivalent time scales (e.g., seconds). Furthermore, results using this algorithm may be more representative of dynamics of arterial blood gases, which vary continuously over time.

However, breath-by-breath rhythms intrinsic to respiratory control may also exist that are asynchronous with respect to absolute time (31, 32). Inasmuch as some feedback to the control system is discrete by nature (e.g., breath-by-breath neuromechanical feedback), algorithm 1 may be preferred, particularly when high-frequency oscillations in a variable are seen (i.e., C¹ < 0 at small time lags). Furthermore, if C^1,2 is to be used for building a nonlinear autoregressive model for analyzing breath-by-breath dynamics, it is advantageous to use the discrete estimate.

Relation of C^1,2 to Return Maps

Local vs. Global Stability of Ventilatory Control

For apparent reasons, many analysts would categorize the ventilatory record of Fig. 9 as nonstationary; certainly, an estimate of C^1,2 over the 200-min interval would be statistically unstable (APPENDIX B). Still, the analytic question remains: To what degree are the observed variations autonomous? Measurements of other physiological variables would suggest a high degree of determinism: patterns 2 and 3 coincide with reduced heart rate and oxymetric arterial saturation, as well as lower mean ventilation (local C⁰). Cycle times of the CSR oscillations increase during intervals 2 and 3, as might be expected with lower cardiac output. Whereas cardiorespiratory interactions are clearly evident, interdependence of these systems with mechanisms controlling sleep state are likely involved as well. Such physiological couplings would result in a synthesized system having a potentially very complex geometry and high number of independent degrees of freedom. As discussed below, if geometry of a control system is such that dX_t/dX_t  1 at is near unity, its dynamics become structurally unstable; i.e., small perturbations can result in large changes in the qualitative behavior of its outputs (6), similar to that seen in Fig. 9. Thus, whereas local analysis characterizes system dynamics within a neighborhood of stationary equilibrium point, global analysis is concerned with domains farther from , i.e., nonequilibrium (or multiequilibria) dynamics. For linear systems, the two perspectives are equivalent, inasmuch as system response is invariant with respect to the magnitude of perturbation. For nonlinear systems, there generally exist thresholds beyond which estimates of C^1,2 are inadequate to describe the dynamics. One approach toward global analysis is a piecewise comparison of parameters (C^0,1,2) estimated locally with respect to time, as in Fig. 9; a statistically more rigorous approach is to estimate higher-order functionals (C^3,4) over the whole sample interval, as estimates of Cⁿ reflect the sensitivity of local estimates of C^n  1 to local variations in C⁰.

AUTOREGRESSIVE SYSTEMS WITH CORRELATED, NON-GAUSSIAN INPUTS

Simulations of previous sections employed autoregressive models with parameters ( and b) that were independent of a symmetrically distributed "white" noise input (). This section addresses interpretations of C^1,2 for systems with a slowly varying equilibrium () or linear stability (a), i.e., when  is itself a linearly autocorrelated variable ("colored" noise) having additive/multiplicative interactions with state variables (X). In this case, the input C¹ is not an impulse at lag zero but, instead, decays linearly to zero. A geometric explanation may also be helpful as to why an asymmetric input to a linear system does not result in a false-positive C² conclusion that X_t is from a nonlinear system. The system equations can be generalized as X_t = h(X_t  1,X_t  2,...,X_t  D,_t,_t  1,...,_t  M) + _t. Orthogonal (e.g., linear) cases can be decomposed as h = f(X_t  1,...,X_t  D) + g(_t  1,...,_t  M), where the poles/peaks of the system transfer function are determined by the autoregressive component f, with an impulse response consisting of an infinite exponential decay (Fig. 2); the "moving-average" component (g) is a finite impulse response filter that determines the transfer function's zeros/valleys. Unless measurements of or presumptions on  are made, f and g cannot be discriminated using C^1,2 analysis of X_t; the "system" is h = f(g) driven by a zero-mean white noise that is orthogonal to f and g. Historical perspectives on Wiener filtering and Wold's decomposition theorem for autoregressive moving-average (ARMA) systems can be found elsewhere (20, 22).

For the simulations of Fig. 10, a linear first-order moving average variable was first generated: _t = 0.5 (u_t + u_t  1), where u is identically distributed white noise, N(0,1). The symmetrically and continuously distributed  was then resampled as follows:  = 0 for || < 0.48,  = 0.1 for  < 0.48,  = 0.15 for  > 0.48. The result is a positively skewed ( = 0.4), linearly correlated signal that varies among three discrete points [0.1,0,0.15], having a probability density of P = [0.25,0.5,0.25]. This piecewise constant  was input to the following five systems where AR is autoregressive and MA is moving average. Figure 10, top, shows the instantaneous response of each map to this three-point input. For the linear cases A and B,  "moves" along the diagonal without altering the slope of the map (dX_t/dX_t  1). Consequently, autoskewness of the output X is statistically zero; the fact that  and X are non-Gaussian has no bearing on this result. For the nonlinear cases, the slope at becomes less positive as the equilibrium is moved up the diagonal for case C (hypoventilatory nonlinearity) and becomes more positive for case D (hyperventilatory); the nonlinear case E is hypoventilatory as well, i.e., the input simultaneously perturbs upward while lowering dX_t/dX_t  1 at . In all cases, C²₁ varies in proportion to d(dX_t/dX_t  1)/, without regard to the (causal) model structure. C² of X_t will also be nonzero where the input is a nonlinear moving-average signal, e.g., as X_t  1 varies in proportion to _t  1, a substitution results in case C. What about multiplicative interactions between X_t and _t? If the summation is linearized (k = 0,1), this case is similar to case E. However, the C² test for linearity would be rejected only for |a|  1 and |1/b| > max{||,|X|}; systems with similar threshold nonlinearities are discussed later.

t  1

t  1

Linear, Time-Varying Approach for Evaluating Autoskewness

t,

t  

If parameter estimates for patient E in Fig. 9 are examined from a similar perspective, it should be apparent that localized estimates of C^1,2 varied with level of I (C⁰). Higher-order functionals (Cⁿ, n > 2) are necessary to approximate the geometries of systems that produce such interdependencies.

METASTABLE SYSTEMS WITH THRESHOLD NONLINEARITIES

Quadratic approximations of systems are not globally stable, in that there exists a perturbation of sufficient magnitude to break the intersection of the estimated map with the diagonal, resulting in a divergence toward ±; for cases C, D, and E, the divergence was transient, because _t was not held beyond the threshold level indefinitely (arrows). Dynamical systems become structurally unstable as the slope of their map along the diagonal approaches unity (6); this section demonstrates how the presence of cubic nonlinearities bound such systems from diverging to infinity. Their analysis requires the next highest-order functional (C³), which is referred to as the autokurtosis function, inasmuch as it quantifies the degree to which dynamics perturb system density (P_X) toward a bimodal distribution. Estimation and statistical testing of C^3,4,... are consistent with that of C² (APPENDIX A).

Simulations of previous sections presumed that system dynamics took place locally to a single equilibrium. For linear and weakly nonlinear control systems, this is a valid assumption; a straight line with slope <1 can intersect the diagonal only at a single point, . Systems having stronger nonlinearities can also exhibit unstable nonequilibrium dynamics such as periodic and chaotic oscillations (discussed later). However, an alternative type of instability must also be considered in which multiple equilibrium states emerge, each of which can be locally (but not globally) stable. Although such "metastable" respiratory control has been proposed to occur during sleep (37), multiequilibria may also unfold as components of respiratory failure. A first-order cubic map will be used to illustrate these topics