## Abstract

Physiological rhythms, including respiration, exhibit endogenous variability associated with health, and deviations from this are associated with disease. Specific changes in the linear and nonlinear sources of breathing variability have not been investigated. In this study, we used information theory-based techniques, combined with surrogate data testing, to quantify and characterize the vagal-dependent nonlinear pattern variability in urethane-anesthetized, spontaneously breathing adult rats. Surrogate data sets preserved the amplitude distribution and linear correlations of the original data set, but nonlinear correlation structure in the data was removed. Differences in mutual information and sample entropy between original and surrogate data sets indicated the presence of deterministic nonlinear or stochastic non-Gaussian variability. With vagi intact (*n* = 11), the respiratory cycle exhibited significant nonlinear behavior in templates of points separated by time delays ranging from one sample to one cycle length. After vagotomy (*n* = 6), even though nonlinear variability was reduced significantly, nonlinear properties were still evident at various time delays. Nonlinear deterministic variability did not change further after subsequent bilateral microinjection of MK-801, an *N*-methyl-d-aspartate receptor antagonist, in the Kölliker-Fuse nuclei. Reversing the sequence (*n* = 5), blocking *N*-methyl-d-aspartate receptors bilaterally in the dorsolateral pons significantly decreased nonlinear variability in the respiratory pattern, even with the vagi intact, and subsequent vagotomy did not change nonlinear variability. Thus both vagal and dorsolateral pontine influences contribute to nonlinear respiratory pattern variability. Furthermore, breathing dynamics of the intact system are mutually dependent on vagal and pontine sources of nonlinear complexity. Understanding the structure and modulation of variability provides insight into disease effects on respiratory patterning.

- nonlinear dynamics
- vagus
- Kölliker-Fuse nucleus
- entropy
- breathing pattern variability

recent research has demonstrated that the human respiratory pattern has characteristics of a nonlinear dynamic system (3–5, 57). Moreover, these characteristics may be considered statistical biomarkers because they are sensitive, both to physiological state (4, 5, 50) and to the presence of several clinical pathologies, including obstructive sleep apnea syndrome, Cheyne-Stokes respiration, extubation failure after mechanical ventilation, asthma, carotid atherosclerosis, and panic disorder in humans (6, 26, 34, 37, 54–56). Taken together, the promise of these previous studies and the concept that cardiorespiratory patterns underlie physiological homeokinesis provided the motivation to develop an understanding of respiratory pattern variability and the network mechanisms that contribute to its genesis.

The breathing pattern is generated by the respiratory central pattern generator (rCPG) in the brain stem and is modulated by peripheral sensory inputs (for review, see Ref. 17). Mechanosensory inputs from the airway and lungs are transduced by pulmonary stretch receptors and transmitted to the nuclei tractus solitarii (nTS) by vagal afferent fibers. Second-order nTS neurons that receive this input project to pontomedullary components of the cardiorespiratory network and affect respiration via mechanisms that control phase switching and duration (35, 41). Removal of this mechanosensory feedback loop with vagotomy affects neural activity within the pontomedullary respiratory network, as well as interactions between pontomedullary respiratory groups (10). In particular, vagotomy increases not only respiratory-modulated pontine activity, but also pontine modulation of the respiratory pattern, such that the pons serves as an “internal vagus” controlling phase-switching mechanisms in a manner similar to vagal feedback (10, 41, 42, 49). Furthermore, blockade of *N*-methyl-d-aspartate receptor (NMDAR) currents in the Kölliker-Fuse nuclei (KFn) after vagotomy induces apneusis in several species (20, 23, 33), suggesting that pontine NMDAR-dependent transmission is critical for the control of inspiratory timing.

Vagotomy affects respiratory patterning, not only by prolonging respiratory phase durations, but also by reducing pattern variability (28). The reduction in pattern variability was measured as a decrease in the coefficient of variation of the phase durations. However, the coefficient of variation does not quantify the nonlinear properties of respiratory variability: specifically, variability that results from deterministic nonlinear relationships between a given point in the respiratory pattern and subsequent points. Nevertheless, we view these analyses as complementary. The coefficient of variation of respiratory phase durations measures variability of respiratory timing about a mean, irrespective of position in the time series. In contrast, our analysis of respiratory pattern variability measures changes in the variability of the instantaneous motor output and depends on the temporal sequence of points. Physiologically, analysis of pattern variability may provide insight into the interplay of peripheral inputs to the pontomedullary cardiorespiratory network that leads to the emergence of nonlinear dynamics, such as the breath-to-breath interdependence observed in the respiratory pattern (15, 22, 32, 36).

In studies that examined nonlinear properties of the breathing pattern, Sammon and Bruce (43) demonstrated that vagotomy reduces the correlation dimension of the respiratory pattern. Correlation dimension is an information theory-based index of pattern variability that quantifies the number of degrees of freedom required to capture the dynamics of the system on an attractor. These authors concluded from both experimental studies and model simulations that the respiratory pattern of the anesthetized, spontaneously breathing rat can be represented by a nonlinear limit-cycle attractor, which can be driven to exhibit more complex dynamics associated with higher nonlinear variability via vagal-dependent, mechanosensory feedback. In the present study, we sought to identify or characterize the determinants of respiratory pattern variability using mutual information, sample entropy (SampEn) and correlation dimension, in combination with iterated amplitude-adjusted Fourier transform (iAAFT) surrogate data sets. The additional metrics were chosen because of their previous application and potential utility in a clinical setting. Our approach should be viewed as complementary to the work of Sammon and Bruce (43).

