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

Symmetrical Case (b = 0)

Symmetry Breaking (b  0)

These simulations demonstrate that the autoskewness function discriminates only those nonlinearities that break system symmetry and alter the mean and skewness of its output. Thus a finding of |C²| > 0 is a sufficient, but not necessary, condition for concluding that a (stationary) signal is nonlinearly autocorrelated; the same can be said for C³. It should be evident from Fig. 11 why a positive cubic nonlinearity would result in a system that oscillates to ± in response to small perturbations; because of their reliance on bilateral negative feedback, C³ tends to be negative for homeostatic control systems. A finding of |C³| > 0 does not necessarily imply that the system is "metastable" (although the converse is true), only that the curvature is asymmetric with respect to C⁰. If cubic memory is significant (M³  0), it does imply that a model possessing a threshold nonlinearity is the best representation of the system at this order of approximation.

Structural Stability of Ventilatory Dynamics

The metastable geometry of Fig. 11 consisted of two locally stable equilibria (zero-dimensional attractors) with diametric curvature separated by an unstable (repellor). The same concept applies to systems possessing attractors/repellors of higher dimension (e.g., 1-dimensional curves). The ventilatory control of patient E in Fig. 9 may be an example of this, where dynamics alternated between Cheyne-Stokes patterns having C² functions with different signs.

PERIODIC AND CHAOTIC OSCILLATIONS

Simulations of previous sections relied on the presence of a time-varying input to perturb the systems from equilibrium to identify components of its impulse response; Cⁿ would become singular for _t = 0, as the deterministic component of the density function (P_X) was a zero-dimensional fixed point (). Once systems enter their periodic/chaotic domains, the attractor dimension is positive and often non-Euclidian (i.e., a noninteger fractal), resulting in oscillations that are self-sustaining (6); estimates of Cⁿ primarily reflect the system attractor rather than its response to perturbation. Whereas higher-order functionals are generally needed to distinguish between attractors, the geometric principles of autocorrelation, skewness, and kurtosis are consistent with those already discussed.

System 2 is useful for generating a broad range of attractors having very complex distributions. The system of Fig. 13 has an elliptical cycle for a₁ = +0.5 and a parabolic density function; however, as X_t narrows toward the diagonal at high values, P_X is negatively autoskewed. As a₁  +1, X_t approaches an orbit that intersects  = 0; this type of Silnikov attractor has been observed experimentally in rats and is associated with nonlinear feedback from vagal mechanoreceptors (31, 32). Amplitude and cycle time of the larger orbit vary according to the proximity with which X_t cycles near  = 0; i.e., higher density at X = 0 in the present correlates with a longer tail in P_X in the future (which is concisely positive autoskewness).

t  1

t  1

t  1

t  3

t  3

t  1

For a₁ > 0, trajectories follow a continuous path in the phase (X_t vs. X_t  1) and the C^1,2 and C^3,4 planes; for a₁ < 0, the trajectory jumps discontinously, e.g., between two closed loops for a₁ = 0.6, 12 points for a₁ = 0.8, and about a Henon-type attractor as a₁  1 (6). Whereas mappings from the phase to the Cⁿ space are generally consistent (i.e., closed curves to closed curves, discrete points to discrete points), we do not know whether attractors that are fractal within the phase space (a = 0.96) are also fractal in the parameter space. Discrimination between the period-12 (a₁ = 0.8) and chaotic (a₁ = 0.96) oscillations is poor using Cⁿ estimates, because their distributions are so close; in the presence of even a small amount of noise, a large sample size would be needed to statistically distinguish between them.

What type of attractors result in higher-order memories, i.e., M^3,4,... > 0? Recall that C¹ is undefined for X_t = constant, and C² is singular for all binary processes. The proximity with which an attractor approaches such states determines how strongly C^3,4 are perturbed; C^3,4 deviate farthest from 0, where X_t approaches a saddle point (a₁ = 0.96) and a saddle cycle (a₁ = 0.6) in Fig. 13. Higher-order nonlinearites are generally associated with very-low-frequency components that switch X_t between states over time (APPENDIX B). The order of system memory (highest n where Mⁿ > 0) is related to the number and symmetry of deterministic modes that can be deconvolved from P_X (e.g., peaks for equilibria, parabolas for cycles). In general, odd functionals discriminate the modes; even Cⁿ quantify their relative symmetry (APPENDIX A).

DISCUSSION

