Evaluating different approaches to identify a three parameter gas exchange model

1 Introduction

Mathematical models have become a valid tool to simulate physiological reactions of patients treated in intensive care units [1]. Adapted to the individual patient properties, they might be used to predict the outcome of changes in the therapy settings. For example, models of gas exchange or respiratory mechanics can be fit to mechanically ventilated patients. The fitted models can then be used to predict the effect of changes in ventilation settings such as breath rate, PEEP or FiO₂. Using those predictions in a decision support system might enable a clinician to calculate optimized settings to fulfil given therapy goals [2], [3], [4].

In order to be applicable at the bedside, such models need to be identifiable both structurally and practically. In order to identify a unique model parameter set that represents the recorded patient data, the parameter set needs to be obtainable from the data available at the bedside. While structural identifiability may be proven through mathematic tools, practical identifiability has to be tested specifically for the intended application. Practical identifiability considers the data available at the bedside and the expected range of patient physiologies [5].

Nonlinear mathematical models that are not identifiable through linear regression need to be identified using gradient based methods. The performance of those methods strongly depends on the selection of appropriate initial estimates. Thus, using a direct approach with arbitrary initial estimates might lead to incorrect parameter values leading to poor predictions. A previously published approach using easily identifiable models of lower order to obtain good initial guesses for the identification of models of higher order has shown satisfying results in models of respiratory mechanics [6]. The approach exploits model hierarchies, in particular, a family of models that are related to each other. In those hierarchies, models of higher order are extensions or modifications of models of lower order in the family.

The presented work aims at incorporating an identification approach that is related to the hierarchical approach to identify a three-parameter model of gas exchange. The hierarchical approach is compared with the standard approach of directly identifying the model with arbitrary initial guesses and a global search algorithm.

2 Methods

2.1 Model

The model of gas exchange is based on the three parameter model presented by Karbing et al. [7]. The model consists of two alveolar compartments which simulate the air distribution among differently ventilated parts of the lung. Additionally, it comprises a shunt to simulate a part of the venous blood not being oxygenated but being mixed directly with the oxygenated arterial blood. The non-shunted blood is distributed differently among the two alveolar compartments to model various ventilation (V˙) to perfusion (Q) ratios (V˙/Q). The three model parameters defining model behaviour thus are the fraction of shunted blood (fs), the fraction of inhaled air into one of the two alveolar compartments (fA) and the fraction of blood being distributed to that compartment (fQ). End-tidal gas fractions are thus defined as:

(1)F⁢e⁢tx=(1−f⁢A)⋅F⁢Ax,1+f⁢A⋅F⁢Ax,2

Index x denotes O₂ and CO₂ here. FA_x are the alveolar gas fractions. Venous concentrations are derived from:

(2)C⁢vx,1=C⁢cx,1−V˙x,1/(Q⋅(1−f⁢s)⋅(1−f⁢Q))

(3)C⁢vx,2=C⁢cx,2−V˙x,2/(Q⋅(1−f⁢s)⋅f⁢Q)

Capillary gas concentrations Cc_x are calculated from FA_x using the gas dissociation equations [8], [9]. Q denotes blood flow, V˙_x are oxygen consumption and CO₂ production. They are defined as:

(4)V˙x,1=(1−f⁢A)⋅V˙A⋅(F⁢ix−F⁢Ax,1)

(5)V˙x,2=f⁢A⋅V˙A⋅(F⁢ix−F⁢Ax,2)

Fi_x are the inspired gas fractions. Arterial gas concentrations are then calculated from:

(6)C⁢ax=C⁢cx,1⋅(1−f⁢s)⋅(1−f⁢Q)+C⁢cx,2⋅(1−f⁢s)⋅f⁢Q+C⁢vx⋅f⁢s

Parameter fs has a range of 0–0.5, while fA and fQ can range between 0.1 and 0.9. Model inputs are inspired oxygen fraction, air flow and Fet_x. Figure 1 shows a schematic representation of the model.

Figure 1

Three parameter gas exchange model.

2.3 Identification algorithms

Three different identification approaches were tested: direct approach, global search, and hierarchical approach. The objective function W to be minimized by all tested approaches was:

(7)W=∑(P⁢m⁢O2−P⁢a⁢O2)2+(P⁢m⁢C⁢O2−P⁢a⁢C⁢O2)2

Direct approach: The direct approach starts at an arbitrary initial guess and uses a gradient based method to find the global minimum of the objective function. In the presented work, the Nelder-Mead Simplex-Search method [10] was used that is incorporated in MATLAB (R2015a, The Mathworks, Natick, MA, USA) as fminsearch. To evaluate the influence of the selected initial estimate, 100 constellations of initial estimates were randomly generated within a defined range. The number of successful identifications, i.e. if the identified parameter values were within 1% of the true value, was counted. These evaluations were repeated for a total of 10 extending ranges around the original values.

