ABSTRACT Title of Dissertation: AUTOMATED MEDICATION INFUSION SYSTEM DESIGN Xin Jin, Doctor of Philosophy, 2019 Dissertation directed by: Associate Professor Jin-Oh Hahn, Department of Mechanical Engineering Automated infusion of medications will be increasingly deployed in patient care as a means to deliver high-quality and continuous monitoring and therapy, and also to alleviate the excessive workload imposed on the clinicians. Therefore, a well-designed automated medication infusion system is an attractive alternative to today?s manual treatment requiring caregiver?s interventions. However, it also presents numerous challenges: 1) Significant inter- and intra-subject variability; 2) Complexity of medication infusion model; 3) Complexity of interaction of multiple medications; 4) Difficulty in coordination of medical targets. In this dissertation, a well-designed automated medication infusion was designed to address the various challenges: First, considering the large degree of individual subject variability, an adaptive controller was pursued instead of non-adaptive controllers which might be difficult to offer decent behavior for all subjects. Since the infusion model of single drug is highly nonlinear and complex, a low-order single-input single-output (SISO) model was proposed and a SISO semi-adaptive control approach which only adapt can adapt model parameters having a large impact on the model?s fidelity was designed. Secondly, the complicated interaction of multiple medications makes the adaptive controller for two medications even more difficult to design. So a model for two interacting dose responses was constructed and linearized at one operation point. Then the SISO semi-adaptive controller was extended to a two-input two-output case. However, this controller is only designed at one operating point. Therefore, based on two models associated with two distinct operating regimes, a two-model switching control technique was developed and combined with the semi-adaptive controller. Thirdly, we presented a coordinate mechanism to deal with the medical targets setting problem. In real clinical scenarios, the reference targets are empirically specified by caregivers, which are not always achievable in all subjects. Therefore, our proposed coordinate mechanism can recursively adjusts the reference targets based on the estimated dose-response relationship of subjects. Lastly, we conducted SISO control experiments on some pigs. Based on the results, large transport delay was observed in the medication infusion of one pig. Therefore, we incorporated transport delay in the controller design. AUTOMATED MEDICATION INFUSION SYSTEM DESIGN by Xin Jin Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park, in partial fulfillment Of the requirements for the degree of [Doctor of Philosophy] [2019] Advisory Committee: Associate Professor Hahn, Jin-Oh, Chair Professor Azarm, Shapour Professor Balachandran, Balakumar Associate Professor Chopra, Nikhil Associate Professor Sanner, Robert ? Copyright by [Xin Jin] [2019] Table of Contents Table of Contents ................................................................................................................ ii List of Tables ..................................................................................................................... iv List of Figures ..................................................................................................................... v Chapter 1: Introduction ....................................................................................................... 1 1.1 Medication Infusion ............................................................................................ 1 1.2 Benefits of Automated Medication Infusion....................................................... 2 1.3 Challenges ........................................................................................................... 3 1.4 Dissertation Contributions .................................................................................. 6 1.5 Outline................................................................................................................. 7 Chapter 2: Low-Order Dose-Response Model ................................................................... 9 2.1. Direct Dynamic Dose-Response Model .............................................................. 9 2.2. Model Identification and Evaluation: Remifentanil Test .................................. 10 2.3. Model Identification and Evaluation: Propofol Test ........................................ 16 2.4. Model Identification and Evaluation: Vasopressor Test ................................... 18 2.5. Sensitivity Analysis .......................................................................................... 20 Chapter 3: SISO Semi-Adaptive Control Design ............................................................. 25 3.1. Model-Reference Adaptive Control (MRAC) .................................................. 25 3.2. Composite Adaptive Control (CAC) ................................................................ 27 3.3. In-Silico Implementation and Simulation ......................................................... 29 Chapter 4: Two-Input Two-Output Semi-Adaptive Control Design ................................ 34 4.1. Dose-Response Model for Two Interacting Medications ................................. 34 4.2. Control Design .................................................................................................. 38 4.3. In-Silico Implementation and Simulation ......................................................... 41 Chapter 5: Switching Two-Input Two-Output Semi-Adaptive Control Design ............... 45 5.1. Control-Oriented Dose-Response Models ........................................................ 45 5.2. Semi-Adaptive Switching Control .................................................................... 48 5.3. In-Silico Evaluation .......................................................................................... 57 5.4. Results and Discussion ..................................................................................... 61 Chapter 6: Coordinated Two-Input Two-Output Semi-Adaptive Control Design ........... 69 6.1. Control Architecture ......................................................................................... 69 6.2. Coordinated Control via Recursive Reference Adjustment .............................. 70 6.3. In-Silico Implementation and Simulation ......................................................... 75 Chapter 7: Pilot Animal Experiments and Further Improvements for the Controller Design ............................................................................................................................... 83 7.1. Pilot Animal Experiment .................................................................................. 83 7.2. Control Design Model: Direct Dynamic Dose-Response Model ..................... 85 7.3. Parametric Sensitivity Analysis ........................................................................ 86 7.4. Sensitivity-Based Semi-Adaptive Pole Placement Control Design .................. 91 7.5. Dose-Response Model Parametrization ............................................................ 94 7.6. Semi-Adaptive Pole Placement Control ........................................................... 95 7.7. In-Silico Implementation and Simulation ......................................................... 99 7.8. Results and Discussion ................................................................................... 100 ii Chapter 8: Conclusions and Future Work ....................................................................... 105 iii List of Tables Table 1: System identification results for direct dynamic dose-response models. Table 2: System identification results for PKPD model. Table 3: Performance of adaptive and non-adaptive controllers. Table 4: Root-mean-squared normalized command tracking errors (RMSNE?s) associated with switching and non-switching controllers, and switching times associated with switching controllers (mean (SD)). Table 5: Command tracking performance of semi-adaptive switching and non-switching controllers with respect to operating regime (mean (SD)). Table 6: Root-mean-squared errors (RMSEs) between the responses associated with the original and all the Pad?-approximated dose-response models (M1~M6). iv List of Figures Figure 1: Classical PKPD (a) versus direct dynamic dose-response (b) models. Figure 2: Measured versus model-predicted RR responses of 24 subjects (one- compartment DDDR model). Figure 3: Measured versus model-predicted CO responses of propofol for 5 pigs (one- compartment DDDR model). Figure 4: Measured versus model-predicted HR responses of 4 pigs of vasopressor for 4 pigs (one-compartment DDDR model). Figure 5: (a) Remifentanil infusion dose profile for in-silico sensitivity analysis. (b) The resulting RR response of the nominal direct dynamic dose-response model. Figure 6: The impact of (a) ? and (b) ?50 on the model?s response. Figure 7: True RR response simulated by direct dynamic dose-response model and its prediction by semi-adaptive direct dynamic dose-response model. Figure 8: Closed-loop controlled RR responses and remifentanil infusion rates associated with (a) MRAC, (b) CAC and (c) non-adaptive control (commanded RR set point = 15 bpm). Figure 9: Representative in-silico results associated with MRAC and CAC. Figure 10: Distribution of parameter estimation errors in the steady state. (a) ??. (b) ?50. Figure 11: Dose-response model for two interacting medications, consisting of a low- order mixing model and a response surface model. Figure 12: In-silico simulation testing results associated with the semi-adaptive and non- adaptive controls in the 20 in-silico subjects with achievable reference targets. v Figure 13: Architecture of semi-adaptive switching control approach. Figure 14: Representative in-silico simulation examples. Blue vertical line indicates the instant at which mode switching occurred. (a) An in-silico subject in the propofol- dominant regime at the set point. (b) An in-silico subject in the remifentanil-dominant regime at the set point. M1: propofol. M2: remifentanil. Figure 15: Parameter Estimation RMSNE?s associated with the semi-adaptive switching versus non-switching control. (a) 28 in-silico subjects in the propofol-dominant regime at the set point. (b) 22 in-silico subjects in the remifentanil-dominant regime at the set point. Figure 16: Mean percent improvement in command tracking performance offered by switching control with respect to operating regime. (a) In-silico subjects subject to propofol-dominant regime. (b) In-silico subjects subject to remifentanil-dominant regime. Figure 17: The architecture of coordinated semi-adaptive control. Figure 18: The set of achievable reference targets specified by the constraints on the input magnitudes and deviation from the original reference targets. Figure 19: The evolution of adjusted reference targets and system outputs associated with the coordinated semi-adaptive control in a in-silico subject with unachievable reference targets. (a) Time plots. (b) Phase plot. Figure 20: Distribution of overall reference tracking and steady-state errors as well as total medication use associated with semi-adaptive control in the presence/absence of coordination control. vi Figure 21: Dependence of the behavior of coordinated semi-adaptive control on cost function weights. Figure 22: Experimental results. (a) 3 pigs with small dose-response delay. (b) 1 pig with large dose-response delay. Figure 23: Root-mean-squared errors (RMSEs) between nominal versus perturbed clinical effect responses to perturbation of dose-response model parameters in 30 in-silico subjects in a case study of regulating a cardiovascular variable (cardiac output (CO)) with a sedative (propofol). Figure 24: Sensitivity-based semi-adaptive pole placement control (semi-APPC) for closed-loop control of medication infusion. The semi-APPC is built upon (1) Pad? approximation of transport delay, (2) novel linear model parameterization, (3) recursive model parameter adaptation, and (4) its integration into pole placement control. Figure 25: Root-mean-squared error (RMSE), median absolute percentage error (MDAPE), and wobble metrics computed for semi-APPC, semi-APPC with nominal transport delay, semi-APPC with worst-case transport delay, and population-based PID control associated with 5 min time constant. For each target cardiac output (CO), each bar denotes from the left semi-APPC, semi-APPC with nominal transport delay, semi- APPC with worst-case transport delay, and population-based PID control. Figure 26: Representative examples of the responses of 30 in-silico subjects to target cardiac output (CO) of 2.7 lpm with 5 min time constant associated with semi-APPC, semi-APPC with nominal transport delay, semi-APPC with worst-case transport delay, and population-based PID control. vii Figure 27: Steady-state error variability of 30 in-silico subjects associated with semi- APPC, semi-APPC with nominal transport delay, semi-APPC with worst-case transport delay, and population-based PID control viii Chapter 1: Introduction In this chapter, the mediation medication and the advantages of automated medication infusion is introduced first. Then the challenges of computerized medication infusion system is discussed. Moreover, the approaches to address them is presented. In the end, the outline of the dissertation is presented. 1.1 Medication Infusion Medication infusion means that a drug is administered intravenously, but the term also may refer to situations where drugs are provided through other non-oral routes, such as intramuscular injections and epidural routes into the membranes surrounding the spinal cord. Typically, the medications sent directly into vein using a needle or tube intravenously. With intravenous administration, a thin plastic tube called an intravenous catheter is inserted into your vein. The catheter allows your healthcare provider to give you multiple safe doses of medication without needing to poke you with a needle each time. Medication infusion has a wide range of application including anesthesia, infections, chemotherapy, dehydration, gastrointestinal diseases or disorders which prevent normal functioning of the gastrointestinal system, and more. To achieving high-quality medication infusions, infusion pharmacies need to provide an extensive array of professional services?patient assessment and admission, education and training, care planning and coordination, care management by clinical infusion pharmacists, trouble- shooting and treatment plan oversight, and much more. 1 Take the current anesthesia care for an example. Anesthesiologists administer the doses based on the average patient. They first infuse the initial dose of anesthetics, observe the response, and adjust the dose accordingly. The performance highly depends on the anesthesiologist?s understanding of both the pharmacokinetics (PK) and the pharmacodynamics (PD) of the drugs in use, as well as the possible drugs interaction. The anesthesiologist hence acts as a feedback controller [1]. Therefore, it is natural to think about if automatic medication infusion is able to be used for research in clinical studies. 1.2 Benefits of Automated Medication Infusion Studies investigating the efficacy of closed-loop medication infusion relative to manual medication infusion suggest that the former performs as good as, or sometimes even outperforms, the latter in sticking to the specified target clinical endpoint in anesthesia care [2]?[4], vasopressor therapy [5] and fluid resuscitation [6]?[8]. Computerized medication infusion may also reduce physiological variability during treatment relative to manual therapy [7], [8], prevent overdosing [9], [10], and facilitate medication weaning [11]. Noting that caregivers are extremely overloaded and frequently distracted by multiple tasks [12]?[14], a well-designed automated medication infusion system is an attractive alternative to today?s manual treatment requiring caregiver?s interventions. 2 1.3 Challenges To develop an automated medication infusion system, several challenges are faced as summarized below. 1.1.1. Model Uncertainty: One of the key challenges in automated medication infusion is the signification inter- and intra-subject variability [7], [15], [16]. However, some existing studies only focused on non-adaptive control approaches [17]?[23]. Though suggested to be robust, these controllers have only either been tested in limited subject populations, and are probably susceptible to uncertainties due to a large inter-individual variability were they to be tested over a wide range of diverse subjects, or designed for the worst case situation, which may make them sluggish. Therefore, a model based adaptive controller should be exploited to deal with the large uncertainty. 1.1.2. Complexity and Nonlinearity of Dose Response Model: Another major source of difficulty results from the complexity of the subjects? pharmacokinetics-pharmacodynamics (PKPD) responses to the drug delivery. Even if for a SISO classical PKPD model, it involves a number of individual-specific parameters that cannot be easily tuned with the dose profiles anticipated in real clinical settings (i.e., dose profiles are typically too simple to faithfully determine all the individual-specific parameters in classical PKPD models). Further, the PKPD model usually involves nonlinear PD models (such as the Hill equation) that cannot be adapted easily using standard linear parameter estimation techniques. Due to these challenges, some of existing adaptive control work are only based on black-box models [24], [25]. But no 3 consensus has been established as to the best black-box model or even as to if a particular black-box model structure is universally applicable to wide-ranging subject populations. Some others are based on the PKPD models. To design an adaptive control, they have employed techniques such as linearizing the Hill equation [26], fixing either PK or PD model while adapting the other [27], [28], using nonlinear filtering techniques (such as the extended Kalman filter) to adapt the nonlinear Hill equation while constraining the PKPD model to reduce the number of parameters to be tuned [29]. However, all these methods have some disadvantages. They may only works at a small region of operation or the parameters they fix actually have large impact on the model?s response. So a better approach which can avoid these disadvantages needs to be proposed. 1.1.3. Complexity of Interaction of Multiple Medications Generally, multiple medications must be infused during care of critically ill subjects to achieve a multitude of treatment goals. However, the interaction of multiple medications is also nonlinear [30]?[34] which makes the adaptive controller even more difficult to design. So there is only limited volume of work reported on the closed-loop control for infusion of multiple medications to track multiple reference targets [35]?[39]. The limited work are usually non-model based [35], [36] or black-box model based [37], [38]. However, it is desirable to design controllers based on physical models which may have better performance. Although [39] utilized a linear approximation of PKPD model around the maintenance values and then applied model based predictive controller for it, it may not have very good performance over a large region of operation. Therefore, a physical model based controller which can work over a wide region should be proposed despite of complexity of the MIMO medication infusion models. 4 1.1.4. Difficulty in Coordination of Medical Targets In addition to closed-loop control design, the coordination of reference targets also presents a challenge. In real clinical scenarios, reference targets are empirically specified by caregivers, e.g., based on population norms and caregivers? experience. These ad-hoc reference targets are not always achievable in all subjects due to the inter-individual variability in dose-response relationships and the bounds on medication dose to ensure subject safety. In fact, inappropriate coordination of reference targets that cannot be achieved in a subject may potentially harm the subject via over-/under-dosing. A critical challenge is that it is not possible to specify achievable reference targets for a subject before the treatment, since the dose-response relationship of the subject is typically not known a priori. Therefore, targets must be recursively adjusted by estimating the subject?s dose-response relationship on-line while respecting the caregivers? therapeutic intent. Existing well-known techniques for the adjustment of reference targets are reference governors and their variants, add-on schemes used to avoid the violation of state and input constraints in a closed-loop control system by adjusting the reference targets during the transients [40]?[42]. Yet, these techniques are not appropriate to address the challenge at hand for at least two reasons. Most importantly, the primary goal of a reference governor is to keep the adjusted reference targets as close as possible to the original ones, and ideally, ultimately bring them back to the original ones. In contrast, our goal is to adjust inherently unachievable reference targets to achievable ones so that the system (the subject in our case) can converge to the adjusted reference targets, while accounting for additional safety-critical considerations such as minimizing the total medication doses. Furthermore, a reference governor can only reduce the reference targets, while our goal is to adjust (i.e., increase or decrease) reference targets according 5 to the sensitivity of each subject to medications and their synergy. Therefore, novel approaches for coordinated adjustment of reference targets may be beneficial in providing closed-loop medication infusion control systems with resilience to inappropriate reference targets specified by caregivers, ultimately improving the safety of subjects receiving medication treatments. 1.4 Dissertation Contributions In an effort to cope with these challenges, we first exploited a low-order dose-response model and presented a SISO semi-adaptive control approach to the model. The rationale underlying this approach was to design a controller that can adapt model parameters having a large impact on the model?s fidelity while fixing the remaining parameters at nominal values to linear parameterize the model because standard adaptive techniques are based on linear parameterized models. Moreover, the SISO approach was extended the two-input two-output case. Because the two-input two-output model is nonlinear, we then linearized the multivariable dose- response model at one operating point to enable the controller design. However, given that the coupling effect between multiple medications frequently exhibits complex behaviors due to the regime-dependent inter-medication synergy, the controller based on one single model maybe not suffice for a trustworthy design. Therefore, we linearized the multivariable model at two distinct operating points and design an adaptive controller for each model. Two locally linearized modes represent the operating regime in which one 6 medication is used as primary therapy for controlling the clinical responses. Then a two- model switching control technique was developed to choose the suitable controller. In addition, because the reference targets have to be adjusted based on the responses of the subjects, a coordinate mechanism was then proposed to recursively adjusts the reference targets based on the estimated dose-response relationship of a subject to ensure that they can be achieved by the subject and minimize the medication use at the same time. In the end, we implemented our SISO controller on pigs to validate the proposed approach. The semi-adaptive controller performed well in 3 pigs while it performed marginally in 1 pig. Our retrospective analysis indicated that dynamic dose-response delay, which is not modeled explicitly in the controller design process, may play an important role in determining the performance of the controller in each pig. Therefore, we accommodated the dose-response delay in the SISO controller design. 1.5 Outline This proposal is organized as follows. Chapter 2 describes a low-order model named the direct dynamic dose-response model and its system identification and sensitivity analysis results. Chapter 3 elaborates on the design of SISO semi-adaptive controllers. Chapter 4 extends the SISO controller to a two-input two-output semi-adaptive controller. Chapter 5 combines the two-model switching techniques to MIMO controller in Chapter 4. Chapter 6 describes the coordinated semi-adaptive controller which combines a coordination 7 mechanism with the MIMO controller in Chapter 4. Chapter 7 shows the pig experiment results of SISO controller and the method to improve the SISO controller. Chapter 8 concludes the dissertation with future directions. 8 Chapter 2: Low-Order Dose-Response Model In this chapter, the SISO model for medication infusion is introduced. Since classical adaptive controller can be usually applied to linearly parametrized models, so we need to propose a model which is accurate enough to represent the dynamics of medication infusion and suitable for adaptive control design at the same time. To achieve this goal, a low-order dose-response model was first introduced and validated using three sets of data for different drugs. However, this model is still nonlinear which does not lend itself to standard adaptive control design techniques. Therefore, sensitivity analysis was performed to examine the impact of each parameter. Based on it, we can fix the parameter which has the least impact on model?s response to linearly parametrized the model. 2.1. Direct Dynamic Dose-Response Model A direct dynamic dose-response model proposed by Hahn et al. [43] was used to reproduce the effect of medication infusion. In the classical PKPD model (Fig. 1a), intermediate state variables (plasma and effect site concentrations of medication, ?? and ??) as well as input (medication dose, ?) and output (? normalized by its baseline, ?0) are directly related to physical quantities. As such, PK and PD can be modeled using the measurements of ?, ??, and ?. However, this is not possible in routine procedures where ?? cannot be measured in real time. The direct dynamic dose-response model eliminates the reliance of the PKPD model on ?? by deriving a direct relationship between ? and ? (Fig. 1b). Compared with the PKPD model, ?????(?) is replaced by a unity-gain transfer function ?????(?), output of which represents a hypothetical medication dose at the site 9 of action, ??. In addition, the Hill equation model in the direct dynamic dose-response model relates ?? to ? in contrast to the PKPD model in which the same Hill equation relates ?? to ? . In this regard, the parameter ?50 can be interpreted as the medication infusion rate associated with 50% depression of ? from its baseline. In sum, the direct dynamic dose-response model can be tuned to each individual solely based on the input- output data, even in the absence of ?? measurements. Note that although the PKPD model may also be tuned in the same way using ? and ?, the model thus derived may not always be faithful since the dose profiles anticipated in real clinical settings are not rich enough to tune a number of parameters involved in the PKPD model. (a) (b) Figure 1: Classical PKPD (a) versus direct dynamic dose-response (b) models. 2.2. Model Identification and Evaluation: Remifentanil Test Multiple medications including remifentanil, propofol and vasopressor were utilized to validate the direct dynamic dose-response model. At first, the remifentanil dose and respiratory rate data collected from 24 pediatric patients undergoing dental restoration [44] were used to derive and compare classical PKPD and direct dynamic dose-response models. The mixed effects modeling analysis was conducted using the NLMEFIT routine available in the MATLAB Statistics Toolbox [45] for system identification as follows. 10 1) Direct Dynamic Dose-Response Model: In this model, the hypothetical remifentanil infusion rate ?? at the site of action is assumed to be related to its intravenous counterpart by ?? ??(?) = ?????(?)?(?) = ?(?) (1) ? + ?? where ?? is the equilibration rate constant and ? is the Laplace operator. The relation between ?? and RR is modeled by the Hill equation model: ??? (?) ??(?) = ??0 [1 ? ] (2) ??50 + ? ? ? (?) where RR0 is the baseline RR in the absence of drug effect and ? is a cooperativity constant. Each of the unknowns is assumed to be composed of a nominal value and a random effect due to the inter-individual variability: ?? = ????? + ? ????1, ?50 = ?50 + ?2, ? = ??? + ?3, ?? = ???????0 0 + ?4 (3) where ? ? represents the nominal value while ?? , ? = 1,? ,4 are zero-mean, normally distributed random variables that represent the inter-individual variability. The nominal value and the variance associated with each ?? are identified simultaneously. Noting the fast and ultra-short acting nature of remifentanil, a zero-compartment model (?? = 0) and a one-compartment model are considered as candidate models for ?????(?). 2) PKPD Model: In the PKPD model, ?? is predicted using a population-based PK model of remifentanil [46]. The effect site remifentanil concentration ?? is related to ?? by ??0 ??(?) = ? (?) (4) ? + ? ??0 11 where ??0 is the equilibration rate constant. The nonlinear Hill equation model is used to relate ?? to ?? , where ?50 denotes the effect site remifentanil concentration corresponding to 50% depression of RR from the baseline and ? is a cooperativity constant: ??? (?) ??(?) = ??0 [1 ? ] (5) ??50 + ? ? ? (?) Similarly to the direct dynamic dose-response model, the nominal value and the variance associated with each ?? (? = 1,? ,4 corresponds to ??0, ?50, ?, and RR0, respectively) are identified simultaneously. 3) Prediction Capability Analysis: The predictive accuracy of the direct dynamic dose- response versus PKPD models was compared in terms of the root-mean-squared errors (RMSE) in individual prediction and the Akaike?s Information Criterion (AIC), as computed by the NLMEFIT routine in the MATLAB Statistics Toolbox. Table 1 and Table 2 summarize the system identification results for the direct dynamic dose-response and PKPD models, respectively. Fig. 2 presents measured versus model- predicted (by the direct dynamic dose-response models) RR responses of 24 subjects. The mixed effects modeling analysis resulted in direct dynamic dose-response and PKPD models with satisfactory individual prediction capability. The model parameters identified for the Hill equation model in both the zero-compartment and one- compartment direct dynamic dose-response models were consistent and comparable to each other (Table 1). This suggests that the Hill equation model was reliably identified in both models to reproduce the steady-state RR response to remifentanil infusion. In fact, 12 the consistency in the Hill equation models associated with these direct dynamic dose- response models is not surprising if the distinct roles of the transfer function model ?????(?) and the Hill equation model in Fig. 1 are taken into account: the former represents the dynamic transients related to the distribution of remifentanil whereas the latter stands for the relationship between the remifentanil dose (??) and the RR response in the steady state. Table 1: System identification results for direct dynamic dose-response models. NV: nominal value. s.e.: standard error. ?IIV: inter-individual variability. RMSE: root mean squared error. AIC: Akaike?s Information Criterion. Zero-Compartment Model One-Compartment Model ?? ?50 ? RR0 ?? ?50 ? RR0 [min-1] [mcg/kg/min] [-] [bpm] [min-1] [mcg/kg/min] [-] [bpm] NV (s.e.) - 0.10 (0.02) 2.09 (1.92) 27.7 (1.99) 0.21 (0.03) 0.08 (0.01) 3.12 (0.24) 26.9 (1.69) ?IIV - 0.04 1.78 8.60 0.16 0.03 2.38 7.87 RMSE 2.15 m-1 1.62 m-1 AIC 1.12?105 0.98?105 Table 2: System identification results for PKPD model. Rigby-Jones PK model [46] was used to predict the plasma concentration response to remifentanil infusion. NV: nominal value. s.e.: standard error. ?IIV: inter-individual variability. RMSE: root mean squared error. AIC: Akaike?s Information Criterion. Effect Compartment PD Model ??0 ??0 [min-1] [min-1] NV (s.e.) 3.80 (2.59) NV (s.e.) 3.80 (2.59) NV (s.e.) ?IIV 5.55 ?IIV 5.55 ?IIV RMSE 1.86 m-1 AIC 1.05X105 13 Figure 2: Measured versus model-predicted RR responses of 24 subjects (one- compartment DDDR model). Circles: RR measurements [m-1]. Solid lines: RR predictions [m-1]. Dashed lines: remifentanil infusion rate [X102 mcg/kg/min]. Compared with its zero-compartment counterpart, a large inter-individual variability was observed for ? in the one-compartment direct dynamic dose-response model (Table 1). The increase in the variability associated with ? in the one-compartment model may be attributed to its interaction with ??. Indeed, the effect of ?? and ? on the RR response is counteracting to each other: in the low infusion rate regime (which is the case of this study), an increase in ?? and a decrease in ? both results in an increase in the speed of response. As a result, it is not trivial to faithfully identify ?? and ? simultaneously. This is also supported by a strong negative covariance between the two in the variance- covariance matrix of the random effects (not shown). In contrast to ?, the inter-individual variability associated with ?50 was small for both models. This implies that ?50 has a large 14 impact on the model?s response and thus may be effectively identified. Also, a modest positive covariance was observed between ? and ?50, meaning that their influence on the model?s response is synergistic: an increase in both yields a decrease in the speed of response. In sum, this pattern illustrates how the random effects are incorporated into their respective nominal values in identifying the individualized direct dynamic dose- response models based on the RR response of each subject. For instance, for a subject whose response is not as sensitive to remifentanil as the population average, the random effects are determined so that 1) ? < ??? ?? ?, 2) ? > ???, and 3) ? ????50 > ?50, which fulfills 1) sluggish equilibration between plasma and effect site, 2) extended dead zone in the Hill equation model to suppress RR responses at low ?? values, and 3) high ?? threshold in the Hill equation model to prevent RR response in low ?? values, altogether decreasing the sensitivity of an individual?s RR response to remifentanil. Compared to the PKPD model, the one-compartment direct dynamic dose-response model exhibited an improved RMSE (by 13 %) as well as the Akaike?s Information Criterion (AIC) (Table 1 and Table 2). On the other hand, the zero-compartment direct dynamic dose-response model was inferior to the PKPD model, both in terms of RMSE and AIC. This suggests that the plasma-effect site equilibration dynamics of remifentanil must be considered in the direct dynamic dose-response model regardless of the ultra- short acting nature of remifentanil. It is also noted that the two-compartment direct dynamic dose-response model did not offer a large incremental improvement in the predictive capability of the direct dynamic dose-response model (not shown). In sum, the smaller AIC value of the direct dynamic dose-response model than the PKPD model 15 indicates its superior trade-off between predictive capability and model parsimony over the PKPD model. The deteriorated performance of the PKPD model, despite its structural complexity, is attributable to its use of population-predicted ?? in that the fidelity of the overall PKPD model and the prediction of the PD response are adversely affected when the discrepancy between predicted versus underlying ?? is unacceptably large, as illustrated in a previous study [47]. 2.3. Model Identification and Evaluation: Propofol Test To validate the model can capture the dynamics of other medications, we also tested propofol. The propofol dose and cardiac output (CO) data from 5 pigs were collected under the protocol approved by the Institutional Animal Care and Use Committee (IACUC) at University of Maryland. Pigs were chosen for this study for several reasons. They are the accepted standard model for studies of cardiovascular physiology due to the high degree of similarity between human and pig cardiovascular systems. In addition, they are of adequate size to accept the instrumentation to be used and are hardy enough to tolerate the dose escalations to be performed while under anesthesia. Each pig received propofol at several distinct infusion rates under general anesthesia and mechanical ventilation to elicit a wide range of CO response. The CO was measured per second. During experiments, RR was perturbed on purposely in some time intervals to examine its effect on CO. Because the perturbation of RR also had effect on CO, these data were dropped and not considered in the optimization. So there are some gaps in the 16 CO measurements as shown in Fig. 3. Before analysis, we pre-processed the CO measurements with a 25-point median filter to remove noise and outliners. We performed system identification of the dose-response model by fitting the model to the experimental data associated with each pig by the MATLAB Optimization Toolbox to minimize the prediction error. As shown in Fig. 3, the DDDR model can capture the response in general which implies its superior predictive capability. Figure 3: Measured versus model-predicted CO responses of propofol for 5 pigs (one- compartment DDDR model). 17 2.4. Model Identification and Evaluation: Vasopressor Test In contrast to the above medications, medications such as vasopressors and inotropes exhibit excitatory dose-dependent effects. Similar to the DDDR model, we can reduce classical ???? model to derive a similar low order dose-response model: ?? ??(?) = ?????(?)?(?) = ?(?) (6) ? + ?? where ?(?) is the intravenous medication infusion rate, ??(?) the infusion rate at the (hypothetical) effect site, ?? the time constant associated with the distribution of medication between blood and effect site. The relation between ?? and HR (heart rate) is modeled by the Hill equation model: ??? (?) ??(?) = ??0 + ???? (7) ??50 + ? ? ? (?) where ???? is maximum possible increase of HR, ? a cooperativity constant, and HR0 the baseline heart rate (i.e., in the absence of medication infusion). Different from the model (5), it has an extra unknown parameter ???? and therefore cannot be readily incorporated into control design. In fact, it is extremely difficult to determine the upper limit of response in a subject in real clinical settings due to subject safety and ethical considerations. To address this challenge, my lab mate Junxi developed a new dose- response model by a nonlinear transformation: ?? ??(?) = ?????(?)?(?) = ?(?) (8a) ? + ?? 1 ??50(?) ??(?) = ??0 [1 + ln ] 1 ? ? (8b) 2ln ( ) ?50 + 2?? (?)3 18 To validate and analyze the proposed dose-response model, Junxi utilized experimental data collected from 4 pigs under the protocol approved by the IACUC at University of North Carolina. Each pig received vasopressor at several distinct infusion rates under general anesthesia and mechanical ventilation to elicit a wide range of HR response. The HR was measured per second and pre-processed with a 50-point median filter to remove noise and outliners. He then performed system identification of the new dose-response model (8) to minimize the prediction errors. The results in Fig. 4 indicate the effectiveness of the proposed model. Figure 4: Measured versus model-predicted HR responses of 4 pigs of vasopressor for 4 pigs (one-compartment DDDR model). 19 2.5. Sensitivity Analysis As shown in the three different medications, the strength of the direct dynamic dose- response model relative to the classical PKPD model is that the former is more parsimonious than the latter. Utilizing the low-order direct dynamic dose-response model ????1? ? (?)for remifentanil ????(?) = ???0?50(?) 2 (?????(?) + ???(?)) as an example, it (?? ?50+?? (?)) is still nonlinear and does not lead to standard adaptive control design techniques based on linear model parameterization. To linearly parameterize the direct dynamic dose- response model, we proposed to derive a semi-adaptive model in which only the parameters having a large impact on the model?s fidelity are adapted while the remaining parameters are fixed at nominal values as determined by the system identification of remifentanil in Chapter 2.2. For this purpose, a sensitivity analysis was performed to examine the importance of the direct dynamic dose-response model parameters ??, ?50 and ? on the response of the model as follows (note that ??0 is not included in the analysis since it is a measured quantity). First, we perturbed the parameters of the direct dynamic dose-response model from their nominal values (as derived from the system identification in Chapter 2.2). With one parameter fixed at its nominal value as well as 25 % and 50 % away from the nominal value in both positive and negative directions (thus five values), the remaining two parameters were widely varied. The models were subject to an escalating dose of remifentanil with 5 distinct infusion rates that can yield a wide range of RR response in the nominal model. The infusion dose profile used for sensitivity analysis and the resulting response of the nominal direct dynamic dose-response model are shown in Fig.5 20 (i.e., the one characterized by the nominal values in Table 1 for the one-compartment model). Noting that the nominal value of ?50 is 0.08 mcg/kg/min, this dose profile could excite the direct dynamic dose-response model over a wide input-output range. 0.2 25 0.15 20 0.1 15 0.05 10 0 5 0 5 10 0 5 10 Time [min] Time [min] Figure 5: (a) Remifentanil infusion dose profile for in-silico sensitivity analysis. (b) The resulting RR response of the nominal direct dynamic dose-response model. Then the difference between the RR responses from the nominal model and the perturbed models were quantified as RMSE. The RMSE contours across the five different values of the fixed parameter were examined to assess the sensitivity of the model response on that parameter. Fig. 6 shows the contour plots of the discrepancy between the nominal and perturbed RR responses across different values of (a) ? and (b) ?50 (nominal value and ?50 % perturbations away from nominal value; ?25 % perturbations are not shown for brevity). It can be clearly observed that the change in the trend of contours with respect to ? is small (Fig. 6a), whereas the contour undergoes a large change with respect to ?50 (Fig. 6b). The contours associated with different values of ?? exhibited a change in trend larger than ? but smaller than ?50 (not shown). In sum, the sensitivity analysis revealed 21 Infusion Rate [mcg/kg/min] RR [bpm] that ? has negligible impact on the model?s response relative to ?50 and ??, meaning that it is difficult to faithfully identify ? and it is reasonable to fix it at its nominal value. (a) (b) Figure 6: The impact of (a) ? and (b) ?50 on the model?s response. Second, based on the results of the sensitivity analysis, the cooperativity constant ? in the direct dynamic dose-response model was fixed at its nominal value ??? determined by the mixed-effects modeling. The goodness of this "semi-adaptive" model was evaluated by 1) generating 100 random models and the corresponding simulated RR responses to the escalating dose of remifentanil via 100 uniformly random model parameter perturbations {??,?, ?50,?, ??}, ? = 1,? ,100. 2) identifying the semi-adaptive model for each of these random models that minimizes the prediction error norm: 22 {?? ??,?, ?50,?} = arg min ???(?; ??,?, ?50,?, ??) ? ??(?; ??,???, ??? 50,? , ???)? 2 (9) ??,?,?50,??? where {???,?, ? ? 50,?} are the optimal parameters of the semi-adaptive direct dynamic dose- response model corresponding to the ? -th random model, ??(?; ??,?, ?50,?, ??) is the simulated RR response of the ?-th random model, ??(?; ? ??,???, ?50,?, ???) is the RR response of the ?-th random model predicted by the semi-adaptive model with the parameters {? ??,???, ?50,?} , and finally, 3) evaluating the discrepancy between the responses of the original versus semi-adaptive models in terms of the distribution of the root-mean- squared prediction error (RMSPE), where ? is the number of data samples: 1 ????? = ???(?; ??,?, ?50,?, ??) ? ??(?; ? ??,???, ?50,?, ???)? (10) ?? 2 The RMSPE associated with the 100 random models tested were only 0.68 ? 0.45 m-1 (mean ? SD ) with the best-case and worst-case RMSPE of 0.02 m-1 and 1.55 m-1, respectively. Fig. 7 presents the (a) best-case and (b) worst-case fits between true (simulated) versus predicted (by the semi-adaptive model derived from (6)) RR responses among the 100 randomly generated direct dynamic dose response models. Therefore, the semi-adaptive direct dynamic dose-response model with ? fixed at ??? = 3.12 was capable of fitting the randomly generated direct dynamic dose-response models with wide- ranging parameter values (all ranging up to ?50 % perturbation from the respective nominal values listed in Table 1) faithfully. 23 25 25 True RR Response Predicted RR Response 20 20 15 15 10 10 True RR Response Predicted RR Response 5 5 0 2 4 6 8 10 0 2 4 6 8 10 Time [min] Time [min] (a) Best fit (RMSPE = 0.02 m-1) (b) Worst fit (RMSPE = 1.55 m-1) Figure 7: True RR response simulated by direct dynamic dose-response model and its prediction by semi-adaptive direct dynamic dose-response model. Therefore, ? can be fixed to its nominal value. Then the semi-adaptive direct dynamic dose-response model can be transformed to a linearly parameterized model as follows. First, (2) can be rewritten as: ? ??0 ? ?? ? = ? (9) ? 50 ? ? ?50? ?? ? ?? ??? where ? = ? 0 . Then, (1) can be rewritten as follows in the time domain where ? = ?? ?: ?? ??? = ???? + ? (10) ?50 ? which can be linearly parameterized with respect to ?? and ? , and thus, ?? and ?50 can ?50 be adapted online based on the measured dose-response data. 24 RR [bpm] RR [bpm] Chapter 3: SISO Semi-Adaptive Control Design In this chapter, a SISO semi-adaptive controller is introduced based on the linearly parameterized model in chapter 2. The ?semi-adaptation? implies one parameter is fixed at its nominal value to derive the linearly parametrized model. Two adaptive control methods are presented in this chapter: direct model-reference adaptive control (MRAC) and composite adaptive control (CAC). Comparative simulation studies were conducted to demonstrate the advantages of the proposed methods over non-adaptive method. The tracking and estimation results of two proposed adaptive controllers were also compared. 3.1. Model-Reference Adaptive Control (MRAC) The desired performance of the closed-loop controlled dose-response system is specified by a 1st-order reference model: ???? = ????? + ??? (11) where ?? = ?? are positive constants, and ? is a bounded external reference signal. The objective is to formulate adaptive control laws to force ? to converge to ?? asymptotically: lim ?(?) ? ??(?) = 0. ??? Consider the following control law: ? = ????? + ????? ? ?(? ? ??) = ????? + ????? ? ?? (12) where ???? and ???? are variable feedback gains, ? is a constant feedback gain, and ? = ? ? ? ? ?? is tracking error. Note that (12) with ???? = ? 50 ? = ?? and ???? = ?? = (?? ? ? 50 ? ? ) ? ?? 25 leads to a perfect matching of the plant dynamics to the reference model under the full knowledge of the plant model parameters: ?? ?? ?50 ?50 ??? = ???? + ? = ???? + [?? ? + (?? ? ??) ? ? ??] ?50 ?50 ?? ?? (13) = ???? + ??? In reality, the plant model parameters ?? and ?50 are unknown. The difference between the ideal (?? and ??) versus estimated (???? and ????) controller parameters are defined as follows: ?50 ???? ? ?? ???? ? ??? ? [ ? ??? ] = (14) ? ?50???? ? (?? ? ? ) [ ? ?? ] Then the dynamics of the tracking error when the plant (10) is subject to the control law (12) becomes: ?? ??? = ??? ? ???? = ???? + (? ? + ? ? ? ??) + ? ? ? ? ? ? ??? ??? ? ? ?50 (15) ?? ?? ?? ?? = ?? ??? + (????? + ?????) ? ?? = ?? ? ? ?? + ??? ? ?50 ? ? 50 ?50 ?50 where ? = [? ?]?. To derive the adaptive laws for the controller parameters ???? and ????, ? 0 consider the following Lyapunov function candidate, where ? = [ 1 ] and ? , ? > 0 0 ? 1 22 are the adaptive gains: 1 1 ?? ? = ?2 + ??????1??? (16) 2 2 ?50 The time derivative of the Lyapunov function candidate is: ?? ?? ?1 ? ?? ????? = ???? + ??? ? ???? = ? (???? + ??? ?? ? ??) + ??????1???? (17) ?50 ?50 ?50 ?50 26 Therefore, the time derivative of the Lyapunov function candidate can be made negative semi-definite by designing the adaptive laws in (18): ????? ???? = [ ] = ???? (18) ????? Indeed: ?? ??? = ?? 2 2?? ? ?? ? 0 (19) ?50 Since ? is positive definite and ??? is negative semi-definite, ? is bounded. Thus, the dose- response system (10) with the adaptive control law (12) and (18) is globally stable and the signals ? , ???? and ???? are bounded. The global asymptotic convergence of ? is guaranteed by the Barbalat?s Lemma: due to the boundedness of ?, ???? and ????, ??? in (15) is bounded. This in turn means that ??? is also bounded, and as a consequence, ??? is uniformly continuous. Hence, lim ??? = 0, and according to (19), lim ?(?) = 0. ??? ??? 3.2. Composite Adaptive Control (CAC) CAC is motivated by the recognition that the information on the plant model parameters are reflected in both the error between commanded versus actual plant outputs (i.e., tracking error) as well as the error between actual versus model-derived plant outputs (i.e., prediction error). Thus, CAC strive to estimate plant model parameters from both errors. In general, CAC exhibits superior control performance and faster parameter convergence relative to direct MRAC [48]. To derive the CAC law, the plant dynamics is rewritten as follows: ?? ??? + ??? = ?(?? ? ??)? + ? (20) ?50 27 1 Filter (20) using the low-pass filter to avoid differentiation yields: ?+?? 1 ?50 ?50 (?? ? ??) ?(?) ? ?(?) = ?(?) + ?(?) ? + ?? ?? ?? ? + ?? (21) ?? ?? = ?(?) + ?(?) = ??? ?? ? + ?? ? ? ? where ? is the Laplace operator, ? = [ ] while ?? and ?? are defined in (14). ?? ?+?? The adaptive law for the parametric model (21) based on the gradient algorithm [49] is given by: ????? [ ] = ????? (22) ????? where ? = (????? ? ?) is the prediction error and ? > 0. The control law for CAC is given by (12) together with the adaptive law incorporating both (18) and (22): ????? [ ] = ???? ? ???? (23) ????? The stability of the CAC-based closed-loop dose-response dynamics can be shown by using the Lyapunov function candidate (16). The time derivative ??? of ? is given by: ? 2 ? ?? ?? ?? ??? = ???? ? ?? 2 ? ??????????? = ???? 2 ? ??2 ? ??2 ? 0 (24) ?50 ?50 ?50 ?50 Since ? is positive definite and ??? is negative semi-definite, ? is bounded. Thus, the dose- response system (10) with the adaptive control law (12) and (18) is globally stable and the signals ?, ?, ???? and ???? are bounded. The global asymptotic convergence of ? and ? is guaranteed by the Barbalat?s Lemma: ??? in (15) and ?? = ?????? ? ??????? are bounded due to the boundedness of ?, ?, ???? and ????. This in turn means that ??? is also bounded, and as a 28 consequence, ??? is uniformly continuous. Hence, lim ??? = 0 , and according to (24), ??? lim ?(?) = 0 and lim ?(?) = 0. ??? ??? 3.3. In-Silico Implementation and Simulation To test the validity of the adaptive controllers designed above, in-silico simulation was conducted under the following conditions: 1) remifentanil concentration is 10 mcg/ml; 2) maximum infusion rate is 1.2 l/h; 3) RR0 is set at 25 bpm; 4) target reference RR is set at 15 bpm; 5) ?? = ?? are set to achieve a 5 % settling time of 3.74 min; and 6) sampling/control update interval is 100 ms. The MRAC and CAC were applied to the direct dynamic dose-response models derived for the 24 in-silico subjects (Chapter 2.2) in setting. The initial parameter estimates were set at their nominal values. The validity and performance of the adaptive controllers were examined by quantifying the distribution of the response characteristics including the settling time, the error norm between ? and ?? , and the steady-state ?? regulation error. For the purpose of comparison, a non-adaptive controller, derived by setting the model parameters ?? and ?50 to their nominal values (Table 1) and the adaptive gains to zero in the MRAC and CAC, was also applied to the models of the 24 in-silico subjects. Then, the performance of the non-adaptive controller was evaluated relative to its adaptive counterparts using the same response characteristics. Table 3 summarizes the performance of the adaptive and non-adaptive controllers on the 24 in-silico subjects represented by the one-compartment direct dynamic dose-response model (Chapter 2.2). Overall, both MRAC and CAC boasted adequate performance in 29 regulating the RR response at the commanded set point. In particular, the adaptive controllers outperformed the non-adaptive controller in terms of the consistency of the RR response with respect to the target response dictated by the reference model (11) in both transient and steady states. The settling time across the 24 in-silico subjects associated with the adaptive controllers was very close to the reference value of 3.74 m with a small interquartile range (IQR), which contradicts against the non-adaptive controller whose settling time deviates largely from the reference value with a large IQR. The superior capability of the adaptive controllers in following the reference model relative to the non-adaptive controller is also clearly seen from the RMSE between ?(?) and ??(?) during the transient state (0 ? ? ? 5 min ). Further, the steady-state error associated with the adaptive controllers was persistently close to zero across all the 24 in- silico subjects, whereas the non-adaptive controller resulted in variable and non-zero steady-state error. It is noted that the absolute magnitude of steady-state RR errors observed in the in-silico simulation was not too large. However, it may be exacerbated in real scenarios by finite RR resolution due to quantization. Table 3: Performance of adaptive and non-adaptive controllers (median (IQR)). 5% Settling Time [min] RMSE between ? and ?? [m -1] Steady-State Error [m-1] 100ms/median(IQR) (Reference = 3.74 min) (0 min ? t ? 5 min) (Reference = 0 m-1) MRAC 3.75 (0.03) 0 0.18 (0.55) 1 0.00 (0.00) 1 CAC 3.74 (0.26) 1 0.18 (0.46) 1 0.00 (0.00) 1 Non-Adaptive Control 3.82 (0.69) ? 0.34 (0.47) ? 0.11 (0.13) ? ?: p<0.05 (paired t-test). CAC and MRAC showed comparable control performance in terms of the response characteristics listed in Table 3. However, CAC modestly outperformed MRAC in terms 30 of parameter estimation (see Fig. 9 for a representative example). Indeed, among the 24 ? in-silico subjects, CAC resulted in more accurate estimation of ?? = ? 50 ? and ? =? ?? ? (? ? ? ) 50? ? than MRAC in 15 subjects (note that the accuracy was assessed with ?? respect to {?? ??,?, ?50,?} derived from (9) for all 24 subjects, i.e., ? = 1,? ,24). Thus, it can be concluded that the adaptation drive contributed by the prediction error, in addition to the one contributed by the tracking error, is in general beneficial in adapting the dose- response model to each subject. However, the absolute accuracy of the estimated parameters was not superb, which may be attributed, at least in part, to the constant target set point that lacks in persistent excitation property [48]. MRAC Composite Adaptive Control Non-Adaptive Control 26 26 26 Reference Model RR Reference Model RR Reference Model RR 24 Target RR 24 Target RR 24 Target RR RR (24 Subjects) RR (24 Subjects) RR (24 Subjects) 22 22 22 20 20 20 18 18 18 16 16 16 14 14 14 0 5 10 0 5 10 0 5 10 4 4 4 3 3 3 2 2 2 1 1 1 0 0 0 0 5 10 0 5 10 0 5 10 Time [min] Time [min] Time [min] Figure 8: Closed-loop controlled RR responses and remifentanil infusion rates associated with (a) MRAC, (b) CAC and (c) non-adaptive control (RR set point = 15 bpm). 31 RR [bpm] Infusion Rate [ml/min] Figure 9: Representative in-silico results associated with MRAC and CAC. 10 10 5 5 0 0 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 Normalized k Error Normalized I Error e 50 Figure 10: Distribution of parameter estimation errors in the steady state. (a) ??. (b) ?50. 32 Histogram In-depth scrutiny of the estimated parameters revealed that ?50 may be accurately estimated regardless of the inaccuracy in estimating ?? and ??. The distribution of the errors associated with ?? and ?50 across the in-silico subjects (normalized by their respective true values) are shown in Fig. 9, which demonstrate that overall ?50 was estimated more reliably than ?? in the steady state. This observation may be explained in two ways. First, the results of the sensitivity analysis in Chapter 2.5 show that the response of the direct dynamic dose-response model is largely affected by ?50 (as evidenced by Fig. 6b). Thus, it intuitively makes sense to anticipate that ?50 can be faithfully estimated. Second, and perhaps more obviously, (13) indicates that ?50 is nothing but a ratio between ? and ? in the steady state: ?? ?(?) ? ??(?) ??? = ???? + ? ? = ?(?) ? ? = ? (25) ?50 ?(?) ?? 50 0 ? ??(?) Thus, ?50 can be readily estimated accurately from the remifentanil infusion rate ? RR response data in the steady state. The knowledge of ?50 may be valuable in securing subject safety in the closed-loop control of remifentanil by preventing an infusion rate that is excessively higher than subject-specific therapeutic range. 33 Chapter 4: Two-Input Two-Output Semi-Adaptive Control Design In this chapter, the SISO semi-adaptive controller is extended to a two-input two-output case. Different from the SISO design, the coordination of targets also presents a challenge. This is because the reference targets are not always achievable in all subjects due to the inter-individual variability in dose-response relationships and the bounds on medication dose to ensure subject safety. Because we only want to focus controller design in this chapter and Chapter 5, so we assume the initial set targets are achievable for in-silico subjects for simplicity in these two chapters. A more general case in which the initial set targets are not achievable will be presented in Chapter 6. In the following, we will first introduce a dose-response model for two interacting medications consists of a low-order mixing model in chapter 2 and a response surface model. Then the control formation for it will be presented. Comparative simulation studies were also done to demonstrate the advantages of the proposed method over non- adaptive method. 4.1. Dose-Response Model for Two Interacting Medications A dose-response model for two interacting medications was derived by combining a low- order mixing model in chapter 2 and a response surface model reported in Minto et al. [33] (Fig. 11). First, the low-order mixing model represents the relationship between the intravenous infusion rate and the (hypothetical) infusion rate at the site of action of a medication: 34 ???1 = ???1?1 + ??1?1 (26) ???2 = ???2?2 + ??2?2 where ?1 and ?2 are the intravenous infusion rates associated with the two medications ?1 and ?2, ?1 and ?2 are the corresponding infusion rates at the sites of action, and ??1 and ??2 are the equilibration constants. Second, the response surface model relates the two infusion rates to the system outputs (i.e., clinical endpoints): ? ?1 ? 1 ? + 2 ( 50,11 ?50,12 )1 ? ? ? + ? ?2 ?1 1 1 1 1 ?50,11 ?1 = ?10 1 ? ? ? ? , ?1 = 1 2 1 ?1 ?2 ? + ? + 1 + ( 50,11 50,12 ?50,11 ?50,12 1 ? ? ? + ? ?2 ) 1 1 1 [ 1 ] (27) ? ? ?1 2 ? + 2 ( 50,21 ?50,22 )1 ? ? 2 ?1 2 ?2 + ?2?2 ?50,21 ?2 = ?20 1 ? ? ? ? , ?2 = 1 2 2 ?1 ?2 +? + 1 + ( 50,21 ?50,22 ?50,21 ?50,22 1 ? ? ? + ? ?2 ) 2 2 2 [ 2 ] where ?1 and ?2 are the system outputs, ?10 and ?20 are the baseline system outputs before the infusion starts, ?50,?? , ?, ? = 1,2 are the infusion rate of medication ?? associated with 50 % change in the system output ? , and ?1 and ?2 the cooperativity constants. The functions ?1 and ?2 denote the relative dominance of the infusion rates associated with the two medications (?1 = ?2 = 1 if only ?1 is infused and ?1 = ?2 = 0 if only ?2 is infused). The parameters ?1 and ?2 represent the degree of synergistic interaction between the medications associated with the outputs ?1 and ?2 (0 < ?1, ?2 < 4, where 0 and 4 correspond to zero and maximum interactions, respectively). 35 Dose-Response Synergy [Eq. (2)] u1 Mixing [Eq. (1)] ? ? Ieq,1 yx 1 + 2 1 (Medication 1) (Medication 1) 1 ?50,11 ?50,12 (Output 1) 1 ? ? 21?1 + ?1?1 Ieq,1 Dose-Response Synergy [Eq. (2)] u2 Mixing [Eq. (1)] ? ? Ieq,2 yx 1 + 2 2 (Medication 2) (Medication 2) 2 ?50,21 ?50,22 (Output 2) 1 ? ?2?2 + ?2? 2 2 Ieq,2 Figure 11: Dose-response model for two interacting medications, consisting of a low- order mixing model and a response surface model. To derive a control-oriented input-output model relating ?1 and ?2 directly to ?1 and ?2, we transform (27) as follows: ? ? 1 1 2 ? +?0 ? ?? ?? ?50,?1 ?50,?2 ?? ? ( ) = 2 , ? = 1,2 (28) ?? 1 ? ????(?1, ?2) + ???? (?1, ?2) which exhibits a large nonlinearity. To facilitate the control design, linearizing (28) around an operating point (?10, ?20) yields: ??? ??? ?? ? ??(?10, ?20) + | (? ? ? ) + | (? ? ? ) ?? 1 10 2 201 ( ???10,? 220) (?10,?20) ? 2 2 2 ( 10 ? + 20 ?10 ?) [( ( ) 20 ?10 ?20 ?50,?1 ?50,?2 ? ) + 1 + ?? (? ) + 2(1 ? ??) ( ) ( )]50,?1 50,?2 ?50,?1 ?50,?2 ?1 = 2 ? 2 ? 2 ?50,?1 (( 10 ) + ( 20 ? ) + (2 ? ? ) ( 10 ? ) ( 20 )) ? ? (29) 50,?1 ?50,?2 ?50,?1 ?50,?2 ? ? 2 ? 2 ? 210 20 ? ?(? + ? ) [ (1 + ??) ( 10 20 ? ) + (? ) + 2(1 ? ??) ( 10 ) ( 20 )] 50,?1 50,?2 50,?1 50,?2 ?50,?1 ?50,?2 ?2 + 2 ? 2 ? 2 ? ? ? 50,?2 (( 10? ) + ( 20 ? ) + (2 ? ? ) ( 10 ? ? ) ( 20 )) 50,?1 50,?2 50,?1 ?50,?2 36 Assuming ?1 and ?2 are primary and secondary medications, respectively, (29) reduces to the following at (?10, ?20) = (?10, 0): ?1 ?2 ?1 ?2 ?1 ?2 ?1 = + (1 + ?1) = + ?1 ? + ?50,11 ?50,12 ?50,11 ?50,12 ?11 ?12 (30) ?1 ?2 ?1 ?2 ?1 ?2 ?2 = + (1 + ?2) = + ?2 ? + ?50,21 ?50,22 ?50,21 ?50,22 ?21 ?22 Therefore, ?1 and ?2 can be expressed by ?1 and ?2 as follows: ?11?21 ?1 = (? ? ? ? ? ) ?12?21 ? ?11? 12 1 22 2 22 (31) ?12?22 ?2 = (? ? ? ? ? ) ?11?22 ? ? ? 11 1 21 2 12 21 Now, differentiating (30) in time yields: 1 1 ???1 = (???1?1 + ??1?1) + (??? ? ?2 ?2 + ??2?2) 11 12 (32) 1 1 ???2 = (???1?1 + ??1?1) + (???2?2 + ??2?2) ?21 ?22 Finally, substituting (31) into (32) yields the following input-output model between ?1 and ?2 versus ?1 and ?2 for control design, where ? ? ?12?21 ? ?11?22: 1 ??1 1 ???1 = (???1?12?21 + ??2?11?22)?1 + ?1 + (?? ? ? ?1 ?22?21 ? ??2?21?22)?2 11 ??2 + ? ? 2 12 (33) 1 ??1 1 ???2 = (???1?12?11 + ??2?11?12)?1 + ?1 + (??1?22?11 ? ??2?21?12)?? ?21 ? 2 ??2 + ? ? 2 22 which leads to the following linearly parameterized model: 37 ??? ?1 ?1 ??? = [ 1] = ? [? ] + ? [? ] (34) ???2 2 2 ??1 ??2 1 ?? ? = [ ?1 ?12?21 + ??2?11?22 ??1?22?21 ? ??2?21?22 ? ? where ] and ? = [ 11 12]. ? ???1?12?11 + ??2?11?12 ??1?22?11 ? ? ? ? ??2 21 12 ?1 ??2 ?21 ?22 The model in (34) is formulated in terms of the transformed (by (28)) outputs ?? rather than the actual ones ??, ? = 1,2. To derive ?1 and ?2 from ?1 and ?2, ?1 and ?2 must be known. In chapter 2, we have shown that ?? makes relatively small influence on the system output compared with ??? and ?50,?? ( ?, ? = 1,2 ) and may thus be fixed at a nominal value [50]. With ?1 and ?2 specified a priori, the model (34) is a linearly parameterized control design model with the elements of the matrices ? and ? as unknowns to be adapted on-line. As the proposed control approach adapts only a subset of the plant model parameters (i.e., parameters other than ?1 and ?2), it is semi-adaptive rather than fully adaptive. 4.2. Control Design Consider the following 1st-order reference model specifying the ideal endpoint responses of the subject to a reference target: ???? = ???? + ??? (35) ? 0 where ?? = ?[ ? ] < ? is a negative definite matrix, ? = ?? 0 ? ? ? , and ? is a ? ?1 bounded reference target: ? = [? ] . Then, the objective is to formulate an adaptive 2 control law that can guide ? to ?? asymptotically: lim ?(?) ? ??(?) = ?. Consider the ??? following model reference adaptive control (MRAC) law: 38 ?1 ?1 ?1 ?1 ? ?1? ? = [? ] = ???? [? ] + ???? [? ] + ? [? ? ? ] = ????? + ????? + ?? (36) 2 2 2 2 2? ? 0 where ???? and ???? are variable feedback gains, ? = ? [ 1 ] < ? a negative definite 0 ?2 constant feedback gain, ?1 and ?2 the reference targets associated with the system outputs, and ? the tracking error. Note that the MRAC law (36) equipped with ? ????? = ?? = ? ?? and ? ????? = ?? = ? (?? ? ?) results in perfect matching of the plant dynamics (34) to the reference model (35): ??? = ?? + ?? = ?? + ?[???(?? ? ?)? + ? ????? + ??] (37) = ??? + ??? + ??? However, the plant model parameters ?? and ?? are not known a priori in reality. Defining ? ? [?? ??] and ??? ? [???? ????] , the difference between the true versus estimated parameters is defined as follows: ??? = ??? ? ? = [???? ????] ? [???? ? ?? ???? ? ??] (38) Then, the dynamics of the tracking error when the plant (34) is subject to the MRAC law (36) becomes: ??? = ??? ? ???? = (? + ?????)? + ?????? + ??? ? ???? ? ??? = (? + ????? ? ??)? + (????? ? ??)? + ??? + ??? (39) = ?????? + ?????? + ??? + ??? ? = ?[???? ????] [ ] + ??? + ??? = ????? + ??? + ??? ? ? where ? = [ ]. ? 39 Assuming that all the leading principal minors of ? are positive, there exists a positive ? ? ?1 ? ? ?2?1 0 1 ? 1 definite matrix ? = [ ] > ? such that ?? = [ 11 ?12 0 ? ? ] > ? holds. Noting 2 ? ?1 ? 2 ? ?2 ? 221 ?22 that all the elements in ? assume positive values due to the physical meaning of the parameters ??? and ??? , ?, ? = 1,2 , there exist ?1 and ?2 making ?? symmetric, i.e., ? ? satisfying ? ?2 = ? ?11 2 . Hence, there exists a matrix ? that satisfies ?? = ? ??. Finally, ?12 ?21 consider the Lyapunov function ? = ???? + ??[??????? ? ?? ??? ? ] and the adaptation law ????? = ?? ???? , where ??[?] is the trace of the argument and ?? is a positive definite adaptation gain matrix. Then, the time derivative of the Lyapunov function candidate reduces to the following: ??? = ?????? + ?????? + ?? [????????????? ? ? + ??????????????? ] = 2?????? + 2?? [??????? ? ?? ???? ? ] = 2???(????? + ? ? + ???) ? 2??[??????? ?? ? ] (40) = 2???????? + 2?????? + 2? ????? ? 2????????? = 2?????? + 2? ????? ? 0 Since ? is positive definite and ??? is negative semi-definite, ? is bounded. Hence, the plant model (34) with the MRAC law (36) and the adaptation law ????? = ????? ?? ? is globally stable, and accordingly, ?, ???? and ???? are bounded. Then, ??? is bounded from (39), and, ??? is also bounded from (40). As a consequence, ??? is uniformly continuous. Therefore, lim ??? = 0, and lim ?(?) = ?. ??? ??? 40 To prevent the drift in ???? and ???? while ?(?) ? ? (which frequently occurs in case the system is regulated at constant set points), we employed the dead zone [48], a scheme to stop parameter adaptation when ?(?) ? ?: ? ????? ?, |?1(?)| > ????? = { ? or |?2(?)| > ?1 ?2 (41) ?, otherwise 4.3. In-Silico Implementation and Simulation To evaluate the proposed semi-adaptive control approach, we considered an example scenario in which cardiac output (?1; ?1(0)=2.0 lpm) and respiratory rate (?2; ?2(0)=15 bpm) are regulated via infusion of propofol (?1) and remifentanil (?2) in an in-silico simulation setting. We used the nonlinear dose-response model in (26)-(28) to simulate the dose-response relationships of a in-silico subject (note that our control design was performed based on the linear model (34), and the controller is subject to the structural uncertainty due to the modeling error originating from linearization). We specifically simulated pediatric patients based on the available dose-response data, by first setting the ranges of the model parameters and then creating a cohort of random subjects. For the parameters associated with the dose-response relationships for propofol, we derived the ranges of ??1, ?50,11, and ?1 by analyzing the experimental data we obtained from swine subjects with body weight ranging 25-30 kg in Chapter 2.3. Then, we derived the ranges of the remaining parameters associated with propofol by assuming that (1) these ranges translate to pediatric patients of comparable body weight and (2) ?50,11 ? ?50,21 [51]. The selected ranges for the parameters were: 0.02 ? ??1 ? 0.10 min -1, 1 ? ?1 ? 5, 0.2 ? ?50,11 ? 0.6 mg/kg/min, and 0.2 ? ?50,21 ? 0.6 mg/kg/min. 41 For the parameters associated with the dose-response relationships for remifentanil, we derived the ranges of ??2 , ?50,22 , and ?2 in Chapter 2.2, while we assumed that the influence of remifentanil on cardiac output is relatively small (translating to a large ? -150,12). The selected ranges for the parameters were: 0.10 ? ??2 ? 0.50 min , 1 ? ?2 ? 5, 0.12 ? ?50,12 ? 0.36 mcg/kg/min, and 0.04 ? ?50,22 ? 0.12 mcg/kg/min. For the parameters associated with the inter-medication interaction (i.e., ?1 and ?2), we simply selected a wide range to simulate diverse inter-medication interaction: 1.5 ? ?1, ?2 ? 3. Based on the ranges of the model parameters selected above and assuming that all the model parameters exhibit uniform distributions within the selected ranges, we created in- silico subjects by selecting a random value for each parameter from the respective range. Because we assume the initial targets are achievable in this chapter, 20 in-silico subjects who can achieve the reference targets specified are selected for the in-silico simulation. In summary, the parameters associated with the semi-adaptive control were as follows: 0.8 0 5 0 ?? = [ ] = ??? , ? = [ ] , ? = ? , ? = 0.01 , ? = 0.01 , ? =0 0.8 0 5 ? 4?4 ?1 ?2 1? 0.8mg/kg/min, ?2? = 0.07 mcg/kg/min. In the in-silico evaluation, nominal cardiac output and respiratory rate before medication infusion were set at 3.0 lpm and 25 bpm. The initial model parameter values in the controller were set at the respective average values from all the in-silico subjects. For 42 control computation, a sampling rate of 1 Hz was used. Considering that higher infusion rates (e.g., bolus infusion) are required for reference target tracking during the initial transients, we allowed higher propofol (4 mg/kg/min) and remifentanil (0.36 mcg/kg/min) infusion rates during the first 10 min after the control action started. Then the performance of semi-adaptive control was compared with the non-adaptive control. Fig. 12 presents the in-silico simulation testing results associated with the semi-adaptive and non-adaptive controllers in the 20 in-silico subjects with achievable reference targets. The semi-adaptive controller was superior to its non-adaptive counterpart both during the transient (0 ? ? ? 5 min) and steady state. On the average, the semi-adaptive controller could reduce the transient reference tracking errors associated with cardiac output and respiratory rate by 15% and 21%, respectively (in terms of root-mean-squared error (RMSE)); and as well, the steady-state errors (in the absolute sense) by 5% and 83%, respectively. In addition, the variances associated with all these errors were consistently smaller in semi-adaptive than non-adaptive controller. Considering that the main rationale underpinning the development of semi-adaptive control is to ensure robust performance against large inter-individual variability in dose-response relationships, the efficacy of the semi-adaptive control can be deemed satisfactory if not excellent. 43 Semi-Adaptive Control (MRAC) Non-Adaptive Control 3 25 3 25 2.5 20 2.5 20 2 15 2 15 1.5 10 1.5 10 0 10 20 30 0 10 20 30 0 10 20 30 0 10 20 30 4 0.2 4 0.1 3 0.15 3 2 0.1 2 0.05 1 0.05 1 0 0 0 0 0 10 20 30 0 10 20 30 0 10 20 30 0 10 20 30 TTiimmee [[mmiinn]] TTiimmee [[mmiinn]] TTiimmee [[mmiinn]] TTiimmee [[mmiinn]] Figure 12: In-silico simulation testing results associated with the semi-adaptive and non- adaptive controls in the 20 in-silico subjects with achievable reference targets. 44 Propofol [mg/kg/min] Cardiac Output [lpm] Remifentanil [mcg/kg/min] Respiratory Rate [bpm] Propofol [mg/kg/min] Cardiac Output [lpm] Remifentanil [mcg/kg/min] Respiratory Rate [bpm] Chapter 5: Switching Two-Input Two-Output Semi-Adaptive Control Design In Chapter 4, the model is only based on one operating point. However, the coupling effect between multiple medications frequently exhibits complex behaviors due to the regime-dependent inter-medication synergy, the model based on one operating regime is not suffice for trustworthy controller design. Therefore, a two-model switching control technique was developed based on two dose-response models associated with two distinct operating regimes in this chapter. In the following, a global nonlinear dose-response model applicable to two interacting medications derived from the global model is first presented. Then the semi-adaptive switching control approach to the model is described. In the end, the results obtained from in-silico testing of the semi-adaptive switching control are discussed. 5.1. Control-Oriented Dose-Response Models For multiple medication infusion, a confounding factor is the bi-phasic behavior elicited by the inter-medication synergy terms, which makes the input-output relationship ? dependent on the operating regime (i.e., on the magnitudes of 1 ? and 2 relative to ?50,?1 ?50,?2 ? ? each other). More specifically, if ?1 is used as the primary therapy ( 1 > 2 ), an ?50,?1 ?50,?2 increase in the co-infusion of ?2 relative to ?1 will effectively depress ?1 and ?2 by improving the inter-medication synergy (i.e., the denominator in ???,?, ? = 1,2 decreases, which adds up to an increase in the numerator in ???,?, ? = 1,2 caused by an increase in 45 ? ? ?2). In contrast, if ?2 is used as the primary therapy ( 1 < 2 ), a further increase in ?50,?1 ?50,?2 the infusion of ?2 relative to ?1 will not be effective in depressing ?1 and ?2 due to the weakening of the inter-medication synergy (i.e., the denominator in ???,? , ? = 1,2 increases, which cancels an increase in the numerator in ???,? , ? = 1,2 caused by an increase in ?2 ). Therefore, we proposed to linear the model at 2 different operating points. First, same as the model in Chapter 4, assume that ?1 and ?2 are secondary and primary medications, respectively. Reducing Eq. (29) at (?10, ?20) = (0, ?20) , differentiating it in time, and expressing ?1 and ?2 by ?1 and ?2 yields: 1 [ ] [ ] [ ] [ ] ??1 ???1 = (?? ? 1 ? 1?1 12 21 + ? 1 1 ? ?2 ?11?22) ?1 + ?[ ? 1 ] 1 11 1 [1] [1] [1] [1] ??2+ (??1?22?21 ? ??2?21?22) ?? 2 + ? [1] 2 ?12 (41) 1 [1] [1] [1] [1] ??1 ???2 = (???1?12?11 + ??2?11?12) ?1 + ?? [1] 1?21 1 [ ] [ + (? ? 1 ? 1 ] [ ] [ ] ??2 ?1 22 11 ? ? 1 ?2?21? 1 ? 12 ) ?2 + ? [ ] 2 ? 122 [1] ?50,11 [1] [1] ?50,21 [1] [1] [1] [1] [1] where ?11 = , ?12 = ?50,12 , ?21 = , ?22 = ?50,22 , and ? ? ?12?21 ? ?11?22 , ?1 ?2 which leads to the following linearly parameterized model: ??? ?1 ?1 ??? = [ 1] = ? ??? 1 [? ] + ?1 [2 2 ? ] (42) 2 [ ] [ ] [ ?? ? 1 ? 1 + ? ? 1 ] [1] [1] [? ? ? ? 1 ] [ ] [ ] 1 ?1 ?2 ?1 ? ??2? 1 ? 1 where ? 12 21 11 22 22 21 21 221 = [ ] and ? =? [1] [1] [?? ? ? + ? ? 1 ] [ ] [ ] [ ] [ ] [ ] 1 ?1 12 11 ?2 11? 1 12 ??1? 1 22? 1 1 1 11 ? ??2?21?12 ??1 ??2 [1] [1] ? [ 11 ?12 ? ]. ?1 ??2 [1] [1] ?21 ?22 46 Second, assume that ?1 and ?2 are primary and secondary medications, respectively. Reducing Eq. (29) at (?10, ?20) = (?10, 0), differentiating it in time, and expressing ?1 and ?2 by ?1 and ?2 likewise yields: ??? ?1 ?1 ??? = [ 1] = ? [ ] + ? [ ??? 2 ? 22 2 ? ] (43) 2 [2] [2] [2] [2] [2] [2] [2] [2] 1 ???1?12?21 + ??2?11?22 ??1?22?21 ? ??2? ?where ? = [ 21 222 ] and ? =? [2] [2] [2] [2] [2] [2] [2] [2] 2???1?12?11 + ??2?11?12 ??1?22?11 ? ??2?21?12 ??1 ??2 [2] [2] ? ? [ ] [ ] ? [ ] [ ] ? [ 11 12] 2 2 50,12 2 2 50,22? ? , ? = ? , ? = , ? = ? , and ? = . ?1 ?2 11 50,11 12 ? 21 50,21 221 ?2 [2] [2] ?21 ?22 It is noted that the control-oriented models in Eq. (42) and Eq. (43) have linear dynamics with state and input coupling. The controllability of these locally linearized models can be easily determined based on the matrix pairs (?1, ?1) and (?2, ?2), respectively. In case inter-medication coupling is not too large (which is in fact a requisite assumption to streamline control design; see Remark 2 in Section 5.2) so that the off-diagonal terms are small compared to the diagonal terms in the matrices ??, ? = 1,2 in Eq. (42) and Eq. (43), the matrices ?? , ? = 1,2 have full rank, and therefore, the linearized models are controllable. Considering that ?1 and ?2 are related to ?1 and ?2 via monotonic transformation in Eq. (28), ?1 and ?2 are likewise controllable. The observability is trivial: given that both ?1 and ?2 are available for measurement, the linearized models are observable. It is also noted that the structure of the control-oriented dose-response models corresponding to the ?1-dominant dose-response model (called the ?1 model hereafter) 47 in Eq. (43) and the ?2-dominant dose-response model (called the ?2 model hereafter) in Eq. (42) is identical, but not the definition of the parameters therein. This illustrates a few notable points on the role of switching control: (1) its most beneficial role may be to assist fast adaptation of dose-response relationship by providing a one-time instantaneous updating of model parameters, based on their definition (specifically, how ?1 and ?2 are [1] [2] incorporated into ??? and ??? , ?, ? = 1,2), to facilitate online parameter estimation in case the operating regime and the dose-response model used in the controller are inconsistent; and (2) the role of switching may be redundant if the dose-response relationship can be readily tuned to the operating regime with the parameter adaptation alone. Given these points, it may be reasonable to anticipate that the instances in which switching control is the most effective are those associated with control set point changes invoking a large perturbation in the operating regime that necessitates the switching of the dose-response model. 5.2. Semi-Adaptive Switching Control To cope with the nonlinearity and inter-individual variability inherent in the dose- response relationship, we investigated a semi-adaptive switching control approach to the infusion of two interacting mediations. The key components of the approach are the switching control strategy and multivariable semi-adaptive control. Considering that a single linearized dose-response model representing a particular operating regime (e.g., ?2 model in Eq. (42) or ?1 model in Eq. (43)) may not well represent the global nonlinear dose-response behavior, it may be a reasonable idea to 48 equip the controller with an ability to switch its mode depending on the operating regimes. In this chapter, we investigated an approach in which the controller can switch its mode between the ?1 -dominant and ?2 -dominant regimes. Considering that the variables ?1 and ?2 in Eq. (27) indicate the relative dominance of the medications ?1 and ?2 on the process outputs (i.e., clinical responses) ?1 and ?2 , the overall relative ? +? dominance of the medication was determined by 1 2: (1) ?2 dominates the process in 2 ? +? case 0 ? 1 2 < 0.5, and a closed-loop controller designed with Eq. (42) must be used 2 to control the dose-response dynamics; and similarly, (2) ?1 dominates the process in ? +? case 0.5 ? 1 2 < 1, and a closed-loop controller designed with Eq. (43) must be used 2 to control the dose-response dynamics. Given the above switching strategy, the controller must estimate ?1 and ?2 to determine the correct mode to operate and perform mode switching. According to Eq. (27), a pre- requisite to estimate ?1 and ?2 is the estimation of ?50,??, ?, ? = 1,2. Given the matrices ?? and ?? in Eq. (42) and Eq. (43) as estimated by the semi-adaptive control technique (as described later in this section; here, ? = 1 if the current control mode is ?2-dominant mode, or 2 otherwise), our approach to estimate ?50,??, ?, ? = 1,2 is as follows. First, ??1 [?] and ??2 are computed from the matrix ?? as its eigenvalues. Second, ??? , ?, ? = 1,2 are [ ] computed from each element in ??. Third, ?50,??, ?, ? = 1,2 are computed from ? ? ?? , ?, ? = 1,2 with the knowledge of the parameters ??, ? = 1,2 (see Remark 1). 49 Remark 1: It is not practically feasible to estimate both ?? and ?50,??, ?, ? = 1,2 from the [?] knowledge of ??? , ?, ? = 1,2. To cope with this challenge, it was assumed here that a priori (yet imperfect) knowledge of ?? , ? = 1,2 is available to compute ?50,?? , ?, ? = 1,2 [?] from ??? , ?, ? = 1,2. To make up for the inaccuracy due to this assumption, the plant (i.e., dose-response) models in Eq. (42) and Eq. (43) were rewritten as follows to take the errors induced by this assumption explicitly into account in the course of control design: ???1 ?1 ?1??? = [ ] = ?10 [? ] + ?10 [? ] + ???? ? (44) 2 2 2 ???1 ?1 ?1??? = [ ] = ?20 [? ] + ?20 [? ] + ?? (45) ???2 2 2 where the matrices ??0 and ??0, ? = 1,2 are the matrices ?? and ?? in Eq. (42) and Eq. (43) constructed with the true ??1 , ??2 , and ?50,?? , ?, ? = 1,2 together with ?? , ? = 1,2 known a priori, while the vectors ?? and ?? are the errors due to the mismatch between Eq. (42) and Eq. (43) versus Eq. (44) and Eq. (45). In essence, Eq. (44) and Eq. (45) instead of Eq. (42) and Eq. (43) were employed in the control design process as described later in this section. The control switching strategy and semi-adaptive control described above were integrated into a semi-adaptive switching controller based on the Lyapunov stability analysis, by proving that the stability of the two semi-adaptive controllers designed for Eq. (44) and Eq. (45) respectively can be established with a common Lyapunov function [52], [53]. Same as Chapter 4, consider the following 1st-order reference model specifying the ideal clinical responses of the subject to a reference target: 50 ??? ??? = [ ?1? ] = ???? + ??? (46) ????2 ? 0 where ?? = ?[ ? ] is a negative definite matrix, ?? = ??? , and ? a bounded 0 ?? ?1 reference target: ? = [? ] . The control objective is to guide ? to ?? asymptotically: 2 lim ?(?) = lim ?(?) ? ??(?) = ?, while at the same time estimating the subject?s dose- ??? ??? response relationship in terms of ??0 and ??0 , ? = 1,2 . For the plant in Eq. (44), consider the following adaptive control law: ?11 ?1 ?1 ? = ?? = [? ] = ????? [? ] + ????? [? ] + ??? = ?????? + ?????? + ??? (47) 12 2 2 where ????? and ????? are variable feedback gains, and ??? a nonlinear feedback term to guarantee the robustness of the controller against the model uncertainty ??. Note that ????? and ????? the estimations of ? ?? ?? = ????? and ? ?? ?? = ??? (?? ? ???) . Defining the difference between the true versus estimated parameters as follows: ???? = ???? ? ?? = [????? ?????] ? [????? ? ??? ????? ? ???] (48) where ? ? [??? ?? ??] and ???? ? [????? ?????], then the dynamics of the tracking error when the plant in Eq. (44) is subject to the adaptive control law ? in Eq. (47) is given by: ??? = ??? ? ???? = ???? + ???? + ?? ? ???? ? ??? = ???? + ???(?????? + ?????? + ???) + ?? ? ??? ? ??? + ??? = (??? + ???????? ? ??)? + (???????? ? ??)? + ?????? + ?? + ??? (49) = ????????? + ????????? + ?????? + ?? + ??? ? = ???[????? ?????] [ ] + ?? ?? ??? + ?? + ??? = ???????? + ?????? + ?? + ??? 51 ? where ? = [ ]. Assuming that all the leading principal minors of ??? are positive (see ? 1 0 Remark 2), there exists a positive definite matrix ? = [ ] > ? such that ?? 0 ? 10 = ??1 ??2 [1] [1] ? [ 11 ?12 ? ? ] > ? holds. Because all the elements in ??? are positive due to the physical ? ?1 ? ?2[1] [1] ?21 ?22 meaning of the parameters ??? and ??? , ?, ? = 1,2, a ? can be found that makes ???? ? ? symmetric, i.e., satisfying ?2[1] = ? ?1 [1]. Thus, there exists a matrix ? that satisfies ???? = ?12 ?21 ???. Consider the Lyapunov function ? = ???? + ?? [?????? ??????? ?? ?] and adaptation law ????? ?? = ????? , where ??[?] is the trace of the argument and ?? is a positive definite adaptation gain matrix. Then, the time derivative of the Lyapunov function is calculated as follows: ??? = ?????? + ?????? + ?? [??????? ?? ??? ? ? ??? + ?? ? ?? ???? ? ??? ? ??? ] = 2?????? + 2?? [?????? ?? ???? ? ? ??? ] = 2???(? ? ? + ? ? + ? + ? ?) ? 2??[?? ????? (50) ?? ??? ?? ?? ? ? ??? ] = 2??????????? + 2? ???????? + 2? ???? + 2? ????? ? 2? ????????? = 2???? ? ?????? + 2? ??? + 2? ???? ?11 If the element-wise upper bound of ??? is given by ?? ? [? ], ??? in Eq. (50) can be 12 |??| made negative semi-definite by designing ??? = ?? ? =?0???? ? [?1sgn(??? 1 ) ?2sgn(?2)] ??, where ?0 is a value smaller than the smallest eigenvalue of ? ?0??? ??10. Indeed: 52 |??| |??| ???? ? ? ? ??? ?? = ?? ????? ? ? ?? ??? ? = ?|? |? ?0???? ? 0 ?0???? ? ? ? (51) ??? = 2????????? + 2? ???? + 2? ??? ??? ? ?2|? |?? + 2|? ?|? + 2??? ???? ? 0 Similarly, for the plant in Eq. (45), consider the following adaptive control law: ? 0 ?21 ? = ?? = [? ] = ?????? + ?????? + [ 1]? 0 ?? (52) 22 ? where ????? and ????? are the estimations of ? = ? ?? ?? ???? and ? ?? ?? = ??? (?? ? ???). By 1 0 noticing that ??? = ??? and ? ??? = ??? [ ] (see Eq. (42) and Eq. (43)), ??? = 0 ? ? 0 ? 0 [ 1] ??? and ??? = [ 1] ??? . Defining the difference between the true versus 0 0 ? ? estimated parameters as follows: ? 0 ???? = ???? ? ?? = [????? ?????] ? [????? ? ??? ????? ? ???] = [ 1]???? (53) 0 ? where ? ? [??? ?? ??] and ???? ? [????? ?????], then the dynamics of the tracking error when the plant in Eq. (45) is subject to the adaptive control law ? in Eq. (52) is given by: ??? = ??? ? ???? = ???? + ???? + ?? ? ???? ? ??? ? 0 = ? 1??? + ??? (?????? + ?????? + [ ]???) + ?? ? ?0 ? ? ? ??? + ??? ? (54) = (??? + ???????? ? ??)? + (???????? ? ??)? + ?????? + ?? + ??? = ????????? + ????????? + ?????? + ?? + ??? 53 ? = ???[????? ?????] [ ] + ?? ?? ??? + ?? + ??? = ???????? + ?????? + ?? + ??? ? 0 Then, with the same Lyapunov function and adaptation law ?????? = ????? ? [ 1], the 0 ? time derivative of the Lyapunov function is calculated as follows: ??? = ?????? + ?????? + ?? [??? ??? ??? ? ?? ???? ???? ? ??? + ??????? ??? ] 1 0 = 2?????? + 2?? [?? ??????????? ? ? [? ] ? ?] 0 ? (55) = 2???(???????? + ?????? + ?? + ???) ? 2??[?? ? ? ????? ? ] = 2??????????? + 2? ??? ? ?????? + 2? ??? + 2? ?? ? ? 2? ???? ?????? = 2???? ? ?????? + 2? ??? + 2? ???? ?21 Similarly to Eq. (51), if the element-wise upper bound of ??? is given by ?? ? [? ], ??? 22 in Eq. (55) can likewise be made negative semi-definite by designing ??? = |??| [? ?? ? = ?? 1 sgn(?1) ?2sgn(?2)] ? . ?0???? ? ?0???? ? In this way, ??? ? 0 can be established for both the plants in Eq. (44) and Eq. (45), and thus, the switching semi-adaptive controller is globally stable in both ?1-dominant and ?2-dominant regimes. Because ? is positive definite and ??? is negative semi-definite, ? is bounded. Hence, ?, ????? and ?????, and also ????? and ????? are bounded. Then, ??? is bounded from Eq. (49) and Eq. (54), and therefore, ??? is also bounded by Eq. (50) and Eq. (55). 54 Thus, ??? is uniformly continuous. Therefore, lim ??? = 0 according to the Barbalat?s ??? Lemma [48] and lim ?(?) = ?. ??? Remark 2: The assumption that all the leading principal minors of ??? , ? = 1,2 are positive may not be readily satisfied in all medication infusion problems. However, the assumption is valid in case the inter-medication coupling is not too large (so that the off- diagonal terms in ???, ? = 1,2 are not too large). Finally, the following practical modifications were made when implementing the semi- adaptive switching controller. First, a hysteresis was incorporated into the switching strategy in order to avoid chattering in the switching process. It was implemented so that ? +? the controller switches from ?? in Eq. (47) to ?? in Eq. (52) if 1 2 > 0.6, and likewise 2 ? +? switches from ?? to ?? if 1 2 < 0.4 . Note that the hysteresis does not impact the 2 system stability; considering that control switching is intended to promote fast adaptation of the controller through instantaneous parameter updating, delayed switching due to the hysteresis may only slow down the adaptation of the controller to individual subjects and operating regimes (and may modestly deteriorate the control performance), but it does not deteriorate the system stability. Second, the nonlinear feedback term ??? = [?1sgn(?1) ?2sgn(??? 2 )] ? ??, ? = 1,2 in the control law, which has a discontinuity at ? = ?, ?0??? was modified into a piecewise continuous nonlinear feedback ??? = ? ? [? ?( 11 ) ?2?( 2)] ?1? ??? 2?? ?? , ? = 1,2 in order to avoid chattering in the control input [48], ?0??? where ?(?) is the saturation function given by: 55 ? , |?| ? 1 ?(?) = { (56) 1 , otherwise and ????, ? = 1,2 specifies the width of the saturation. Note that this approximation, widely used in robust control (e.g., sliding mode control [48]) to suppress the adverse impact of uncertainties while avoiding control chattering, still keeps the tracking error bounded (i.e., |?1| < ?1? and |?2| < ?2? ), although it does not guarantee the convergence of the tracking error to zero. Third, upper and lower bounds were imposed on the control input ?: 0 ? ?? ? ????, ? = 1,2, where ???? is the upper bound associated with the medication ??. Specifically, the lower bounds were set to zero due to physical constraint (that medications can only be infused but not drawn), while upper bounds were set to ensure subject safety against over-dosing. Fourth, a dead zone scheme [48] was implemented in the adaptation law to prevent the drift in the process model parameters ????? and ?????, ? = 1,2 in the neighborhood of ? = ? (which frequently occurs in case the system is regulated at constant set points due to the lack of persistent excitation) while maintaining the system stability and error boundedness, by de-activating the parameter adaptation when ? is small: ? ? ?? ? ??? , |?1(?)| > ?? or |?2(?)| > ?? = { 1 ?2??? , ? = 1,2 (57) ? , otherwise In sum, Fig. 13 schematically illustrates the overall control architecture. 56 Fig. 13: Architecture of semi-adaptive switching control approach. 5.3. In-Silico Evaluation An array of in-silico evaluation was conducted to examine the semi-adaptive switching control approach to the infusion of two interacting medications. For this purpose, an example scenario was considered, in which cardiac output (CO) and respiratory rate (RR) were controlled through the infusion of a sedative propofol ( ?1 ) and an opioid remifentanil (?2). In the in-silico evaluation, the nonlinear dose-response model in Eq. (26)-(27) was used to simulate the process (i.e., in-silico subject?s dose-response dynamics), while the control design was conducted based on its linearization in Eq. (44) and Eq. (45). Hence, the structural uncertainty due to the modeling error was considered in the in-silico evaluation. By exploiting available dose-response data, synthetic pediatric subjects were simulated by first setting the appropriate ranges of the model parameters and then creating a cohort of random subjects. Details follow. For the parameters associated with the dose-response relationships of propofol, the ranges of ??1 , ?50,11 , and ?1 were obtained from the experimental data collected from swine subjects with body weight ranging 25-30 kg and by assuming that these ranges translate 57 to the pediatric subjects of comparable body weight used in our previous work [54] to derive the parameters associated with remifentanil (based on the hemodynamic similarity between swine and humans). The range of the ?50,21 was set by assuming that ?50,11 ? ?50,21 [51]. The ranges thus selected were: 0.03 ? ??1 ? 0.09 min -1, 1.5 ? ?1 ? 4.5 , 0.2 ? ?50,11 ? 0.6 mg/kg/min, and 0.2 ? ?50,21 ? 0.6 mg/kg/min. For the parameters associated with the dose-response relationships of remifentanil, the ranges of ??2, ?50,22, and ?2 were obtained from our previous work in pediatric subjects [54]. The range of ?50,12 was set by assuming that the influence of remifentanil on CO is relatively small compared with its effect on RR (thus a large ?50,12). The ranges thus selected were: 0.15 ? ??2 ? 0.45 min -1, 1.5 ? ?2 ? 4.5 , 0.1 ? ?50,12 ? 0.3 mcg/kg/min, and 0.04 ? ?50,22 ? 0.12 mcg/kg/min. Finally, a wide range of inter-medication interaction parameters (i.e., ?1 and ?2 ) was considered to simulate diverse inter-medication interaction: 0.25 ? ?1, ?2 ? 2.75. During the in-silico simulation, nominal CO and RR before the initiation of medication infusion were set at 3.0 lpm and 25 bpm. The initial process parameters ????, ? = 1,2 were set at the average values from the respective parameter ranges selected above. The adaptation gain matrix of ?? = ?4?4 was chosen by trial and error so that the speed of adaptation is maximized while oscillatory adaptation is prevented in all the in-silico subjects investigated in this study. A sampling rate of 1 Hz was used to simulate the process and control actions. The upper bounds of propofol (?1?) and remifentanil (?2?) infusion rates were set to 0.8 mg/kg/min and 0.16 mcg/kg/min, respectively. But, considering that higher infusion rates (e.g., bolus infusion) are typically used for 58 reference target tracking during the initial transients, higher propofol (4 mg/kg/min) and remifentanil (0.8 mcg/kg/min) infusion rates were allowed during the first 10 min after the control action started. To investigate the efficacy and limitations of the semi-adaptive switching control approach, two in-silico evaluations were performed. First, to examine the role and benefit of switching control relative to non-switching control (i.e., the remifentanil-dominant control in Eq. (47) and the propofol-dominant control in Eq. (52)), a series of in-silico simulations were performed in which the semi-adaptive switching and non-switching controllers guided in-silico subjects with diverse dose-response profiles to a target clinical set point. For this purpose, 50 in-silico subjects were generated through a uniformly random selection of parameters in the dose-response models in Eq. (26)-(27) from the respective range described earlier in this section, and the controllers performed the task of guiding the in-silico subjects from the nominal CO and RR to the CO (?1) and RR (?2) set points of 2.4 lpm and 15 bpm, respectively. The performance of the semi- adaptive switching versus non-switching controllers was compared based on the command tracking errors. In addition, to test the hypothesis that the most beneficial role of switching control may be to assist fast adaptation of dose-response relationship by providing one-time instantaneous updating of model parameters in case the operating regime and the dose-response model used in the controller are inconsistent, the accuracy of the estimated process parameters (i.e., ????? and ????? , ? = 1,2) associated with the switching versus non-switching controllers were computed and compared to each other. 59 Second, to examine the limitation of switching control relative to non-switching control by quantifying its benefit with respect to the operating regime, a series of in-silico simulations were performed in which the semi-adaptive switching and non-switching controllers guided in-silico subjects to a set of target clinical set points corresponding to different levels of medication dominance. For this purpose, 40 in-silico subjects were generated through a uniformly random selection of parameters in the dose-response models in Eq. (26)-(27) from the respective range, and the controllers performed the task of guiding the in-silico subjects from the nominal CO and RR to a set of CO (?1) and RR (?2) set points corresponding to a range of medication dominance levels (in terms of ?1+?2). Specifically, the RR set point was fixed at 15 bpm. In 20 in-silico subjects, the 2 semi-adaptive switching and non-switching controllers were tasked to achieve CO targets ? +? corresponding to 1 2 values between 0.6 and 1 in addition to the RR set point of 15 2 bpm. Then, the performance of the switching versus non-switching controllers was ? +? comparatively examined with respect to 1 2 based on the command tracking errors. 2 Considering that the operating regime was propofol-dominant and the control law in Eq. (52) must be used, of particular interest was to compare the switching control starting from remifentanil-dominant mode to the remifentanil-dominant non-switching control in Eq. (47). In the remaining 20 in-silico subjects, the semi-adaptive switching and non- ? +? switching controllers were tasked to achieve CO targets corresponding to 1 2 values 2 between 0 and 0.4 in addition to the RR set point of 15 bpm. Then, the performance of the switching versus non-switching controllers was likewise examined based on the command tracking errors. Considering that the operating regime was remifentanil- 60 dominant and the control law in Eq. (47) must be used, of particular interest was to compare the switching control starting from propofol-dominant mode to the propofol- dominant non-switching control in Eq. (52). 5.4. Results and Discussion Within the 50 in-silico subjects, the target CO and RR set points of 2.4 lpm and 15 bpm resulted in 28 and 22 simulated subjects operating in the propofol- and remifentanil- dominant regimes, respectively. Therefore, the in-silico subjects generated in this study were adequate in evaluating the semi-adaptive switching control in a wide range of dose- response profiles. Overall, the switching control offered added benefit relative to non-switching control. Table 4 shows the root-mean-squared normalized errors (RMSNE?s; in terms of mean and standard deviation) in command tracking associated with the propofol-dominant non- switching, remifentanil-dominant non-switching, and switching controllers. RMSNE was 1 ? ?? 1 ? ?? computed as ? 1 ?1? + ? 2 ?2? , 0 ? ? ? 15 min, where ??1 and ??? ? ?? ? ?2 are ?1 2 ?2 2 the reference trajectories for ?1 (CO) and ?2 (RR) that can be obtained from Eq. (28) and Eq. (46). On the average, the switching control reduced the RMSNE by 22% and 26% compared to the propofol-dominant and remifentanil-dominant non-switching control if it started from the former, and by 23% and 27% if it started from the latter. It appeared that the benefit offered by the switching control was related to the correct and prompt mode switching, which assisted the controller with fast adaptation to the in-silico subject?s correct dose-response relationship. Table 4 shows the time at which mode switching 61 occurred in the switching controllers. Within the 50 in-silico simulations, there was no incidence of incorrect mode switching, and mode switching, when required, occurred very quickly on the average (within 0.92 min after enforcing the set point). Fig. 14 shows 2 representative in-silico simulation examples. In Fig. 14(a), the in-silico subject was subject to propofol-dominant regime at the set point. Therefore, the propofol- dominant non-switching controller performed the best (RMSNE = 0.009). Not surprisingly, the switching controller starting from propofol-dominant mode showed the same performance to the propofol-dominant non-switching controller without any inappropriate mode switching (not shown). As expected, the remifentannil-dominant non- switching controller performed poorly in comparison to its propofol-dominant counterpart (RMSNE = 0.032). In contrast, the switching controller starting from remifentanil-dominant mode could promptly switch to the propofol-dominant mode (at 0.48 min) to minimize the degradation in performance due to incorrect control mode (RMSNE = 0.011). Similarly, in Fig. 14(b) the in-silico subject was subject to remifentanil-dominant regime at the set point. Therefore, the remifentanil-dominant non- switching controller (RMSNE = 0.008) exhibited superior performance to its propofol- dominant counterpart (RMSNE = 0.032). In this case as well, the switching controller starting from propofol-dominant mode could promptly switch to the remifentanil- dominant mode (at 1.45 min) to minimize the degradation in performance (RMSNE = 0.025). 62 In sum, the above in-silico results suggest that the semi-adaptive switching control offers added benefit to its non-switching counterpart by virtue of its ability to operate in a correct mode best suited to the given operating regime. Table 4: Root-mean-squared normalized command tracking errors (RMSNE?s) associated with switching and non-switching controllers, and switching times associated with switching controllers (mean (SD)). Control RMSNE Switching Time [min] Propofol-Dominant Non-Switching Control 0.019 (0.012) N/A Remifentanil-Dominant Non-Switching Control 0.020 (0.017) N/A Switching Control, Stating from Propofol-Dominant 0.015 (0.009) 0.77 (0.55) Mode Switching Control, Stating from Remifentanil- 0.014 (0.012) 1.03 (1.63) Dominant Mode The results also suggested that the added value offered by the switching control may be due to its superior parameter estimation accuracy to the non-switching control. Fig. 15 shows the parameter estimation errors associated with the semi-adaptive switching versus non-switching controllers. The parameter estimation error was computed as the sum of 8 element-wise RMSNE?s between the actual (? computed from Eq. (29)) versus estimated (????(?) , ? = 1,2, 0 ? ? ? 15 min) controller parameters. In 28 in-silico subjects who resided in the propofol-dominant regime at the target set point, the parameter estimation RMSNE?s associated with the semi-adaptive switching control and remifentanil- dominant non-switching control were 13.5 and 24.2 on the average, respectively (Fig. 63 15(a)). Likewise, in 22 in-silico subjects who resided in the remifentanil-dominant regime at the target set point, the parameter estimation RMSNE?s associated with the semi-adaptive switching control and propofol-dominant non-switching control were 16.0 and 25.6 on the average, respectively (Fig. 15(b)). Thus, switching control could reduce the parameter estimation errors by a large amount. Importantly, the semi-adaptive control exhibited smaller parameter estimation errors than non-switching control in most in-silico subjects used in this study: Fig. 15, which presents the cumulative distribution of the parameter estimation errors, shows that switching control persistently outperforms its non-switching counterpart in terms of parameter estimation accuracy. Hence, it can be deduced that the instantaneous parameter updating offered by the switching control may improve the command tracking performance by improving the model adaptation accuracy. 64 Fig. 14: Representative in-silico simulation examples. Blue vertical line indicates the instant at which mode switching occurred. (a) An in-silico subject in the propofol- dominant regime at the set point. (b) An in-silico subject in the remifentanil-dominant regime at the set point. M1: propofol. M2: remifentanil. 65 Fig. 15: Parameter Estimation RMSNE?s associated with the semi-adaptive switching versus non-switching control. (a) 28 in-silico subjects in the propofol-dominant regime at the set point. (b) 22 in-silico subjects in the remifentanil-dominant regime at the set point. The above results show that the switching control consistently outperformed non- switching control via its ability to provide instantaneous model parameter updating. Yet, its benefit may not be large if the dose-response relationship can be readily tuned to the operating regime with the parameter adaptation alone. In Table 5, the command tracking performance of semi-adaptive switching versus non-switching controllers are compared with respect to the operating regime. Table 5(a) shows that, in subjects who were subject to propofol-dominant regime at the target set point, the difference in the command tracking performance between the switching control (which started from remifentanil- dominant mode) versus remifentanil-dominant non-switching control increased as the ? +? operating regime became more propofol-dominant (i.e., as 1 2 became larger). 2 Likewise, Table 5(b) shows that the same trend was observed for the switching control (which started from propofol-dominant mode) versus propofol-dominant non-switching 66 control in subjects who were subject to remifentanil-dominant regime at the target set point. Hence, the effect of switching was large when the gap between the dose-response relationship estimated by the controller versus the subject?s actual dose-response relationship was large (so that instantaneous model parameter updating could expedite the fast adaptation of dose-response relationship), but small otherwise (in which the gap may have been easily handled by the adaptation law alone). Indeed, Fig. 16 shows that the degree of improvement in command tracking performance (in terms of the RMSNE) ranged from 6% to nearly 35% depending on the operating regime of the subject. This suggests that it may be prudent to selectively activate mode switching (i.e., activate it only in case it can improve the controller?s command tracking performance) to maximize its efficacy. Fig. 16: Mean percent improvement in command tracking performance offered by switching control with respect to operating regime. (a) In-silico subjects subject to propofol-dominant regime. (b) In-silico subjects subject to remifentanil-dominant regime. 67 Table 5: Command tracking performance of semi-adaptive switching and non-switching controllers with respect to operating regime (mean (SD)). (a) In-silico subjects subject to propofol-dominant regime (?1 + ?2)?2 0.6 0.7 0.8 0.9 1.0 0.016 0.017 0.019 0.026 0.065 Switching (0.005) (0.006) (0.008) (0.014) (0.044) Remifentanil-Dominant Non- 0.017 0.021 0.026 0.038 0.094 Switching (0.005) (0.008) (0.013) (0.022) (0.062) (b) In-silico subjects subject to remifentanil-dominant regime (?1 + ?2)?2 0.0 0.1 0.2 0.3 0.4 0.015 0.013 0.013 0.015 0.017 Switching (0.009) (0.007) (0.007) (0.008) (0.009) Propofol-Dominant Non- 0.023 0.020 0.019 0.018 0.020 Switching (0.014) (0.011) (0.010) (0.011) (0.013) 68 Chapter 6: Coordinated Two-Input Two-Output Semi-Adaptive Control Design In this chapter, the problem if the initial set targets are not achievable is addressed. A coordinated two-input two-output semi-adaptive controller which consists of a coordinated control and the controller in chapter 4 is introduced. In the following, the architecture of the controller is discussed first. Then the coordinated control is presented. Comparative simulation studies were also done to demonstrate the advantages of the proposed method over the method in chapter 4. Moreover, the detailed behavior of the coordinated controller was examined. 6.1. Control Architecture The proposed control approach consists of an upper-level coordination controller and a low-level semi-adaptive controller (Fig. 17). The coordination controller recursively adjusts the reference targets, based on the dose-response relationship of a subject estimated by the semi-adaptive controller and constraints imposed on the medication use, in order to ensure that the reference targets thus derived can be safely achieved by the subject. For given reference targets provided by the coordination controller, the semi- adaptive controller in chapter 4 is utilized to compute and execute the requisite medication infusion rates to guide the subject towards the reference targets while estimating the subject?s dose-response relationship on-line and providing it to the coordination controller. 69 Figure 17: The architecture of coordinated semi-adaptive control. 6.2. Coordinated Control via Recursive Reference Adjustment The primary role of the coordinated control is to keep the reference targets achievable by the system driven by the semi-adaptive control outlined in Chapter 4. A coordinated control scheme based on a recursive adjustment law for the reference targets was developed by considering constraints on (1) the input magnitudes (i.e., bounds on the medication infusion rates) and energy (i.e., total medication use) as well as (2) the degree of discrepancy between the original (i.e., specified by the caregiver) versus adjusted reference targets. The set of achievable reference targets can be determined by the bounds imposed on the input as well as the discrepancy between the originally specified versus adjusted reference targets as follows. First, in case of medication infusion, the elements of the input ? must be positive while limited by an upper bound to ensure subject safety: 0 ? ?? ? ????, ? = 1,2. Hence, considering (34), the reference targets must satisfy the following inequalities in the steady state: 70 ?1 ?1 ? 0 ? ? = [? ] = ? ??? [? ] ? [ 1?] (58) 2 2 ?2? Referring to (33), (58) yields the following pair of inequalities in terms of ?1 and ?2: ?11?12?21 ?11?21?22 ?22 ?22 ??1? 0 ? ? ? 1 ? ?2 ? ?1? ? ?? ? 2 ? ?1 ? ?2 + 12 ?12 ?12?11?21 (59) ?12?21?22 ?11?12?22 ?21 ??2? ?21 0 ? ?2 ? ?1 ? ?2? ? ? ? ? ?? ? ? 2 ? ? ? 1 ? ? ? 2 11 11 12 22 11 which altogether yields a parallelogram in the (?1, ?2) space as shown in Fig. 14 (note that the parallelogram is guaranteed to exist, because the slopes of the constraints in (59) ? ? satisfy 21 > 22 given that all the leading principal minors of ? are positive). Second, the ?11 ?12 discrepancy between the original (i.e., caregiver-specified) versus adjusted reference targets may be limited to respect the expertise of the caregiver, which results in the following inequalities in terms of ?1 and ?2 (note that these limits may be specified by the caregiver in advance): ??0 ? ?? ? ?? ? ??0 + ?? , ? = 1,2 (60) ? ? where ??0 is the originally specified value for ?? , ? = 1,2. These constraints altogether yield a rectangle in the (?1, ?2) space as shown in Fig. 18. Finally, the set of achievable reference targets is determined as the intersection between the parallelogram and rectangle, as shown in Fig. 18. If the current pair of reference targets (?1, ?2) is not achievable, one (and perhaps the easiest) way to make it achievable is to adjust it towards the nearest achievable pair of reference targets. From Fig. 18, the nearest achievable pair can be found as the intersection between the side of the parallelogram having the smallest distance from the 71 r1 ?22 ??1? 2?? ?1 = ? + 2 ? 212 ?12?11?21 A D Achievable Targets 2?? 1 B ?21?1 = ?2 ?11 C Optimal ?21 ??2? Regime ?1 = ?2 ? ?11 ?11?12?22 ?22 ?1 = ? ? 212 r2 Figure 18: The set of achievable reference targets specified by the constraints on the input magnitudes and deviation from the original reference targets. current pair of reference targets and the line perpendicular to the side that passes through the current pair of reference targets (e.g., the pair ?A? is adjusted towards ?B?). However, this approach is not ideal in terms of input energy (i.e., medication use), because the adjusted reference targets will be on a side of the parallelogram, at which ?? = 0 or ?? = ????, ? = 1,2 and the synergistic interaction between the two medications is not exploited to minimize the medication use. Hence, it may be sensible to penalize the medication use in determining the direction of adjustment of the reference targets in order to minimize the medication use. The optimal direction of adjustment, obtained by considering both the smallest distance to the set of achievable reference targets and input energy, may not be strictly perpendicular to the side of the parallelogram. Regardless, considering that all the sides of the parallelogram have positive slopes, that ?1 and ?2 must be adjusted in the opposite direction (i.e., if one is increased, the other must be decreased) is a viable 72 requirement in adjusting the reference targets in order to effectively reach the achievable reference targets (e.g., adjusting the pair ?A? towards ?C? is more effective than adjusting towards ?D?). Hence, this can be enforced in the course of reference targets adjustment. Summarizing the above considerations, we developed a coordinated control scheme based on a recursive adjustment law for the reference targets using the model predictive control formalism as follows. Consider a simple recursive adjustment law for ?1 and ?2 given by: ? (?) ? (?) ???(?) = [ 1? ] = ? [ 1? ] = ???(?) (61) ?2?(?) ?2(?) ?? 0 where ?? = [ 1 ] is a positive definite adaptation gain matrix and ?(?) the 0 ??2 adjustment policy that is to be designed. Then, the dynamics of the plant (35), reference model (36), MRAC law (37), and reference target adjustment law (61) combined all together is given in the discrete-time domain by: ?(? + 1) = ?(?) + ??{(???(?) + ???(?)????(?) + ???(?)?)?(?) ? ???(?)???(?) + ???(?)????(?)?(?)} ? ? ? ??? ? ?(?)??? (? + 1) = [ ] = ??? (?) ? ?? {?? [ ] [?(?) ? ??(?)] ?} (62) ? ? ?(?)??? ??(? + 1) = ??(?) + ??{????(?) + ???(?)} ?(? + 1) = ?(?) + ?????(?) ?? ?? with ???(?) = ?????? and ???(?) = ?? ? ?????? ???? . To reconcile constraints on the input magnitudes (59) and energy (i.e., total medication use) as well as the degree of 73 discrepancy between the original versus adjusted reference target (60), ?(?), ? ? ? ? ? + ?? is determined by minimizing the following cost function, where ?? and ?? denote the control and prediction horizons, respectively: ?+?? 2 2 ?(?) = ? [?1(?1(?) ? ?10(?)) + ? 2 2(?2(?) ? ?20(?)) + ?3?1(?) ?=? (63) + ?4? 2 2(?)] subject to the linear inequality constraints (59) and (60), and the following nonlinear inequality constraints to ensure that ?1 and ?2 are adjusted in the opposite direction: ?2(?) ?1(?)[?2(?) ? ?????1(?)] ? 0???? ? ? ???? ? { ,?1(?) ?1(?)[?2(?) ? ?????1(?)] ? 0 (64) ? ? ? ? ? + ?? where ???? < ???? < 0. At each sampling time step, the coordinated control problem is solved to yield ?(?), ? ? ? ? ? + ?? that minimizes (63) subject to the dynamics (62) with the dose-response model parameters ???? and ???? adapted in the previous sampling time step as well as the constraints (59), (60), and (64). Then, ?(?) is applied to the system to adjust the reference targets, which are tracked by the semi-adaptive controller while the subject?s dose- response model parameters are adapted. This process is repeated at each sampling time step to guide the subject outputs to the (recursively adjusted) reference targets. To prevent the drift in ?1 and ?2, we employed a dead-zone scheme to the coordinated controller, so that the adjustment of the reference targets are made only when ?1 and ?2 are sufficiently large: 74 ??(?) + ???? ??(?), |?1(?)| > ?? or |?2(?)| > ?? ? (? + 1) = { ? 1 2? (65) ??(?), otherwise which, together with the dead-zone scheme employed in the adaptation law (41), helps to avoid unnecessary drift in the reference targets. 6.3. In-Silico Implementation and Simulation To examine the performance of the coordinated semi-adaptive controller, we utilized the same settings in Chapter 4.3. 20 in-silico subjects which cannot achieve the initial targets are selected for simulation. The weights in the cost function (63) were selected based on the in-silico simulation result in an average in-silico subject, so that (1) the cost associated with the discrepancy between the original versus adjusted reference targets is 1,000 times larger than the cost associated with the total medication use; (2) the costs associated with the two reference targets are equal; and (3) the costs associated with the two medications are equal. In summary, the parameters associated with the coordinated semi-adaptive control were 0.8 0 5 0 as follows: ?? = [ ] = ??? , ? = [ ], ?? = ?4?4 , ?? = 0.01, ?? = 0.01, 0 0.8 0 5 1 2 60 0 ?? = [ ] , ?? = 2 , ?? = 20 , ?1 = 1000 , ?2 = 825.5 , ?3 = 1 , ?0 60 4 = 40 , ?1? = 0.8mg/kg/min, ?2? = 0.07 mcg/kg/min, ?? = 0.3, ?? = 0.3, ???? = ?2, ???? = ?0.5, 1 2 ?? = 0.02, ?? = 0.02. 1 2 75 2 key aspects related to the performance of the coordinated semi-adaptive control were tested. First, we examined the performance of coordinated semi-adaptive control compared with the same semi-adaptive control without reference target adjustment. For this purpose, we conducted in-silico evaluation of coordinated versus non-coordinated semi-adaptive controllers using the 20 in-silico subjects with unachievable reference target. Second, we examined the detailed behavior of the coordinated controller, especially the way it reconciles the constraints on (1) the infusion rate limits and total medication use versus (2) the degree of discrepancy between the original versus adjusted reference targets, as the values of the weights in the cost function (63) vary. Fig. 19 presents a representative example of the evolution of adjusted reference targets and system outputs associated with the coordinated semi-adaptive control in an in-silico subject with unachievable reference targets. The result demonstrates that the initially specified unachievable reference targets are adjusted towards the optimal regime (shown in Fig. 18), and the semi-adaptive control guides the in-silico subject to these new achievable reference targets as time evolves. In contrast, the system outputs converges to a point on the parallelogram in case the coordination control is not employed, meaning at least one input is saturated with large steady-state error(s). This desired behavior was consistently observed in all in-silico subjects. 76 Figure 19: The evolution of adjusted reference targets and system outputs associated with the coordinated semi-adaptive control in a in-silico subject with unachievable reference targets. (a) Time plots. (b) Phase plot. 77 Overall, the coordinated controller could largely reduce the overall reference tracking errors associated with cardiac output (60%) and respiratory rate (96%) over the entire in- silico simulation ( 0 ? ? ? 30 min) as well as steady-state error associated with respiratory rate (99%; in the absolute sense), compared with when the coordinated controller was not employed to adjust the reference targets (Fig. 20). In addition, the coordinated controller could also decrease propofol use while increasing remifentanil use, by adjusting the reference targets so that inter-medication synergy can be exploited more effectively (Fig. 16). The steady-state error associated with cardiac output was relatively large with the coordinated controller. This was due to the dead zone implemented for the parameter adaptation law in (41). In the absence of coordinated controller, the adaptation law was persistently enabled to drive cardiac output error to zero, because the tracking error associated with respiratory rate was large (since its reference target could not be reached) and the dead zone was never activated. On the other hand, the coordinated controller activated the dead zone by adjusting the reference targets, thereby allowing the tracking errors to become small. Noting that the dead zone for cardiac output in (41) was 0.03 lpm and that the steady-state cardiac output error was consistently smaller than 0.03 lpm, these errors were deemed appropriate. It is worthwhile to scrutinize the behavior of the system inputs and outputs pertinent to the particular in-silico simulation conducted in this section to glean more insights on the coordinated controller. From our in-silico simulation, a large decrease in the respiratory rate error and a large increase in the remifentanil use were consistently observed. 78 However, the changes associated with the cardiac output error and propofol were relatively small. We speculate that these observations may be interpreted as follows. First, the dose-response model parameters we employed in the in-silico simulation dictate that cardiac output is primarily influenced by propofol while respiratory rate is influenced by both propofol and remifentanil. Second, the propofol infusion rate required to achieve the original cardiac output target (?1(0)=2.0 lpm) in the steady state leads to steady-state respiratory rate lower than the original respiratory rate target (?2(0)=15 bpm) on the average. Hence, the coordinated controller increases the reference target for cardiac output (as realized by a decrease in ?1; Fig. 19) and decreases the reference target for respiratory rate (as realized by an increase in ?2; Fig. 19) in order to decrease the required propofol infusion rate and increase the required remifentanil infusion rate while exploiting the synergy between the two. In this way, the coordinated controller fulfills the desired objectives of driving unachievable reference targets to achievable ones, and at the same time, minimizing the total medication use. It is noted that the proposed control approach exhibited consistent performance regardless of the location of the original reference targets in our in-silico simulation testing. For the dependence of coordinated semi-adaptive control on cost function weights, the coordinated controller exhibited intuitively relevant behaviors as the weights in the cost function (63) were varied on the average, as described below (Fig. 20). 79 0.1 12 10 20 8 15 0.05 6 10 4 2 5 0 0.5 1 1.5 0 0Coordinated Non-Coo2rdinate2d.5 0.5Coord1inated1.N5on-Coo2rdinat2e.d5 0.5 Coord1inated 1N.5on-Coo2rdinate2d.5 0.1 12 0.15 0.08 10 8 0.06 0.1 6 0.04 4 0.05 0.02 2 0 0 0 0.5 Coord1inated 1N.5on-Coo2rdinate2d.5 0.5Coord1inated1N.5on-Coo2rdinat2e.d5 0.5Coord1inated 1N.5on-Coo2rdinate2d.5 Figure 20: Distribution of overall reference tracking and steady-state errors as well as total medication use associated with semi-adaptive control in the presence/absence of coordination control. (a) {r3,r4}=k?(r30,r40} (b) r1=k?r10 (c) r3=k?r30 k k k Figure 21: Dependence of the behavior of coordinated semi-adaptive control on cost function weights. The left and right axes apply to squares and circles, respectively. 80 Propofol Use [mg/kg] |r1-r1(0)| [lpm] Steady-State Error: CO [lpm] RMSE: CO [lpm] Remifentanil Use [mcg/kg] |r2-r2(0)| [bpm] Steady-State Error: RR [bpm] RMSE: RR [bpm] Propofol Use [mg/kg] |r1-r1(0)| [lpm] Remifentanil Use [mcg/kg] |r2-r2(0)| [bpm] Remifentanil Use [mcg/kg] Propofol Use [mg/kg] Propofol Use [mg/kg] |r1-r1(0)| [lpm] Remifentanil Use [mcg/kg] |r2-r2(0)| [bpm] First, when the penalty on the control energy terms is increased relative to the discrepancy between the original versus adjusted reference targets (by simultaneously altering ?3 and ?4 ; Fig. 21(a)), the coordinated controller decreases propofol and remifentanil uses. This leads to a decrease in ?1 and ?2 via an increase in the reference targets for cardiac output and respiratory rate (note that ?1 and ?2 are specified in terms of ?1 and ?2). Then, noting that ?1 < ?1(0) and ?2 > ?2(0) under the nominal weights in our in-silico simulation testing (see ?A? and ?B? in Fig. 5), |?1 ? ?1(0)| increases while |?2 ? ?2(0)| decreases. Second, when the penalty on the discrepancy between the original versus adjusted cardiac output reference target is increased relative to other penalties (by altering ?1; Fig. 21(b)), the coordinated controller decreases |?1 ? ?1(0)| . To still drive the reference targets towards the optimal regime (Fig. 3) while limiting |?1 ? ?1(0)|, it increases |?2 ? ?2(0)|. Since ?1 < ?1(0) and ?2 > ?2(0) under the nominal weights in our in-silico simulation testing, this leads to a decrease in cardiac output and respiratory rate reference targets, and accordingly, an increase in propofol and remifentanil uses. The opposite behavior is observed when ?2 is altered instead (not shown). Third, when the penalty on the propofol use is increased relative to the other penalties (by altering ?3; Fig. 21(c)), the coordinated controller decreases propofol use while increases remifentanil use. This leads to a decrease in ?1 , and since ?1 < ?1(0), an increase in |?1 ? ?1(0)|. In addition, a decrease propofol use, which is relatively large compared with the increase in remifentanil use, results in the corresponding decrease in ?2 and 81 (since ?2 > ?2(0)) |?2 ? ?2(0)|. The opposite behavior is observed when ?4 is altered instead, though the extent is relatively small (not shown). The abovementioned intuitive behaviors of the proposed coordinated semi-adaptive control approach may be useful in tuning the weights in the cost function (63) to customize its performance when applied to other medication infusion problems. 82 Chapter 7: Pilot Animal Experiments and Further Improvements for the Controller Design To verify the proposed controller, the experiment of SISO controller in Chapter 3 was performed. However, the performance was not as good as the simulation results. After carefully examination of the experiment results, we observed that transport delay existed in the dose-response relationship, which was not considered in the controller design and led to the deterioration of controller performance. Therefore, a semi-adaptive controller which incorporates transport delay is developed in this chapter. In this chapter, the experiment results are first presented and discussed. Then a novel dynamic dose-response model for control design is presented. Moreover, parametric sensitivity analysis of the dynamic dose-response model is performed to assess the influence of transport delay on the model output in comparison with the other dose- response model parameters. Then the indirect APPC design for medication infusion problems based on the Pad? approximation of transport delay is developed. In the end, the results are presented and discussed. 7.1. Pilot Animal Experiment To examine the performance of the semi-adaptive controller, we first performed experiments on 5 pigs. Because estimation of controller parameters were very large for some pigs, we augmented parametric projection and dead zone in the controller to avoid this. Based on input-output data of 5 pigs, we performed system identification and created in-silico pigs of them. Then we simulated our proposed controllers on the 5 in- 83 silico pigs with different controller gains. Based on the best controller gains from simulations, we performed the semi-adaptive control experiments in another 4 pigs. In the experiments, CO was derived from the aortic flow waveform measurement (14 PAU and T402, Transonic Systems, Ithaca, NY) by first low-pass filtering the aortic flow waveform and then computing its average over 5-second window. Computed CO was then transmitted to a laptop equipped with the semi-adaptive controller with parametric projection and dead zone algorithms through an A/D converter (NI-USB6002, National Instruments, Austin, TX). The infusion rate computed by the controller was transmitted to a digital infusion pump (Fusion 100, Chemyx, Stafford, TX) via serial communication. Then, the pump delivered the commanded infusion rate to the pig. A sampling rate of 100 Hz was used in all the computations, while the control command was updated every 5 seconds. The semi-adaptive controller performed well in 3 pigs while it performed marginally in 1 pig. Our retrospective analysis indicated that dynamic dose-response delay, which is not modeled explicitly in the controller design process, may play an important role in determining the performance of the controller in each pig. Indeed, the semi-adaptive controller showed good performance in pigs with small dose-response delay, but it suffers from limited performance in pigs with large dose-response delay (Fig. 22). The results illustrate that incorporating dynamic dose-response delay into controller design must be pursued. 84 Fig. 22: Experimental results. (a) 3 pigs with small dose-response delay. (b) 1 pig with large dose-response delay. 7.2. Control Design Model: Direct Dynamic Dose-Response Model To design the controller for medication infusion problems with non-negligible transport delay, a novel dynamic dose-response model is developed. The model consists of a low- order lumped parameter latency model (66a) and a modified Hill equation model (66b): ????(?) = ?????(?) + ???(? ? ?) ? ?[??(?)] (66a) 100 ? ( ? 1) ?? (?) ?(?) = ?? (1 ? ? ) = ?[??(?)] (66b) ? 100 ? ?? + ( ? 1) ?? (?)? where ?(?) is the intravenous medication infusion rate, ??(?) is the medication infusion rate at the site of action, ? is the transport delay between medication infusion and the onset of the intended clinical effect, ?? is the rate constant associated with the distribution of medication from intravenous site to the site of action, ?(?) and ?? are clinical effect and its nominal (i.e., in the absence of medication infusion) value, ? is the 85 target percentage depression of ?(?) from ?? (meaning that target clinical effect ?(?) is 100?? specified by ?(?) = ??), ?? is the infusion rate required to maintain ?(?) = ?(?) in 100 the steady state, and ? is the cooperativity constant. There are a few advantages associated with this dose-response model compared with conventional PKPD models. First, intended to primarily capture the direct relationship between dose and the resulting clinical effects, its dynamic component is much simpler than multi-compartmental PKPD models. Second, it explicitly incorporates transport delay which can be exploited to represent PD and/or subject monitor delays. Third and most importantly, it facilitates sensitivity-based semi-adaptive control design [55], in which model parameters exerting high sensitivity on the model?s behavior are adapted while those with low sensitivity are fixed at nominal values (see section 7.3 for details). 7.3. Parametric Sensitivity Analysis In this section, analytical and numerical parametric sensitivity analyses are performed to demonstrate that (1) the proposed dose-response model allows us to minimize the influence of ? with appropriate choice of ?; and (2) the transport delay L needs to be adapted for control efficacy. Details follow. The proposed dose-response model (1) de-sensitizes the cooperativity constant ? in the steady state when ?(?) = ?(?), or equivalently, when ??(?) = ??, making it possible to fix it at a nominal value in the control design process without compromising controller robustness against its uncertainty. To illustrate, regard (1a) and (1b) as the state and 86 output equations for the process dynamics. Then, the following sensitivity functions can be derived: ?? ?? ?? ?? ?? ? ?? (?) = ?? (?) + [ ] ? ?? (?) ?? ??? ?? ??? ?? ??(? ? ?) = ????? (?) + [???(?) + ?(? ? ?) 0 0 ?? ] ? ? ?(? ? ?) ?? ?? ?? ?? ?? (67) ??(?) = ?? (?) + [ ] ???(?) ? ??? ?? ??? ?? 100 ? ??1 ? ( ? 1) ???? (?)? ?? ?? ?? ??= ??? ?? (?)2 + [ ] ? 100 ??? ?? ??? ? ? ?? (?? + ( ? ? 1) ?? (?)) ?? (?) ?? (?) ?? (?) ?? (?) ??(?) ??(?) ??(?) ??(?) where ?? (?) ? [ ? ? ? ? ] , ??(?) ? [ ] , and ? ??? ?? ??? ?? ??? ?? ??? ?? 100 ? ? 100 ??1 ? ?H ?H ?H ?H ( ?1)I?Ie(t)(log(I?)?log(Ie(t))) ?( ?1)I? Ie(t)? ? [ ] = [0 yo 2 yo 2 0] . ?ke ?? ?I? ?L ? 100 ? ? 100 ?(I?+( ?1)Ie(t)) (I?+( ?1)I (t))? ? e Then, the closed-form formula for ??(?) is given by: 100 ??1 ? ? ( ( ) ? ?? ? ? 1) ?? ?? ? ? 2 ? ? ????(???(?) + ?(? ? ?))?? ? 100 ? 0(?? + ( ? 1) ?? (?)) ? 100 ? ? ( ? ? 1) ???? (?)(???(??) ? ???(? ? (?))) ? ? 2 ? 100 ? (?? + ( ? ? 1) ?? (?)) ??(?) = ? (68) 100 ? ??1? ( ? ? 1) ???? (?) ?? 2 ? 100 ? (?? + ( ? ? 1) ?? (?)) 100 ??1 ? ? ( ? 1) ? ?? ?? (?) ??(? ? ?) ?? ? ? ?????? 2 (??? )?? ? 100 ?(? ? ?) ? 0 (?? + ( ? 1) ?? (?)) [ ? ] 87 ??(?) Hence, , which is the second element of ??(?) in (68), becomes zero when ??(?) =?? ??. The implication of the zero sensitivity of ?(?) on ? for ??(?) = ?? is that the influence of ? on the process behavior can be minimized by setting ? according to the target clinical effect, that is, equal to the percentage difference between ?? and ?(?), i.e., ? = ????(?) ? 100. Indeed, with this choice of ?, it can be readily shown that the process ?? behavior becomes insensitive to ? under steady-state target tracking condition (?(?) = ?(?)). To confirm the above analytical sensitivity analysis as well as to assess the significance of transport delay ? relative to the other parameters in the dose-response model, a numerical parametric perturbation analysis was performed using the regulation of a cardiovascular variable cardiac output (CO) with a sedative propofol as a case scenario. In this analysis, a fine-tuned population-based PID controller was used to regulate CO via propofol infusion in 30 randomly created in-silico subjects for nominal closed-loop response. Then, the dose-response model parameters (including ?? , ?? , ?, and ?) were perturbed, one at a time, by ?25% and ?50% (thus, 4 perturbations per each parameter), and the perturbed in-silico subjects were simulated with the same PID controller for perturbed closed-loop responses. Then, the difference between nominal and perturbed responses were quantified in terms of root-mean-squared errors (RMSEs) between the two up to the settling time of the desired response. Finally, a total of 120 RMSEs associated with each parameter was aggregated to compute mean and standard deviation (SD). We used the following ranges of the dose-response model parameters obtained from our prior work [56] and our unpublished in-house experimental data to create 30 in- 88 Fig. 23: Root-mean-squared errors (RMSEs) between nominal versus perturbed clinical effect responses to perturbation of dose-response model parameters in 30 in-silico subjects in a case study of regulating a cardiovascular variable (cardiac output (CO)) with a sedative (propofol). 89 silico subjects: 0.02 ? ?? ? 0.10 min-1, 0.2 ? ?? ? 0.6 mg/kg/min, 1 ? ? ? 5 , and 50 ? ? ? 100 s. In all in-silico simulations, ?0 was set at 3.0 lpm. A multitude of desired responses, in terms of both magnitude and rate, was employed in the analysis: 2.1, 2.4, 2.7 lpm with the time constants of 5, 2.5, and 1.7 min. For each of these 9 desired responses, the aggregated RMSEs associated with the dose-response model parameters were compared to assess if the perturbation in the transport delay makes a large impact on the closed-loop control performance relative to the other dose-response model parameters, in order to determine if the adaptation of transport delay is necessary. Results from the parametric sensitivity analysis indicated that the impact of transport delay on the closed-loop clinical effect response (as measured by the RMSE between nominal and perturbed responses up to the settling time of the desired response) was larger than cooperativity constant but not as large as ?? and ??. Figure 23 presents the RMSEs between nominal and perturbed responses associated with each model parameter, which offers several key observations. First, all in all ?? exerted the largest impact on the model?s clinical effect response, followed by ?? , and then ? , and finally ? . Second, although ? was not the most crucial parameter in (1) in terms of average sensitivity, it exhibited a very large variability in sensitivity relative to its average counterpart (i.e., large coefficient of variation) in comparison with the remaining parameters in (1). This is attributed to a large deviation of the closed-loop controlled clinical effect response from its nominal counterpart when L assumes very large values, since a large L reduces the stability margin associated with the population-based PID control (which does not accommodate the increase in ? ). Hence, it deemed reasonable to investigate the 90 advantage of its adaptation by comparing an adaptive control equipped with the capability of adapting ? versus an adaptive control with ? fixed at a nominal value. Third, although the dose-response model (66) by construction exhibits zero sensitivity to ? when ??(?) = ?? if ? is specified in accordance with the target clinical effect (i.e., if ? is set equal to the percentage difference between the baseline and target clinical responses), its sensitivity to ? is not zero during transients. Hence, the finding that ? exerts the smallest influence on the dose-response model?s clinical response even during transients suggests that all in all the model?s sensitivity to ? is the smallest and that ? may indeed be fixed at a nominal value, thereby supporting the validity of the transformed dose-response (6) under the assumption of a predetermined ?. 7.4. Sensitivity-Based Semi-Adaptive Pole Placement Control Design Our sensitivity-based semi-APPC is built upon (1) Pad? approximation of transport delay, (2) linear model parameterization, (3) recursive model parameter adaptation, and (4) its integration into pole placement control. Figure 24 shows the scheme of the sensitivity- based semi-APPC. Details follow. For the sake of adaptive control design, the transport delay in the dose-response model (66) was simplified into the Pad? approximation. It is obvious that the higher the order of approximation, the more accurate the approximation is. But, a high-order approximation complicates control design by increasing the order of process dynamics. Besides, a relatively low-order approximation may still be appropriate for adaptive control, because 91 the adaptation of transport delay parameter may mitigate the approximation error and still produce correct phase delay [57]. Fig. 24: Sensitivity-based semi-adaptive pole placement control (semi-APPC) for closed- loop control of medication infusion. The semi-APPC is built upon (1) Pad? approximation of transport delay, (2) novel linear model parameterization, (3) recursive model parameter adaptation, and (4) its integration into pole placement control. To determine the Pad? approximation relevant to the problem at hand, 6 candidate approximations in (69) were considered: ?? M1: ??(?) = ?(?) ? + ?? ?? 1 M2: ??(?) = ?(?) ? + ?? 1 + ?? ?? 2 ? ?? M3: ??(?) = ?(?) (69) ? + ?? 2 + ?? ?? 6 ? 2?? M4: ??(?) = ?(?) ? + ?? 6 + 4?? + (??)2 ?? 12 ? 6?? + (??) 2 M5: ??(?) = ?(?) ? + ?? 12 + 6?? + (??)2 92 ? 2? 60 ? 24?? + 3(??) M6: ??(?) = ?(?) ? + ?? 60 + 36?? + 9(??)2 + (??)3 Then, the dose-response model (66) with these candidate approximations equipped with nominal parameter values (associated with 30 in-silico subjects) used in the parametric sensitivity analysis in Chapter 7.3 was in-silico simulated with an escalated propofol dose with 5 infusion rate levels designed to elicit a wide range of CO response. The CO responses associated with each of the 6 candidate approximations were then compared with the response of the original dose-response model without Pad? approximation. Specifically, RMSEs between the responses associated with the original and all the Pad?- approximated models were computed, and the optimal order of the Pad? approximation for the problem was determined based on the trend of the RMSE with respect to the order. Table 6 shows the RMSE between the responses associated with the original (i.e., (66)) and all the Pad?-approximated (i.e., (66) with the M1~M6) dose-response models. Clearly, M1 (which does not account for transport delay) suffers from the largest error, whereas the models incorporating the approximation of transport delay can largely reduce the error. The amount of reduction in error becomes trivial beyond M3. Hence, we used M3 for semi-APPC design in this study. 93 Table 6: Root-mean-squared errors (RMSEs) between the responses associated with the original and all the Pad?-approximated dose-response models (M1~M6). M1 M2 M3 M4 M5 M6 RMSE [lpm] 0.093 0.009 0.004 0.002 0.002 0.001 (mean (SD)) (0.039) (0.004) (0.002) (0.001) (0.001) (0.001) 7.5. Dose-Response Model Parametrization For the sake of control design, the dose-response model (66) was parameterized as follows. First, by fixing ? to a nominal value ???, (66b) can be written for ??(?) as follows: ??? ?0 ? ?(?) ??(?) = ?? ? ? ? ?(?) 100 ? (70) ( ( )? ? 1) ? ? ??? ? ??(?) where ?(?) = 0? 100 . Then, (66a) can be written in terms of ?(?) as follows: ( ?1)?(?) ? 1 ???(?) ?? ?? 1 ?? ???(?) = = ? ?(?) + ?(? ? ?) ? ?(?) = ?????(?) (71) ?? ?? ?? ?? ?? (? + ??) Using the Pad? approximation M3 in (69), ?(?) in (71) can be written as follows: 2 ??(?) 1 ?? (? ? ) ?(?) = ?(?) = ? ? ?(?) (72) ??(?) ?? 2(? + ??) (? + ?) For a given the desired clinical effect ?(?) , the reference model was specified by ? ??(?) = ? ?(?) with ?? > 0. Note that ?(? + ??)?? = 0 if ?(?) is constant. ?+?? 94 7.6. Semi-Adaptive Pole Placement Control We aim to design an adaptive controller for (72). Since (72) is non-minimum phase, the application of direct adaptive control techniques (e.g., direct model reference adaptive control) is not trivial [58]. Therefore, we pursued indirect APPC due to its ability to handle non-minimum phase plants [58]. Details follow. The adaptive law for on-line parameter estimation was derived from the standard recursive gradient algorithm [58]. The input-output model (72) can be re-written as ??(?)?(?) = ??(?)?(?), or equivalently: 2 2? ? 2? [?2 ? ? ? + (?? + ) ? + ] ?(?) = (? ? + )?(?) (73) ? ? ?? ??? Rearranging into linear parametric form: 2 2? ? 2? ?2 ? ? ? ?(?) = ?(?? + ) ??(?) ? ?(?) ? ??(?) + ?(?) ? ? ?? ??? ??(?) (74) ? 2? 2 2? ? ? ? ?(?) = [? (?? + ) ] ?? ??? ? ? ???(?) [ ??(?) ] 1 Multiplying a stable low-pass filter with ?0 > 0 to avoid differentiation and (?+? )20 ?2 defining the output as ?(?) = ?(?)2 leads to the following linear parametric model: (?+?0) ?2 ?(?) = ?(?) = ???(?) (75) (? + ? )20 95 ? ?(?)(? + ? )2 0 1 ?(?) ?? 2?? 2 2?? (? + ? )2 = [? (?? + ) ] 0 ? ? ? ? ? ? ? ? ? ?(?)(? + ? 2 0 ) 1? ?(?) [ (? + ? 20) ] ? ?? 2?? 2 2? ? ? where ? ? [?1 ?2 ?3 ?4] = [? (?? + ) ] . Then, the following ?? ??? ? ? adaptive law can be derived from the recursive gradient algorithm [58]: ?(?) ? ??(?)?(?) ???(?) = ??(?)?(?) = ??(?) (76) 1 + ??(?)?(?) To robustify the adaptive law by preventing the drift of the parameter estimates when ?(?) ? 0, (76) was augmented by the following dead zone scheme: ??(?)?(?), |?(?)| > ?0 ???(?) = { (77) 0, |?(?)| ? ?0 Here, we derive standard pole placement control law for the dose-response model (72) incorporating the Pad? approximation, and implement the APPC by combining the pole placement control law thus derived with the adaptive law (76)-(77) by leveraging the certainty equivalence principle. Let ?? and ?? the degrees of ??(?) (?? = 2) and ? ? ? ??(?) (?? = 2), respectively. The standard pole placement control law is given by ? [58]: ??(?)?(?)?(?) = ??(?)[?(?) ? ??(?)] = ??(?)?(?) (78) where ??(?) = ?(? + ??) is the internal model associated with ??(?), ?(?) is a monic polynomial of degree ?? ? 1 , and ?(?) is a polynomial of degree ?? + ? ? 1 . ? ? ?? 96 Hence, ?(?) = ? + ?0 and ?(?) = ?3? 3 + ? ?22 + ?1? + ?0, where ?0, ?3, ?2, ?1, and ?0 are unknown coefficients to be determined via pole placement. Substituting ?(?) in (72) by (78) yields the following closed-loop transfer function between ?(?) and ??(?): ??(?)?(?) ?(?) = ? (?) (79) ?(?)??(?)??(?) + ?(?)??(?) ? Hence, the characteristic equation ?(?)??(?)? th ?(?) + ?(?)??(?) is a 5 -order polynomial. The objective of pole placement control design is to select ?(?) and ?(?) so as to design the characteristic equation: ?(?)??(?)??(?) + ?(?)??(?) = ? ?(?) (80) where ??(?) is a desired 5th-order polynomial. Denoting ??(?) = ?5 + ?4 ? ??=0 ?? ? , the solutions ?(?) and ?(?) to (70) can be found by solving the following algebraic equation: ??? (0 4 ? ?3 + ? ? ) ? ?3 ( ) ?3 ? ?4 + ???3 ?? ?? ? ,? ?2 = ?2 ? ?4?? (81) ? ? ? ?1 ? ? 1 [? ?0] [ ?0 ] where ?? ? ,? is the Sylvester matrix associated with ? 4 ? ? ? ? (?)??(?) = ? + (? + ? )?3 + (? + ? ? )?23 ? 4 ? 3 + ?4??? and ??(?) = ?1? + ?2 given by: 1 (?3 + ??) (?4 + ? ? ?3) ?4?? 0 ?1 ?2 0 0 0 ?? ? ,? = 0 ? ? 0 0 (82) ? ? ? 1 2 0 0 ?1 ?2 0 [ 0 0 0 ?1 ?2] 97 Since ??(?) (especially the value of ??) can be chosen so that ??(?)??(?) and ??(?) are coprime (as long as ? > 0), the Sylvester matrix (72) has full rank. Therefore, the unknown polynomial coefficients of ?(?) and ?(?) in (71) can be determined by: ?0 ? ? 4 ? (?3 + ??) ? ?? ? (? 3 ?1 3 4 + ???3) ? = (?? ) ?? 2 ????,?? 2 ? ?4?? (83) ?1 ? ?1 [?0] [ ??0 ] Hence, ?(?) and ?(?) to yield a desired closed-loop characteristic polynomial ??(?) can be determined if the dose-response model parameters ?? , ? = 1,2,3,4 are known. In APPC, these parameters are provided by the adaptive law (76)-(77). In sum, APPC is realized by the pole placement control law (78) with ?(?) and ?(?) determined by (83) with the aid of the adaptive law (76)-(77). The stability of the resulting APPC can be established using available procedures [58]. One last consideration is concerned with the stability of control law. Specifically, ?(?) = ?(?)?(?) ? and ??(?)?(?) is not guaranteed to be Hurwitz. Noting that ??(?)?(?) ??(?)?(?)?(?) + ?(?)?(?) = 0 from (78), an alternative realization of ?(?) can be obtained as follows: 1 ?(?) = ?(?) ? [? (?)?(?)?(?) + ?(?)?(?)] ?(?) ? (84) 1 ?(?) = [?(?) ? ??(?)?(?)]?(?) ? ?(?) ?(?) ?(?) where ?(?) is a monic Hurwitz polynomial of degree 3 for proper filtering of ?(?) and ?(?). 98 7.7. In-Silico Implementation and Simulation The controller was in-silico evaluated using a case study of regulating a cardiovascular variable (CO) with a sedative (propofol), under a wide range of transport delay and pharmacological profiles as well as desired responses. For the sake of in-silico evaluation, we used the 30 randomly created in-silico subjects in Chapter 7.3. In all in-silico subjects, the baseline CO (y0) was set at 3.0 lpm, while a multitude of desired responses, in terms of both magnitude and rate, was employed: 2.1, 2.4, 2.7 lpm with the time constants of 5, 2.5, and 1.7 min. The parameters associated with the indirect semi-APPC were chosen empirically as follows: ? = ? ?4?4 , ? (?) = (? + ?)5 with ? > 0, ?(?) = (? + ?)3 with ? = ?0, and ?0 = 0.02. Since the goal of this case study was to examine the efficacy of APPC for medication infusion problems with transport delay, these parameters were not rigorously optimized. The values of ? (and thus ?0 as well) and ? were likewise tuned for each desired response. The semi-APPC was then implemented in the discrete-time domain using the zero-order hold method, in which a sampling interval of 5 sec was used for control computation. To ensure subject safety against over-dosing, a target-dependent (i.e., ?-dependent) upper bound of infusion rate was augmented to the control law (0.8 mg/kg/min for 2.7 lpm target, 1.2 mg/kg/min for 2.4 lpm target, and 1.6 mg/kg/min for 2.1 lpm target). In this way, a total of 270 in- silico simulations was performed and analyzed for semi-APPC. To investigate the significance of on-line transport delay adaptation, the semi-APPC was compared with controllers without explicit account for transport delay. First, semi-APPC 99 with transport delay L fixed at nominal (i.e., 75 sec) and worst-case (i.e., 100 sec) values were simulated in the same 30 in-silico subjects. Second, the population-based PID controller used in Chapter 7.3 was also simulated in the same 30 in-silico subjects. Then, the performance of these controllers was quantitatively compared for transient and steady-state behaviors. 7.8. Results and Discussion The objective of this chapter was to investigate adaptive control design for closed-loop medication infusion problems with on-line adaptation of transport delay originating from PD and subject monitors. We used the Pad? approximation and indirect APPC techniques to streamline control design. The results provided below suggest that adaptation of transport delay may benefit in minimizing the variability in closed-loop response against a wide range of dose-response variability. Details follow. Comparing semi-APPC, semi-APPC with nominal ?, semi-APPC with worst-case ?, and population-based PID control, semi-APPC outperformed all the other controllers. Figure 25 presents RMSE, median absolute percentage error (MDAPE) [59], and wobble [59] metrics computed for all four controllers associated with 5 min time constant (i.e., am = 0.2 s-1; results for other time constants (a = 0.4 s-1 and a = 0.6 s-1m m ) exhibited similar trends). Figure 26 shows representative examples of the responses of 30 in-silico subjects associated with all four controllers under target CO of 2.7 lpm with 5 min time constant. Figure 27 shows the variability of steady-state set point tracking errors (as measured by the standard deviation of the set point tracking errors between settling time and settling 100 time + 10 min of the desired response) associated with all four controllers. The semi- APPC was in particular superior to all the other controllers in terms of RMSE (by 16%, 12%, and 35% against semi-APPC with nominal L, semi-APPC with worst-case L, and population-based PID control), MDAPE (by 17%, 18%, and 42% against semi-APPC with nominal L, semi-APPC with worst-case L, and population-based PID control), and wobble (by 16%, 16%, and 30% against semi-APPC with nominal L, semi-APPC with worst-case L, and population-based PID control). The semi-APPC, semi-APPC with nominal L, and semi-APPC with worst-case L exhibited comparable median percentage error (MDPE) and divergence performance, which was superior to population-based PID control (not shown). The semi-APPC also outperformed all the other controllers in terms of the consistency in steady-state response: semi-APPC consistently exhibited smaller steady-state set point tracking error variability than all the other controllers (by 28%, 30%, and 48% against semi-APPC with nominal L, semi-APPC with worst-case L, and population-based PID control) under all desired clinical effect response. The population- based PID control suffered from the largest average variability and variability thereof, indicating that neglecting transport delay in closed-loop medication infusion may yield drastic degradation in performance in case dose-response delay is not negligibly small. Although semi-APPC with nominal and worst-case L were superior to population-based PID control, they fell short of semi-APPC in both transient and steady-state response characteristics, implying that the use of a population-based transport delay is not ideal in closed-loop medication infusion control. 101 Fig. 25: Root-mean-squared error (RMSE), median absolute percentage error (MDAPE), and wobble metrics computed for semi-APPC, semi-APPC with nominal transport delay, semi-APPC with worst-case transport delay, and population-based PID control associated with 5 min time constant. For each target cardiac output (CO), each bar denotes from the left semi-APPC, semi-APPC with nominal transport delay, semi-APPC with worst-case transport delay, and population-based PID control. 102 Fig. 26: Representative examples of the responses of 30 in-silico subjects to target cardiac output (CO) of 2.7 lpm with 5 min time constant associated with semi-APPC, semi-APPC with nominal transport delay, semi-APPC with worst-case transport delay, and population-based PID control. All in all, the results suggest that (1) transport delay may need to be accounted for in medication infusion control problems, and that (2) transport delay may even need to be adapted to secure robustness in control performance despite uncertainty and variability in dose-response relationship. 103 Fig. 27: Steady-state error variability of 30 in-silico subjects associated with semi-APPC, semi-APPC with nominal transport delay, semi-APPC with worst-case transport delay, and population-based PID control. 104 Chapter 8: Conclusions and Future Work In the dissertation, model and control of medication infusion are investigated. The unique contributions of our work are given as follows: (1) The classical medication infusion model is very complex which results in the difficulty of controller design. We proposed a control-oriented model for multiple medication infusion. After comparison of classical model and our proposed model, our model showed an improved fitting error and AIC which indicates it has better predictive capability and model parsimony. (2) Due to complexity of traditional medication infusion model, most of the existing work of MIMO medication infusion controller are not model-based. Even if for model-based controllers, some of them made some unrealistic assumptions such as all parameters of model are known. Based on our proposed model, we designed a semi-adaptive controller for multiple infusion of interacting medications. To the best of our knowledge, our proposed controller is the first work based on rigorous Lyapunov stability proof. The validity and performance of the approach was also examined in simulations. (3) Based on our proposed controller, we further improved its performance by combining switching controller. It was shown that the switching control could overall provide added benefit to improve command tracking performance via more accurate estimation of subject?s dose-response relationship. 105 (4) Furthermore, targets of multiple medication infusion should be coordinated. However, it was not considered in other controller design. Our work is first one which designed an coordination algorithm to optimize the targets. We also demonstrated that the proposed approach could achieve consistent reference target tracking in the presence of large inter- individual variability in dose-response relationships, and adjust unachievable reference targets to achievable ones while minimizing the medication use by exploiting the inter- medication synergy. (5) In the end, to verify our proposed controller, pilot pig experiments were performed and the results were examined. Because transport delay were observed in the model, so we investigated a semi-adaptive indirect pole placement control for infusion of medications with large transport delay in dose-response dynamics. The benefits of it were shown from comparison with PID controller. The following work need to be conducted in the future: (1) Although we have designed the switching controller, coordination algorithm and controller which incorporates transport delay, we need to combine them together. The combined controller can deal with the existing chanllenges we found at the same time. (2) Our proposed controller was based on a two-input two-output model. This is not always the case in real clinical scenarios. Therefore, we need expand our proposed controller to enable infusion of more than two medications and control of more than two endpoints. 106 (3) Because most of my work was focused on controller design and simulations, we only implemented our SISO controller in some preliminary experiments. In the future, we need to experimentally evaluate our MIMO controller more rigorously. 107 Bibliography [1] D. SILVA, ?System identification and control for general anesthesia based on parsimonious Wiener models,? 2012. [2] N. Liu et al., ?Closed-loop coadministration of propofol and remifentanil guided by bispectral index: A randomized multicenter study,? Anesth. Analg., vol. 112, no. 3, pp. 546?557, 2011. [3] T. M. Hemmerling, S. Charabati, C. Zaouter, C. Minardi, and P. a. Mathieu, ?A randomized controlled trial demonstrates that a novel closed-loop propofol system performs better hypnosis control than manual administration,? Can. J. Anesth., vol. 57, no. 8, pp. 725?735, 2010. [4] T. J. Sieber, C. W. Frei, M. Derighetti, P. Feigenwinter, D. Leibundgut, and A. M. Zbinden, ?Model-based automatic feedback control versus human control of end- tidal isoflurane concentration using low-flow anaesthesia.,? Br. J. Anaesth., vol. 85, no. 6, pp. 818?825, 2000. [5] B. L. Sng, H. S. Tan, and a. T. H. Sia, ?Closed-loop double-vasopressor automated system vs manual bolus vasopressor to treat hypotension during spinal anaesthesia for caesarean section: A randomised controlled trial,? Anaesthesia, vol. 69, no. 1, pp. 37?45, 2014. [6] J. Salinas et al., ?Closed loop and decision assist resuscitation of burn patients,? J. Trauma Acute Care Surg., vol. 64, no. 4, pp. S321?S332, 2008. [7] J. Rinehart, C. Lee, C. Canales, A. Kong, Z. Kain, and M. Cannesson, ?Closed- loop fluid administration compared to anesthesiologist management for hemodynamic optimization and resuscitation during surgery: An in vivo study,? 108 Anesth. Analg., vol. 117, no. 5, pp. 1119?1129, 2013. [8] J. Rinehart, E. Chung, C. Canales, and M. Cannesson, ?Intraoperative stroke volume optimization using stroke volume, arterial pressure, and heart rate: Closed- loop (learning intravenous resuscitator) versus anesthesiologists,? J. Cardiothorac. Vasc. Anesth., vol. 26, no. 5, pp. 933?939, 2012. [9] G. C. Kramer et al., ?Closed-loop control of fluid therapy for treatment of hypovolemia.,? J. Trauma, vol. 64, no. 4, pp. S333?S341, 2008. [10] N. F. Chaisson, R. a Kirschner, D. J. Deyo, J. A. Lopez, D. S. Prough, and G. C. Kramer, ?Near-infrared spectroscopy-guided closed-loop resuscitation of hemorrhage.,? J. Trauma, vol. 54, no. 5, pp. S183?S192, 2003. [11] M. Merouani et al., ?Norepinephrine weaning in septic shock patients by closed loop control based on fuzzy logic.,? Crit. Care, vol. 12, no. 6, p. R155, 2008. [12] M. Verdon, P. Merlani, T. Perneger, and B. Ricou, ?Burnout in a surgical ICU team,? Intensive Care Med., vol. 34, no. 1, pp. 152?156, 2008. [13] X.-C. Zhang et al., ?Job burnout among critical care nurses from 14 adult intensive care units in northeastern China: a cross-sectional survey,? BMJ Open, vol. 4, pp. e004813?e004813, 2014. [14] N. Embriaco, L. Papazian, N. Kentish-Barnes, F. Pochard, and E. Azoulay, ?Burnout syndrome among critical care healthcare workers,? Curr. Opin. Crit. Care, vol. 13, no. 5, pp. 482?488, 2007. [15] A. Gentilini et al., ?Multitasked closed-loop control in anesthesia,? IEEE Eng. Med. Biol. Mag., vol. 20, no. 1, pp. 39?53, 2001. [16] G. A. Dumont, A. Martinez, and J. M. Ansermino, ?Robust control of depth of 109 anesthesia,? Int. J. Adapt. Control Signal Process., vol. 23, no. 5, pp. 435?454, 2009. [17] S. Yelneedi, L. Samavedham, and G. P. Rangaiah, ?Advanced control strategies for the regulation of hypnosis with propofol,? Ind. Eng. Chem. Res., vol. 48, pp. 3880?3897, 2009. [18] J. O. Hahn, G. a. Dumont, and J. M. Ansermino, ?Robust closed-loop control of hypnosis with propofol using WAV CNS index as the controlled variable,? Biomed. Signal Process. Control, vol. 7, no. 5, pp. 517?524, 2012. [19] A. L. G. Caruso, T. W. Bouillon, P. M. Schumacher, E. Zanderigo, and M. Morari, ?Control of drug administration during monitored anesthesia care,? IEEE Trans. Autom. Sci. Eng., vol. 6, no. 2, pp. 256?264, 2009. [20] F. N. Nogueira, T. Mendon?a, and P. Rocha, ?Controlling the depth of anesthesia by a novel positive control strategy,? Comput. Methods Programs Biomed., vol. 114, no. 3, pp. e87?e97, 2014. [21] A. Gentilini, C. Schaniel, M. Morari, C. Bieniok, R. Wymann, and T. Schnider, ?A new paradigm for the closed-loop intraoperative administration of analgesics in humans,? IEEE Trans. Biomed. Eng., vol. 49, no. 4, pp. 289?299, 2002. [22] C. M. Ionescu, R. D. Keyser, B. C. Torrico, T. D. Smet, M. M. Struys, and J. E. Normey-Rico, ?Robust predictive control strategy applied for propofol dosing using BIS as a controlled variable during anesthesia,? IEEE Trans. Biomed. Eng., vol. 55, no. 9, pp. 2161?2170, 2008. [23] K. van Heusden et al., ?Design and Clinical Evaluation of Robust PID Control of Propofol Anesthesia in Children,? IEEE Trans. Control Syst. Technol., vol. 22, no. 110 2, pp. 491?501, 2014. [24] W. M. Haddad, J. M. Bailey, T. Hayakawa, and N. Hovakimyan, ?Neural Network Adaptive Output Feedback Control for Intensive Care Unit Sedation and Intraoperative Anesthesia,? IEEE Trans. Neural Netw., vol. 18, no. 4, pp. 1049? 1066, 2007. [25] J. A. M?ndez, S. Torres, J. A. Reboso, and H. Reboso, ?Adaptive computer control of anesthesia in humans.,? Comput. Methods Biomech. Biomed. Engin., vol. 12, no. 6, pp. 727?734, 2009. [26] K. Volyanskyy, W. M. Haddad, and J. M. Bailey, ?Adaptive disturbance rejection control for compartmental systems with application to intraoperative anesthesia influenced by hemorrhage and hemodilution effects,? Int. J. Adapt. Control Signal Process., vol. 23, pp. 1?29, 2009. [27] Y. Sawaguchi, E. Furutani, G. Shirakami, M. Araki, and K. Fukuda, ?A model- predictive hypnosis control system under total intravenous anesthesia,? IEEE Trans. Biomed. Eng., vol. 55, no. 3, pp. 874?887, 2008. [28] K. Soltesz, J. O. Hahn, T. H?gglund, G. a. Dumont, and J. M. Ansermino, ?Individualized closed-loop control of propofol anesthesia: A preliminary study,? Biomed. Signal Process. Control, vol. 8, no. 6, pp. 500?508, 2013. [29] M. M. Da Silva, T. Wigren, and T. Mendonca, ?Nonlinear identification of a minimal neuromuscular blockade model in anesthesia,? IEEE Trans. Control Syst. Technol., vol. 20, no. 1, pp. 181?188, 2012. [30] D. J. F. Nieuwenhuijs et al., ?Response surface modeling of remifentanil-propofol interaction on cardiorespiratory control and bispectral index.,? Anesthesiology, vol. 111 98, no. 2, pp. 312?322, 2003. [31] M. J. Mertens, E. Olofsen, F. H. Engbers, A. G. Burm, J. G. Bovil, and Lj. Vuyk, ?Propofol Reduces Perioperative Remifentanil Requirements in a Synergistic Manner,? J. Am. Soc. Anesthesiol., vol. 99, no. 2, pp. 347?359, 2003. [32] S. E. Kern, G. Xie, J. L. White, and T. D. Egan, ?A Response Surface Analysis of Propofol?Remifentanil Pharmacodynamic Interaction in Volunteers,? J. Am. Soc. Anesthesiol., vol. 100, no. 6, pp. 1373?1381, 2004. [33] C. F. Minto, T. W. Schnider, T. G. Short, K. M. Gregg, A. Gentilini, and S. L. Shafer, ?Response Surface Model for Anesthetic Drug Interactions,? Anesthesiology, vol. 92, pp. 1603?1616, 2000. [34] T. W. Bouillon et al., ?Pharmacodynamic interaction between propofol and remifentanil regarding hypnosis, tolerance of laryngoscopy, bispectral index, and electroencephalographic approximate entropy.,? Anesthesiology, vol. 100, no. 6, pp. 1353?1372, 2004. [35] Y. C. Chou, M. F. Abbod, J. S. Shieh, and C. Y. Hsu, ?Multivariable fuzzy logic/self-organizing for anesthesia control,? J. Med. Biol. Eng., vol. 30, no. 5, pp. 297?306, 2010. [36] Q. Lu and M. Mahfouf, ?Multivariable self-organizing fuzzy logic control using dynamic performance index and linguistic compensators,? Eng. Appl. Artif. Intell., vol. 25, no. 8, pp. 1537?1547, 2012. [37] H. H. Lin, C. L. Beck, and M. J. Bloom, ?On the use of multivariable piecewise- linear models for predicting human response to anesthesia,? IEEE Trans. Biomed. Eng., vol. 51, no. 11, pp. 1876?1887, 2004. 112 [38] H. H. Lin, C. Beck, and M. Bloom, ?Multivariable LPV control of anesthesia delivery during surgery,? Proc. Am. Control Conf., no. 1, pp. 825?831, 2008. [39] C. M. Ionescu, I. Nascu, and R. De Keyser, ?Lessons learned from closed loops in engineering: towards a multivariable approach regulating depth of anaesthesia,? J. Clin. Monit. Comput., vol. 28, no. 6, pp. 537?546, 2014. [40] I. Kolmanovsky, E. Garone, and S. Di Cairano, ?Reference and command governors: A tutorial on their theory and automotive applications,? Proc. Am. Control Conf., pp. 226?241, 2014. [41] E. Gilbert, I. Kolmanovsky, and K. T. Tan, ?Discrete?time reference governors and the nonlinear control of systems with state and control constraints,? Int. J. Robust Nonlinear Control, vol. 5, no. October 1994, pp. 487?504, 1995. [42] E. G. Gilbert and C. J. Ong, ?Constrained linear systems with hard constraints and disturbances: An extended command governor with large domain of attraction,? Automatica, vol. 47, no. 2, pp. 334?340, 2011. [43] J. O. Hahn, G. A. Dumont, and J. M. Ansermino, ?A direct dynamic dose-response model of propofol for individualized anesthesia care,? IEEE Trans. Biomed. Eng., vol. 59, no. 2, pp. 571?578, 2012. [44] J. M. Ansermino, P. Brooks, D. Rosen, C. a. Vandebeek, and C. Reichert, ?Spontaneous ventilation with remifentanil in children,? Paediatr. Anaesth., vol. 15, no. 2, pp. 115?121, 2005. [45] MATLAB Statistics Toolbox User?s Guide. Natick, MA: MathWorks: Natick, MA: MathWorks, 2010. [46] a. E. Rigby-Jones et al., ?Remifentanil-midazolam sedation for paediatric patients 113 receiving mechanical ventilation after cardiac surgery,? Br. J. Anaesth., vol. 99, no. 2, pp. 252?261, 2007. [47] S. Khosravi, J. O. Hahn, G. a. Dumont, and J. M. Ansermino, ?A monitor- decoupled pharmacodynamic model of propofol in children using state entropy as clinical endpoint,? IEEE Trans. Biomed. Eng., vol. 59, no. 3, pp. 736?743, 2012. [48] J. J. Slotine and W. Li, Applied Nonlinear Control. Englewood Cliffs: Prentice- Hall, 1991. [49] P. a. Ioannou, J. Sun, and P. List, Robust Adaptive Control. 1996. [50] X. Jin, C. S. Kim, G. A. Dumont, J. M. Ansermino, and J. O. Hahn, ?A semi- adaptive control approach to closed-loop medication infusion,? Int. J. Adapt. Control Signal Process., 2016. [51] A. Koroglu, H. Teksan, O. Sagir, A. Yucel, H. I. Toprak, and O. M. Ersoy, ?A comparison of the sedative, hemodynamic, and respiratory effects of dexmedetomidine and propofol in children undergoing magnetic resonance imaging,? Anesth. Analg., vol. 103, no. 1, pp. 63?67, 2006. [52] D. Liberzon, Switching in systems and control. New York: Springer Science & Business Media, 2003. [53] D. Liberzon and A. S. Morse, ?Basic problems in stability and design of switched systems,? IEEE Control Syst., vol. 19, no. 5, pp. 59?70, 1999. [54] X. Jin, C. Kim, G. A. Dumont, J. M. Ansermino, and J. Hahn, ?A semi-adaptive control approach to closed-loop medication infusion,? Int. J. Adapt. Control Signal Process., vol. 31, no. 2, pp. 240?254, 2017. [55] X. Jin, C. S. Kim, G. A. Dumont, J. M. Ansermino, and J. O. Hahn, ?A semi- 114 adaptive control approach to closed-loop medication infusion,? Int. J. Adapt. Control Signal Process., vol. 31, pp. 240?254, 2017. [56] X. Jin, C.-S. Kim, S. T. Shipley, G. A. Dumont, and J.-O. Hahn, ?Coordinated semi-adaptive closed-loop control for infusion of two interacting medications,? Int. J. Adapt. Control Signal Process., vol. 32, no. 1, pp. 134?146, 2018. [57] M. Krstic and A. Banaszuk, ?Multivariable adaptive control of instabilities arising in jet engines,? Control Eng. Pract., vol. 14, no. 7, pp. 833?842, 2006. [58] P. Ioannou and J. Sun, Robust Adaptive Control. Mineola, NY, USA: Dover Publications, 2012. [59] J. R. Varvel, D. L. Donoho, and S. L. Shafer, ?Measuring the predictive performance of computer-controlled infusion pumps,? J. Pharmacokinet. Biopharm., vol. 20, no. 1, pp. 63?94, 1992. 115