To assess the role of mechanosensory feedback and its pontine processing, we quantify variability in the breathing pattern with vagi intact, after vagotomy, after dorsolateral (DL) pontine MK-801 microinjection, and after both vagotomy and MK-801. We use autocorrelation and coefficient of variation of inspiratory burst statistics to quantify linear stochastic variability of the respiratory pattern. Furthermore, we employ two statistical techniques: mutual information, a measure of statistical dependence in a time series (see Fig. 1); and SampEn, a measure of pattern predictability that arises from both linear and nonlinear correlations (see Fig. 2). To distinguish between linear and nonlinear sources of variability using these techniques, we compared the mutual information or SampEn of the original data set with that of iAAFT surrogate data sets that preserve the amplitude distribution and autocorrelations of the original data, while eliminating any nonlinear correlations (for review, see Ref. 45). We conduct this analysis by comparing points separated by multiple time delays (τ) from neighboring points, up to and including those separated by one cycle length to quantify the components of respiratory pattern variability on the basis of their characteristic time scale. Finally, to confirm the presence of nonlinear deterministic dynamics within the time series, we estimated the correlation dimension. We hypothesize that *1*) mechanosensory input of vagal origin contributes to both linear and nonlinear variability in respiratory motor activity; and *2*) signaling via NMDARs in the KFn is critical for the expression of vagal-dependent nonlinear variability.

## METHODS

Experiments were performed under protocols approved by Case Western Reserve University's Institutional Animal Care and Use Committee and were performed with strict adherence to all American Association for Accreditation of Laboratory Animal Care International, National Institutes of Health, and National Research Council guidelines.

### Surgical Preparation

Adult male rats (Sprague-Dawley, obtained from Zivic Miller, Pittsburgh, PA, *n* = 11, 350–375 g) were anesthetized with urethane (1.2 g/kg ip). Anesthetic depth was assessed regularly by testing the withdrawal reflex: an absence of blood pressure increase and heart rate response to a paw pinch. Supplemental doses of urethane (15% of original dose, intravenous) were administered as necessary. Once the animal was at a surgical plane of anesthesia, we exposed the vagi bilaterally and the trachea through a midline incision. We tied sutures loosely around the vagal nerve bundles and included the aortic depressor nerve within the loop.

Body temperature was monitored with a rectal probe and maintained between 38 and 39°C with a heating pad. We cannulated the following: *1*) a femoral artery to monitor blood pressure and to sample arterial blood; *2*) a femoral vein to administer fluids and supplemental anesthetic; and *3*) the trachea to maintain a patent airway in the anesthetized animal. Animals breathed spontaneously throughout the experiment. Diaphragmatic electromyographic (DiaEMG) activity was recorded as an index of the respiratory pattern. Bipolar Teflon-coated tungsten wire electrodes with bared tips were inserted percutaneously into the diaphragm, adapting the technique of Basmajian and Stecko (2a). The wire leads were twisted, and distance between the two electrodes was minimized to attenuate noise and artifacts related to the electrocardiogram. DiaEMG activity was filtered (band-pass 3 Hz to 1 kHz) and amplified (AC preamplifier; model P511K, Grass Instruments). The raw DiaEMG activities were stored in a computer sampled at 2 kHz via an analog-to-digital sampling system (Power1401, CED). To obtain an “integrated” signal, the data were full-wave rectified and smoothed (time constant, 100 ms).

### Experimental Protocol

Before initiating data collection, an arterial blood sample was collected for analysis (ABL5, Radiometer, Copenhagen, Denmark). The postsurgical acidic pH was corrected with an appropriate volume of 8.4% (wt/vol) sodium bicarbonate solution. Analysis of arterial blood gases was conducted at later time points during the experiment to ensure that respiratory drive was stable (7.43 ± 0.02, arterial pH; 186 ± 7 Torr, arterial Po_{2}; 38.8 ± 1.4 Torr, arterial Pco_{2}; and 25.2 ± 1.4 mmol/l, arterial bicarbonate). For baseline recording of DiaEMG, room air was supplemented with 100% O_{2} (2–3 l/m). To assess the independent and joint contributions of vagal feedback and NMDAergic transmission in the KFn to nonlinearity of respiratory motor output, we divided the experiments into two groups. In *protocol 1* (*n* = 6), after recording baseline activity, we transected the vagi before proceeding with bilateral microinjection of MK-801 into the DL pons. In *protocol 2* (*n* = 5), we microinjected MK-801 bilaterally before vagal transection.

### Microinjection of MK-801 into the KFn

