A Comparative Data-Based Modeling Study on Respiratory CO2 Gas Exchange during Mechanical Ventilation

The goal of this study is to derive a minimally complex but credible model of respiratory CO2 gas exchange that may be used in systematic design and pilot testing of closed-loop end-tidal CO2 controllers in mechanical ventilation. We first derived a candidate model that captures the essential mechanisms involved in the respiratory CO2 gas exchange process. Then, we simplified the candidate model to derive two lower-order candidate models. We compared these candidate models for predictive capability and reliability using experimental data collected from 25 pediatric subjects undergoing dynamically varying mechanical ventilation during surgical procedures. A two-compartment model equipped with transport delay to account for CO2 delivery between the lungs and the tissues showed modest but statistically significant improvement in predictive capability over the same model without transport delay. Aggregating the lungs and the tissues into a single compartment further degraded the predictive fidelity of the model. In addition, the model equipped with transport delay demonstrated superior reliability to the one without transport delay. Further, the respiratory parameters derived from the model equipped with transport delay, but not the one without transport delay, were physiologically plausible. The results suggest that gas transport between the lungs and the tissues must be taken into account to accurately reproduce the respiratory CO2 gas exchange process under conditions of wide-ranging and dynamically varying mechanical ventilation conditions.

The goal of this study is to derive a minimally complex but credible model of respiratory CO2 gas exchange that may be used in systematic design and pilot testing of closed-loop end-tidal CO2 controllers in mechanical ventilation. We first derived a candidate model that captures the essential mechanisms involved in the respiratory CO2 gas exchange process. Then, we simplified the candidate model to derive two lower-order candidate models. We compared these candidate models for predictive capability and reliability using experimental data collected from 25 pediatric subjects undergoing dynamically varying mechanical ventilation during surgical procedures. A two-compartment model equipped with transport delay to account for CO2 delivery between the lungs and the tissues showed modest but statistically significant improvement in predictive capability over the same model without transport delay. Aggregating the lungs and the tissues into a single compartment further degraded the predictive fidelity of the model. In addition, the model equipped with transport delay demonstrated superior reliability to the one without transport delay. Further, the respiratory parameters derived from the model equipped with transport delay, but not the one without transport delay, were physiologically plausible. The results suggest that gas transport between the lungs and the tissues must be taken into account to accurately reproduce the respiratory CO2 gas exchange process under conditions of wide-ranging and dynamically varying mechanical ventilation conditions. Keywords: respiratory cO2 gas exchange, data-based modeling, closed-loop mechanical ventilation control inTrODUcTiOn It is anticipated that autonomous closed-loop controlled mechanical ventilators will be increasingly used in the future to enhance the safety of mechanical ventilation. Automation will enable standardized treatment protocols as well as fill the gap between the increasing demands versus the limited number of respiratory experts. First, it is estimated that the failure to use recommended respiratory interventions in the intensive care unit (ICU) results in 170,000 preventable deaths per year in US (Pronovost et al., 2004). Thus, closed-loop-controlled mechanical ventilators can be an attractive option to translate the research knowledge into clinical practice via automatic knowledge transfer, which can potentially reduce errors, inappropriate interventions, and heterogeneity of knowledge and practice. Second, the number of prolonged mechanical ventilations (>96 h) was projected to have more than doubled in 2020 compared with that in 2000 (Zilberberg et al., 2008). It was also forecasted that the shortage of medical personnel with expertise in mechanical ventilation will begin in 2007 and will worsen thereafter (Angus et al., 2000). In this regard, autonomous and closed-loop-controlled mechanical ventilators can be a viable solution to guarantee quality care without exacerbating the medical personnel's workload. Despite their promising potential, closed-loop mechanical ventilation controllers have not yet penetrated, at least as much as anticipated, into clinical practice. This may be attributed to many reasons [such as conservatism and inertia (Cabana et al., 1999;Dent and Goldberg, 1999;Rubenfeld et al., 2004)], but a critical challenge may have been the concerns raised on the safety of these autonomous systems. Indeed, proving the validity and safety of closed-loop-controlled mechanical ventilators in preclinical testing is not trivial (see the limited number of subjects used in the testing of closed-loop mechanical ventilation controllers reported in the literature summarized in Table 1). In fact, even the new modes available in the commercialized mechanical ventilators [such as PAV (Younes, 2002), NAVA (Sinderby et al., 1999), ASV , and SmartCare™ (Dojat et al., 1992)] have not been extensively evaluated (Verbrugghe and Jorens, 2011;Cordioli et al., 2013;Rose et al., 2013).
Mathematical models of physiologic system have been promoted as viable alternative to preclinical testing. Indeed, computational models may potentially be widely used to examine and assess the safety of a range of medical devices and systems, as suggested by the guidance document recently released by the U.S. Food and Drug Administration (FDA) on the reporting of computational modeling studies in the medical device submissions (CDRH and FDA, 2014). To rigorously assess the safety, efficacy and robustness of closed-loop mechanical ventilation controllers, a mathematical model of respiratory physiology is required that possesses several desired characteristics: (1) it must reproduce a realistic respiratory response to mechanical ventilation with accuracy for trustworthy assessment of mechanical ventilation controllers, (2) it must be low-order with minimum number of parameters for tractability in controller design and testing, and (3) it must be physiologically transparent to streamline the interpretation of the testing results. However, models of respiratory physiology reported in today's literature do not fulfill all these characteristics. Most importantly, existing models are typically too complex to be used in Monte-Carlo simulation-based testing (Chiari et al., 1997;Anderson et al., 2003;Wolf and Garner, 2007;Cheng et al., 2010).
In pursuit of the ultimate goal of realizing a computational model of respiratory physiology applicable to the design and testing of closed-loop mechanical ventilation controllers, the goal of this study is to derive a minimally complex but credible model of respiratory CO2 gas exchange that may be used in systematic design and pilot testing of closed-loop end-tidal CO2 controllers in mechanical ventilation. We first derived a candidate model that captures the essential mechanisms involved in the respiratory CO2 gas exchange process. Then, we simplified the candidate model to derive two lower-order candidate models. We compared these candidate models for predictive capability and reliability using experimental data collected from 25 pediatric subjects undergoing dynamically varying mechanical ventilation during surgical procedures.
This paper is organized as follows. Section 2 outlines the model structure, its data-based modeling procedure, and the comparative analysis approaches. Section 3 describes the experimental data and data analysis methods. Section 4 presents and discusses the results. Section 5 concludes the paper with future directions. resPiraTOrY cO 2 gas eXchange MODel Model structure selection Respiratory CO2 gas exchange process during mechanical ventilation involves the lungs, arteries, and veins as well as body tissues (Grodins et al., 1954;Khoo et al., 1982;Batzel and Tran, 2000) ( Figure 1). CO2 gas is produced in the tissues as a consequence of metabolism and is subsequently transported to the lungs via venous blood (whose CO2 tension is thus relatively high). CO2 gas is then excreted in the lungs by alveolar ventilation. The arterial blood, whose CO2 tension is relatively lower than the venous blood, is delivered back to the tissues to collect CO2 from them.
In this study, we model the lungs and the tissues as homogeneous compartments, while approximate the delivery of CO2 gas between them via arteries and veins as pure transport delays. The CO2 tension in the lungs (alveolar CO2 tension) is determined by the balance between CO2 supply from the venous blood and CO2 removal by the alveoli. The amount of CO2 gas supplied from the venous blood to the lungs depends on the difference between alveolar ( ) C ACO 2 versus mixed venous ( ) C vCO 2 CO2 concentrations as well as cardiac output (Q), while the amount of excreted CO2 gas is determined by alveolar CO2 tension ( ) P ACO 2 and alveolar ventilation ( )  V A . Therefore, the rate of change in P ACO 2 can be written as follows: where VL is the effective lung volume, P ICO 2 is the CO2 tension in the inspired air, τ1 is the transport delay between the tissues and the lungs, and λ = 863 mmHg lBTPS/lSTPD. The CO2 tension in the tissues ( ) C vCO 2 is determined by the production of CO2 due to the body's metabolic activity ( )  V CO 2 and its removal from the tissues via arterial blood. The amount of CO2 gas removed from the tissues depends on the difference between C ACO 2 versus C vCO 2 as well as Q. Therefore, the rate of change in C vCO 2 can be written as follows: where VB is effective tissue volume in the body and τ2 is the transport delay between the heart and the tissues. By assuming that CO2 concentration and tension are related to each other via the Henry's law (Khoo et al., 1982;Lumb, 2005): where α = 0.0065 lSTPD/mmHg l, β = 0.244 lSTPD/l, Eqs 1 and 2 are reduced to the following with P ACO 2 and P vCO 2 as state variables: and noting that P ICO = 0 2 , yields the following state variable representation that dictates the respiratory gas exchange process during mechanical ventilation: and the transport delays are the unknowns to be derived from the data. It is noted that Q and  V CO 2 are regarded as constants (see Limitation and Future Perspectives for details).
We intend to conduct data-based modeling using minute ventilation ( )  V and end-tidal CO2 tension ( ) P etCO 2 data. To this aim, we assume that (1)   V V » A (i.e., negligible dead space ventilation; see Limitation and Future Perspectives for details) and (2) P P etCO ACO 2 2 » (Burton, 1966;Coles et al., 1973;Ritchie et al., 1987;Williams and Babb, 1997;Brunner, 2002) (see Limitation and Future Perspectives for details). Since Eq. 5 cannot be identified solely based on  V and P etCO 2 because x2 is not accessible, it is re-formulated into the following regression on ẍ1: Note that only τ, but not τ1 and τ2 individually, appears in Eq. 6 because only x1 is measured. Both P etCO 2 and  V data are given by discrete-time sequences. Thus, the regression Eq. 6 is discretized using the forward difference approximation (Nakamura, 1993) where Ts is sampling interval and τ is restricted to positive integers to be compatible with the discretization. Since the regression vector in Eq. 7 can be constructed solely based on x1 and u, Θ  θ θ θ θ { } and τ may be identified by fitting Eq. 6 to P etCO 2 and  V . We considered two avenues in deriving lower-order models from Eq. 4. First, noting that transport delays were not taken into account in many previous studies [see, e.g., Melo et al. (1993), Olofsen et al. (2010), Wang et al. (2010), and Karbing et al. (2011) versus Grodins et al. (1967), Khoo et al. (1982), Batzel and Tran (2000), Beda et al. (2010), and Dunn and Whiteley (2010)], a lower-order model was derived by neglecting τ in Eq. 4. This results in the following regression model (called M2 hereafter): where Ω 2 Second, Eq. 8 can be further simplified by aggregating the compartments associated with the lungs and the tissues. If we assume that gas exchange and metabolism occur in a single compartment representing the lungs and the tissues altogether, the following regression model results from Eq. 4 by setting x1 = x2 and τ1 = τ2 = 0, and then combining the two equations in Eq. 4 (called M1 hereafter): where Ω 1 1

