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

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.

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, ℝ^𝑛 denotes the set of 𝑛×1 real column vectors, and ℝ^𝑛×𝑚, denotes the set of 𝑛×𝑚 real matrices. For 𝑥∈ℝ^𝑛 we write 𝑥≥≥0 (resp., 𝑥≫0) to indicate that every component of 𝑥 is nonnegative (resp., positive). In this case, we say that 𝑥 is nonnegative or positive, respectively. Likewise, 𝐴∈ℝ^𝑛×𝑚 is nonnegative or positive if every entry of 𝐴 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, ℝ^𝑛₊ and ℝ^𝑛₊ denote the nonnegative and positive orthants of ℝ^𝑛, that is, if 𝑥∈ℝ^𝑛, then 𝑥∈ℝ^𝑛₊ and 𝑥∈ℝ^𝑛₊ are equivalent, respectively, to 𝑥≥≥0 and 𝑥≫0. Finally, we write (⋅)^𝑇 to denote transpose, (⋅)′ to denote Frèchet derivative, and 𝛿𝑥 to denote the first variation of the function 𝑥.

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 Figure 1. In this model, the lungs are represented as a single lung unit with nonlinear compliance 𝑐(𝑥) connected to a pressure source by an airway unit with resistance (to airflow) of 𝑅. At time 𝑡=0, a driving pressure 𝑝_in(𝑡) is applied to the opening of the parent airway, where 𝑝_in(𝑡) is generated by the respiratory muscles or a mechanical ventilator. This pressure is applied over the time interval 0≤𝑡≤𝑇_in, which is the inspiratory part of the breathing cycle. At time 𝑡=𝑇_in, the applied airway pressure is released and expiration takes place passively, that is, the external pressure is only the atmospheric pressure 𝑝_ex(𝑡) during the time interval 𝑇_in≤𝑡≤𝑇_in+𝑇_ex, where 𝑇_ex is the duration of expiration.

Figure 1 

Single-compartment lung model.

The state equation for inspiration (inflation of lung) is given by𝑅_in(𝑡)+1𝑥̇𝑐_in(𝑥)𝑥(𝑡)=𝑝_in(𝑡),𝑥(0)=𝑥ⁱⁿ₀,0≤𝑡≤𝑇_in, (1) where 𝑥(𝑡)∈ℝ, 𝑡≥0, is the lung volume, 𝑅_in∈ℝ is the resistance to airflow during the inspiration period, 𝑐_in∶ℝ→ℝ₊ is a nonlinear function defining the lung compliance at inspiration, and 𝑥ⁱⁿ₀∈ℝ₊ 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 𝑝_in(𝑡),0≤𝑡≤𝑇_in, 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𝑅_ex(𝑡)+1𝑥̇𝑐_ex(𝑥)𝑥(𝑡)=𝑝_ex(𝑡),𝑥𝑇_in=𝑥^ex₀,𝑇_in≤𝑡≤𝑇_in+𝑇_ex, (2) where 𝑥(𝑡)∈ℝ, 𝑡≥0, is the lung volume, 𝑅_ex∈ℝ is the resistance to air flow during the expiration period, 𝑐_ex∶ℝ→ℝ₊ is a nonlinear function defining the lung compliance at expiration, and 𝑥^ex₀∈ℝ₊ is the lung volume at the start of expiration.

Next, we develop the state equations for inspiration and expiration for a 2^𝑛-compartment model, where 𝑛≥0. In this model, the lungs are represented as 2^𝑛 lung units which are connected to the pressure source by 𝑛 generations of airway units, where each airway is divided into two airways of the subsequent generation leading to 2^𝑛 compartments (see Figure 2 for a four-compartment model).

Figure 2 

Four-compartment lung model.

Let 𝑥_𝑖,𝑖=1,2,…,2^𝑛, denote the lung volume in the 𝑖th compartment, let 𝑐ⁱⁿ_𝑖(𝑥_𝑖)(resp.,𝑐^ex_𝑖(𝑥_𝑖)), 𝑖=1,2,…,2^𝑛, denote the compliance at inspiration (resp., expiration) of each compartment as a nonlinear function of the volume of 𝑖th compartment, and let 𝑅ⁱⁿ_𝑗,𝑖 (resp., 𝑅^ex_𝑗,𝑖), 𝑖=1,2,…,2^𝑗,𝑗=0,…,𝑛, denote the resistance (to air flow) of the 𝑖th airway in the 𝑗th generation during the inspiration (resp., expiration) period with 𝑅ⁱⁿ₀₁ (resp., 𝑅^ex₀₁) denoting the inspiration (resp., expiration) of the parent (i.e., 0th generation) airway.