Bilateral craniotomies were made to expose the surface of the cerebellum for KFn microinjections. Subsequent local administration of glutamate (10 mM), MK-801 (15 mM), and fluorescent microspheres was made by pressure microinjection (Veradyne). Injection volumes were measured by observing the drop in the meniscus of a triple-barreled micropipette (tip diameter, 50 μm) using a surgical scope equipped with a calibrated reticule. After the baseline DiaEMG signal was recorded, the respiratory-responsive sites within the DL pons were mapped using the apneic response observed upon glutamatergic stimulation (20–40 nl; Fig. 3*A*). Effective coordinates were found at 0–0.2 mm caudal to lambda, 2.9–3.2 mm lateral to midline, and 5–6 mm below the cerebral surface. After identification of the KFn, MK-801 (100 nl) and fluorescent microspheres (100 nl) were injected. The contralateral DL pons was injected with MK-801 without mapping the response to glutamate to minimize the interval before and after injections. Recordings were then made in the presence of MK-801 before or after vagotomy. After the experiment, animals were perfused with heparinized saline followed by 4% (wt/vol) paraformaldehyde/PBS. Brains were removed and postfixed in 4% paraformaldehyde/PBS containing 20% sucrose. Coronal slices (50 μm) of the brain stem were cut on a freezing microtome (Leica Jung Histoslide 2000) and mounted for post hoc identification of the injection sites (Fig. 3, *B* and *C*).

### Data Analysis

To test our hypotheses, we analyzed at least three stationary epochs, each 100 s in duration, of the integrated DiaEMG signal with vagi intact, after vagotomy, in the presence of MK-801 alone, and after vagotomy in the presence of MK-801. We chose to analyze this integrated DiaEMG, as this signal has been shown to be physiologically equivalent to recordings of respiratory flow (13, 27). Furthermore, application of the present analytic approach to the raw DiaEMG signal would reflect the fast time scale nonlinear dynamics associated with motoneuron spiking. Instead, we are interested in the pattern variability associated with the slow dynamics of the rCPG that is revealed in the integrated DiaEMG signal. We applied Spike2 software (CED, Cambridge, UK) for inspiratory burst detection and calculated the mean, standard deviation, and coefficient of variation of phase durations, and inspiratory burst statistics, including peak height (Pk Ht), root mean square amplitude and area. The linear statistical relationship between sampled points in the pattern was measured by the autocorrelation function for time delays/lags (τ) from neighboring points to points separated by one cycle length, as well as by the coefficient of variation of the amplitude distribution (Fig. 4). We calculated mutual information and SampEn, together with surrogate data testing (Figs. 1, 2, 5⇓⇓–8), to investigate nonlinear variability in the pattern.

#### Surrogate data testing.

To test for the presence of nonlinear dynamics, we compared computations using the original data set with those obtained from the surrogate data sets (*n* = 19 for each epoch) that were generated using the iAAFT method (25, 52). The iAAFT surrogate data sets were constructed such that both the autocorrelation function (linear correlations) and the amplitude distribution of the original data set were preserved [Figs. 5 and 6 (25, 45)]. Thus the nonlinear correlations were eliminated in the surrogate data sets, while maintaining the linear correlation structure and amplitude distribution of the data. A representative tracing of a surrogate data set before and after vagotomy is shown in Fig. 5*B*. Subsequent nonlinear analyses (mutual information, SampEn, and correlation dimension) were applied to both the original and surrogate data sets (Figs. 5⇓–7). Although we interpret differences between the surrogates and original data sets as indicating the contribution of deterministic nonlinear variability to these dynamic measurements, these differences could also arise from a stochastic non-Gaussian dynamic system.

#### Autocorrelation.

Autocorrelation functions for each epoch were computed across a range of τ to characterize linear variability in the data set using available MatLab functions. Autocorrelation quantifies the strength of linear correlations in the time series on the basis of their characteristic time scale τ. Because the respiratory pattern is highly periodic, autocorrelation, mutual information, and SampEn analyses all showed characteristic relationships with the cycle length. This led us to report the τ as a function of cycle length to facilitate comparison across animals with different breathing frequencies.

#### Mutual information.

Mutual information is a measure of statistical dependence in the data set (19, 47) and quantifies how knowledge of a coordinate *x*(*t*) reduces the amount of uncertainty associated with a time-shifted coordinate *x*(*t* + τ) (Fig. 1*A*). Mutual information (MI) is defined as
*P*[*x*(*t*), *x*(*t* + τ)] is the joint probability of *x*(*t*) and *x*(*t* + τ) (this probability can be obtained from the joint probability distribution function, Fig. 1*B*, *bottom*), and *P*[*x*(*t*)] and *P*[*x*(*t* + τ)] are the marginal probabilities of *x*(*t*) and *x*(*t* + τ) for the original or time-shifted time series, respectively. In practice, the computation begins by constructing a return map (Fig. 1*B*, *top*) from the original and time-shifted time series (Fig. 1*A*). The joint probability distribution function (Fig. 1*B*, *bottom*) is generated by dividing the area of the return map into bins to determine the relative frequency that a point in the time series is in in a given bin. We used a uniform number of bins (40 bins), where the size of each bin was scaled by the variance of the data set. The individual probability distributions are determined by taking the sum of each column or row of the joint probability distribution for the *x*(*t*) or *x*(*t* + τ) time series, respectively. The individual probability distributions of the original and time-shifted time series data sets are identical. Finally, mutual information, a statistic of the joint probability distribution, can be directly computed using *Eq. 1*. Because the breathing pattern is highly periodic (autocorrelation for 1 cycle length is 0.99 ± 0.00, *n* = 10 epochs), we chose to compute mutual information for τ from one sample (neighboring points separated by 5 ms) to one cycle length to generate the mutual information functions before and after vagotomy in Figs. 5 and 6. Additionally, because mutual information measures statistical dependence from both linear and nonlinear sources, we compute mutual information of both the original and surrogate data sets to test for the presence of nonlinear dependence in the time series (e.g., if mutual information is greater in the original data set relative to the average mutual information of the surrogate data sets, then the analysis suggests that the breathing pattern contains nonlinear dependence for that τ). We used the mean difference between the mutual information of the original and surrogate data sets across all τ up to one cycle length as an index of nonlinear dependence (NLdi) in the time series.