Model identification
Since Eq. 7 is a pseudo-linear regression model due to the dependence of the regression vector ψ3 on τ, it cannot be solved via standard linear least-squares method. In addition, its parameter vector Ω3 contains the products of unknowns θ1θ3, θ2θ3, and θ2θ4, which presents a challenge in deriving θi, i = 1, …, 4 from Ω3 uniquely. Further, standard least-squares method is not ideal in maximizing the predictive capability since it minimizes one-step-ahead prediction error instead of pure (i.e., infinite-step-ahead) prediction error. To cope with these challenges, we identified Θ and τ by solving the following optimization problem: where {Θ * ,τ * } is the set of optimal model parameters, N is the number of data samples used to solve Eq. 10, x1(k) is the P etCO 2 data at sample time k, while x k Once Θ is determined, the physical respiratory parameters can be derived as follows. First, VL is derived as = α θ . Thus, VL, VB, Q, and  V CO 2 can be uniquely derived from Θ. The parameters in M1 and M2 can be determined by solving optimization problems similar to Eq. 10 in order to maximize predictive capability, so that an objective comparison can be made on the three candidate models.

comparative Model analysis
The three candidate models were compared in terms of predictive capability and reliability. First, the predictive capability was measured in terms of the root-mean-squared error (RMSE) between note that τ * = 0 for M1 and M2): Second, the reliability was measured in terms of the asymptotic variance (Ljung, 1999), which represents the expected parametric variance estimated by the prediction error and the parametric sensitivity: is the prediction error variance associated with  ε k ( ), and SN(Θ * )is the inverse sensitivity covariance matrix. Finally, the Akaike's Information Criterion (AIC) (Burnham and Anderson, 1998) was used to compare the overall quality of the candidate models by rewarding the predictive capability and penalizing the model complexity simultaneously: where K is the number of parameters in the model. The first term in Eq. 14 denotes the predictive capability (i.e., the goodness of fit), while the second term reflects the model complexity. The third term is to compensate for finite sample size (Burnham and Anderson, 1998).