As in the single-compartment model, we assume that a pressure of 𝑝_in(𝑡),𝑡≥0, is generated (by the inspiratory muscles) or applied (by a mechanical ventilator) during inspiration. Now, the state equations for inspiration are given by𝑅ⁱⁿ_𝑛,𝑖𝑥̇_𝑖(𝑡)+1𝑐ⁱⁿ_𝑖𝑥_𝑖(𝑡)𝑥_𝑖(𝑡)+^𝑛−1_𝑗=0𝑅ⁱⁿ_{𝑗,𝑘𝑗}^{𝑘𝑗2𝑛−𝑗}_{𝑙=(𝑘𝑗−1)2𝑛−𝑗+1}𝑥̇_𝑙(𝑡)=𝑝_in(𝑡),𝑥_𝑖(0)=𝑥ⁱⁿ_𝑖0,0≤𝑡≤𝑇_in,𝑖=1,2,…,2^𝑛, (3) where 𝑐ⁱⁿ_𝑖(𝑥_𝑖),𝑖=1,2,…,2^𝑛, are nonlinear functions of 𝑥_𝑖, 𝑖=1,2,…,2^𝑛, given by [22]𝑐ⁱⁿ_𝑖𝑥_𝑖≜𝑎ⁱⁿ_𝑖1+𝑏ⁱⁿ_𝑖1𝑥_𝑖,if0≤𝑥_𝑖≤𝑥ⁱⁿ_𝑖1,𝑎ⁱⁿ_𝑖2,if𝑥ⁱⁿ_𝑖1≤𝑥_𝑖≤𝑥ⁱⁿ_𝑖2,𝑎ⁱⁿ_𝑖3+𝑏ⁱⁿ_𝑖3𝑥_𝑖,if𝑥ⁱⁿ_𝑖2≤𝑥_𝑖≤𝑉_𝑇𝑖,𝑖=1,…,2^𝑛, (4) where 𝑎ⁱⁿ_𝑖𝑗, 𝑗=1,2,3, and 𝑏ⁱⁿ_𝑖𝑗, 𝑗=1,3, are model parameters with 𝑏ⁱⁿ_𝑖1>0 and 𝑏ⁱⁿ_𝑖3<0, 𝑥ⁱⁿ_𝑖𝑗, 𝑗=1,2, are volume ranges wherein the compliance is constant, 𝑉_𝑇𝑖 denotes tidal volume, and𝑘_𝑗=𝑘_𝑗+1−12+1,𝑗=0,…,𝑛−1,𝑘_𝑛=𝑖, (5) where ⌊𝑞⌋ denotes the floor function which gives the largest integer less than or equal to the positive number 𝑞. 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 a 2^𝑛-compartment model, consider the four-compartment model shown in Figure 2 corresponding to a two-generation lung model. Let 𝑥_𝑖,𝑖=1,2,3,4, denote the compartmental volumes. Now, the pressure (1/𝑐ⁱⁿ_𝑖(𝑥_𝑖(𝑡)))𝑥_𝑖(𝑡) due to the compliance in 𝑖th 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 𝑖th compartment. In particular, for 𝑖=3 (see Figure 2),1𝑐ⁱⁿ₃𝑥₃(𝑡)𝑥₃(𝑡)=𝑝_in(𝑡)−𝑅ⁱⁿ_0,1𝑥̇₁(𝑡)+𝑥̇₂(𝑡)+𝑥̇₃(𝑡)+𝑥̇₄(𝑡)−𝑅ⁱⁿ_1,2𝑥̇₃(𝑡)+𝑥̇₄(𝑡)−𝑅ⁱⁿ_2,3𝑥̇₃(𝑡), (6)

