## Abstract

Cyclic ventilatory instabilities are widely attributed to an increase in the sensitivity or loop gain of the chemoreflex feedback loop controlling ventilation. A major limitation in the conventional characterization of this feedback loop is the need for labor-intensive methodologies. To overcome this limitation, we developed a method based on trivariate autoregressive modeling using ventilation, end-tidal Pco_{2} and Po_{2}; this method provides for estimation of the overall “loop gain” of the respiratory control system and its components, chemoreflex gain and plant gain. Our method was applied to recordings of spontaneous breathing in 15 anesthetized, tracheostomized, newborn lambs before and after administration of domperidone (a dopamine D_{2}-receptor antagonist that increases carotid body sensitivity). We quantified the known increase in hypoxic ventilatory sensitivity in response to domperidone; controller gain for O_{2} increased from 0.06 (0.03, 0.09) l·min^{−1}·mmHg^{−1} to 0.09 (0.08, 0.13) l·min^{−1}·mmHg^{−1}; median (interquartile-range). We also report that domperidone increased the loop gain of the control system more than twofold [0.14 (0.12, 0.22) to 0.40 (0.15, 0.57)]. We observed no significant changes in CO_{2} controller gain, or plant gains for O_{2} and CO_{2}. Furthermore, our estimate of the cycle duration of periodic breathing compared favorably with that observed experimentally [measured: 7.5 (7.2, 9.1) vs. predicted: 7.9 (7.0, 9.2) breaths]. Our results demonstrate that model-based analysis of spontaneous breathing can *1*) characterize the dynamics of the respiratory control system, and *2*) provide a simple tool for elucidating an individual's propensity for ventilatory instability, in turn allowing potential therapies to be directed at the underlying mechanisms.

- periodic breathing
- apnea
- chemoreflex
- loop gain
- lung

ventilatory control instabilities in the form of cyclic apnea or periodic breathing are important in a variety of pathological conditions, including Cheyne-Stokes breathing in congestive heart failure (25, 59, 6), obstructive sleep apnea in adult humans (28, 60, 54), as well as in adults at altitude (31, 5, 11, 10) and in neonates (49, 53, 58). Although the specific mechanisms underlying each condition may vary, ventilatory instability can result from increases in the ventilatory sensitivity to hypoxia/hypercapnia (controller gain), in the efficiency of gas exchange (plant gain), or in the circulatory delay between lungs and chemoreceptors (28, 5, 7). The concept of loop gain has been increasingly emphasized because it integrates each of these physiological components into a single value that describes the stability of the respiratory control system (see Fig. 1*A*). A high loop gain generally describes a system that is intrinsically less stable, whereas a low loop gain describes a more stable system.

The role of loop gain in understanding respiratory control (and in particular, its link to instability) has been recognized for many years (12, 28). However, the clinical potential of loop gain has been limited by the fact that it depends on complex dynamic interactions of central and peripheral chemoreceptor control loops, on circulatory delays, and on the efficiency of CO_{2} and O_{2} exchange in the lung. These individual factors are cumbersome and impractical to measure directly (21). Indirect measures for assessing the risk of periodic breathing include proportional assist ventilation (38), pressure support ventilation (17), and pseudorandom binary stimulation using CO_{2} (22). All existing techniques to characterize the loop gain of the respiratory control system require interventions that alter breathing patterns and blood gases during either quiet wakefulness or sleep. By contrast, in our study the same goal is achieved by taking advantage of spontaneous fluctuations in ventilation V̇e and in respiratory gas tensions Pco_{2} and Po_{2} (23, 9).

The current study assessed spontaneous fluctuations in ventilation, end-tidal CO_{2}, and end-tidal O_{2} using a trivariate autoregressive modeling technique to derive estimates of controller and plant gains for the feedback loops controlling both CO_{2} and O_{2} and thus the overall loop gain of the control system. We tested our method using data from spontaneously breathing anesthetized lambs, before and after administration of domperidone, a drug previously shown to increase the dynamic hypercapnic and hypoxic ventilatory responses, i.e., the associated controller gains (19). Crucially, if our method has validity, it must yield controller gains that are comparable to published estimates of dynamic hypercapnic and hypoxic ventilatory responses in the lamb. Further, the method would need to reveal an increase in overall loop gain, driven predominantly by an increase in the O_{2}-specific loop gain, following domperidone administration. In addition to passing these tests, we show that spontaneous fluctuations in ventilation are determined by the loop gain of the respiratory control system (29, 41, 16) and that our analysis closely predicts the cycle duration of experimentally induced periodic breathing in lambs. Thus the concept of loop gain provides a unifying framework that links the natural variability in ventilation during spontaneous breathing to the cycle duration of experimentally induced periodic breathing.

## METHODS

Our model quantifies the inter-relationships among spontaneous fluctuations in breath-to-breath values of ventilation, end-tidal CO_{2}, and end-tidal O_{2} by characterizing the “transfer path” between each pair of variables (see Fig. 1*B*). The transfer path between O_{2} and ventilation defines the controller gain for O_{2} and similarly for CO_{2}. Likewise, the O_{2} and CO_{2} plant gains are equal to the transfer path functions between ventilation and O_{2} and CO_{2}, respectively. The product of controller and plant gains for O_{2}, and similarly for CO_{2}, are the respective loop gains (LG_{O2} and LG_{CO2}). Thus our model quantifies the two individual feedback loops, one for the feedback control of O_{2} and one for control of CO_{2}, The overall loop gain of the system is equal to the sum of LG_{O2} and LG_{CO2} (as demonstrated in appendix b). We note that loop gain is a frequency dependent quantity (28). In the current study, where duration is represented by the number of breaths, cycle frequency is therefore described in units of cycles per breath. For example, an oscillation with cycle duration (or period) of 10 breaths corresponds to frequency of 0.1 cycles/breath.

