An enhanced respiratory mechanics model based on double-exponential and fractional calculus

We address mathematical modelling of respiratory mechanics and put forward a model based on double-exponential and fractional calculus for parameter estimation, model simulation, and evaluation based on actual data. Our model has been implemented on a publicly available executable code with adjustable parameters, making it suitable for different applications. Our analysis represents the first application of fractional calculus and double-exponential modelling to respiratory mechanics, and allows us to propose a hybrid model fitting experimental data in different ventilation modes. Furthermore, our model can be used to study the mechanical features of the respiratory system, improve the safety of ventilation techniques, reduce ventilation damages, and provide strong support for fast and adaptive determination of ventilation parameters.

1 Introduction

With the aging population and the increase in chronic respiratory diseases (Divo et al., 2014; Soriano et al., 2020; Xie et al., 2020), the demand for mechanical ventilation treatment has largely increased (Wunsch, 2020). At the same time, the outbreak of COVID-19 in recent years has further boosted the demand for ventilators. As a consequence, the analysis of more accurate respiratory mechanics models has become extremely important, as it can be used to predict and optimize parameter settings in mechanical ventilation therapy. Further, this kind of models are required to dynamically adjust parameters according to the respiratory characteristics of different patients, in order to improve treatment outcomes and reduce ventilation injuries.

In the field of mechanical ventilation, the respiratory system is often modelled as a linear dynamical system described by linear differential equations. Examples are provided by the single-compartment model (Bates, 2009; Ghafarian et al., 2016; Gertler, 2021) corresponding to a first-order linear differential system, the two-compartment model (Bates, 2009; Carvalho and Zin, 2011; Ghafarian et al., 2016) corresponding to a second-order differential system, and the viscoelastic model (Ganzert et al., 2009). Over the last few years, different types of respiratory mechanics models have been put forward, including nonlinear models involving elastance and resistance of the respiratory system (Chiew et al., 2015; Redmond et al., 2017; Ellwein Fix et al., 2018; Marconi and De Lazzari, 2020). These studies generally provide a good theoretical modelling and may be used to simulate the respiratory system. However, they generally involve idealized elementary function description of ventilation pressure (Ellwein Fix et al., 2018; Marconi and De Lazzari, 2020; Dincel, 2021), and ignore the inherent randomness and background noise/error of real respiratory processes. Other studies are based on experimental data, but produce results in disagreement with actual data. For example, the simulations of Albanese et al. (2016), lead to a minimum CV(RMSE) (stands for the coefficient of variations of the root mean square errors) volume curve equal to 0.079, which is an order of magnitude higher than the simulation results of all models in this paper. At the same time, these studies did not provide executable program code and corresponding original data.

In this paper, we propose an accurate respiratory mechanics model based on previous research and studies based on actual measurement data. We consider Bama pigs for actual measurements, since they provide conditions close to the healthy state of the human respiratory system, and assess results of our models against actual measurement data. By combining fractional calculus with a double-exponential model, we are able to describe the dynamical behaviour of the coefficients in the differential equations, and to put forward a hybrid model consistent with the experimental in wide range of ventilation modes, and provided. We also provide a public executable code package for parameter estimation, simulation methods, and optimized parameter search.

2 Methods

2.1 Respiratory mechanics models

2.1.1 Linear models

First order linear differential models with constant coefficients (linear single-compartment model). According to previous studies on respiratory mechanics models (Bates, 2009; Carvalho and Zin, 2011; Hess, 2014; Ghafarian et al., 2016; Gertler, 2021), the classical linear single-compartment model corresponds to the following differential equation:

$P\left(t\right)=R\frac{dV\left(t\right)}{dt}+EV\left(t\right)+P_0(1)$

for the pressure (applied power) to the lung P(t), where V(t) is the air volume entering the lung (Tidal volume), E and R denote the elastance and resistance of respiratory system respectively. The offset pressure P₀ represents the positive end-expiratory pressure (PEEP). The compliance C of respiratory system is the reciprocal of elastance E, that is C = 1/E.

The electrical circuit that represents the linear single-compartment model consists of a resistance R, a capacitor C (1/E) and a generator P(t) (Figure 1A). The corresponding impedance is (Carvalho and Zin, 2011; Li, 2011):

$Z=R-\frac{j}{\omega C}(2)$

where $\omega$ is the angular frequency in radians/second, and the symbol j represents the unit imaginary number.

FIGURE 1

FIGURE 1. Electrical analogue of the respiratory system. (A) Scheme corresponding to Eq. 1, (B) scheme corresponding to Eq. 10, (C) scheme corresponding to Eqs 14, (D) different impedance elements in series and in parallel corresponding to Eq. 16. The resistance R of respiratory system in Eq. 1 is equivalent to the resistance in the circuit, and the compliance C of respiratory system is equivalent to the capacitor in the circuit. The current in the circuit represents the flow in the airway.

Moreover, based on the electrical analogous (a), if we assume the initial value V (0) = 0 and P₀ = 0 (or define a new pressure ${\widehat{P}\left(t\right)=P\left(t\right)-P}_0$) for simplicity, we can obtain the open-loop transfer function between the pressure P(t) and volume V(t):

$G\left(s\right)=\frac{V\left(s\right)}{P\left(s\right)}=\frac{C}{CRs+1}=\frac{1}{Rs+E}(3)$

The inverse Laplace transform of transfer function G(s) is given by:

${G\left(t\right)=L}^{-1}\left(G\left(s\right)\right)=L^{-1}\left(\frac{\frac{1}{R}}{s+\frac{1}{CR}}\right)=\frac{1}{R}{e}^{-\frac{t}{CR}}=\frac{1}{R}{e}^{-\frac{tE}{R}}(4)$

If P(t) is a unit step signal, this will be the solution in time domain.

If we express the pressure P(t) as a linear combination of elementary functions:

$P\left(t\right)={\displaystyle \sum _{i=0}^n}{k_{i}P}_i\left(t\right)(5)$

Such as sin (at), $e^{-at}$, or polynomial functions, and the k_i are constant coefficients, then we can theoretically obtain the general solution of V(t) in time domain as follows

$V\left(t\right)=L^{-1}\left(G\left(s\right)P\left(s\right)\right)=L^{-1}\left(G\left(s\right)\right)*L^{-1}\left(P\left(s\right)\right)={\int }_0^{t}G\left(t-\tau \right)P\left(\tau \right)d\tau (6)$

where “*” denotes convolution. The result is the same as that obtained by directly solving the equation by Laplace transform. Based on the superposition principle of linear differential equation, the final solution of V(t) is given by

$V\left(t\right)={\displaystyle \sum _{i=0}^n}{\int }_0^t{k}_{i}G\left(t-\tau \right)P_i\left(\tau \right)d\tau (7)$

For example, if ${\mathrm{P}\left(\mathrm{t}\right)=\mathrm{P}}_{constant}\sin \left(\text{at}\right),\mathrm{t}\in \left[0,\mathrm{\pi }\right]$ , Eq. 1 can be written as: $R\frac{dV\left(t\right)}{dt}+EV\left(t\right)={\mathrm{P}}_{constant}\sin \left(\text{at}\right)-P_0$ , with general solution $V\left(t\right)=K{e}^{-\frac{t}{CR}}+\frac{PC}{a^2{C}^2{R}^2+1}\left(\mathit{sin}\left(at\right)-aCRcos\left(at\right)\right)-C{P}_0$, where K is a constant determined by the initial conditions and P_constant is also a constant that determines the magnitude of pressure. Other simple examples are illustrated in the Supplementary Material Section S1.

Although it is theoretically possible to represent any P(t) based on the Fourier transform, this method is not particularly useful due to the computational complexity and the lack of accuracy.

Second order models with constant coefficients. Assuming that the airway is an approximately cylindrical structure of length l and radius r and using Newton’s second law (Carvalho and Zin, 2011; Hess, 2014) we have:

$\pi r^2{P}_i=P_{i}S=F=ma=\pi r^{2}l\rho \left(\frac{d^{2}V}{d{t}^2}\right)(8)$

that is:

$P_i=l\rho \frac{d^{2}V}{d{t}^2}=I\frac{d^{2}V}{d{t}^2}(9)$