MeThODs experimental Data
The P etCO 2 and  V data used in this study were extracted from a large and anonymized physiologic database, which was constructed out of a standardized set of clinical data collected as part of a larger investigation which was approved by the Children's and Women's Health Centre of British Columbia. We selected 25 anonymous pediatric subjects receiving pressure-controlled mechanical ventilation during a surgical procedure. These subjects were selected based on the evidence of a large change in the ventilation setting made by the caregivers, manifested by their variability in respiratory rate (RR), peak inspiratory pressure, and I:E ratio (the ratio between inspiratory time versus expiratory time), rendering the corresponding data appropriate for data-based modeling analysis. The data were collected as trend values every 5 s using a standard respiratory module (M-CAiOVX, Datex-Ohmeda, Finland) and then saved via a custom-built, centrally located data collection software.

Data analysis
The P etCO 2 and  V data thus recorded were used to solve the optimization problems (e.g., Eq. 10 for M3) to derive the candidate models for each subject. A constrained optimization routine in MATLAB's Optimization Toolbox was used. Multiple initial parameter values were employed to assure that the solution converged to the global minimum. The model-predicted P etCO 2 was produced by applying u V =  and the requisite initial conditions on x P 1 2 = etCO to Eq. 11. In deriving M3, the optimization problem was repetitively solved while the value of τ was varied via exhaustive search within a domain that was specified a priori based on the physiologically relevant values reported in the literature [see, e.g., Khoo et al. (1982) and Batzel and Tran (2000)]. In this way, multiple candidate solution sets (i.e., {Θ * ,τ * }) were derived for different τ. Then, the optimal solution set was determined as the one associated with the minimum J value.
Once the optimal data-based models (M1, M2, and M3) corresponding to the 25 subjects were derived and the corresponding model-predicted P etCO 2 produced, the candidate models associated with each subject were compared with each other using RMSE, asymptotic variance and AIC, which were computed by Eqs 12-14, respectively. The number of data samples that was used to compute RMSE, asymptotic variance, and AIC was dependent on subjects, since the amount of data available from each subject was different (see Table 2). The significance in difference among the candidate models was assessed as follows. First, the difference in predictive capability was analyzed using the paired t-test. Second, the difference in reliability associated with each element in Θ * * * * *  θ θ θ θ 1 2 3 4 , , , { } was analyzed by applying the paired t-test to the asymptotic variance corresponding to each element in Θ * . The difference was regarded as significant in case p < 0.05. Third, in comparing the candidate models by AIC, the frequency in which M1, M2, and M3 attained the minimum AIC values was counted across the 25 subjects in order to assess how many times each of the candidate models was suggested as the best model among them.

