On the practical identifiability of a two-parameter model of pulmonary gas exchange

Successful application of mechanical ventilation as a life-saving therapy implies appropriate ventilator settings. Decision making is based on clinicians’ knowledge, but can be enhanced by mathematical models that determine the individual patient state by calculating parameters that are not directly measurable. Evaluation of models may support the clinician to reach a defined treatment goal. Bedside applicability of mathematical models for decision support requires a robust identification of the model parameters with a minimum of measuring effort. The influence of appropriate data selection on the identification of a two-parameter model of pulmonary gas exchange was analyzed. The model considers a shunt as well as ventilation-perfusion-mismatch to simulate a variety of pathologic pulmonary gas exchange states, i.e. different severities of pulmonary impairment. Synthetic patient data were generated by model simulation. To incorporate more realistic effects of measurement errors, the simulated data were corrupted with additive noise. In addition, real patient data retrieved from a patient data management system were used retrospectively to confirm the obtained findings. The model was identified to a wide range of different FiO2 settings. Just one single measurement was used for parameter identification. Subsequently prediction performance was obtained by comparing the identified model predicted oxygen level in arterial blood either to exact data taken from simulations or patients measurements. Structural identifiability of the model using one single measurement for the identification process could be demonstrated. Minimum prediction error of blood oxygenation depends on blood gas level at the time of system identification i.e. the measurement situation. For severe pulmonary impairment, higher FiO2 settings were required to achieve a better prediction capability compared to less impaired pulmonary states. Plausibility analysis with real patient data could confirm this finding. Dependent on patients’ pulmonary state, the influence of ventilator settings (here FiO2) on model identification of the gas exchange model could be demonstrated. To maximize prediction accuracy i.e. to find the best individualized model with as few data as possible, best ranges of FiO2-settings for parameter identification were obtained. A less effort identification process, which depends on the pulmonary state, can be deduced from the results of this identifiability analysis.


Background
Mechanical ventilation is a life-saving intervention in intensive care, maintaining pulmonary function in critically ill patients. Appropriate ventilator settings need to be found by the clinician to ensure both sufficient oxygenation and carbon dioxide removal. Target values for arterial partial pressures of oxygen (PaO 2 ) and carbon dioxide (PaCO 2 ) can be reached by changing inspired oxygen fraction (FiO 2 ) and minute ventilation (MV). Removal of CO 2 and therefore PaCO 2 in the patient is mainly regulated by adjusting MV. In critically ill patients, e.g. patients suffering from acute respiratory distress syndrome (ARDS), high levels of FiO 2 and appropriate PEEP are usually necessary to ensure sufficient oxygenation. Finding the appropriate FiO 2 setting follows a trial-anderror approach that may not only be tedious but also exposes the patient to the potential risk of hypoxia and hyperoxia [1][2][3][4]. Pulse oximetry allows a continuous measurement of peripheral oxygen saturation (SpO 2 ), however this method has limitations in sensitivity and accuracy due to calibration assumptions, optical interference, and signal artifact [5]. Therefore, an invasive blood gas analysis is required at the end of each trial to evaluate the individual effect of a change in FiO 2 accurately.
In mechanical ventilation therapy, both the risk of ventilator induced lung injury (VILI) and the effort to find adequate settings may be reduced if medical decision support would provide recommendations on how to adjust a patient's settings to reach a prescribed treatment goal. In general, decision support can be divided into knowledge based (KDSS) and model based systems (MDSS). KDSS builds on rules of typical i.e. average patient behavior to represent reactions to changes of ventilator settings. In contrast, MDSS that are adapted to patient specific physiologic properties can simulate the individual reaction to changes in therapy settings. Using the inverse model, MDSS therefore suggests individualized ventilator settings by evaluating the approximated physiology of the patient. Parameters of a model contain compact information about the individual patient state and dynamics once they are quantified in a parameter identification process (PIP). Parameter identification requires information from patient measurements often obtained during certain clinical maneuvers. Success and robustness of the PIP strongly depends on the properties of the model to reflect the required dynamics of the patient, the signal quality and amount of data available at the bedside. As the identified parameters are used for forward calculation of the model equations, they directly influence prediction performance of the model. While using multiple and even redundant measurements helps compensating noise induced errors, measurement efforts and applying the necessary maneuvers should not interfere with clinical processes. Thus, those measurements should be kept to the minimum necessary to ensure a robust PIP.
Models of pulmonary gas exchange are able to predict the effect of FiO 2 and MV on PaO 2 and PaCO 2 in the patient. One-parameter models [6,7] usually only consider shunt, i.e. the amount of venous blood that is mixed with the oxygenated blood, to describe a patient's oxygenation status and to predict the effect of an increase of FiO 2 on PaO 2 . However, using only one parameter to describe gas exchange impairments fails at low FiO 2 when mismatches between alveolar ventilation (V) and perfusion (Q) occur. Several studies [8,9] have come to the conclusion that PaO 2 /FiO 2 ratio, usually used to categorize lung impairment, changes with FiO 2 . Thus, besides shunt, mathematical models of gas exchange should either include a parameter to describe oxygen diffusion limitation [10,11] or a parameter to characterize V /Q mismatch [9,[12][13][14]. Latter have shown to reproduce measurements at different oxygenation levels with sufficient accuracy compared to more complex models or MIGET measurements [15]. A two-parameter model of pulmonary gas exchange including shunt and V /Q mismatch has previously been published by Kjaergaard et al. [12]. Karbing et al. [16] evaluated this model with data of severely ill intensive care patients. The model has been found to be identifiable with four pulse oxymetric (SpO 2 ) measurements at different levels of FiO 2 and one blood gas analysis (BGA) providing PaO 2 and PaCO 2 together with the acid-base parameters pH, base excess (BE) and the hemoglobin concentration (cHb) as well as the end-tidal gas fractions of oxygen (FetO 2 ) and carbon dioxide (FetCO 2 ). Although systems have been built to perform the necessary measurements in 10-15 min [17], lowering the number of measurements required for identification and therefore minimizing the required time and effort is highly relevant. Therefore we investigated if the number of measurements that are necessary to identify the model can be reduced to one FiO 2 -setting. Additionally, we evaluated the influence of the chosen level of FiO 2 during model identification.