#### SampEn.

SampEn is a measure of self-similarity in a time series and provides a robust computational method for quantifying predictability in short time series data (40). Briefly, for every point in the time series, a template of *m* points (each separated by a given τ) was constructed (Fig. 2*A*, red and green squares within the black box, *m* = 2 and τ = 10). The time series was then searched to identify possible template matches (Fig. 2*A*, each series of a red and green point). A pair of points was considered a template match if its member points were within a tolerance *r* of the corresponding point in the original template. Next, the template of *m* points was extended to include an additional point (Fig. 2*A*, addition of a black square within the black box). The time series was then searched for matches of this extended template of *m* + 1 points (Fig. 2*A*, each series of a red, green, and black point). This searching algorithm was repeated for all possible starting points to compute the total number of *m*-point matches and (*m* + 1)-point matches. Informally, SampEn is a permutation of the ratio of the number of matches for the (*m* + 1)-point templates (Fig. 2*A*, green boxes) to the number of matches for the *m*-point templates (Fig. 2*A*, red and green boxes), excluding self-matches. Formally, SampEn is computed using the following equation
*A ^{m}*(

*r*,τ) is the probability that two (

*m*+ 1)-point sequences match,

*B*(

^{m}*r*,τ) is the probability that two

*m*-point sequences match,

*N*is the number of points in the epoch,

*m*defines the template length,

*r*is the tolerance (a fraction of the standard deviation of the amplitude envelope of the epoch), and τ is the time interval/delay between points in a template or matching sequence. Thus a higher SampEn indicates either a fewer number of (

*m*+ 1)-point template matches or a greater number of

*m*-point template matches (i.e., a lower SampEn indicates increased temporal pattern regularity in the time series and, as such, greater predictability and lesser complexity).