### Experimental Setup and Protocol

Using our model, we quantified the impact of the administration of domperidone, a dopamine D_{2}-receptor antagonist that increases carotid body sensitivity to O_{2} and CO_{2} (19). The impact on the control system transfer paths, and thus on loop gain, was assessed using direct measurement of the respiratory variables associated with spontaneous breathing in newborn lambs. Animal data were obtained from previous experiments performed at Monash University, Australia, by a number of the present authors with other collaborators (19). All surgical and experimental procedures conformed to the guidelines of the National Health and Medical Research Council of Australia and had the approval of the Standing Committee in Ethics in Animal Experimentation of Monash University.

The surgical preparation of the animals, and the variables measured and derived in this study, have been described in our previous publication (19). Briefly, 15 newborn lambs (aged 10–20 days) were given a loading dose of ketamine hydrochloride (5 mg/kg Ketamil 100 mg/ml; ILEUM Veterinary Products) via a nonocclusive catheter (Becton-Dickinson; Intracath 19GA) inserted percutaneously into the left jugular vein. A bolus of α-chloralose (80 mg/kg) was then given, followed by a continuous infusion at 20 mg·kg^{−1}·h^{−1}. The adequacy of anesthesia was checked by regular stimulation of the inner and outer canthus of the eye and by monitoring heart rate and blood pressure. If needed, supplemental doses of α-chloralose were delivered. Once the animals were anesthetized, measurements of respiratory flow, end-tidal Pet_{CO2}, and end-tidal Po_{2} (Pet_{CO2} and Pet_{O2}) were made. Respiratory flow was measured using a Hans Rudolph (3500A) pneumotachograph. Gas concentrations were measured using a Morgan 901 MK2 CO_{2} analyzer and an Ametek S-3A/I O_{2} analyzer. Nonocclusive catheters (Datamasters, Victoria, Australia) placed in the right carotid artery and jugular vein were used for blood gas sampling. Flow and gas tension signals were acquired at a sampling frequency of 400 Hz. In this study, we used 5–10 min of spontaneous breathing data (V̇e, Pet_{CO2}, and Pet_{O2}) before and ≥10 min following the intravenous administration of domperidone.

In our previous study, we experimentally measured the ventilatory sensitivity to both O_{2} and CO_{2} (the respective controller gains for O_{2} and CO_{2}) by manipulating inspired O_{2} and CO_{2} under baseline conditions and again under the influence of domperidone (19). Additionally, the risk for ventilatory instability under both conditions was characterized using an established model of hyperventilation-induced periodic breathing (13, 56) (see Fig. 2). All lambs were euthanized at the end of the experiment using an overdose of anesthetic (Lethabarb; sodium pentobarbitone, 150 mg/kg; Virbac, Sydney, Australia).

### Preprocessing

All data were low-pass filtered at a cut-off frequency of 40 Hz. Respiratory flow recordings were corrected for drift and were integrated breath-by-breath to obtain a continuous tidal volume signal. The start and end of inspiration were determined by detecting local maxima and minima of the tidal volume signal and were confirmed visually. Minute ventilation (V̇e) was calculated for each breath as V_{T}/T_{tot}, where V_{T} is the tidal volume and T_{tot} is the duration of that breath. To minimize the influence of noise, Pet_{CO2} and Pet_{O2} values were calculated by taking the mean of the last quarter of the relevant waveform during expiration. The breath-to-breath time series data were further high-pass filtered to remove the steady-state baseline and any oscillation in the signals slower than 50 breaths/cycle (using a 7^{th} order Butterworth digital filter with cutoff frequency of 0.02 cycles/breath). The resulting breath-to-breath V̇e, Pet_{CO2}, and Pet_{O2} time series, representing deviations from a steady-state baseline during spontaneous breathing, were used for subsequent analysis (see Fig. A1).

### Trivariate Autoregressive Modeling

Our model represents the three key variables (V̇e, Pet_{CO2}, and Pet_{O2}) as linear functions of their history (hence the label “autoregression”) and random fluctuations. The deterministic component represents the chemical control, including chemoreflexes and gas exchange, while the stochastic component describes how external random fluctuations propagate through the system. The parameters that characterize this system can be utilized to describe the pairwise interactions between the model components (42, 27, 30, 2) and to derive the frequency-domain and stability characteristics of the underlying system. A trivariate model with maximal lag (or memory) of *p* breaths that describes the interactions among the three modeled ventilatory measurements V̇e, Pet_{CO2}, and Pet_{O2} can be represented by the following matrix equation:
*x*(*n*) comprises V̇e, Pet_{CO2}, and Pet_{O2} at breath *n*; *x*(*n*) lists the values of these three variables at the *kth* previous breath; the matrices *a*(*k*) for *k* = 1,…, *p* represent the static gains that relate *x*(*n-k*) to *x*(*n*); *w*(*n*) represents the variations in V̇e, Pet_{CO2}, and Pet_{O2} that are not explained by the chemical control system properties and are therefore considered to be the result of external stochastic disturbances to the system, i.e., noise. *Equation 1* states that the values of V̇e, Pet_{CO2}, and Pet_{O2} at any breath *n* are *a* linear functions of their *p* previous values, plus an independent random term *w*. Each variable is therefore broken down into four components, incorporating the additive influence of its own history, the histories of the other two variables, and noise *w*. Furthermore, the model is “breath invariant”: the matrix coefficients *a*(*k*) do not depend on the breath number *n*. We assume that the individual elements of *w* are uncorrelated with each other at any breath, and across different breaths; we also assume they have mean value zero and constant variances σ_{V̇e}^{2}, σ_{PCO2}^{2}, and σ_{Po2}^{2}, respectively.