or, equivalently,𝑅ⁱⁿ_2,3𝑥̇₃(𝑡)+𝑅ⁱⁿ_1,2𝑥̇₃(𝑡)+𝑥̇₄(𝑡)+𝑅ⁱⁿ_0,1𝑥̇₁(𝑡)+𝑥̇₂(𝑡)+𝑥̇₃(𝑡)+𝑥̇₄(𝑡+1)𝑐ⁱⁿ₃𝑥₃(𝑡𝑥)³(𝑡)=𝑝_in(𝑡). (7)

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 𝑝_ex(𝑡),𝑡≥0. Following an identical procedure as in the inspiration case, we obtain the state equation for expiration as𝑅^ex_𝑛,𝑖𝑥̇_𝑖(𝑡)+^𝑛−1_𝑗=0𝑅^ex_{𝑗,𝑘𝑗}^{𝑘𝑗2𝑛−𝑗}_{𝑘𝑙=𝑗2−1𝑛−𝑗+1}𝑥̇_𝑙(𝑡)+1𝑐^ex_𝑖𝑥_𝑖(𝑡)𝑥_𝑖(𝑡)=𝑝_ex(𝑡),𝑥_𝑖𝑇_in=𝑥^ex_𝑖0,𝑇_in≤𝑡≤𝑇_ex+𝑇_in,𝑖=1,2,…,2^𝑛, (8) where𝑐^ex_𝑖𝑥_𝑖≜𝑎^ex_𝑖1+𝑏^ex_𝑖1𝑥_𝑖,𝑎^ex_𝑖2,𝑎^ex_𝑖3+𝑏^ex_𝑖3𝑥_𝑖,𝑖𝑓0≤𝑥_𝑖≤𝑥^ex_𝑖1,𝑖𝑓𝑥^ex_𝑖1≤𝑥_𝑖≤𝑥^ex_𝑖2,𝑖𝑓𝑥^ex_𝑖2≤𝑥_𝑖≤𝑉_𝑇𝑖,𝑖=1,…,2^𝑛, (9)𝑎^ex_𝑖𝑗, 𝑗=1,2,3, and 𝑏^ex_𝑖𝑗, 𝑗=1,3, are model parameters with 𝑏^ex_𝑖1>0 and 𝑏^ex_𝑖3<0, 𝑥^ex_𝑖𝑗, 𝑗=1,2, are volume ranges wherein the compliance is constant, and 𝑘_𝑗 is given by (5).

Next, we provide a smooth (i.e., C^∞) characterization of the nonlinear compliance using sigmoidal functions [23]. Specifically, for inspiration, 𝑐ⁱⁿ_𝑖(𝑥_𝑖) can be approximated as𝑐ⁱⁿ_𝑖𝑥_𝑖≈𝑎ⁱⁿ_𝑖2𝑆^(𝛽)_𝑎,𝑏𝑥_𝑖−𝑆^(𝛽)_𝑐,𝑑𝑥_𝑖,𝑖=1,…,2^𝑛, (10) where 𝑎=−𝑎ⁱⁿ_𝑖1/𝑏ⁱⁿ_𝑖1, 𝑏=(𝑎ⁱⁿ_𝑖2/𝑏ⁱⁿ_𝑖1)+𝑎, 𝑐=−𝑎ⁱⁿ_𝑖3/𝑏ⁱⁿ_𝑖3, 𝑑=𝑎ⁱⁿ_𝑖2/𝑏ⁱⁿ_𝑖3+𝑐, 𝑆^(𝛽)_𝑎,𝑏(𝑥_𝑖)≜1/(𝑏−𝑎)ln(𝜎_𝑏^(−𝛽)(𝑥_𝑖)/𝜎_𝑎^(−𝛽)(𝑥_𝑖))^1/𝛽 with 𝜎_𝑏^(−𝛽)(𝑥_𝑖)≜1/(1+𝑒^{−𝛽(𝑥𝑖−𝑎)}), and 𝛽>0 is an approximation parameter. Figure 4 shows the smoothed approximation of the piecewise linear compliance function 𝑐ⁱⁿ_𝑖(𝑥_𝑖). A similar approximation holds for 𝑐^ex_𝑖(𝑥_𝑖) 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 𝑥≜[𝑥₁,𝑥₂,…,𝑥_2𝑛]^𝑇, where 𝑥_𝑖 denotes the lung volume of the 𝑖th compartment. Now, the state equations (3) for inspiration can be rewritten as𝑅_in𝑥̇(𝑡)+𝐶_in(𝑥(𝑡))𝑥(𝑡)=𝑝_in(𝑡)𝐞,𝑥(0)=𝑥ⁱⁿ₀,0≤𝑡≤𝑇_in, (11) where 𝐞≜[1,…,1]^𝑇 denotes the one vector of order 2^𝑛, 𝐶_in(𝑥) is a diagonal matrix function given by𝐶_in(𝑥)1≜diag𝑐ⁱⁿ₁𝑥₁1,…,𝑐ⁱⁿ_2𝑛𝑥_2𝑛𝑅_in,≜^𝑛2_𝑗=0^𝑗_𝑘=1𝑅ⁱⁿ_𝑗,𝑘𝑍_𝑗,𝑘𝑍^T_𝑗,𝑘, (12)(13) where 𝑍_𝑗,𝑘∈ℝ^2𝑛 is such that the 𝑙th element of 𝑍_𝑗,𝑘 is 1 for all 𝑙 = (𝑘−1)2^𝑛−𝑗+1,(𝑘−1)2^𝑛−𝑗+2,…,𝑘2^𝑛−𝑗, 𝑘=1,…,2^𝑗, 𝑗=0,1,…,𝑛, and zero elsewhere.

