Respiratory Parameters. Minimization of the Energy Cost of Breathing I. KARPOUZAS MEDIMA T - Insritut Biomkdiwl des Cordeliers. Universitk Pans VI. 15. rue de I’Ecole de Mgdecine. 75770 Pork Cedex 06, France Received I March 1984; rewed
3 June 1985
ABSTRACT We study the mechanics of the ventilatory system by mathematical modeling. WC include in our model the regulatory function of the system, which is the minimization of the energy cost required for the gaseous exchange of oxygen analysis techniques have been used for the mathematical technique,
1.
and carbon dioxide. Numerical model. and a new optimization
called ALIENOR, has been tried for the minimization
of the energy
required
INTRODUCTION
For a long time, the idea of optimization has been applied to the study of biological systems, particularly the respiratory system [l, 21. In considering the thoracopulmonary system, and particularly its gaseous exchange system, the mechanical work of the ventilatory system has been plotted as a function of respiratory frequency. It has been established that, for a given ventilation, there is a frequency which is optimal in the sense that the work is minimal [3, 41. We do not know in detail the regulatory mechanism and how it influences constriction and dilatation at the level of the bronchial structures. The efficiency of the muscular work depends on the speed of elongation and contraction, but for a given ventilatory level, the optimal working conditions are not necessarily realized. Therefore, we propose a new approach to the minimization principle applied to the ventilatory system when the mechanical work and the efficiency p of the system are taken into account simultaneously: J = (1/p)
x (mechanical 5
work)
PO,
where the numerator represents the total energy used by the system. The model consists of a system of nonlinear algebraic equations. MATilEMATICAL
BIOSCIENCES
7R:l-20
1 Elsevier Science Publishing Co., Inc., 1986 52 Vanderbilt Ave.. New York, NY 10017
1
(1986)
0025-5564/86/$03.50
2
I. KARPOUZAS I
0
/ 1 1
B
0
4
_
1 _I__________~_______ v
CT
h t
--
FIG. 1. Diagram active expiration.
of the ventilatory
system:
0
inspiratory
model,
@
model
with
First, we consider the ventilatory rates at which expiration is completely passive. This problem is solved using a method of local variations, the Vignes method [15]. Then the model is refined to take account of active expiration. We have solved this problem with the help of a new global optimization method called ALIENOR [16]. 2.
MECHANICAL
MODEL
CHOSEN
The ventilatory system carries gas between the ambient atmosphere and the organism. We have represented it by a pump in which a motor-driven piston moves in a cylinder. The ventilatory muscles are represented by a single muscle (Figure 1). Active inspiration is insured by moving the piston from y:, (the relaxation volume of the system) to V,, (th e maximum volume of the system). The eventual variations of volume from V, to VR (the residual volume) require the
RESPIRATORY
3
PARAMETERS
intervention of active expiration. The ventilatory system works as an alternating pump with a frequency v and an intensity 6’. Let VT be the tidal volume necessary for each respiratory cycle; VA, the alveolar volume, which is that part of the ventilation volume which is effective; and V,, the volume of the dead space which does not participate in the gaseous exchanges. It follows that
r;=vrv, VT = VD + VA.
3.
TREATMENT
OF EXPERIMENTAL
DATA
The distensibility of the system is a function of the volume V [5]. The pressure P,, which must be applied to maintain the volume V is taken as a cubic polynomial: PC= IV3 + JV2 + KV + L. For m different points (v], PC,), i = 1,. . , m, we form an m x 4 system with four unknowns (I, J, K, L): IV; + JV; + KV, + L = Pc,l, IV; + JV,2 + KV, + L = P<,2, (3.1) IV,’ + JV; + KV,, + L = PC,,, Taking
A=(%,)
with
a, =F4-J,
i=l,...,
B= (b,)
with
6, = PC,,
j =1,2,3,4,
m,
a system of the type AX= B is obtained, where XT = (I, J, K, L). Let R be the resistance of the bronchial and pulmonary passages to the flow [5]. Similarly, we can write the equation
R=-
aV+ b v+c
4
I. KARPOUZAS
For m different
points (U;, R,), i = 1,.
, m, we form an m X 3 system
aV,+h-cR,=V,R,, aV,+h-cR2=V,R2,
(3.2) aV,, + h - CR,, = Vn,R,,
of the type AX = B with A = (a,, ) and B = (h, ), where a,, = y
ar2 =
3
1,
or3 = -
c, = V;R,
R,,
with three unknowns a, b, and c. Considering AX= B to be a linear I x n algebraic system. At least one solution exists if I >, n. There is a unique solution of row(A) = n. Equations (3.1) and (3.2) have been solved using a least-squares method [6]. Neither a numerical solution nor an analytical method was used because of experimental errors. We want to minimize
6(X)= ,=l
We look for X* = min{ 6( X)]X E R” }. Since X* is a minimum, = 0, so that’
j&x*)
=o,
r=l
grad 6(X*)
,...,n,
r
that is,
y,$(X*)=25,[(AX*),-B,]
r
“‘;x”‘)‘=O,
r=l,...,n,
I
,=l
where ~(X*)=~&Z,[(A~)-B,]A,~=O. r ,=l
‘The residue generated
6 is generally
by the columns
of A.
not zero. It is zero if and only if R belongs
to the subspace
RESPIRATORY
PARAMETERS
Thus
2 a,(AX*),A,,= ,=l5 a,A,,B,,
t=l k
F (a,A,,A,,)
/=1
X,* = 2 aiA,,B,
r=l
for
r=l,...,n
r=l
Let
Z,. =
4, = i aA4,~ I=1
2 a,A,,B,, 1=1
where a, = 1, i = 1,. . , m, is a weighting coefficient. If the square matrix D,, is regular, then there is a solution X,* for the system. It has been solved using the Gauss method with maximum column pivots [7]. 4.
FIRST
MODEL
With the help of a mathematical modeling technique, and taking the experimental data into account, we have proposed a system of nonlinear algebraic equations as a model [8, 91: PC=
where Vi = V-t- VT/m,
Iv; + Jr-: + KV, + L,
with m a constant
(1)
to be chosen, and
aV,+h R=-
v,+c
The third equation is a purely numerical relation in which the variations Vrl with the ventilation are taken into account: v,=0.1v,+10-4. The fourth equation
of
(3)
represents the gaseous airflow: (4)
The fifth equation relates the oxygen flow consumed by the organism ( l&,) to the ventilator-y flow ( f) and the alveolar flow ( PA’,>:
(5)
I. KARPOUZAS
6
where E is a constant
which represents
the used fraction of oxygen, and
ri= V,V.
(6)
We use the Nelson relation [19], which represents the known variations of the muscular strength F of any skeletal muscle as a function of its speed of elongation and contraction. With parameters for the ventilatory system and for active inspiration (k’ > VO),we have
Let t be the duration
of the active muscular
intervention.
Then
t = 0.3T +0.4. We minimize the following regulation and for different energy levels
equation
with respect to V and k,
wherepmax is constant. 4.1.
DEFINITION
The following (P) Find
OF THE
PROBLEM
problem [lo-121 must be solved mathematically:
S* E A such that min
J( S*) =
J(S)
SEACR~
with
or the problem
with penalty function:
(Q) Find (S*, e) E R2 J(S*,
X
E) = min SGRZ
Ri
such that
;
[(min{ A’,-X,,O})’
I=1
+(min{
X, - n, ,O})‘]
II
RESPIRATORY
7
PARAMETERS
with Xi =V, X,=c, 10-3, n,=0.1x10-3. Let ( r’oz, .) be a subdivision
N,=6~10-~,
N2 =2x10m3,
n, = 2.4~
of the interval being studied: r’,,., = voio,,O+ ?+r
r=O
at M,
/2=3x10-6,
so that
Ijo2.0 = 4x10-6,
l&M
=4x10-j
At step r, we know Sr,in. We have to solve (Q,)
Find (ST, E) E R* X R: such that J,(S*,e)
= min J(S,,e). S,E R=
We can find a solution
of Q, for each given ventilatory level po2,T (r = 0 at M). In [13] and [14], theorems about the existence of solutions of the problems posed in (P), (Q), and (Q,) and for the convergence of the method are given. Problem (Q,) can be treated by the method of local variations (Vignes method) [15]. 4.2.
RESULTS
The theoretical curves obtained (Figures 2, 3,4, 5) are not at variance with the physiological variations observed when the energy consumption Qo, increases at a moderate rate, that is, when there is no intense muscular exercise. In this first, very simplified approach, many facts have been neglected. Thus, we cannot expect to obtain theoretical results which exactly reproduce the physiological data. In particular, the variations of the ventilatory frequency as a function of the energy level parallel those of VT more closely than seems to be appropriate. 5.
SECOND
MODEL
In the first model we supposed pm, constant. Thus, its value does not occur in the optimization of J. Here, pmax is a function of V, and three
1. KARPOUZAS
m
d ii
i
RESPIRATORY
PARAMETERS
10
I. KARPOUZAS
RESPIRATORY PARAMETERS
11
12
I. KARPOUZAS TABLE
li,,
V (m’)
(m’/sec)
1
P (m3/sec)
J, minimum
2.29954232
x 10 - ’
1.85877428
x 10
’
1.50175642x10-
’
7x10
h
2.29541636
x 10 - ’
2.92480835
x 10
4
2.03361731
x 10
’
1x10
5
2.28409946
x 10.
’
3.88385386
x 10
4
2.51139471
x 10
’
2.27876156
x 10 - ’
4.79770733
x 10
4
2.9680648
x 10 ~’
2.27398923
x 10
’
5.66196140
X 10
’
3.40770716
x 10
2.26790186
x 10
’
6.5570575X
x 10
4
3.83509333
x 10
3
7.410485
x10
4
4.25376535
x 10
’
4x10mh
1.3x1o-5 1.6x10
’
1.Yx1o-i
’
2.2x10
5
2.26246217
x 10
’
2.5x10
5
2.25843082
x 10
’
8.20698569
x 10
4
4.66540031
x 10
’
2.x x 10
5
2.25252664
x 10
’
9.05278183
x 10
4
5.07245055
x 10
’
3.1x10
s
2.24781923
x 10
’
9.85419963
x 10
’
5.4753144
x 10
3
2.25400559
x 10
’
1.07493064X
10
’
5.88310033
x 10
x
2.24193589
x 10
’
1.13473113
x 10
3
6.27271225
x 10
’
3.4x10-5 3.7x10
different
5
hypotheses
have been proposed:
(a) pm=.1 = p(V) reaches its maximum for the I$ relaxation (Figure 6). (b) pmax.z = p(V) as above, but with less decrease during inspiration (V > V,) (Figure 7). (c) Pmax.3 = P ( V) is like P,,,~~.~, but its maximum is not reached for V = VO of relaxation (Figure 8). Polynomial interpolation within segments was used for this second model, because we wanted to interpolate only a few points on a certain graph. To elaborate this model (with active expiration), we have used all the equations of the first model (some of them modified) and some new equations: (a) The mean value of IP, 1:
with P, (u) = Iu3 + Ju2 + Ku + L (already calculated (b) The mean value of [RI:
in the first model).
RESPIRATORY
PARAMETERS
13
I. KARPOUZAS
RESPIRATORY
15
PARAMETERS
b
i
2-
16
I. KARPOUZAS
with R ( u) = (uu + b)/( u + c) (already calculated in the first model). (c) v,=0.1xV,+10~? (d) PA- PB=Rli+fC21j2. (e) u = (l/V”)(VVo,/Z). (f) ri=V,Xv. (g) From the Nelson relation [19], we have for active inspiration ( V > 6) )
and for active expiration F -= 4,
( V < VO) O.lV, t(H-
i
VI) i
.[I.,,(~ (t
is
the duration
-’
1+
-O.+
of the active muscular
,,,u,:yY,)] t = T/2).
intervention,
Finally,
we minimize with respect to V and ?, for three different and for different energy levels: of P,,m,
TABLE
Minimization
I&, (m3/sec) 7x10mh
2.2956062X X 10
’
2.29180176
X 10
’
2.28738402
x 10 - ’
2.32524736
X 10 - ’
1x10
5
1.6~10
’
7, + F< p4 “0,
I; (m’/sec)
V (ml) 2.2997699
1.3x1o-5
2
1 41 of J, = ~ xFxPz. “Jill
4x10mh
x 10 - ’
choices
1.X3045564X
J, minimum 10
‘l
1.5X370663
x 10
x 10
4
2.15365024
x 10 - ’
3.X515,1729 x 10
4
2.68176796
X 10 - ’
4.78177003
x 10
4
3.18729805
x 10
’
6.16970091
x 10
4
3.68639892
x 10
’
2.90949251
’
2.32188633
x 10
’
6.94885387
x 10
4
4.13351188
x 10
2.2 x 10
s
2.31563695
x 10 - ’
7.81293407
x 10
4
4.57516672
x 10 -
3 ’
2.5 x 10
’
2.31203X2
x10
’
8.58725017
x 10
4
5.01279471
x 10
’
2.31836198
x 10 -
’
9.46821698
X 10
4
5.4514009Xx
1.9x1o-5
2.Xx1o-5
10 - 3
3.1 x 10
5
2.3141772
x10
’
1.02639599
x 10
’
5.X77379
x10
3.4x10
5
2.30303901
x 10 - ’
1.08402466
x 10
’
6.31579803
x 10 ~’
3.7 x 10
5
2,31115599x10-”
1.16869872
x 10
’
6.72497018
x 10
3
’
RESPIRATORY
PARAMETERS
17 TABLE
Minimization
Vo, (mj/sec)
3 F, + P,
1 5, of J, = Pl.¶llZUxfxw
pa vo 1
Q (m’/sec)
V(m’)
J, (minimum)
4x10-h
2.29989048
x 10
’
1.81500039
x 10 ~’
2.22117128
x 10
3
7x10mh
2.34692882
x 10
’
3.27838212
x 10
4
3.01179715
x 10
3
lxlo~s
’
2.34020471
x 10
’
4.26212941
x 10
4
3.70269494
x 10
1.3 x 10 -s
2.34369566
x 10
’
5.17333835
x 10
4
4.37273415
x lo-’
1.6~10~~
2.3393659
x 10
’
6.03519838
x 10 - 4
5.00379237
x 10
’
1.9x10m5
2.33671971
x 10
3
6.81295906
x 10
4
5.62405485
x 10
’
2.2 x 10
5
2.3339151
x 10
’
7.59170216
x 10
4
6.23707801
x 10 ~’
2.5 x 10
5
2.36593099
x 10
3
8.94922015
x 10
4
6.79902393
x 10
’
2.8 x 10
5
2.37249631
x 10
’
9.8342186
x10
4
7.35656483
x 10
~’
3.1x10~s
2.37299403
x 10
3
1.04860321
x 10 ~’
10
’
3.4x10ms
2.36921457
x 10
’
1.12637408
x 10
’
8.41299352
x 10
3
3.7 x 10 - 5
2.37024181
x 10 -- ’
1.19006171
x 10
’
8.94007364
x 10
3
7.8X466229x
using
the global optimization method ALIENOR. The physiological bounds are as follows: The unknowns V, c, and the parameter PO’,,must belong to the respective intervals
5.1.
ALIENOR
METHOD
We have previously established a mathematical transformation which reduces the several variables of a function to a single one. This transformation is called ALIENOR [16, 171. The basic principle of ALIENOR is a special property of the Archimedean spiral: The curve corresponding to the polar equation R = e/a lies at a maximum distance 27r/a from any point in the plane (the distance being measured along the radius vector with 0 positive). If B takes positive and negative values, a double spiral is obtained. It lies at a maximum distance ?r/a from any point in the plane. Let
elcosel C=1.36+
~
50
and
elsinel D=O.l+p
50
.
If C d 6, let X = C; otherwise X = 1.36 + C - int( C)+ 3.64 X RND(~). If D G 2, let Y = D; otherwise Y = 0.1+ D - int( D)+ 0.9 X RND(9). Finally let V = X x10P3 and ri= Y x~O-~.
18
LKARPOUZAS
Then an optimization problem in terms of a single variable is obtained. A relative minimum was calculated using the probability technique. The next step is to explore along the curve from 0 = 0 looking for successively lower minima. The lowest minimum obtained was accepted as the required approximation, and was in fact close to the absolute minimum. The absolute minimum was then found by a subroutine using a local-variation method
P81. 5.2. RESULTS
The results obtained are compatible with the physiological variations during moderate muscular exercise, even though the intervention of the expiratory muscles was taken into account in the global energy expenditure of the ventilator-y system. 5.3. CONCLUSION
The energy minimization of the respiratory muscle has been taken as a regulation function determining, for a given I’,,, the ventilatory profile (V, V,, v) under conditions of moderate gaseous flow corresponding to moderate muscular exercise. The proposed regulation function is not incompatible with the physiological data for such exercise. Many methods are used in the determination of the local minima of the functions with several variables, but there is no effective method for finding a global optimum. ALIENOR is a very simple method which can reduce the minimization of a function of several variables to the search for the minimum of a function of only one. NOMENCLATURE A.
UNKNOWNS
V=volume of air V, = alveolar volume li= gaseous airflow v = respiratory frequency Vr = tidal volume VD= deadspace volume PC= compliance pressure Pa = alveolar pressure R = resistance of the bronchial and pulmonary passages to the airflow F=muscu!ar strength of any skeletal muscle at elongation I and initial contraction speed v p = the efficiency of the system
RESPIRATORY B
19
PARAMETERS
PARAMETER
oxygen flow
l&
=
C.
CONSTANTS
3 = 0.05 = used fraction of oxygen F,=lOO N =muscular strength of any skeletal muscle at elongation I, and contraction speed 0 H= 14.14 x 10 3 m3 (first model) or 12.96 X 10 _ 3 m3 (second model) = total volume of the fictitious heterogeneous “muscular and gaseous” cylinder P,=O = reference
pressure k,=2x107 kg me4 set = constant of the system V,=2.4X10P3 m3 = resting thoracopulmonary volume V,,=6~10-~ m3 = maximum volume of the system Pm, = 0.20 = constant of the system VR=1.36X10-3 m3 = residual volume REFERENCES 1
T. Julia Apter. Biosystems modeling, and .I. H. Milsum, Eds.). Chapter 5.
2
J. Mead, Control and respiratory frequency, J. Appl. Ph.vsiol. 15(3):325-336 (1960). J. Ph. Derenne, P. T. Macklem, and Ch. Roussos, State of the art, Part I, The respiratory muscles: Mechanics, control and pathophysiology. Amer. Reo. Resp. Dir. 118:119-133 (1978). R. Flandrois, J. Brune, and T. Wiesendanger, in Ph.wiologie Humame. Vol. 2, Simep Ed., 1976. R. M. Chemiack. Eprewes Fonctionnelles Respirutoires, Thtorie et Prutique. Doin. Paris. 1980. Ci. Ribier. Amelioration du residu dans la resolution de systemes lineaires au sens des moindres car&s, Mathematiques a l’usage du calculateur, 20, C.N.R.S., Feb. 1967. M. Laporte and J. Vignes, Algorithmes Numhiques, Antr!vse et Mise en Oeuore - T. I, Arithmkrrques .Qst>mes Linbires. Editions Technip, Paris, 1974. C. Fromageot. Organisation fonctionnelle de la mecanique ventilatoire pour optimisation de l’energie musculaire utilisee, diplame d’etudes et de recherches en biologie humaine, Univ. Rene Descartes, Facultt de Medecine de Paris-ouest, 1983.
3
4 5 6 7 8
in Biomedicd
Engineermg
Swrems
(M. Clynes
20 9
10 11 12 13
1. KARPOUZAS I. Karpouzas, Etude de la mecanique ventilatoire par application du principe de la minimisation de l’energie par des methodes mathematiques et numeriques. These de 3eme cycle, Univ. de Paris VI. 1983. J. L. Lions, Cows d’Anu!vse NumPrique. Ecole Polytechnique, Paris, 1973. B. Pchnitchny and Y. Daniline, MtSthodes Nun&yues dms /es Prohlhnm d ‘E~tremun~. Editions de Moscou. 1977. D. Wismer and R. Chattergy, Inrroductm to Nonlrneur Optimization, North Holland, Amsterdam, 1978. P. Berdot. Y. Cherruault, G. Korobenlik, P. Loridan. G. N&en, and F. Tonnelier, Contribution a la resolution de problemes d’optimisation par des methodes directcs, Compte Rendu. C.N.R.S., Paris, Dec. 1969.
14
Y. Cherruault, Recherches Biomathematiques. methodes fusion OFFILIB, 4X rue Gay Lussac, 75005 Paris, 1983.
15 16
J. Vignes, Etude et mise en oeuvre d’algorithmes de recherche d’une fonction a plusieurs variables, these de Doctorat d’Etat. Faculte des Sciences de Paris. 1969. Y. Cherruault and A. Guillez. Une methode pour la recherche du minimum global
17
d’une fonctionnelle, C. R. Acud. Sci. Paris ScG. I, 296, 24 Jan. 19X3. I. Karpouzas, Resolution numerique du principe d’optimisation applique a la mecanique yentilatoire. Cordeliers.
1X 19
presented at Seminaire 8 Nov. 1983.
de biomathematiques.
et exemples,
Institut
C.I.M.P.A.
Biomedical
Dif-
des
Y. Cherruault, Biomathematiques. presented at COB. Que sais-je?. P.U.F., Paris, 19X3. P. Nelson. Loglque des Neurotles er du SvstPme Neroeux, Ed. Doin. Paris, 197X.