Method for modelling respiratory system based on fractional calculus

Abstract

A method for modelling a respiratory system based on fractional calculus, comprising a fractional mechanical model, including an optimal fractional order and various coefficient parameters, based on a

Description

[0001] The present disclosure relates to a respiratory system model, in particular to a method for modeling a respiratory system based on fractional calculus. BACKGROUND

[0002] In the study of mechanical ventilation of a respiratory system, lung tissue has a certain power-law characteristics, memory properties and path-dependent properties, the lung tissue is not an ideal viscoelastic material. Therefore, it is difficult for the classical lumped parameter model based on integer differential equations to describe and reflect these characteristics well. Classical calculus is to obtain integer differential and integral of functions. If non-integer differential and integral are calculated, the field of fractional calculus is involved. At present, there is no method to model the mechanics of the respiratory system by using the fractional calculus based on the measured data in any respiratory mode in the prior art. Therefore, there is a lack of a modeling method based on fractional calculus in the field of modeling the mechanics of the respiratory system at present. SUMMARY

[0003] The present disclosure provides a method for modeling a respiratory system based on fractional calculus, which shows a better fitting degree than a classical model in the measured data.

[0004] In order to solve the above technical problems, the technical scheme of the present disclosure is as follows.

[0005] The present disclosure relates to a method for modeling a respiratory system based on fractional calculus, where:

[0006] expressing a fractional mechanical model as:

[0007] p = EV+ RV' + W(K(0) + W(V(t)) + - + kmD^(K(t)) + Po

[0008] where represents calculating cq -th fractional derivative of V(t), D“2(U(0) represents calculating a2-th fractional derivative of V(t), D“m(V(t)) represents calculating am -th fractional derivative of V(t), m is a number of fractional terms of the model, V is a respiratory volume, V is a first-order derivative of V, which is a gas flow rate, E is an elasticity of the respiratory system, R is an airway resistance, Po is a trans-pulmonary pressure, a2, -, am are fractional orders, and klt k2,..., km are coefficient parameters of the fractional mechanical model; 2 = [ yt y\ l yt p

[0009] obtaining fitting results of various coefficients by using [ | / anj iteratively searching an optimal order, where

[0010] km p . 0 V, V’ ... v; 1 v2 K v} ... v™ 1

[0011] , X is a data matrix of n rows and m+3 columns,

[0012] P is an n-dimensional air pressure vector, which consists of values of a respiratory circuit acquired at n points in time.

[0013] In some embodiments, in the method for modeling the respiratory system based on fractional calculus, a fractional mechanical model of a single fractional derivative term is expressed as: P = EV RV+P()

[0014] In some embodiments, in the method for modeling the respiratory system based on fractional calculus, a fractional mechanical model of two fractional derivative terms is expressed as:

[0015] P = EV + RV' + k^ (^(0) + k:\);’'(V(t)) P.

[0016] In some embodiments, in the method for modeling the respiratory system based on fractional calculus, in the process of iteratively searching the optimal order, parameters to be optimized are a series of orders a2, a^, all coefficient parameters of model are obtained upon specifying a set of orders and in turn a sum of squared residuals or a root mean square error is obtained, and the respective coefficient parameters are continuously updated following a change trend of the sum of squared residuals or the root mean square error so as to minimize the sum of squared residuals or the root mean square error.

[0017] The present disclosure has the following beneficial effects.

[0018] In comparison with a general integer calculus model, the method for modeling the respiratory system based on fractional calculus according to the present disclosure can reflect power-law characteristics and memory properties of the respiratory system. It is easier for the method to obtain a better fitting degree than the classical integer calculus model, and no specific respiratory model is required for data acquisition. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] In order to explain the embodiments of the present disclosure or the technical schemes in the prior art more clearly, the drawings that need to be used in the detailed description or the prior art will be briefly introduced hereinafter.

[0020] FIG. 1 is a simulation result of a fractional mechanical model of a method for modeling a respiratory system based on fractional calculus according to the present disclosure. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0021] In order to make the above objects, technical scheme and advantages of the present disclosure more obvious, the technical scheme in the embodiment of the present disclosure will be described clearly and completely with reference to the drawings in the embodiment of the present disclosure hereinafter.

