# Nonlinear identification algorithm for online and offline study of pulmonary
mechanical ventilation

##  1 INTRODUCTION

##  2 Introduction

From the point of view of intensive medicine, mechanical ventilation therapy
is a procedure used to maintain breathing temporarily, for the time necessary
until the patient recovers their functional capacity and can resume
spontaneous ventilation [1]. The ventilatory strategies are based on
calculated parameters and help patient monitoring against adverse effects of
mechanical ventilation [2], [3].

A patient connected to a breathing machine offers the opportunity to compute
various physiological parameters and to evaluate their modifications to
changes in the breathing support strategy [4]. However, ventilators in
intensive care units usually provide to the clinicians a rough estimation of
two key physiological parameters: the airway resistance and the compliance of
the breathing system [5]. Standard equipment typically computes them
considering no more than five points of the total data available from the
respiratory cycles: the Plateau Pressure, the Peak Inspiratory Pressure (PIP),
the Peak Inspiratory Flow (PIF), the Total Inspiration Volume, and the
Positive End-Expiratory Pressure (PEEP) [6]. A few devices make use of the
whole cycle data, such as Hamilton ventilators, which obtain values for the
two parameters through least squares fitting based on a linear model of the
system [7].

In this sense, model-based methods which use all the data from the cycle
measured by the ventilator fit the respiratory behaviour better, giving
clinicians a deeper physiological picture of the patient’s state [8], [9],
[10], [11], [12]. Such models can be used to optimise the ventilator settings,
i.e., to reduce harmful effects of the ventilation therapy or to keep the
patient within a “zone of comfort” during mechanical ventilation [13].
Nevertheless, most of those models characterise the dynamics of the pulmonary
system as linear, simplifying their description at the expense of accuracy
[14], [15].

In this paper, a more accurate nonlinear model is considered to better
represent the respiratory system of a sedated patient under assisted
ventilation, especially when there is pulmonary disease, as the physiology of
their lungs may change remarkably and cause them to breath in the nonlinear
regions [16]. Thus, an algorithm based on nonlinear identification techniques
is proposed to estimate the model parameters, which uses a hierarchical
identification process. Firstly, the parameters of a linear version of the
model are estimated and then, those parameters are used for the initialisation
and identification of the nonlinear model. In this way, with the flow and the
pressure data measured in the patient’s mouth, the algorithm computes a set of
parameters of the nonlinear model for each respiratory cycle.

This work comprises two stages. In the first one, a simulated patient is
initially used to assess the algorithm efficiency, considering nonlinear
sigmoid models. Then, an even more realistic nonlinear model which includes
the occurrence of local hystereses is implemented, to evaluate the capability
of the algorithm to get information from more realistic data. After the
successful in-silico validation, in the second stage the proposed algorithm is
applied to real respiratory systems of sedated patients with a breathing
disease, particularly COVID-19. To this end, measurements of flow and pressure
obtained from PEEP titration manoeuvres are utilised, since during these
manoeuvres the patient is ventilated in a wide range of breathing regions.

The main contributions of this work are listed as follows:

  * •

An algorithm that carries out an identification process to obtain the
nonlinear model parameters of a patient under assisted ventilation is
developed.

  * •

The algorithm can determine the respiratory region where the patient is
ventilated, i.e. Atelectasis, Linear, and Overdistension Regions.

  * •

The proposed parameter identification algorithm is validated on extensive
simulated data (including noise corrupted data), which uses a sigmoidal with
hysteresis behaviour to describe patient compliance.

  * •

The algorithm is tested using actual patients’ data, and the results confirm
its potential value in critical care.

The remainder of this article is organised as follows. Section 3 introduces
respiratory physiology to understand the mathematical models used in this
work, which are explained in Section 4. Section 5 presents the developed
algorithm used to obtain the parameter values of the models. Section 6
presents the results obtained using simulated and actual patients under
assisted ventilation. Finally, Section 7 concludes the paper, highlighting the
main findings and describing future research directions.

##  3 Respiratory physiology

To understand the models that describe the behaviour of the pulmonary system,
some respiratory physiology concepts are explained in this section.

Air moves due to a pressure gradient. If external pressure (pressure in the
mouth) is higher than internal pressure (pressure in respiratory system),
entry of air into the pulmonary system occurs (inhalation). By contrast, if
external pressure is less than the internal, the system expels air
(exhalation). These changes in pressure can occur spontaneously (carried out
by respiratory muscles) or artificially (performed by a mechanical
ventilator).

###  3.1 Spontaneous ventilation

Considering atmospheric pressure as the base reference level, the respiratory
system achieves the entry of air from variations of the internal pressure,
generated by the respiratory muscles. In spontaneous ventilation, this
pressure variation is responsible for the instantaneous value of the air
volume moving in or out of the lungs during the respiratory process.

Defining the instantaneous volume of the lung at time t as $V_{L}(t)$, it can
be observed in Fig. 1 that the total amount of air that enters with each
breath (known as Tidal Volume $V_{T}$) is considerably less than the total
volume that can be forcibly inhaled from a maximal inspiration following a
maximal exhalation, called Vital Capacity. Defining the Residual Volume (RV)
as the amount of air that remains in a person’s lungs after a maximum
exhalation, the Total Lung Capacity (TLC) is the total lung volume,
considering the Vital Capacity plus the Residual Volume, as it is schematised
in Fig. 1 [17]. Finally, after a normal passive exhalation, the amount of air
remaining in the lungs is called Functional Residual Capacity (FRC), which is
a volume greater than the Residual Volume. The term $V(t)$ is used to denote
$V_{L}(t)-FRC$.

Figure 1:  Characterisation of volumes in the system.

###  3.2 Mechanical ventilation

In the event that the patient is unable to breathe on their own, either
totally or partially, mechanical ventilation can be used to maintaining their
good oxygenation and monitoring their lung mechanics.