resUlTs anD DiscUssiOn
A parsimonious and credible model of respiratory mechanics and physiology may expedite the design, analysis, testing, and deployment of closed-loop mechanical ventilation controllers. As an initial step toward such a computational model, the goal of this study was to derive a minimally complex and credible respiratory CO2 gas exchange model applicable to the design and pilot testing of closed-loop end-tidal CO2 controllers during mechanical ventilation. The novelty of this study is that we have rigorously compared respiratory CO2 gas exchange models with different degrees of complexity. We have evaluated the reliability and physiologic transparency of the respiratory parameters as well as the model accuracy in reproducing the relationship between minute ventilation and end-tidal CO2, which may offer unique value to the data-based computational modeling of respiratory physiology.

Data
On the average, 24.4 min length of data was available from each subject ( Table 2). The ventilation settings were largely altered to yield considerable changes in P etCO 2 and  V, as evidenced by Table 2 and Figure 2. In particular, Figure 2 shows that more than 30% change in  V (from its within-subject mean value) was made in 24 subjects, resulting in more than 29% change in P etCO 2 in those subjects (again, from its within-subject mean value).

Predictive capability and reliability
Overall, the two-compartment models (M2 and M3) exhibited superior predictive accuracy than their one-compartment counterpart (M1). Compared with M1, M2, and M3 achieved 35 and 45% reductions in RMSE, respectively, which were statistically significant ( Table 3). A representative model-predicted P etCO 2 in response to  V is shown in Figure 3. In this subject, both RR (7-13 1/min) and tidal volume (TV; 31.2-56.9 ml/min) were varied to yield a large change in  V and P etCO 2 . This example indicates that M2 and M3 can reproduce the relationship between P etCO 2 and  V with high fidelity, whereas M1 clearly exhibits deficiency. Due to its largely degraded predictive capability compared with      M2 and M3, the remaining investigation was devoted to the comparison between M2 and M3.
The asymptotic variance (expressed as asymptotic standard variation) associated with θ 1 * , θ 2 * , and θ 3 * was reasonably small in both M2 and M3, which suggests that these parameters were adequately identified. Considering that θ 1 * , θ 2 * , and θ 3 * are made up of V L * , V B * , and Q * (which are closely related to the time constants associated with x1 and x2), reliable determination of these parameters may be attributed to the dynamic, widely varying  V that resulted in informative transient components in the P etCO 2 response. On the other hand, the asymptotic variance associated with θ 4 * was relatively large for both M2 and M3. This may be at least in part attributed to the underlying assumption in these models that  V CO = α θ where V B * can be reliably derived from θ 2 * and θ 3 * , uncertainty in θ 4 * may arise from  V CO * 2 . This interpretation is in fact supported by the examination of SN(Θ * ) in Eq. 13. For all the 25 subjects, we found that the smallest eigenvalue of the sensitivity covariance matrix S N − ( ) 1 Θ * associated with M2 and M3 was consistently aligned with the direction of θ 4 * , which means that P etCO 2 is more sensitive to θ 1 * , θ 2 * , and θ 3 * than θ 4 * . The reason for the small sensitivity of P etCO 2 to θ 4 * may be because the corresponding regressor element is constant (i.e., 1), whereas the rest of the elements in ψ2 and ψ3 are widely varying. It may thus be concluded that the relatively large asymptotic variance associated with θ 4 * is caused by the lack of excitation in its direction by the data, which is essentially due to the assumption that  V CO * 2 is constant. The analysis of AIC suggested that M3 is overall superior to M1 and M2, while M2 is superior to M1. The AIC results indicate that transport delay plays a crucial role in reproducing the respiratory gas exchange process under dynamically varying mechanical ventilation conditions.