Similarly, the state equation (8) for expiration can be rewritten as𝑅_ex𝑥̇(𝑡)+𝐶_ex(𝑥(𝑡))𝑥(𝑡)=𝑝_ex(𝑡)𝐞,𝑥𝑇_in=𝑥^ex₀,𝑇_in≤𝑡≤𝑇_ex+𝑇_in, (14) where𝐶_ex(𝑥)1≜diag𝑐^ex₁𝑥₁1,…,𝑐^ex_2𝑛𝑥_2𝑛,𝑅_ex≜^𝑛2_𝑗=0^𝑗_𝑘=1𝑅^ex_𝑗,𝑘𝑍_𝑗,𝑘𝑍^T_𝑗,𝑘. (15)(16) Finally, it follows from [6, Proposition 4.1] that 𝑅_in and 𝑅_ex are positive definite and, hence, 𝑅_in and 𝑅_ex 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 𝑃𝑉, where 𝑃 is the pressure driving inflation and 𝑉 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 state 𝑥ⁱⁿ₀ 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 solution 𝑥^∗(𝑡),0≤𝑡≤𝑇_in, 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 𝑥^∗(𝑡),0≤𝑡≤𝑇_in, be given by the solution to the minimization problem 𝒥_in(𝑥)=^𝑇in₀𝑥̈^𝑇(𝑡)𝑥̈(𝑡)+𝛼₁𝑝_in(𝑡)𝐞^𝑇𝑥̇(𝑡)d𝑡,𝛼₁≥0, (17) subject to the natural boundary conditions 𝑥(0)=𝑉₀,𝑥̇(0)=0, (18)𝑥𝑇_in=𝑉₀+𝑉_𝑇,𝑥̇𝑇_in=0, (19) where 𝑉₀∈ℝ^2𝑛 is the end expiratory volume and 𝑉_𝑇∈ℝ^2𝑛 is the tidal volume. If 𝛼₁>0, then 𝑥^∗(𝑡),0≤𝑡≤𝑇_in, is given by 𝑥^∗(𝑡)=𝑑₁+𝑑₂𝑡+exp√𝛼₁𝑅^1/2_in𝑡𝑑₃+exp−√𝛼₁𝑅^1/2_in𝑡𝑑₄,𝑡≥0, (20) and if 𝛼₁=0, then 𝑥^∗(𝑡)=𝑑₁+𝑑₂𝑡+𝑑₃𝑡²+𝑑₄𝑡³,𝑡≥0, (21) where 𝑑₁, 𝑑₂, 𝑑₃, and 𝑑₄∈ℝ^2𝑛 are constant vectors determined by the boundary conditions (18) and (19), and 𝑅^1/2_in denotes the (unique) positive-definite square root of 𝑅_in.