In practice, each epoch was decomposed into a series of contiguous 10-s windows, each composed of 2,000 data points. The SampEn of each window was computed across the epoch. Traditionally, SampEn is computed with the following parameters: *m* = 1 or 2 points, *r* = 0.2 × SD without incorporating a τ, e.g., τ = 1 sample (30, 38, 39). However, oscillatory signals tend to have a high degree of linear correlation at a τ of 1, so SampEn may not adequately capture nonlinear contributions to pattern complexity at this time scale. To investigate complexity on different characteristic time scales, we computed SampEn across multiple τ as suggested in Ref. 25 from neighboring points (5 ms) to those separated by one cycle length (Fig. 2*B*). In addition, we computed the uncertainty associated with the SampEn computation according to the following equation: *U* = SD_{data}·*t*_{B − 1, 0.975}/*B*, where SD_{data} is the standard deviation of the time series, *t*_{B − 1, 0.975} is the upper 2.5th percentile of a *t*-distribution with *B* − 1 degrees of freedom, and *B* is the number of template matches (40). Only stationary epochs were included in this analysis, and SampEn and its uncertainty were averaged across all windows for a given τ. For each τ, a binomial test was used to determine whether the SampEn of the surrogates was statistically different (α = 0.05) from the SampEn of the original data set for all windows in a given epoch. When the SampEn of the original data was statistically different from the SampEn of the surrogate data sets, we used the distance between the SampEn of the surrogates and that of the original data as a measure of the amount of nonlinear complexity present in the epoch for that τ. If the SampEn of the surrogates was statistically indistinguishable from the SampEn of the original data set, the contribution of nonlinear determinants of complexity was taken as 0 bits for that τ. We used the mean of these differences across all τ as a global index of nonlinear complexity (NLci) irrespective of time scale.

#### Correlation dimension.

Correlation dimension quantifies the active degrees of freedom necessary to replicate the dynamics of the system in a reconstructed phase space (attractor). We reconstructed the nonlinear system using the time-delayed embedding technique, with the τ chosen as the first minimum of the mutual information function, as suggested by Fraser and Swinney (19). Correlation dimension was computed using the Grassberger-Procaccia algorithm for embedding dimensions, *m* = 3–15 (21). In all cases, we used a Theiler correction window of 40, which was always greater than the τ used for the reconstruction to correct for autocorrelation effects (51). To confirm the findings obtained with mutual information and SampEn analyses, we computed the correlation dimension of both the original and surrogate data sets (Fig. 7).

All data were expressed as means ± SE, and, unless stated otherwise, we applied a one-way ANOVA with repeated measures to determine the statistical significance of differences. If significant, we used a Student-Newman-Keuls post hoc test to identify specific differences.

## RESULTS

To test our hypotheses, we analyzed stationary epochs of respiratory motor activity with vagi intact, after vagotomy alone, after microinjection of MK-801 alone, and after both vagotomy and MK-801 microinjection. Representative tracings with vagi intact (*left*), after vagotomy (*middle*), and after vagotomy and MK-801 microinjection (*right*) are shown in Fig. 5*A*. The KFn was identified physiologically by observing the expiratory prolongation in response to local glutamate application (Fig. 3*A*). The injections sites were found predominantly in the KFn, a critical component of the pontine respiratory group. In most cases, injection volumes were sufficient to perfuse the full rostrocaudal extent of the KFn.

Vagotomy alone was associated with a significant increase in both inspiratory (Ti) and expiratory duration (Te) by 37 and 211%, respectively (Ti with vagi intact, 0.30 ± 0.01 s vs. after vagotomy, 0.41 ± 0.01 s, *P* < 0.01; Te with vagi intact, 0.38 ± 0.03 s vs. after vagotomy, 0.80 ± 0.04 s, *P* < 0.01; Fig. 4). With respect to patterning of respiratory motor activity, vagotomy significantly increased Pk Ht by 1.13-fold (*P* < 0.05) and inspiratory burst area by 1.62-fold (*P* < 0.01). Vagotomy tended to decrease the coefficient of variation of respiratory phase durations, but this tendency was not significant for the group.

Bilateral microinjection of MK-801 into the KFn with vagi intact resulted in a mild, but significant increase in Ti, with little change in the Te (Ti with vagi intact, 0.30 ± 0.01 s vs. after MK-801, 0.40 ± 0.05 s, *P* < 0.01). The coefficient of variation of Ti and Te tended to increase after MK-801 microinjection. However, these changes were not found to be significant for the group. MK-801 administration was also associated with a significant 1.38-fold increase in inspiratory burst area (*P* < 0.05).

The combination of vagotomy and MK-801 resulted in an apneustic breathing pattern, characterized by significant increases in both Ti and Te (Ti with vagi intact, 0.30 ± 0.01 s vs. after vagotomy and MK-801, 0.54 ± 0.03 s, *P* < 0.01; Te with vagi intact, 0.38 ± 0.03 s vs. after vagotomy and MK-801, 0.61 ± 0.05 s, *P* < 0.01). Moreover, the failure of the inspiratory off-switch after vagotomy and MK-801 was highlighted by the additional significant increase in Ti from either vagotomy or MK-801 alone (*P* < 0.01 in both cases). The combination of MK-801 and vagotomy also induced a tendency for an increase in the coefficient of variation of respiratory phase durations. Finally, with vagotomy and bilateral MK-801, Pk Ht, root mean square amplitude, and area increased significantly by 1.22-fold (*P* < 0.01), 1.23-fold (*P* < 0.01), and 2.18-fold (*P* < 0.01), respectively. In all states analyzed, the coefficient of variation of all inspiratory burst statistics measured did not change significantly across treatments (Fig. 4*D*). Together, these results are consistent with the reported literature on the effects of vagotomy and blockade of DL pontine NMDAergic transmission (20, 23, 33).

### Analysis with iAAFT-surrogate Data Testing Demonstrates the Presence of Nonlinear Variability During Spontaneous Anesthetized Breathing

Representative examples of our analytic approach for *protocols 1* and *2* are shown in Figs. 5 and 6, respectively. Our approach to quantify nonlinear variability critically depends on the ability of the iAAFT surrogate data sets to maintain both the autocorrelation structure and amplitude distribution of the original time series to adequately test the null hypothesis that the time series is obtained from a linear Gaussian stochastic system. Representative traces from a surrogate time series are shown below their original counterparts in Fig. 5, *A* and *B*. From the traces, it is apparent that the amplitude distribution and general periodicity of the motor pattern is preserved in the surrogate data sets. Furthermore, autocorrelation functions are shown for both the original (solid trace) and surrogate (shaded trace) data sets in Figs. 5*C* and 6*C*. As expected, the autocorrelation functions of the surrogate data sets were not significantly different from that of the original data. Thus differences between metrics applied to surrogate vs. original data sets reveal the contribution of nonlinear or non-Gaussian variability to these measures.

With vagi intact, at τ = 1, a point at which linear correlation is maximal (Fig. 5*C*, *left*), the average mutual information of the surrogate data sets was not significantly different from the original data, suggesting that linear Gaussian stochastic correlation alone accounts for the relationship between adjacent points. By contrast, as the τ increased from one sample (5 ms) to one cycle length, the amount of linear correlation in the time series decreases (Fig. 5*D*, *left*), and differences between the mutual information of the original and surrogate data sets indicated the presence of significant nonlinear dependence in the time series over many time scales (Fig. 5*D*, *left*). Similar results were observed using SampEn as the statistical metric for comparison of the original and surrogate data sets. With vagi intact, at very short τ, the SampEn of the original data was not found to be significantly different from that of the surrogate data sets (Fig. 5*E*, *left*). As the τ was increased, differences between the SampEn of the surrogate and original data sets emerged, suggesting the presence of significant nonlinear complexity. Together, these observations demonstrate that analysis of the system dynamics at small τ leads to a conclusion that the dynamics can be adequately explained by linear Gaussian stochastic variability alone. However, nonlinear breathing variability is revealed at higher τ, including those that are optimal for the reconstruction of the phase space using the time-delayed embedding technique (19).

### Vagal Feedback Contributes to Nonlinear Respiratory Dynamics

With vagotomy alone, the difference between the mutual information of the surrogate and original data sets decreased significantly at many τ, suggesting that vagotomy decreased nonlinear dependence in the breathing pattern, regardless of time scale (Fig. 5*D*, *middle*). For the group, vagotomy alone significantly decreased the NLdi by 29% (with vagi intact, 0.31 ± 0.05 bits vs. after vagotomy, 0.22 ± 0.05 bits, *P* < 0.01; Fig. 8*A*). With vagotomy, the difference in SampEn of the surrogate and original data sets also decreased significantly over many time scales (Fig. 5*E*, *middle*). For the group, the NLci fell significantly by 26% (with vagi intact, 0.31 ± 0.02 bits vs. after vagotomy, 0.23 ± 0.04 bits, *P* < 0.01; Fig. 8*B*). Together, these observations support the conclusion that vagal feedback contributes to the nonlinear respiratory dynamics observed in the intact animal.

### Blockade of DL Pontine NMDAergic Transmission Decreases Nonlinear Breathing Dynamics

Representative autocorrelation, mutual information, and SampEn curves are shown before and after MK-801 microinjection in Fig. 6. With MK-801 microinjection alone, we observed a decrease in the difference between the mutual information of the original and surrogate data sets, suggesting a loss of nonlinear dependence (Fig. 6*B*). Across all τ, the NLdi decreased significantly by 35% (with vagi intact, 0.31 ± 0.05 bits vs. with MK-801, 0.20 ± 0.04 bits, *P* < 0.01; Fig. 8*A*). The magnitude of the difference between the SampEn of the surrogate and original data sets also fell with MK-801 microinjection (Fig. 6*C*). For the group, the NLci decreased significantly by 42% (with vagi intact, 0.31 ± 0.02 bits vs. with MK-801, 0.18 ± 0.03 bits, *P* < 0.01; Fig. 8*B*). Together, these results support the conclusion that, even with intact vagal feedback, DL pontine NMDAergic transmission is also required for the expression of eupneic nonlinear respiratory dynamics.

### Vagal Feedback and KFn Are Mutually Dependent Sources of Nonlinear Respiratory Dynamics

With both vagotomy and MK-801 microinjection, the degree of nonlinear dependence in the breathing pattern was significantly different from baseline (with vagi intact, 0.31 ± 0.05 bits vs. with vagotomy and MK-801, 0.20 ± 0.03 bits, *P* < 0.01; Fig. 8*A*). The NLci was also significantly less than baseline (with vagi intact, 0.31 ± 0.02 bits vs. with vagotomy and MK-801, 0.21 ± 0.02 bits, *P* < 0.01; Fig. 8*B*). However, the NLdi and NLci with vagotomy and MK-801 were not significantly different from the NLdi and NLci of either the vagotomy-alone or MK-801-alone groups, suggesting that vagotomy occluded any further MK-801-dependent changes in nonlinear variability in *protocol 1* and vice versa in *protocol 2*.

### Correlation Dimension Estimates Confirmed the Existence of Nonlinear Respiratory Dynamics

While the surrogate data testing above suggested the existence of nonlinear dynamics within the rCPG, our analysis can only conclude that the dynamics of the respiratory time series cannot be explained by a linear Gaussian stochastic dynamic system. To provide further evidence that our analysis with mutual information and SampEn with iAAFT surrogate testing revealed the presence of nonlinear deterministic variability, we computed the correlation dimension in conjunction with iAAFT surrogate data testing with vagi intact, after vagotomy, after MK-801 microinjection, and after both vagotomy and MK-801. The correlation dimension was taken to be the value for which the slopes of the correlation integrals converge for at least 2 log lengths. The correlation dimension estimates the number of degrees of freedom required to reconstruct the dynamics on the attractor. Noninteger values of correlation dimension suggest a fractal characteristic of the attractor and are indicative of the presence of nonlinear system dynamics. Representative correlation integrals from *protocol 1* are plotted in Fig. 7*A*. In all conditions analyzed, the respiratory motor activity had a correlation dimension estimate of ∼1.2 (Fig. 8*C*). To further confirm that these noninteger estimates of correlation dimension were accurate indications of the existence of nonlinear dynamics within the breathing pattern, we also attempted to estimate the correlation dimension of the iAAFT surrogate data sets (Fig. 7*B*). As expected, there was no meaningful convergence of the slopes of correlation integrals because the surrogate data sets were adequately represented by a linear Gaussian stochastic system. Together, these observations support the conclusion that anesthetized, spontaneous breathing dynamics arise from a nonlinear oscillator.

## DISCUSSION

Both vagal-dependent and vagal-independent nonlinear determinants of variability were evident from a thorough analytic investigation using information theory-based measures of variability in conjunction with iAAFT-surrogate data analysis. Nonlinear respiratory pattern variability was attenuated to a similar degree by bilateral vagotomy or bilateral MK-801 microinjection in the KFn. When applied in conjunction, these perturbations did not result in any additional changes in nonlinear dependence or complexity, irrespective of sequence, suggesting that these circuit elements are mutually involved in the genesis of mechanosensory-dependent nonlinear breathing variability.

### Limitations of the Analytic Approach

The surrogate data sets used in this investigation are the principal feature of our analytic approach that distinguishes it from previous studies on nonlinear respiratory dynamics. Surrogate data sets were generated using the iAAFT method, whereby surrogate data sets preserve both the amplitude distribution and autocorrelation functions of the original data (25, 52). Such surrogate data sets provide a means of hypothesis testing for the presence of nonlinear structure within the time series using metrics such as mutual information and SampEn (45), which are sensitive to both linear and nonlinear dynamics. Significant differences between these metrics when applied to surrogate and experimental data sets indicate the presence of either stochastic non-Gaussian or deterministic nonlinear variability. To further validate that the observed differences between the SampEn or mutual information of the surrogate and original data sets represented nonlinear structure, we examined the convergence of correlation integrals for both the original time series and its iAAFT surrogate data sets (Fig. 7). Correlation integrals showed convergence for the original but not surrogate time series, confirming the existence of nonlinear deterministic structure in the time series. However, the correlation dimension (∼1.2) was not dependent on vagal feedback or intact NMDAergic transmission in the KFn. The failure of our correlation dimension estimates to decrease with vagotomy contrasts with previous estimates (43). These authors, however, eliminated two-thirds of samples contributing to expiration from the time series collected after vagotomy to avoid biasing the analysis by the change in the respiratory period from baseline. This selective downsampling may account for the differences observed in the correlation dimension estimates of the present study. Other differences between the two studies include the type of respiratory pattern recording (DiaEMG vs. inductance plethysmography) and the animal's posture during the experiment. Nonetheless, our use of iAAFT surrogates, in conjunction with correlation dimension estimates, did confirm the presence of nonlinear variability in both states. Moreover, the ability of mutual information and SampEn, combined with iAAFT surrogate data testing to distinguish differences in nonlinear variability, despite differences in the respiratory period, highlights the tractability of these metrics for application in a clinical setting compared with analysis using correlation dimension.

### Sensory Inputs Influence Nonlinear Respiratory Dynamics

Maintenance of pulmonary ventilation is critical for blood-gas (oxygen and carbon dioxide) homeostasis. As such, one finds that the network responsible for generating the motor patterns associated with this behavior is relatively insensitive to perturbation and usually operates within a regime of limit-cycle dynamics (43). However, ventilation can be coupled with autonomic rhythms for efficient and optimized gas exchange. This coupling may be accomplished through the interconnection between the respiratory and cardiovascular control systems and/or through sensory afferent feedback loops within the network. From a computational perspective, feedback loops have the potential to generate constrained nonlinear deterministic phenomena (53). Furthermore, in a phenomenological model of vagal feedback on respiration in rats, Sammon and Bruce (43) demonstrated, using bifurcation analysis, that an increase in the pulmonary stretch receptor feedback delay parameter above one-eighth of the respiratory period was sufficient to replicate the increase in complexity associated with the presence of vagal feedback in their experimental studies. Thus we hypothesize that nonlinear dynamics within the CPG arise from the interactions of the peripheral inputs feeding into the network. However, the network mechanism by which this feedback circuit leads to a bifurcation of the complexity in the motor output remains to be described.

Recent work has specifically examined the role of various sensory inputs in generating nonlinear respiratory dynamics in humans and animal models. With respect to vagal mechanosensory feedback, while our work and that of Sammon and Bruce (43) supports a role for vagal-dependent mechanosensory feedback, Akay (2) found that vagotomy did not result in any change in respiratory complexity in neonatal and juvenile piglets. These authors used approximate entropy to quantify respiratory complexity and did not investigate the effect of τ on computational results. Restricting analysis to a short time scale (τ = 1) may not capture the nonlinear dynamics of time series data that have strong linear correlations (autocorrelation) or limit-cycle behavior that masks the nonlinear dynamics operating over longer time scales (25). Accordingly, we found that nonlinear complexity or dependence in the breathing pattern does not exist at a τ of one sample (5 ms with a sampling frequency of 200 Hz). Thus we conclude that vagal mechanosensory feedback contributes to nonlinear respiratory dynamics when properly quantified.

With respect to chemosensory feedbacks, Fiamma and colleagues (18) found that, with hypercapnia, respiratory pattern complexity increased using the largest Lyapunov exponent and Kolgorov-Sinai entropy and the noise limit to quantify pattern complexity. While the authors also did not use a τ in their computations, the data were downsampled to 5 Hz before analysis. Compared with our results, this is well within a regime in which nonlinear respiratory dynamics exist, but analysis with surrogate data testing would have improved the conclusiveness of their findings. In contrast, Akay (1) did not find any influence of hypercapnia on respiratory complexity in neonatal and juvenile piglets. Again, this group's approach using approximate entropy without a τ might explain their conclusion, even though nonlinear variability was present in the respiratory pattern.

Peripheral chemosensory feedbacks have also been shown to influence nonlinear breathing pattern variability in animal models. Restricting their analysis to inspiratory bursts, Chen and colleagues found that inspiratory pattern complexity decreased during NaCN-induced gasping in rats in vivo (7). The authors' use of multiscale entropy is similar to our own in that it quantifies nonlinear variability over many time scales or τ values. However, the exact contribution of nonlinear vs. stochastic influences remains unclear without surrogate data testing. Nonetheless, the authors also demonstrated that activation of the pre-Bötzinger region by dl-homocysteic acid microinjection mimicked the gasping-induced reduction in complexity, suggesting that reorganization of the rCPG in response to strong chemoafferent activation can account for changes in respiratory dynamics. Taken together, these observations suggest that chemosensory and mechanosensory feedbacks contribute to nonlinear respiratory dynamics observed during eupnea.

### Pontine Networks Also Contribute to Nonlinear Respiratory Dynamics

The microinjection studies presented in this paper demonstrated that NMDAR blockade could selectively decrease the nonlinear components of pattern variability, despite both the trend for increased variability in phase durations and the presence of intact vagal feedback. Moreover, when vagotomy occurred after MK-801 microinjection, there was no additional change in nonlinear pattern dependence or complexity, suggesting that vagal-dependent nonlinear variability is also dependent on the presence of intact DL pontine synaptic activity.

Neurons in the DL pons are well situated to be a central neural correlate for the nonlinear component of respiratory pattern variability. The pons is not essential for inspiratory rhythmogenesis, but modulates inspiratory and expiratory timing. Yet how the pons modulates timing remains a complex question. Even though anatomical studies show intense projections to the medullary ventrolateral respiratory column, the pons appears to have a diffuse tonic influence on neurons in the column because few pontine neurons influence activity of medullary respiratory neurons directly, as revealed by cross-correlation histograms (46). Furthermore, most DL pontine neurons are only weakly modulated by the respiration discharging in tonic or phase-spanning fashion (16). The pons may also influence respiratory pattern by modulating afferent input. While the initial synapse of baro-, chemo- and mechanosensory feedback occurs within the nTS, secondary fibers project to the pons, and, reciprocally, pontine neurons can gate peripheral feedback in the nTS (11). Finally, pontine neurons can show various types of plasticity in response to peripheral feedback (12, 24, 48). The respiratory system displays memory-like behavior whereby perturbations of peripheral inputs have effects that last for many breaths (8, 14, 22, 32, 48). The memory-like behaviors observed after stimulation of the vagal or superior laryngeal nerves depend on an intact pons (36). Taken together, these data support a role for pontine neurons in providing a modulatory, coupling influence to the pontomedullary respiratory network that, by virtue of intrinsic synaptic delays and recurrent connectivity, may contribute significantly to the nonlinear component of breathing variability.

### Nonlinear Respiratory Pattern Variability as a Biomarker of Physiological Status

While many studies have demonstrated that nonlinear dynamics of the heart rate time series are indicative of pathophysiology due to the sensitivity of these measures to the balance between parasympathetic and sympathetic tone, a growing literature supports the application of similar analytic approaches to the respiratory pattern as a biomarker of pathological status. Veiga and colleagues (54) demonstrated that the approximate entropy of airflow recordings from asthmatic patients decreases with the severity of the disease. Kaikamakis and colleagues (26) demonstrated that the severity of obstructive sleep apnea syndrome could be accurately distinguished using an analysis of approximate entropy, detrended fluctuation analysis, and largest Lyapunov exponent of respiratory flow signals. Moreover, such approaches can be combined with attractor neural networks to generate more sophisticated algorithms for pathology detection. For instance, Weinreich and colleagues (55) were able to distinguish obstructive sleep apnea syndrome and Cheynes-Stokes from normal respiration with a sensitivity of 90% and specificity of 91% using spectral entropy as one of the test criteria. Similar analytic techniques are also relevant for weaning from mechanical ventilation. Recent findings from Schmidt and colleagues (44) and from Papaioannou and colleagues (37) both demonstrated that the complexity of the respiratory airflow pattern increases during spontaneous breathing trials, which preceded successful weaning from mechanical ventilation. Moreover, in analyzing the interbreath interval time series from patients during spontaneous breathing trials, White and colleagues (56) similarly observed that extubation failure was associated with decreased respiratory pattern complexity. Additionally, Letellier and colleagues (31) have shown that recurrence plots can be used as an indicator of patient-ventilator asynchronisms that contribute to patient discomfort. As such, recurrence plots may be useful to tune ventilator parameters in the intensive care unit. Because control of respiratory and cardiovascular physiology is coupled, the utility of nonlinear respiratory variability as a biomarker is not limited to respiratory pathologies. Recently, Mangin and colleagues (34) have shown that analysis of respiratory pattern variability is decreased in the presence of carotid atheroscelerosis and returns after the occlusion was removed. In an animal model of stroke, Koo and colleagues (29) demonstrated that respiratory pattern complexity is increased after stroke. Together, these observations provide ample motivation to investigate the underlying mechanisms behind nonlinear respiratory variability in health and disease.

## GRANTS

We also gratefully acknowledge that this work was supported by National Institutes of Health (NIH) Grants HL-080318, NIH HL-087377 (T. E. Dick), HL-007913 (R. R. Dhingra), the Veterans Affairs Research Service (F. J. Jacono), and NIH NS-057815 (I. A. Rybak).

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

## ACKNOWLEDGMENTS

The authors thank Dr. Kendall Morris for insightful comments on this manuscript.