Physiologic Transparency
The respiratory parameters derived for M2 were largely different from those for M3 with statistical significance (see Table 4). Compared with the respiratory parameter values for adults reported in the literature, the respiratory parameters derived for M3 appeared to be more physiologically plausible than those associated with M2. In particular, the median values of V L * , V B * , Q * , and  V CO * 2 derived for M3 were approximately 72, 31, 22, and 46% of the typical adult values, which appears to be at least qualitatively reasonable considering that the data used in this study were acquired from pediatric subjects. 1 In contrast, the respiratory parameter values derived for M2 were unacceptably small to be physiologically realistic. Further investigation of the models showed that the predictive accuracy of M3 drastically deteriorated when its transport delay was set to 0 (much worse than M2), 1 The weights of the pediatric subjects were not available in the anonymized data. However, the median size of the subjects may be estimated by comparing the values of  V and TV in which suggests that the parameters derived for M2 were actually optimized reasonably well to minimize the prediction error in the absence of transport delay. Therefore, the above observation is an additional evidence to support that wide-ranging changes in mechanical ventilation settings result in emphasized transient responses in P etCO 2 that may not be captured without transport delay.
The range of the respiratory parameters derived for M3 were wide; the interquartile ranges (IQRs) of V L * , V B * , Q * , and  V CO * 2 were approximately 143, 175, 213, and 92% of the respective median values. The results suggest that the subjects were associated with a wide-ranging respiratory variability, though the validity of the absolute values of the respiratory parameters cannot be established. The fact that our model could reproduce the responses of all these subjects without any a priori knowledge of these subjects is notable, because its governing Eq. 5 was derived by explicitly considering physiologic principles rather than to simply fit the data. In sum, the ability to derive subject-specific respiratory parameters suited to each subject is the strength of our model.

