This was a retrospective study using VitalDB, an open dataset containing intraoperative biosignal data and perioperative clinical information from Seoul National University Hospital (SNUH), a tertiary-level hospital in South Korea [21]. The dataset encompassed 6388 cases of noncardiac surgery, with an average of 2.8 million data points per case, collected from August 2016 to June 2017. Vital signs were recorded during surgery, while pertinent clinical information was retrospectively retrieved from the electronic medical records system.

This study utilized a publicly available, deidentified dataset (VitalDB). The collection and use of this dataset for research purposes were approved by the Institutional Review Board (IRB) of SNUH (H-1408-101-605). The requirement for obtaining individual informed consent was waived by the IRB because the dataset contained only deidentified data. All data were anonymized before inclusion in the database, ensuring patient privacy and confidentiality. No compensation was provided to participants as part of this secondary data analysis.

Case selection was performed based on the following inclusion and exclusion criteria. Surgical cases involving patients meeting the following criteria were eligible: (1) aged 18-80 years, (2) mechanically ventilated under general anesthesia using an oral airway, (3) had ASA physical status classification grade 5 or lower, (4) had ETCO₂ monitoring records, and (5) had at least 1 or more PaCO₂ measurement records obtained during surgery. Each PaCO₂ measurement obtained via ABGA was defined as an event point. PaCO₂ measurements before or after surgery were excluded and not considered event points.

Each event point was defined by a trigger time and an observation window for the prediction model that estimates PaCO₂ based on noninvasive parameters (Figure 1). The trigger time is defined as the time at which an estimation of PaCO₂ is performed based on intraoperative biosignals and perioperative clinical information. The observation window is the time frame during which the biosignals were extracted for model training.

PaCO₂ is measured via ABGA, which evaluates arterial gas concentrations such as carbon dioxide and oxygen. Usually, a small volume of blood is extracted from the radial artery using a syringe and a thin needle. During surgery, the anesthesiologists draw arterial blood via an A-line. The temporal point of the ABGA records might be compromised as a result of the transit time of the specimen to the analyzer or manual data entry delays. The timestamp of ABGA records frequently denotes the moment that analysis results are entered into the database instead of precisely indicating the moment at which blood is drawn from the patient.

To increase the precision of the ABGA timestamp, in this study, we utilized mean arterial pressure (MAP) recordings available in VitalDB. We hypothesized that the ABGA procedure (blood draw via A-line) induces a rapid increase in the MAP due to transient line occlusion or a patient physiological reaction. On the basis of this hypothesis, we employed a z score–based outlier detection approach. Abnormal periods surpassing a z score of 3 in MAP records were identified by observing the 20-minute interval that preceded the originally recorded ABGA timestamp. The nearest such abnormal period was defined as the estimated actual timepoint at which ABGA was performed. The trigger time for feature extraction was set at 60 seconds prior to the newly calculated ABGA timestamp to capture the patient’s physiological state immediately preceding the likely blood draw event. In this phase, the PaCO₂ values that did not have identifiable MAP surge points meeting these criteria within the preceding 20-minute window were excluded because they were not considered reliably timed event points. Of the 6311 potentially relevant PaCO₂ measurements considered for this timestamp adjustment, 720 (approximately 11.4%) were excluded due to the absence of a detectable MAP surge. The median (IQR) difference between the newly estimated ABGA timepoint based on MAP surge and the original database timestamp for the remaining measurements was 34.00 (22.00‐54.00) seconds, with a mean standard deviation (SD) of 47.29 (57.53) seconds (range 0-936.00 seconds), indicating that the original database timestamp often did not precisely reflect the physiological event of blood sampling. While this pressure-surge detection method provides a systematic approach to approximate the blood sampling time based on physiological responses, we acknowledge that its precise accuracy compared to the true sampling time was not formally validated against a gold standard timestamp in this study.