Gas exchange model with V /Q-mismatch and shunt
The mathematical model of human pulmonary gas exchange consists of two alveolar compartments that are perfused and ventilated and one shunt compartment that is perfused but not ventilated. The alveolar compartments are separated into a compartment with high V /Q-ratio and a compartment with low V /Q-ratio. This allows the consideration and simulation of limitations in gas exchange for both oxygen (O 2 ) and carbon dioxide (CO 2 ) concentrations in blood. Shunt, i.e. the fraction of venous blood not participating in gas exchange, is quantified by model parameter f s multiplied with blood flow Q. 90 % of the non-shunted blood ((1 − f s )*Q) is distributed to the low V /Q compartment, 10 % of the non-shunted blood is delivered to the high V /Q compartment. Model parameter f A represents the fraction of alveolar ventilation V̇A that reaches the low V /Q compartment. Figure 1 shows the model structure of the pulmonary gas exchange model.
The model assumes equilibrium in blood gas concentrations as well as constant alveolar ventilation and perfusion without separating ventilation into phases of inspiration and expiration. Model inputs are FiO 2 as well as end-tidal blood gas fractions of oxygen (FetO 2 ) and carbon dioxide (FetCO 2 ), respectively. Inspired carbon dioxide is set to 0. Tidal volume V tid and respiratory frequency f R are assumed to be constant during simulation and are provided as additional model inputs. Model outputs are the resulting arterial blood gas parameters PaO 2 and PaCO 2 .
Alveolar ventilation V̇A is calculated from f R and the difference between tidal volume V tid and the anatomic dead space volume V ds : FetO 2 and FetCO 2 are composed of alveolar gas fractions FAO 2 and FACO 2 in both compartments, such that Index x represents O 2 and CO 2 in Eq. 2 and in all following equations. Index 1 refers to the alveolar compartment with high V /Q, while index 2 denotes the low V /Q compartment. Oxygen consumption V̇O 2 and carbon dioxide production V̇C O2 are derived from alveolar air flow to each of the compartments and the difference between inspired and alveolar gas fractions as described in Eqs. (3) and (4): Capillary blood gas concentrations Cc x are derived from alveolar gas fractions using O 2 and CO 2 dissociation curves [18,19] (T-temperature): Venous blood gas concentrations are then calculated as described in Eqs. (6) and (7).
Cardiac output Q is measured at the bedside or estimated from the patient's body surface area. Fig. 1 Schematic representation of the gas exchange model. The model consists of two alveolar compartments, one with low and one with high ratio between ventilation V and perfusion Q respectively. Blood flow is distributed among the shunt f s and the two alveolar compartments. The low V /Q-compartment is perfused with a fixed fraction f Q of 90 % of the non-shunted blood and ventilated with fraction f A . FiO 2 describes the fraction of oxygen in inspired air. PaO 2 and PaCO 2 are arterial partial pressures of oxygen and carbon dioxide respectively FAO 2 and FACO 2 are solved numerically for a given measured FetO 2 and FetCO 2 with the condition that venous concentration in both compartments is equal. Finally, arterial blood gas concentrations CaO 2 and CaCO 2 are calculated as: Arterial partial pressures of oxygen and carbon dioxide are then calculated from the reversed dissociation curves:

Model simulation
Forward calculation of the model equations is termed as model simulation. The flowchart of the model simulation process M Θ is depicted in Fig. 2 on the left. Vector Ψ summarizes physical constants for measurements needed for model simulation: Here, minute ventilation MV is calculated as: Vector Θ includes model parameters f s and f A as well as Ψ: Output values PaO 2,sim and PaCO 2,sim are calculated depending on Θ and FiO 2 :

Model identification
Parameters that need to be identified in the presented model are shunt fraction (f s ) and the fraction of alveolar ventilation that is distributed to the alveolar compartment with low V /Q-ratio (f A ). The process of model identification (M a ) is shown on the right of  Identified parameters are f s * and f A PaO 2,meas and PaCO 2,meas are the measured blood gases obtained at a specific condition described by α. They are used to determine f s * and f A * that best reproduce the measurements in the forward model: Parameter identification was performed by minimizing the sum of the squared error (SSE) between measured (meas) and predicted (pred) partial blood gas pressures in arterial blood: The weighting factor of 3 for PaCO 2 was chosen to avoid imbalanced influence of PaO 2 data on the identification process, as dimension of PaCO 2 is approximately three times smaller than PaO 2 . Minimization of the above described objective function was carried out using fminsearchbnd in MATLAB (R2012a, The Mathworks, Natick, MA, USA). fminsearchbnd is distributed under the BSD license and is based on fminsearch, the MATLAB function that employs the Nelder-Mead simplex search method [20]. According to [21] a shunt of 50 % and above leads to increases in FiO 2 having no effect on PaO 2 . Additionally, f A values above 0.9 lead to a swap of the high V /Q-compartment with the low V /Q-compartment, essentially mirroring V /Q-values of f A below 0.9. Thus parameter constraints for {f s , f A } were set as nonnegative lower boundaries LB = {0, 0} and as upper boundaries UB = {0.5, 0.9}.
Structural identifiability of the model using multiple measuring points was shown in a previous report [22]. Initial fs was set to 0.2, which showed to lead to the global minimum of the objective function in all test cases. Initial f A was set arbitrarily to 0.5 to start the minimization process at a certain initial mismatch between ventilation and perfusion. Constant patient state for the time of model prediction is assumed.

Model prediction
Prediction of blood gas levels depending on FiO 2 is done using forward calculation of