[0022] The method for modeling the respiratory system based on fractional calculus of the present disclosure can quickly obtain a corresponding fractional mechanical model, including the optimal fractional order and various coefficient parameters, based on a pressure value, a flow rate value and an accumulated gas volume acquired in real time from a respiratory circuit in any ventilation mode.

[0023] The method for modeling the respiratory system based on fractional calculus according to the present disclosure includes steps SI to S3.

[0024] In step SI, according to the fractional mechanical model, the measured data are collated into the corresponding term matrix-form data, as input.

[0025] The fractional mechanical model has the general form:

[0026] p = EV+ RV' + (E(0) + k.D? (^t)) + - + kmD^ + P()

[0027] where • ' ' ” represents calculating ax-th fractional derivative of V(t), such as 0.23, 2.8, etc.; D^(K(0) represents calculating a2-th fractional derivative of V(t); D“"(V(t)) represents calculating am -th fractional derivative of V(t), m is a number of fractional-order terms of the model; Visa respiratory volume (in unit of L); V is a first-order derivative of V, that is, a gas flow rate; E is an elasticity of the respiratory system, R is an airway resistance, Po is a trans-pulmonary pressure, a1,a2,---,am are fractional orders, and kltk2,...,km are coefficient parameters of the fractional mechanical model.

[0028] In practical application, the model can also be simplified with a few fractional-order terms to reflect the power-law characteristics and memory properties, such as:

[0029] P = EV + RV' + kxDp(V(t)) + Po

[0030] or

[0031] P = EV + RV' + (7( / )) + kJ) J (V(t)) + Po

[0032] For a model with m fractional derivative terms (plus 3 non-fractional derivative terms, amounting to m+3 terms in total), the present disclosure acquires the pressure value, the flow rate value and the accumulated gas volume value of the respiratory circuit in real time. The fractional differential values are obtained according to the number m of fractional-order terms of the preset model and the corresponding form, and a specified set of fractional orders {a-y, a2, am}, the fractional differential values are collated into corresponding matrix-form data as the input of the algorithm. If values of the respiratory circuit are acquired at n points in time, the input data includes an n-dimensional air pressure vector P and a data matrix X with n rows and m+3 columns:

[0033] P = < 1 jrm

[0034] For example, column Vm of the above formula (the second column from the right) denotes values of am -th fractional derivative of the gas volume V(t), and the same is true for other corresponding fractional derivative terms.

[0035] Note: if there is only one of the flow rate value and the accumulated gas volume value due to a limitation of sampling conditions, the missing other value can be obtained by appropriate numerical differentiation and numerical integration methods.

[0036] In step S2, model coefficients of a specified order are calculated.

[0037] For the fractional model with m fractional derivatives, after a set of initial fractional orders { a2,--> am } are specified, the matrix calculation is carried out according to the following formula, so as to obtain the fitting results of various coefficients.

[0038] J = [XrX] 'XrP

[0039] The vector A including various coefficients is as follows:

[0040] A

[0041] The calculation method proposed in the above formula is derived from an extreme point where the partial derivative for each coefficient is 0 when the sum of squared residuals (SSR) is taken as a minimum value. By this step, the parameters of the specified fractional model have been obtained, and the fractional order can be further optimized.

[0042] In step S3, the optimal order is iteratively searched based on the above method.

[0043] After the coefficients of the specified fractional model are obtained, the coefficients and the measured data can be substituted into the model for simulation to calculate the sum of squared residuals (SSR) and the root mean square error (RMSE) of the corresponding model. In order to optimize an optimal series of orders based on data, the present disclosure preferably selects the Particle Swarm Optimization (PSO) among various machine learning algorithms in consideration of computational efficiency. In comparison with genetic algorithms, simulated annealing and other similar algorithms, although the particle swarm optimization may fall into the local optimal solution, the particle swarm optimization converges faster and pays more attention to timeliness in practical application.

[0044] In the iterative optimization process, parameters to be optimized are a series of orders {a1, a2, Each time a set of orders is specified, all the coefficient parameters of the model can be obtained in step S2, and in turn the SSR / RMSE can be obtained. The parameters can be continuously updated following the change trend of the SSR or the RMSE. The model optimization goal is to minimize the SSR or the RMSE. A model better than the initial state can be obtained after iterative calculation.

[0045] Note: after calculating the model coefficients of the specified order in step S2, in step S3, other similar machine learning algorithms can be selected to further optimize the fractional order. The selection of the machine learning algorithm is not the core of the method. After the same type of methods is replaced or changed, its overall framework still belongs to the method.

[0046] The calculation amount of the method is higher than that of the classical single-compartment model in the respiratory mechanics modeling. The more calculation amount of tens of seconds can obtain a better fitting effect in the measured data.

[0047] Hereinafter, with two typical models, that is, the classical single-compartment model P = EV+RV'+Po and the classical second-order model P = EV+RV'+IV"+P0 in the respiratory mechanics modeling as reference, two examples are supplemented to show the calculation results, and the parameters of the two models of the following Formula (1) and Formula (2) are estimated by using the method in the present disclosure:

[0048] P = EV + RV' + kxD?(V(t}) + P0 (1)

[0049] P = EV + RV' + k^^(y{t^ + k^^{V{t^ + P0 (2)

[0050] The parameter estimation is based on the measured data of animal experiments (measured by Drager Savina 300 ventilator). The data of the pressure P (in unit of mbar), the flow rate V (in unit of L / s) and the respiratory volume V (in unit of L) every 10ms are recorded as one row in chronological order.

[0051] Based on an example python+matlab code main.m (other *.py and *.m are python and matlab subfunctions called by the main code, and *.slx is a model file), the calculation results of parameters of three model after reading the measured data are shown in Table 1.

[0052] Table 1: Estimation Results of Parameters of Each Model

[0053] Classical model Step S2: estimating parameters Single-compartment model E=26.1130, R=7.1716, P0=3.3972, RMSE=10.699, runtime=0.3087(s) Second-order model E=26.1135, R=7.1716, 1=0.0085, P0=3.3972, RMSE= 10.799, run_time=0.6519(s) Fractional model 10 particles and 10 iterations 30 particles and 30 iterations Single-term fractional order E=26.3609, R=6.7338, P0=3.3755, RMSE=8.5794, kl=0.3638, alphal=1.1486, run_time=126.6592(s) E=26.3610, R=6.7345, P0=3.3755, RMSE=8.5842, kl=0.3632, alphal=1.1490, runUime=1342.9390(s) Two-term fractional order E=26.2672, R=7.0210, P0=3.3822, RMSE=9.5220, kl=0.0722, alphal=1.4069, k2=0.0470, alpha2=1.0559, runJ:ime=134.5450(s) E=26.1952, R=6.7190, P0=3.3898, RMSE=8.7406, kl=0.0023, alphal=1.5625, k2=0.4209, alpha2=1.0192, run_time=1029.3 016(s)

[0054] Four decimal places are reserved, the RMSE is converted into the unit of milliliter, and run time stands for running time (running on an ordinary PC with CPU i7-1165G7 @ 2.80GHz and a memory of 16GB).

[0055] It can be seen from the result that the method of the present disclosure is practical and effective, the fitting effect of the fractional model is better, and the RMSE value is lower. The fitting results of the measured data (a solid line) and the flow V (a dotted line) of Formula (1) are shown in FIG. 1, which are the results with 10 particles over 10 iterations.

[0056] In addition, taking into account the calculation time consumption and the fitting degree comprehensively, the Formula (1) can be generally considered by default based on PSO with 10 particles over 10 iterations, so that a stable optimization result can be basically obtained and the calculation is faster and the result is better. With the increase of the number of fractional differential terms and the increase of the number of the orders to be estimated, machine learning algorithms need more calculation amount to obtain stable results. For the PSO method, Formula (2) needs more particles and iterations or narrows the search range based on pre-calculated features to obtain better results.

[0057] In addition, one of various implementations of numerical calculation of fractional calculus is supplemented as a reference. There are different definitions of fractional calculus: Grunwald-Letnikov definition, Riemann-Liouville definition, Caputo definition, Erdelyi-Kober definition, etc. Each of the definitions has different application scenarios and numerical calculation methods. Taking the G-L definition as a possible calculation example, the a-th derivative of the given function f(t) is expressed as:

[0058] D / / ( / ) = limo—¥~h~  + 0............................._ Z / ) (3) ha. nrx                        J ’ j

[0059] In the above formula, h is a sufficiently small time step, r represents Gamma function, that is, T(z) = . This definition assumes that the value of the function f(t) is 0 when t<to, involving all function values from time to. It can be considered that the fractional derivative has memory properties and is suitable for the differential of a>0 and the integral of a<0.

[0060] For the case of numerical calculation, when the step h is small enough, there is the following formula:

[0061] 0- / (,) = / f ¢ / (,- ,0) (4). 0 h j=o

[0062] Moreover, the coefficient (Oj in the above Formula (4) has the following recursive formula which is convenient for numerical calculation:

[0063] = 1,^ = [1- — J = (5)-

[0064] The Gamma function has the property F(z+1) = zF(z) according to the integration by parts method in calculus, then: (-iyrGz + 1)

[0065] = r(j + - j +!) =--r(j)r(q< — j + 2) 0^ (-l)Hr(g + l) r(j + l)T(a -j + 1) r(j)T(a^j + 2) H + l)T(a-j +1) = _ a-j + 1 = j _ 0 + 1 jr(j)r(o-j + i) j j

[0066] In comparison with a general integer calculus model, the method for modeling the respiratory system based on fractional calculus according to the present disclosure can reflect power-law characteristics and memory properties of the respiratory system. It is easier for the method to obtain a better fitting degree than the classical integer calculus model, and no specific respiratory model is required for data acquisition.

[0067] The above embodiments are only specific embodiments of the present disclosure, which are used to illustrate the technical scheme of the present disclosure, rather than limiting the technical scheme. The scope of protection of the present disclosure is not limited thereto. Although the present disclosure has been described in detail with reference to the above embodiments, it should be understood by those skilled in the art that any person familiar with the technical field can still modify or improve the technical scheme described in the above embodiments or make equivalent substitutions of some technical features within the technical scope disclosed by the present disclosure. However, these modifications, changes or substitutions do not make the essence of the corresponding technical scheme deviate from the spirit and scope of the technical scheme of the embodiments of the present disclosure, and should be included in the scope of protection of the present disclosure. Therefore, the scope of protection of the present disclosure should be based on the scope of protection of the claims.

Claims

1. A method for modeling a respiratory system based on fractional calculus, wherein: expressing a fractional mechanical model as:P = EV + RV + kJX (F( / )) + k2D^ (F(0) + • • • + km(V(i» + Powherein D^(F(t)) represents calculating -th fractional derivative of V(t), D“2(F( / )) represents calculating a2 "th fractional derivative of V(t), D^(V(t)) represents calculating am -th fractional derivative of V(t), m is a number of fractional terms of the model, Visa respiratory volume, V is a first-order derivative of V, which is a gas flow rate, E is an elasticity of the respiratory system, R is an airway resistance, PQ is a trans-pulmonary pressure, a1( are fractional orders, and k1, k2,..., km are coefficient parameters of thefractional mechanical model;obtaining fitting results of various coefficients by using A = [XT X]x XT P , and iteratively searching an optimal order, whereinV[ V / ... Vf lhv; V / ... V?1 1 Iy V1 ... vm 1Jn n ■■■ n , X is a data matrix of n rows and m+3P is an n-dimensional air pressure vector, which consists of values of a respiratory circuit acquired at n points in time.

2. The method according to claim 1, wherein a fractional mechanical model of a single fractional derivative term is expressed as:P = EV + RV + kxD? (F(0) + Po.

3. The method according to claim 1, wherein a fractional mechanical model of two fractional derivative terms is expressed as:P = EV + RV' + tai(7( / )) + O (V(t)) + Po.

4. The method according to claim 1, wherein in a process of iteratively searching the optimal order, parameters to be optimized are a series of orders a2, •••, am}, all coefficient parameters of the model are obtained upon specifying a set of orders and in turn a sum of squared residuals or a root mean square error is obtained, and the respective coefficient parameters are continuously updated following a change trend of the sum of squared residuals or the root mean square error so as to minimize the sum of squared residuals or the root mean square error.

Images