# Optimal Determination of Respiratory Airflow Patterns Using a Nonlinear Multicompartment Model for a Lung Mechanics System

Authors: Hancao Li; Wassim M. Haddad
Journal: Computational and Mathematical Methods in Medicine
 (2012)

---

## Abstract

We develop optimal respiratory airflow patterns using a nonlinear multicompartment model for a lung mechanics system. Specifically, we use classical calculus of variations minimization techniques to derive an optimal airflow pattern for inspiratory and expiratory breathing cycles. The physiological interpretation of the optimality criteria used involves the minimization of work of breathing and lung volume acceleration for the inspiratory phase, and the minimization of the elastic potential energy and rapid airflow rate changes for the expiratory phase. Finally, we numerically integrate the resulting nonlinear two-point boundary value problems to determine the optimal airflow patterns over the inspiratory and expiratory breathing cycles.


---

## Body

## 1. Introduction

Respiratory failure, the inadequate exchange of carbon dioxide and oxygen by the lungs, is a common clinical problem in critical care medicine, and patients with respiratory failure frequently require support with mechanical ventilation while the underlying cause is identified and treated. The goal of mechanical ventilation is to ensure adequate ventilation, which involves a magnitude of gas exchange that leads to the desired blood level of carbon dioxide, and adequate oxygenation, which involves a blood concentration of oxygen that will ensure organ function. Achieving these goals is complicated by the fact that mechanical ventilation can actually cause acute lung injury, either by inflating the lungs to excessive volumes or by using excessive pressures to inflate the lungs. The challenge to mechanical ventilation is to produce the desired blood levels of carbon dioxide and oxygen without causing further acute lung injury.With the increasing availability of microchip technology, it has been possible to design partially automated mechanical ventilators with control algorithms for providing volume or pressure control [1–5]. More sophisticated fully automated model reference adaptive control algorithms for mechanical ventilation have also been recently developed [6, 7]. These algorithms require a reference model for identifying a clinically plausible breathing pattern. However, the respiratory lung models that have been presented in the medical and scientific literature have typically assumed homogenous lung function. For example, in analogy to a simple electrical circuit, the most common model has assumed that the lungs can be viewed as a single-compartment characterized by its compliance (the ratio of compartment volume to pressure) and the resistance to air flow into the compartment [8–10]. While a few investigators have considered two compartment models, reflecting the fact that there are two lungs (right and left), there has been little interest in more detailed models [11–13].Early work on the optimality of respiratory control mechanisms using simple homogenous lung models dealt with the frequency of breathing. In particular, the authors in [14–17] predicted the frequency of breathing by using a minimum work-rate criterion. This work involves a static optimization problem and assumes that the airflow pattern is a fixed sinusoidal function. The authors in [17, 18] developed optimality criteria for the prediction of the respiratory airflow pattern with fixed inspiratory and expiratory phases of a breathing cycle. These results were extended in [19] by considering a two-level hierarchical model for the control of breathing, in which the higher-level criterion determines values for the overall control variables of the optimal airflow pattern derived from the lower-level criteria, and the lower-level criteria determine the airflow pattern with the respiratory parameters chosen by minimizing the higher-level criterion.Although the problem for identifying optimal respiratory patterns has been addressed in the literature (see [14, 16–19] and the references therein), the models on which these respiratory control mechanisms have been identified are predicated on a single compartment lung model with constant respiratory parameters. However, the lungs, especially diseased lungs, are heterogeneous, both functionally and anatomically, and are comprised of many subunits, or compartments, that differ in their capacities for gas exchange. Realistic models should take this heterogeneity into account. In addition, the resistance to gas flow and the compliance of the lung units are not constant but rather vary with lung volume. This is particularly true for compliance. While more sophisticated models entail greater complexity, since the models are readily presented in the context of dynamical systems theory, sophisticated mathematical tools can be applied to their analysis. Compartmental lung models are described by a state vector, whose components are the volumes of the individual compartments.A key question that arises in the consideration of nonlinear multicompartment models is whether or not experimental data support a complex model. This question can be addressed by considering an analogy to pharmacokinetics. Specifically, the earliest pharmacokinetic models were typically one-compartment models. This reflected the challenges of sampling and drug assay. These models were adequate for quantifying drug disposition on a long time scale. For example, simple one-compartment models were adequate in describing the total clearance or volume of distribution. However, for even open-loop control of drug concentrations, the one compartment model was inadequate. More complex models (two- and three-compartment models) were needed that accounted for distribution as well as elimination processes (see [20] and the references therein).Similarly, for adaptive control of mechanical ventilation, that is, more advanced controller architectures than simple volume- or pressure-controlled ventilation, more elaborate models are needed, especially when accounting for nonlinear compliance and resistance and lung heterogeneity [6]. In the case of pharmacokinetics, the control algorithm can only be as complex as the data supports. This is also true for control of mechanical ventilation. Flow and pressure patterns in the airway are not simple waveforms, although clinicians to date have modeled them as such. There is considerable information embedded in these waveforms. The purpose of our work in this paper is to provide a mathematically rigorous and general framework developing optimal determination of respiratory airflow patterns using a nonlinear multicompartment model for a lung mechanics system. It is a an easy task to simplify this framework to be congruent with the granularity of the data. The reverse process, however, is not possible without the development of a general framework.In this paper, we extend the work of [17, 18] to develop optimal respiratory airflow patterns using a nonlinear multicompartment model for a lung mechanics system. (The usage of the word optimal throughout the paper refers to an optimal solution of the calculus of variations problems addressed in the paper and not an optimal breathing pattern in the sense of respiratory physiology.) First, we extend the linear multicompartment lung model given in [6] to address system model nonlinearities. Then, we extend the performance functionals developed in [17, 18] for the inspiratory and expiratory breathing cycles to derive an optimal airflow pattern using classical calculus of variations techniques. In particular, the physiological interpretation of the optimality criteria involves the minimization of work of breathing and lung volume acceleration for the inspiratory breathing phase, and the minimization of the elastic potential energy and rapid airflow rate changes for the expiratory breathing phase. Finally, we numerically integrate the resulting nonlinear two-point boundary value problems to determine the optimal airflow patterns over the inspiratory and expiratory breathing cycles.The notation used in this paper is fairly standard. Specifically,ℝn denotes the set of n×1 real column vectors, and ℝn×m, denotes the set of n×m real matrices. For x∈ℝn we write x≥≥0 (resp., x≫0) to indicate that every component of x is nonnegative (resp., positive). In this case, we say that x is nonnegative or positive, respectively. Likewise, A∈ℝn×m is nonnegative or positive if every entry of A is nonnegative or positive. (In this paper, it is important to distinguish between a square nonnegative (resp., positive) matrix and a nonnegative-definite (resp., positive-definite) matrix.) Furthermore, ℝ̅+n and ℝ+n denote the nonnegative and positive orthants of ℝn, that is, if x∈ℝn, then x∈ℝ̅+n and x∈ℝ+n are equivalent, respectively, to x≥≥0 and x≫0. Finally, we write (·)T to denote transpose, (·)′ to denote Frèchet derivative, and δx to denote the first variation of the function x.