The pressure imposed by the ventilator at the patient’s mouth, $P(t)$, is
composed by a variable component, $P_{v}(t)$, mounted on a programmable
constant pressure, the PEEP (i.e., $P(t)=P_{v}(t)+PEEP$). $P_{v}(t)$
establishes the aforementioned variable volume $V(t)$, ensuring a predefined
respiratory cycle set by the clinician. On the other hand, the PEEP is a
ventilator parameter, also configurable by the clinician, used to attain that
at the end of expiration the lung does not return to atmospheric pressure,
hence preventing the collapse of alveoli every cycle [18] [19]. The PEEP,
adequately set, can decrease the respiratory work and improve the patient’s
oxygenation, since it may increase the FRC of the patient [20].

Most of the ventilation machines can work in one of the following two types of
operation. In the first one, assisted ventilation, the machine assists the
patient after the detection of an inspiratory effort and so, the patient can
control the timing and size of the breath. In the second type, controlled
ventilation, a patient effort is not required, and the ventilator starts
automatically every respiratory cycle with a predefined period [21]. When in
this operation type, the patient needs to be sedated to avoid asynchronies.

The breathing pattern produced by the ventilator requires the setting of
control parameters such as the size of the breath, how fast and how often air
is brought in and let out, and how much effort, if any, the patient must exert
to signal the ventilator to start a breath. Numerous ‘ventilation modes’ have
been developed to make ventilators produce breathing patterns that coordinate
the machine’s activity with the needs of the patient [21]. The best known are:

  * •

Volume control ventilation (VCV). In this ventilation mode, the physician sets
the volume of air that will enter the patient in each respiratory cycle. The
ventilator will then establish a flow profile for inspiration, and the
pressure in the mouth will depend on the patient’s respiratory mechanics.

  * •

Pressure control ventilation (PCV). The physician sets the peak value of
pressure that the ventilator will reach in each respiration cycle. In this
case, the amount of air volume entering the lungs is not fixed, but instead
depends on the respiratory mechanics of the corresponding patient.

###  3.3 Hysteresis and Pressure-Volume curves

The Quasi-Static Pressure-Volume (Q-S P-V) curve represents the limits of
pressure and volume inside of which respiratory cycle paths can take place
[22] (see Fig. 2). It covers a range from FRC to TLC and its contour can be
empirically obtained by measuring pressure in the mouth while maintaining
lengthy low and constant inspiration and expiration flows, so that the
pressure drops in the airways are negligible. Under these conditions, $P(t)$
equals the alveolar pressure $P_{A}$. [23] [24].

A physiological phenomenon that can be observed in a patient is the
hysteresis. Hysteresis is defined as different volume paths upon application
and withdrawal of a certain pressure, and reflects an imperfect elastic
response to deformation [25] [26].

Additionally, two characteristic inflection points can be appreciated in the
curve, denoted as LIP (i.e. Low Inflection Point) and UIP (i.e. Upper
Inflection Point) [27]. They define the three main respiratory regions of the
patient: (a) the Atelectasis Region: from the LIP zone downwards, the lung is
poorly recruited and many alveoli must be opening and closing with each
respiratory cycle [28]; (b) the Linear Region: between the LIP and the UIP, is
the desired breathing region where alveoli are highly recruited plus their
elasticity is at peak performance; and (c) the Overdistension Region: from the
UIP zone upwards, the lung is functioning near its elastic limit and
overdistension of some alveoli units can take place, causing injury to the
patient [29].

Figure 2:  In dashed lines, the Q-S P-V curve measured from FRC to TLC is
shown. In solid lines, illustrative respiratory cycles curves starting from
different equilibrium points $Eq_{i}$, $i=1,2,3$, showing the various curve
shapes that can be observed depending on the region where the patient is
ventilated. $P_{A}$ is the alveolar pressure, $V_{L}$ is the total volume and
$V_{PEEP_{i}}$ is the volume due to the ventilator $PEEP_{i}$.

Some illustrative respiratory hysteresis cycles are depicted in each one of
the respiratory regions in Fig. 2. Such respiratory cycles start on different
initial equilibrium points ($Eq_{i}$, $i=1,2,3$ in Fig. 2), which depend on
the patient’s FRC plus the extra air volume due to the ventilator PEEP [30].
From a given equilibrium, Volume $V(t)$ evolves in accordance with the
variable component ($P_{v}(t)$) of the ventilator pressure. A diagram of a
generic respiratory cycle hysteresis which starts from an $Eq_{i}$ is shown in
Fig. 3.

Figure 3:  Hysteresis in the lungs-chest wall system in a typical respiratory
cycle of a patient. The arrow pointing upwards corresponds to the inspiration
path, and the arrow pointing downwards to the expiration path.

##  4 Respiratory system modelling

In this section the equation of motion that relates pressures in the
respiratory system with the volume of air moving in and out of lungs-chest
wall system is presented. Subsequently, the corresponding linear and nonlinear
models are introduced, whose parameters will be estimated by the proposed
identification algorithm.

###  4.1 The equation of motion

The breathing dynamics of a sedated patient can be described by relating the
pressure applied in the patient’s mouth with the volume and flow in the
respiratory system.

The elastic properties are determined by the relationship between lung volume
and pressure in the lungs-chest wall system. The resistive properties are
represented by the airway resistance. Such relations, starting from a given
$Eq_{i}$, are shown in Eq. (1), commonly known as Equation of Motion [31], and
its electrical equivalent model is represented in Fig. 4:

$P_{v}(t)=\dot{V}(t)R_{aw}+P_{c}(V(t))$ (1)

