Abstract
Liang, PeiJi, Jaideep J. Pandit, and Peter A. Robbins. Statistical properties of breathtobreath variations in ventilation at constant Pet _{CO2} and Pet _{O2} in humans. J. Appl. Physiol. 81(5): 2274–2286, 1996.—The purpose of this study was to provide a statistical description of the breathtobreath variations in ventilation during steady breathing in both rest and during light exercise, with the endtidal gases controlled by using an endtidal forcing system. Sixty data sets were studied, only one of which was white (i.e., did not show autocorrelation). Three simple autoregressive moving average (ARMA) models, i.e., AR_{1}, AR_{2}, and AR_{1}MA_{1}, and one simple statespace model were fitted to the data and resulted in white residuals in 15, 31, 46, and 48 out of 60 occasions, respectively. Evolutionary spectral analysis revealed that only 13 data sets had a constant power spectrum, although 50 were uniformly modulated. An autoregressive estimate of variance could be used to “demodulate” the data in most cases, but the results were not significantly affected by fitting the model to the demodulated data. The results indicate that 1) both simple ARMA models and a simple statespace model can describe the autocorrelation present; 2) variations in spectral power were present in the data that cannot be described by these models; and3) these variations were often due to a uniform modulation and did not significantly affect the coefficients for the models. For these kinds of data, a heteroscedastic form of statespace model provides an attractive theoretical structure for the noise processes.
 timeseries models
 statespace model
 constant power spectrum
 uniform modulation
 heteroscedastic form
 endtidal partial pressure of oxygen
 endtidal partial pressure of carbon dioxide