where ρ is the gas density, and $I=l\rho$ . Due to the relatively low gas density, the inertia coefficient I (Inertia) is generally very small compared to other parameters with the same dimensions in respiratory mechanics models. Taking airway resistance and elasticity into account, the pressure P and volume V are related by the following Eq. 10.

$P\left(t\right)=I\frac{d^{2}V\left(t\right)}{d{t}^2}+R\frac{dV\left(t\right)}{dt}+EV\left(t\right)+P_0(10)$

This is a second order linear differential equation, which corresponds to the electrical analogue shown in Figure 1B with impedance:

$Z=R+j\omega I-\frac{j}{\omega C}(11)$

where I represents the total inductance of the circuit.

If we assume the initial conditions V (0) = 0 and V’ (0) = 0 (the first-order derivative of V) and P₀ = 0 (or define new pressure ${\widehat{P}\left(t\right)=P\left(t\right)-P}_0$) for simplicity, we can also obtain the transfer function, which can be written as the sum of two fractions:

$G\left(s\right)=\frac{V\left(s\right)}{P\left(s\right)}=\frac{1}{I{s}^2+Rs+E}=\frac{c_1}{s+a}+\frac{c_2}{s+b}(12)$

where a, b, c₁, c₂ are constant and the real part of a and b are positive, that is Re(a) > 0, Re(b) > 0 (satisfying the stability constraint conditions) (Zhong, 1988; Hu, 2007).

The inverse Laplace transform of the transfer function in Eq. 12 is given by:

$G\left(t\right)=L^{-1}\left(G\left(s\right)\right)=L^{-1}\left(\frac{c_1}{s+a}+\frac{c_2}{s+b}\right)=c_1{e}^{-at}+c_2{e}^{-bt}(13)$

Assuming that P(t) is a step signal, the general solution of Eq. 10 has the same form as Eq. 14. On the other hand, if we have single-compartment models connected in series or in parallel, we obtain the two-compartment model proposed by Bates (Bates, 2009), whose dynamics is similar to that governed by Eq. 10 (Figure 1C), i.e.,

$P\left(t\right)+k\frac{dP\left(t\right)}{dt}=I\frac{d^{2}V\left(t\right)}{d{t}^2}+R\frac{dV\left(t\right)}{dt}+EV\left(t\right)(14)$

where k is a constant. The transfer function of Eq. 14 is similar to that reported in Eq. 12:

$G\left(s\right)=\frac{V\left(s\right)}{P\left(s\right)}=\frac{1+ks}{I{s}^2+Rs+E}=\frac{k_1}{s+a}+\frac{k_2}{s+b}(15)$

In other words, Eqs 10, 14 has the same solutions.

According to the superposition principle, we can also obtain the general solution of second order linear differential models based on Eq. 7. Simple examples are given in the Supplementary Material Section 2.

High-order linear differential models. The term $\mathit{I}\frac{{\mathit{d}}^{\mathbf{2}}\mathit{V}\left(\mathit{t}\right)}{\mathit{d}{\mathit{t}}^{\mathbf{2}}}$ in Eq. 10 corresponds to an acceleration, however, in a respiratory mechanical model (Carvalho and Zin, 2011; Eager et al., 2016; Dincel, 2021) we may also consider higher derivatives (like jerk and snap), such that more general models may be obtained, in the following general form:

$b_{0}P\left(t\right){+b}_1\frac{dP\left(t\right)}{dt}+\dots +b_m\frac{d^{m}P\left(t\right)}{d{t}^m}=a_{0}V\left(t\right)+a_1\frac{dV\left(t\right)}{dt}+a_2\frac{d^{2}V\left(t\right)}{d{t}^2}+\dots +a_n\frac{d^{n}V\left(t\right)}{d{t}^n}(16)$

where n > m.

Assuming the initial values V (0), V’ (0), V’’ (0), … = 0, we obtain the transfer function as

$G\left(s\right)=\frac{V\left(s\right)}{P\left(s\right)}=\frac{b_m{s}^m+\dots +b_2{s}^2+b_{1}s+b_0}{a_n{s}^n+\dots +a_2{s}^2+a_{1}s+a_0}(17)$

In a stable respiratory system, the system should satisfy Routh-Hurwitz criterion. The inverse Laplace transform of transfer function G(s) is as follows:

${G\left(t\right)=L}^{-1}\left(G\left(s\right)\right)={\displaystyle \sum _{i=0}^n}{k_{i}e}^{{\mathrm{\lambda }}_{i}t}(18)$

where ${\text{Re}\left(\mathrm{\lambda }\right.}_i\left)<0\right.$.

If P(t) is a unit step signal, this represents the solution of Eq. 16.

Using again the electrical analogy, Eq. 16 corresponds to put more components in series or in parallel (in a circuit, components connected in parallel are equivalent to adding their transfer functions in the s domain, while components connected in series are equivalent to multiplying their transfer functions in the s domain). The overall circuit has an equivalent impedance Z_total:

$Z_{total}=R_{total}+j\omega I_{total}-\frac{j}{\omega C_{total}}(19)$

where R_total, I_total, C_total represent the overall equivalent resistance, capacitor, inductance respectively in the circuit (Figure 1D). Eq. 19 and Eq. 11 have indeed a similar form.

However, high-order linear differential models are rarely encountered and applied in practical applications. We suggest caution in using high-order models to describe the respiratory system, primarily based on the following reasons:

1. Since Eq. 19 and Eq. 11 have a similar form, we can approximate the effect of high-order differential models by using first and second order differential models.

2. Higher-order differential models are difficult to stabilize in practical applications, as there are always disturbances present in real-world situations. We illustrate some examples in the Supplementary Material Section 3.

3. There are many parameters in high-order differential models, which typically require more measured data to be estimated. In the following sections, we introduce fractional calculus models with fewer parameters and similar accuracy.

4. In nondegenerate mechanics systems, (the Hessian matrix is full rank), based the Hamilton’s Principle and Ostrogradsky’s theorem, terms of order three and higher in higher-order differential models do not exist. More details on thus point may be found in the Supplementary Material Section 4.

Overall, given the above analysis, we are not going to employ high-order differential models in processing measured data.

2.1.2 Fractional calculus models

If there is uncertainty about the appropriate order of the differential terms in the respiratory mechanics model and on the actual weight of each order, we may consider introducing fractional-order terms involving different orders of derivatives into the model. This would also allow us to take into account the memory-dependent viscoelasticity and heterogeneity features of the lung tissue (Bates, 2009; Ganzert et al., 2009; Chen et al., 2010; Carvalho and Zin, 2011; Ionescu and Kelly, 2017; Xue and Bai, 2023), which cannot be captured by classical integer-order differential equations-based lumped-parameter models. Here, we propose a fractional calculus model to solve this problem, which has fewer parameters to be estimated compared to classical models.

Fractional calculus may be introduced in different and slightly inequivalent ways, involving Grünwald-Letnikov, Riemann–Liouville, and Caputo fractional derivatives. Here, we primarily consider the Grünwald-Letnikov definition that is convenient for numerical calculations. The α-order derivative of a given function f(t) is expressed as (Xue and Bai, 2023):

${\mathcal{D}}_{t_0}^{\alpha }f\left(t\right)=\underset{h\to 0}{\lim }\frac{1}{h^{\alpha }}{\displaystyle \sum _{j=0}^{\frac{t-t_0}{h}}}\frac{{\left(-1\right)}^j\mathrm{\Gamma }\left(\mathrm{\alpha }+1\right)}{\mathrm{\Gamma }\left(\mathrm{j}+1\right)\mathrm{\Gamma }\left(\mathrm{\alpha }-\mathrm{j}+1\right)}f\left(t-jh\right)(20)$

In the above formula, h is a sufficiently small time step, and $\mathrm{\Gamma }$ denotes the Gamma function, that is $\mathrm{\Gamma }\left(\mathrm{z}\right)={\int }_0^{\infty }{e}^{-t}{t}^{z-1}dt$。This definition assumes that the value of the function f(t) is 0 when t ≤ t₀, in agreement with the assumption that any physical quantity is zero before the initial time t0. Fractional order derivatives are able to capture memory properties (applicable to differentiation where α > 0 and to integration where α < 0). More about numerical calculations of fractional order calculus can be found in the Supplementary Material Section 5.