In this study, we fixed the analysis window at 125 breaths and utilized a least-square-error based identification technique (34, 35) to identify the model parameters (see appendix d for a discussion on model order and window size selection).

#### Calculation of controller, plant, and loop gain.

The transfer path functions of the control system (controller and plant gains for CO_{2} and O_{2}) describe the direct effect of one variable in the system on another (27). These functions are derived from the model coefficients *a _{ij}* (see appendix a). Each such function is frequency dependent (in units of cycles per breath), and is characterized by a magnitude and phase. We denote by T

_{j→i}(

*f*) the frequency-dependent characteristic of the direct pathway connecting the

*jth*signal to the

*ith*signal (see Fig. 1

*B*). For instance, | T

_{Pco2→V̇e}(

*f*)| (units of l·min

^{−1}·mmHg

^{−1}) describes the magnitude of the change in the component of V̇e of frequency

*f*per change in the corresponding frequency component of Pa

_{CO2}(controller gain for CO

_{2}at frequency

*f*). Similarly, | T

_{V̇e→Pco2}(

*f*)| characterizes the magnitude of the change in the component of Pa

_{CO2}at frequency

*f*per change in the corresponding frequency component of V̇e (plant gain for CO

_{2}at frequency

*f*). The CO

_{2}and O

_{2}loop gains [denoted by LG

_{CO2}(

*f*) and LG

_{O2}(

*f*), respectively, or with frequency dependence understood] are defined as the product of their respective plant and controller gains (i.e., LG

_{CO2}= T

_{V̇e→Pco2}T

_{Pco2→V̇e}= and LG

_{O2}= T

_{V̇e→Po2}T

_{Po2→V̇e}). Finally, the overall loop gain (LG) turns out to be the sum of CO

_{2}and O

_{2}loop gains: LG = LG

_{CO2}+ LG

_{O2}(see appendix b for the demonstration).

#### Impact of external disturbances on ventilatory variability: role of loop gain.

Our model-based characterization of the ventilatory control system allows for the quantification of the degree to which fluctuations in ventilation result from external fluctuations in Pco_{2}, Po_{2}, or V̇e, namely the “noise” component of these variables that is not explained by the chemical control system (i.e., the *w* term in *Eq. 1*). By characterizing the size of the fluctuations in ventilation with respect to the external noise, we can interpret how external disturbances in Pco_{2}, Po_{2}, or V̇e propagate through the feedback loops, resulting in the emergence of oscillatory patterns in ventilation. This is achieved by using the fluctuation transfer function (appendix b), denoted by H_{Wj→1}(*f*), whereby random disturbances in signal *j* contribute to the power in signal *i*. For example, |H_{WPCO2→V̇E}|^{2} [units of (l·min^{−1}·mmHg^{−1})^{2}] quantifies the power (i.e., time average of the squared magnitude) of fluctuations in V̇e at a given frequency contributed per unit power in the external random disturbance in *w*_{Pco2} at the same frequency.

The dependence of H_{Wj→1} on the overall loop gain and the transfer path functions within the various system loops can be described explicitly (see appendix b). For instance, we show that the ventilatory response to noise in CO_{2} is given by H_{WPCO2→V̇E} = T_{Pco2→Po2}T_{Pco2→V̇e}/(1 − LG). As expected, a smaller noise-buffering effect for CO_{2} (greater T_{Pco2→Po2}) or a greater controller gain for CO_{2} (T_{Pco2→V̇e}) contributes to a greater impact of noise in CO_{2} on ventilation. Notably, the finding that H_{WPCO2→V̇E} is dependent on the overall loop gain (LG = LG_{CO2} + LG_{O2}), and specifically on how close the loop gain is to 1, demonstrates that a given extrinsic disturbance in CO_{2} results in greater ventilatory variability in a higher gain system, even if overall loop gain is dominated by the O_{2} control loop. Thus, if the overall loop gain increases with domperidone, one expects from the analysis that H_{WV̇E→V̇E}, H_{WPCO2→V̇E} and H_{WPO2→V̇E} will increase in magnitude and will exhibit a maximal amplitude at the frequency at which periodic breathing is expected to occur (i.e., the natural frequency of the system). To obtain the natural frequency of the system, we therefore calculated the frequency at the peak amplitude of the H matrix within the medium frequency range (see appendix b for details) and compared it to the measured cycle duration of periodic breathing in individual animals, both before and after domperidone.

#### Signal power as a measure of variability.

The power spectrum of a signal *i* (denoted P* _{i,i}*) describes the frequency dependent variability of the signal, where a prominent peak in the power spectrum indicates the existence of a dominant oscillatory component. Similarly, the cross-power between two signals