limitation and Future Perspectives
This study has a number of limitations that need to be addressed in follow-up studies.
First, the lungs of the subjects analyzed in this study were mostly normal without any lung disease. Therefore, the validity of our model in subjects with respiratory pathophysiology has yet to be evaluated. It may be that our model needs to be expanded to include more physiologic mechanisms relevant to accurately reproduce the CO2 gas exchange phenomena in the presence of lung diseases. Future work must make improvements to our model to develop and validate minimally complex lumpedparameter models of respiratory CO2 gas exchange that can be used for both normal and pathologic conditions. Second, the subjects examined in this study were subject to different anesthesia protocols and surgical procedures. This may not alter the conclusion regarding the validity of our model, because the model could yield physiologically plausible respiratory parameter values that accurately reproduced the data in all the subjects examined in this study. However, the values of the respiratory parameters associated with each subject may have been modestly affected in different ways to capture the influence of different anesthesia protocols and surgical procedures. To solidly establish the validity of our model, it must be tested against data collected under strictly standardized protocols.
Third, several key mechanisms in the model were simplified. For example, CO2 production and cardiac output were assumed constant in the model presented in this study. Noting that the subjects were undergoing surgical procedures, it may be reasonable to assume that CO2 production and cardiac output were stable due to the effect of anesthesia and the end-tidal CO2 was varied mainly by the change in the mechanical ventilator settings. However, it is possible that CO2 production and cardiac output may have changed, which may in turn have had an influence on the model parameters. The validity of the model presented in this study may be limited in subjects undergoing a large change in CO2 production and/or cardiac output. The models in this study also assumed negligible dead space ventilation and pulmonary shunt because accurate subject-specific estimation of dead space ventilation and pulmonary shunt solely based on minute ventilation and end-tidal CO2 data is in general challenging. Though the subjects examined in this study were all healthy, even healthy subjects have non-zero dead space, which may further increase under anesthesia despite the mechanical ventilation. Thus, this assumption may have had an influence on the values of the model parameters. In addition, this assumption may be far less justified in subjects with lung disease, limiting the applicability of our model. Thus, our model must be improved to explicitly incorporate dead space ventilation and pulmonary shunt.
Fourth, the models presented in this study are continuous ventilation models. As such, the tidal nature of breathing was not explicitly incorporated in the models. It has been suggested that tidal ventilation models explicitly incorporating the discontinuity in breathing due to inspiratory and expiratory phases are more desirable than continuous ventilation models when attempting to obtain a rigorous and comprehensive understanding of physiology and pathophysiology in respiratory gas exchange during mechanical ventilation [see, e.g., Hahn and Farmery (2003)] and the references therein]. The model presented in this work was able to reproduce the end-tidal CO2 response to minute ventilation. On the other hand, the validity and utility of this model beyond end-tidal CO2 (e.g., the relationship between x2 and mixed venous CO2 concentration) still need to be investigated. In this regard, care must be taken in using this model in analyzing responses to mechanical ventilation other than end-tidal CO2. At the same time, efforts to extend our continuous ventilation model to include tidal nature of breathing may also be rewarding. Finally, the model presented in this study did not incorporate oxygen (O2) gas exchange. Considering that the primary objective of mechanical ventilation therapy is to achieve adequate arterial oxygenation, a truly viable computational model of respiratory physiology for mechanical ventilation beyond closed-loop endtidal CO2 control must incorporate O2 dynamics. Future work must be conducted to expand our CO2 gas exchange model to enable the reproduction of both O2 and CO2 responses to mechanical ventilation.

cOnclUsiOn anD FUTUre WOrK
In our effort to derive a minimally complex, high-fidelity, and transparent data-based respiratory physiology model applicable to the design and pilot testing of closed-loop end-tidal CO2 controllers, we conducted a comparative data-based modeling analysis of the respiratory CO2 gas exchange process. It was shown that a two-compartment model with transport delay was able to accurately reproduce the end-tidal CO2 tension response to dynamically varying minute ventilation challenge in human subjects. Future work will include further investigation of the model as well as the design of model-based closed-loop end-tidal CO2 controllers.

aUThOr cOnTriBUTiOns
All authors listed have made substantial, direct, and intellectual contribution to the work and approved it for publication.

FUnDing
This material is based upon work supported in part by the Office of Naval Research (ONR) under Grant No. N000141410591. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the ONR.