Iterative integral parameter identification of a respiratory mechanics model

Abstract

Background

Patient-specific respiratory mechanics models can support the evaluation of optimal lung protective ventilator settings during ventilation therapy. Clinical application requires that the individual’s model parameter values must be identified with information available at the bedside. Multiple linear regression or gradient-based parameter identification methods are highly sensitive to noise and initial parameter estimates. Thus, they are difficult to apply at the bedside to support therapeutic decisions.

Methods

An iterative integral parameter identification method is applied to a second order respiratory mechanics model. The method is compared to the commonly used regression methods and error-mapping approaches using simulated and clinical data. The clinical potential of the method was evaluated on data from 13 Acute Respiratory Distress Syndrome (ARDS) patients.

Results

The iterative integral method converged to error minima 350 times faster than the Simplex Search Method using simulation data sets and 50 times faster using clinical data sets. Established regression methods reported erroneous results due to sensitivity to noise. In contrast, the iterative integral method was effective independent of initial parameter estimations, and converged successfully in each case tested.

Conclusion

These investigations reveal that the iterative integral method is beneficial with respect to computing time, operator independence and robustness, and thus applicable at the bedside for this clinical application.

Introduction

Mechanical ventilation (MV) therapy in intensive care units (ICU) carries the risk of severe additional complications for the patient due to non-optimal ventilator settings [1]. To reduce patient risk, optimal patient-specific ventilator settings must be found. Mathematical models of respiratory mechanics effectively predict the outcome of specific ventilator settings and thus support clinical evaluation of lung protective settings [2]. Prediction quality depends on how well model parameters correspond to the true patient properties and on the accuracy of the model representing the individual’s lung physiology. To assure clinical applicability at the bedside, especially in closed loop settings [2, 3], reliable and immediate parameter identification is required. However, information for analyzing respiratory mechanics directly at the bedside is generally restricted to measurements of airway pressure and flow rate. Therefore, any model must enable easy and fast parameterization in terms of the available clinical data. Thus, models should capture all necessary dynamics and be as simple as possible. Hence, lumped parameter models are a common approach for optimizing ventilator settings [4].

Typical parameter identification methods minimize the least square error (LSE) between measured samples and model simulations. For single compartment models of respiratory mechanics, multiple linear regression is an established straightforward approach [5, 6]. Multiple linear regression is also applicable to higher order linear models, but the quality of results can be distorted by noise [6]. Alternative parameter identification methods for higher-order and nonlinear models are iterative error-mapping methods [7–11], such as the Simplex-Search Method [12] or Levenberg-Marquardt Algorithm [13, 14]. These methods require initial guesses for the model parameters and iterate towards LSE values by approaching the minimum on the error surface. However, with an increasing number of degrees of freedom, several parameter constellations can appear as additional possible solutions (local minima) [8]. Hence, such methods can report non-optimal parameter values and are highly dependent on initial parameter estimates. To reduce the probability of convergence to local minima, various global search strategies, such as random search, simulated annealing or genetic algorithms can be used [15]. However, these methods are time consuming and the quality of the solution is sensitive to the parameterization of the algorithm [16].

This research presents an iterative integral method (IIM) [17, 18] for respiratory mechanics parameter identification, and compares the performance against established parameter identification methods. The method was originally developed for parameter identification of glucose-insulin models [19]. In the context of the glucose-insulin models, the IIM is comparatively simple to apply, requires comparatively minimal computing time and does not require the estimation of initial values [17, 19]. The IIM is tested with simulation data to demonstrate its robustness and efficiency, and clinical study data is used to demonstrate the clinical potential.

Model & methodology

Viscoelastic model

Viscoelastic Models (VEM) of respiratory mechanics are established models and assume that the tissues comprising the walls of the alveolar compartment are viscoelastic, rather than simply elastic [20]. The VEM used in this research is an established model of respiratory mechanics whose applicability with respect to healthy and lungs with Acute Respiratory Distress Syndrome (ARDS) has been demonstrated [21, 22]. Ganzert et al. [10] showed the different trends of the VEM parameters in healthy and ARDS subjects, which might be used as indicator for ARDS development.

The analogous electrical circuit for the VEM is shown in Figure 1 and the mathematical description is presented in state-space representation in Eq. (1):