Furthermore, according to Eq. 20, any fractional order derivative term also involves some higher order derivative features.

Based on the fractional order differentiation proposed above, we modify Eqs 1, 20, thus obtaining the following fractional order derivative respiratory mechanics model:

$P=EV\left(t\right)+R{V}^{\prime }\left(t\right)+a{\mathcal{D}}_{t_0}^{\alpha }V\left(t\right)+P_0(21)$

where ${\mathcal{D}}_{t_0}^{\alpha }V\left(t\right)$ is the α-order derivative of the $V\left(t\right)$, with a fractional order differential term α > 0. When α = 1, it takes the form of Eq. 10. In the subsequent sections, we mainly focus on cases with 0<α < 1.

The above formula contains a single fractional order derivative term. Analogously, we may also introduce models containing more fractional terms, such as the generalization of Eq. 22 containing 2 fractional order derivative terms, i.e.,

$P=EV\left(t\right)+R{V}^{\prime }\left(t\right)+a{\mathcal{D}}_{t_0}^{\alpha }V\left(t\right)+b{\mathcal{D}}_{t_0}^{\beta }V\left(t\right)+P_0(22)$

As a matter of fact, upon considering actual measurement data, one sees that more fractional order terms do not provide better results than a single one. Parameters estimation and model simulation results for Eq. 21, 22 are reported in the Supplementary Material Section 6.

2.1.3 Non-linear models

The models mentioned earlier are all cases where the coefficients of the differential equations are constants. In the actual respiratory system, the elasticity E and airway resistance R may be variables rather than constants (Van Drunen et al., 2014; Chiew et al., 2015; Redmond et al., 2017; Ellwein Fix et al., 2018; Marconi and De Lazzari, 2020). When the elasticity E changes with V(t), and the resistance R changes with V′(t), based on Eq. 1, pressure P(t) and volume V(t) satisfy the following equation:

$P\left(t\right)=Ef\left(V\left(t\right)\right)V\left(t\right)+Rh\left(V^{\prime }\left(t\right)\right)V^{\prime }\left(t\right)+P_0(23)$

where V′(t) represents the first-order derivative of V with respect to t.

Specifically, when only the elasticity E varies with V(t) or only the resistance R varies with V′(t), Eq. 23 can be evolved into Eq. 24 and Eq. 25.

$P\left(t\right)=Ef\left(V\left(t\right)\right)V\left(t\right)+R{V}^{\prime }\left(t\right)+P_0(24)$

$P\left(t\right)=EV\left(t\right)+Rh\left(V^{\prime }\left(t\right)\right)V^{\prime }\left(t\right)+P_0(25)$

Nonlinear models containing polynomial terms. In order to describe physical situations in which the elasticity E is monotonically (and strongly) positively correlated with V(t) and the resistance R is monotonically strongly positively correlated with V′(t), we should consider polynomial functions f and h in Eq. 23. To a first approximation, we may safely ignore the different inflection points of the polynomials and focus on the highest order term that mostly affects the behaviour of the functions. This amounts to consider f of the form of ${\mathit{E}}_{\mathit{o}\mathit{f}\mathit{f}\mathit{s}\mathit{e}\mathit{t}}+\mathit{\alpha }{\mathit{V}\left(\mathit{t}\right)}^{\mathit{\beta }}$, where E_offset > 0, and the constants α and β capture the behaviour of the elasticity E as a function of V(t). The form of the function h is analogue, i.e., ${\mathit{R}}_{\mathit{o}\mathit{f}\mathit{f}\mathit{s}\mathit{e}\mathit{t}}+\mathit{\alpha }{{\mathit{V}}^{\prime }\left(\mathit{t}\right)}^{\mathit{\beta }}$. Overall, we have

$P\left(t\right)=E_0\left(E_{offset}+{\alpha }_E{V\left(t\right)}^{{\beta }_E}\right)V\left(t\right)+R_0\left(R_{offset}+{\alpha }_R{V^{\prime \left(t\right)}}^{{\beta }_R}\right)V^{\prime }\left(t\right)+P_0(26)$

Then Eq. 23 can be written as

$P\left(t\right)=EV\left(t\right)+R{V}^{\prime }\left(t\right)+a{V\left(t\right)}^x+b{V^{\prime \left(t\right)}}^y+P_0(27)$

where a, b, x and y are constant.

Similarly, when only the elasticity E changes with V(t) in a polynomial form or only the resistance R changes with V′(t) in a polynomial form, Eq. 27 can be written as

$P\left(t\right)=EV\left(t\right)+R{V}^{\prime }\left(t\right)+a{V\left(t\right)}^x+P_0(28)$

$P\left(t\right)=EV\left(t\right)+R{V}^{\prime }\left(t\right)+b{V^{\prime \left(t\right)}}^y+P_0(29)$

As a matter of fact, there exists some nonlinear models with analytic solutions in terms of elementary functions. For instance, for the model $\left(t\right)=E_0\left(E_{offset}+kV\left(t\right)\right)V\left(t\right)+R{V}^{\prime }\left(t\right)+P_0=EV\left(t\right)+R{V}^{\prime }\left(t\right)+a{V\left(t\right)}^2+P_0$ , if we assume the P(t) is a step signal with constant value P_c, then its analytic general solution is $V\left(t\right)=\frac{-\mathrm{E}+\sqrt{-{\mathrm{E}}^2+4a\left({\mathrm{P}}_0-{\mathrm{P}}_c\right)}\tan \left[\frac{\sqrt{-{\mathrm{E}}^2+4a\left({\mathrm{P}}_0-{\mathrm{P}}_c\right)}\left(-t+kR\right)}{2R}\right]}{2a}$, where k is the constant coefficient determined by the initial conditions.

Non-linear models containing exponential terms. In some cases, the respiratory system elasticity E may show a step response property^10,12. Combined with the results presented in the following sections, this corresponds to the model described by Eq. 24. If function f is an exponential function, double-exponential function and triple exponential function, P and V satisfy one of the following equations.

$P\left(t\right)=EV\left(t\right)+a{e}^{xV\left(t\right)}V\left(t\right)+R{V}^{\prime }\left(t\right)+P_0(30)$

$P\left(t\right)=EV\left(t\right)+a{e^e}^{xV\left(t\right)}V\left(t\right)+R{V}^{\prime }\left(t\right)+P_0(31)$

$P\left(t\right)=EV\left(t\right)+a{e^{e^e}}^{xV\left(t\right)}V\left(t\right)+R{V}^{\prime }\left(t\right)+P_0(32)$

The growth in single exponential function is fast, but not as fast as the factorial function. The double-exponential function grows faster than the factorial function (see the Supplementary Material Section 7). The growth rate of the Ackermann’s function is much faster than the double-exponential function. As it may intuitively be expected, we will show that without using piecewise functions the faster is the growth, the better is the fitting of the measured data (for a step response).

2.1.4 Hybrid model

Considering the step-like features of the respiratory system elasticity, as well as the memory-dependent viscoelasticity and heterogeneity characteristics of lung tissue, we may combine Eq. 21 with Eq. 30 to obtain the hybrid model of Eq. 33, which itself incorporates different properties:

$P\left(t\right)=EV\left(t\right)+R{V}^{\prime }\left(t\right)+a{\mathcal{D}}_{t_0}^{\alpha }V\left(t\right)+\mathrm{b}{e^e}^{\beta V\left(t\right)}+P_0(33)$

As we will see in the following Sections, this hybrid model provides the best fit for measured data in most situations.

2.2 Data collection

We use a Dräger Savina 300 to measure airway pressure P(t), flow speed V’(t) and respiratory volume V(t) in the actual respiration of three Bama pigs. The first Bama pig (sample ID as case1) weighs 28 kg, and data are taken for three ventilation modes, namely, pressure control - assist control (PC-AC), volume control - assist control (VC-AC), and volume control - synchronized intermittent mandatory ventilation (VC-SIMV). The second one (sample ID as case2) weighs 31 kg and the third one (sample ID as case3) weighs 24.3 kg, and data are taken for the VC-AC ventilation mode. In all measurements, pressure P(t), flow speed V′(t), and respiratory volume V(t) data are collected every 10 ms for at least 6 min. In the pressure control ventilation mode, set peak breathing pressure = 16 mbar, plateau pressure = 15 mbar and PEEP = 3 mbar. In the volume control ventilation mode, the tidal volume is set to 0.3 L.