Data
We have employed both simulated and recorded real patient data to evaluate how well the described model is identifiable with data obtained at one single FiO 2 level.
Simulated data The same two-parameter model of pulmonary gas exchange was used to create experimental data. It allows calculating the impact of noise in the data because the correct results for parameter identification are known a priori. Twelve classes of patient data sets have been generated, that differ in the parameter combinations of f s and f A . Those have been chosen to represent different stages of pulmonary disease. Model parameters used for data generation, resulting FiO 2 /PaO 2 -ratios and the classifications of pulmonary impairment [24] are listed in Table 1. BGA and physiological standard .
values of an adult man were used for data generation. These physical constants are listed and explained in Table 2.
More formally, we define twelve patient classes by the model parameters For each of the twelve patient classes, 1000 simulated measurements equidistant between FiO 2 of 21 % and 100 % were determined. Depending on FiO 2 settings, model simulation led to Measuring PaO 2 and PaCO 2 via blood gas samples drawn from the arterial line is the current gold standard in clinical practice [25,26], while measuring arterial oxygen saturation via pulse oximetry is accurate within ±2 % of the true value [27]. Thus, to   values ranging from 21 to 100 % in steps of 5 %: Real patient data: Two real patient data sets were used for the plausibility check of the results obtained from the theoretical analysis. Real patient data including at least four blood gas measurements at different FiO 2 settings in mandatory ventilation mode were retrieved from a patient data management system of the university medical centre in Kiel [7]. Two data sets, a mild (Pat R1) and a critically ill patient (Pat R2), met those demands. The recorded levels of FiO 2 were applied on a therapeutical basis, thus not systematically in the context of a clinical trial. The data sets included measurements of PaO 2 , PaCO 2 , f R , Q, V tid , V ds , cHb, pH, T and FetCO 2 at each of the applied FiO 2 levels. Patient data did not include FetO 2 measurements, thus FetO 2 was approximated from: Here, the respiratory quotient RQ was assumed to be 0.8. As with the simulated data sets, initial conditions of {fs, fA} for parameter identification were set to {0.2, 0.5}.

Analysis of structural identifiability
To verify structural identifiability of a model system, uniqueness of the solution of parameter identification has to be proven. The simplicity of the two-parameter model allows a numerical calculation and two-dimensional visualization of the objective function. The SSE is calculated and plotted for different parameter combinations to visualize the contour of the error surface. A single global minimum of the objective function indicates structural identifiability of the model.
Structural identifiability of the model using one measurement point at one FiO 2 level was analyzed with the synthetic as well as the two real patient data sets. The error surfaces (SSE) were plotted as a function of model parameters f s and f A with a resolution of 90 × 90.