$$\begin{array}{l}\left[\begin{array}{c}\hfill {\dot{p}}_{C1}\hfill \\ \hfill {\dot{p}}_{C2}\hfill \end{array}\right]=\left[\begin{array}{cc}\hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1/\left(R_2{C}_2\right)\hfill \end{array}\right]\left[\begin{array}{c}\hfill p_{C1}\hfill \\ \hfill p_{C2}\hfill \end{array}\right]+\left[\begin{array}{c}\hfill 1/C_1\hfill \\ \hfill 1/C_2\hfill \end{array}\right]\dot{V}\hfill \\ p_{\mathit{aw}}=\left[\begin{array}{cc}\hfill 1\hfill & \hfill 1\hfill \end{array}\right]\left[\begin{array}{c}\hfill p_{C1}\hfill \\ \hfill p_{C2}\hfill \end{array}\right]+R_1\dot{V}\hfill \end{array}$$

(1)

Figure 1

Viscoelastic Model of Respiratory Mechanics. R ₁ denotes the airway resistances and C ₁ the static compliance of the respiratory system. R ₂ and C ₂ are the resistance and the compliance of the viscoelastic component. The respiratory airflow $\dot{V}$ represents the input and the airway pressure p _aw the output of the model.

The flow rate $\dot{V}$ [mL/s] defines the model input and the airway pressure p _aw [cmH₂O] defines the model output. Airway resistance is described by R ₁ [cmH₂O·s/mL] and the static lung compliance by C ₁ [mL/cmH₂O]. R ₂ [cmH₂O·s/mL] and C ₂ [mL/cmH₂O] account for the viscoelasticity described by the viscoelastic time constant R ₂ ·C ₂ . p _C1 and p _C2 [cmH₂O] are the internal state variables being the pressure in the compliant static and viscoelastic lung compartment.

The four parameters (R ₁, C ₁ , R ₂, C ₂) correspond to the unknown patient-specific parameters:

$$X_{\mathit{VEM}}:=\left\{R_1,C_1,R_2,C_2\right\}$$

(2)

A fundamental prerequisite for successful parameter identification is a-priori structural identifiability of the parametric model [23]. This necessary criterion states that under ideal conditions of noise-free observations and error-free model structure, the unknown parameters of the model can be uniquely recovered from the measured input–output variables. The VEM is globally identifiable [8] using the DAISY (Differential Algebra for Identifiability Systems) computer program to automatically check for a-priori structural identifiability [23].

Data

First, the parameter identification methods were compared in-silico using synthetic data sets generated by simulating parameterized VEMs. Second, the methods were tested using clinical data to assess their clinical applicability.

Simulation data

A flow profile corresponding to the Volume-Controlled-Ventilation (VCV) mode was applied to parameterized VEMs. Inflation was simulated by a constant flow of 500 mL/s over 1 s, followed by an end-inspiratory pause (EIP) of 4 s corresponding to the closure of the expiratory valve as observed in clinical data [24].

100 different sets of model parameters were generated for a Monte-Carlo Analysis. The model parameters were randomly selected from an evenly distributed range between 0.5 and 1.5 times a realistic physiological parameter set (R ₁  = 0.010 cmH₂O·s/mL, C ₁  = 30.00 mL/cmH₂O, R ₂  = 0.020 cmH₂O·s/mL, C ₂  = 80.00 mL/cmH₂O).

After simulating each VEM, 10 different sets of white noise (5% of given pressure sample) were added in order to impose noisy measurement error. An example of the simulated pressure response is shown in Figure 2. It consists of an increase during inspiration and an exponential drop during the EIP.

Figure 2

Predefined flow profile and simulated pressure response of the VEM. The simulation consists of an inspiration part (1 s), followed by an end-inspiratory pause (4 s).

Clinical data

Measurement sets from N = 13 mechanically ventilated patients were selected from a previous ARDS– Study [10, 25], where SCASS-Maneuvers (Static Compliance Automated Single Step) were performed. The SCASS-Maneuver consists of airway occlusions within the inspiration phase of a breathing cycle that are initiated when a randomized inspiration volume is reached. The measurement set consisted of flow rate ($\dot{V}$) and airway opening pressure (p _aw ) signals sampled at 125 Hz. For each patient, two breathing cycles were extracted. During occlusion, the airway pressure reached a quasi-static equilibrium exponentially.

