Limit Cycle Stability Analysis of a Multi-Compartment Model for a Pressure-Limited Respirator and Lu

Limit Cycle Stability Analysis of a Multi-Compartment Model for a Pressure-Limited Respirator and Lung Mechanics SystemVijaySekhar Chellaboina∗,Wassim M.Haddad†,James M.Bailey‡,Hancao Li∗∗Mechanical,Aerospace,and Biomedical Engineering,The University of Tennessee,Knoxville,TN37996-2210†School of Aerospace Engineering,Georgia Institute of Technology,Atlanta,GA30332-0150‡Department of Anesthesiology,Northeast Georgia Medical Center,Gainesville,GA30503Abstract—Acute respiratory failure due to infection,trauma, or major surgery is one of the most common problems encoun-tered in intensive care units and mechanical ventilation is the mainstay of supportive therapy for such patients.In this paper, we develop a general mathematical model for the dynamic be-havior of a multi-compartment respiratory system in response to an arbitrary applied inspiratory pressure.Specifically,we use compartmental dynamical system theory to model and analyze the dynamics of a pressure-limited respirator and lung mechanics system,and show that the periodic orbit generated by this system is globally asymptotically stable.Furthermore, we show that the individual compartmental volumes,and hence the total lung volume,converge to steady-state end-inspiratory and end-expiratory values.I.I NTRODUCTIONAcute respiratory failure due to infection,trauma,and major surgery is one of the most common problems encoun-tered in intensive care units and mechanical ventilation is the mainstay of supportive therapy for such patients.Numer-ous mathematical models of respiratory function have been developed in the hope of better understanding pulmonary function and the process of mechanical ventilation[1–5]. However,the 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 airflow into the compartment[1],[2], [5].While a few investigators have considered two compart-ment models,reflecting the fact that there are two lungs (right and left),there has been little interest in more detailed models[6–8].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.While more sophisti-cated 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 partmental lung models are described by a state vector,whose components are the volumes of the individual compartments.One interesting and important question is the stability,in the sense of dynamical systems theory,of the model.This research was supported in part by NSF under Grants ECS-0551947and ECS-0601311,and AFOSR under Grant FA9550-06-1-0240.For a simple one compartment model,it is easy to demon-strate that the model exhibits an asymptotically stable limit cycle behavior.And indeed,in clinical practice it appears that the total lung volume converges to the steady-state end-inspiratory and end-expiratory values after the institution of mechanical ventilation.However,a more subtle question for a multi-compartment lung model is whether the volumes in the individual compartments could be unstable,even when the total volume of the lung(the sum of all the compartment volumes)converges to a steady-state value. That is,is it possible that individual compartment volumes oscillate or even demonstrate chaotic behavior while the total lung volume is stable?This question has interesting clinical implications as there is also heterogeneity in the amount of bloodflowing to individual subunits of the lung.If there is significant disparity in the ratio of ventilation(reflected in the compartment volume)to bloodflow,gas exchange is impaired,resulting in decreases in the oxygen or increases in the carbon diox-ide content of blood,which is a serious clinical problem. Instability of the compartment volumes could be reflected in unstable measures of basic pulmonary function,such as oxygen or carbon dioxide levels in the blood.In this paper, we develop a generalized multi-compartment lung model and subsequently analyze its stability properties.Specifically,we use compartmental dynamical system theory to model and analyze the dynamics of a pressure-limited respirator and lung mechanics system,and show that the periodic orbit generated by this system is globally asymptotically stable. Furthermore,we show that the individual compartmental volumes,and hence the total lung volume,converge to steady-state end-inspiratory and end-expiratory values.II.N OTATION AND M ATHEMATICAL P RELIMINARIES In this section,we introduce notation,several definitions, and some key results that are necessary for developing the main results of this paper.Specifically,for x∈R 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∈R n×m is nonnegative or positive if every entry of A is nonnegative or positive,respectively,which is written as A≥≥0or A>>0,respectively.Furthermore,for A∈R n×n we write A≥0(resp.,A>0)to indicate that A is a nonnegative-definite(resp.,positive-definite)matrix.In addition,(·)T denotes transpose and(·)−1denotes inverse.LetR n+and x∈R n+are equivalent, respectively,to x≥≥0and x>>0.Finally,let e n∈R n denote the ones vector of order n,that is,e n=[1,···,1]T; if the order of e n is clear from context we simply write e for e n.Proceedings of the 2007 American Control ConferenceMarriott Marquis Hotel at Times SquareNew York City, USA, July 11-13, 2007WeC18.4The following definitions introduce the notions of essen-tially nonnegative,compartmental,and strictly ultrametric matrices.Definition 2.1([9],[10]):Let A ∈R n ×n .A is essen-tially nonnegative if A (i,j )≥0,i,j =1,...,n ,i =j .A is compartmental if A is essentially nonnegative and A T e ≤≤0.Definition 2.2([11]):Let A ∈R n ×n be such that A ≥≥0.A is strictly ultrametric if A is symmetric,A (i,i )>max {A (i,k ):k =1,...,n,k =i },i =1,...,n ,and A (i,j )≥min {A (i,k ),A (k,j )},k =1,...,n ,i,j =1,...,n,i =j .The following lemmas and propositions are key in estab-lishing the main results of the paper.Lemma 2.1([9],[10]):Let A ∈R n ×n .Then A is es-sentially nonnegative if and only if e At is nonnegative for all t ≥0.Proposition 2.1:The following statements hold:i )Let λ1,λ2≥0be such that λ1+λ2>0and let A 1,A 2∈R n ×n be strictly ultrametric.Then λ1A 1+λ2A 2is strictly ultrametric.ii )Let x ∈R n be such that x i =0or 1,i =1,...,n,and let P ∈R n ×n be a positive diagonal matrix.Then P +xx T is a strictly ultrametric matrix.Lemma 2.2([11]):Let A ∈R n ×n be such that A ≥≥0.If A is strictly ultrametric,then −A −1is essentially nonnegative and A −1e ≥≥0.Proposition 2.2:Let A ∈R n ×n and assume that there exists an n ×n matrix P >0such thatA T P +P A <0.(1)Then e A TP e A <P.Remark 2.1:It is well known that A is Hurwitz if and only if e A is Schur.Hence,it follows from Proposition 2.2that the Lyapunov function V (x )=x T P x can be used to establish the stability of both A and e A .In this paper,we consider the nonlinear periodic dynamical system˙x (t )=f (t,x (t )),x (0)=x 0,t ∈I x 0,(2)where x (t )∈D ⊆R n ,t ∈I x 0,is the system state vector,D is an open set,f :[0,∞)×D →R n satisfies f (t,x )=f (t +T,x ),x ∈D ,t ≥0,for some T >0,and I x 0=[0,τx 0),0<τx 0≤∞,is the maximal interval of existence for the solution x (·)of (2).A function x :I x 0→D is said to be a solution to (2)on the interval I x 0⊆[0,∞)with initial condition x (0)=x 0if x (t )satisfies (2)for all t ∈I x 0.It is assumed that f (·,·)is such that the solution to (2)is unique for every initial condition in D and jointly continuous in t and x 0.A sufficient condition ensuring this is Lipschitz continuity of f (t,·):D →R n for all t ∈[0,t 1]and continuity of f (·,x ):[0,t 1]→R n for all x ∈D .Here,we assume that all solutions to (2)are bounded over I x 0,and hence,by the Peano-Cauchy theorem can be extended to infinity.Next,we introduce the notions of periodic solutions and periodic orbits for (2).For the next definition,we denote the solution x (·)to (2)with initial conditon x 0∈D by s (t,x 0).11Notethat since (2)is a time-varying dynamical system it is typical todenote its solution as ˆs (t,t 0,x 0)to indicate the dependence on both the initial time t 0and the initial state x 0.In this paper,we assume that t 0=0and define s (t,x 0)△=ˆs (t,0,x 0).Definition 2.3:A solution s (t,x 0)of (2)is periodic ifthere exists a finite time T >0such that s (t +T,x 0)=s (t,x 0)for all t ≥0.A set O ⊂D is a periodic orbit of (2)if O ={x ∈D :x =s (t,x 0),0≤t ≤T }for some periodic solution s (t,x 0)of (2).Next,we introduce the notions of Lyapunov and asymp-totic stability of a periodic orbit of the nonlinear dynamical system (2).For this definition,dist (p,M )denotes the small-est distance from a point p to any point in the set M ,that is,dist (p,M ) inf x ∈M p −x .Definition 2.4:A periodic orbit O of (2)is Lyapunov stable if,for all ε>0,there exists δ=δ(ε)>0such that if dist (x 0,O )<δ,then dist (s (t,x 0),O )<ε,t ≥0.A periodic orbit O is asymptotically stable if O is Lyapunov stable and there exists ε>0such that if dist (x 0,O )<ε,then dist (s (t,x 0),O )→0as t →∞.Finally,in this paper,we develop a multi-compartment lung model based on a directed tree architecture.The fol-lowing definitions are necessary for the main results of this paper.Definition 2.5([12]):A weighted directed graph G is a triple (V ,E ,W ),where V ={v 1,v 2,...,v N }is the set of vertices ,E ={e 1,e 2,...,e M }⊆V ×V is the set of edges ,and W ∈R N ×N is the weighted adjacency matrix .Every edge e l ∈E corresponds to an ordered pair of vertices (v i ,v j )∈V ×V ,where v i and v j are the initial and terminal vertices of the edge e l .In this case,e l is incident into v j and incident out of v i .The adjacency matrix W is such that W (i,j )>0,i,j =1,...,N,if (v i ,v j )∈E ,and W (i,j )=0otherwise.The in-degree d i (v i )of v i is the number of edges incident into v i and the out-degree d o (v j )of v j is the number of edges incident out of v j .A directed path from v i 1to v i k is a set of distinct vertices {v i 1,v i 2,...,v i k }such that (v i j ,v i j +1)∈E ,j =1,...,k −1.A vertex v i is a root of G if,for every v j =v i ,there exist directed paths from v i to v j .G is connected if,for every pair of v i ,v j ∈V ,there exists v k ∈V such that there are directed paths from v k to v i and v k to v j .A vertex v i ∈V is a leaf of G if d o (v i )=0.Definition 2.6([12]):A weighted directed graph G is a weighted directed tree if G is connected and there exists a vertex v i ∈V such that d i (v i )=0and d i (v j )=1,v j ∈V \{v i }.Remark 2.2:Note that if G is a weighted directed tree,then there exists exactly one root v i ∈V and exactly one directed path from v i to v j for all v j ∈V \{v i }.See [12]for details.III.C OMPARTMENTAL M ODELING OF L UNG D YNAMICS :D ICHOTOMY A RCHITECTURE In this section,we develop a general mathematical model for the dynamic behavior of a multi-compartment respiratory system in response to an arbitrary applied inspiratory pres-sure.Here,we assume that the bronchial tree has a dichotomy architecture [13],that is,in every generation each airway unit branches in two airway units of the subsequent generation.First,however,we start by considering a single-compartment lung model as shown in Figure 1.In this model,the lungs are represented as a single lung unit with compliance c connected to a pressure source by an airway unit with resistance (to air flow)of R .At time t =0,an arbitrary pressure p in (t )is applied to the opening of the parent airway,where p in (t )is determined by the mechanical ventilator.A typical choice for p in (t )is p in (t )=αt +β,where αand βare positive constants.This pressure is applied to the airway opening over'cx (t )=p in (t ),x (0)=x in 0,0≤t ≤T in ,(3)where x (t )∈R ,t ≥0,is the lung volume,R in ∈R is theresistance to air flow during the inspiration period,x in 0∈R is the lung volume at the start of the inspiration and serves as the system initial condition.We assume that expiration is passive (due to elastic stretch of lung unit).During the expiration process,the state equation is given by R ex ˙x (t )+1c ix i (t )+n −1 j =0R inj,k jk j 2n −jl =(k j −1)2n −j +1˙x l (t )=p in (t ),x i (0)=x in i 0,0≤t ≤T in ,i =1,2,...,2n ,(5)'&%$'&%$c 1c 2c 3c 4x 1x 2x 3x 4R in2,1R in 2,2R in2,3R in2,4R in1,1R in1,2R in0,1p appFig.2.Four-compartment lung modelwhere k j =⌊k j +1−1c i x i (t )due to the compliancein i -th compartment will be equal to the difference between the external pressure applied and the resistance to air flow at every airway in the path leading from the pressure source to the i -th compartment.In particular,for i =3(see Figure 2),1c 3x 3(t )=p in (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 is p ex (t ).Following an identical procedure as in the inspiration case,we obtain the state equation for expiration asR ex n,i ˙x i (t )+n −1 j =0R exj,k jk j 2n −jl =(k j −1)2n −j +1˙x l (t )+1IV.S TATE S PACE M ULTI-C OMPARTMENT L UNG M ODEL In this section,we rewrite the state equations(5)and(7) for inspiration and expiration,respectively,as a switched dynamical system.To describe the dynamics of the multi-compartment lung model in terms of a state space model, define the state vector x [x1,x2,···,x2n]T,where x i denotes the lung volume of the i-th compartment.Now,the state equation(5)for inspiration can be rewritten asR in˙x(t)+Cx(t)=p in(t)e,x(0)=x in0,0≤t≤T in,(8) where C diag[1c2n]andR innj=02j k=1R in j,k Z j,k Z T j,k,(9)where Z j,k∈R2n is such that the l-th element of Z j,k is1 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(7)for expiration can be rewritten asR ex˙x(t)+Cx(t)=p ex(t)e,x(T in)=x ex0,T in≤t≤T,(10) where T T ex+T in andR exnj=02j k=1R ex j,k Z j,k Z T j,k.(11)Note that if R in and R ex are invertible,then(8)and(10) can be equivalently written as˙x(t)=A in x(t)+B in p in(t),x(0)=x in0,0≤t≤T in,(12)˙x(t)=A ex x(t)+B ex p ex(t),x(T in)=x ex0,T in≤t≤T ex+T in,(13) where A in△=−R−1in C,B in△=R−1in e,A ex△=−R−1ex C,and B ex R−1ex e.The following proposition states and proves several impor-tant properties of R in,R ex,A in,and A ex that are essential for the main results of this paper.Proposition4.1:Consider the dynamical system(8)and(10).Then the following statements hold:i)R in>0and R ex>0.ii)A T in C+CA in<0.iii)A T ex C+CA ex<0.iv)R in and R ex are strictly ultrametric.v)A in and A ex are compartmental and Hurwitz,andB in≥≥0and B ex≥≥0.Proof:The proof is a direct consequence of Proposition 2.1and Lemma2.2.R n+,the solution x(t),t≥0,to(14)satisfies x(t)≥≥0,t≥0.Theorem4.1:Consider the switched dynamical system (14)where x in0≥≥0.Then x(t)≥≥0,t≥0,where x(t) denotes the solution to(14).Proof:The proof is a direct consequence of Proposition 4.1and Lemma2.1.whereΓei△=ΓexΓin.Similarly,it can be shown thatx ex1=Γie x ex0+Γinδ+θ,(26) whereΓie△=ΓinΓex.More generally,x in m+1=Γei x in m+Γexθ+δ,m=0,1,...,(27) x ex m+1=Γie x ex m+Γinδ+θ,m=0,1,. (28)The following proposition states and proves two key properties forΓei andΓie which are useful in characterizing a periodic orbit for the switched dynamical system G. Proposition5.1:The following statements hold:i)ΓT ex CΓex<C andΓT in CΓin<C.ii)ΓT ei CΓei<C andΓT ie CΓie<C.Proof:The proof is a direct consequence of Proposition 4.1.R n+,the following statements hold:i)lim m→∞x in m=ˆx in and lim m→∞x ex m=ˆx ex.ii)For every t∈[0,T in],limm→∞x(t+mT)=e A in tˆx in+ t0e A in(t−τ)B in p in(τ)dτ, and for every t∈[T in,T],limm→∞x(t+mT+T in)=e A ex tˆx ex+ t0e A ex(t−τ)B ex p ex(τ+T in)dτ.Proof:The proof is a direct consequence of Proposition 5.1.R n+:x=s(t,ˆx),where s(t,ˆx)is the solution to(14)}.(29)With x in0=ˆx note that x in m=ˆx,m=1,2,...,or, equivalently,x(mT)=ˆx,m=1,2,...,which implies that Oˆx is a periodic orbit of(14).The following theorem presents the main result of this paper.Theorem5.1:Consider the switched dynamical system G given by(14).Let x(t)and y(t),t≥0,denote the solutions to(14)with initial conditions x(0)∈Remark5.2:Note that Theorem 5.1is valid for ar-bitrary nonnegative functions(possibly discontinuous) p in(t)and p ex(t)as long as T in0e−A in t B in p in(t)d t and T T in e−A ex t B ex p ex(t)d t arefinite.In the case,where p in(t)=αt+βand p ex(t)=γfor some positive constantsα,β,and γ,θandδare given byθ=A−2in[(αI+βA in)(e A in T in−I)−αA in T in]B in,δ=γA−1ex(e A ex T ex−I)B ex.Remark5.3:Note that the error dynamics e(t),t≥0, given by(30)is a switched dynamical system where each of the switched systems is a linear dynamical system,and V(e)=e T Ce is a common Lyapunov function for both linear systems.VI.A G ENERAL T REE S TRUCTURE M ODELIn this section,we extend the model presented in Sections III–V to the case where the bronchial tree has a general tree architecture[14–16].The general tree structure includes the regular and irregular dichotomy[13].Specifically,let the bronchial tree be represented by a weighted directed tree G=(V,E,R),where each vertex corresponds to a branching point of an airway unit or the terminal compartment(alve-olus)of the lung.In this case,the trachea corresponds to the root v1of the tree and all the alveoli correspond to the leaves of the tree.Every edge,(v l,v m)∈E corresponds to an airway unit and R(l,m),the weight of the edge,corresponds to the resistance of the airway unit;we use R(l,m)=R in l,m and R(l,m)=R ex l,m for resistance during inspiration and expiration,respectively.Let L {v i∈V:v i is a leaf of G}and let the number of leaves of G(or,equivalently,compartments of the lung)be n so that L={v i1,v i2,...,v in},where i k∈{1,2,...,N},k=1,2,...,n,and N is the number of vertices of the graph.To develop the dynamical model for the inspiration process,let c k,k=1,2,...,n,denote the compliance of each compartment,and let x k,k=1,2,...,n, denote the lung volume in the k-th compartment so that the state equations for inspiration are given by1Next,let x [x 1,...,x n ]T so that (31)can be written as R in ˙x (t )+Cx (t )=p in (t )e ,x (0)=x in 0,0≤t ≤T in ,where C diag[1c n ]andR in =(v l ,v m )∈ER in l,m Z l,m Z Tl,m ,(34)where Z l,m ∈R n is such that the k -th eletment of Z l,m is1if v i k ∈L l,m and 0otherwise.An identical procedure yields the state equations for expiration given by R ex ˙x (t )+Cx (t )=p ex (t )e ,x (T in )=x ex 0,T in ≤t ≤T,(35)whereR ex =(v l ,v m )∈ER ex l,m Z l,m Z Tl,m.(36)Note that it can be easily shown that R in >0and R ex >0and it follows from (34),(36),and Proposition 2.1that R in and R ex are strictly ultrametric.Hence,for a general tree structure model all of the results of Sections IV and V are valid with R in and R ex given by (34)and (36),respectively.v 1p appR in1,2v 2R in 2,3v 3R in 3,5v 5x 1,c 1R in 3,6v 6x 2,c 2R in3,7v 7x 3,c 3R in 2,4v 4R in4,8v 8x 4,c 4R in4,9v 9x 5,c 5Fig.3.Five-compartment tree structure modelTo illustrate the general tree structure lung model,consider the five-compartment model shown in Figure 3.Here,the bronchial tree is represented by a weighted directed tree G =(V ,E ,R )consisting of nine nodes V ={v 1,v 2,...,v 9}and eight edges E ={(v 1,v 2),(v 2,v 3),(v 2,v 4),(v 3,v 5),(v 3,v 6),(v 3,v 7),(v 4,v 8),(v 4,v 9)}.In this case,the set of leaves L ={v 5,v 6,...,v 9}corresponds to the five compartments of the lung.Let v i k =v k +4,k =1,...,5.Now,the pressure 1c 3x 3(t )=p in (t )−R in1,2[˙x 1(t )+˙x 2(t )+˙x 3(t )+˙x 4(t )+˙x 5(t )]−R in 2,3[˙x 1(t )+˙x 2(t )+˙x 3(t )]−R in 3,7˙x3(t ),or,equivalently,1。