Evaluation of quality of fit
Besides quantity, quality of measurements used for model identification is essential for the accuracy of parameter identification. To verify practical identifiability, the influence of measuring errors on identification behavior of the model system was evaluated using the 1000 virtual measurements (PaO

Verification of results with patient data
Two real patient data sets were used to confirm the findings obtained with simulated data. Identification was conducted at each FiO 2 value that was recorded in the particular patient. Predictive performance was evaluated by comparing measured and predicted PaO 2 and PaCO 2 at all four recorded FiO 2 levels:

Results
Visualizing the objective function Figure 3 visualizes the contour of the objective function evaluated at one single measuring point. Figure 3a shows synthetic patient data, while Fig. 3b is devoted to real patient data. The contour lines (SSE) are scaled logarithmically to improve visibility of the minimum. For all analyzed data sets, a single global minimum in the error surface as a function of the model parameters f s and f A could be detected. The parameter combination leading to the global minimum was in agreement with the parameters used for data generation. Within accuracy of numerical representation (double), the SSE value was zero at PaCO ts 2,pred i − PaCO ts 2,sim i .

Prediction of PaCO 2 and PaO 2 depending on FiO 2
Over all tested data sets, mean PaCO 2 was 2.4 % (1.3 mmHg) of the true value with a standard deviation of ±1.6 % (±0.8 mmHg). Figure 4a shows mean deviation of PaO 2 with respect to α. PaO 2 prediction was less accurate when identification data were recorded at low FiO 2 levels, especially for data sets representing pulmonary impairment. For data set 1 (healthy lung) minimal PaO 2 was achieved for 40 % < FiO 2 < 50 % whereas model identification with data set 12 (severely impaired lung) shows the smallest PaO 2 for α = (100 %, Ψ). In the minima, the model was able to reproduce PaO 2 of all of the simulated patient data sets with a mean deviation of less than 2.5 % (<2.5 mmHg) of the true value with a standard deviation of less than 3 % (<2 mmHg). In Fig. 4a, the minima are marked with vertical lines and the respective patient numbers.  Fig. 4 Left clustered mean deviation of PaO 2 over FiO 2 for simulated data. Deviation of prediction of PaO 2 depends on FiO 2 range used for model identification. Broken lines show respective minima for the different data sets (numbered). Minimum is located at a higher FiO 2 range for data representing a higher pulmonary distress. Right mean of ΔPaO 2 over FiO 2 for real patient data. Deviation of prediction of PaO 2 varies with the FiO 2 level used for model identification. The location of the minimum depends on patients' pulmonary state. Mean deviation of ΔPaO 2 was less than 10 % at the minimum for both data sets Mean deviation of ΔPaO 2 for real patient data sets is shown in Fig. 4b. Predictions with model parameters being identified at low and high FiO 2 show higher deviations from measured values than for identification at medium FiO 2 levels. Mean deviation of ΔPaO 2 was less than 10 % (8 % or 5.8 mmHg for Patient R1 and 6 % or 4.6 mmHg for Patient R2) at the minimum. The best performance of PaO 2 prediction was found for α = (40 %, Ψ) (Pat R1) and α = (80 %, Ψ) (Pat R2) respectively. Figure 5 summarizes the FiO 2 levels leading to a minimum of mean deviation of PaO 2 for simulated data sets (colored markers). 32 additional simulated data sets (black markers) were generated to illustrate the relation between optimal FiO 2 and pulmonary impairment more precisely. Location of the minimum was shifted to a higher FiO 2 cluster with increasing pulmonary impairment, i.e. decreasing PaO 2 /FiO 2 -ratio.

Discussion
Using mathematical models for decision support in clinical practice requires high level of safety and accuracy of the model predictions. Furthermore, measuring effort for model identification, i.e. identification of patient specific parameters, should be kept to a minimum.
Identification of the two-parameter model of human pulmonary gas exchange requires data from arterial blood gas analysis and photoplethysmographic saturation measurement. The model has previously been applied with four measurements at different levels of FiO 2 . To reduce the effort required in clinical practice to identify gas exchange models it was investigated if a reduced number of FiO 2 levels for data acquisition is sufficient.
Structural identifiability of the model applying a single identification data point was examined in this study with both simulated and real patient data sets. One single global minimum of the objective function is an indication for structural identifiability of the model. Without noise in the data, one single blood gas measurement is sufficient for robust model identification.  Table 1 (numbered). Additional data were generated to visualize the curve progression and to confirm the findings (black markers). Minimum of PaO 2 mean deviation is shifted to a higher FiO 2 cluster for increasing severity of pulmonary impairment (lower PaO 2 /FiO 2 -ratio) The effect of a change in FiO 2 on the concentration of carbon dioxide in arterial blood is negligible. PaCO 2 data could be reproduced by the model with high accuracy. Prediction error of CO 2 data was below the noise level of ±5 % for all data sets. However, measuring errors may decrease the accuracy of parameter identification and therefore model prediction of PaO 2 , especially if identification is based on only one measurement. Results show that the gas exchange model with shunt and V /Q-mismatch is able to fit both the synthetic and the real patient data with good accuracy, as already presented in former work [16]. Oxygenation data of all data sets could be reproduced by the model with a mean deviation below 10 % in spite of measuring errors in the identification data. However, when identifying the model with noisy data, the FiO 2 setting influences the prediction accuracy of PaO 2 . This influence was therefore examined to find a guideline how to choose an appropriate FiO 2 level for data acquisition in the identification process.
The identification processes of both simulated and real patient data sets representing a variety of different disease states showed surprisingly similar results. It could be pointed out that accuracy of model prediction of blood gas concentration is related to the FiO 2 setting when recording identification data. The optimum FiO 2 level depends on the level of pulmonary impairment whereupon FiO 2 should be increased with increasing severity of pulmonary impairment. In severely ill patients, oxygenation of arterial blood is inhibited, thus higher FiO 2 levels are required in order to achieve adequate PaO 2 levels. Figure 6a shows the mean deviation of PaO 2 with respect to PaO j 2,meas . It could be observed that a minimal prediction error is achieved for PaO j 2,meas in the range of 150-200 mmHg for the entire simulated data sets. Identification at higher levels leads to a small increase of mean deviation. However, mean deviation of PaO 2 was found to be still below 5 % for all tested patient cases. Identifying the model with PaO 2 levels of less than 100 mmHg is potentially leading to an increase in both mean and standard deviation of model prediction. Because of the high severity of pulmonary impairment, data sets 7 and 8 do not achieve a PaO 2 of 100 mmHg, even when a FiO 2 of 100 % is applied. Patients with such high pulmonary impairment have to be ventilated using the highest FiO 2 possible to achieve a sufficient oxygenation. Fig. 6 Left clustered mean deviation of PaO 2 over PaO j 2,meas for simulated data. Minimum in prediction error of PaO 2 data is in the range of 150-200 mmHg for data sets 1-7. Because of the high severity of pulmonary impairment, data sets 8-12 do not achieve this oxygenation range, even when a FiO 2 of 100 % is applied. Right mean of ΔPaO 2 over PaO real 2,meas . Real patient data tested in our study confirm the curve progression of the study with the simulated data. The best prediction performance is shown for identification in the PaO 2 range of 70-80 mmHg Mean deviation of ΔPaO 2 over PaO real 2,meas is shown in Fig. 6b. Both curves representing real data confirm the results of the analysis with simulated data, but minimum mean deviations were found for a PaO 2,meas of 73 mmHg and 81 mmHg respectively.
When only one measuring point is used for model identification, success of parameter identification obviously depends on the level of PaO 2 at the particular time of the measurement. Low PaO 2 levels bear the risk of large influence of measuring errors in identification data. In low PaO 2 levels, even small changes in PaO 2 used for identification may lead to an overestimation of shunt fraction f s and therefore to a smaller increase of PaO 2 for higher FiO 2 levels compared to correct data. In severely ill patients, this effect is more prominent than in patients with less pulmonary distress.
PaO 2 levels suggested to be optimal for identification of the model may be not achievable in patients with severe pulmonary impairment. Here, identification at an FiO 2 of 100 % was shown to achieve predictions with highest accuracy.
In the presented work, we have investigated a FiO 2 range 21-100 %. To separate effects of shunt from those caused by low V /Q-ratio, subatmospheric oxygen levels should be considered as well. However, the intended use of the applied model and the presented routine of identifying the model parameters with only one blood gas measurement are in a critical care environment where such oxygen levels are not applied.
Karbing et al. [16] have previously presented a three-parameter extension of the model used in this work which uses an adjustable distribution of non-shunted blood among the alveolar compartments. This model shows to be superior in terms of reproducing PaCO 2 especially in patients with V /Q-ratios below 0.27. The presented work however focuses mainly on finding the optimal calibration point of FiO 2 , which has only a minor effect on PaCO 2 . Nevertheless, investigating the structural identifiability and the possibility of using only one measurement set to also identify the three-parameter model should be considered in future work.
The model of gas exchange applied in the presented study includes assumptions such as steady state conditions of blood gases and constant alveolar ventilation and perfusion. Thus, tidal breathing as is the reality in humans is not considered. Those assumptions present shortcomings compared to the reality those models try to reproduce [28]. Several models including tidal breathing have been presented in the past [28][29][30] however continuous measurement of blood gases in combination with continuous measurements of gases in inspired and expired air is not routinely available at the bedside at this moment. Thus clinical applications are currently limited to models assuming continuous ventilation and perfusion as well as equilibrated blood gases. Still, Karbing et al. [16] have shown previously that the applied model is well capable of reproducing patient data for a wide range of lung impairments and that the model can thus be used in a clinical environment as a prediction tool.
Model based decision support in clinical practice implies that the mathematical model is identifiable with a minimum of measuring effort. We could show that the two-parameter gas exchange model with shunt and V /Q-mismatch is structural identifiable with only one single blood gas measurement. Using only one single measurement, possible measuring errors are not averaged. However, simulation results show that model based prediction of blood gases for different FiO 2 is possible with a mean prediction error below 10 % for a maximum measuring error of 5 %. Simultaneously, we could determine