The study was approved by the local ethics committees of the participating university hospitals. Informed consent was obtained from patients or their legally authorized representative. The detailed description of the experimental setup is shown in [25].

Parameter identification

Five parameter identification methods were used to identify the patient-specific VEM parameters Eq. (2) using the simulation and clinical data. These methods also provide references for the iterative integral method. To compare computational efficiency, all methods were tested on the same Desktop PC (Intel Core 2 Duo, 2.80 GHz), where computation time was measured.

Multiple linear regression (MLR)

Given a linear model, MLR is an obvious approach for parameter identification [5, 6]. MLR requires the input–output relation of the model, eliminating the state variables. This relation is derived from the state space description of the VEM by calculating the derivative of the output equation from Eq. (1) and inserting the definition of the state variables:

$${\dot{p}}_{\mathit{aw}}={\dot{p}}_{C1}+{\dot{p}}_{C2}+R_1\ddot{V}=\frac{1}{C_1}\dot{V}-\frac{1}{R_2{C}_2}{p}_{C2}+\frac{1}{C_2}\dot{V}+R_1\ddot{V}$$

(3)

To eliminate p _C2 in Eq. (3), the output of Eq. (1) is rewritten in terms of p _C2 :

$$p_{C2}=p_{\mathit{aw}}-p_{C1}-R_1\dot{V}=p_{\mathit{aw}}-\frac{1}{C_1}V-R_1\dot{V}$$

(4)

Eq. (4) is inserted in Eq. (3) leading to the input–output relation of the model:

$${\dot{p}}_{\mathit{aw}}=-\frac{1}{R_2{C}_2}{p}_{\mathit{aw}}+R_1\ddot{V}+\left(\frac{1}{C_1}+\frac{1}{C_2}+\frac{R_1}{R_2{C}_2}\right)\dot{V}+\frac{1}{R_2{C}_1{C}_2}V$$

(5)

Effectively, Eq. (6) re-arranges the state-space equation Eq. (1) to place the model variables as functions of measured values (p _aw and $\dot{V}$) and their derivatives or integrals.

$$p_{\mathit{aw}}=-R_2{C}_2{\dot{p}}_{\mathit{aw}}+R_1{R}_2{C}_2\ddot{V}+\left(\frac{R_2{C}_2+R_2{C}_1+R_1{C}_1}{C_1}\right)\dot{V}+\frac{1}{C_1}V$$

(6)

MLR requires a sum of independent variables scaled by a multiplicative factor. Therefore the coefficients of Eq. (6) can be represented by new variables:

$$p_{\mathit{aw}}=A{\dot{p}}_{\mathit{aw}}+B\ddot{V}+C\ddot{V}+DV$$

(7)

For parameter identification according to the Least-Square-Error (LSE) principle, Eq. (7) was arranged in matrix form:

$$\begin{array}{ll}\left[\begin{array}{cccc}\hfill {\dot{p}}_{\mathit{aw}}\left(1\right)\hfill & \hfill \ddot{V}\left(1\right)\hfill & \hfill \dot{V}\left(1\right)\hfill & \hfill V\left(1\right)\hfill \\ \hfill {\dot{p}}_{\mathit{aw}}\left(2\right)\hfill & \hfill \ddot{V}\left(2\right)\hfill & \hfill \dot{V}\left(2\right)\hfill & \hfill V\left(2\right)\hfill \\ \hfill {\dot{p}}_{\mathit{aw}}\left(3\right)\hfill & \hfill \ddot{V}\left(N\right)\hfill & \hfill \dot{V}\left(N\right)\hfill & \hfill V\left(N\right)\hfill \end{array}\right]\hfill & \left[\begin{array}{c}\hfill A\hfill \\ \hfill B\hfill \\ \hfill C\hfill \\ \hfill D\hfill \end{array}\right]\hfill \end{array}=\left[\begin{array}{c}\hfill p_{\mathit{aw}}\left(1\right)\hfill \\ \hfill p_{\mathit{aw}}\left(2\right)\hfill \\ \hfill p_{\mathit{aw}}\left(N\right)\hfill \end{array}\right]$$

(8)