where the time derivative $\dot{V}(t)$ of $V(t)$ corresponds to the airflow
$F(t)$ into the lungs. $R_{aw}$ represents the airflow resistance of the
airways and $P_{c}(V(t))$ is the variable pressure in the lungs-chest wall
system, which depends on $V(t)$ and on the elastic properties of the region
where the patient is being ventilated. To simplify notation, the time
dependency will not be written in the following.

Figure 4:  Electrical equivalent model of the respiratory system.

In the electrical equivalent model shown in Fig. 4, the capacitor
characterised by $P_{c}(V)$ corresponds to the lung-chest wall compliance.

A widely accepted characterisation of the lung-chest wall system in the
complete volume range from the FRC to the TCL corresponds to a sigmoid curve
[23]. However, when describing normal either healthy or unhealthy breathing
cycles, simpler ‘local’ descriptions can be used, as the sigmoid curve can be
well approximated via piecewise continuous polynomial functions depending on
the $V$ and the respiratory region.

A linear model for $P_{c}(V)$ provides a very accurate approximation in the
Linear Region, but a coarse one outside of it. Thus, in this work, a quadratic
(nonlinear) model for $P_{c}(V)$ is used for identification, to be able to
obtain an accurate approximation also in the Atelectasis and Overdistension
Regions. Both the Linear and Nonlinear identification models are presented in
sections 4.2, and 4.3, respectively.

###  4.2 Respiratory System Linear Model (LM)

The brief physiological explanation of the respiratory system in subsection
3.3 helps understanding the complexity of accurately modelling its dynamics.
However, there are simple models that still provide the most significant
information. In that sense, the linear model is the most widespread model
among clinicians. According to it, the relationship between $P_{c}$ and $V$ is
characterised as:

$P_{c}(V)=\frac{V}{C}=VE$ (2)

where $C$, the compliance of the respiratory system, is a constant which
depends on the region of the P-V curve where breathing is taking place and, of
course, on the health condition of the patient. The inverse of the compliance
is known as elastance ($E=1/C$). The value of $C$ is maximum when the patient
is ventilated in the Linear Region.

Thus, using the linear description of Eq. (2) and the equation of motion, Eq.
(1), the dynamics of the linear model can be expressed as:

$\dot{V}=-\frac{1}{R_{aw}C}V+\frac{P_{v}}{R_{aw}}$ (3)

It can be observed that two parameters characterise the patient’s respiratory
system in the region where breathing is occurring:

$\theta_{l}=[C,R_{aw}]$ (4)

They will be estimated via the identification algorithm developed and
presented in this work, by using the airflow $F=\dot{V}$ and pressure $P_{v}$
measured on the patient.

###  4.3 Respiratory System Nonlinear Model (NLM)

A better characterisation of the relation between $P_{c}$ and $V$ can be
obtained by using a quadratic description:

$P_{c}(V)=V(a_{1}+a_{2}V)$ (5)

where the elastance varies with $V$. Eq. (5) provides improved fit over that
of the linear model, particularly when the lungs are ventilated in the areas
near the inflection points (UIP or LIP) [31]. Additionally, the concavity of
the curve obtained in Eq. (5) can provide information related to the patient’s
ventilation region, i.e. Overdistension, Linear or Atelectasis.

In this case, the nonlinear state equation that the estimation algorithm will
use is:

$\dot{V}=-\frac{1}{R_{aw}}(a_{1}V+a_{2}V^{2})+\frac{P_{v}}{R_{aw}}$ (6)

and the three parameter vector to be estimated is:

$\theta_{nl}=[a_{1},a_{2},R_{aw}]$ (7)

##  5 Nonlinear identification algorithm

In this section, the proposed systematised identification process is
presented. Beginning with the estimation of the linear model, it then uses
these parameters to initialise and estimate those of the nonlinear quadratic
model.

The algorithm core is based on a grey-box modelling approach, since the
mathematical models of the respiratory system are available (Eq. (2) and Eq.
(5)) [32] [33]. The parameters vector values of both models are obtained using
the Levenberg-Marquardt algorithm, which is an iterative method for nonlinear
least squares problems, widely adopted in a broad spectrum of disciplines [34]
[35].

The algorithm presented in this work consists of three stages:

  * •

Signals acquisition, where pressure, flow and volume signals are measured.

  * •

Parameter estimation, where the signals are used to obtain the nonlinear set
of parameters, one for each respiratory cycle.

  * •

Fit evaluation, where the output of the estimated models is compared with the
real output of the system. The model is saved or discarded depending on a
normalised root-mean-square-error criterion.

###  5.1 Signal acquisition

The algorithm uses the pressure and flow signals measured in the mouth of a
totally sedated patient. These data have been obtained with the acquisition
device FluxMed® GrE with a sampling frequency of 256 Hz.

The signal pressure, $P_{v}$, is measured using a differential pressure
sensor, which provides the variable pressure at the mouth.

On the other hand, the flow signal $F$ is measured through a fixed-orifice
pneumotachograph, in which the pressure difference on both sides of the
orifice is proportional to the flow value squared.

$V$ is obtained indirectly, using trapezoidal integration of the flow signal.
This numerical integration is performed for each respiratory cycle.

To close this stage, an auxiliary variable $Aux$ is set equal to zero, when
the end of cycle is detected. This variable will be used, during subsequent
stages, to prevent a possible loop when the estimation fails.

###  5.2 Parameter estimation

The proposed parameter estimation process can be summarised as follows (see
the flow chart in Fig. 5 for further clarification):

Figure 5:  Proposed algorithm flow chart.

  * •

Linear model initialisation. The algorithm computes initial values for both
linear parameter, $C$ and $R_{aw}$, using data from the first cycle of the
registered signals, as indicated:

$\displaystyle C=\frac{V_{Tinsp}}{P_{Plat}-PEEP}$ (8) $\displaystyle
R_{aw}=\frac{PIP-P_{Plat}}{PIF}$ (9)

