Abstract
Liang, PeiJi, Daphne A. Bascom, and Peter A. Robbins.Extended models of the ventilatory response to sustained isocapnic hypoxia in humans. J. Appl. Physiol. 82(2): 667–677, 1997.—The purpose of this study was to examine extensions of a model of hypoxic ventilatory decline (HVD) in humans. In the original model (model I) devised by R. Painter, S. Khamnei, and P. Robbins (J. Appl. Physiol. 74: 2007–2015, 1993), HVD is modeled entirely by a modulation of peripheral chemoreflex sensitivity. In the first extension (model II), a more complicated dynamic is used for the change in peripheral chemoreflex sensitivity. In the second extension (model III), HVD is modeled as a combination of both the mechanism of Painter et al. and a component that is independent of peripheral chemoreflex sensitivity. In all cases, a parallel noise structure was incorporated to describe the stochastic properties of the ventilatory behavior to remove the correlation of the residuals. Data came from six subjects from a study by D. A. Bascom, J. J. Pandit, I. D. Clement, and P. A. Robbins (Respir. Physiol.88: 299–312, 1992). For model II, there was a significant improvement in fit for two out of six subjects. The reasons for this were not entirely clear. For model III, the fit was again significantly improved in two subjects, but in this case the subjects were those who had the most marked undershoot and recovery of ventilation at the relief of hypoxia. In these two subjects, the chemoreflexindependent component contributed ∼50% to total HVD. In the other four subjects, the chemoreflexindependent component contributed ∼10% to total HVD. It is concluded that in some subjects, but not in others, there may be a component of HVD that is independent of peripheral chemoreflex sensitivity.
 sustained isocapnic hypoxia
 hypoxic ventilatory decline
 peripheral chemoreflex sensitivity