Values of A-D are identified applying the Moore-Penrose pseudo-division for over-defined matrix equations. Patient-specific R ₁ , C ₁ , R ₂ , C ₂ values are then uniquely derived from the identified values of A-D.

Integral Method (IM)

The IM is similar to MLR in terms of identification steps. However, integrals are used to improve robustness to noise [19]. Hence, Eq. (5) was integrated assuming p _aw (0) = 0:

$$p_{\mathit{aw}}=-\frac{1}{R_2{C}_2}\int p_{\mathit{aw}}dt+R_1\dot{V}+\left(\frac{1}{C_1}+\frac{1}{C_2}+\frac{R_1}{R_2{C}_2}\right)V+\frac{1}{R_2{C}_1{C}_2}\int Vdt$$

(9)

The coefficients of Eq. (9) are represented by new variables:

$$A=-\frac{1}{R_2{C}_2}$$

(10a)

$$B=R_1$$

(10b)

$$C=\frac{1}{C_1}+\frac{1}{C_2}+\frac{R_1}{R_2{C}_2}$$

(10c)

$$D=\frac{1}{R_2{C}_1{C}_2}$$

(10d)

Yielding:

$$p_{\mathit{aw}}=A\int p_{\mathit{aw}}dt+B\dot{V}+CV+D\int Vdt$$

(11)

Incorporating Eq. (11) into an over-defined matrix system yields:

$$\left[\begin{array}{cccc}\hfill {\int }_0^1{p}_{\mathit{aw}}dt\hfill & \hfill \dot{V}\left(1\right)\hfill & \hfill V\left(1\right)\hfill & \hfill {\int }_0^{1}Vdt\hfill \\ \hfill {\int }_0^2{p}_{\mathit{aw}}dt\hfill & \hfill \dot{V}\left(2\right)\hfill & \hfill V\left(2\right)\hfill & \hfill {\int }_0^{2}Vdt\hfill \\ \hfill {\int }_0^N{p}_{\mathit{aw}}dt\hfill & \hfill \dot{V}\left(N\right)\hfill & \hfill V\left(N\right)\hfill & \hfill {\int }_0^{N}Vdt\hfill \end{array}\right]\left[\begin{array}{c}\hfill A\hfill \\ \hfill B\hfill \\ \hfill C\hfill \\ \hfill D\hfill \end{array}\right]=\left[\begin{array}{c}\hfill p_{\mathit{aw}},\left(1\right)\hfill \\ \hfill p_{\mathit{aw}}\left(2\right)\hfill \\ \hfill p_{\mathit{aw}}\left(N\right)\hfill \end{array}\right]$$

(12)

Patient-specific parameters R ₁ , C ₁ , R ₂ , C ₂ are uniquely regained using the identified values of A-D:

$$R_1=B$$

(13a)

$$R_2{C}_2=-\frac{1}{A}$$

(13b)

$$C_1=\frac{1}{D{R}_2{C}_2}$$

(13c)

$$C_2=\frac{R_2{C}_1{C}_2}{C{R}_2{C}_1{C}_2-R_2{C}_2-R_1{C}_1}$$

(13d)

Iterative Integral Method (IIM)

The IIM uses much the same method as the IM. However, the identified A-D values are used to re-simulate p _aw according to Eq. (11). Resimulated p _aw is then used to update the left-hand-side (LHS) integrals of Eq. (12). Subsequently, Eq. (12) is solved again yielding updated values of A-D. This process is repeated until convergence. The sum of squared error (SSE) between simulated and measured p _aw was calculated after every iteration. The relative termination tolerance (TolFUN) for changes in the function value (SSE) was set to 10^-4 This value is in accordance to the default setting of the proprietary gradient-based search methods in MATLAB.

Proprietary Levenberg-Marquardt Algorithm (LMA) and Simplex Search Method (SSM)

Gradient-Based Methods such as the Levenberg-Marquardt Algorithm (LMA) (MALTAB command: lsqnonlin) and the Simplex-Search Method (SSM) (MALTAB command: fminsearch) were applied as a references to the IIM using following default convergence criteria: terminate the optimization if the relative change in parameters is smaller than 10^-4 (MATLAB: TolX value) or when the relative change in the SSE value is smaller than 10^-4 (MATLAB: TolFUN value).

Gradient-Based Methods require initial estimates for the corresponding parameters. Initial estimates close to the global minimum increase the probability of successful parameter identification.