*i*and

*j*(denoted P

*) is a frequency-dependent function that quantifies their cross-correlation properties. From our model, analytical forms for P*

_{i,j}*and P*

_{i,i}*are derived (and is given in*

_{i,j}*Eqs. C1*and

*C2*of appendix c). An empirical estimate of the power spectrum of a signal, for comparison with the analytical expression, can be obtained by using the periodogram averaging technique, computed via the fast Fourier transform of the signal (16). A close match between the model-based and the periodogram-based spectra indicates that the model provides an adequate description of the periodic components of the data.

### Data Analysis and Statistics

To facilitate comparison between control and domperidone conditions, we considered three distinct frequency bands: low frequency (LF) oscillations with periods 16–50 breaths/cycle, medium frequency (MF) oscillations of 5–15 breaths/cycle, and high frequency (HF) oscillations from 2–4 breaths/cycle. Within each frequency band, the magnitude of the derived quantities (controller gains, plant gains, loop gains, P* _{i,j}*, etc.) were averaged, and the resulting indexes compared between the control and domperidone studies, using the Wilcoxon signed rank test with significance level

*P*< 0.05. Data are presented as medians and interquartile ranges. Since the MF band spans the range of cycle durations of periodic breathing observed experimentally (22), we focus on the values of controller and plant gain, and the overall loop gain, measured in this band. For the purpose of predicting the predisposition to ventilatory instability, we consider loop gain magnitude in the MF band to be an indicator of nearness to instability, with a value >1 suggesting an unstable system, and a value <1 suggesting a stable system. This indicator attempts to capture in a simple way, suited to our purpose, the essence of the Nyquist criterion (28).

## RESULTS

### Respiratory Variables and Experimentally Derived System Properties

The average baseline respiratory variables derived from the breath-by-breath time-series are shown in Table 1. The administration of domperidone resulted in a slight increase in V̇e, mainly due to a decrease in T_{i} (*P* < 0.05) and a near significant decrease in T_{tot} (*P* = 0.09). These changes were accompanied by a decrease in Pet_{CO2} (*P* < 0.05) and an increase in Pet_{O2} (*P* < 0.05).

### Trivariate Analysis

#### Controller, plant, and loop gain.

Transfer path analysis revealed a significant increase in the magnitude of O_{2} chemosensitivity (T_{Po2→V̇e}) across the LF and MF band following administration of domperidone (see Fig. 3). Figure 4 shows the comparison of the control and domperidone studies for the average gain magnitudes within the MF range for each subject. We found a significant increase in the median value of the average MF band gains in T_{Po2→V̇e}. Although the median value of the average MF band gain in T_{Pco2→V̇e} was increased with domperidone, this failed to reach statistical significance. The complete analysis results for the CO_{2}- and O_{2}-specific loop gains, and the overall LG, are presented in Table 2; both the LG_{O2} magnitude and the overall LG magnitude increased significantly in the LF and MF regions. Notably, the LG_{O2} magnitude under both control and domperidone conditions was greater than LG_{CO2} magnitude, and the overall LG was comparable to the LG_{O2}. Moreover, a twofold increase in the overall LG magnitude with domperidone was observed in conjunction with the greater propensity towards periodic breathing (1/15 vs. 13/15 animals exhibited periodic breathing; see Table 2). Figure 5*A* shows that model-based CO_{2} and O_{2} controller sensitivities in the MF band compared favorably with the experimentally measured ventilatory sensitivity to CO_{2} and O_{2}, respectively (*P*< 0.05).

#### Impact of external disturbances on ventilatory variability: role of loop gain.

Figure 6, *A–C*, shows the average fluctuation transfer function (H) spectra that quantify the influence of external disturbances on V̇e. Note the increase in the amplitude of the spectra and emergence of a sharp peak in the MF band after administration of domperidone. Figure 5*B* is a plot of the cycle duration of experimentally induced periodic breathing vs. the cycle duration of these dominant MF band peaks, showing good agreement between the predicted and experimentally measured cycle durations [predicted: 7.9 (7.0, 9.2), measured: 7.5 (7.2, 9.1) breaths]. Moreover, Fig. 5*C* demonstrates a significant increase in the number of experimentally observed periodic breathing cycles with an increase in the overall LG (*r*^{2} = 0.50; *P* < 0.05).

#### Power and cross-power.

Figure 7 shows the model-based power spectra of the individual signals overlaid on their periodogram-based power spectra. Before the administration of domperidone, the V̇e spectrum (Fig. 7*A*) shows a peak around the cycle duration of 5.5 breaths/cycle, although the signal power is relatively low. After the administration of domperidone (Figs. 7, *D–F*), all three signals exhibit marked periodicity around the cycle duration of eight breaths/cycle. Moreover, under domperidone (Fig. 7*D*) the model-based and the emperical power spectra of V̇e closely match.

## DISCUSSION