During the measurement process, the original HL7 data are first obtained through serial communication, and then the pressure P, flow speed V′, and respiratory volume V data at each moment are further parsed and extracted. The data at the same moment are placed in the same row, with each row divided into three columns: pressure P, flow speed V’, and respiratory volume V, separated by tabs. The data are arranged in chronological order. Finally, the unit of P is converted to mbar, the flow rate is in L/s, and the respiratory volume is in L, which are written into the txt file.

For cross-validation at different time periods, two data files of 150 s duration (15,000 data points, about 50 breathing cycles) are extracted from each breathing mode. After sorting, the data files are respectively named as case1PCAC1.txt, case1VCAC1.txt, case1VCSIMV1.txt, case2VCAC1.txt, case3VCAC1.txt and case1PCAC2.txt, case1VCAC2.txt, case1VCSIMV2.txt, case2VCAC2.txt case3VCAC2.txt.

2.3 Parameter estimation

Parameter estimation for all the above-mentioned ordinary differential equations with constant coefficients can be performed using the following procedure (Bates, 2009; Redmond et al., 2017; Avilés-Rojas and Hurtado, 2022). The first observation is that all the models analyzed in this paper may be recasted as

$P\left(t\right)=a_1{f}_1\left(t\right)+a_2{f}_2\left(t\right)+\dots +a_m{f}_m\left(t\right)+P_0(34)$

where $\left\{f_1\left(t\right),f_2\left(t\right)\dots f_n\left(t\right)\right\}$ denotes any linear or nonlinear term of the differential equation, for example, $f_n\left(t\right)$ may represent ${\left(V^{\prime }\left(t\right)\right)}^2$, ${\mathcal{D}}_{t_0}^{\alpha }V\left(t\right)$, $e^{xV\left(t\right)}$, etc., while $\overrightarrow{A}=\left\{a_1,a_2\dots a_{\mathrm{n}}\right\}$ is the corresponding set of coefficients, i.e., the quantities to be estimated.

The values of pressure P(t), flow speed V′(t), and respiratory volume V(t) collected in real time are organized in vector and matrix form and used as input for the algorithm. If the collected values of the respiratory system are at n time points $\left\{t_1,t_2\dots t_{\mathrm{n}}\right\}$, then the input data are the n-dimensional pressure vector P(t), and the data matrix X with n rows and m+1 columns:

$\overrightarrow{P}=\left\{\begin{array}{c}\left.{P(t}_1\right)\\ \left.{P(t}_2\right)\\ \begin{array}{c}⋮\\ \left.{P(t}_n\right)\end{array}\end{array}\right\}(35)$

$X=\left\{\begin{array}{c}\begin{array}{ccc}{f}_1\left(t_1\right)& \dots & \begin{array}{cc}{f}_m\left(t_1\right)& 1\end{array}\end{array}\\ \begin{array}{ccc}{f}_1\left(t_2\right)& \dots & \begin{array}{cc}{f}_m\left(t_2\right)& 1\end{array}\end{array}\\ \begin{array}{c}⋮\\ \begin{array}{ccc}{f}_1\left(t_{\mathrm{n}}\right)& \dots & \begin{array}{cc}{f}_m\left(t_{\mathrm{n}}\right)& 1\end{array}\end{array}\end{array}\end{array}\right\}(36)$

Then, the coefficients can be estimated from

$\overrightarrow{A}={\left[X^{T}X\right]}^{-1}{X}^T\overrightarrow{P}(37)$

In the above Eq. 37, $X^T$ represents matrix transposition, and ${\left[*\right]}^{-1}$ represents matrix inversion. We estimate parameters by minimizing the sum of squared residuals (SSR). At the extreme point, the partial derivatives with respect to each coefficient are zero, and this provides constraints to obtain their values. In order to determine the fractional-order in Eq. 21, we employ Particle Swarm Optimization (PSO) algorithm in machine learning (Kennedy and Eberhart, 1995; Gad, 2022; Shami et al., 2022; Xie et al., 2022; Yang et al., 2022), searching for the stable optimal solution with the lowest overall SSR. In this way, we also find the coefficient x for the exponential Eqs 30–32, and the hybrid Eq. 33. For all the models, suitable values of the number of particles N and iterations M may be found, which ensure stable optimal solutions. In general, for Eq. 21, M = N = 5 is enough, whereas for the exponential model we need M = N = 10, and for the hybrid model, M = N = 20. Once the values of α and x are determined, their corresponding constant coefficients can be calculated according to Eq. 37.

To ensure computational accuracy, we use the different packages in Python to calculate the fractional order derivatives (Adams, 2019), and call Python code within MATLAB while optimizing the SSR value of the models. In addition, since the exponential Eqs 30–32, is prone to produce a singular matrix, we perform address Eq. 37 using a pseudo-inverse matrix to reduce the distortion and deviation of numerical results.

More details about faster empirical estimation methods of the parameters for fractional calculus and exponential terms may be found in the Supplementary Material Section 8.

2.4 Numerical simulations

After obtaining the estimated values of the parameters, the stability of the models is assessed using the Routh-Hurwitz stability criterion, and then the subsequent simulations are carried out carefully considering their actual physical meaning (Hu, 2007; Wang et al., 2007). Generally, we accept estimated values of the parameters only if they are positive. On the one hand, this is required by the Routh-Hurwitz stability criterion, dictating that our models are stable only if all coefficients are positive. On the other hand, from a practical physical perspective, most of the variables of interest, e.g., the elasticity E of lung tissue are positive numbers.

After the parameters meet the above conditions, we use MATLAB/Simulink and Python to complete the calculations and construct the Simulink model for all the models considered in this work. The open-source vfoderiv3 module of Simulink (Scherer et al., 2011; Sierociuk et al., 2015; Sierociuk, 2023) is used for the simulation of fractional calculus. The Simulink model for the other models may be found in the Supplementary Material Section 9.

2.5 Evaluation criteria

The comparison of models with the same number of parameters to be estimated is made in terms of the SSR and the root mean square error (RMSE). Models with lower values of these two indicators should be preferred. The definition is given by

$\text{SSR}={\displaystyle \sum _{i=1}^n}{\left(V_i-\widehat{V_i}\right)}^2(38)$

$\text{RMSE}=\sqrt{\frac{SSR}{n}}(39)$

where $V_i$ represents the actual respiratory volume value from measured data at time i (the i th data point), and $\widehat{V_i}$ represents the value from model simulations. The actual measurement period contains a total of n data points (an overall duration of 150 s corresponds to n = 15,000). The comparison between different models is made using the Bayesian Information Criterion (BIC) which considers the number of parameters. If $\text{residuals}\sim N\left(0,{\sigma }^2\right)$, we ignore the irrelevant constants and use the equivalent BIC value, i.e., Eq. 40, for model comparison and evaluation (Sheng et al., 2020).

$\text{BIC}=\mathrm{n}*\ln \left(\text{SSR}\right)+\mathrm{k}*\ln \left(\mathrm{n}\right)(40)$

In the above formula, k represents the number of model parameters. The process of deriving Eq. 40 based on the original BIC definition and the assumption of residual distribution may be found in the Supplementary Material Section 10.

In addition, we also use the Pearson correlation coefficient between the measured and simulated values of the respiratory volume as an additional criterion for optimization (the closer the correlation coefficient is to 1, the better).

2.6 Computation details