Proof. First, note that 𝑝_in(𝑡)𝐞,0≤𝑡≤𝑇_in, in (17) can be eliminated using the state equation (11). Thus, the integrand of the performance criterion (17) can be written as 𝐿_in(𝑥(𝑡),𝑥̇(𝑡),𝑥̈(𝑡))=𝑥̈^𝑇(𝑡)𝑥̈(𝑡)+𝛼₁𝑅_in𝑥̇(𝑡)+𝐶_in(𝑥(𝑡))𝑥(𝑡)^𝑇𝑥̇(𝑡)=𝑥̈^𝑇(𝑡)𝑥̈(𝑡)+𝛼₁𝑥̇^𝑇(𝑡)𝑅_in𝑥̇(𝑡)+𝑥^𝑇(𝑡)𝐶_in(𝑥)(𝑡,𝛼𝑥̇)₁≥0. (22) The first variation of the performance criterion 𝒥_in(𝑥) is given by 𝛿𝒥_in𝑥^∗,𝛿𝑥=^𝑇in₀𝛿𝐿_in𝑥^∗(𝑡),𝑥̇^∗(𝑡),𝑥̈^∗(𝑡)=d𝑡_𝑇in0𝜕𝐿_in(𝑡)+𝜕𝑥𝛿𝑥𝜕𝐿_in(𝑡)+𝜕𝑥̇𝛿𝑥̇𝜕𝐿_in(𝑡=𝜕𝑥̈𝛿𝑥̈)d𝑡𝜕𝐿_in𝜕𝑥̈𝛿𝑥̇+𝜕𝐿_in−d𝜕𝑥̇d𝑡²²𝜕𝐿_in𝜕𝑥̈𝛿𝑥^𝑇in₀+^𝑇in₀𝜕𝐿_in−d𝜕𝑥d𝑡𝜕𝐿_in+d𝜕𝑥̇d𝑡𝜕𝐿_in(𝑡)𝜕𝑥̈𝛿𝑥d𝑡. (23) Using the boundary conditions (18) and (19), it follows that 𝛿𝑥(0)=𝛿𝑥(𝑇_in)=𝛿𝑥̇(0)=𝛿𝑥̇(𝑇_in)=0. Now, since 𝑇_in is fixed, it follows from the fundamental theorem of the calculus of variations that the variation of 𝒥_in(𝑥) must vanish on 𝑥^∗; that is, the extremals optimizing the performance criterion 𝒥_in(𝑥) satisfy the Euler-Lagrange equation 𝜕𝐿_in𝜕𝑥^𝑇−dd𝑡𝜕𝐿_in𝜕𝑥̇^𝑇+d²d𝑡²𝜕𝐿_in𝜕𝑥̈^𝑇=0. (24)
Next, using 𝐶_in(𝑥) given by (12), 𝜕𝐿_in𝜕𝑥^𝑇=𝛼₁𝐶_in(𝑥(𝑡))𝑥̇(𝑡)+𝛼₁𝐶^′_in(𝑥(𝑡))𝑋̇(𝑡)𝑥(𝑡),0≤𝑡≤𝑇_in,𝜕𝐿_in𝜕𝑥̇^𝑇=2𝛼₁𝑅_in𝑥̇(𝑡)+𝛼₁𝐶_in(𝑥(𝑡))𝑥(𝑡),0≤𝑡≤𝑇_in,𝜕𝐿_in𝜕𝑥̈^𝑇=2𝑥̈(𝑡),0≤𝑡≤𝑇_in, (25) where 𝐶^′_in(𝑥(𝑡))≜diag[(𝜕/𝜕𝑥_𝑖)(1/(𝑐ⁱⁿ_𝑖(𝑥_𝑖(𝑡))))] and 𝑋̇(𝑡)≜diag[𝑥̇_𝑖(𝑡)],𝑖=1,…,2^𝑛. Thus, (24) yields the fourth-order differential equation 𝑥⁽⁴⁾(𝑡)−𝛼₁𝑅_in𝑥⁽²⁾(𝑡)=0,0≤𝑡≤𝑇_in, (26) where 𝑥^(𝑛)(𝑡)≜(d^𝑛𝑥(𝑡)/d𝑡^𝑛), with boundary conditions given in (18) and (19). Now, using standard analysis techniques, the solution 𝑥(𝑡),0≤𝑡≤𝑇_in, to (26) satisfies (20) if 𝛼₁>0 and (21) if 𝛼₁=0.

Remark 2. The vectors 𝑑₁, 𝑑₂, 𝑑₃, and d₄ in Theorem 1 can be uniquely determined using the four boundary conditions given by (18) and (19). Specifically, if 𝛼₁=0, it can be shown that 𝑑₁=𝑉₀, 𝑑₂=0, 𝑑₃=(3/𝑇²_in)𝑉_𝑇, and 𝑑₄=−(2/𝑇³_in)𝑉_𝑇. Hence, in this case, 𝑥̇(𝑡) = 𝑑₂+2𝑑₃𝑡+3𝑑₄𝑡² = (6𝑡/(𝑇²_in))𝑉_𝑇(1−(𝑡/𝑇_in))≥≥0, 0≤𝑡≤𝑇_in, which implies that the solution 𝑥^∗(𝑡),0≤𝑡≤𝑇_in, to (26) is increasing during inspiration, and hence, 𝑉_0𝑖≤𝑥^∗_𝑖(𝑡)≤𝑉_0𝑖+𝑉_𝑇𝑖,𝑖=1,…,2^𝑛, where 𝑉_0𝑖,𝑥_𝑖(𝑡) and 𝑉_𝑇𝑖 are the 𝑖th components of 𝑉₀,𝑥(𝑡), and 𝑉_𝑇, respectively. A similar result holds for the case where 𝛼₁>0.
Next, we give the optimal solution 𝑥^∗(𝑡),𝑇_in≤𝑡≤𝑇_in+𝑇_ex, for the expiratory airflow breathing pattern.