The major finding of our study is that there is sufficient information in the natural fluctuations in V̇e, Pco_{2}, and Po_{2} during spontaneous breathing to *1*) estimate, through the use of a physiologically based dynamic model, the magnitude of the dynamic hypercapnic and hypoxic ventilatory responses, and *2*) predict the propensity to periodic breathing (associated with high loop gain), and the corresponding cycle duration of periodic breathing, should it occur. The power of our analytical method is that it permits quantitative determination of controller and plant gains for both CO_{2} and O_{2}, and thus overall loop gain of the respiratory control system, without the need for externally applied perturbations of the respiratory system. Our results are consistent with our earlier finding in newborn lambs: that administration of domperidone increases the loop gain, predominantly by an increase in the O_{2} controller gain, making the respiratory control system more prone to instability (19). For the first time to our knowledge, we demonstrate quantitatively that the controller gains for O_{2} and CO_{2} in the lamb have comparable magnitudes, but a higher plant gain for O_{2} results in the overall loop gain being dominated by the O_{2} feedback path rather than that for CO_{2}. This finding is consistent with experimental observation in the same animals: that stabilizing the O_{2} feedback path during periodic breathing, but not the CO_{2} feedback path, stabilizes ventilation (20). By virtue of being based entirely on noninvasive measures of spontaneous breathing, our approach represents a new methodology to *1*) assess the risk of ventilatory instability in vulnerable populations, *2*) identify the primary causes of instability in individual patients and design appropriate treatment, and *3*) follow patient progress in response to therapeutic interventions.

### Comparison with Other Approaches: Estimated vs. Measured Hypoxic/Hypercapnic Responses

In our previous work, by presenting different concentrations of inhaled CO_{2} and O_{2}, we determined the corresponding controller gains under control conditions and following domperidone administration. We observed that CO_{2} response to three breaths of CO_{2} increased with domperidone (19). However, in the current study, we found no statistically significant increase in CO_{2} controller gain. This disparity is likely to be due to the large perturbation used previously (14% inspired CO_{2}) whereas the current method relies on spontaneous oscillations in CO_{2} to calculate CO_{2} controller gain. Nevertheless, the MF band CO_{2} and O_{2} transfer path gains were significantly correlated with the experimental (traditional) estimates of controller gain (see Fig. 5*A*) indicating our trivariate analysis sufficiently captured the measured changes. Interestingly, the experimental eupneic controller gain for O_{2} was not significantly different between control and domperidone (data not shown); despite the difference previously described at hypoxic levels (19). Importantly, our novel method was sufficiently sensitive to detect a significant increase in MF band O_{2}-controller gain during spontaneous breathing with the administration of domperidone (see Fig. 4). As such, our method accurately captures changes in sensitivities of the respiratory controller to O_{2} from spontaneous breathing records, obviating the need to perform labor-intensive and technically challenging experiments.

### Comparison with Other Approaches: Independent Laboratories

The link between respiratory variability during spontaneous breathing and the chemoreflex feedback system in adult human subjects has been previously investigated (29, 41). Of note is the work of Van den Aardweg and Karemaker (16), who observed that mean inspiratory flow time series in healthy human subjects contains oscillatory components with a cycle duration of ∼10 breaths/cycle. Based on a simplified model of the CO_{2} chemoreflex feedback loop, they concluded that the correlations in breath-to-breath values of the mean inspiratory flow may reflect the response of the chemoreflex feedback system to (uncorrelated) noisy disturbances in the partial pressure of arterial carbon dioxide (Pa_{CO2}). More recently, approaches using spontaneous changes in ventilation and CO_{2} have been developed, with one study (40) demonstrating that opioid administration (remifentanil) substantially increases plant gain but reduces controller gain. These studies assume that fluctuations in O_{2} have negligible impact on ventilation, and have not shown a direct link between spontaneous breathing variability and the occurrence or cycle duration of periodic breathing. The current work goes beyond these previous studies in two important ways. First, the use of a trivariate analysis incorporates the influence of Po_{2} fluctuations as well as the Pco_{2} on changes in V̇e. As such, our method has general appliability since it does not make any prior assumptions regarding the dominance of one loop vs. the other. Second, our study goes beyond a phenomenological description of the ventilatory variability and examines the mechanism responsible for the amplitude and frequency characteristics of the observed ventilatory periodicity.

Our study confirms the view (41, 26) that fluctuations in breath-to-breath values of V̇e provide insight into ventilatory instability. Thus increasing the gain of the system with domperidone causes an increase in magnitude of the ventilatory variability per unit of extraneous noise in Po_{2}, Pco_{2}, and ventilation (as reflected in the fluctuation transfer functions for V̇e) within the medium frequency band. We also show explicitly that the closer LG is to 1, the greater the ventilatory variability that will result from any given amplitude of external perturbation (see *Eq. B2–B5* in appendix b). Specifically, small random perturbation in ventilation (i.e., noise) of amplitude A, will become approximately augmented by 1/(1 − LG) to exhibit fluctuations of amplitude A/(1 − LG); for example, if LG = 0.9 then an underlying variability with a mean amplitude of 0.5 l/min will be amplified 10-fold to exhibit oscillations of mean amplitude of 5 l/min. Hence, while periodic central apneas will consistently occur in a control system with LG > 1, periodic-like ventilatory variability is also expected to occur in a system with moderately high gain (loop gain <1, but close to 1) in the presence of sufficient underlying ventilatory variability. Importantly, our data demonstrate (Fig. 6) that frequency content of this ventilatory variability reflects the natural frequency of the system, as confirmed by our finding that cycle duration predicted from spontaneous breathing corresponds closely with the cycle duration of experimentally induced periodic breathing. In addition, the experimentally observed number of cycles of periodic breathing increases as the estimated overall LG approaches 1. Although a large LG is only one of the several factors contributing to an increase in the number of cycles, the number of cycles has been previously used as a surrogate measure of an increased loop gain (57, 19); because the number of cycles is expected to increase as the “damping” of oscillations decreases, i.e., as LG approaches 1. Our results indicate that >50% of the variation in the number of cycles of periodic breathing can be explained by the overall LG estimates during spontaneous breathing. Taken together, our evidence leads us to conclude that the periodic components of spontaneous ventilatory variability and periodic breathing are emergent properties of a high gain respiratory control system amplifying the external disturbances entering the chemoreflex feedback loop.