For the case of identifying simulated data, the parameter values that were used for simulation are known and were used as initial values to provide the fastest, easiest solution, and a conservative comparison.

For the case of clinical data, patient outcomes are initially unknown. Thus, convenient estimates for initial values are important. To decrease the potentially deleterious influence of poorly estimated initial values for gradient-based methods, a hierarchical method was developed by Schranz et al. [8]. This approach uses identification of simple respiratory models to provide prior knowledge. It is preferable for the simpler models to be identified with linear regression to avoid the influence of initial values. The identified parameter values are then used as initial values to support the identification of more complex models. In the VEM case, the 1^st order model (FOM) of respiratory mechanics (Figure 3) and Equation (14) was first identified by MLR leading to estimates of airway resistance (R) and respiratory compliance (C).

$$p_{\mathit{aw}}=R\dot{V}+\frac{1}{C}\int Vdt$$

(14)

Figure 3

1^st Order Model (FOM). R denotes the airway resistance and C static compliance of the respiratory system. The respiratory airflow $\dot{V}$ represents the input and the airway pressure p _aw the output of the model.

Additionally, the time constant τ of the exponential pressure drop during the zero-flow phase (Figure 2) was estimated by extracting the corresponding pressure interval fitted by an exponential function. Given an error-free model structure, this estimated time constant corresponds to the viscoelastic time constant R ₂ ·C ₂ . With the information of R, C and τ, initial values can be provided for gradient-based parameter identification of the VEM. This method improved the robustness of the subsequent parameter identification significantly by reducing the potential for poorly estimated initial values to result in local minima identification [8]. It is thus a more effective and conservative method comparison for the IIM.

Results

Parameter identification using simulation data

The identified parameter values from each parameter identification method were equal in the noise free synthetic data set and were effectively equal to the simulation values (R = 1.000). However, adding noise to the simulation data causes the true global minimum to slightly shift away from the original parameter values. In this case, the solutions reported by the SSM were used to define the new error minimum. The Monte-Carlo analysis revealed that the IIM, IM, SSM and LMA reported the same solution error minimum in every data set. The MLR didn’t locate a reasonable error minimum in any single case. The resulting parameter values of the MLR deviate by up to 100% from the values at the reported error minima. This outcome led to a median SSE 380% larger than the SSE of the error minima. The IIM always reached convergence after its first iteration and was thus equivalent to the IM in this case. The required median [interquartile] computing time of the IM and IIM was 0.8 [0.7-0.9] ms. In comparison, the gradient-based methods require 108.3 [101.6-115.7] ms for the SSM and 15.1 [15.0 – 15.5] ms for LMA.

Parameter identification using clinical data

All results using the clinical data are summarized in Table 1. Within all clinical data sets, the SSM converged to a minimal SSE with physiologically plausible parameter values. These values are considered as error minimum. The SSM required an average of 251 ms per data set. The LMA found the same minima in 17 data sets, but located local non-physiological minima with partly negative parameter values and higher SSE values in 9 out of 26 (35%) data sets. MLR identification resulted in non-physiological parameter values and significantly higher SSE values in every (100%) data set. IM identification reported non-physiological parameter values that were either negative or significantly higher than the error minimum in 8 of 26 (31%) data sets with higher SSE values than the SSM. Parameter identification using the IIM resulted in similar minima as the SSM in every data set.

Table 1 Median and minimal-, maximal values of the identified parameters with SSE and required computing time

Table 2 shows the IIM convergence in an example data set where the IM reported erroneous results. Starting from the erroneous IM result (0^th iteration), the SSE and non-physiological parameter values converge step by step towards the error minimum approaching physiologically plausible parameter values.

Table 2 Parameter convergence of the Iterative Integral Method (IIM) compared to the Simplex-Search Method (SSM)

The median required computing time in Table 1 for the IIM equals 5.02 ms reaching the minimum error faster than the SSM by a factor of 50. Example pressure responses for parameter sets identified by IIM, IM and MLR are shown together with the residuals in Figure 4. The residuals of the pressure response produced by IM parameter values are higher than the response of IIM parameter values. The pressure responses of the parameter sets identified by SSM and LMA are congruent with the solution of the IIM. Therefore these plots are neglected in Figure 4.

Figure 4