Theorem 3. Consider the nonlinear system model for expiration given by (14). Let the optimal solution 𝑥^∗(𝑡),𝑇_in≤𝑡≤𝑇_in+𝑇_ex, be given by the solution to the minimization problem 𝒥_ex(𝑥)=^{𝑇in+𝑇ex}_𝑇in𝑥̈^𝑇(𝑡)𝑥̈(𝑡)+𝛼₂𝑝²_ex(𝑡)𝐞^𝑇𝐞d𝑡,𝛼₂≥0, (27) subject to the natural boundary conditions 𝑥𝑇_in=𝑉₀+𝑉_𝑇,𝑥̇𝑇_in=0, (28)𝑥𝑇_in+𝑇_ex=𝑉₀,𝑥̇𝑇_in+𝑇_ex=0. (29) If 𝛼₂>0, then 𝑥^∗(𝑡),𝑇_in≤𝑡≤𝑇_in+𝑇_ex, satisfies 𝑥⁽⁴⁾(𝑡)−𝛼₂𝑅²_ex𝑥⁽²⁾(𝑡)+𝛼₂𝐶²_ex(𝑥)𝑥(𝑡)+𝛼₂𝐶_ex(𝑥)𝑅_ex𝑥̇(𝑡)−𝑅_ex𝐶_ex(𝑥)(𝑡)(𝑡)𝐶𝑥̇+𝑋^′_ex(𝑥)𝑅_ex𝑥̇(𝑡)−𝑅_ex𝐶^′_ex(𝑥)𝑋(𝑡)(𝑡)(𝑡)𝐶𝑥̇+𝑋^′_ex(𝑥)𝐶_ex(𝑥)𝑥(𝑡)=0, (30) where 𝑋(𝑡)≜diag[𝑥_𝑖(𝑡)] and 𝐶^′_ex(𝑥)≜diag[(𝜕/𝜕𝑥_𝑖)(1/𝑐^ex_𝑖(𝑥_𝑖))],𝑖=1,…,2^𝑛, and if 𝛼₂=0, then 𝑥^∗(𝑡)=𝑑₁+𝑑₂𝑡+𝑑₃𝑡²+𝑑₄𝑡³,𝑡≥0, (31) where 𝑑₁, 𝑑₂, 𝑑₃, and 𝑑₄∈ℝ^2𝑛 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 𝐿_ex(𝑥(𝑡),𝑥̇(𝑡),𝑥̈(𝑡))=𝑥̈^𝑇(𝑡)𝑥̈(𝑡)+𝛼₂𝑝_ex(𝑡)𝐞^𝑇𝑝_ex(𝑡)𝐞=𝑥̈^𝑇(𝑡)𝑥̈(𝑡)+𝛼₂𝑅_ex𝑥̇(𝑡)+𝐶_ex(𝑥(𝑡))𝑥(𝑡)^𝑇×𝑅_ex𝑥̇(𝑡)+𝐶_ex(𝑥(𝑡))𝑥(𝑡)=𝑥̈^𝑇(𝑡)𝑥̈(𝑡)+𝛼₂𝑥̇^𝑇(𝑡)𝑅²_ex𝑥̇(𝑡)+𝑥^𝑇(𝑡)×𝐶²_ex(𝑥(𝑡))𝑥(𝑡)+2𝑥̇^𝑇(𝑡)𝑅_ex𝐶_ex(𝑥(𝑡𝑥(𝑡,𝛼)))₂>0. (32) Thus, the variation of 𝒥_ex(𝑥) on an extremal solution gives 𝛿𝒥_ex𝑥^∗,𝛿𝑥=^{𝑇in+𝑇ex}_𝑇in𝛿𝐿_ex𝑥^∗(𝑡),𝑥̇^∗(𝑡),𝑥̈^∗(𝑡)=d𝑡_{𝑇in+𝑇ex𝑇in}𝜕𝐿_ex(𝑡)+𝜕𝑥𝛿𝑥𝜕𝐿_ex(𝑡)+𝜕𝑥̇𝛿𝑥̇𝜕𝐿_ex(𝑡=𝜕𝑥̈𝛿𝑥̈)d𝑡𝜕𝐿_ex𝜕𝑥̈𝛿𝑥̇+𝜕𝐿_ex−d𝜕𝑥̇d𝑡𝜕𝐿_ex𝜕𝑥̈𝛿𝑥^{𝑇in+𝑇ex}_𝑇in+^𝑇ex₀𝜕𝐿_ex−d𝜕𝑥d𝑡𝜕𝐿_ex+d𝜕𝑥₂̇d𝑡²𝜕𝐿_ex(𝑡)𝜕𝑥̈𝛿𝑥d𝑡=0. (33) Using the boundary conditions (28) and (29), it follows that 𝛿𝑥(𝑇_in)=𝛿𝑥(𝑇_in+𝑇_ex)=𝛿𝑥̇(𝑇_in)=𝛿𝑥̇(𝑇_in+𝑇_ex)=0. Hence, the extremals optimizing the performance criterion 𝒥_ex(𝑥) satisfy the Euler-Lagrange equation 𝜕𝐿_ex𝜕𝑥^𝑇−dd𝑡𝜕𝐿_ex𝜕𝑥̇^𝑇+d²d𝑡²𝜕𝐿_ex𝜕𝑥̈^𝑇=0. (34)
Now, using 𝐶_ex(𝑥) given by (15), 𝜕𝐿_ex𝜕𝑥^𝑇=𝛼₂2𝐶²_ex(𝑥(𝑡))𝑥(𝑡)+2𝐶_ex(𝑥(𝑡))𝑅_ex(𝑡)(𝑡)𝐶𝑥̇+2𝑋^′_ex(𝑥(𝑡))𝑅_ex(𝑡)(𝑡)𝐶𝑥̇+2𝑋^′_ex(𝑥(𝑡))𝐶_ex(𝑥(𝑡𝑥(𝑡,𝑇)))_in≤𝑡≤𝑇_in+𝑇_ex,𝜕𝐿_ex𝜕𝑥̇^𝑇=𝛼₂2𝑅²_ex𝑥̇(𝑡)+2𝑅_ex𝐶_ex(𝑥(𝑡𝑥(𝑡,𝑇)))_in≤𝑡≤𝑇_in+𝑇_ex,𝜕𝐿_ex𝜕𝑥̈^𝑇=2𝑥̈(𝑡),𝑇_in≤𝑡≤𝑇_in+𝑇_ex, (35) which yields (30). Finally, in the case where 𝛼₂=0, (30) collapses to 𝑥⁽⁴⁾(𝑡)=0,𝑇_in≤𝑡≤𝑇_in+𝑇_ex, which satisfies (31).