where $V_{Tinsp}$ is the maximum tidal air volume that enters during the
inspiration semicycle, and $P_{Plat}$ is the pressure value when $F=0$ (in the
change from inspiration to expiration).

  * •

LM identification for NLM initialisation. The algorithm finds the optimal
values of $\theta_{l}$ for the linear model through identification, aiming to
minimise the quadratic cost function:

$M_{N}(\theta_{l})=\sum_{k=1}^{n}(\hat{V}-V)^{2}$ (10)

where $n$ is the number of samples of the respiratory cycle and
$\hat{V}=f(P_{v},\theta_{l})$ is the tidal volume computed by the model (at
this step, the linear one) with $P_{v}$ as input signal.

  * •

Nonlinear quadratic model initialisation. The initial values for the nonlinear
model identification are set based on the linear parameters estimated in the
above step. So, the previous value of $R_{aw}$ is directly utilised for
initialisation, while the initial values of parameters $a_{1}$ and $a_{2}$ are
obtained from Eq. (6), which can be rewritten as:

$P_{v}-FR_{aw}=(a_{1}V+a_{2}V^{2}).$ (11)

Then, the initial values of the nonlinear description of $P_{c}$, $a_{1}$ and
$a_{2}$, are computed by least squares adjustment.

  * •

Nonlinear quadratic model estimation. With the $\theta_{nl}$ vector properly
initialised, the algorithm starts an identification process to obtain the
nonlinear model (Eq. (10) is minimised using the output of the NLM,
$\hat{V}=f(P_{B_{v}},\theta_{nl})$).

For subsequent cycles, the algorithm initialises the identification process
with the parameter values of the NLM computed in the previous respiratory
cycle.

###  5.3 Fit evaluation

In this stage, the estimated nonlinear model for each respiratory cycle is
validated using the normalised root-mean-square-error ($NRMSE_{\%}$)
criterion:

$NRMSE_{\%}=100\left(1-\frac{||V-\hat{V}||}{||V-\bar{V}||}\right)$ (12)

where $\bar{V}$ is the mean value of $V$. In this work, the validation stage
uses a threshold to decide whether to keep or discard the estimated model.

If the result is higher than this threshold, the algorithm saves the estimated
nonlinear model parameters. Then, it executes the Signal Acquisition Block,
subsequently initialising the next nonlinear model estimation. Note that
different approaches can be used to design the threshold. In the present case,
the threshold was set equal to the $NMRSE_{\%}$ obtained from the standard
linear model as computed by the ventilator (according to (8) - (9)).

On the other hand, if the result is lower than the threshold, the estimation
process was unsuccessful. In this case, if $Aux=0$, the algorithm goes to
Linear Based Initialisation for a second chance with different initial
condition. By contrast, if such second chance fails ($Aux=1$), the algorithm
discards the current respiratory cycle and continues with the next one, by the
Signal Acquisition Block.

In the specific case that the estimation fails for the first cycle, the
algorithm discards it instead of returning to Linear Based Initialisation.

##  6 Results from simulated and real patients

This section presents results obtained with simulated and real patients data.
The first subsection is intended to validate the effectiveness of the proposed
algorithm for the estimation of pulmonary system parameters. To this end, the
nonlinear quadratic identification algorithm is thoroughly assessed using
multiple signals generated by simulation and representative results are
displayed. Once the suitability of the estimation algorithm is proved, it is
used, in the second subsection, to characterise actual pulmonary parameters
computed using data obtained from real sedated patients with COVID-19 under
different conditions of mechanical ventilation.

###  6.1 Algorithm validation with simulated patients

Patients are simulated to obtain airflow and pressure signals to evaluate,
analyse and validate the algorithm. The main objectives are, firstly, to
confirm the superiority of the quadratic (nonlinear) identification over the
linear one, to locally approximate a simulated patient characterised with a
sigmoidal model, particularly in the Atelectasis and Overdistension Regions.
Secondly, to draw conclusions regarding the degree of information that the
estimator can provide when dealing with models with hysteresis, which render a
more accurate representation of real patients.

$\boldsymbol{(a)}$ For the first objective, the proposed identification
algorithm is tested using data acquired from a sigmoidal-modelled patient
without hysteresis, i.e., the P-V curve of the simulated patient was
programmed according to the sigmoid curve equation [27]:

$V=a+\frac{b}{1+e^{-(P_{A}-c)/d}}$ (13)

where a, b, c, d are constants, and $V$ goes from 0 to TLC-FRC.

As input signal for the simulation trials, actual pressure measured from a
sedated patient was utilised (see 10 illustrative $P_{v}$ cycles in Fig. 6).
Simulations were carried out to obtain the corresponding flow signal $F$ at
different PEEP values, near the LIP with PEEP=4 $cmH_{2}O$, in the Linear
region with PEEP=13 $cmH_{2}O$ and near the UIP with PEEP=22 $cmH_{2}O$.

Figure 6:  Pressure from actual sedated patient, used as input to the
simulator.

The algorithm was run for every simulation, so that linear and nonlinear
estimated models were obtained. Figure 7 shows the $P_{c}(V)$ of both LM in
black, NLM in red, superimposed on the P-V curve of the simulated patient (in
blue). It can be appreciated that the quadratic model fits better than the
linear model when the patient is ventilated in LIP ($94.40\%$ vs $99.01\%$)
and UIP region ($89.10\%$ vs $97.28\%$). As expected, for the patient
ventilated in the linear region, both models achieve similar good fit
($99.28\%$ LM, $99.60\%$ NLM).

Figure 7:  P-V curve of the sigmoid-modelled simulated patient (blue). The
$P_{c}(V)$ estimated by the algorithm (LM in black, NLM in red) are also
shown.