## 2. A Nonlinear Multicompartment Modelfor Respiratory Dynamics

In this section, we extend the linear multicompartment lung model of [6] to develop a nonlinear model for the dynamic behavior of a multicompartment respiratory system in response to an arbitrary applied inspiratory pressure. Here, we assume that the bronchial tree has a dichotomy architecture [21]; that is, in every generation each airway unit branches into two airway units of the subsequent generation. In addition, we assume that the lung compliance is a nonlinear function of lung volume.First, for simplicity of exposition, we consider a single-compartment lung model as shown in Figure1. In this model, the lungs are represented as a single lung unit with nonlinear compliance c(x) connected to a pressure source by an airway unit with resistance (to airflow) of R. At time t=0, a driving pressure pin(t) is applied to the opening of the parent airway, where pin(t) is generated by the respiratory muscles or a mechanical ventilator. This pressure is applied over the time interval 0≤t≤Tin, which is the inspiratory part of the breathing cycle. At time t=Tin, the applied airway pressure is released and expiration takes place passively, that is, the external pressure is only the atmospheric pressure pex(t) during the time interval Tin≤t≤Tin+Tex, where Tex is the duration of expiration.Figure 1
Single-compartment lung model.The state equation for inspiration (inflation of lung) is given by(1)Rinẋ(t)+1cin(x)x(t)=pin(t),x(0)=x0in,0≤t≤Tin,
where x(t)∈ℝ, t≥0, is the lung volume, Rin∈ℝ is the resistance to airflow during the inspiration period, cin:ℝ→ℝ+ is a nonlinear function defining the lung compliance at inspiration, and x0in∈ℝ+ is the lung volume at the start of the inspiration and serves as the system initial condition. Equation (1) is simply a pressure balance equation where the driving pressure pin(t),0≤t≤Tin, applied to the compartment is proportional to the volume of the compartment via the compliance and the rate of change of the compartmental volume via the resistance. We assume that expiration is passive due to the elastic stretch of the lung unit. During the expiration process, the state equation is given by(2)Rexẋ(t)+1cex(x)x(t)=pex(t),x(Tin)=x0ex,Tin≤t≤Tin+Tex,
where x(t)∈ℝ, t≥0, is the lung volume, Rex∈ℝ is the resistance to air flow during the expiration period, cex:ℝ→ℝ+ is a nonlinear function defining the lung compliance at expiration, and x0ex∈ℝ+ is the lung volume at the start of expiration.Next, we develop the state equations for inspiration and expiration for a2n-compartment model, where n≥0. In this model, the lungs are represented as 2n lung units which are connected to the pressure source by n generations of airway units, where each airway is divided into two airways of the subsequent generation leading to 2n compartments (see Figure 2 for a four-compartment model).Figure 2
Four-compartment lung model.Letxi,i=1,2,…,2n, denote the lung volume in the ith compartment, let ciin(xi)(resp.,ciex(xi)), i=1,2,…,2n, denote the compliance at inspiration (resp., expiration) of each compartment as a nonlinear function of the volume of ith compartment, and let Rj,iin (resp., Rj,iex), i=1,2,…,2j,j=0,…,n, denote the resistance (to air flow) of the ith airway in the jth generation during the inspiration (resp., expiration) period with R01in (resp., R01ex) denoting the inspiration (resp., expiration) of the parent (i.e., 0th generation) airway.As in the single-compartment model, we assume that a pressure ofpin(t),t≥0, is generated (by the inspiratory muscles) or applied (by a mechanical ventilator) during inspiration. Now, the state equations for inspiration are given by(3)Rn,iinẋi(t)+1ciin(xi(t))xi(t)+∑j=0n-1Rj,kjin∑l=(kj-1)2n-j+1kj2n-jẋl(t)=pin(t),xi(0)=xi0in,0≤t≤Tin,i=1,2,…,2n,
where ciin(xi),i=1,2,…,2n, are nonlinear functions of xi, i=1,2,…,2n, given by [22](4)ciin(xi)≜{ai1in+bi1inxi,if0≤xi≤xi1in,ai2in,ifxi1in≤xi≤xi2in,ai3in+bi3inxi,ifxi2in≤xi≤VTi,i=1,…,2n,
where aijin, j=1,2,3, and bijin, j=1,3, are model parameters with bi1in>0 and bi3in<0, xijin, j=1,2, are volume ranges wherein the compliance is constant, VTi denotes tidal volume, and(5)kj=⌊kj+1-12⌋+1,j=0,…,n-1,kn=i,
where ⌊q⌋ denotes the floor function which gives the largest integer less than or equal to the positive number q. Figure 3 shows a typical piecewise linear compliance function for inspiration. A similar compliance representation holds for expiration and is also shown in Figure 3.Figure 3
Typical inspiration and expiration compliance functions as function of compartmental volumes.To further elucidate the inspiration state equation for a2n-compartment model, consider the four-compartment model shown in Figure 2 corresponding to a two-generation lung model. Let xi,i=1,2,3,4, denote the compartmental volumes. Now, the pressure (1/ciin(xi(t)))xi(t) due to the compliance in ith compartment will be equal to the difference between the driving pressure and the resistance to air flow at every airway in the path leading from the pressure source to the ith compartment. In particular, for i=3 (see Figure 2),(6)1c3in(x3(t))x3(t)=pin(t)-R0,1in[ẋ1(t)+ẋ2(t)+ẋ3(t)+ẋ4(t)]-R1,2in[ẋ3(t)+ẋ4(t)]-R2,3inẋ3(t),or, equivalently,(7)R2,3inẋ3(t)+R1,2in[ẋ3(t)+ẋ4(t)]+R0,1in[ẋ1(t)+ẋ2(t)+ẋ3(t)+ẋ4(t)]+1c3in(x3(t))x3(t)=pin(t).Next, we consider the state equation for the expiration process. As in the single-compartment model, we assume that the expiration process is passive and the external pressure applied ispex(t),t≥0. Following an identical procedure as in the inspiration case, we obtain the state equation for expiration as(8)Rn,iexẋi(t)+∑j=0n-1Rj,kjex∑l=(kj-1)2n-j+1kj2n-jẋl(t)+1ciex(xi(t))xi(t)=pex(t),xi(Tin)=xi0ex,Tin≤t≤Tex+Tin,i=1,2,…,2n,
where(9)ciex(xi)≜{ai1ex+bi1exxi,ai2ex,ai3ex+bi3exxi,if0≤xi≤xi1ex,ifxi1ex≤xi≤xi2ex,ifxi2ex≤xi≤VTi,i=1,…,2n,aijex, j=1,2,3, and bijex, j=1,3, are model parameters with bi1ex>0 and bi3ex<0, xijex, j=1,2, are volume ranges wherein the compliance is constant, and kj is given by (5).Next, we provide a smooth (i.e.,C∞) characterization of the nonlinear compliance using sigmoidal functions [23]. Specifically, for inspiration, ciin(xi) can be approximated as(10)ciin(xi)≈ai2in(Sa,b(β)(xi)-Sc,d(β)(xi)),i=1,…,2n,
where a=-ai1in/bi1in, b=(ai2in/bi1in)+a, c=-ai3in/bi3in, d=ai2in/bi3in+c, Sa,b(β)(xi)≜1/(b-a)ln(σb(-β)(xi)/σa(-β)(xi))1/β with σb(-β)(xi)≜1/(1+e-β(xi-a)), and β>0 is an approximation parameter. Figure 4 shows the smoothed approximation of the piecewise linear compliance function ciin(xi). A similar approximation holds for ciex(xi) and is also shown in Figure 4.Figure 4
Original and the smoothed compliance functions,β=30.Finally, we rewrite the state equations (3) and (8) for inspiration and expiration, respectively, in vector-matrix state space form. Specifically, define the state vector x≜[x1,x2,…,x2n]T, where xi denotes the lung volume of the ith compartment. Now, the state equations (3) for inspiration can be rewritten as(11)Rinẋ(t)+Cin(x(t))x(t)=pin(t)e,x(0)=x0in,0≤t≤Tin,
where e≜[1,…,1]T denotes the one vector of order 2n, Cin(x) is a diagonal matrix function given by(12)Cin(x)≜diag⁡[1c1in(x1),…,1c2nin(x2n)],(13)Rin≜∑j=0n∑k=12jRj,kinZj,kZj,kT,
where Zj,k∈ℝ2n is such that the lth element of Zj,k is 1 for all l = (k-1)2n-j+1,(k-1)2n-j+2,…,k2n-j, k=1,…,2j, j=0,1,…,n, and zero elsewhere.Similarly, the state equation (8) for expiration can be rewritten as(14)Rexẋ(t)+Cex(x(t))x(t)=pex(t)e,x(Tin)=x0ex,Tin≤t≤Tex+Tin,
where(15)Cex(x)≜diag⁡[1c1ex(x1),…,1c2nex(x2n)],(16)Rex≜∑j=0n∑k=12jRj,kexZj,kZj,kT.
Finally, it follows from [6, Proposition 4.1] that Rin and Rex are positive definite and, hence, Rin and Rex are invertible matrices.


## 3. Optimal Determination of Inspiratory and Expiratory Airflow in Breathing

In this section, we use the respiratory dynamical system characterized by (11) and (14) to develop an optimal model for predicting airflow patterns in breathing. The optimization criteria used allows for the minimization of oxygen expenditure of the respiratory muscles as well as rapid changes in the lung volume flow rate. The oxygen consumption of the lung muscles is mainly due to the work carried out by the respiratory muscles to overcome the resistive forces and stretch the lung and chest wall. In [24], this work is defined as PV, where P is the pressure driving inflation and V is the lung unit volume. The efficiency of gas exchange in the lungs is related to the volume acceleration, since rapid changes in lung volume can cause discomfort and inefficacy of muscular contraction and control. Moreover, high-volume acceleration can result in overexpansion of the lung resulting in lung tissue rupture as well as excessive work of breathing with subsequent ventilatory muscle fatigue.In the ensuing discussion, we assume that the inspiration process starts from a given initial statex0in and is followed by the expiration process where its initial state will be the final state of the inspiration. An inspiration followed by an expiration is called a single breathing cycle. Furthermore, we assume that each breathing cycle is followed by another breathing cycle where the initial condition for the latter breathing cycle is the final state of the former breathing cycle. Since the respiratory process is periodic, we need only focus on one breathing cycle.The next result gives the optimal solutionx*(t),0≤t≤Tin, for the inspiratory airflow breathing pattern using classical calculus of variations techniques.Theorem 1.
Consider the nonlinear system model for inspiration given by (11). Let the optimal air volume x*(t),0≤t≤Tin, be given by the solution to the minimization problem
(17)Jin(x)=∫0Tin[ẍT(t)ẍ(t)+α1pin(t)eTẋ(t)]dt,α1≥0,
subject to the natural boundary conditions
(18)x(0)=V0,ẋ(0)=0,(19)x(Tin)=V0+VT,ẋ(Tin)=0,
where V0∈ℝ2n is the end expiratory volume and VT∈ℝ2n is the tidal volume. If α1>0, then x*(t),0≤t≤Tin, is given by
(20)x*(t)=d1+d2t+exp⁡(α1Rin1/2t)d3+exp⁡(-α1Rin1/2t)d4,t≥0,
and if α1=0, then
(21)x*(t)=d1+d2t+d3t2+d4t3,t≥0,
where d1, d2, d3, and d4∈ℝ2n are constant vectors determined by the boundary conditions (18) and (19), and Rin1/2 denotes the (unique) positive-definite square root of Rin.Proof.
First, note thatpin(t)e,0≤t≤Tin, in (17) can be eliminated using the state equation (11). Thus, the integrand of the performance criterion (17) can be written as
(22)Lin(x(t),ẋ(t),ẍ(t))=ẍT(t)ẍ(t)+α1[Rinẋ(t)+Cin(x(t))x(t)]Tẋ(t)=ẍT(t)ẍ(t)+α1[ẋT(t)Rinẋ(t)+xT(t)Cin(x)ẋ(t)],α1≥0.
The first variation of the performance criterion 𝒥in(x) is given by
(23)δJin(x*,δx)=∫0TinδLin(x*(t),ẋ*(t),ẍ*(t))dt=∫0Tin{(∂Lin∂x)δx(t)+(∂Lin∂ẋ)δẋ(t)+(∂Lin∂ẍ)δẍ(t)}dt=[∂Lin∂ẍδẋ+(∂Lin∂ẋ-d2dt2∂Lin∂ẍ)δx]0Tin+∫0Tin{(∂Lin∂x)-ddt(∂Lin∂ẋ)+ddt(∂Lin∂ẍ)}δx(t)dt.
Using the boundary conditions (18) and (19), it follows that δx(0)=δx(Tin)=δẋ(0)=δẋ(Tin)=0. Now, since Tin is fixed, it follows from the fundamental theorem of the calculus of variations that the variation of 𝒥in(x) must vanish on x*; that is, the extremals optimizing the performance criterion 𝒥in(x) satisfy the Euler-Lagrange equation
(24)(∂Lin∂x)T-ddt(∂Lin∂ẋ)T+d2dt2(∂Lin∂ẍ)T=0.
Next, usingCin(x) given by (12),
(25)(∂Lin∂x)T=α1Cin(x(t))ẋ(t)+α1Cin′(x(t))Ẋ(t)x(t),0≤t≤Tin,(∂Lin∂ẋ)T=2α1Rinẋ(t)+α1Cin(x(t))x(t),0≤t≤Tin,(∂Lin∂ẍ)T=2ẍ(t),0≤t≤Tin,
where Cin′(x(t))≜diag[(∂/∂xi)(1/(ciin(xi(t))))] and Ẋ(t)≜diag[ẋi(t)],i=1,…,2n. Thus, (24) yields the fourth-order differential equation
(26)x(4)(t)-α1Rinx(2)(t)=0,0≤t≤Tin,
where x(n)(t)≜(dnx(t)/dtn), with boundary conditions given in (18) and (19). Now, using standard analysis techniques, the solution x(t),0≤t≤Tin, to (26) satisfies (20) if α1>0 and (21) if α1=0.Remark 2.
The vectorsd1, d2, d3, and d4 in Theorem 1 can be uniquely determined using the four boundary conditions given by (18) and (19). Specifically, if α1=0, it can be shown that d1=V0, d2=0,  d3=(3/Tin2)VT,  and d4=-(2/Tin3)VT. Hence, in this case, ẋ(t) = d2+2d3t+3d4t2 = (6t/(Tin2))VT(1-(t/Tin))≥≥0, 0≤t≤Tin, which implies that the solution x*(t),0≤t≤Tin, to (26) is increasing during inspiration, and hence, V0i≤xi*(t)≤V0i+VTi,i=1,…,2n, where V0i,xi(t) and VTi are the ith components of V0,x(t), and VT, respectively. A similar result holds for the case where α1>0.
Next, we give the optimal solutionx*(t),Tin≤t≤Tin+Tex, for the expiratory airflow breathing pattern.Theorem 3.
Consider the nonlinear system model for expiration given by (14). Let the optimal solution x*(t),Tin≤t≤Tin+Tex, be given by the solution to the minimization problem
(27)Jex(x)=∫TinTin+Tex[ẍT(t)ẍ(t)+α2pex2(t)eTe]dt,α2≥0,
subject to the natural boundary conditions
(28)x(Tin)=V0+VT,ẋ(Tin)=0,(29)x(Tin+Tex)=V0,ẋ(Tin+Tex)=0.
If α2>0, then x*(t),Tin≤t≤Tin+Tex, satisfies
(30)x(4)(t)-α2Rex2x(2)(t)+α2Cex2(x)x(t)+α2[Cex(x)Rexẋ(t)-RexCex(x)ẋ(t)+X(t)Cex′(x)Rexẋ(t)-RexCex′(x)X(t)ẋ(t)+X(t)Cex′(x)Cex(x)x(t)]=0,
where X(t)≜diag[xi(t)] and Cex′(x)≜diag[(∂/∂xi)(1/ciex(xi))],i=1,…,2n, and if α2=0, then
(31)x*(t)=d1+d2t+d3t2+d4t3,t≥0,
where d1, d2, d3, and d4∈ℝ2n are constant vectors determined by the four boundary conditions (28) and (29).Proof.
Using (14), the integrand of the performance criterion (27) can be written as
(32)Lex(x(t),ẋ(t),ẍ(t))=ẍT(t)ẍ(t)+α2(pex(t)e)T(pex(t)e)=ẍT(t)ẍ(t)+α2[Rexẋ(t)+Cex(x(t))x(t)]T×[Rexẋ(t)+Cex(x(t))x(t)]=ẍT(t)ẍ(t)+α2[ẋT(t)Rex2ẋ(t)+xT(t)×Cex2(x(t))x(t)+2ẋT(t)RexCex(x(t))x(t)],α2>0.
Thus, the variation of 𝒥ex(x) on an extremal solution gives
(33)δJex(x*,δx)=∫TinTin+TexδLex(x*(t),ẋ*(t),ẍ*(t))dt=∫TinTin+Tex{(∂Lex∂x)δx(t)+(∂Lex∂ẋ)δẋ(t)+(∂Lex∂ẍ)δẍ(t)}dt=[∂Lex∂ẍδẋ+(∂Lex∂ẋ-ddt∂Lex∂ẍ)δx]TinTin+Tex+∫0Tex{d2dt2(∂Lex∂x)-ddt(∂Lex∂ẋ)+d2dt2(∂Lex∂ẍ)d2dt2}δx(t)dt=0.
Using the boundary conditions (28) and (29), it follows that δx(Tin)=δx(Tin+Tex)=δẋ(Tin)=δẋ(Tin+Tex)=0. Hence, the extremals optimizing the performance criterion 𝒥ex(x) satisfy the Euler-Lagrange equation
(34)(∂Lex∂x)T-ddt(∂Lex∂ẋ)T+d2dt2(∂Lex∂ẍ)T=0.
Now, usingCex(x) given by (15),
(35)(∂Lex∂x)T=α2[2Cex2(x(t))x(t)+2Cex(x(t))Rexẋ(t)+2X(t)Cex′(x(t))Rexẋ(t)+2X(t)Cex′(x(t))Cex(x(t))x(t)],Tin≤t≤Tin+Tex,(∂Lex∂ẋ)T=α2[2Rex2ẋ(t)+2RexCex(x(t))x(t)],Tin≤t≤Tin+Tex,(∂Lex∂ẍ)T=2ẍ(t),Tin≤t≤Tin+Tex,
which yields (30). Finally, in the case where α2=0, (30) collapses to x(4)(t)=0,Tin≤t≤Tin+Tex, which satisfies (31).Remark 4.
In the case whereα2=0, the vectors d1, d2, d3, and d4 in Theorem 3 can be uniquely determined using the four boundary conditions (28) and (29). In particular, d1 = V0+VT+3βTin2TexVT+2βTin3VT,  d2 = -β(6Tin2VT+6TexTinVT), d3 = β(3TexVT+6TinVT), and d4=-2βVT,  where β=1/(3Tex3+12Tex2Tin+12TexTin2+4Tin3). Hence,  ẋ(t)=d2+2d3t+3d4t2=-6βVTt(Tin+Tex-t)-6βVTt(t-Tin)≤≤0, Tin≤t≤Tin+Tex, which implies that the solution x*(t),Tin≤t≤Tin+Tex, is decreasing during expiration, and hence, V0i≤xi*(t)≤V0i+VTi, i=1,…,2n. The case where α2>0 involves the solution to (30), and hence, we have been unable to show that x*(t),Tin≤t≤Tin+Tex, is decreasing during expiration analytically. However, this has been verified numerically.Remark 5.
If optimal solutions to Theorems1 and 3 exist, then the optimal respiratory airflow patterns and their corresponding driving pressures can be computed using the lung mechanics model developed in Section 2. The input signal to this model can then be used as the driving pressure of a mechanical ventilator needed to achieve the optimal respiratory airflow pattern.
The physiological interpretations of the performance criteria for inspiration and expiration used in Theorems1 and 3 are slightly different. In particular, the inspiratory criterion 𝒥in(x) involves a weighted sum of squares of the lung volume acceleration and the mechanical work performed by the inspiratory muscles. Alternatively, during the expiratory phase, the respiratory muscles remain active in the beginning of expiration since they continue their action by opposing expiration and hence consume oxygen thereby performing negative work. Thus, mechanical work alone is not a satisfactory criterion for describing control of breathing at rest. As in [25], we assume that oxygen consumption of expiration correlates with the integral square of the driving pressure. This assumption is supported in [26] which shows that an index of average respiratory pressure can predict the total oxygen cost of breathing. Hence, instead of mechanical work, we use the integral square of the applied pressure in the expiratory criterion 𝒥ex(x), which corresponds to minimizing the mean standard potential energy in the lung.
It can be seen that the optimal solutionsx*(t),t≥0, depend on the variables Tin, Tex, V0, and VT through the boundary conditions. Moreover, the nonlinearities in (30) are due to nonlinearities in the lung compliance Cex(x), which make analytical solutions to (30) difficult to obtain. It is interesting to note that although the optimal solutions x*(t),Tin≤t≤Tin+Tex, to (30) during the expiration phase depend on the nonlinear compliance of Cex(x), the optimal solutions x*(t),0≤t≤Tin, to (26) during the inspiration phase are independent of the nonlinear system compliance Cin(x). In the case where n=0 (i.e., a single-lung-compartment model), x(t)∈ℝ, Rex∈ℝ, and Cex(x)=Cex are constants, (30) reduces to
(36)x(4)(t)-α2Rex2x(2)(t)+α2Cex2x(t)=0.
This case is extensively discussed in [25] wherein the authors characterize four different solutions to (36) corresponding to α2=0, 0<α2<4Cex2/Rex4, α2=4Cex2/Rex4, and α2>4Cex2/Rex4.


## 4. Numerical Determination of Optimal Volume Trajectories

The optimal volume trajectories formulated in Section3 result in two-point nonlinear boundary-value problems. Numerical methods for solving such problems include shooting methods [27] and steepest descent methods [28]. In this section, we use the collocation method implemented by bvp4c in MATLAB [29] to numerically integrate the differential equations (26) and (30) to obtain the optimal volume trajectory x*(t),t≥0.For our simulations, we first consider a two-compartment lung model and use the values for the lung compliance found in [22]. In particular, we set ai1in=0.018ℓ/cm H2O, bi1in=0.0233, ai2in=0.025ℓ/cm H2O, ai3in=0.2532ℓ/cm H2O, bi3in=-0.01, xi1in=0.3ℓ, xi2in=0.48ℓ, ai1ex=0.02ℓ/cm H2O, bi1ex=0.078, ai2ex=0.038ℓ/cm H2O, ai3ex=0.1025ℓ/cm H2O, bi3ex=-0.15, xi1ex=0.23ℓ, xi2ex=0.43ℓ, and i=1,2. Here, we assume that the bronchial tree has a dichotomy structure (see Section 2). The airway resistance varies with the branch generation, and typical values can be found in [30]. Furthermore, the expiratory resistance will be higher than the inspiratory resistance by a factor 2 to 3. Here, we assume that the factor is 2.5.For our simulation, we assume that the inspiration timeTin=2sec and the expiration time Tex=3sec. The two weighting parameters α1 and α2 differ from person to person. Nominal values for the weighting parameters are α1=2.0l/sec3cmH2O and α2=0.1l2/sec4cmH2O, which correspond to spontaneous breathing at rest [25]. Figure 5 shows the optimal air volume eTx*(t),t≥0, and the optimal airflow rate eTẋ*(t),t≥0, given by the two-point nonlinear boundary-value problems (24) and (34). Note that the airflow curve for inspiration is symmetric, since the nonlinearities in Cin(x) do not appear in (26). However, x*(t),t≥0, obtained using (30) during expiration involves Cex(x), and hence, the airflow curve is asymmetric.Figure 5
Volume and airflow rate patterns for the total lung compartments.Figure6 shows the sensitivity of the optimal volume and airflow rate patterns to changes in the parameters α1 and α2. As can be seen from the figure, the inspiratory airflow rate is symmetric and the maximum value of the airflow rate decreases as a function of increasing α1. Furthermore, the asymmetric pattern of the expiratory airflow rate reflects the fact that the minimum value becomes steeper with increasing α2. Specifically, if we set the weighting parameter α2=0, it follows from (30) that the airflow curve for the expiration is given by a parabolic arc. The airflow patterns in Figure 6 exhibit typical respiratory characteristics observed in spontaneous breathing, that is, the inspiratory airflow rate curve is relatively flat and the expiratory airflow rate waveform is asymmetric with an initial trough, and quite similar to “real” airflow signals [31].Figure 6
Volume and airflow rate patterns for differentα1’s and α2’s.Figure7 shows the driving pressure generated by the respiratory muscles using the optimal air volume eTx*(t),t≥0. Figure 8 compares the optimal air volume trajectory eTx*(t),t≥0, with a nonoptimal air volume trajectory eTx(t),t≥0, generated by the linear pressure pin(t)=20t+5cm H2O, t∈[0,Tin], and pex(t)=0cm H2O, t∈[Tin,Tin+Tex], [6]. Note that eTx*(t),t≥0, switches between the end expiratory level eTV0=0.2l and the tidal volume eTVT=1.2l. Figure 9 shows the phase portrait of the optimal trajectories x1*(t) and x2*(t) and suboptimal trajectories x1(t) and x2(t). Note that both sets of trajectories asymptotically converge to a limit cycle, with the optimal solutions satisfying the boundary conditions given in (18), (19), (28), and (29). Figure 10 compares the value of the total performance criterion generated by the optimal air volume with the value of the total performance criterion generated by the nonoptimal air volume.Figure 7
Pressure generated by optimal solution.Figure 8
Optimal volumeeTx*(t) and nonoptimal volume eTx(t) versus time.Figure 9
Phase portrait forx1*(t) versus x2*(t) and x1(t) versus x2(t).Figure 10
Performance criterion comparison versus time.Finally, Figure11 shows the optimal air volume trajectories for a four-compartment model with each air volume trajectory satisfying the boundary conditions given in (18), (19), (28), and (29). For this simulation, the compliance parameters are taken to be identical to those used for the two-compartment model with i=1,2,3,4, and the values for airway resistances are generated using the results of [30].Figure 11
Optimal volumex*(t) versus time for a four-compartmental model.


## 5. Conclusion and Directions for Future Work

In this paper, we developed an optimal respiratory air flow pattern using a nonlinear multicompartment model for a lung mechanics system. The determination of the optimal air volume trajectories is derived using classical calculus of variations techniques and involves optimization criteria that account for oxygen expenditure of the respiratory lung muscles, lung volume acceleration, and elastic potential energy of the lung. Future work will include the development of multivariable and adaptive control algorithms that will utilize these models within a model reference control architecture for fully automating mechanical ventilation to ensure adequate ventilation and oxygenation for critical care patients in intensive care units.Since sedation in intensive care units is often administered to prevent the patient from fighting the ventilator, it seems plausible to use respiratory parameters as a performance variable for closed-loop control. Calculation of patient work of breathing requires measurement of a patient-generated pressure/volume loop or work of breathing. Since work of breathing can be measured using a commercially available esophageal balloon [32], work of breathing can serve as a performance variable for closed-loop control of sedation. Furthermore, patient-ventilator dyssynchrony can be identified by analysis of pressure/flow wave forms [33].Closed-loop control algorithms can use either work of breathing as measured by an esophageal balloon or patient respiratory rate as a performance variable for closed-loop control of sedation. The need for optimal control algorithms is necessary for achieving a target performance value while satisfying certain constraints. For example, we could seek to design a control algorithm that seeks to minimize the patient respiratory rate (above the set ventilator rate) but which does not result in hypotension. This requires the development of a constrained optimal control framework that seeks to minimize a given performance measure (e.g., patient respiratory rate) within a class of fixed-architecture controllers satisfying internal controller constraints (e.g., controller order, control signal nonnegativity, etc.) as well as system constraints (e.g., blood pressure, system state nonnegativity, etc.). The results in the present paper can serve as a starting point for developing multivariable controllers for mechanical ventilation of critically ill patients.



---

*Source: 165946-2012-06-08.xml*