Remark 4. In the case where 𝛼₂=0, the vectors 𝑑₁, 𝑑₂, 𝑑₃, and 𝑑₄ in Theorem 3 can be uniquely determined using the four boundary conditions (28) and (29). In particular, 𝑑₁ = 𝑉₀+𝑉_𝑇+3𝛽𝑇²_in𝑇_ex𝑉_𝑇+2𝛽𝑇³_in𝑉_𝑇, 𝑑₂ = −𝛽(6𝑇²_in𝑉_𝑇+6𝑇_ex𝑇_in𝑉_𝑇), 𝑑₃ = 𝛽(3𝑇_ex𝑉_𝑇+6𝑇_in𝑉_𝑇), and 𝑑₄=−2𝛽𝑉_𝑇, where 𝛽=1/(3𝑇³_ex+12𝑇²_ex𝑇_in+12𝑇_ex𝑇²_in+4𝑇³_in). Hence, 𝑥̇(𝑡)=𝑑₂+2𝑑₃𝑡+3𝑑₄𝑡²=−6𝛽𝑉_𝑇𝑡(𝑇_in+𝑇_ex−𝑡)−6𝛽𝑉_𝑇𝑡(𝑡−𝑇_in)≤≤0, 𝑇_in≤𝑡≤𝑇_in+𝑇_ex, which implies that the solution 𝑥^∗(𝑡),𝑇_in≤𝑡≤𝑇_in+𝑇_ex, is decreasing during expiration, and hence, 𝑉_0𝑖≤𝑥^∗_𝑖(𝑡)≤𝑉_0𝑖+𝑉_𝑇𝑖, 𝑖=1,…,2^𝑛. The case where 𝛼₂>0 involves the solution to (30), and hence, we have been unable to show that 𝑥^∗(𝑡),𝑇_in≤𝑡≤𝑇_in+𝑇_ex, is decreasing during expiration analytically. However, this has been verified numerically.

Remark 5. If optimal solutions to Theorems 1 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 Theorems 1 and 3 are slightly different. In particular, the inspiratory criterion 𝒥_in(𝑥) 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(𝑥), which corresponds to minimizing the mean standard potential energy in the lung.
It can be seen that the optimal solutions 𝑥^∗(𝑡),𝑡≥0, depend on the variables 𝑇_in, 𝑇_ex, 𝑉₀, and 𝑉_𝑇 through the boundary conditions. Moreover, the nonlinearities in (30) are due to nonlinearities in the lung compliance 𝐶_ex(𝑥), which make analytical solutions to (30) difficult to obtain. It is interesting to note that although the optimal solutions 𝑥^∗(𝑡),𝑇_in≤𝑡≤𝑇_in+𝑇_ex, to (30) during the expiration phase depend on the nonlinear compliance of 𝐶_ex(𝑥), the optimal solutions 𝑥^∗(𝑡),0≤𝑡≤𝑇_in, to (26) during the inspiration phase are independent of the nonlinear system compliance 𝐶_in(𝑥). In the case where 𝑛=0 (i.e., a single-lung-compartment model), 𝑥(𝑡)∈ℝ, 𝑅_ex∈ℝ, and 𝐶_ex(𝑥)=𝐶_ex are constants, (30) reduces to 𝑥⁽⁴⁾(𝑡)−𝛼₂𝑅²_ex𝑥⁽²⁾(𝑡)+𝛼₂𝐶²_ex𝑥(𝑡)=0. (36) This case is extensively discussed in [25] wherein the authors characterize four different solutions to (36) corresponding to 𝛼₂=0, 0<𝛼₂<4𝐶²_ex/𝑅⁴_ex, 𝛼₂=4𝐶²_ex/𝑅⁴_ex, and 𝛼₂>4𝐶²_ex/𝑅⁴_ex.