The executable code and actual respiratory data can be found in the Supplementary Material and our Gitee repository (https://gitee.com/lizwtest/mechanical_ventilation).

For models that do not contain exponential or fractional order differential terms, such as Eq 1 and Eq. 2 and Eq. 10, and (27–29), the measured respiratory mechanical data of the respiratory system can be substituted into Eq. 37 based on Section 2.3. This results in calculated values with practical physical characteristics such as respiratory system elasticity E, airway resistance R, P0 (PEEP), and so on.

For the calculation of models containing double exponential or fractional order differential terms, taking the double exponential model of Eq. 31 as an example, when a suitable value of the exponential parameter x is given in the model, the respiratory system mechanics data is substituted into Eq. 37 to obtain the elasticity “E,” airway resistance “R” and coefficient “a” which can reflect the weight of the step mechanical characteristic. At this time, the SSR value can be obtained by Eq. 38 as the score of the PSO function based on the predicted respiratory volume $\widehat{V}$ (t) and the measured volume V(t), and then the velocity vector of N particles is continuously updated for M rounds to find a better value of x (corresponding to the position of the particle), thus obtaining a lower SSR value corresponding to the optimized parameter x. The main flow of the proposed Algorithm 1 is shown in the pseudo-code. The NCP in the pseudo-code stands for non-coefficient parameters. P_SSR_best, G_SSR_best, P_best, G_best, represent the particle best SSR value, the global best SSR value, the particle best location, and the global best location, respectively.

Algorithm 1. Pseudo-code for models with PSO algorithm.

input:

Data -- Respiratory mechanical data P(t), V’ (t), and V(t)

D -- Number of NCPs

M -- Max number of iterations

N -- Swarm size

LB -- Lower boundary of the search space

UB -- Upper boundary of the search space

output:

G_best -- the best positions (NCPs) found so far

(1) Initialization;

(2) for each particle i = 1 to N do

(3) for each NCP j = 1 to D do

(4) Initialize particle velocity and position randomly within [LB, UB]^D;

(5) end for

(6) Substituting Data into Eq. 37 to determine all the model parameters;

(7) Compute predicted respiratory volume $\widehat{V}$ (t) based on respiratory mechanical model;

(8) Substituting $\widehat{V}$ (t) and V(t) to Eq. 38 to initialize P_SSR_best to the initial value;

(9) end for

(10) Initialize G_best and G_SSR_best based on the lowest P_SSR_best;

(11) while iteration < max generation M

(12) for each particle i = 1 to N do

(13) Update the velocity and position of particle i

(14) if the position of particle i exceeds the boundary [LB, UB]^D then

(15) the position of particle i is set to the boundary value;

(16) Substituting Data into Eq. 37 to determine all the model parameters;

(17) Compute predicted respiratory volume $\widehat{V}$ (t) based on respiratory mechanical model;

(18) Substituting $\widehat{V}$ (t) and V(t) to Eq. 38 to calculate temporary SSR (fitness) value;

(19) if P_SSR_best > temporary SSR value then

(20) P_SSR_best = temporary SSR value;

(21) P_best = location of particle i

(22) if G_SSR_best > P_SSR_best then

(23) G_SSR_best = P_SSR_best;

(24) G_best = P_best;

(25) end for

(26) end while

Furthermore, for the Grünwald-Letnikov discretization calculation formula of the open-source vfoderiv3 module in Simulink, you can refer to the following Eqs 41, 42 or the instruction file and related literature of this module (Scherer et al., 2011; Sierociuk et al., 2015; Sierociuk, 2023; Xue and Bai, 2023). When the step size h is sufficiently small and the start time involved in the calculation is t₀, the calculation of fractional calculus under the Grünwald-Letnikov definition is as follows:

${\mathcal{D}}_{t_0}^{\alpha }f\left(t\right)=\frac{1}{h^{\alpha }}{\displaystyle \sum _{j=0}^{\frac{t-t_0}{h}}}{\omega }_{j}f\left(t-jh\right)(41)$

where the coefficient ${\omega }_j$ in the above formula (41) has the following recursive formula that is convenient for numerical calculation:

${\omega }_0=1,{\omega }_j=\left(1-\frac{\alpha +1}{j}\right){\omega }_{j-1},j=1,2,\cdots (42)$

For instance, applying it to calculate the α fractional derivative of the respiratory volume V(t) at a certain moment t, is to substitute the measured values of V (t₀) to V(t) into Eq. 41 for calculation. The coefficients ${\omega }_j$ of Eq. 41 are determined by recursive calculation from Eq. 42.

3 Results

3.1 Our results indicate that the parameter estimation and model simulation methods proposed in this paper are fast and effective for all the considered classes of models

Four representative different models have been employed to analyze 150 s of PC-AC mode data (corresponding to case1PCAC1.txt, 15,000 data points) for parameter estimation and model simulation. Results for the different models are summarized in Table 1 along with their mean computation time obtained from three runs (all the calculations are retained to four decimal places).

TABLE 1

TABLE 1. Parameter estimation and simulation results of four representative models.

As it is apparent from Table 1, the fractional model has the best performance, although the computation time is slightly higher than the other models. For all the models, parameter estimation and model simulation and evaluation may be completed within seconds. Figure 2 shows that models reproduce accurately the measured respiratory volumes V(t) data.

FIGURE 2

FIGURE 2. Comparison of experimental and simulated V-T curves for different representative models. Models 1 to 4 correspond to P = EV + RV’ + P0, P = EV + RV’ + IV'’ + P0, P = EV + RV’ + P0 + a*V^2, and P = EV + RV’ + a*D1.1 (V) + P0 respectively. (A) 10 s simulation results for each model, (B) zoom of the orange region in (A).

Looking at the first row of Table 1 as an example, we see that the values of respiratory system elastance E = 26.1130 cmH2O/L and respiratory system resistance R = 7.1716 cmH2O*s/L obtained from Eq. 37 both fall within the known range of normal measurements (Protti et al., 2011; Henderson et al., 2017).

In addition, compared to the study by Albanese et al. (2016), our modelling provides an improvement in accuracy by an order of magnitude in terms of the CV(RMSE). In particular, their model achieved an optimal CV(RMSE) = RMSE/average (V(t)) = 0.079 for the volume V(t). While just using the most basic Eq. 1 we find CV(RMSE) = 0.0089 (the results of other models are even better). This demonstrates the accuracy and reliability of our modelling and parameter estimation methods.

3.2 The double-exponential model can be chosen as a representative of exponential models

As mentioned above, we have proposed the exponential Eqs 30–32 to describe data measured in three different ventilation modes: PC-AC, VC-AC, VC-SIMV (corresponding to data files case1PCAC1.txt, case1VCAC1.txt, case1VCSIMV1.txt). The x parameter in the exponentials of Eqs 30–32 is optimized using PSO to obtain the most suitable value, and the criterion for optimization is the SSR value (lower SSR value is better). The number of particles N and the number of iterations M are taken as 5, 10, 20, respectively. The SSR and computation time of each model are summarized as shown in Figure 3.

FIGURE 3

FIGURE 3. Comparison of SSR and computation time for Eqs 30–32 in three respiratory modes. The left y-axis in the Figure denotes time (in seconds), corresponding to the line graph, while the right y-axis represents the SSR values, corresponding to the bar graph. The x-axis is the particle number N (equal to the number of iterations M). The purple line and bar denote results for Eq. 30, green is for Eq. 31, and blue is for Eq. 32. (A) reports results for the PC-AC ventilation mode, (B) for the VC-AC ventilation mode and (C) for the VC-SIMV ventilation mode.

Overall, the computational results improve (lower SSR values) as the exponent increases (i.e., at a faster growth rate), suggesting that the respiratory system elasticity E has certain step characteristics (the exponential term corresponding to the component separated from the elasticity E), and the closer is to the step characteristics, the better are the results.

Looking at the dashed line in the graph, we conclude that Eq. 32 shows a significant improvement compared to Eq. 30 (with a significant difference of p-value = 0.0001009, and an average difference rate of 12.91% in the three ventilation modes). However, there are no significant differences between Eq. 32 and 31 (with a non-significant difference of p-value = 0.8271, and an average difference rate of 1.05% in the three ventilation modes). Further increasing the exponential order results only in an increase of the computational time. We thus recommend using the double exponential model as a representative of this type of models. The particle number N and iteration rounds M for this type of model PSO optimization are both set to 10, which allows us to achieve saturation in calculations.

3.3 Exponential models provide superior results compared to the polynomial models in volume-controlled ventilation mode

In order to investigate and find out a suitable nonlinear model to describe the actual respiratory data, we simulate the dynamics of Eq. 27, with the exponent ranging from 1 to 11, and select the best model based on its SSR value. When x = 1 Eq. 27 is equivalent to Eq. 29, whereas for y = 1 Eq. 27 is equivalent to Eq. 28. When both x = 1 and y = 1, Eq. 27 reduces to (1). We use three different ventilation modes, PC-AC, VC-AC and VC-SIMV, to estimate the parameters of each model using measured data (corresponding to data files case1PCAC1.txt, case1VCAC1.txt, case1VCSIMV1.txt), and only the parameter values that meet the Routh-Hurwitz criterion and the actual physical meaning are taken for subsequent simulation calculations of SSR values. The results may be summarized as follows:

Figure 4 shows that the exponent y only weakly influences the overall SSR. To reduce the undetermined parameters Eq. 28 is thus taken as the more suitable model for the discussion in the subsequent sections.

FIGURE 4

FIGURE 4. SSR of nonlinear polynomial models. (A–C) illustrate the behaviour of the SSR as a function of the exponent x of the aV(t)^x term in Eq. 27. The fitting functions corresponding to the mean of all curves in the 3 ventilation modes are shown in Figures 5A–C, respectively (retaining 2 decimal places). (D–F) illustrate the behaviour of the SSR as a function of the exponent y of the bV′(t)^y term in Eq. 27. The numbers in x1, x2. and y1, y2. in the legend represent the values corresponding to x and y in Eq. 27.

The SSR decreases with the exponent x in the term aV(t)^x with an approximately exponential behaviour. We take the mean and the logarithm of the mean of the SSR values of all curves in Figures 4A–C respectively, and through linear regression analysis, we find that the logarithm of the mean of the SSR is more significantly correlated to the independent variable x (higher coefficient of determination r2, see Supplementary Material Section 11). This suggests that the SSR of Eq. 28 converges slowly and other may be more convenient in terms of convergence and model fitting.

To simplify and more clearly illustrate our results, we further take the best results obtained by using polynomial Eqs 27–29 of orders 1–11, and compare them with results obtained from Eq. 31. Comparison is made in terms of SSR for the three different ventilation modes, PC-AC, VC-AC, VC-SIMV, for parameter estimation and simulation (corresponding to data files case1PCAC1.txt, case1VCAC1.txt and case1VCSIMV1.txt). For the double exponential Eq. 31, the exponent x is optimized by PSO for a number of particles N = 10 and iterations M = 10.

As can be seen from the results in Table 2, in the Volume-Controlled ventilation mode (VC-AC, VC-SIMV), the double exponential model Eq. 31 has significantly superior BIC and SSR, obtained in a computational time of the same order of magnitude compared with the other models. According to Eq. 26, the exponential term $a{e^e}^{xV\left(t\right)}V\left(t\right)$ in Eq. 31 and the aV(t)^x term in Eq. 28 are physically related to the elasticity E (component of respiratory elasticity E), and the results of this Section indicates that the elasticity E of the respiratory system has a step-like characteristic, which is more suitable to be modeled and described by exponential or similar fast-growing functions.

TABLE 2

TABLE 2. Comparison of double-exponential model and polynomial models.

3.4 The fractional calculus model provides superior results compared to integer-order differential models in Pressure-Controlled ventilation mode

We further compare results for the Pressure-Controlled ventilation mode, introducing the fractional calculus term into the respiratory mechanical model.

Before comparing Eq. 21 with other different models, the undetermined fractional order in the fractional calculus model of Eq. 21 should be estimated using the PSO algorithm to obtain the optimal solution. First, we use PSO algorithm to determine for which values of the particle number N and iteration rounds M one obtains a stable optimal solution for Eq. 21. We use data files from three different breathing patterns: case1PCAC1.txt, case1VCAC1.txt, case1VCSIMV1.txt, and compare SSR for N and M ranging from 3 to 30.

The results are shown in Figure 5. Using M = N = 5, PSO already reaches a stable optimal solution, and increasing N and M provides only little improvement on the final result. For N = M = 5, the maximum difference rate from the optimal SSR (N = 30 and M = 30) in the three ventilation modes is 0.14% for the PC-AC mode. Hence, for the PSO algorithm optimization of Eq. 21 we consistently set N = M = 5.

FIGURE 5

FIGURE 5. SSR values as obtained from different PSO parameters in the three ventilation modes for Eq. 21. The black line is for the PC-AC mode, green for VC-AC mode, and blue for VC-SIMV mode. The horizontal axis represents the values of M and N in the PSO algorithm, for example, the position of x = 5 denotes N = 5 and M = 5. The SSR value data points calculated from different M and N values are marked with red “x” in the Figure.

After determining that Eq. 21 has reached a stable optimal solution based on the PSO algorithm with particle number N = 5 and iteration number M = 5, we can further compare the single-compartment model Eq. 1, fractional calculus model Eq. 21, polynomial model Eqs 27–29, and double-exponential model Eq. 31 based on the actual measured data from the three ventilation modes case1PCAC1.txt, case1VCAC1.txt, case1VCSIMV1.txt. We calculate their respective BIC values, SSR values, and computation time. The PSO algorithm is used to optimize the fractional order of Eq. 21 and the exponential parameter of Eq. 31, with the goal of minimizing SSR. In the optimization calculation based on the PSO algorithm, for Eq. 21, the particle number N = 5 and the iteration number M = 5 are taken, while for Eq. 31, according to the previous results, N = M = 10 is taken to achieve saturation. To simplify and clearly present our results, the best results from the polynomial models up to order 11 are compared, and the results are shown in Table 3.

TABLE 3

TABLE 3. Comparison of fractional calculus model and integer-order differential models.

From the above results, it is apparent that the fractional calculus model performs better (lower BIC and SSR values) compared to the single-compartment Eq. 1 for all the ventilation modes. The running time is generally within half a minute on our computing platform. The additional computation time of the fractional order model is mostly due to the iterative search, required to find the suitable order, whereas the computational time required for saturation is shorter compared to the exponential Eq. 31. In particular, for the Pressure-Controlled ventilation mode (PC-AC), the fractional calculus model of Eq. 21 is superior to other types of models.

3.5 Significant improvement of respiratory model fitting based on hybrid model

As illustrated in the previous Sections, Eq. 21 is optimal for the Pressure-Controlled ventilation mode and Eq. 31 for the Volume-Controlled ventilation mode. In this Section, we investigate whether the hybrid model in Eq. 33, which has features of both models may achieve better results (lower BIC and SSR values) for all the ventilation modes.

Since the fractional order α and the exponential term parameter β in the hybrid Eq. 33 are obtained by PSO algorithm, we first calculate the number of particles N and the number of iterations M required to achieve stable optimal solution based on data from the three different ventilation modes, case1PCAC1.txt, case1VCAC1.txt and case1VCSIMV1.txt.

From Figure 6, we see that the hybrid Eq. 33 reaches stable optimal solutions for all ventilation modes when N = M = 20, and that there is no statistical difference compared with the optimal SSR (N = 50 and M = 50) (p-value = 0.1835, average difference rate 0.0087% in three ventilation modes). In the following, we use N = M = 20 to optimize the parameters of the hybrid model Eq. 33.

FIGURE 6

FIGURE 6. SSR values obtained from different PSO parameters in three ventilation modes for Eq. 33. The horizontal axis denotes the values of M = N in the PSO algorithm, and the SSR value data points calculated from different M and N values are marked with red “x” in the Figure.

Next, we compare the BIC and computational time of the single-compartment model Eq. 1, fractional calculus model Eq. 21, polynomial model Eqs 27–29, double-exponential model Eq. 31 and hybrid model Eq. 33 using the same data (case1PCAC1.txt, case1VCAC1.txt, case1VCSIMV1.txt). For the polynomial model Eqs 27–29, we take the best results from orders 1–11. For Eqs 21, 31, 33, which need to be optimized based on the PSO algorithm, we take their corresponding stable optimal solutions, i.e., N = M = 5, N = M = 10 and N = M = 20, respectively. The results are shown in Figure 7A. After parameter estimation, we further use the relationship between pressure P(t) and respiratory volume V(t) in the data of case1PCAC2.txt, case1VCAC2.txt, case1VCSIMV2.txt for cross-validation (i.e., the pressure in the second time period is used as the input, and the model parameters estimated in the first time period are used for model simulation. The simulated volume $\widehat{V}\left(t\right)$ results are compared with the actual measured respiratory volume V(t) in the second time period to calculate SSR and BIC values), and the results of their BIC, SSR values and computational time correspond to Figure 7B.

FIGURE 7

FIGURE 7. BIC and computational time for the three respiratory modes. The left axis in the figure represents BIC value (bar graph), while the right axis represents time (line graph). The horizontal axis represents three different ventilation modes [black for Eq. 1, green for (Eq. 21), cyan for (Eqs 27–29), blue for (Eq. 31), and purple for (Eq. 33)]. (A) The results based on data from the first time period (case1PCAC1.txt, case1VCAC1.txt, case1VCSIMV1.txt). (B) The cross-validation results based on data from the second time period (case1PCAC2.txt, case1VCAC2.txt, case1VCSIMV2.txt), with the parameter estimation based on data from the first time period.

Results in Figure 7A show that our proposed hybrid model Eq. 33 is the optimal model for all ventilation modes, and that its BIC value is significantly better than other models. When applied to actual measured data from different time periods of the same mode after estimating the model parameters (Figure 7B), similar results are obtained (the behaviour of the BIC for each model is similar), that is, the hybrid model is the optimal model within a given time period and for cross-validation in different time periods.

In terms of computation time, the optimal hybrid model also takes the longest time for parameter estimation stage. The single-compartment model has the worst BIC value and the shortest calculation time (average 0.26 s). Since the hybrid model involves two parameters to be estimated, the main overhead in computation is the optimization of parameters by the PSO algorithm. However, if the personalized characteristic parameters are known, the single-round running time for each model is on the same order of magnitude (Figure 7B), i.e., less than 1 s on our computing platform.

In addition, in order to verify the consistency of the results of each model in different experimental animals, we validated results on measured data from another two Bama pig for the VC-AC mode. First, based on the measured data file case2VCAC1.txt, we estimated the parameters and calculated the BIC and SSR values to evaluate the results of each model. Meanwhile, the model parameters obtained from case2VCAC1.txt were used for cross-validation of the relationship between pressure P(t) and respiratory volume V(t) in the measured data of case2VCAC2.txt. Moreover, the validation for case3 data (case3VCAC1.txt and case3VCAC2.txt) is the same. The calculated results of BIC and SSR values are summarized in Table 4.

TABLE 4

TABLE 4. Validation of results based on VC-AC mode case2 and case3 data.

The validation results are consistent with the results in Figure 7, indicating that the hybrid model Eq. 33 can achieve significantly better results for different ventilation modes. Compared with the basic single-compartment model Eq. 1, the SSR value of hybrid model Eq. 33 is nearly halved and the difference in BIC value is even more significant. In addition, the results of Table 4 also re-verify that if personalized characteristic parameters are obtained based on several respiratory cycle data, they can be used for simulations under the same ventilation mode, and even predict the respiratory volume V(t) curve based on the specified pressure P(t), without having to recalculate the parameter values, with a single round of calculation that can be completed in seconds.

In summary, if there is sufficient computing power or quick calculations is not required, the hybrid model is recommended as the optimal model to universally describe different ventilation modes. When quick calculations are required, the fastest model is the basic Eq. 1. For the same individual, the model parameters that conform to its personalized characteristics can also be pre-calculated based on measured data (Figure 7B), and respiratory volume V(t) can be directly calculated based on pressure P(t) in the same conditions.

4 Discussion

We have addressed effective modelling of respiratory mechanics and put forward a set of methods for rapid parameter estimation based on measured data (each model round calculation is less than 1 s), model simulation and evaluation system, also providing the corresponding executable code (https://gitee.com/lizwtest/mechanical_ventilation) and measured data. Using our approach, different models and parameters can be selected and tailored to different applications to further explore the mechanical characteristics of the respiratory system, improve the ventilation mode of the ventilator, and lay the foundation for the implementation of rapid adaptive parameter ventilation.

Some models analyzed in the article have time-domain general solutions, and so when P(t) in the respiratory system is very close to some elementary functions (or multiple superpositions), the volume V(t) curve (or multiple solutions superposition in linear systems) can be obtained directly through the time-domain general solution. However, considering certain autonomous randomness of individual’s breathing and the presence of background, more realistic simulations may be obtained by substituting each point in time order.

In addition, the actual elasticity E of respiratory systems shows a step characteristic. In Volume-Controlled ventilation modes, the pressure P(t) is not fixed, and decreases as it approaches the expected tidal volume, causing a sudden change in the pressure P(t) and respiratory volume V(t) curve. In Pressure-Controlled ventilation modes, the pressure is fixed, which reflects the memory-dependent viscoelasticity and heterogeneity characteristics of lung tissue, whereas in volume-controlled ventilation modes, the pressure P(t) has not a fixed value. The pressure P(t) decreases when it is close to the expected tidal volume, and if it reaches the thoracic/lung volume, the elasticity E of the respiratory system will significantly increase (a larger pressure increment ΔP is required to produce the same volume increment ΔV). This corresponds to a sudden change in the pressure P(t) and volume V(t) curve, and the volume V(t) curve is more likely to reflect the plateau feature. In Pressure-Controlled ventilation modes, the pressure is fixed to reflect the memory-dependent viscoelasticity and heterogeneity characteristics of lung tissue.

Finally, using the systematic approach proposed in this paper, we have proposed and validated an optimal hybrid model useful for various ventilation modes. To this aim, we have for the first time applied fractional calculus and double exponential terms to the modelling of respiratory mechanics. Compared to other models, the hybrid one shows significantly better BIC and SSR results for different ventilation modes, i.e., it better describes the actual pressure P(t) and respiratory volume V(t) relationship. This implies that we can more accurately describe and simulate the characteristics of respiratory mechanics. Based on this model, we can further study the mechanical characteristics of the respiratory system, deepen our understanding of the impact of mechanical ventilation on patients, and optimize the ventilation mode of the ventilator to improve ventilation effects and reduce ventilation damage.

Furthermore, after pre-calculating personalized parameters using several respiratory cycles of different individuals, we can also implement rapid adaptive parameter ventilation based on the hybrid model. Traditional ventilation mode parameters are generally fixed and cannot be dynamically adjusted according to the patient’s real-time respiratory characteristics. However, in the hybrid and other models and methods mentioned in this article, we can dynamically adapt the model and estimate the corresponding model parameters according to the actual data measured at different time intervals of individuals. This makes the overall approach more personalized and adaptable to the actual situation of the patient, thereby achieving better results.

5 Conclusion

We have proposed a system for parameter estimation, model simulation, and evaluation based on actual measured data, and after carefully comparing different models, we have put forward an optimal hybrid model valid for various ventilation modes. It may not only help doctors to better understand the mechanical characteristics of the respiratory system, but also improve the ventilation mode of the ventilator, enhance the effect of respiratory therapy, and provide strong support for the realization of rapid adaptive parameter ventilation, which is a relevant feature to improve the safety, effectiveness, and personalization level of mechanical ventilation.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Ethics statement

The studies involving animal participants were reviewed and approved by the Chinese PLA General Hospital ethics committee (reference: 2021-NO.-772). The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

ZL: Conceptualization, Data curation, Formal Analysis, Methodology, Software, Visualization, Writing–original draft. YP: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Validation, Writing–review and editing. YW: Resources, Supervision, Writing–review and editing. QT: Funding acquisition, Project administration, Resources, Supervision, Writing–review and editing.

Funding

The authors declare financial support was received for the research, authorship, and/or publication of this article. This work was supported by the Key Technologies Research and Development Program of China under grant no. 2020YFC0122202.

Acknowledgments

The authors would like to express their gratitude to EditSprings (https://www.editsprings.cn) for the expert linguistic services provided.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2023.1273645/full#supplementary-material

References

Adams M. (2019). differint: a python package for numerical fractional calculus. arXiv preprint arXiv:1912.05303.

Google Scholar

Albanese A., Cheng L., Ursino M., Chbat N. W. (2016). An integrated mathematical model of the human cardiopulmonary system: model development. Am. J. Physiology-Heart Circulatory Physiology 310 (7), H899–H921. doi:10.1152/ajpheart.00230.2014

PubMed Abstract | CrossRef Full Text | Google Scholar

Avilés-Rojas N., Hurtado D. E. (2022). Whole-lung finite-element models for mechanical ventilation and respiratory research applications. Front. Physiology 13, 2060. doi:10.3389/fphys.2022.984286

CrossRef Full Text | Google Scholar

Bates J. H. (2009). Lung mechanics: an inverse modeling approach. Cambridge University Press.

Google Scholar

Carvalho A. R., Zin W. A. (2011). Respiratory system dynamical mechanical properties: modeling in time and frequency domain. Biophys. Rev. 3, 71–84. doi:10.1007/s12551-011-0048-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen W., Sun H., Li X. (2010). Fractional order derivative modeling of mechanics and engineering problems. Beijing: Science Press.

Google Scholar

Chiew Y. S., Pretty C., Docherty P. D., Lambermont B., Shaw G. M., Desaive T., et al. (2015). Time-varying respiratory system elastance: a physiological model for patients who are spontaneously breathing. PLoS ONE 10 (1), e0114847. doi:10.1371/journal.pone.0114847

PubMed Abstract | CrossRef Full Text | Google Scholar

Dincel E. (2021). Advanced mechanical ventilation modes: design and computer simulations. Comput. Methods Biomech. Biomed. Engin. 24 (6), 673–686. doi:10.1080/10255842.2020.1845319

PubMed Abstract | CrossRef Full Text | Google Scholar

Divo M. J., Martinez C. H., Mannino D. M. (2014). Ageing and the epidemiology of multimorbidity. Eur. Respir. J. 44 (4), 1055–1068. doi:10.1183/09031936.00059814

PubMed Abstract | CrossRef Full Text | Google Scholar

Eager D., Pendrill A.-M., Reistad N. (2016). Beyond velocity and acceleration: jerk, snap and higher derivatives. Eur. J. Phys. 37 (6), 065008. doi:10.1088/0143-0807/37/6/065008

CrossRef Full Text | Google Scholar

Ellwein Fix L., Khoury J., Moores R. R., Linkous L., Brandes M., Rozycki H. J. (2018). Theoretical open-loop model of respiratory mechanics in the extremely preterm infant. PLoS ONE 13 (6), e0198425. doi:10.1371/journal.pone.0198425

PubMed Abstract | CrossRef Full Text | Google Scholar

Gad A. G. (2022). Particle swarm optimization algorithm and its applications: a systematic review. Archives Comput. Methods Eng. 29 (5), 2531–2561. doi:10.1007/s11831-021-09694-4

CrossRef Full Text | Google Scholar

Ganzert S., Möller K., Steinmann D., Schumann S., Guttmann J. (2009). Pressure-dependent stress relaxation in acute respiratory distress syndrome and healthy lungs: an investigation based on a viscoelastic model. Crit. Care 13, 1999–R210. doi:10.1186/cc8203

PubMed Abstract | CrossRef Full Text | Google Scholar

Gertler R. (2021). Respiratory mechanics. Anesthesiol. Clin. 39 (3), 415–440. doi:10.1016/j.anclin.2021.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghafarian P., Jamaati H., Hashemian S. M. (2016). A review on human respiratory modeling. Tanaffos 15 (2), 61–69.

PubMed Abstract | Google Scholar

Henderson W. R., Chen L., Amato M. B., Brochard L. J. (2017). Fifty years of research in ARDS. Respiratory mechanics in acute respiratory distress syndrome. Am. J. Respir. Crit. Care. Med. 196 (7), 822–833. doi:10.1164/rccm.201612-2495CI

PubMed Abstract | CrossRef Full Text | Google Scholar

Hess D. R. (2014). Respiratory mechanics in mechanically ventilated patients. Respir. Care 59 (11), 1773–1794. doi:10.4187/respcare.03410

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu S. (2007). Principles of automatic control. Beijing: Science Press.

Google Scholar

Ionescu C., Kelly J. F. (2017). Fractional calculus for respiratory mechanics: power law impedance, viscoelasticity, and tissue heterogeneity. Chaos, Solit. Fractals 102, 433–440. doi:10.1016/j.chaos.2017.03.054

CrossRef Full Text | Google Scholar

Kennedy J., Eberhart R. (1995). Particle swarm optimization. Proc. ICNN’95 - Int. Conf. Neural Netw. 4, 1942–1948. doi:10.1109/ICNN.1995.488968

CrossRef Full Text | Google Scholar

Li N. (2011). Analog circuits. Beijing: Tsinghua University Press.

Google Scholar

Marconi S., De Lazzari C. (2020). In silico study of airway/lung mechanics in normal human breathing. Math. Comput. Simul. 177, 603–624. doi:10.1016/j.matcom.2020.05.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Protti A., Cressoni M., Santini A., Langer T., Mietto C., Febres D., et al. (2011). Lung stress and strain during mechanical ventilation: any safe threshold? Am. J. Respir. Crit. Care. Med. 183 (10), 1354–1362. doi:10.1164/rccm.201010-1757OC

PubMed Abstract | CrossRef Full Text | Google Scholar

Redmond D. P., Kim K. T., Morton S. E., Howe S. L., Chiew Y. S., Chase J. G. (2017). A variable resistance respiratory mechanics model * *the authors acknowledge the support of the health research council (HRC) of New Zealand. IFAC-PapersOnLine 50 (1), 6660–6665. doi:10.1016/j.ifacol.2017.08.1533

CrossRef Full Text | Google Scholar

Scherer R., Kalla S. L., Tang Y., Huang J. (2011). The Grünwald–Letnikov method for fractional differential equations. Comput. Math. Appl. 62 (3), 902–917. doi:10.1016/j.camwa.2011.03.054

CrossRef Full Text | Google Scholar

Shami T. M., El-Saleh A. A., Alswaitti M., Al-Tashi Q., Summakieh M. A., Mirjalili S. (2022). Particle swarm optimization: a comprehensive survey. IEEE Access 10, 10031–10061. doi:10.1109/ACCESS.2022.3142859

CrossRef Full Text | Google Scholar

Sheng Z., Xie S., Pan C. (2020). Probability theory and mathematical statistics. Beijing: Higher Education Press.

Google Scholar

Sierociuk D. (2023). Fractional variable order derivative Simulink toolkit, Available online: https://www.mathworks.com/matlabcentral/fileexchange/38801-fractional-variable-order-derivative-simulink-toolkit.

Google Scholar

Sierociuk D., Malesza W., Macias M. (2015). Derivation, interpretation, and analog modelling of fractional variable order derivative definition. Appl. Math. Model. 39 (13), 3876–3888. doi:10.1016/j.apm.2014.12.009

CrossRef Full Text | Google Scholar

Soriano J. B., Kendrick P. J., Paulson K. R., Gupta V., Abrams E. M., Adedoyin R. A., et al. (2020). Prevalence and attributable health burden of chronic respiratory diseases, 1990–2017: a systematic analysis for the Global Burden of Disease Study 2017. Lancet Respir. Med. 8 (6), 585–596. doi:10.1016/S2213-2600(20)30105-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Drunen E. J., Chiew Y. S., Pretty C., Shaw G. M., Lambermont B., Janssen N., et al. (2014). Visualisation of time-varying respiratory system elastance in experimental ARDS animal models. BMC Pulm. Med. 14 (1), 33–39. doi:10.1186/1471-2466-14-33

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang G., Zhou Z., Zhu S., Wang S. (2007). Ordinary differential equations. Beijing: Higher Education Press.

Google Scholar

Wunsch H. (2020). Mechanical ventilation in COVID-19: interpreting the current epidemiology. Am. J. Respir. Crit. Care. Med. 202 (1), 1–4. doi:10.1164/rccm.202004-1385ED

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie M., Liu X., Cao X., Guo M., Li X. (2020). Trends in prevalence and incidence of chronic respiratory diseases from 1990 to 2017. Respir. Res. 21 (1), 49–13. doi:10.1186/s12931-020-1291-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie Z., Huang X., Liu W. (2022). Subpopulation particle swarm optimization with a hybrid mutation strategy. Comput. Intell. Neurosci. 2022, 9599417. doi:10.1155/2022/9599417

PubMed Abstract | CrossRef Full Text | Google Scholar

Xue D., Bai L. (2023). Fractional calculus: numerical algorithms and implementation. Beijing: Tsinghua University Press.

Google Scholar

Yang M., Liu Y., Yang J. (2022). A hybrid multi-objective particle swarm optimization with central control strategy. Comput. Intell. Neurosci. 2022, 1522096. doi:10.1155/2022/1522096

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong Y. (1988). Complex analysis. Beijing: Higher Education Press.

Google Scholar