Nonlinear variants (Cⁿ) of the autocorrelation function C¹ were introduced as analytic tools for identifying nonlinearities within the output of a control system. With primary focus on autoskewness (C²), the analysis was tested using simulations of several nonlinear systems driven by Gaussian and non-Gaussian distributed noise. Analysis of ventilatory patterns of patients with various degrees of cardiorespiratory dysfunction demonstrated the necessity of quantifying the linear (n = 1) and nonlinear components when analyzing such signals. Where statistical testing indicated the presence of significant autocorrelation, the hypothesis of linearity was rejected at the 99% confidence level, indicating the presence of significant nonlinearities at various time delays. Observations that greater persistence in C¹ coincides with proportionally greater C² memory are consistent with the notion that, as the set point becomes less stable in the local sense (elevated M¹), quadratic nonlinearities have greater influence on the ventilatory dynamics (elevated M²). In general, elevated Mⁿ indicates greater persistence in C^n + 1, with the sequence M¹,M²,...,M^N converging to zero (although not necessarily monotonically). We view estimation of Cⁿ as a critical first step in time-series analysis, because these autofunctionals have implications for so many statistical properties of the system.

Stochastic vs. Deterministic Components of Density Functions

t  

1,...,n  1

Any Cⁿ test that concludes that a time series is linearly/nonlinearly autocorrelated can be checked using a surrogate X_t generated by randomizing the sequence of the original X_t. The density function of the surrogate X is identical to the original P_X, so that identical basis vectors are generated for estimating Cⁿ. However, the surrogate Cⁿ should be memoryless for all n.

Similarly, any Cⁿ test (n > 1) that detects the presence of nonlinear autocorrelation can be verified by comparing results found for the original X_t with those for surrogate data generated by phase randomization (12). Briefly, this algorithm consists of 1) calculating the fast Fourier transform (FFT) of X_t, 2) randomizing its phase components while holding its magnitude constant, 3) making the FFT self-adjoint, and 4) calculating the inverse transform back to the time domain. The result is a real-valued time series (X_t) that has the same autocorrelation function as the original signal but with nonlinear components that are now uncorrelated. Inasmuch as the inverse FFT produces Gaussian-distributed variables, the surrogate P_X is generally different from the original P_X. Any statistical test of X_t that rejects the null hypothesis of linearity should then result in the opposite conclusion when applied to X_t. When surrogate data were generated for the ventilatory patterns analyzed in Figs. 7, 9, and 12 (which were concluded to be nonlinearly correlated), surrogate C^2,3 were statistically zero.

If Cⁿ (n > 1) is nonzero for X_t generated by either of the surrogate methods described above, it may be indicated that the original signal X_t was statistically unstable relative to the length of its sample interval (APPENDIX B).

Dynamic Changes in Mean, Skewness, and Coefficient of Determination

From linear statistics, the concept of variance having random and nonrandom components is familiar to most systems analysts, with the coefficient of determination (r²) reflecting the correlated fraction of the signal variance (25, 29). For nonlinear systems with nonzero second derivatives, r² depends also on the SNR, with system memory generally increasing with level of noise (23). Less familiar, however, is the concept that the mean and skewness of the output signal can also be viewed in terms of deterministic and random components. Whereas the system equilibrium is the primary determinant of the mean output, deviations in the mean from occur in systems possessing asymmetrically distributed nonlinear autocorrelation.

Thus, for patients with strong symmetry-breaking nonlinearities in their ventilatory control, an injection of noise within the system can induce relative changes in mean I (as well as in r² and ). Although the presence of such asymmetries can be detected using autoskewness analysis, we are unaware of a way to reliably quantify the magnitude of dynamic hypo- or hyperventilation (in l/min) using estimates of C^1,2.

Relation of Cⁿ to Return Maps

"Can Periodic Breathing Have Advantages for Oxygenation?"

Structural Instabilities vs. Time-Varying Dynamics

Pulse Amplitude vs. Pulse Width Modulation for Kernel Estimates

t + 

t + 

t  s

t + 

t + 

t + 

t  s

t  s

t  s

Differential Response to Small- vs. Large-Impulse Perturbations

Other Nonlinear Analyses

Summary

ACKNOWLEDGEMENTS

This research was funded by The Whitaker Foundation.

FOOTNOTES

Address for reprint requests: M. Sammon, Div. of Pulmonary and Critical Care Medicine, UMass Medical Center, 55 Lake Ave., N., Worcester, MA 01655.

Received 22 July 1996; accepted in final form 5 May 1997.

APPENDIX A

t  

n  1

APPENDIX B

t  

t  

t  

t  1

t  1

t  5

t  1

t  5

t  1

t  1

t  5

t  5

t  1

t  5

REFERENCES