the ventilatory response to sustained hypoxia is biphasic: there is an initial rapid increase in ventilation, which is then followed by a slow decline (7, 18). This slow decline is known variously as hypoxic ventilatory decline or depression (HVD) or hypoxic ventilatory “rolloff.”
In anesthetized animals, HVD may be independent of the peripheral chemoreflex sensitivity, since 1) in anesthetized cats, the peripheral chemoreflex sensitivity to CO_{2} is unchanged during HVD induced by central hypoxia (16); 2) in anesthetized cats, the rapid ventilatory response at the onset and relief of hypoxia is almost symmetrical (4); and 3) respiratory depression during hypoxia has been a consistent finding in anesthetized animals that have undergone peripheral chemodenervation (12).
In conscious humans, however, HVD may arise predominantly through an alteration of peripheral chemoreflex sensitivity. The evidence for this is as follows: 1) the rapid ventilatory response at the relief of sustained hypoxia is much smaller than the rapid response at the onset of the period of hypoxia (6, 7, 9); 2) the reintroduction of hypoxia after a brief relief from sustained hypoxia will cause a ventilatory response with a smaller magnitude compared with the first exposure to hypoxia (3, 8); and 3) brief exposures to additional hypoxia during the development of HVD elicit progressively smaller ventilatory responses as the HVD develops (2).
On the basis of the above findings, Painter et al. (13) developed a mathematical model to describe HVD during sustained hypoxia in humans in which HVD arises as the result of a progressive decline in peripheral chemoreflex sensitivity. The ventilatory response to sustained hypoxia and the asymmetrical magnitudes of the ventilatory on and offtransients at the onset and relief of hypoxia are both described well by this model. However, as noted by Painter et al., this model fails to describe very well the ventilatory behavior after the step out of hypoxia in cases in which there may be some undershoot and subsequent recovery of the ventilatory response. One possible cause is that the assumptions about the dynamics are too simple. Another possible cause of the undershoot and recovery sometimes observed at the relief of hypoxia might be the existence of an additional component of HVD in conscious humans, which is independent of peripheral chemoreflex sensitivity, as has been described in the studies relating to anesthetized animals.
The purpose of the present study was to explore extensions to the model of Painter et al. First, we wished to explore whether the dynamics of the model were too simple. Second, we wished to explore whether there is an additional component for HVD that is not related to the peripheral chemoreflex sensitivity. In particular, we wished to determine 1) whether either extension would improve the fit of the model to the data at the relief of hypoxia and 2) in the case of the model with the additional component for HVD, what proportion of HVD would be attributed to the component that is independent of the peripheral chemoreflex sensitivity.
In the original model of Painter et al., the parameter estimation was performed with a simple least squares fitting technique (this corresponds to a measurementonly noise model). With respiratory data, use of this technique generally results in residuals that are correlated, and this can make statistical comparisons between different models difficult. To overcome this, a parallel noise structure was used in conjunction with the deterministic models to remove this correlation. The noise model was fitted by using a Kalman filter algorithm.
METHODS
Data
Data were taken from a previous study performed in our laboratory (3). Each experiment lasted 40 min. Endtidal Po _{2}(Pet _{O2}) was held at 100 Torr for the first 10 min, at 50 Torr for the next 20 min, at 100 Torr for the next 5 min, and at 50 Torr for the final 5 min. Endtidal Pco _{2}(Pet _{CO2}) was kept constant throughout the experiment at 2–3 Torr above the subject’s natural value. Minute ventilation (V˙e) was measured breath by breath. Data were available from six subjects, with three to six repeats of the protocol on each subject.
Both the strength of the stimulus (Pet _{O2} = 50 Torr) and background Pet _{CO2} were the same as for the data used by Painter et al. (13). However, the data differ in the sense that the protocol used to generate the data modeled in the present study was prolonged to allow for the reintroduction of hypoxia to provide a second ontransient. In the data modeled by Painter et al., there was no second ontransient. This second ontransient may help to discriminate between a component of HVD that arises through modulation of peripheral chemoreflex sensitivity and a component that does not act via such a mechanism.
Models
Model of Painter et al. (13) (model I).
The model of Painter et al. can be written in the form
In the original formulation, Painter et al. (13) used a parameter k_{h}, whereas in the present formulation the parameter g_{100} has been used. The relationship between these two parameters is given by
There were two reasons for preferring our revised parametrization of the model by Painter et al. First, g_{100}, the hypoxic sensitivity after exposure to hyperoxia, is a more intuitive parameter than k_{h} in a physiological sense. Second, g_{h} and k_{h} appeared to be highly correlated in the study of Painter et al., whereas there appeared to be much less correlation between g_{h} and g_{h} + k_{h} ( g_{100}, see Table 1 in Ref. 13).
Assuming that both saturation S and gain g
_{p}remain constant throughout the breath, these differential equations can be solved to provide a set of difference equations that can be written as
Model incorporating secondorder dynamics for g_{p} (model II).
In the model of Painter et al. (13), a firstorder dynamic was chosen for g
_{p}. To explore whether this choice of dynamic was too simple, a more complicated dynamic was chosen to determine whether such a change could provide a significant improvement in the fit. The particular dynamic chosen was secondorder, and, in a sense, it represents the next level of sophistication that could be introduced. The model can be written as
Assuming, as before, that saturation and g_{p} remain constant throughout a breath, the difference equations become
Model incorporating additional component for HVD (model III).
The original model of Painter et al. (13) describes the whole of the hypoxic ventilatory decline as arising from a modulation of peripheral reflex sensitivity. The object of this extended model is to introduce a component of HVD (V˙_{d}) that may be independent of the peripheral chemoreflex. This model can be written as
The time delay for this component is assumed to be the same as for the sensitivityrelated component d_{p}.
Fitting Technique
Simple least squares parameter estimations.
The approach to estimating parameter values adopted by Painter et al. (13) was to use a simple least squares fitting technique. With this technique, the objective function to be minimized is
Parameter estimation by using a Kalman filter.
One problem associated with estimating the parameters of the deterministic models is that successive breaths are not independent (14). On the basis of a previous study of breathing during steady levels of stimulation (10), the stochastic behavior was modeled by using a parallel noise process of the form
and the prediction and update of a modified variance (10) are given by
where
is the Kalman gain at the (n + 1)th breath.
The physiological parameters of the deterministic models along with the two parameters of the stochastic model ( f, R_{v}/R_{w}) can now be estimated by minimizing the sum of the squared prediction errors of the model, which is given by
Minimization process.
The simple least squares objective function (Eq. 20 ) and the objective function for the Kalman filter (Eq. 29 ) were minimized by using a standard Numerical Algorithms Group (Oxford, UK) subroutine (E04FDF) that finds an unconstrained minimum of a sum of squared residuals. The fitting procedure was repeated over a range of fixed values for d_{p} for each data set. The value for d_{p} corresponding with the lowest sum of squared error was regarded as the optimal value for d_{p} for the data set.
In the study from which the data are taken (3), the authors report that the speed of the initial ventilatory response to the rapid changes of Pet _{O2} appeared to be different at the different stages in the protocol. Thus, instead of fitting a single value for the fast time constant T_{p}, separate values for the time constant were fitted for the onset, relief, and the second ontransient of the hypoxic stimulus (T_{p1}, T_{p2}, and T_{p3}, respectively).
Goodness of fit.
The improvement of the fit with the extended models, as compared with the original model, was assessed using an Fratio test on the sum of squared prediction errors (1). The F statistic is given by
Parameter values for each subject.
Overall parameter values for each subject for each model were obtained by using a technique employed by Swanson and Bellville (15) in their study of modeling the ventilatory responses to CO_{2}. First, an idealized stimulus of 5min euoxia (Pet
_{O2} = 100 Torr), 20min hypoxia (
Sensitivity Analysis
To determine whether the parameters of the model were theoretically identifiable, a sensitivity analysis was undertaken.
Model data were generated for the input function associated with the protocol. The more complex extended model (model III) was used with a set of physiologically reasonable parameter values (actually, the parameter values determined for subject 758) and a breath duration of 3 s. The model parameters were then estimated repeatedly for the data set, with the exception that one of the parameters would be fixed, in turn, to the value used to generate the data and then to values ranging from 80 to 120% of this value (except for k_{p}, where a range of 0.0–0.02 was used, since its value is close to zero). The sum of squared residuals could then be plotted against the value of the fixed parameter. If the parameter is not identifiable, then the plot should be flat, whereas if the value of the parameter affects the fit, then the plot should have a definite minimum (zero) at the parameter value used to generate the data.
RESULTS
Sensitivity Analysis
The results for the sensitivity analysis are shown in Fig.1. For all parameters, a minimum (zero) for the sum of squared residuals can be observed at the parameter value used to generate the data set. This indicates that all the parameters are theoretically identifiable in the model with the input function used to generate the data. However, the sensitivity of the sum of squared residuals to variations in parameter value varies considerably for the different parameters, and in some cases the sensitivity plots are quite asymmetric. The importance or lack of it of these features is difficult to gauge.
Fitting Technique
The effect on the residuals of incorporating the noise model in the Kalman filter technique is illustrated in Fig.2. A large degree of autocorrelation is seen in the residuals from the simple least squares fitting process, which appears to have been almost totally removed by incorporating the parallel noise structure. Results from the portmanteau test demonstrated that all 31 sets of residuals from the simple least squares fitting algorithm were significantly nonwhite, with Pvalues close to zero, whereas only 5 out of 31 sets of residuals from the model using a Kalman filter algorithm were found to be nonwhite. These five sets of nonwhite residuals were still much less correlated than the ones from the simple least squares fitting technique.
Average parameter values from both the simple least squares fitting technique and the Kalman filter algorithm are shown for the model by Painter et al. (13) in Table 1. In general, the parameter estimates from the two techniques appear quite similar. However, for the time constants associated with the rapid ontransients (T_{p1}, T_{p3}), the Kalman filter algorithm consistently yielded lower values than did the simple least squares fitting procedure. For the time constant for the hypoxic ventilatory decline T_{h}, the reverse was true, with the values from the Kalman filter algorithm being consistently higher than those from the simple least squares fitting procedure. Consistent effects across all the six subjects are significant at the level P < 0.05 (sign test) and suggest some element of bias in estimating these values with a simple least squares fitting routine.
Model Comparison
Figure 3 shows one example of the respiratory data together with the model fits for the three models. (Model fits for the model by Painter et al. (13) when using the simple least squares fitting procedure are not shown, as the appearance was essentially identical to the appearance of the results for the deterministic component with the use of the Kalman filter algorithm.) Figure 3, left, shows the fit for model I developed by Painter et al. (13), Fig. 3, center, shows the fit formodel II, which has secondorder dynamics for g_{p}, and Fig. 3, right, shows the fit for model III, which incorporates the additional component of HVD, independent of chemoreflex sensitivity. Figure 3, bottom left (which shows the deterministic component of the output for model I) illustrates clearly for this particular data set that the simpler model of Painter et al. can describe most of the experimental record reasonably well, but that the model does not reflect particularly well the transient undershoot in ventilation at the relief of the stimulus. This aspect of the fit does not seem to be improved by model II, but the introduction of the extra chemoreflexinsensitive component of model III greatly improves the fit to this part of the record.
Figure 4 shows average data and model fits for each subject for model I, Fig.5 shows these data for model II,and Fig. 6 shows these data for model III. These were calculated from the individual data sets and model fits by interpolating the data every 3 s and then averaging across all the individual responses for each subject. The average of the residuals sequences for each subject was calculated in the same way, and the 95% confidence intervals (calculated as ±1.96 SE) were plotted along with the averaged ventilatory output. For much of the response, the outputs of all three models are very similar. However, for the two subjects (subjects 758, 796) who show the greatest undershoot and recovery at the relief of hypoxia, model III appears to describe this process more accurately than do the other two models.
When model II is compared with model I, Fratio tests on the sum of squared prediction errors revealed that there was no significant improvement in fit for four out of the six subjects but there was an improvement in fit for the other two (subjects 766, 824). Careful inspection of Figs. 4 and 5 suggests that there may have been some improvement in the fit during the ontransients, but no improvement was observed during the recovery from hypoxia, where visual inspection suggests the fit is least good.
When model III is compared with model I, Fratio tests revealed that there was, again, no significant improvement in fit for four out of the six subjects but that for the two other subjects (subjects 758, 796) model III fitted significantly better than model I. Visual inspection of Figs. 4 and 6suggests that for these two subjects the undershoot and recovery in ventilation were most marked and that model III did produce an apparent improvement in fit to this portion of the response.
A whiteness test (portmanteau test) on the prediction errors showed that 26 out of 31 sequences could be accepted as white for model I, 27 out of 31 sequences for model II, and 28 out of 31 sequences for model III. Thus fitting a stochastic component to the data makes the use of an Fratio test possible without the complications of trying to correct the degrees of freedom for correlation present within the residuals.
Overall parameter values for each subject for each model are shown in Table 2. For all three models, the time constants for the rapid ventilatory response to the onset and the relief of hypoxia (T_{p1}, T_{p2}, and T_{p3}) are faster than the time constant for the slow change in the peripheral chemoreflex sensitivity T_{h}. In most cases, the time constant for the rapid ventilatory response at the relief of hypoxia (T_{p2}) appears faster than those for the rapid ventilatory response at the onset of hypoxia (T_{p1} and T_{p3}). However, this effect does not reach statistical significance overall.
For model II, in five out of six subjects, the parameter values for T_{I} and T_{h} are of the same order of magnitude (for three subjects, the values for T_{I} and T_{h} are the same). In only one subject (subject 796) are the values radically different. In this one subject, the very low value for T_{I} compared with T_{h} means that the dynamic for g _{p} is very close to first order, but for the other five subjects this seems not to be true. This alteration in dynamic has had consistent (albeit often small) effects on the other parameter values of the model. In particular, when model II is compared with model I, V˙_{c} is consistently increased and g_{100} and g_{h}consistently decreased across all the subjects.
Turning to the parameters for model III, some effect was attributed to the extra component of the model by the fitting process in five out of six subjects. The values estimated for both g_{d} and T_{d} were quite variable among the different subjects. Table 3 shows the percentage contributions to HVD of the component acting via modulation of the peripheral chemoreflex sensitivityV˙_{p} and the component that is independent of the peripheral chemoreflex (V˙_{d}) for each subject. These figures were obtained by calculating the steadystate values for V˙_{p} andV˙_{d} under conditions of a sustained steady Pet _{O2} of 100 Torr and under conditions of a sustained steady Pet _{O2} of 50 Torr. For V˙_{p}, this was done by setting the derivatives of Eqs. 2 and 3 to zero and obtaining an explicit expression for V˙_{p}, and, forV˙_{d}, by setting the derivative of Eq.18 to zero and obtaining an explicit expression forV˙_{d}. The averaged contributions calculated across all the six subjects were 73.8% for the component acting via the peripheral chemoreflex and 26.2% for the peripheral chemoreflexinsensitive component.
When we compared parameter values across all models, one striking feature we noticed was the consistent nature of the differences in value for k_{p}. For five out of six subjects, the value for k_{p} for model II was smaller than for model I, and for five out of six subjects, the value for k_{p}was smaller for model III compared with model II. Formodel I, a positive value for k_{p} was obtained for all six subjects; for model II, a positive value was obtained for five out of six subjects; but for model III, a positive value was only obtained for one out of the six subjects. Physiologically, k_{p} represents the component of peripheral chemoreflex drive that is present in high oxygen, and this suggests that model I attributes a higher proportion of the total ventilation to peripheral chemoreflex drive in euoxia than doesmodel II, and that model II attributes a higher proportion than does model III. This is confirmed by inspection of Table 3, where the contributions of each component to ventilation have been calculated for each of the models in euoxia. The mean contributions to ventilation by the peripheral chemoreflex drive are significantly different among models, being highest with model I and lowest with model III (analysis of variance,P < 0.05). In the case of model III, the reduction of total ventilation by V˙_{d} is small in euoxia, as would be expected by the fact that the saturation is close to 1.
DISCUSSION
Fitting Technique
In the study by Painter et al. (13), no model of the noise structure was employed, whereas in the present study both a simple least squares fitting procedure and a simple statespace model for a parallel noise process were explored. A suitable noise model allows the inherent correlation between successive breaths to be included in the model and results in residuals that are overall white or close to white. This is useful when comparing the overall quality of fit between different models. This approach has been adopted by others for modeling the ventilatory response to CO_{2} (5). The ability of the particular statespace model chosen to describe ventilatory variability has been demonstrated previously (10), and this result is supported by the observation that 26–28 out of 31 sets of residuals could be accepted as white on the basis of a portmanteau test in the present study, whereas all the residuals’ sequences from the simple least squares fitting technique were highly correlated.
It is of some interest to determine whether the values for the stochastic component of the model ( f, R_{v}/R_{w}) were similar to those obtained from a previous study under conditions of steady chemical stimulation (10). From this previous study, the mean ± SE values for f and R_{v}/R_{w} were 0.69 ± 0.06 and 0.95 ± 0.13, respectively, for the resting data and 0.81 ± 0.03 and 0.44 ± 0.11, respectively, for the data during exercise. The mean values for f and R_{v}/R_{w} in the present study were 0.76 ± 0.04 and 1.55 ± 0.41, respectively, formodel I; 0.75 ± 0.03 and 1.85 ± 0.53, respectively, formodel II; and 0.72 ± 0.03 and 1.74 ± 0.53, respectively, for model III. Thus the mean values for f in the present study lay between the values for rest and exercise from the previous study, whereas the values for R_{v}/R_{w}appeared to be a little higher. When the values obtained from the fits for the three different models of the present study were compared, analysis of variance showed that there was no significant difference in the estimated value of R_{v}/R_{w}between the three models (P > 0.05), but there was a significant difference between the estimated values of f for the different models (P < 0.05). One interpretation of this is that when the model structure varies, the residual correlation that has to be explained by the stochastic process ( f ) varies too, but not necessarily the estimate for R_{v}/R_{w}.
Model Comparison
Comparisons among the models may be drawn in a number of ways. First, the question may be asked whether the extended models produce any statistically significant improvement in the fit, compared with the original model of Painter et al. (13). This can be determined both by examining the overall sum of squared residuals and by inspection of the confidence intervals for the ensembleaveraged residuals over the time course of the protocol. Second, some physiological judgment may be exercised in terms of whether the estimated parameter values are likely to be consistent with other known features of the respiratory system.
In statistical terms, for model II, two subjects (subjects 766, 824) showed a significant reduction in the sum of squared residuals, as compared with model I. These subjects do not appear to have any particularly distinguishing features in terms of their response to the hypoxic protocol. In addition, a direct visual comparison of the model fit between model II and model I did not suggest much by way of difference between them. Overall, it is difficult to be certain exactly why these two subjects show a significant improvement in fit, but we feel that it is most likely to be related to differences during and immediately after the onset of hypoxia. The model does not appear to improve the fit to the period immediately following the relief from hypoxia, which was one of the specific purposes of the study.
For model III, again only two subjects (subjects 758, 796) showed a significant reduction in the total sum of squared residuals when compared with model I. These two subjects consistently showed an apparent undershoot and then recovery of ventilation at the relief of hypoxia in all the individual experimental repeats of the protocol. In other subjects, an undershoot followed by recovery was either not observed or appeared inconsistently in some, but not all, of the repeats. Inspection of the plots of the model fit and ensembleaveraged residuals over time shows that it is during the period following the relief from hypoxia that the models are most distinct. We conclude that the use of model III may significantly improve the fit in some subjects, who show a reasonable degree of undershoot and recovery of ventilation after the relief of hypoxia, but that this is not the case for all subjects.
The two subjects in whom the fit of model III was a significant improvement over the original model of Painter et al. (13) had the two largest magnitudes for the gain term g_{d} of the additional chemoreflexindependent component. In these two subjects, slightly over 50% of HVD (mean 56%) was attributed to the peripheral chemoreflexindependent component. In the other four subjects, the value was very much smaller, with a mean value of 11%. Again, this emphasizes the importance of intersubject differences.
For the two subjects with a substantial contribution to ventilation from the peripheral chemoreflexindependent compartment, the time constants associated with this compartment were rather fast (22 and 50 s) and are consistent with the comments of Bascom et al. (3) that the baseline ventilation appears to recover more quickly than does the peripheral chemoreflex sensitivity. These values do not really correlate well with the sorts of values previously obtained for the development of HVD in anesthetized animals [onset of HVD 149.7 ± 8.5 (SE) s; offset of HVD 105.5 ± 10.1 s (17)]. Thus, on the basis of data from these two subjects, it would appear that this peripheral chemoreflexinsensitive component cannot really be equated simply with a chemoreflexinsensitive form of HVD observed in anesthetized animals, as postulated in the introductory part of this paper. However, for the other three subjects in whom some component of the ventilatory response was attributed to the peripheral chemoreflexindependent term, the time constants were all somewhat longer than those associated with HVD in anesthetized animals. Overall, we feel that it would be unwarranted to draw any very firm conclusions from these values.
One notable feature that differs among the models is the parameter estimates determined for k_{p}. Physiologically, k_{p} is a parameter relating to the contribution of the peripheral chemoreflex to ventilation in hyperoxia when the arterial blood is fully saturated. We consider that the positive values for k_{p} in the original model of Painter et al. (13) may be necessary to fit the slightly lower level of euoxic ventilation seen after the relief of hypoxia, when compared with the original baseline ventilation. In model III, this feature of the data can be modeled via the additional term V˙_{d}. Which of these approaches in a physiological sense is correct is difficult to determine. Although there is some evidence to suggest that the peripheral chemoreflex is inactivated by hyperoxia in humans (11), in which case model III might represent the physiology better, it is also possible that positive values of k_{p} might simply represent some nonlinearity of the response of the peripheral chemoreflex to saturation in this range.
Acknowledgments
This work was supported by The Wellcome Trust. P.J. Liang was supported by a Run Run Shaw studentship.
Footnotes

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