Subsequently, to further assess the proposed algorithm under non-ideal
conditions, uniformly distributed noise was added to both $P_{v}$ and $F$
signals, to account for measurement and digitalisation errors introduced by
the sensors and the ADCs of the FluxMed® GrE [36]. For these tests,
simulations were conducted considering three patients, each one modelled with
different sigmoidal P-V curves and $R_{aw}$ values and selected PEEP values
were used to cover the whole ventilation range.

The fits obtained with the estimated models are presented in Fig. 8, where the
blue bars correspond to the NLM, and the green ones to the LM. It can be seen
that even with high levels of noise, the estimated quadratic models have
really good and higher fit values than the linear ones in the LIP and UIP
Regions, while in Linear Region both models have similar good results.

Note that the noise considered in this tests significantly exceeds the
measurement and digitalisation errors bounds provided by the manufacturer of
the FluxMed® GrE, which are $3\%$ for the $P_{v}$ signal and $3.5\%$ for the
$F$ signal. Consequently, these results confirm the proficiency of the
proposed nonlinear quadratic based algorithm for the identification in the
different regions of sigmoidal-modelled patients, even under non-ideal
conditions.

Figure 8:  Fit results obtained by the algorithm using input signals with
different percentages of added noise (column indicators at the top)

. The simulated patients were ventilated in LIP, linear and UIP regions. Bars
in blue correspond to the NLM and green bars to the LM.

$\boldsymbol{(b)}$ Next, the second objective of evaluating the algorithm
capability to draw information from more realistic data is addressed. To that
effect, the nonlinear quadratic identification algorithm was tested with data
obtained from 10 simulations of a patient with nonlinear-modelled lungs that
present hysteresis in the $P_{c}(V)$ curve, using the two ventilation modes
(see section 3.2).

Several breathing zones were considered, covering the range from FRC to TLC,
with representative hysteresis curves for each of the three respiratory
regions. The aforementioned pressure $P_{v}$ from an actual sedated patient
(Fig. 6) was also used as input signal in each simulation, to obtain the
corresponding flow signals. The simulated curves are presented in dashed red
lines in Fig. 9. The hysteretic characteristics were constructed based on the
descriptions provided in reference [30]. For their design, it was ensured that
both $V$ and $P_{c}$ fell within actual acceptable breathing ranges.

Then, the corresponding pairs $P_{v}$-$F$ were fed to the algorithm, and the
estimation process was conducted. The resultant estimated $P_{c}(V)$ curves
are depicted in black in Fig. 9, where it is worth noting that they are mostly
confined inside their corresponding modelled hysteresis curve, and their
corresponding estimated parameters are shown in Table 1.

(a) Simulated Patient ventilated in volume control mode. (b) Simulated Patient
ventilated in pressure control mode.

Figure 9: In blue, P-V curve programmed in the simulator. The simulated respiratory cycles (dashed red) can be compared with the $P_{c}-V$ curves estimated by the NLM of the proposed algorithm (in black). Sim # | $a_{1}$ | $a_{2}$ | $a_{1}V_{T}$ | $a_{2}V_{T}^{2}$
---|---|---|---|---
1 | 0.026 | 1.45e-4 | 6.7 | 9.7
2 | 0.025 | 1.84e-5 | 7.4 | 4.1
3 | 0.020 | 1.38e-6 | 9.6 | 0.3
4 | 0.052 | -5.09e-5 | 19.1 | -6.8
5 | 0.084 | -1.02e-4 | 26.9 | -10.5
6 | 0.036 | 9.7e-5 | 8.6 | 5.6
7 | 0.028 | 2.80e-5 | 7.6 | 2.1
8 | 0.015 | 3.36e-6 | 8.4 | 1.1
9 | 0.040 | -7.29e-6 | 14.1 | -0.8
10 | 0.056 | -4.83e-5 | 14.6 | -3.3
Table 1: Identification results: second and third columns show the $P_{c}(V)$
estimated parameters (in $cmH_{2}O/ml$ and $cmH_{2}O/ml^{2}$, respectively).
Fourth and fifth columns show the ranges for the linear and quadratic pressure
terms of $P_{c}$, respectively (in $cmH_{2}O$).

It can be observed that the concavity of the estimated curves has a direct
relation with the respiratory region where the patient is ventilated. For
instance, when the estimated parameter $a_{2}$ is negligible, it means that
the patient is ventilated in the Linear Region, while $a_{1}$ gives
information about the pulmonary compliance. On the other hand, the more
negative $a_{2}$, the more into the Atelectasis Region is the patient
breathing. Conversely, the more positive $a_{2}$, the higher the risk of
alveolar overdistension. Moreover, the physicians could use the estimated
parameters together with empirical information of the patient (for instance,
from titration manoeuvre data) to have a better knowledge of their condition.
In this way, they could establish quantitative ranges for $a_{1}$ and $a_{2}$
that more precisely indicate the respiratory region where the ventilation is
conducted. It seems likely that further analysis could provide a relation
between those parameters and the patient’s Stress Index, whose value
quantitatively represents whether or not the patient is being ventilated in
the linear region [37].

With the above satisfactory results at hand, the proposed nonlinear
identification algorithm proves to be a powerful tool to extract information
from the ventilated patient. Therefore, it will be used next with real
patients data.

###  6.2 Algorithm applied to actual patients data

In this section, the algorithm is tested with input-output signals ($P$-$F$)
from sedated patients under assisted ventilation, measured during a PEEP
titration manoeuvre.

Such manoeuvre is intended to recruit alveoli, aiming to ventilate the patient
in the linear region, where both the alveolar collapse and the overdistension
are minimum. In this way, the pulmonary homogenisation would be maintained,
protecting the lungs from mechanical ventilation injuries [38] [39] [40]. It
consists of increasing the PEEP value in stages until reaching a maximum
admissible pressure and then decreasing it, enabling to detect the highest
compliance region in the process.