Global search: A global search approach starts at multiple initial estimates (trial points) and evaluates the convergence and the final value of the objective function for each of those estimates. The trial points that lead to the best scores (in terms of the objective function and violation of given constraints) are then used to create new trial points through deterministic combination [11]. The algorithm terminates after a given number of trial points have been evaluated. In the presented work we used the global search algorithm implemented in MATLAB.

Hierarchical identification: Hierarchical identification uses identification results of models of lower order to calculate good initial estimates for the identification of a model of higher order. The models used need to be related to each other in terms of a hierarchy. We have previously created a hierarchical family of gas exchange models shown in Figure 1 [12]. Specifically, a simple shunt model [3] is used to calculate the shunt from a blood gas measurement. The calculated shunt (fs) is then used as an initial estimate in the identification of a two parameter model that is derived from the three parameter model. Here, fQ is fixed to a certain value (fQi = 0.1, 0.2, 0.3,… 0.9). Parameters fs and fA of the two parameter model are then identified for the specific fQi. The calculated fs from the simple shunt model along with the combination of fA and fQi where the objective function resulted in the lowest value are finally used as initial estimates for the identification of the three parameter model. Figure 2 shows an overview of the described algorithm. The identification of the two and the three parameter model was done with fminsearch.

Figure 2

The hierarchical identification approach. The simple shunt model is used to calculate a valid initial guess (fs0) for the identification of the two parameter model. There, fA is identified for nine different fQi. The best combination of fA and fQi along with fs0 is then used as initial values for the identification of the three parameter model.

Global search and the hierarchical approach do not require initial estimates, thus in each patient only one identification was done. To allow a fair comparison, the maximum number of function evaluations for the direct approach was set to 10,000, while each of the hierarchical steps (nine iterations with fixed fQi and one final identification of the three parameter model) was allowed to call the function 1000 times. The global search was limited to a computing time of 500 s per run. Computing time for the tested approaches was recorded on an i7-4770 CPU (4x3.4GHz) with 12GB RAM.

3 Results

Figure 3 shows the averaged hit rate over all 12 virtual patients for the three parameter identification approaches. The hit rate was defined as the number of identifications that resulted in values within 1% of the true value (X_true) for all three parameters, with respect to the range (n) from which the initial estimates (X) were randomly drawn (X_true ± n). Hit rates of parameter identification with direct approach showed a decrease from 68% to 47% when the n-Range is increased from 0.05 to 0.5. The averaged deviation between the identified parameter values and the true values was 2.5% at an n-Range of 0.05 and increased to 17.6% when drawing initial estimates from X_true ± 0.5. Identification with global search showed a hit rate of 0% in all of the tested patients. In contrast, the hierarchical approach found the correct values (hit rate 100%) in all 12 patients.

Figure 3

Average hit rates for all 12 patients. Hit rate denotes the number of identifications that resulted in parameter values within 1% of the true value for all three model parameters. Hit rates are shown with respect to the range from which the initial estimates X were randomly drawn.

Averaged computing time for the direct approach was 39.4 s, 476.7 s for the global search approach, and 191.9 s for the hierarchical approach.

4 Discussion

Mathematical models can help in gaining a better understanding of the human physiology. When they are adapted to an individual human using various measurements, they are able to reproduce physiological reactions to changes in therapeutic settings. Thus they may be used to predict those reactions and can therefore be implemented in a decision support system that aims to optimize therapeutic outcomes to achieve goals defined by a clinician. To be applicable in a real clinical setting, those models need to be robustly identifiable from the data available at the bedside. Not only do incorrect parameter values influence the prediction of the reaction to changes in the therapy but they may also yield an incorrect picture of the current disease state of a patient.

The results show that when using a direct approach to identify the non-linear three-parameter gas exchange model, the outcome is highly dependent on the quality of the initial estimates. In particular, the direct approach had a higher hit rate when the initial values were close to the real values (Figure 3). Greater initial value ranges yielded lower hit rates for the direct approach parameter identification.

The proposed hierarchical approach eliminates the influence of the initial estimates by using models related to the three parameter models but of less complexity to compute adequate initial estimates. It was therefore able to determine better initial value estimations for the higher order three-parameter model. The hierarchical approach was thus able to find the true parameter values in all patients.

The evaluation of global search showed that this approach was not applicable to the presented problem, at least not using the applied settings, i.e. the maximum allowed computing time. The longer the computing time in a global search approach is set, the more trial points can be tested, thus with an infinite number of trial points and an unlimited computing time, the global search will should be able to find the correct parameter values. However, in a clinical setting parameter identification must be possible on a regular basis to allow the model to be adapted to changes in the patient’s physiology. Thus computing time needs to be low. The results show that the global search was allowed twice the computational time needed by the hierarchical approach and still yielded a higher deviation between identified and true values.

Using artificial patient data to test algorithms is a valid tool during development and initial evaluations. Since the true parameter values are known, they can be easily compared to the values identified by the evaluated algorithms. However to test the clinical applicability, real patient data must be employed for evaluation. Tests with data drawn from real patients are therefore planned for the future.