A total of 19 variables were selected based on previous literature and domain knowledge. These variables were classified as follows: clinical information, including age, sex, height, weight, surgical approach, surgery type, and preoperative pulmonary function test (PFT) results; and intraoperative biosignals, including body temperature (BT), heart rate, percutaneous oxygen saturation (SpO₂), minute ventilation from the ventilator (MV), positive end-expiratory pressure (PEEP), peak inspiratory pressure (PIP), plateau pressure (PPLAT), mean airway pressure (MAWP), RR based on capnography, tidal volume (TV), fraction of inspired oxygen (FiO₂), and ETCO₂.

During the 60-second observation window before the trigger time, intraoperative biosignals were extracted. Extreme outliers were eliminated using the IQR method, removing data points outside 3 times of the IQR below the first quartile or above the third quartile. Aberrant points deemed theoretically unacceptable were also removed based on clinical expertise. Median values for each biosignal in the observation windows were then calculated and used as features for the prediction model, along with the corresponding surgical patients’ electronic medical records. This approach simplifies the high-frequency biosignal data into static features for point-in-time PaCO₂ estimation but does not explicitly model temporal dependencies within the observation window or directly track dynamic changes over longer periods.

Furthermore, additional features were established by utilizing feature engineering techniques [20]. The variables for PFT results were reportedly classified into nine classes in the original dataset description; however, due to potential inconsistencies or inadequate representation of minority classes, we regenerated the preoperative PFT as an indicator variable, categorizing it into three grades: obstructive, restrictive, and mixed type [22]. Another feature engineering technique employed in this study involved generating interaction features that comprise combinations of two or more existing variables. The incorporation of interaction features, representing known physiological indices, primarily aimed to improve model interpretability by including recognized clinical parameters, although tree-based models like CatBoost can implicitly capture interactions.

Several interaction features, including TV per kilogram of ideal body weight (TV/IBW), oxygen saturation to FiO2 ratio (SpO₂/FiO₂), PEEP to FiO₂ ratio (PEEP/FiO₂), compliance of the respiratory system (CRS), and rapid shallow breathing index (RSBI), were generated using domain knowledge. Each of these interaction features brings a dimension of clinical relevance that reflects the interaction between multiple aspects of patient respiratory mechanics and ventilator settings. Additionally, this study utilized APCONET, which is an external application programming interface (API) that is able to estimate CO using arterial pressure waveforms as an input feature [23]. Detailed information, such as the unit of data, recording device, and source of data, for all the selected variables is available in Multimedia Appendix 1.

Event points containing missing values in any of the selected predictor variables were removed prior to model training. The listwise deletion approach reduced the number of event points from an initial 5591 (derived from cases meeting inclusion criteria before handling missing data for specific event points) to the final 4651 event points used for analysis. The entire dataset was split into training, validation, and testing sets using the nested cross-validation (7 outer folds, 6 inner folds). This configuration was chosen based on the dataset size and common practices aiming for stable performance estimation [24]. All partitioning procedures were conducted at the surgical case level. For example, all data from a single surgical case belonged exclusively to either the training/validation set or the testing set within a given outer fold. This approach was designed to maintain similar data distributions across subsets while ensuring patient independence between training and testing.

Categorical variables were initially processed using OneHotEncoder for compatibility with some preliminary models tested. Continuous variables were scaled using RobustScaler, which scales data using IQR, making it robust to outliers.

To identify a suitable regression model for this task, several ML algorithms were evaluated in preliminary experiments (see Multimedia Appendix 2). Based on its superior performance in these experiments, CatBoostRegressor was selected for final model development and evaluation. The CatBoost model is a robust and effective library for gradient boosting on decision trees. It is particularly adept at handling categorical features natively using its built-in encoder, which was employed in this study for the CatBoost model, thus not requiring the prior OneHotEncoding for these features when training the final CatBoost model [25].

To develop the final predictive model for PaCO₂, training sets were utilized for training the CatBoostRegressor model, and the corresponding validation set was used for hyperparameter tuning within each outer fold. The training sets were used to train the CatBoostRegressor model, and the validation set was utilized to tune the hyperparameters (see Multimedia Appendix 3). During hyperparameter optimization, we utilized Optuna, a framework designed to streamline the hyperparameter tuning process [26]. This framework facilitates the search for the optimal hyperparameter space configuration for a given model in an efficient manner. Subsequently, the performance of the final model was evaluated using held-out testing sets from the seven outer-fold cross-validation methods on the testing sets.