### Methodological Considerations

The model presented here is based on a number of simplifying assumptions. First, it is assumed that Pet_{CO2} and Pet_{O2} are accurate proxies for their corresponding arterial values (14). However, arterial and end-tidal Po_{2} can differ substantially when gas exchange is not ideal. On the basis that Pet_{CO2} and Pet_{O2} are strongly correlated with their respective arterial values, our measures of loop gain are not expected to be affected unduly by this assumption. However, a constant ratio of arterial and end-tidal Po_{2}, as observed dynamically in the lamb (55), is expected to yield an underestimated controller gain and an overestimated plant gain compared with the true arterial values.

Our second assumption was that the current values of V̇e are linearly dependent on their own history and the previous values of Pet_{CO2} and Pet_{O2} and also that Pet_{CO2} and Pet_{O2} are only dependent on their own respective histories along with the previous values of V̇e. Thus we assumed a negligible synergistic interaction between Pco_{2} and Po_{2} at the controller or the plant, which we considered valid for the case of small purturbations in Po_{2} and Pco_{2} around resting levels. However, such an effect may be important in the case of large disturbances in the human adult, when a large hypoxic dip would be expected to increase the controller gain for CO_{2}. For our application, such Po_{2}-Pco_{2} interaction was deemed negligible, since the controller responses to Pco_{2} and Po_{2} are additive in the lamb (20). Third, although the calculated values for CO_{2} controller gain in the mid-frequency range could be a combination of both peripheral and central chemosensitivity to CO_{2}, it is likely that the relatively long brain tissue wash-in/wash-out time constants for CO_{2} (39, 4, 43, 50) limit the central chemoreceptors from contributing significantly to rapid or transient changes resulting from breath-to-breath variability.

From a model-fitting prospective, nonstationarities may arise as a result of nonlinearly interacting processes (48, 46) at time scales corresponding to fluctuations in cardiac output, respiratory rate, sleep-state or arousal-related changes in respiratory mechanics and the chemical control system, behavioral factors, and those induced by interventions. In the current study, we addressed this issue by using a relatively short analysis window to minimize the impact of nonstationarities. Moreover, a strong collinearity (or high correlation) between modeled variables may be a confounding factor in determining the true contribution of CO_{2} and O_{2} to ventilation. In our data, the maximum correlation coefficient between CO_{2} and O_{2} (over the range of −4 to 4 breath lags) was −0.48 (−0.38 −0.64) under baseline conditions and −0.61 (−0.41 −0.78) after domperidone administration. This intermediate degree of correlation between breath-by-breath fluctuations in CO_{2} and O_{2} is expected from respiratory system physiology, given the different buffering capacities of the body for CO_{2} and O_{2}.

Reliability of the estimated model parameters relies on the existence of adequate perturbations in the frequency bands of interest so that every possible pathway is sufficiently excited. Notably, the spontaneous variability was sufficient in our study to achieve a close correspondence between the spectral characteristics of our model output and the experimental data (Fig. 7), although the correspondence appeared to improve as ventilatory variability was increased following domperidone administration. Under conditions of insufficient natural fluctuations, it may be beneficial to subject the system to external perturbations; in principle, our autoregressive method is also applicable to the pseudorandom binary sequence CO_{2} stimulation technique (22).

We explicitly assume that the pathways in Fig. 1*B* are the dominant ones but are mindful that the respiratory control system could be affected by fluctuations in heart rate, blood pressure, cerebral blood flow, and sympathetic activity. For instance, Henry et al. (24) reported a possible coupling between the chemoreflex and baroreflex mechanism, in which case extra pathways representing the interactions among heart rate, blood pressure, and respiration need to be incorporated. However, this particular effect is likely to be minimal in the current work, as the lamb data we used came from a study in which mean arterial blood pressure was stable (19). Regarding other confounding factors, we simply note that the concordance between our estimates and the experimental measurements (in both the cycle duration of periodic breathing and in the independent assessments of hypoxic and hypercapnic respiratory responses) strongly suggests that the model we propose captures the essentials of the linear feedback mechanisms of the respiratory control system.

### Utility of Our Analytic Approach

Our noninvasive method for quantifying the dynamic hypercapnic and hypoxic ventilatory responses and overall loop gain potentially has widespread clinical applicability. Increased loop gain as a phenotypic trait in obstructive sleep apnea is predictive of the response of patients to treatment with supplemental oxygen (54). An elevated hypoxic ventilatory response at high altitude is associated with the presence of unstable breathing (31). Moreover, heightened dynamic hypercapnic responses delineate heart failure patients with periodic breathing in the form of Cheyne-Stokes respiration from those patients that do not exhibit ventilatory instability (51, 21, 52). Importantly, in heart failure patients, elevated hypoxic and/or hypercapnic sensitivity (45) as well as the presence of Cheyne-Stokes respiration itself (32) are important predictors of mortality. Thus model-based characterization of the respiratory control system using spontaneous breathing provides a simple noninvasive diagnostic tool, allowing potential therapies to be directed at the dominant feedback loop (i.e., CO_{2} vs. O_{2} loop) responsible for the instability.