Three case studies are considered, shown in Fig. 10. In the first two cases,
the $P$ signals are measurements from the same patient, under PEEP titration
manoeuvres conducted two days apart. Case 3 corresponds to a different
patient, also under a PEEP titration recruitment manoeuvre. Note that to
facilitate results interpretation, particularly to visualise the correlation
between Fig.10 and Fig.12, the signals were distinguished by coloured
sections, where grey corresponds to the highest level of PEEP in the
manoeuvre, and black to the lowest one.

Figure 10: $P$ signal measured on two patients during PEEP titration
manoeuvres. Cases 1 and 2 correspond to the same patient, two days apart.

Fig. 11 comparatively displays the fit obtained with both identification
models (LM in blue and NLM in red). For the three cases under study, the NLM
demonstrates a better fit than the LM when the patient is ventilated with high
PEEP levels (corresponding to the Overdistension Region or its immediate
environs). Conversely, when the PEEP decreases and the patient is ventilated
in the Linear Region, the fits of both models are similar.

In this context, Case 3 is particularly illustrative. It can be appreciated
that a considerable worse fit is obtained with the LM. From this information,
it can be inferred that the PEEP titration manoeuvre was conducted mainly in
the Overdistension Region, rather than in the Linear one.

Figure 11: Fit results obtained with the linear model (in blue) and the
nonlinear model (in red), using the NRMSE criterion.

The P-V curves results provided by the proposed estimation algorithm are
depicted in Fig. 12, subfigure (b), while in subfigure (a) the results using
the linear model of the respiratory system are shown. For the sake of clarity,
the curves obtained with the ascendant PEEP series are presented in the left-
hand side column and those obtained with the descendant PEEP series are in the
right-hand side column.

(a) Results using the linear model. (b) Results using the quadratic model
(NLM).

Figure 12: $P_{c}-V$ curves identified by the algorithm in each respiratory
cycle of patients measured during PEEP titration manoeuvres. Cases 1 and 2
correspond to Patient #1 and Case 3 to Patient #2. The left-hand side column
shows the curves estimated with the ascendant PEEP series, and the right-hand
side the curves estimated with the descendant PEEP series.

Note that the slopes of the lines obtained with the linear model (Fig. 12(a))
give information about an “average” compliance for the different regions where
the patient is ventilated, it being expected that the highest corresponds to
the Linear Region. This ”average” compliance is typically used by the
physicians, but it is incomplete if compared to the information available from
the proposed quadratic identification algorithm. In effect, the slopes of the
linear identification on their own, in certain cases, do not suffice to
determine if the patient is outside the Linear Region. The information
regarding curvature is of decisive importance.

The best way to illustrate this situation is with Case 3. A standard PEEP
titration manoeuvre was performed on Patient #2. As a result, a series of
straight lines were obtained with the LM (Fig. 12(a) bottom). An analysis of
their slopes (relatively low, if it is assumed the titration has covered a
great deal of the Linear Region) might lead to the conclusion that the patient
suffers from a severe reduction of compliance. The actual situation is quite
different and it can be easily detected from the curves in Fig. 12(b), Case 3.
As it was previously mentioned, the curvatures reveal that the titration
manoeuvre was conducted mainly in the Overdistension Region.