4. Numerical Determination of Optimal Volume Trajectories

The optimal volume trajectories formulated in Section 3 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 𝑥^∗(𝑡),𝑡≥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 𝑎ⁱⁿ_𝑖1=0.018ℓ/cm H₂O, 𝑏ⁱⁿ_𝑖1=0.0233, 𝑎ⁱⁿ_𝑖2=0.025ℓ/cm H₂O, 𝑎ⁱⁿ_𝑖3=0.2532ℓ/cm H₂O, 𝑏ⁱⁿ_𝑖3=−0.01, 𝑥ⁱⁿ_𝑖1=0.3ℓ, 𝑥ⁱⁿ_𝑖2=0.48ℓ, 𝑎^ex_𝑖1=0.02ℓ/cm H₂O, 𝑏^ex_𝑖1=0.078, 𝑎^ex_𝑖2=0.038ℓ/cm H₂O, 𝑎^ex_𝑖3=0.1025ℓ/cm H₂O, 𝑏^ex_𝑖3=−0.15, 𝑥^ex_𝑖1=0.23ℓ, 𝑥^ex_𝑖2=0.43ℓ, and 𝑖=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 time 𝑇_in=2sec and the expiration time 𝑇_ex=3sec. The two weighting parameters 𝛼₁ and 𝛼₂ differ from person to person. Nominal values for the weighting parameters are 𝛼₁=2.0𝑙/sec³cmH₂O and 𝛼₂=0.1𝑙²/sec⁴cmH₂O, which correspond to spontaneous breathing at rest [25]. Figure 5 shows the optimal air volume 𝐞^𝑇𝑥^∗(𝑡),𝑡≥0, and the optimal airflow rate 𝐞^𝑇𝑥̇^∗(𝑡),𝑡≥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 𝐶_in(𝑥) do not appear in (26). However, 𝑥^∗(𝑡),𝑡≥0, obtained using (30) during expiration involves 𝐶_ex(𝑥), and hence, the airflow curve is asymmetric.

Figure 5 

Volume and airflow rate patterns for the total lung compartments.

Figure 6 shows the sensitivity of the optimal volume and airflow rate patterns to changes in the parameters 𝛼₁ and 𝛼₂. 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 𝛼₁. Furthermore, the asymmetric pattern of the expiratory airflow rate reflects the fact that the minimum value becomes steeper with increasing 𝛼₂. Specifically, if we set the weighting parameter 𝛼₂=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.

Figure 7 shows the driving pressure generated by the respiratory muscles using the optimal air volume 𝐞^𝑇𝑥^∗(𝑡),𝑡≥0. Figure 8 compares the optimal air volume trajectory 𝐞^𝑇𝑥^∗(𝑡),𝑡≥0, with a nonoptimal air volume trajectory 𝐞^𝑇𝑥(𝑡),𝑡≥0, generated by the linear pressure 𝑝_in(𝑡)=20𝑡+5cm H₂O, 𝑡∈[0,𝑇_in], and 𝑝_ex(𝑡)=0cm H₂O, 𝑡∈[𝑇_in,𝑇_in+𝑇_ex], [6]. Note that 𝐞^𝑇𝑥^∗(𝑡),𝑡≥0, switches between the end expiratory level 𝐞^𝑇𝑉₀=0.2𝑙 and the tidal volume 𝐞^𝑇𝑉_𝑇=1.2𝑙. Figure 9 shows the phase portrait of the optimal trajectories 𝑥^∗₁(𝑡) and 𝑥^∗₂(𝑡) and suboptimal trajectories 𝑥₁(𝑡) and 𝑥₂(𝑡). 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 volume 𝐞𝑇𝑥∗(𝑡) and nonoptimal volume 𝐞𝑇𝑥(𝑡) versus time.

Figure 9 

Phase portrait for 𝑥∗1(𝑡) versus 𝑥∗2(𝑡) and 𝑥1(𝑡) versus 𝑥2(𝑡).

Figure 10 

Performance criterion comparison versus time.

Finally, Figure 11 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 𝑖=1,2,3,4, and the values for airway resistances are generated using the results of [30].

Figure 11 

Optimal volume 𝑥∗(𝑡) 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.

Acknowledgments

This research was supported in part by the US Army Medical Research and Material Command under Grant 08108002 and the QNRF under NPRP Grant 4-187-2-060.