### Conclusions

This study shows that it is possible to characterize and quantify the dominant feedback loops within the respiratory control system from measurements of fluctuations in ventilation and respiratory gas tensions during spontaneous ventilation. Consequences of this approach include the ability to assess risk of ventilatory instabilities, cycle duration of periodic breathing, and hypoxic and hypercapnic ventilatory responses. The major advantages of our method are that it requires minimal subject intervention and it makes no demand on subjects to comply with complex breathing protocols. Through overall loop gain, our method promises new ways to assess the potential for respiratory instabilities, especially important in cases of congestive heart failure or obstructive sleep apnea. Importantly, identification of the mechanisms responsible for an increase in loop gain in such individuals would enable clinicians to target particular therapies on an individualized basis and to monitor patient response to such therapies.

## GRANTS

This work was supported by the American Heart Association Grant 0840159N and National Heart, Lung, and Blood Institute Grants R01-HL-73146, HL-085188-01A2, HL-090897-01A2, and K24-HL-093218-01A1 and Training Grant T32-HL07901. B. A. Edwards is the recipient of the Thoracic Society of Australia and New Zealand/Allen and Hanburys Respiratory Research Fellowship.

## DISCLOSURES

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

## ACKNOWLEDGMENTS

The content of this document is solely the responsibility of the authors and does not necessarily represent the official views of the American Heart Association or National Institutes of Health.

## APPENDIX

In this section, we describe the method of transfer path analysis, which yields frequency-domain expressions (transfer path functions) for the gains of the chemoreflex controllers, plants, and feedback loop for CO_{2} and O_{2}. We also explicitly demonstrate the dependence of fluctuation transfer functions on the overall loop gain and the transfer path functions within the various system loops. Furthermore, we derive analytical expressions for the power spectra of the modeled variables. Finally, we conclude with a discussion of model order selection and data segmentation.

## APPENDIX A: TRANSFER PATH ANALYSIS