In this study, estimating PaCO₂ based on noninvasive parameters was approached as a regression predictive modeling task. To assess the effectiveness of the ML-based model, we conducted a comparison with two baseline methods that employed the ETCO₂ value. One method was a simple offset model (ETCO₂+5 mm Hg), while the other method involved utilizing linear regression with ETCO₂ measurements as the sole predictor.

Model performance was assessed using two commonly employed metrics: mean absolute error (MAE) and root mean squared error (RMSE). Additionally, we evaluated model performance across different conditions by establishing subgroups. The subgroups were established based on PaCO₂ levels as follows: hypocapnic (<35 mm Hg), normocapnic (35-45 mm Hg), and hypercapnic (>45 mm Hg) cases. A Bland-Altman plot was used to calculate the limits of agreement and analyze the agreement between the real and predictive values; additionally, the intraclass correlation coefficient (ICC) was utilized to evaluate the relative reliability and consistency of the model compared to actual measurements [27,28].

To assess the clinical utility of the predictive model, the percentage of estimation errors for PaCO₂ was computed by calculating the differences between the real and predicted values. Disparities below 5 mm Hg were deemed highly acceptable, those between 5 and 10 mm Hg were moderately acceptable, and any value exceeding 10 mm Hg was considered unacceptable. The establishment of this threshold was based on prior research and Clinical Laboratory Improvement Amendments recommendations [29].

Furthermore, model interpretability was analyzed using the Shapley additive explanation (SHAP) value [30]. The SHAP values, derived from an additive feature attribution model, succinctly illustrate the impact of the input variables on the model outputs, enhancing the understanding of the model’s decision-making process.

The performance metrics were evaluated by averaging the results from the testing sets across the outer folds and computing a CI of 95%. All experimental and data analysis procedures were conducted in the Python 3.10.12 environment.

This study developed an ML-based model capable of estimating PaCO₂ in mechanically ventilated patients under general anesthesia with greater accuracy than traditional ETCO₂-based methods. Utilizing noninvasive physiological parameters and clinical information, the CatBoost model demonstrated strong overall performance, achieving an MAE of 2.38 mm Hg, an RMSE of 3.26 mm Hg, and an ICC of 0.87, indicating excellent agreement with arterial measurements. Critically, the model significantly increased the proportion of clinically highly acceptable predictions (error<±5 mm Hg) to 90.02%, comparable to 80.43% for a linear regression baseline, and reduced unacceptable errors (>±10 mm Hg) to 1.20% from 2.64%. The model’s superiority was consistent across hypocapnic, normocapnic, and hypercapnic subgroups.

Interpretability analysis using SHAP identified ETCO₂ as the most influential feature, as anticipated. Beyond the dominant contribution of ETCO₂, other parameters such as BT, SpO₂/FiO₂ ratio, age, sex, and CRS were found to be important for refining PaCO₂ estimations. For instance, the model tended to predict higher PaCO₂ values with lower BT, lower SpO₂/FiO₂, lower CRS, and higher age, suggesting it learned complex physiological relationships. These findings highlight the value of a multiparameter approach to capture variability not explained by ETCO₂ alone.

The limitations of relying solely on ETCO₂ for PaCO₂ estimation are well documented. While ETCO₂ provides some insights, the PaCO₂-ETCO₂ gradient is known to be variable and influenced by numerous physiological factors, often exceeding the commonly cited 3-5 mm Hg range in healthy individuals [11,16,31]. Our study corroborates this, finding a median gradient of 7 mm Hg (average 7.57 mm Hg, range −16 to 34 mm Hg), underscoring the unreliability of a fixed gradient assumption. Conventional statistical models, such as multivariable linear regression, offer some improvement but are often constrained by linearity assumptions and may not fully capture the complex, nonlinear interactions inherent in physiological systems [32,33].