Simulation Results. Exemplary measurement set of pressure and flow with the simulated pressure responses and residuals of the VEM. VEM identified by the Iterative Integral Method (IIM), Multiple Linear Regression (MLR) and Integral Method (IM). Note, even if the IIM and the IM reported different parameters (Table 2), the curves are almost congruent.

Figure 5 shows the correlation plots of the reported IIM and SSM parameter values and confirms that the same solutions were reached in each data set. As the termination of the method is determined by convergence tolerances, minimal deviations occur in each parameter, leading to correlations lower than 1.000, ranging from R = 0.986 – 1.000. The correlation for the MLR, IM and LMA to the SSM are R = −0.209 –0.992, R = −0.448 – 0.999, R = 0.158 – 0.999, respectively.

Figure 5

Correlation Plots. Correlation of the VEM parameter values identified by the Simplex-Search Method (SSM) and by the Iterative Integral Method (IIM) of 26 data sets.

Discussion

The noise-free in-silico analysis showed that all identification methods where capable of identifying the parent parameter values via the simulated model response. This result confirms that the model is structurally identifiable [8, 23], and the identification methods were accurately implemented. When limited white noise was added to the simulated profile all tested methods, except MLR successfully located a similar error minimum with identified parameter values close to their parent values. MLR resulted in significantly higher SSE values in every single data set and parameter values that were highly removed from their parent values. These erroneous results are caused by the noise-amplifying effects of derivative terms in Equation (6). The differentiation terms in Equation (6) are due to the second-order characteristics of the underlying two-compartment model. However, for some simpler models (e.g. linear and non-linear one-compartment models) MLR can be a reasonable and legitimate method for parameter identification [5, 8]. However, this analysis indicates that MLR should not be applied for parameter identification of this VEM or such other higher-order models.

While the IM utilizes much the same methodology as the MLR, the outcomes of the noisy in-silico analysis and the clinical analysis show the merit of integrating the governing equation. The integrals in the IM low-pass filter the measured signals, yielding significant improvement in robustness over MLR for noisy simulation data. However, the IM located some non-physiological parameter values in the clinical data. In total, there were negative parameter values in 8 of 26 (31%) data sets. The erroneous parameter constellations resulted in higher residuals than reported by SSM and IIM methods, but still produced reasonable simulation results. These unacceptable parameter values results are in accordance with the findings of Sato et al [6], who found improvement in integral-based VEM identification over MLR, but still obtained negative parameter values. The IM failure to consistently produce accurate parameter values might be due to un-modeled effects in the p _aw (t) term and the subsequent amplification of error via the comparatively highly parameterized model. Less parameterized models have been robust to such estimation of the simulated profile and successfully identified with the IM [19, 26, 27].

In contrast to the IM, the IIM enabled successful parameter identification in every case tested. Table 2 highlights a case wherein the IM (effectively the first cycle of IIM) yielded a highly un-physiological value for C ₂ . In this case, the IIM converged to a lower error state with more plausible physiological values. By repeatedly evaluating the matrix equation Eq. (12) the p _aw (t) term in the left hand side is updated, whereas the p _aw (t) on the right hand side remains the measured data. Thus, the linear least square error minimization process optimizes the simulated p _aw (t) against the measured values directly. Figure 5 shows that the IIM and the SSM located the same parameter values and Table 1 shows that the methods found the lowest overall error values.

The LMA is a popular parameter identification method for linear and non-linear problems [9, 28, 29]. LMA parameter identification in clinical data was 3.5 times faster than the SSM due to its higher order, but still 14 times slower than the IIM. However, the LMA appears more sensitive to initial values than the SSM. With the same hierarchically derived initial values, the LMA reported local minima with comparatively high SSE values in 9 out of 26 (35%) data sets. Based on this study, the LMA with hierarchical support is not suitable for reliable straightforward VEM parameter identification.

The IIM and SSM located very similar minimal error in every case. In contrast, MLR consistently failed in the presence of measurement error and noise. The IM was unable to report the error minima located by IIM and SSM, perhaps due to un-modeled effects in p _aw . LMA occasionally converged successfully, but located local minima in a significant number of cases. Hence, only the IIM and SSM could be considered as reliable methods for VEM parameter identification. The robustness of the SSM was dependent on the initial values, which were provided here by a hierarchical method. In contrast, the IIM was independent of initial values, which is a significant reduction in complexity and is thus more robust in application.