*Equation 1* is a discrete convolution; in the frequency domain, this becomes simply multiplicative. Define the Fourier transform of a time signal by Q(*f*) = Σ* _{k}q*(

*k*)exp(−2π

*k*indexes the time domain, and

*f*indexes the frequency domain. Our convention is to denote time domain variables in lowercase, and transformed variables in uppercase. The Fourier transform of

*Eq. 1*is thus

*f*will be omitted for notational simplicity. Note that the summation over time in defining A(

*f*) = Σ

_{k = 1}

^{p}=

*a*(

*k*)exp(−2π

*fk*) includes only the first

*p*points, where

*p*is the maximal lag or memory in the model of

*Eq. 1*. Indexes (1, 2, 3) denote (V̇e, Pco

_{2}, and Po

_{2}).

*Eq. A1*can be written as

_{Pco2,Po2}and A

_{Po2,Pco2}are set to zero (see

*Choice of Model Order*). Figure 1

*B*is a graphical representation of

*Eq. A2*. The transfer path function from the

*jth*signal to the

*ith*signal is denoted by T

*and is given by*

_{j→i}*Eq. A3*is a consequence of the assumption that the current values of the time series depend on their previous values and as such it captures the notion of a memory in the system (Fig. A1).

## APPENDIX B: FLUCTUATION TRANSFER FUNCTIONS

Solving for X in *Eq. A1* yields:
^{−1} is a 3 × 3 matrix of functions relating the fluctuations in the measured variables to the sources of fluctuation defined in *Eq. 1*. Taking A_{Pco2,Po2} = A_{Po2,Pco2} = 0 implies T_{Pco2,Po2} = T_{Po2,Pco2} = 0, in which case the components of H pertinent to the influence of external disturbances on V̇e are given by
(B2)
(B3)
(B4) where,
(B5) *Equations B2–B5* show clearly that increased loop gain amplifies the influence of the noise terms *w*_{V̇e}, *w*_{Po2}, and *w*_{Po2}. [Note that Mason's rule (37) generalizes *Eqs. B2–B4* to an arbitrary number of variables in a feedback system.] From *Eq. B5*, we draw the important conclusion that the overall loop gain, for the specific model structure shown in Fig. 1*B*, is simply the sum of the individual CO_{2} and O_{2} loop gains.

The maximum possible amplification achievable by any combination of external disturbances entering the chemoreflex feedback loop is known as the matrix 2-norm of the H matrix [denoted by ∣∣H(*f*)∣∣_{2}]. It is frequency dependent, which implies that the system may selectively amplify certain frequencies in the input disturbances. The location of these frequencies is determined by the poles of the H matrix, which determine the natural frequencies of the system, that is, frequencies where even small periodic driving forces can produce large amplitude oscillations. For the system under study, the poles are given simply by the roots of the denominator 1 − LG (*Eq. B5*). In general, the system poles are complex numbers of the form *re*^{iθ}, with magnitude *r* and angle θ = 2π*f*. We therefore calculated the cycle duration (i.e., inverse of the frequency) of the largest magnitude pole in the MF band for comparison with the cycle duration of posthyperventilation periodic breathing in individual animals, before and after domperidone.

## APPENDIX C: POWER

The power spectrum of the vector time-series *x*(*k*) is defined in the frequency domain by (47):
^{T}] is the spectral density matrix associated with the fluctuations, and E is the expectation operator. The components of P are estimated by simple averages, being unbiased estimators of the expectations. Under the assumption that the fluctuations are independent (i.e., uncorrelated) across the variables and also from breath-to-breath on each variable, the matrix Σ has zero off-diagonal elements and constant (i.e., not f-dependent) diagonal elements (σ_{1}^{2}, σ_{2}^{2}, σ_{3}^{2}) = (σ_{V̇e}^{2}, σ_{PETCO2}^{2},σ_{PETO2}^{2}). The diagonal elements of P give the power spectra of the corresponding variables. For instance, the power spectrum of the fluctuations in ventilation is given by
(C2) This is a weighted combination of the powers due to each individual source of fluctuation, with weights given by the squared magnitudes of the fluctuation transfer functions. Similarly, the off-diagonal elements of the power spectral matrix P give the cross-power term between each pair of variables in the ventilation, CO_{2}, and O_{2} triple.

## APPENDIX D: MODEL ORDER SELECTION AND DATA SEGMENTATION

#### Choice of Model Order

In autoregressive modeling, the optimal model order is typically estimated from the data through evaluation of the Akaike's information criterion (AIC) (1). The AIC selects a model that optimizes the trade-off between a parsimonous model (lowest model order) and one that provides the best predictive power (minimizing prediction error). The use of AIC is based on the assumption that no prior information concerning the model structure is available. Note that the proposed model is physiologically structured, that is, it is based on our current understanding of the physiology of respiratory control system. Dependence of the key variables on their own history captures the notion of a memory in the system. In the case of ventilation this memory is neural (3); while in the case of the blood gases, this memory may reflect the buffering effect of the functional residual capacity and lung tissue (16).

In this work, the orders of the autoregressive terms associated with the dynamics of the plant (that is, the terms *a*_{Pco2,V̇e}, *a*_{Pco2,Po2}, *a*_{Po2,V̇e}, and *a*_{Po2,Po2} in *Eq. 1*) were set equal to 1 in accordance with earlier studies (33, 22, 16). Moreover, the terms *a*_{Pco2,Po2}(*k*) and *a*_{Po2,Pco2}(*k*) were set to zero for all *k*, since we assumed a negligible interdependence of Pet_{CO2} and Pet_{O2}. Such an assumption vastly simplifies the model analysis. This is justified because our model deals with small (1^{st} order) fluctuations about nominal resting values; the interaction term is second order and can therefore have only a minor effect on the measured ventilatory control variables. Although these choices of model parameters for the plant resulted in 9 out of 15 unwhite residuals (that is, model residuals with statistically significant correlation with respect to the variable being estimated for lags 1 or 2), the identified plant dynamics is physiologically sound. Increasing the model order by up to two additional lags eliminated the issue of whiteness of residuals but made the results harder to interpret, so those results are not discussed in this work.

To achieve a minimal model capturing the chemoreflex dynamics, the model order associated with the autoregressive terms *a*_{V̇e,V̇e}, *a*_{V̇e,Pco2}, and *a*_{V̇e,Po2} in *Eq. 1* was increased from 2 to 10 breaths, and the lowest AIC order that resulted in white residuals was selected (34). The model orders varied between two and five time lags across subjects. From a physiological point of view, it takes at least two to three breaths for a change in blood gases at the alveolar level to reach the peripheral chemoreceptors. Therefore, to have a model that is physiologically sound and to facilitate intersubject comparison and group averaging, we set the model order equal to 4 for all the subjects and across different protocols. The model order utilized here is well within the bounds of the orders reported in univariate analysis of the chemoreflexes in Khoo and Marmarelis (29), also by Ghazanshahi and Khoo (22), who used autoregressive and moving average of order 1 and one to three pure lags representing the convective delay for pulmonary capillary blood to reach the carotid chemoreceptors (thus considering signal history up to 4 lags), and the order reported by Clement and Robbins (15), who estimated a median lag of two breaths to carotid body chemoreceptors. Figure *D1* shows the averages and SDs of the model parameter values between each pair of variables.

#### Window Size Selection

The influence of window size on model fitting was studies by varying the window size (75, 100, 125, and 150 breaths), and a window size of 125 breaths was selected based on the smallest AIC and loss function. The loss function is defined as the determinant of the residual covariance matrix Σ. The AIC is related to the loss function, model order (*p*), and window size (*N*) via the relationship AIC = log[det(Σ)] + 2*p/N*. In general, it is recommended to have several times as many data points as the number of model parameters (here we have 16 parameters to fit, so 125/16∼8). On the other hand, the window size has to be kept small enough to assume local stationarity of the joint time series (18). The choice of 125 breaths as the window size resulted in 1–5 windows of spontaneous breathing data across different subjects. For each subject the window with the largest power in the medium frequency band of V̇e (assessed via periodogram-based power spectral estimation) was selected for model fitting purposes. Our rationale behind this choice of data segment selection was to enhance model identification by choosing a segment of data with sufficient variability (see *Methodological Considerations*) and to characterize the system at its most variable (or unstable) observed state during spontaneous breathing.

- Copyright © 2011 the American Physiological Society