The following provides a more detailed analysis of the results obtained with
the proposed nonlinear identification algorithm. Focusing the attention in Fig
12(b), in Cases 1 and 2 (Patient#1), right column, a slight change of
concavity can be observed (concave left, i.e. negative $a_{2}$, for the lowest
values of PEEP). This means that the patient is entering into the Atelectasis
Region. Therefore, it is not desirable to go further below that minimum PEEP
with this patient. On the other hand, in Case 3 (Patient#2) the concavity has
not yet changed ($a_{2}$ is positive) for the lowest PEEP value. With this
result, the physician can infer that there is still room for PEEP reduction in
Patient#2’s titration manoeuvre.

For its part, information of the top end of the PEEP titration manoeuvre can
be easily obtained from the left-hand side column of Fig. 12(b). For instance,
the pronounced curvature of Patient#2 during the highest PEEP is a clear
indicator that the patient is ventilated in the Overdistension Region.
Consequently, the maximum PEEP applied during its titration manoeuvre is
possibly excessive. This is the kind of online feedback that clinicians can
use to make the decision whether or not to increase the PEEP, since in this
case the patient could generate a lung injury due to the overdistension of
alveolar units.

The aforementioned is not the only information about the lungs condition that
can be extracted from the algorithm results. The slopes of the curves in Fig.
12(b) provide an insight on the patient’s compliance. Moreover, they are
inversely related to parameter $a_{1}$, corresponding this parameter to the
lung’s elastance in the Linear Region.

Therefore, the steeper the slope (lower $a_{1}$), the higher the compliance.
For instance, comparing Patient#1 against Patient#2, it is apparent that the
lung of the latter presents a lower compliance (or expandability), requiring
more pressure to move the same amount of tidal volume.

Such comparative analysis can also help to assess the evolution of a single
patient over time. For this example, let’s compare Case 1 and Case 2,
recalling that those results were obtained from Patient#1’s data collected
with a difference of two days. Examining the curves in Fig. 12(b), it can be
appreciated that the patient’s condition is relatively stable. However, from
direct observation of the slopes, hints of minor decline seems to be
manifesting. This can be better confirmed by plotting parameters $a_{1}$ and
$a_{2}$ against the PEEP in Fig. 13. Effectively, a growing trend in $a_{1}$
is discernible, implying a compliance reduction. Additionally, comparative
inspection of parameter $a_{2}$ shows that, two days later, Patient#1 is
approaching the Atelectasis Region with higher values of PEEP (note that in
Case 2, the sign of $a_{2}$ is negative for PEEPs below 15 $cmH_{2}O$).

Figure 13: Comparison of parameters $a_{1}$ (above) and $a_{2}$ (below)
estimated for Cases $1$ and $2$, according to the corresponding PEEP value.

Of course, for conclusive medical actions, more complementary data and
clinical evidence need to be taken into consideration, but these examples
illustrate how the information provided by the proposed algorithm may assist
health professionals with decision-making tasks, treatment and diagnosis.

As an additional remark, comparing the employed techniques, based on the
Levenberg-Marquardt algorithm, with the well-known least square fitting (LSF)
techniques, our approach achieves higher fit to the data, but at the expense
of more processing time. Depending on the microprocessor and programming, this
extra time may cause that one cycle every other is not identified. This is not
significant, as it is not an uncommon loss in respiratory monitoring.

##  7 Conclusions

An algorithm capable of estimating, online and offline, the internal
parameters of the pulmonary system of patients under assisted ventilation has
been presented. With noninvasive measurements of pressure and flow at the
patient’s mouth and nonlinear identification techniques, the algorithm
estimates the parameters of the quadratic model that best fits the patient in
each respiratory cycle.

Firstly, the performance of the nonlinear identification algorithm was
analysed and validated via simulation. At this stage, the proposed
hierarchical identification algorithm was tested using a nonlinear sigmoidal
modelled patient. Results show that both approaches, the LM and the NLM,
provide excellent fits in the Linear Region. However, an important improvement
of the fit is attained with the proposed NLM based estimation method when the
patient was ventilated outside the linear region. Furthermore, to account for
even more stringent and realistic conditions, the algorithm was subsequently
tested considering noise, resulting in a remarkably better performance of the
NLM approach over the LM in the Atelectasis and Overdistension regions.

Once the algorithm proficiency was established, hysteresis was included in the
simulation models attaining a more realistic characterisation of the patients.
From those tests, it was concluded that the parameters estimated by the
algorithm, $R_{aw}$, $a_{1}$ and $a_{2}$, can provide very useful online
information about the airflow resistance, the lung’s compliance and the region
where the patient is currently ventilated, respectively.

Boosted by those successful results, in a final stage, the proposed algorithm
was used to characterise actual sedated patients with COVID-19 under assisted
ventilation, during PEEP titration manoeuvres. Three case studies were
considered, where two of them correspond to the same patient, in different
days, and the last one to another patient. From the concavity of the curves
obtained for each respiratory cycle, quantified by parameter $a_{2}$, the
algorithm proved to be a useful tool to determine when the patient moves out
from the Linear Region, into the Atelectasis or the Overdistension Regions.
Also, the slope of the curve and its evolution, quantified by parameter
$a_{1}$, supply information to the clinician about the thoracic-pulmonary
compliance and its progress over time.

As future work, it is planned the use of the proposed estimation tool on a
large number of patients, to correlate pathologies, health evolution and other
representative data, which will allow for comparison and a comprehensive
study. It is expected that this will provide quantitative indicators to better
characterise the patient’s status and the ventilation conditions, helping the
physician to better diagnose, to adjust therapy and to improve recovery.

## Competing Interest

The authors have no relevant financial or nonfinancial competing interests to
disclose.

## Founding

This work was carried out with the support of the Facultad de Ingeniería,
Universidad Nacional de La Plata (UNLP project 11-I255), CONICET (PIP
112-2020-010281CO and PIP 112-2020-010281CO), I + D + i Agency (PICT N°
2018-03747), Bunge y Born Foundation - I + D + i Agency (Project COVID-19 #
873) and Hospital I.E.A y C. San Juan de Dios de La Plata.

## References

  * [1] Gattinoni L, Marini JJ, Collino F et al.  The future of mechanical ventilation: lessons from the present and the past.  _Critical Care_ 2017; 21(1): 1–11.
  * [2] Grasselli G, Cattaneo E, Florio G et al.  Mechanical ventilation parameters in critically ill COVID-19 patients: a scoping review.  _Critical Care_ 2021; 25(1): 1–11.
  * [3] Silva PL and Rocco PR.  The basics of respiratory mechanics: ventilator-derived parameters.  _Annals of Translational Medicine_ 2018; 6(19): 376.
  * [4] Borges DL, Nina VJdS, Costa MdAG et al.  Effects of different PEEP levels on respiratory mechanics and oxygenation after coronary artery bypass grafting.  _Brazilian Journal of Cardiovascular Surgery_ 2013; 28: 380–385.
  * [5] Walter JM, Corbridge TC and Singer BD.  Invasive mechanical ventilation.  _Southern Medical Journal_ 2018; 111(12): 746.
  * [6] Baconnier PF, Carry PY, Eberhard A et al.  A computer program for automatic measurement of respiratory mechanics in artificially ventilated patients.  _Computer Methods and Programs in Biomedicine_ 1995; 47(3): 205–220.
  * [7] Iotti G and Braschi A.  _Measurements of Respiratory Mechanics During Mechanical Ventilation._ Rhazuns, Switzerland: Hamilton Medical Scientific Library, 1999.
  * [8] Schranz C, Knöbel C, Kretschmer J et al.  Hierarchical parameter identification in models of respiratory mechanics.  _IEEE Transactions on Biomedical Engineering_ 2011; 58(11): 3234–3241.
  * [9] Redmond DP, Chiew YS, Major V et al.  Evaluation of model-based methods in estimating respiratory mechanics in the presence of variable patient effort.  _Computer Methods and Programs in Biomedicine_ 2019; 171: 67–79.
  * [10] Damanhuri NS, Chiew YS, Othman NA et al.  Assessing respiratory mechanics using pressure reconstruction method in mechanically ventilated spontaneous breathing patient.  _Computer Methods and Programs in Biomedicine_ 2016; 130: 175–185.
  * [11] Chase JG, Preiser JC, Dickson JL et al.  Next-generation, personalised, model-based critical care medicine: a state-of-the art review of in silico virtual patient models, methods, and cohorts, and how to validation them.  _Biomedical Engineering Online_ 2018; 17(1): 1–29.
  * [12] Morton SE, Knopp JL, Chase JG et al.  Optimising mechanical ventilation through model-based methods and automation.  _Annual Reviews in Control_ 2019; 48: 369–382.
  * [13] Kim KT, Morton S, Howe S et al.  Model-based PEEP titration versus standard practice in mechanical ventilation: a randomised controlled trial.  _Trials_ 2020; 21(1): 1–18.
  * [14] Schranz C, Docherty PD, Chiew YS et al.  Iterative integral parameter identification of a respiratory mechanics model.  _Biomedical Engineering Online_ 2012; 11(1): 1–14.
  * [15] Lucangelo U, Bernabé F and Blanch L.  Respiratory mechanics derived from signals in the ventilator circuit.  _Respiratory care_ 2005; 50(1): 55–67.
  * [16] Grinnan DC and Truwit JD.  Clinical review: respiratory mechanics in spontaneous and assisted ventilation.  _Critical Care_ 2005; 9(5): 1–13.
  * [17] West JB.  _Respiratory Physiology: The Essentials_.  Lippincott Williams & Wilkins, 2012.
  * [18] Sahetya SK, Hager DN, Stephens RS et al.  PEEP titration to minimize driving pressure in subjects with ARDS: a prospective physiological study.  _Respiratory Care_ 2020; 65(5): 583–589.
  * [19] van der Zee P and Gommers D.  Recruitment maneuvers and higher PEEP, the so-called open lung concept, in patients with ARDS.  _Critical Care_ 2019; 23(1): 1–7.
  * [20] Dries DJ and Marini JJ.  Finding best PEEP: A little at a time.  _Respiratory Care_ 2020; 65(5): 722–724.
  * [21] Tobin MJ.  Mechanical ventilation.  _New England Journal of Medicine_ 1994; 330(15): 1056–1061.
  * [22] Rimensberger PC, Pristine G, Mullen JBM et al.  Lung recruitment during small tidal volume ventilation allows minimal positive end-expiratory pressure without augmenting lung injury.  _Critical Care Medicine_ 1999; 27(9): 1940–1945.
  * [23] Harris RS.  Pressure-volume curves of the respiratory system.  _Respiratory Care_ 2005; 50(1): 78–99.
  * [24] Lu Q, Vieira SR, Richecoeur J et al.  A simple automated method for measuring pressure–volume curves during mechanical ventilation.  _American Journal of Respiratory and Critical Care Medicine_ 1999; 159(1): 275–282.
  * [25] Escolar Castellón JdD et al.  Lung hysteresis: a morphological view.  _Histology and Histopathology_ 2004; 19(1): 159–166.
  * [26] Fedullo AJ, Jung-Legg Y, Snider GL et al.  Hysteresis ratio: a measure of the mechanical efficiency of fibrotic and emphysematous hamster lung tissue.  _American Review of Respiratory Disease_ 1980; 122(1): 47–52.
  * [27] Venegas JG, Harris RS and Simon BA.  A comprehensive equation for the pulmonary pressure-volume curve.  _Journal of Applied Physiology_ 1998; 84(1): 389–395.
  * [28] Tsuchida S, Engelberts D, Peltekova V et al.  Atelectasis causes alveolar injury in nonatelectatic lung regions.  _American Journal of Respiratory and Critical Care Medicine_ 2006; 174(3): 279–289.
  * [29] González-López A, García-Prieto E, Batalla-Solís E et al.  Lung strain and biological response in mechanically ventilated patients.  _Intensive Care Medicine_ 2012; 38(2): 240–247.
  * [30] Muñiz-Albaiceta G.  Pressure-volume curves in acute pulmonary injury.  _Medicina Intensiva_ 2009; 33(5): 243–250.
  * [31] Bates JH.  _Lung Mechanics: an Inverse Modeling Approach_.  Cambridge University Press, 2009.
  * [32] Li Y, O’Neill Z, Zhang L et al.  Grey-box modeling and application for building energy simulations - A critical review.  _Renewable and Sustainable Energy Reviews_ 2021; 146: 111174.
  * [33] Estrada-Flores S, Merts I, De Ketelaere B et al.  Development and validation of “grey-box” models for refrigeration applications: A review of key concepts.  _International Journal of Refrigeration_ 2006; 29(6): 931–946.
  * [34] Gavin HP.  The Levenberg-Marquardt algorithm for nonlinear least squares curve-fitting problems.  _Department of Civil and Environmental Engineering, Duke University_ 2019; 19.
  * [35] Lourakis MI et al.  A brief description of the Levenberg-Marquardt algorithm implemented by levmar.  _Foundation of Research and Technology_ 2005; 4(1): 1–6.
  * [36] Widrow B, Kollar I and Liu MC.  Statistical theory of quantization.  _IEEE Transactions on Instrumentation and Measurement_ 1996; 45(2): 353–361.
  * [37] Ranieri VM, Zhang H, Mascia L et al.  Pressure–time curve predicts minimally injurious ventilatory strategy in an isolated rat lung model.  _The Journal of the American Society of Anesthesiologists_ 2000; 93(5): 1320–1328.
  * [38] Hess DR.  Recruitment maneuvers and PEEP titration.  _Respiratory Care_ 2015; 60(11): 1688–1704.
  * [39] Costa EL, Borges JB, Melo A et al.  Bedside estimation of recruitable alveolar collapse and hyperdistension by electrical impedance tomography.  In _Applied Physiology in Intensive Care Medicine 1_. Springer, 2012\.  pp. 165–170.
  * [40] Borges JB, Okamoto VN, Matos GF et al.  Reversibility of lung collapse and hypoxemia in early acute respiratory distress syndrome.  _American Journal of Respiratory and Critical Care Medicine_ 2006; 174(3): 268–278.