in the early 1960s, Priban (17) demonstrated statistically that, in resting human subjects, successive breaths were correlated. Since this study, a number of investigators have used a variety of models to describe the breathtobreath variation in breathing (1, 3, 8, 9, 11, 13, 15). These studies cover a wide range of different conditions, including studies of breathing in animals as well as in humans, studies during sleep and anesthesia as well as during consciousness, and studies in the presence of changing stimulation as well as in the presence of steady stimulation. This makes the results from these studies difficult to compare. Furthermore, many of these studies are somewhat limited in scope. In particular 1) the investigators often have not examined some of the basic statistical properties of the sequences before applying their models to their data;2) the investigators often have studied only one model rather than undertaking a comparison across a range of possible models; and3) the investigators often have not checked the adequacy of the model for describing the correlation once the model had been fitted.
The present study attempts to meet some of the criticisms that can be leveled against these earlier studies. It focuses on steady breathing in conscious humans during both rest and mild exercise under euoxic conditions. The CO_{2} was slightly elevated compared with air breathing so that the technique of endtidal forcing could be used to hold endtidal Pco _{2} and Po _{2} (Pet _{CO2}and Pet _{O2}, respectively) as constant as possible. The reason for this was that we wished to concentrate on the statistical properties of the ventilatory controller in the openloop setting rather than on the properties of the entire closedloop system, which would include the dynamics of the gas stores within the body. Thus the purpose of the study was to provide a statistical description of the breathbybreath ventilatory sequences associated with the respiratory controller in an openloop state in conscious humans.
As well as the description’s own intrinsic interest, it has a particular application for studies of the respiratory control system using the endtidal forcing technique to generate dynamic stimuli in the endtidal gases. These studies often involve fitting one or more models to the ventilatory sequences obtained, and to do this the model should contain elements to describe both the deterministic and the stochastic parts of the process. Therefore, a particular aim of the present study was to develop a suitable structure for the stochastic part of such models.
METHODS
Data
The experimental data came from a previously published study from our laboratory (16). Experiments were performed on five subjects, with two different protocols each lasting 43 min. In the first protocol ( protocol A ), the subject was at rest, and the Pet _{CO2} was held constant at 2–5 Torr above the subject’s natural resting value. In the second protocol ( protocol B ), the subject undertook exercise at 70 W, and the Pet _{CO2} was held at 2–5 Torr above the subject’s natural value during exercise. During both protocols, Pet _{O2} was held at 100 Torr. Each protocol was repeated six times for each subject.
Ventilatory data during steadystate breathing were used for model fitting and statistical analysis. Data collected during the first 50 breaths for protocol A and the first 300 breaths forprotocol B were discarded to remove any initial transients in the response.
ModelIndependent Statistical Analysis
The mean value for ventilation together with any linear trend term were removed from the data before any statistical analysis and model fitting was undertaken. The values for the mean ventilations and the trend terms are shown in Table 1.
Specific periodicities.
Spectral analysis was used to investigate whether there were any specific periodicities in the data. For each data set, a smoothed periodogram was calculated by using a modified Daniell window. The modified Daniell window was obtained by combining three Daniell windows of lengths of 3, 5, and 7. The ratio of the estimated spectral densityĥ(ω) over its real value h(ω) approximately follows a χ^{2} distribution such that {ĥ(ω)/h(ω)} ∼
Truncation of each data set from the same subject and the same protocol to the same length ensures that the estimated spectral densityĥ(ω) over its real value h(ω) will follow a χ^{2} distribution with the same degrees of freedomv for each of the data sets (since this depends entirely on the structure of the window and the length of the data set). Thus the average value for the six estimated spectral densities (for the same subject and the same protocol) over its real value will also follow a χ^{2} distribution, with the degrees of freedom being simply six times that of the original data sets (6v). This enables the 95% confidence intervals for the mean spectral density to be calculated in a similar way to that used for the spectral densities from the individual data sets.
Evolutionary spectral analysis.
Evolutionary spectral analysis is a technique for assessing whether there are changes occurring in the power spectrum of a series over time (18, 20). The principle is to apply a double window in both time and frequency domains to calculate values for evolutionary spectral density corresponding to the selected times and frequencies. First, the estimate (U_{t}) for the spectrum for a central frequency of angular velocity ω_{0} is calculated by using
An example of this form of analysis is shown in Table 2 for subject 797 at rest. The values for logarithmic spectral densities for selected time and frequency points are listed in Table2, whereas the results of variance analysis are listed in Table 3. The first step is to test the “interaction + residual” term [χ^{2} = (sum of squares)/ς^{2}; in this study, ς^{2} = 4h/3T ′ = (4 × 7)/(3 × 100) = 7/75 (see Ref. 20), where h and T ′, selected to be 7 and 100, respectively, are the parameters for the windows of the frequency and time domains]. If the interaction is not significant, we conclude that the sequence is a uniformly modulated process and proceed to examine the “between times” term. The sequence is accepted as having a constant power spectrum provided the between times term is not significant.
Using the particular windows described above and data that had been simulated to correspond closely in structure and length to the experimental data, we found that the analysis was oversensitive. Based on these simulation results, we used a value of χ^{2}corresponding to P value of 0.0025 as significant for this test (see discussion). In the example in Table 3, the value of χ^{2} for interaction + residual results in aP value of 0.0074, which implies that the sequence analyzed can be accepted as uniformly modulated. However, the value of χ^{2} for between times is high, which results in aP value close to zero and means that the sequence cannot be accepted as having a constant power spectrum.
All the statistical calculations were done by using Splus statistical software (Statistical Sciences UK, Oxford, UK).
ModelDependent Statistical Analysis
Two simple linear timeindependent approaches to modeling the sequences were employed. The first is the simple autoregressive moving average (ARMA) model. The second is a simple model in statespace form. Both of these models have a timeindependent structure. Both of the models assume there is no exogenous input. These assumptions are explored further in both results and discussion.
Simple ARMA model fitting.
The general form for an ARMA model is
The coefficients of ARMA models were estimated by using Splus statistical software.
Statespace model fitting.
As an alternative to specifying the noise processes in terms of an ARMA model, the noise processes can also be modeled in statespace form
In this statespace form, provided the system parameter f,R
_{v}, and R
_{w}are all known, the system state and variance (P ) can be predicted from the measurement y(t) through a Kalman filter as follows (2)
However, in the present study, both R
_{v} andR
_{w} are unknown and they can not be identified separately using the scalar observationy(t). Therefore, we define P ′ =P/R
_{w}. The procedure for prediction and update then becomes
The initial values for x(0‖0) andP ′(1‖0) were set arbitrarily at 0 andR _{v}/R _{w}, respectively. The initial transient lasts for only a few steps of iteration.
By minimizing the prediction error ∑_{t} [ y(t) −x̂(t‖t − 1)]^{2}, the system parameter f along with the ratioR _{v}/R _{w}can be estimated. The sequences obtained using the Kalman filter forv(t) and w(t) [v̂(t) = x̂(t + 1‖t + 1) − fx̂(t‖t) andŵ(t) = y(t) −x̂(t‖t)] are such that one is simply a scalar multiple of the other. It should be noted that this procedure provides a maximum likelihood estimate for the parameters and is different from the common procedure of using an extended Kalman filter to estimate parameters on line (14).
The parameters of the statespace model were estimated using a least squares method provided in the Numerical Algorithms Group Fortran library (NAG; Oxford, UK), subroutine E04FDF.
RESULTS
ModelIndependent Statistical Analysis
Specific periodicities.
The smoothed periodicities for the six data sets for one subject and protocol (subject 811, exercise) are shown in Fig.1. These show that the power is concentrated in the low frequencies, suggesting the presence of correlation within the data sets. No particular peaks are apparent that are consistent across the six data sets. Furthermore, in general, a smooth decline in power can be traced within the 95% confidence intervals. These findings were generally true for all the individual data sets, and following this observation the periodograms were averaged to produce a single periodogram for each subject and protocol; the results are shown in Fig. 2. Again, smooth lines for the power can be traced within the 95% confidence intervals, suggesting the absence of marked specific periodicities within the data.
Evolutionary spectral analysis.
Evolutionary spectral analysis was used to determine whether there were changes occurring over time in the power spectrum. The results are shown in Table 4. When the values for the interaction + residual term are examined, they show that 24 of these 30 data sets for rest can be accepted as uniformly modulated but only 7 out of the 24 could be accepted as having a constant power spectrum over time. For exercise, 26 out of the 30 data sets could be treated as uniformly modulated, but only 6 of these 26 as having a constant power spectrum.
ModelDependent Statistical Analysis
The results from the evolutionary spectral analysis suggest that simple linear timeinvariant models are unlikely to describe the data fully. Nevertheless, we start with these models to investigate what aspects of the data the models fail to fit and also to determine whether there are aspects of the data that these models can describe well despite inadequacies elsewhere.
Simple ARMA model fitting.
Three kinds of loworder timeseries models were applied to each individual data set. They were AR_{1}, a single autoregressive term; AR_{2}, two autoregressive terms; and AR_{1}MA_{1}, a single autoregressive term together with a moving average term. The result of each model fitting was tested by assessing the whiteness of the residuals using a portmanteau test. The numbers of data sets with white residuals for each subject and protocol are listed in Table 5. The AR_{1}MA_{1} model is the most satisfactory model for both protocols. This model described 22 of the 30 resting data sets and 24 of 30 exercise data sets with overall white residuals.
Examples of the autocorrelation functions of the original data sequences and the residuals after different model fittings for two typical runs are plotted in Fig. 3. Thetop panel shows an example where only the AR_{1}MA_{1} model is satisfactory in the sense that the residual sequence is white, and the lower panel shows another example, which was not satisfactorily described by any of the three different models. It should be noted that all the three models are successful in removing a considerable degree of the correlation present in the original data sequence and that, even when the models fail to fit, the autocorrelation functions are not too far from white when compared with the original data sequence.
The mean and SD values for each coefficient of the three fitted models for each subject and protocol are listed in Table6. In general, these coefficients are not too different between protocols and subjects, although the values forsubject 807 at rest are unusual. The autocorrelation function for this subject at rest appeared much closer to white than in all other subjects and protocols.
The residuals of the AR_{1}MA_{1} model were tested by using evolutionary spectral analysis to determine whether their power spectra were constant. The results are shown in Table 4. The findings were very similar to those before any model had been fitted. This suggests that the model does not describe this feature of the data.
Statespace model fitting.
As an alternative to the approach of using ARMA models, one statespace model of the noise processes was investigated.
The residuals were tested for whiteness in the same manner as for the residuals from the ARMA model. Overall, 22 out of 30 resting data sets were white and 26 out of 30 exercise data sets were white. The results are listed in Table 5. This is broadly similar to the findings with the AR_{1}MA_{1} model.
The parameters (means ± SD) of statespace model for each subject and protocol are listed in Table 7. Theoretically, the statespace model specified reduces to an AR_{1}MA_{1} model with certain restrictions on the values of the a _{1} and c _{1}coefficients of the AR_{1}MA_{1} model ( and Refs. 10, 14). This means that the values fora _{1} and c _{1} coefficients can be calculated from the statespace model, and these values are also given in Table 7. In most cases, these parameter values were similar to those obtained from the AR_{1}MA_{1} model fitting (Table6).
Finally, the residuals were tested by using evolutionary spectral analysis to determine whether the power spectrum was constant. The results are shown in Table 4 and are very similar to those obtained from the original data sequences, showing that the fitting of the model has not affected this feature of the sequences.
Changes in ventilatory variance.
The residuals from the ARMA and statespace models are mostly overall white, but evolutionary spectral analysis shows that the residual sequence does not have a constant power spectrum throughout. This result, coupled with the finding that many of the sequences may be classified as uniformly modulated by evolutionary spectral analysis, suggests that the variance of the sequence may be slowly changing.
Assuming that the mean value of the data sequence remains constant (equal to zero) and that the only factor causing the variation in power spectrum is a changing variance, the “uniformlymodulated” sequence X(t) can be written as
where X
_{0}(t) is a timeinvariant process and C(t) is the modulating function. We chose to try to construct a modulating functionC(t), using an autoregressive estimate of the variance as follows
Examples for C(t) for each subject for both the rest and exercise protocols are shown in Fig.4. Based on the average values for breath durations, this function may be considered to have a time constant of 66 s for the data obtained at rest and 47 s for the data obtained during exercise.
After demodulating the original sequence by C(t), we obtain the demodulated sequence. Evolutionary spectral analysis on this demodulated sequence was then used to determine whether the variation in power spectrum had been removed by demodulation. The results are shown in Table 4.
The number of data sets that had a constant power spectrum rose from 7 and 6 before demodulation to 24 and 26 after demodulation for resting and exercising data, respectively. It is possible that these figures could rise further with a different choice of demodulation function. Finally, it is worth noting that similar results are obtained if demodulation is applied to the residuals from the statespace model instead of the original data sequences (Table 4).
Suitability of models.
Both the ARMA models and the statespace model used above are linear timeinvariant models. However, the results from evolutionary spectral analysis show that in most sequences the power spectrum is varying, and this may suggest that the process underlying the generation of the sequence is nonstationary and/or nonlinear. This variation is clearly not described by these models. The question thus arises as to whether these models are at all useful for describing the data sequences.
One approach to answering this question is to apply the models to the data after demodulation, which removes, or at least attenuates, the variation in the power spectrum from the sequence. There was little difference between the number of residual data sets accepted as white before demodulation compared with after demodulation (Table 5). Also, the coefficients from the ARMA model fitting on the demodulated data (Table 6, bottom) are similar to the values obtained before demodulation (Table 6, top). Analysis of variance performed on the AR_{1}MA_{1} model coefficients revealed no significant differences between the values obtained with the modulated and demodulated data.
These findings suggest that the autocorrelation in the sequence can usefully be described by a simple model despite the presence of a slowly varying variance. The statespace model is equivalent to the AR_{1}MA_{1} model with certain restrictions, and similar results for whiteness were obtained before and after demodulation (Table 5).
DISCUSSION
The major findings of this study are as follows. 1) Both a simple ARMA model and a simple linear timeinvariant statespace model for the breathtobreath correlations present in the ventilatory data gathered by using an endtidal forcing system to hold endtidal gas tensions constant can produce residuals that, in the majority of cases, are overall white. 2) Neither the simple ARMA models nor the simple linear timeinvariant statespace model fully describe the data because the power spectra associated with the residuals sequences vary with time. 3) In a sizable proportion of cases, this variation in power spectrum can be explained by a process of uniform modulation of the ventilatory data or residuals. 4 ) One chosen modulating function was shown to be effective in reducing these data sequences to ones with constant power spectra in the majority of cases.
Comparison with Previous Studies
In general, it is quite difficult to draw direct comparisons between the results from our study and results from previous studies. First, by using data gathered with the endtidal forcing technique, our study has concentrated on the properties of the respiratory controller, whereas most other studies have examined breathbybreath variations when the chemical feedback loop is intact. Such a feedback loop might be expected to affect the correlational structure very significantly. Second, studies have varied in terms of the background level of stimulation employed (if any), the species employed, and whether the study was conducted with subjects awake, asleep, or under anesthesia.
Specific periodicities.
In our study, the results from spectral analysis failed to reveal any sustained specific periodicities in the data. In conscious humans, by displaying the spectral densities for sequential subsets of the original data series, Jensen (11) failed to find any peaks or troughs that were persistently present throughout the length of the series. This is also in keeping with a study of human subjects during stage 2 sleep (15). These results support the absence of specific periodic components in the data series. However, in the present study, the results of evolutionary spectral analysis generally show that the power spectrum is not constant, and thus none of the above studies can exclude the presence of brief periods of periodicity that are not present throughout the data sequence.
Evolutionary spectral analysis.
Most studies have not considered whether the data sequence retains a constant power spectrum throughout. In conscious humans, Jensen (11) calculated power spectra for sequential subsets of ventilatory data sequences but did not test whether they differed significantly from one another. Ackerson et al. (1) commented that data sequences were often “nonstationary” when assessed by using a “likelihood ratio test,” but did not provide any results to compare with those in this study. Nevertheless, in their study, Ackerson et al. developed an adaptive timeseries model for describing nonstationary ventilatory variations. This model has many similarities to the modulation model proposed in the present study.
ARMA models.
The results from our study suggest that, of AR_{1}, AR_{2}, and AR_{1}MA_{1} model structures, the AR_{1}MA_{1} structure was the most satisfactory for describing the breathbybreath correlation in the data sequence. We have not found another study that has investigated the AR_{1}MA_{1} structure for ventilatory sequences. Jensen (11) reports values in humans for coefficients from multivariate AR_{1} and AR_{2} models for tidal volume, inspired time, and expired time pooled from a range of different experimental conditions. All coefficients for the lagged effects of each variable on itself are positive. Modarreszadeh et al. (15) examined an AR_{2} model for describing ventilatory variability in sleeping humans. They found that the AR_{1} coefficient was significantly positive but that the AR_{2} coefficient did not differ significantly from zero.
In the anesthetized cat, Benchetrit and Bertrand (3) used an “isolated respiratory centre” preparation to examine the statistical properties of the phrenic nerve output. This was done so as to circumvent the influence of the chemical feedback system, rather as we have attempted with the endtidal forcing system in our present study. The study employed a multivariate AR_{1}structure and found a positive autocorrelation coefficient, suggesting that there was indeed breathtobreath correlation arising within the respiratory controller. The firstorder model structure was not always adequate.
Khatib et al. (13) used an AR_{1} model structure to compare the correlation between tidal volumes in anesthetized spontaneously breathing rats, with the correlation between the integrated phrenic neurogram in anesthetized vagotomized artificially ventilated rats. In hyperoxia and hypercapnia, they found positive autoregressive coefficients for the tidal volume data in the spontaneously breathing rats but not for the integrated phrenic neurogram in the artificially ventilated rats. From these results, they conclude that the autocorrelation in tidal volume arises from the presence of noise within the chemical feedback system rather than as a function of the respiratory controller. They suggest that the difference between their results and those of Benchetrit and Bertrand (3) could either be due to species differences or to Benchetrit and Bertrand failing to eliminate experiments in which there was nonstationarity. The results from our study in conscious humans accord much more closely with those of Benchetrit and Bertrand than with those of Khatib et al. (13).
Statespace models.
We are unaware of any studies in which a statespace formulation has been used for the noise processes in the absence of additional deterministic input stimuli. However, there have been studies of different noise models in experiments involving step changes in Pet _{CO2} in both anesthetized cats and in conscious humans (8, 9). These studies do not report any of the parameter values associated with the stochastic part of the model, but the example autocorrelation functions suggest that both a “process” and a “parallel” noise model lead to a considerable whitening of the residuals when compared with the “measurement only” noise model.
Effects of exercise.
We are unaware of any previous studies that have provided a comparison of these models between the states of rest and exercise. In general, the results from the two states were surprisingly similar given the differences in the mean ventilations, the associated differences in respiratory cycle durations, and the somewhat different mental state that may exist between rest and exercise. Inspection of Tables 6 and 7suggests that there may be some increase in the degree of autocorrelation observed in exercise compared with rest, but this did not reach significance when the values of the coefficients were compared for the AR_{1}MA_{1} and statespace models.
Evolutionary Spectral Analysis, Stationarity, and Time and StateDependent Models
Evolutionary spectral analysis.
Evolutionary spectral analysis is a means of assessing whether the power spectrum of a data sequence is constant by estimating the power spectrum at certain specified points in time. This has been applied as a test of whether the underlying data sequence is stationary (18, 20). This section discusses the utility of evolutionary spectral analysis for this purpose, together with what may be inferred about the underlying model structure from this analysis. It also considers the reliability of the analysis for the particular form and length of data sequence that has been used in this study.
A series is (weakly) stationary if the first and second moments are constant over time (i.e., the expected value of x, E(x) = μ and the variance of xVar (x) = ς^{2}). If there are multiple independent realizations of a time series, then stationarity may be tested for by estimating the mean and variance at each point in the sequence. However, in other cases, multiple realizations may not exist, as is the case with our data. In this case, it is necessary to test for stationarity by sampling the series at intervals that are sufficiently far apart for the samples to be considered independent. Evolutionary spectral analysis is one method based on this approach. This leaves the question of how far apart do the samples have to be in order to be sufficiently far apart? The problem is that very slow variations in the structure of a time series may occur (especially with nonlinear series) despite the fact that the series is actually still stationary. Thus it is safer to interpret results from evolutionary spectral analysis in terms of whether there is longterm variation in a sequence than whether the process is truly stationary. Indeed, the autoregressive model used to estimate the slowly changing variance in this study is a timeinvariant model, despite the fact that it is being used to describe “nonstationary” behavior detected by evolutionary spectral analysis.
There is an interesting form of duality between the ideas of nonstationarity and nonlinearity (20). A general stationary statedependent model of order (k, l) may be written as
where X _{t} is the state at timet; e _{t} is a white sequence;x _{t} denotes the state vector,x _{t} = (e _{t − l − 1}, … ,e _{t},X _{t − k + 1}, … ,X _{t})′; φ_{u} and ψ_{u} are AR and MA coefficients, respectively; and μ is local mean.
This form is a stationary nonlinear model. However, for a single realization, it is also possible to think of the statedependent coefficients φ_{u}(x _{t}), ψ_{u}(x _{t}) as timedependent coefficients φ_{u}(t), ψ_{u}(t) via the definitions φ_{u}(t) ≙ φ_{u}(x _{t}) and ψ_{u}(t) ≙ ψ_{u}(x _{t}). In this formulation, the model may be thought of as a linear nonstationary model (19). This duality illustrates clearly the difficulty of distinguishing between linear nonstationary models and stationary nonlinear models when there is only a single realization of the data sequence.
The success, or otherwise, of using evolutionary spectral analysis to classify time series must depend on the type and length of the actual sequences used together with the particular forms of window chosen. To explore this further, 60 AR_{1}MA_{1} data sets were generated, each set matching one of the real data sets for AR_{1} and MA_{1} coefficients and for record length. Evolutionary spectral analysis was then undertaken on these sequences. Using a significance level of P < 0.05 for the χ^{2} test, over 20% of the sequences were incorrectly classified as arising from a process generating sequences without a constant power spectrum. However, modulation of each of the sequences with a modulating function that was the same as that which had been determined for the matching experimental data sequence and then repetition of the evolutionary spectral analysis on the modulated sequences generally resulted in levels of significance for the χ^{2} test that were very much higher than those for the unmodulated sequences. When a value for χ^{2}corresponding to P < 0.0025 was used, it resulted in only 4 of the 60 stationary data sets being misclassified as not having constant power, but only 7 of the 60 modulated data sets were classified as being of constant power. Furthermore, it must be remembered that in a number of cases the modulating function may have been very “mild.” The results from this simulation determined the choice of significance level used for the χ^{2} test in this study.
Adequacy of simple linear timeinvariant models.
A major result of this study is that the correlation between breaths can be modeled by simple linear timeinvariant ARMA or statespace models but that there remains a longerterm variation in the variance of the sequence. One question that naturally arises is to what extent can we separate these two features of the sequence and, in particular, to what degree may a longer term variation in ventilatory variance affect the estimation of the coefficients describing the shorter term correlational structure? A partial answer to this question is found inresults, where the coefficients of the AR_{1}MA_{1} model were estimated both before and after demodulation. No significant differences were found. However, it might be argued that this result arises because both the modulation and the autoregressive coefficients were estimated from the same data sequence.
To overcome this objection, a simulation study was performed. A set of 10 data sequences was constructed with a known underlying AR_{1}MA_{1} structure, with the coefficient values equal to the mean values for each subject and protocol. First, the coefficients were estimated from these simulated data by using the same technique as used for the real data sequences. Second, the simulated data sequences were modulated by using a randomly selected modulating function, one from each subject and protocol. Then the coefficients were reestimated by using these modulated sequences. Finally, the values estimated from both the original and modulated sequences were compared. Comparison of the parameter values between the modulated and unmodulated data sets revealed an average deviation for thea _{1} coefficient of 0.8 ± 8.35 (SD) % and for thec _{1} coefficient of 5.9 ± 20.08%. Analysis of variance revealed no significant differences for either coefficient. These results support the notion that the shorter term correlations in the data sequence can be estimated in the presence of the slower variations in the variance of the sequence.
Influence of CO_{2} on Data Sequences
A major intent of this study was to examine the openloop characteristics of the respiratory controller with chemical levels of stimulation held constant. However, the experimental technique of dynamic endtidal forcing is imperfect, and small fluctuations in Pet _{CO2} and Pet _{O2} remain. Because all experiments were conducted at a mean Pet _{O2} of 100 Torr, it is most unlikely that the variations in Pet _{O2} will have any effect on the ventilatory sequence because the respiratory controller’s sensitivity to changes in Pet _{O2} around this mean level will be extremely low (6). The same may not necessarily be true for Pco _{2}, since the ventilatory sensitivity to changes in Pco _{2} is much higher and, indeed, crosscorrelation analysis reveals a significant effect of Pet _{CO2} on minute ventilation (V˙e) in 30 out of 60 data sets (portmanteau test for cross correlation). Thus one question that arises is whether our results reflect the combined properties of the human respiratory system coupled to our dynamic endtidal forcing apparatus rather than simply the openloop properties of the respiratory controller.
To try to answer this question, a multivariate timeseries approach has been adopted. This was necessary because not only may the Pet _{CO2} affect the ventilation but the breathbybreath fluctuations in ventilation will affect Pet _{CO2}. Many different model structures are possible when using multivariate analysis, and to make progress we decided to focus on one simple model that has some structural basis in terms of the known properties of both the physiological and the dynamic endtidal forcing systems.
The Pet
_{CO2} of a breath is determined (assuming that CO_{2} delivery to the lungs via the blood remains constant) by 1) the Pet
_{CO2} of the previous breath;2) the inspired Pco
_{2} of the current breath; and 3) the ventilation of the current breath. Thus we may write
where h is a general function and Pi
_{CO2} is inspiratory Pco
_{2}. The endtidal forcing system uses a predictioncorrection scheme, with the correction based on an integralproportional controller structure. We may write for this
where subscripts p, d, and m are predicted, desired, and measured gas tensions, respectively, and g_{p} and g_{i} are the proportional and integral feedback gains, respectively.
In our experiments, Pet
_{CO2}
_{ d} is constant, and after a period of time the integral term should be relatively constant. Thus Pi
_{CO2}(n) is really just a linear function of Pet
_{CO2}
_{ m}(n − 1). Hence, we may write Eq. 22
as
where h′ is a general function. For small fluctuations in Pet
_{CO2}, we linearize around the means and write
This now gives a simple expression for the way ventilation may affect CO_{2} in our system.
The Pet
_{CO2} may affect ventilation either through the peripheral chemoreceptors or through the central chemoreceptors or both. A reasonable estimate for the median lag to the peripheral chemoreceptors is two breaths (5). Estimates for the lag for the central chemoreceptors vary, but a value of three to four breaths might be reasonable (7). To keep the model simple, we chose to use a lag of two breaths, which should model a peripheral response and still be capable of modeling at least some of any central response that is present. When this is combined with our AR_{1}MA_{1} model for ventilation and the model derived above for Pet
_{CO2}, we get an entire multivariate model of the form
This model was fitted to our data by using the NAG library routine G13DCF with and without the coefficient a _{12} being fixed at zero, to examine the effects of either incorporating or not incorporating the effects of Pet _{CO2}on ventilation. The exclusion of the effects of Pet _{CO2} on ventilation reduced the average value of coefficient a _{11} by 7% and reduced the average value of coefficient c _{11} by 14%. Both these changes were significant when using analysis of variance. Thus we conclude that when the variations in Pet _{CO2} are ignored, it causes a modest underestimation of the degree of correlation associated with the respiratory controller.
After fitting the multivariate model, the cross correlations between the residuals for ventilation and the Pet _{CO2} were calculated. There was no significant effect of Pco _{2} onV˙e in 51 out of 60 data sets, confirming that the effect of fluctuations in Pet _{CO2} onV˙e had been modeled satisfactorily for most data sequences.
Generalization of StateSpace Model to Incorporate a Time and StateDependent Form
The results have shown that the degree of variation in ventilatory variance observed with a constant level of stimulation is such that it does not affect the fitting of the correlational structure unduly. However, it is not clear whether this would still be the case if the level of stimulation were varying because, in this case, the variation in ventilatory variance is likely to be substantially greater (cf. variances for rest and exercise in Fig. 4). In this sense, the results of the current study might be seen as rather limited when applying models to situations in which the level of stimulation is varying (i.e., there is a deterministic input present, rather than the purely stochastic processes considered in the present study).
As an alternative to the linear timeinvariant statespace model, a more general statespace model can be written as follows
where ξ(t) and η(t) are noise processes with state and timedependent variances
This model allows the variances to vary with amplitude of ventilation and with time through the functiong[x(t), t]. This is known as a heteroscedastic form (10). The advantage is that providedg[x(t), t] is fully determined at time t, then the same Kalman filter can be used to estimate the coefficients much as before, since the term ing drops out while the ratio of the variances is estimated. The only differences are that the prediction and update equations for the variance P should be defined by using R
_{ξ}and R
_{η} instead of R
_{v}and R
_{w}
where
is the Kalman gain at time t + 1, and thatP ′ is defined as P ′ =P/R _{η} for the estimating procedure. The ratio R _{ξ}/R _{η} equals the ratio R _{v}/R _{w}and may be estimated in the estimation procedure.
The question now remains, is this heteroscedastic form likely to be an adequate model for ventilatory variability? Two possible limitations are 1) the same function g is used to modulate variances of both noise sequences, and 2) the functiong has to be fully determined (i.e., nonstochastic) at timet. For the scalar output of ventilation, the first of these limitations is unimportant as ξ(t) and η(t) are not separately identifiable and the sequencesξ̂(t) and η̂(t) are scalar multiples of each other. Consequently, the same function φ(t) may be constructed to modulate bothv(t) and w(t) to produce ξ(t) and η(t), respectively.
The importance of the second of these limitations is less clear theoretically. However, if g(t) is varying sufficiently slowly with respect to the observations, then a good estimation of g(t) may be obtained from prior observations. The observation in this study that many sequences were uniformly modulated, combined with the practical construction of the demodulation function C(t) for the residual sequences, suggests that such a function does exist (or at least a good approximation exists).
For use in parameter estimation, the heteroscedastic noise model is no more complicated than the ordinary form, since the algorithms associated with each are the same. Thus, for a parallel noise structure, the heteroscedastic noise model can be combined with deterministic models in the same way as an ordinary noise model. The deterministic and stochastic components are simply added together. Because the algorithms are the same for both the ordinary and heteroscedastic models, parameter estimates obtained by using these models are the same in each case. The important point is that the heteroscedastic model provides the necessary theoretical extension to deal with noise processes that are not of constant power but that are, nevertheless, uniformly modulated.
In summary, the heteroscedastic structure allows the variation in ventilatory variance to be introduced into the model in a general way, without specifying the precise form that it takes. Although this does not contribute further to understanding the origins of the variation in ventilatory variance, it does provide a useful means to incorporate this feature into overall models of ventilatory control and still provide a structure that is amenable to parameter estimation.
Acknowledgments
P. A. Robbins thanks many colleagues, including Drs. Berkenbosch, DeGoede, and Marriott, and Profs. Swanson, and Ward for all the discussion that has helped to give rise to this paper.
Appendix
Relationship between StateSpace Form and AR_{1}MA_{1} Form
The first section of the derives the appropriate ARMA representation of the statespace model employed in this study. The second section derives the relations between the parameters of the statespace model and those of the ARMA model.
ARMA Representation of StateSpace Model
The statespace model with noise processes and zero input employed in this study may be written as
From Eq. EA2, we may write
Relationship between Parameters of AR_{1}MA_{1}and StateSpace Model
Define z(t) as
Also, from Eqs. EA6
and
EA7
After equating Eqs. EA10
and
EA13
and eliminatingR
_{z}, we get
In addition to the above result, we may also calculate E[z(t)z(t − 1)]. FromEq. EA8
we obtain
Footnotes

This study was supported by the Wellcome Trust. P. J. Liang was supported by a Run Run Shaw Scholarship.

Address for reprints requests: P. A. Robbins, University Laboratory of Physiology, Parks Road, Oxford OX1 3PT, UK.
 Copyright © 1996 the American Physiological Society