In contrast, our ML-based approach effectively models these intricate associations by integrating a wider array of biosignals and clinical data. The ability of ML to learn from these complex patterns without a priori assumptions about relationships has shown promise in various medical prediction tasks [20]. Previous studies have also indicated that factors like surgical techniques and patient positioning can affect the PaCO₂-ETCO₂ gradient [34], further supporting the need for adaptive models like the one developed in this study, which can account for such patient-specific and contextual variability more effectively than simpler methods.

This study possesses several strengths that enhance the credibility and potential impact of its findings. First, the use of VitalDB, a large, publicly available, real-world dataset from a tertiary university hospital, provides a diverse cohort from various noncardiac surgeries, improving the generalizability of our results. Second, model performance was rigorously assessed using nested cross-validation, offering a robust estimate of its predictive capabilities on unseen data. Third, the ML model was benchmarked against two clinically relevant baseline methods, clearly demonstrating its superior accuracy. Fourth, our evaluation encompassed not only standard error metrics (MAE, RMSE) but also reliability (ICC) and clinical utility based on predefined error categories (Table 4), providing a multifaceted view of performance. Fifth, the inclusion of SHAP analysis offers a degree of transparency into the model’s decision-making process, which is crucial for clinical translation. Finally, the exploration of an MAP-based timestamping method, while requiring further validation, represents a novel attempt to address a common challenge in retrospective EMR-based research.

Nevertheless, several limitations should be acknowledged when interpreting the results of this study. First, while the ML model outperformed baselines across all PaCO₂ subgroups, its performance was relatively lower in the hypocapnic and hypercapnic groups compared to the normocapnic group. This may be partly due to data imbalance, as these extreme ranges were less frequently observed. Future work could explore techniques like targeted data augmentation or specialized modeling to address this. Second, our SHAP analysis focused on explaining direct PaCO₂ predictions. As ETCO₂ is inherently a dominant feature, this makes it harder to isolate the specific contributions of other features to the PaCO₂—ETCO₂ gradient. Analyzing this gradient directly would be a valuable future direction. Third, the model relies on point-in-time estimations using median values from a 60-second window, simplifying the rich time-series data available. This approach does not capture temporal trends or predict rapid PaCO₂ changes. Fourth, the MAP-surge-based ABGA timestamping method, though systematically applied, was not formally validated against a gold-standard timing reference. Any imprecision here could introduce noise into the feature-target alignment. Fifth, the listwise deletion approach for handling missing data, which reduced our event points from 5951 to 4651, may have introduced selection bias if the pattern of missingness was not completely at random and reduced the overall sample size available for training. Seventh, estimated CO via an external API did not emerge as a highly influential feature in SHAP plots. This might be attributed to the indirect nature or potential inaccuracies of the CO estimation rather than CO itself lacking physiological relevance. Finally, being a single-institution study, the findings require external validation to ensure generalizability across different settings and populations.

Building on these findings and limitations, several avenues for future research are essential for advancing noninvasive PaCO₂ monitoring. First and foremost, external validation of the ML model in diverse, multicenter clinical settings is crucial to confirm its robustness and general applicability. Second, developing time-series models, such as recurrent neural networks, long short-term memory, and transformers, which can process the continuous stream of biosignals, is a key next step. This could improve accuracy and enable the prediction of PaCO₂ trends and rapid changes. Third, future studies should explicitly investigate the model’s capacity to track longitudinal changes in the PaCO₂-ETCO₂ gradient within individual patients. Exploring the linkage between the predicted PaCO₂-ETCO₂ gradient and critical events like hemodynamic instability could yield clinical value. Fourth, research correlating intraoperative PaCO₂ deviations identified by accurate monitoring with postoperative outcomes would further strengthen the rationale for enhanced continuous monitoring. Finally, the MAP-based timestamping approach warrants further investigation and validation.

This study demonstrated that an ML-based model integrating multiple noninvasive parameters can estimate PaCO₂ with higher accuracy and reliability than traditional ETCO₂-based methods in mechanically ventilated surgical patients. The model shows particular strength in increasing the proportion of highly accurate predictions. While acknowledging the need for further development, particularly in incorporating time-series data and external validation, this work highlights the considerable potential of AI to serve as a valuable supplementary tool for enhancing respiratory monitoring and patient management in the perioperative setting.