Limitations

All models are representations of a true system and thus all models have limitations. In this case, the underlying VEM is a basic model and is limited in its abilities to describe all effects in respiratory mechanics, in large part to capture fundamental physiology and also to ensure it is identifiable with readily available data or simple maneuvers. Pathophysiological effects that are not captured by the VEM include for example alveolar recruitment and over-distension, flow limitation, and edema among leading issues that can affect such models and experiments. Assuming these un-modeled effects are not present in the data, a perfect model agreement could be guaranteed, given a perfect experiment.

However, additional un-controlled effects are always present in clinical data due to measurement noise and the clinical impact of these devices on physiological responses. Secondly, inadequate experimental handling e.g. spontaneous muscle activities caused by insufficient sedation interferes with the assumed experimental hypotheses and model assumptions that neglect this activity. With this variety of un-modeled effects in the data, the IM can report inaccurate results that could be corrected by applying the IIM.

For example, assuming a hypothetical impulse of pressure in the measured data at t = t ₂ due to any of the above described reasons, then the integral of the pressure will be affected from t ₂ to t _end in Eq 8. Hence, the IM will have the aberration in both sides of the equation and the model functions that are most prominent in the period t ₂ will be negatively affected. In contrast, the IIM uses a re-simulated model response and thus, the impulse gets smoothed in the left hand side of Eq 8 over one or more iterations. The error induced by the impulse at t = t ₂ is present in the model simulation contained in the left hand side of Eq 8 from t = 0 to t _end . The error is not limited to a particular period of the response and thus the error does not exaggerate any particular model function. With ongoing iterations the influence of the impulse is reduced, and the results of the IIM approach the solution of the gradient-based methods. Hence, the IIM offers built-in model-based compensation for such errors within its framework, as the iterations enforce the model dynamics onto the data mitigating the impact of un-modeled effects. These effects are presumed small, although it should be noted that this presumes the model itself captures the fundamental mechanics accurately.

Small deviations within the reported parameter values of both methods were found (Figure 5). These deviations are well within clinical tolerances and were unlikely to significantly alter therapeutic choices. The lowest correlation (R = 0.986) of the IIM parameters to the SSM parameters were found in R ₂ , which is, together with C ₂ , predominantly sensitive to the relaxation phase [30]. This phase is most-likely to be disturbed by cardiogenic oscillations. Additionally, in cases with no distinct exponential pressure drops during EIP, information content for R ₂ , C ₂ identification is low, decreasing the convexity in the error plane. Low convexity causes an early termination of the tested identification algorithms leading to minor deviation in the corresponding parameters. Changes in tolerances used would also improve these small differences but are with a variation of less than 1% utilizing a TolFun value being 10 times smaller, relatively small.

Technical and clinical impact

The IIM is potentially applicable to other linear respiratory mechanics models of higher order, such as the Inhomogeneity Model [31], as this model poses a transfer function of the same form as the VEM [32]. Furthermore, the advantage of the IIM of being independent on initial values is beneficial for the principle of hierarchical approaches of parameter identification. A robust identification process of more advanced models like the VEM avoids gradient-based methods within this hierarchical level and shifts these more susceptible methods further down in the model hierarchy. Thus, an implementation of IIM would lead to increased robustness and efficiency in hierarchical approaches of more extensive hierarchical model structures [33].

Applying the IIM to the VEM allows online-monitoring the trends of patient-specific viscoelasticity at the bedside, and thus supports controlling the mechanical ventilator to maintain the viscoelastic properties in a normal range. Secondly, using a frequency dependent model supports the adjustment of ventilation frequency at the lowest impedance. Further clinical studies are required to develop a control scheme based on viscoelastic properties of the individual patient’s lung.

Conclusion

Compared to gradient-based methods, the IIM provides robust parameter identification and operates without the initial value problem. Therefore it can offer more operator independence in model parameter identification.

The VEM investigations confirm that parameter identification via the IIM is also superior to error-mapping methods with respect to efficiency and robustness. The IIM takes advantage of the linear mathematical structure of the model and offers fast computing time and maximal robustness. These features offer potential for the IIM to be implemented in an online tool at the bedside for ventilation management.