UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Automation in clinical anesthesia Bibian, Stephane 2006

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2006-199398.pdf [ 31.31MB ]
Metadata
JSON: 831-1.0065536.json
JSON-LD: 831-1.0065536-ld.json
RDF/XML (Pretty): 831-1.0065536-rdf.xml
RDF/JSON: 831-1.0065536-rdf.json
Turtle: 831-1.0065536-turtle.txt
N-Triples: 831-1.0065536-rdf-ntriples.txt
Original Record: 831-1.0065536-source.json
Full Text
831-1.0065536-fulltext.txt
Citation
831-1.0065536.ris

Full Text

AUTOMATION, IN CLINICAL ANESTHESIA by Stephane Bibian Master of Applied Science, University of British Columbia, Vancouver, Canada, 1999 Bachelor of Science, E.S.I.E.E., Amiens, France, 1997 A thesis submitted in partial fulfillment of the requirements for the degree of DOCTORATE OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Electrical and Computer Engineering) T H E U N I V E R S I T Y OF B R I T I S H C O L U M B I A July 2006 ©Stephane Bibian, 2006 Abstract This thesis investigates the design and performance of a controller for the maintenance of anesthesia during surgery. The controller is designed to be robustly stable for a large population of patients. Even though anesthetic drugs are amongst the most dangerous drugs used in today's clinical setting, anesthesia procedures are known to be very safe. Hence, the impact of automation in anesthesia in terms of patients' safety cannot be clearly established. However, there are a number of significant clinical advantages to be gained by closing the loop: 1. Recent evidences suggest that most patients undergoing anesthesia procedures are overdosed. This is one of the main reasons for patients' discomfort and slow recovery. Literature suggests that closed-loop systems can significantly reduce drug consumption and lessen recovery times, thus improving the patient outcome while reducing drug-associated costs and bed occupancy. 2. Anesthesiologists have access to intravenous agents with fast onset of action and fast metabolism. Using a closed-loop controller would allow for an infusion-type titration that provides smoother tran-sitions, thus avoiding the respiratory and hemodynamic depression observed in a bolus-based manual regimen. 3. Closed-loop controllers are also particularly well-suited for solving complex optimization problems. The profound synergy that exists between intravenous anesthetics and opioids could then be fully exploited. This could be a significant factor contributing to a reduction in drug usage and the improvement of patients' comfort. This project is particularly challenging. In particular: 1. There is no accepted measure of depth of anesthesia. Hence, it is necessary to work at the conceptual and sensor levels in order to define adequate feedback measures. 2. Drug effect modeling suffers from many shortcomings. In particular, published studies are often not in agreement regarding model parameters. 3. Uncertainty of dose/response models is daunting. Measuring this uncertainty is necessary in order to ensure stability of the control design. ii ABSTRACT AND THESIS ORGANIZATION i i i While the anesthesia closed-loop concept has already been investigated in the past, no breakthrough has yet been achieved. We feel it is necessary to investigate the anesthesia system from a control engineering perspective. This thesis is divided into two distinct parts. Part A contains the first 4 chapters and presents a thorough introduction to clinical anesthesia. The main concepts, terminology and issues are covered, including anesthesia monitors and basic pharmacology principles. A review of the prior closed-loop control attempts is presented in Chapter 4. Part B contains the chapters 5 to 8. In these chapters, we investigate a new sensor technology to quantify both cortical and autonomic activity. This technology is used to derive drug effect models, from which uncertainty bounds are derived. Based on this uncertainty analysis, we derive robustly stable controllers achieving clinically adequate performances. Finally, we invite the readers to refer to Chapter 9 for a complete synopsis and summary of this thesis. Table of Contents Abstract i i Table of Contents iv List of Tables ix List of Figures x i List of Acronyms xx Acknowledgements xxii Dedication xxii i P A R T A : The Anesthesia Framework 1 1 Clinical Anesthesia: Terminology, Concepts, and Issues 2 1.1 Clinical Anesthesia: A Key Specialty in Critical Care 2 1.1.1 The Early Days 3 1.1.2 Risks and Outcome in Anesthesia 3 1.2 Modern Concepts 6 1.2.1 The Role of the Anesthesiologist 6 1.2.2 Functional Components of Clinical Anesthesia 7 1.2.3 Modern Balanced Anesthesia 8 1.3 A n Extensive Pharmacopoeia 8 1.3.1 Anesthetics (i.e. hypnotics) 8 1.3.2 Opioids • 10 1.3.3 Neuromuscular Blocking Agents (NMBs) 10 1.3.4 Drugs: Action, Effect and Interaction 11 1.4 Conduct of Anesthesia 12 1.5 Drug Delivery and Monitoring in Anesthesia 14 1.5.1 Drug Delivery Devices 14 1.5.2 Anesthesia Monitors 15 1.6 Summary 16 iv ABSTRACT AND THESIS ORGANIZATION v 2 Quantifying Depth of Anesthesia: a Review 17 2.1 The Use of Surrogate Measures 18 2.2 Quantifying the Depth of Hypnosis: a Review 19 2.2.1 Review of Tools and Techniques 19 2.2.2 Clinical Relevance, Interest and Potentials of the BIS Monitor 26 2.3 Quantifying the Depth of Analgesia 29 2.3.1 Using Physiological Measures 29 2.3.2 Using Heart Rate Variability 30 2.4 Summary 31 3 Modeling Anesthetic Drugs: the Traditional Approach 32 3.1 Pharmacokinetics 32 3.1.1 Principles and Concepts 33 3.1.2 Modeling 34 3.1.3 Literature Review 38 3.1.4 Concluding Remarks 40 3.2 Pharmacodynamics 42 3.2.1 The Dose/Response Relationship: a Steady-State Model 42 3.2.2 The keo Parameter: a Measure of the Effect Dynamics 45 3.2.3 Literature Review 46 3.2.4 Concluding Remarks 48 3.3 Summary 51 4 Closing the Loop in Anesthesia: a Review 53 4.1 Review 53 4.2 Summary 54 PART B: Sensing, Modeling, and Control 58 5 Quantifying Cortical and Autonomic Activity Using Wavelets 59 5.1 The Discrete Wavelet Transform 61 5.1.1 Standard Wavelet Dyadic Decomposition 61 5.1.2 Change in Phase and Wavelets 62 5.2 Estimating the Anesthetic Drug Effect: the W A V C N S 62 5.2.1 Concepts and Derivation 62 5.2.2 Practical Issues and Implementation 73 5.2.3 Clinical Results 76 ABSTRACT AND THESIS ORGANIZATION vi 5.3 Measuring the Autonomic Nervous System (A-NS) State: the W A V A N S Index 82 5.3.1 The Use of Heart Rate Variability 82 5.3.2 The W A V A N S 84 5.3.3 Implementation Issues 85 5.3.4 Case Reports 86 5.4 Dynamic Behavior 91 5.5 Summary 94 6 A System Oriented Approach to Pharmacodynamic modeling 95 6.1 A New System Oriented Approach to P D Modeling 95 6.1.1 Sheiner's Approach 96 6.1.2 Proposed Methodology : 101 6.2 Application to Propofol and W A V C N S Index 103 6.2.1 Clinical Data 103 6.2.2 P D Identification - Traditional Approach 104 6.2.3 P D Identification - New Approach 105 6.2.4 Adequacy of the Identification Data 107 6.2.5 Identification Results 107 6.2.6 Model Validation 109 6.2.7 Nominal Propofol P D Models 112 6.3 Discussion 114 6.3.1 Clinical Adequacy of the Identification Data 114 6.3.2 System Oriented vs. Traditional Approach 116 6.3.3 Dose vs. Response Characteristics 120 6.3.4 Frequency Response Characteristic 121 6.4 Summary 123 7 Managing PKPD Uncertainty 126 7.1 Quantifying Uncertainty 126 7.1.1 The Relative Uncertainty Framework 128 7.1.2 Selection of a Near Optimal Nominal Model 128 7.2 Application to Propofol P K P D Models 134 7.2.1 Origins of P K P D Uncertainty 134 7.2.2 Propofol P K P D Uncertainty Results 138 7.3 Reducing Uncertainty 141 7.3.1 Passive Methods 141 7.3.2 Total and Partial Self Tuning , 144 7.4 Summary 149 ABSTRACT AND THESIS ORGANIZATION vii 8 SISO Control Design 155 8.1 Patient Simulator 156 8.1.1 Simulator Structure 156 8.1.2 Noise Modeling 158 8.1.3 Patient Simulator Output 161 8.1.4 Tests Setup and Performance Requirements 161 8.2 P ID Loop Shaping Design 164 8.2.1 First Design 164 8.2.2 Second Design 166 8.2.3 Third Design 166 8.3 Hoo Control Design 168 8.3.1 Hoo Design Principle 168 8.3.2 Implementation 170 8.3.3 Application to Anesthesia Control 170 8.4 Reducing the Uncertainty 172 8.4.1 Accounting for age 174 8.4.2 Identifying the P K time delay 174 8.5 The Clinicians' Point of View 176 8.5.1 Increased Controller Gain 178 8.5.2 Use of a Smith Predictor Structure 179 8.5.3 Faster Sensing Dynamics 182 8.5.4 Improved Performance: Important Trade-offs 183 8.6 Summary 184 9 Conclusion, Contributions and Recommendations 187 9.1 Significance 187 9.2 Synopsis and Contributions 188 9.2.1 Synopsis 189 9.2.2 Contributions and Implications 190 9.3 Future Work 192 Bibliography 197 A B S T R A C T A N D THESIS ORGANIZATION vi i i A Pharmacopoeia 209 A . l Vapour Anesthetics 209 A.2 Intravenous Anesthetics 209 A . 3 Opioids 210 B Propofol and Remifentanil P K P D 213 B . l Propofol 213 B . l . l Pharmacokinetics 213 B.1.2 Pharmacodynamics 216 B.2 Remifentanil 219 B.2.1 Pharmacokinetics 219 B.2.2 Pharmacodynamics 220 C Propofol P K P D Nominal and Uncertainty Models 222 D Clinical Studies 229 D . l Arthroscopy Observational Study 229 D . l . l Protocol 229 D . l . 2 Results 230 D. l .3 Comments and Unexpected Results 230 D.2 Laryngeal Mask Airway Study 230 D.2.1 Protocol 232 D.2.2 L M A Results 235 D.3 Electro-Convulsive Shock Therapy Study 236 E Survey 242 List of Tables 1.1 Contemporary anesthetic mortality rates (adapted from Brown [5]) 5 3.1 Propofol P K parameter sets from [101] where BW stands for body weight, ven= 1 for venous sampling, ven= 0 for arterial sampling, bol= 1 for bolus administration, and bol= 0 for infusion administration. 39 3.2 Parameter estimates from the N O N M E M analysis ([101]) 39 3.3 Pharmacokinetic parameters of remifentanil as a function of the age and lean body mass (from [103]). 40 4.1 Prior Art: a Literature Survey (Part 1) 55 4.2 Prior Art: a Literature Survey (Part If) 56 6.1 PD models obtained from the proposed approach with E m a 3 : = 0 and Eo=100 108 6.2 Population-normed PD models with age as a linear covariate 114 6.3 Goodness of fit of patient-specific models obtained from the traditional PD approach. 118 6.4 Goodness of fit of patient-specific models obtained from the system-oriented PD approach 118 6.5 Goodness of fit of population-normed models obtained from the traditional PD approach 119 6.6 Goodness of fit of population-normed models obtained from the system-oriented PD approach. . . . 119 7.1 Propofol P K P D Inter-patient Variability 135 7.2 E C T Intra- and Inter-patient Variability 137 8.1 Control Performance Requirements 163 8.2 PID Controller Characteristics 165 8.3 Control Performance - PID Loop Shaping Design 167 8.4 Hoo Controller Characteristics 171 8.5 Control Performance - Hoc Control Design 172 8.6 HQO Controller Characteristics 174 8.7 Control Performance - HQO Control Design based on Age 175 8.8 Control Performance Requirements - Clinicians' Point of View 177 8.9 HQO Controller Characteristics with Increased Gain 179 8.10 Control Performance - Improved HQO Control Designs 180 ix LIST OF TABLES x B . l Propofol P K parameter sets (hybrid form) 214 B.2 Propofol effect site half-life and pharmacodynamic parametric constants 217 B.3 Remifentanil PK parameter sets (hybrid form) 219 B . 4 Remifentanil effect site half-life and pharmacodynamic parametric constants 221 C. l P K P D Propofol nominal models - part I 223 C.2 P K P D Propofol nominal models - part II 224 C.3 P K P D Propofol nominal models - part III 225 C.4 P K P D Propofol uncertainty weights - part I 226 C.5 P K P D Propofol uncertainty weights - part II 227 C. 6 P K P D Propofol uncertainty weights - part III 228 D. l Patients' data and titration. Doses are expressed in \ig (sufentanil, fentanyl, remifentanil) or mg (propofol) 231 D.2 L M A Study population (patients #001 to #049) 233 D.3 L M A Study population (patients #050 to #076) 234 D.4 L M A study population - demographics and dosage summary 235 D.5 E C T Study Population 238 D.6 Thiopental pharmacokinetic models 240 D.7 Thiopental identified pharmacodynamic models 241 List of Figures 2.1 Therapeutic window of the E E G for measuring the opioid effect (from [31]) 18 2.2 Simultaneous records from the left and right motor regions, illustrating the changes in frequency during cyclopropoane anesthesia and recovery, (a) awake (b) breathing cyclopropane for 2 min. (c) 1 min later (d) immediately on substitution of room air for cyclopropane (e) 3 seconds later (f) 2 min. later (g) 5 min. later (h) complete recovery (from [35]) 20 2.3 Awake and anesthetized E E G and their power spectrum 21 2.4 Changes in latencies in a signal can be tracked by the bicoherence index (from [52]) 23 2.5 Bispectral Index Scale and its meaning (from [59]) 24 2.6 Midlatency Auditory Evoked Potentials obtained from a responsive (awake) and non-responsive (asleep) dog subject to tail clamping. The waveforms have been obtained after averaging of 1000 samples (adapted from [64]) 25 2.7 Time of arousal following discontinuation of propofol and remifentanil infusion (from [70]) 27 3.1 Simulation of the time course of the percentage of a thiopental bolus in the central blood pool, viscera, lean tissue and fat as a function of time, assuming no elimination clearance (adapted from [97]). . . 33 3.2 Time course of plasma concentration with two distinct phases: distribution and elimination (adapted from [98]) 35 3.3 2-compartment pharmacokinetic model. A third compartment is often added for improved accuracy (adapted from [99]) 36 3.4 Bolus vs. infusion pharmacokinetic of propofol (using [101]). (a) Impulse response (b) Frequency response 41 3.5 Dose-response relationships differ according to the drug potency, efficacy, slope and subjects' variability. 43 3.6 Drug dose-response curves can be obtained by observing whether a patient responds to a particular stimulus at a given concentration (from [113]) 44 3.7 Different dose/response characteristics (a) On/Off effect (b) The observed endpoint cannot charac-terize the effect of small doses (c) Ideal characteristic 44 3.8 Dose/response model for Alfentanil and using the probability of maximal E E G response (from [31]). 45 3.9 Effect and plasma concentration of fentanyl and alfentanil following rapid intravenous administration (adapted from Stanski [89]) 46 xi LIST OF FIGURES xii 3.10 Individual (dashed lines) and mean (solid lines) Hill non-linear characteristics of published pharma-codynamic models, (a) Propofol study from Kuizenga et al. [118] (b) Remifentanil study from Egan et al. [102] 47 3.11 Drug interaction between propofol and fentanyl (adapted from [120]) 48 3.12 Simulation of propofol and remifentanil pharmacodynamic interactions (adapted from [23]) 49 3.13 Measured and predicted time courses of the BIS following a series of induction sequence administered to a volunteer (from [118]). A pharmacodynamic model obtained based on the first repeat was used to predict the BIS time course. Note the good agreement between measurement and prediction for the first repeat. Results worsen considerably for the two other repeats 50 3.14 Development of acute tolerance to analgesia during short-term constant-rate infusion of remifentanil (adapted from [121]) 51 5.1 Dyadic frequency tiling with L = 3 61 5.2 Changes in latencies in a signal can also be tracked by the wavelet coefficients (from [151]). (a) Time series signal. Note how the signal pattern changes significantly when changing the phase of one of the frequency component, (b) 'Periodogram' of the wavelet coefficients (note how the change in phase is clearly noticeable in the wavelet domain), (c) Classical Fourier periodogram. (d) Time profile of the phase change 63 5.3 The function h quantifies the state of a system by attributing a unique value to each of the operating modes of the system. If the system evolves in a monotonous fashion from state a to state b, it is required that h is also monotonous 64 5.4 The characteristic of the function h depends on the selection of the feature function / . (a) Both hi and hi are valid functions, whereas h$ and h± cannot be used to represent adequately intermediate states, (b) The linearity of a given function h can be assessed by introducing new data sets, (c) Variability is another important aspect. For instance, even though h\ is nearly linear, it fails in characterizing properly the intermediate states T\ and T2. The function hi is preferred here since it allows for the proper discrimination between consecutive states 65 5.5 Normalized E E G signals at different anesthetic depths. Each 1-second epoch is first detrended and then normalized using the RMS amplitude, (a) Awake healthy subject (b) Light R E M sleep (the spikes due to R E M activity were manually removed) (c) General anesthesia (d) Deep anesthetic state (8-waves) (e) Isoelectric state (normalized by the RMS amplitude of the last non-isoelectric epoch). 69 5.6 Selection of the frequency band and wavelet filter for the sampling frequency of 256 S/s. (a) Linearity parameter L. A positive value indicates a concave characteristic, while a negative value indicates convexity. By appropriately choosing the wavelet filter, the characteristic of h can change from concave to convex, (b) Variability parameter 71 5.7 Awake and isoelectric reference PDFs based on / (y A V 71 LIST OF FIGURES xi i i 5.8 Characteristic of the optimal function ^ Y V A V C N S ^ o r t ^ e sampling frequency of 256 S/s. (a) Linearity characteristic, (b) W A V Q N S values for the data sets Ra,Ti,T^, T 3 and R^. Note the larger variability in the intermediate data sets 72 5.9 Effect of ocular artifact de-noising during an awake period, (a) Raw E E G and de-noised E E G . (b) Note how the high frequency information remains unaltered by the de-noising technique, both in terms of amplitude and phase 74 5.10 Characteristic of a 30-second averaging filter and a 2 n d order low pass filter (u>o = 0.02 Hz and C=l). 76 5.11 Time course of the BIS and W A V C N S for Patient #22 and Patient #8 77 5.12 Correlation during steady-state operation between the BIS (v.3.4) and W A V C N S - A total of 13 hours of steady-state data were collected from the arthroscopy study, (a) Correlation density plot (b) Bland-Altman test 78 5.13 Time courses of the BIS (v.3.4) and W A V C N S during induction for both patient groups (each case is synchronized at the LOC). (a) Population showing no reaction to L M A insertion (n=10). (b) Population showing some reaction to the L M A insertion (n=9) 79 5.14 Time delay between the BIS and W A V C N S - (a) Induction (both study groups), (b) Emergence. . . . 80 5.15 Time course of the BIS (v.3.4) and W A V C N S during an episode of burst suppression (patient #4). The corresponding E E G is also plotted for comparison purposes 81 5.16 The HRV signal as a measure of autonomic activity, (a) 3 consecutive QRS complex (the HRV is obtained based on the detection of the R-wave and the measurement of the R-R interval between two consecutive R-waves), (b) Tachogram of a subject performing a relaxation exercise, and a patient reacting to surgical stimulation 83 5.17 Application of the WAV technology to the quantification of autonomic activity using the HRV signal, (a) Representative data sets of HRV signals for different stress levels: Ra: relaxed subject (y-scale: 0.2 s/div), Tj: cold pressure test on volunteer (y-scale: 0.2 s/div), T?: clinical patient undergoing surgical stimulation (y-scale: 0.1 s/div), T3: clinical patient reacting to surgical stress (y-scale: 0.1 s/div). R),: no vagal tone (i.e., no heart rate variability - theoretical state), (b) Result of the best wavelet and frequency band selection. Note that no post-processing filtering is required due to a long analysis epoch 85 5.18 Case report #1 - Female patient, 55 years old, undergoing an hysterectomy procedure. The large gaps in W A V A N S index correspond to periods of time heavily corrupted by electrocautery artifacts (the HRV signal could not be obtained from the E C G signal). Electrocautery also affected the E E G signal and created data loss in the W A V C N S (indicated by vertical bars in the hypnosis index). The 'patient light' events correspond to the anesthesiologist's assessment of the patient during the surgery. Periods of intense electrocautery activity are indicated (it is assumed that they coincide with surgical noxious stimulation) 88 LIST OF FIGURES xiv 5 .19 Case report # 2 - Female patient, 39 years old, undergoing a hysterectomy procedure. The heart rate (top graph) was calculated ad hoc based on the E K G . The blood pressure from the cuff sensor was recorded manually each time this measurement was updated by the anesthesia monitor 8 9 5.20 Case report # 3 - Female patient, 38 years old, undergoing a laparotomy procedure. Due to the nature of the surgery, skin incision was kept to a minimum 9 0 5.21 Experiment design for identifying the sensor dynamics. The E E G signal is composed based on E E G epochs arbitrarily chosen from a database of signals. The input identification signal corresponds to the simulated patient's instantaneous state (direct correspondence with the source E E G signal), while the output identification signal is the W A V C N S 9 2 5.22 Bode plots of the sensor dynamic (analytic and identified models) 9 2 5.23 Accuracy of HQNS and . H A N S for predicting the W A V C N S and W A V A N S time courses calculated for various step changes, (a) The HQNS is able to predict accurately the W A V C N S time course, (b) The HANS is an approximation of the index dynamics. As such, some discrepancies between measured and predicted outputs exist. Also, there is a noticeable index variation in T3. This can be the result of the fact that the patient state was not as stable as in the other states.) 9 3 6.1 Traditional PD model block diagram 9 6 6.2 Example of the 'hysteresis' phenomenon observed between the drug plasma concentration and the corresponding effect, (a) In this example, an output signal is obtained from a first order system, (b) Plotting the output vs. input relationship yields the 'hysteresis' curve described by pharmacologists. By tuning appropriately the rate constant ke0 of a first order model, one can collapse the loop, (c) Under-compensation: the rate constant is too large, (d) Perfect equilibration. Note that the measured vs. predicted characteristic is linear (i.e., in this case, a Hill saturation would not improve further the fit), (e) Over-compensation: the rate constant is too small, which results in an inverted loop. . 9 9 6.3 Effect of a pure time delay, (a) A pure 10 seconds time delay was added to the true system of Figure 6.2. (b) Plotting the output vs. input relationship yields a 'hysteresis' curve similar to that of the previous example (i.e., there is no discriminating feature which can stress out the presence of a delay), (c) Under-compensation. (d) Equilibration. Note, however, that the rate constant is now about 6 times smaller than the true time constant of the system. The first order model therefore mostly captures the time delay dynamics. The model is just a coarse approximation of the true system. Even though the 'hysteresis' loop is reduced, the error between the predicted output and the measured output becomes quite large, (e) Over-compensation 100 6.4 Block diagram of the proposed system oriented approach to PD modeling 101 6.5 PD identification illustrative example (patient #52) . (a) Traditional approach (b) System oriented approach 105 6.6 Block diagram used in stage #3 : the non-linear element is omitted and the sensor filter is used as an input filter 106 LIST OF FIGURES xv 6.7 LTI part of the models, (a) In the traditional approach, the Hill equation is replaced by its static gain, (b) In the system-oriented approach, the sensor dynamics is to be explicitly included 109 6.8 Propofol plasma concentration (bottom plot) and W A V C N S time course (top plot) compared to the outputs obtained from the PD linear models (Patient #65) 110 6.9 Measured vs. predicted W A V C N S time courses for both models (Patient #65). (a) Standard residuals. (b) Prediction errors (assuming an autoregressive model) 110 6.10 Standard residual analysis for model validation, (a) Whiteness test (note that the autocorrelation coefficient at lag 0 is not supposed to be contained within the confidence intervals). Note that the traditional PD model fails this test, indicating the need for a higher degree model, (b) Independence of the residuals with respect to the input signal Cp 113 6.11 Effect of age on the E C B O parameter, (a) Traditional approach (b) New system-oriented approach. . 113 6.12 Comparison between the identified PD parameters from the traditional approach, and the parameters from literature 115 6.13 Error in peak effect (time and value) 117 6.14 Identified W A V C N S VS. Propofol dose-response relationships, (a) Traditional approach (b) New ap-proach. Note that the thick fines represent averages over the corresponding age group 121 6.15 P K P D block diagrams, (a) Traditional approach. Note that the Hill saturation is characterized by both the steepness coefficient 7 and the EC50 parameter, (b) System-oriented approach. Conversely to the traditional approach, the Hill saturation is defined uniquely by its steepness coefficient 7. . . 122 6.16 Linearized P K P D block diagrams. Note that the scaling function was removed since it does not affect the dynamic response of the system, (a) Traditional approach, (b) System-oriented approach. . . . 122 6.17 Frequency response of the PD models derived in Section 6.2.5. (a) Traditional approach (b) System-oriented approach (frequency response also models the sensor dynamics for consistency with the traditional approach) 123 6.18 Comparison of traditional vs. system-oriented frequency response models (all cases averaged). . . . 124 6.19 Comparison of the static dose vs. response relationships obtained for propofol. In their study, Kazama et al. used Gepts P K parameter set (see Table B. l ) . This set was originally derived for infusion administration. As a result, the EC50 found in their study are significantly different than the ones we originally found. However, this difference originates mostly from the choice of the P K set. In order to limit the influence of the P K set, we therefore re-processed all the cases using Gepts P K parameter set. We also added Eo as an identified parameter. We found that Eo is comprised between 0.02 (Gl) and 0.06 (G4). When scaling the effect using the W A V C N S scale, the Eo values translate in a baseline W A V C N S of 98 (Gl) to 94 (G4). These slight differences in the overall awake baseline values were found to be also a function of age 125 7.1 Unstructured uncertainty expressed in a multiplicative framework 128 LIST OF FIGURES xvi 7.2 Quantifying uncertainty, (a) Uncertainty expressed in the frequency domain via Bode plots, (b) Nyquist mapping of the uncertainty at the given frequency u> 129 7.3 Effect of the nominal model Nyquist path with respect to the uncertainty radius: note that nominal models which are close to the uncertainty edges can significantly increase the unstructured uncertainty radius 130 7.4 Relocation of the nominal model Go to the location to reduce the uncertainty radius, (a) ( < ^ M A X ( W ) — < £ > M I N ( W ) ) < 7r (large gain uncertainty relative to phase): in this case the circumscribing circle minimizes the uncertainty disc, (b) ( ^ M A X ( W ) — ^ M I N ( W ) ) < i" (small gain uncertainty relative to phase): in this case, the minimizing uncertainty disc center is located on the center-point of the {Ci , C2} segment, (c) ( ^ M A X ( W ) — V M I N ( ^ ) ) > TT: in this case, the minimizing uncertainty disc center is located on the inner radius of the ring section, directly opposite to the center of the circle which circumscribes the complementary ring section 131 7.5 Linearization of the Hill equation 136 7.6 Intra- and inter-patient variability during thiopental induction, (a) W A V C N S time courses for 2 patients. Each repeat is synchronized with the end of thiopental injection, (b) Identified P K P D gains and P K time delays 137 7.7 Normalized Propofol P K P D uncertainty bounds expressed in the frequency domain. The PK part of the models was derived based on the published P K parameter sets of Schiittler et al. (a) Based on Traditional PD models (b) Based on system-oriented PD models (does not include the sensor dynamics). 138 7.8 Propofol P K P D uncertainty weights for the complete adult population (drug administration: infusion and bolus ; age group: 18-60 yrs old ; W A V C N S range: 80 to 20). (a) Uncertainty weights calculated based on the initial GQ(JU>) nominal model, (b) Uncertainty weights based on an optimized Nyquist path, (c) Near optimum uncertainty weights ||w(ju;)|| (and its realization ||«>(jaj) ||) obtained based on the near optimal realization G Q P ' ( J U ; ) of the nominal model 140 7.9 Nominal model selection in the frequency domain. Note that the rapid gain drop in the high frequency range could not be captured adequately by wnom(ju). Also, some small phase differences between the optimal model and the realization G Q P ' at about w = 1 • 10~3 rad-s - 1 results in 'bumps' which significantly increase the uncertainty gain 141 7.10 Effect of age consideration on P K P D relative uncertainty 142 7.11 Effect of drug administration consideration on P K P D relative uncertainty. Note that the uncertainty weight for bolus administration presents only an academic interest since it is unlikely that a controller be designed to work uniquely in this mode 143 7.12 Effect of operating W A V C N S range reduction on uncertainty 144 LIST OF FIGURES xvii 7.13 Inter- vs. intra-patient uncertainty (thiopental P K P D models). The patient-specific nominal P K P D models used to calculate the uncertainty weights were derived based on the 4 (or 5) P K P D models derived to calculated the uncertainty bounds. In this case, prior knowledge of each patient P K P D characteristics was necessary 147 7.14 Intra-patient upper uncertainty bounds. This bound (in thick dash line) represents the minimum uncertainty weight which must be considered to guarantee stability if a P K P D model can be identified from the induction data (this model is then substituted to the generic nominal model to calculate the controller parameters) 148 7.15 Limiting the PD identification to the sole estimation of the P K time delay can significantly reduce the uncertainty weight in the high frequency region. As expected, there is no reduction in the low frequency band. (Note that the increase in the uncertainty magnitude observed in the [2-10 - 4; 2-10 - 3] rad-s - 1 frequency range may be the result of slight differences in the way the optimization function that calculates the realization of the optimal nominal model was set.) 149 7.16 Selected uncertainty strategy. The corresponding P K P D nominal models, as well as the uncertainty weights, are summarized in the Table C.2 and Table C.5 (models #11 to #14), and Table C.3 and Table C.6 (models #15 to #18). (a) Application: general anesthesia in the peri-operative environment (the controller takes into account the age of the patient and limit the drug administration to infusions only), (b) Application: sedation in the ICU (in addition to accounting for the patient's age and limiting the drug administration to infusions only, the controller is also limited to maintain a sedation level in the 80 to 50 W A V C N S range) 152 7.17 The frequency uncertainty bounds are represented as ring sections in the Nyquist plot. In the clas-sical uncertainty analysis method followed in this chapter, the ring sections represent the system uncertainty. However, when plotting each individual P K P D model, we can notice that the Nyquist uncertainty disk can be significantly reduced. This approach introduces therefore less conservatism than the classical approach. Note that the absence of models in the outer left and right corners of the ring section may reflect the fact that P K P D models are essentially positive systems 153 7.18 Comparison between the uncertainty weight obtained from the classical method (where the uncertainty is first expressed as frequency domain bounds), and the new approach where all the models axe directly expressed into the complex domain 154 8.1 Patient Simulator block diagram 156 8.2 Patient Simulator block - Details level 1 157 8.3 Patient Simulator - Pharmacokinetics block diagram. Depending on the infusion rate, the simulator will either output the plasma concentration corresponding to the bolus or infusion P K models. . . . 157 8.4 Patient Simulator - Pharmacodynamics block diagram 158 8.5 Patient Simulator - WAV Sensor block diagram 158 8.6 W A V C N S noise measurement based on the raw index 159 LIST OF FIGURES xviii 8.7 Measurement noise frequency characteristic 159 8.8 Measurement and modeling noise frequency characteristic 160 8.9 Identification of the noise characteristics. Two IIR filters were derived to capture the spectral nature of the measurement and modeling noises 160 8.10 Patient Simulator - Noise Simulator block diagram. The drug effect E calculated in the Pharmacodynamics block is used to determine which noise profile that needs to be applied 161 8.11 Measured vs. predicted W A V C N S time courses (LMA patient #065). The noise filter T^oiee was switched from 90 s to 170 s. The noise filter T^ o i s e was used during the rest of the time. Note that only 5 min 20 sec of data were collected in this case. As a result of a human operator error, no inhalational anesthetic was provided to the patient after the L M A was placed. The anesthesia record mentions that the patient woke up 7 min 35 seconds after the start of the surgery. The predicted W A V C N S time course does indeed indicate that the patient was above 80 by that time 162 8.12 Controller test protocol 163 8.13 Close loop system with anti-windup and pre-filter 164 8.14 Infusion rates and W A V C N S time courses in the 3 test patients using the classical PID loop shaping design method (first design). Note that patients presenting an aggressive dynamic behavior are unstable. 165 8.15 Infusion rates and W A V C N S time courses in the 3 test patients using the classical PID loop shaping design method (second design). Note the poorly damped oscillations in the high gain patient. This controller is shown to be mathematically unstable 166 8.16 Close loop sensitivity analysis using the augmented nominal model GQpt(jui). The second design is unstable since there exists a frequency range where the complementary sensitivity is greater than the inverse of the uncertainty weight 167 8.17 Infusion rates and W A V C N S time courses in the 3 test patients using the classical PID loop shaping design method (third design). Performances are strongly reduced, in particular for patients presenting a low P K P D gain 168 8.18 Mixed-sensitivity minimization 169 8.19 Hoo controller implementation with anti-windup 170 8.20 Initial and optimized Hoo designs, (a) Closed-loop complementary sensitivity, (b) Control action profile in response to a setpoint change of 1 W A V C N S unit 171 8.21 Close loop sensitivity analysis 172 8.22 Infusion rates and W A V C N S time courses in the 3 test patients. The HQO design robustified significantly the system which is now stable for the whole adult population. As a result, patients presenting low gains are not well compensated 173 . 8.23 Frequency response of the Hoo and PID controllers. Note that the PID has a faster roll-off in high frequencies, as well as a higher gain in low frequencies due to its integrator 173 LIST OF FIGURES xix 8.24 Infusion rates and W A V C N S time courses in the 3 test patients. Accounting for age does result in a significant increase in performance 175 8.25 Infusion rates and W A V C N S time courses in the 3 test patients. The on-line identification of the P K time delay reduces the system uncertainty and allows for a larger control bandwidth. Performances are improved at the expense of an increased complexity in the controller design 176 8.26 Settling time definition. In the standard definition, the settling time is defined as the time when the system settles within 5% of the setpoint. In cases with overshoot higher than 5%, this may result in much longer settling time values. For instance, in this example, one system has a significantly longer settling time. Therefore, if the overshoot is less of a concern, it may be indicated to use the setpoint crossing as settling time, in which case the two responses plotted are almost equivalent 178 8.27 Smith Predictor control structure 181 8.28 Infusion rates and W A V C N S time courses in the 3 test patients. The use of a Smith Predictor allows for a faster settling time at the expense of a larger overshoot. Note that, in this simulation, the pre-filter was removed to reduce the settling time during setpoint change 182 8.29 Sensitivity peaks of the initial Hoo age targeted design, and the design involving an increased controller gain 183 8.30 Effect of measurement/modeling noise on the closed-loop system response during saturation, (a) The age targeted design with increased gain keeps on administrating propofol, which delays the system response, (b) Using a modified anti-windup scheme where the low frequency pole integrator resets for negative values instead of 0 can be a solution 184 8.31 Summary of the controllers performance in terms of settling time and overshoot 186 A . l Context sensitive half times of fentanyl, sufentanil, alfentanil and remifentanil (from [21]) 211 D . l Comparison between the Movers and Non-Movers groups 232 D.2 Average time courses of W A V C N S and A W A V C N S for both groups. A W A V C N S is the rate of change of the W A V C N S index. Note the significant difference between the two groups before the L M A insertion. 236 D.3 Correlation between W A V C N S value 40 seconds after LOC and the administered propofol dose (note that the study population was extended to patients receiving propofol doses outside the protocol range) 237 D. 4 Identification procedure for 4 different cases. The thiopental blood plasma concentration Cp is esti-mated through the 3-compartment PK model proposed by Stanski et al. This model is weight and age dependent, (a) Successful identification, (b) Successful identification despite limited fascicula-tion. (c) Identification unsuccessful due incomplete data (the index does not reach its peak value), (d) Unsuccessful identification due to a large fasciculation disturbance 239 E . l Survey 1 243 E . 2 Survey II 244 List of Acronyms A C R O N Y M D E S C R I P T I O N A C L Anterior Cruciate Ligament A S A American Society of Anesthesiologists A S S R Auditory Steady State Response BIS Bispectral Index Scale B S R Burst Suppression Ratio CNS Central Nervous System D O A Depth of Anesthesia D W T Discrete Wavelet Transform E C G Electrocardiogram E C T Electro-Convulsive Therapy E E G Electroencephalogram E M G Electromyography E O G Electrooculogram E S U Electro-surgical unit E R P Evoked Response Potential F D A Food and Drug Administration F IR Finite Impulse Response H R V Heart Rate Variability I C U Intensive Care Unit IIR Infinite Impulse Response L M A Laryngeal Mask Airway L O C Loss Of Count M A C Minimum Alveolar Concentration X X LIST OF ACRONYMS A C R O N Y M D E S C R I P T I O N Median Frequency Mid-Latency Auditory Evoked Response Neuromuscular Blocking Agent Ocular Artifact Ocular Artifacts Operating Room Pharmacodynamics Probability Density Function Positron Emission Tomography Proportional-Integrative-Derivative control law Pharmacokinetics Pharmacokinetic - Pharmacodynamic Patient Simulator Rapid Eye Movement Root Mean Square Respiratory Sinus Arrhythmia Redundant Wavelet Transform Semilinear Canonical Correlation Spectral Edge Frequency Somatosensory Evoked Potential Smith Predictor Short Time Fourier Transform Stationary Wavelet Transform The University of British Columbia The University of British Columbia Hospital Vancouver General Hospital Wavelet-based Anesthesia Value Wavelet-based Anesthesia Value for Central Nervous System monitoring Wavelet-based Anesthesia Value for Autonomic Nervous System monitoring Wavelet Transform M E F M L A E R N M B O A OAs O R P D P D F P E T P I D P K P K P D P S R E M R M S R S A R W T s e c S E F S E P S P S T F T S W T U B C U B C H V G H W A V W A V C N S W A V C N S W T Acknowledgements As I am concluding this thesis, I would like to express my deepest gratitude to my supervisor, Professor Guy Dumont, for his valuable support, advice, and guidance throughout the course of this research. Being part of the multidisciplinary research team has given me the unique opportunity to familiarize myself to an exiting new field of research. I am also particularly grateful to Dr. Mihai Huzmezan, my co-supervisor, for his constant encouragement and unfaltering support. His optimism and pragmatism have often smoothed over many difficulties and moments of doubt. I am also deeply indebted to Dr. Craig Ries, my co-supervisor, for sharing with me his expertise in anesthesia, and for opening the many doors that have allowed this work to flourish. I am looking forward to our continuing collaboration. I would also like to thank Ms. Tatjana Zikov, who shared with me her scientific expertise in signal processing, and her commitment to anesthesia research. I will always be grateful for the long hours we spent together reviewing data, and sharing in the excitement of new discoveries. Thanks are extended to Dr. Kris Lundqvist, Dr. Henrik Huttunen, Dr. Jason Waechter, and Ms. Lou-Ann Mendoza from the U B C departments of Anesthesia and Electrical and Computer Engineering, for helping prepare and co-ordinate the clinical studies, which provided the basis for most of the results and contributions presented in this thesis. I would also like to acknowledge Cleveland Medical Devices Inc. (Cleveland, O H , U S A ) , and in partic-ular Dr. Mohammad Modarreszadeh and M r . Robert Schmidt, for their interest in carrying forward the development of the NeuroSENSE™ anesthesia monitor, which integrates the sensing technology developed in part in this thesis. Finally, I would like to thank all my friends here in Vancouver, for their support and friendship, and for having made my stay in Vancouver memorable. Stephane Bibian July, 2006 xxii This thesis is dedicated to the memory of my grandfather, Albert Vande Kerckhove, a truly righteous and altruistic man, who will always live in our memories. xxiii PART A: The Anesthesia Framework Chapter 1 introduces the readers to the terminology and concepts in use in the field of clinical anesthesia. A section dedicated to the most commonly administered drugs is presented, as well as information regarding drug delivery and standard of care in monitoring. The second chapter addresses the problem of Depth of Anesthesia (DOA) measurement. This particular topic has received significant attention from both the anesthesia and engineering communities. Even as of today, there is still no definite consensus as to how this measurement should be performed. While there are commercially available monitors of the anesthetic state, current practice does not mandate these monitors as a standard of care. These monitors are often criticized by practitioners, and are the source of debates and controversies. Knowledge of the process is vital in any control design venture. Since online identification of the process is not practical for this particular application, it is necessary to establish proper models relating infusion rates and observed effects, as well as the uncertainty range. Chapter 3 investigates the pharmacokinetic and pharmacodynamic concepts used in the clinical field to model the drug-response relationship. Finally, Chapter 4 lays out the foundation for an anesthesia controller. Prior attempts are reviewed. Based on the experience gained from these attempts, as well as the concepts presented in the previous chapters, we propose a control framework suitable for this particular application. The main requirements and challenges are outlined. 1 Chapter 1 Clinical Anesthesia: Terminology, Concepts, and Issues Industrialized societies have enjoyed the benefits brought by clinical surgery and modern medicine for over a century. Undergoing a surgical procedure has become nowadays a rather common event in one's life. However, this remarkable progress has only been made possible through the development of modern anesthesia practice since its discovery in the mid 19 t h century. While the aim of this chapter is not to present a thorough review of the practice of clinical anesthesia, we feel that it is necessary to introduce our readers with an overview of the concepts now in use in this medical field. The short history presented in the first section will emphasize the tremendous impact that the discovery of anesthesia had in medicine, as well as the numerous developments that have taken place ever since. The concepts and definition of anesthesia, as well as the different anesthetic and opioid agents used in today's practice, will then be discussed. Finally, we will conclude this chapter by briefly reviewing the different devices and sensors most commonly found in today's practice. 1.1 Clinical Anesthesia: A Key Specialty in Critical Care Before the advent of anesthesia, surgical procedures demanded extremely fast execution. Early regional techniques such as nerve compression or the application of cold were used to provide a slight relief from the pain. Decreased cerebral perfusion obtained by compressing the carotid artery was also used to render patients unconscious. Quite clearly, the discovery of inhalation gases that could provoke a state of anesthesia, and thus make invasive surgeries possible, was a major event in the development of modern medicine. Numerous books and textbooks give the historical perspective to clinical anesthesia. We invite interested readers to refer to [1], [2], and [3] for a more in depth look at this medical specialty. 2 CHAPTER 1. CLINICAL ANESTHESIA: TERMINOLOGY, CONCEPTS, AND ISSUES 3 1.1.1 The Early Days While the alteration of the senses using drugs (opium, laudanum, mandragora, etc.) and alcohol was well known since antiquity, it's only in 1840 that the idea of using inhalation gases and vapors, otherwise reserved to entertain the public in village fairs, was put forward by Hickman. The first recorded anesthesia procedure was performed by Crawford Williamson Long using diethyl ether. Having often inhaled ether, he realized that this gas had the property of rendering him particularly insensitive to painful knocks and bruises. He concluded that this could be applicable to surgery, and tried out his theory in March 1842. While he continued using ether for surgeries, he did not publish his results until later. In 1844, Gardner Colton and Horace Wells performed the first dental surgery using nitrous oxide as an anesthetic. However, they failed in convincing the medical community, and it is only later, in October 1846, that Will iam Morton succeeded in demonstrating to his colleagues that ether could be used to deprive patients of their sense during surgery. The success of his 'etherization' technique spread quickly throughout the civilized world. The term 'anesthesia' (lack of esthesia, i.e. sense) was later proposed to Morton by Oliver Wendell Holmes to describe this new phenomenon. While ether has been advocated by many in the early days, other agents such as chloroform, nitrous oxide, and ethyl chloride were investigated. Nowadays, only nitrous oxide is still in use as a supplement to the modern agents developed in the 1960s and 1970s. The development of drugs used for anesthesia procedures will receive a more thorough review in Section 1.3. Parallel to the development of anesthetic agents, the 19 t h century witnessed the invention of many inhalation apparatuses and techniques for the administration of anesthetics. For instance, there was a need for an artificial ventilator, since a major side effect of many anesthetics and opioids is ventilatory depression. The first ventilator was developed in 1896 by Northrup, who had shown that patients who had opium poisoning could be maintained alive through artificial ventilation, using an endotracheal tube and bellows. But only in the 1950s, after the introduction of muscle relaxants and polio epidemics, did the use of ventilators became common in the operating room. Due to the increasing complexity of the administration and management of anesthesia, it became clear that anesthesia and surgery were two complementary specialties that demanded their own practitioners. But it was only in 1935 that the first diploma of anesthesia was offered, and in 1936 that the American Society of Anesthesiologists (ASA) was founded. 1.1.2 Risks and Outcome in Anesthesia 1.1.2.1 Mortality rate Development of clinical anesthesia since the 19 t h century has resulted mostly from concerns about patients' safety. The first recorded death attributed directly to anesthesia was that of Hannah Greener in 1848, after CHAPTER 1. CLINICAL ANESTHESIA: TERMINOLOGY, CONCEPTS, AND ISSUES 4 inhaling chloroform for the removal of an ingrown toenail. The question whether the patient's death resulted from an inadequate anesthesia management or from undesirable side effects due to chloroform was hotly debated. Anesthesia and the use of chloroform being rather new, this case received considerable attention from both the medical community and the public, and was particularly well-publicized. Chloroform was by then a newly introduced drug. Less flammable than ether, chloroform offered an attractive alternative to surgeons. But after more deaths were reported throughout the world, the safety of chloroform as an anesthetic drug was questioned. Even then, it was only in the 1920s that chloroform was replaced by new, more potent, and safer gases. The cause of Hannah Greener's death remains unclear, and 150 years later this case still generates interest [4]. B y the late 19 t h century, the incidence of death due to anesthesia was less than 0.1% [5]. Nowadays, clinical anesthesia is probably one of the safest components of any surgical operation. A 1986 survey [6] revealed that the overall death rate attributable directly to anesthetic drugs was 1:185,056 (see Table 1.1). W i t h approximately 28 million patients undergoing anesthesia and surgery in the United States [7], it is estimated that about 150 patients die each year in the U.S .A . from complications due to anesthesia. This very low mortality rate can be attributed mainly to the following three aspects of the clinical practice: i . First, anesthesiologists select an appropriate combination of drugs and drug dosage according to the patient's age, weight, co-morbid disease, and the type and duration of the operation. In standard practice, anesthesiologists often use several drugs in order to reach a state of balanced anesthesia (refer to Section 1.2.3), thus limiting the potential lethal side effects of each drug. i i . The second aspect concerns the equipment that monitors patients' vital signs and eventually warns the practitioner of possible complications. Modern equipment is fairly sophisticated and includes standard devices such as mass spectrometers, capnographs, pulse oximeters, heart rate and blood pressure monitors, etc. i i i . Finally, education has had a key role in making anesthesia a particularly safe and reliable procedure. Postgraduate training programs in the specialty of anesthesiology are offered by every major medical schools. The content of these programs are usually supervised by a centralized specialty college in most countries. Also, this medical specialty benefits from the publication of numerous clinical research journals (Anesthesia and Analgesia, Anesthesiology, British Journal of Anaesthesia, e tc . . . ) . CHAPTER 1. CLINICAL ANESTHESIA: TERMINOLOGY, CONCEPTS, AND ISSUES 5 S T U D Y D A T E T O T A L M O R T A L I T Y C A S E S R A T E Beecher 1948-52 599,548 1:1,560 a Clifton 1952-62 205,640 1:3,955 a Harrison 1967-76 240,483 1:4,537 b Hatton 1977 190,389 1:2,885 a Lunn 1979 1,147,362 1:6,789 e Eichhorn 1976-85 757,000 1:151,400 d Eichhorn 1985-88 244,000 0 d C E P O D 1986 486,000 1:185,056 e C E P O D 1986 486,000 1:1,185 c a: All operative cases considered in calculation b: Cases included if death occurred in less than 24 hours c: Cases included if some contribution by anesthetic d: Only ASA physical status I and II patients included e: Only deaths directly attributable to anesthetic included CEPOD: Confidential Enquiry into Postoperative Deaths Table 1.1: Contemporary anesthetic mortality rates (adapted from Brown [5]) The method used for estimating the mortality rate depends greatly on the methodology used by each researcher. While deaths during surgery were easily attributed to the anesthetic drugs in the early days, such cases are now thoroughly investigated and documented. Also, the patient physical status is now accounted for. The outcome of any surgery depends indeed whether the patient is healthy (ASA Physical Status I) or in critical condition (ASA V ) . The mortality rate is obviously much higher in A S A V patients (about 10% [5]) who undergo emergency surgery following major trauma than for A S A I or II patients (0.08% and 0.27% respectively) who undergo elective surgery. The A S A (American Society of Anesthesiology) rating was introduced in 1941. Human error is probably the most common cause of death (hypoxic gas mixture, airway obstruction, errors in drug administration, lapses in vigilance, e t c . . ) [8]. According to a 1987 study [9], 75% of anesthetic related deaths are due to the anesthetists' failure to apply life saving knowledge, while only 1.7% of cases involve equipment failure. Also, where the most common reason cited leading to such events used to be overdosing (1960-1969), it is nowadays the inadequate preoperative preparation and patient assessment that leads to errors in anesthesia management. Finally the actual trend of increased efficiency in the operating room (so called production pressures) is an additional factor of anesthesiologists' fatigue that might provoke misjudgement, resulting in reduced patients' safety. 1.1.2.2 In t raopera t ive awareness While patient's safety is still an important issue, new concerns arose from the use of neuromuscular block-ing agents (NMBs) in the early 1940s. These drugs, used to block muscle movement and thus reduce muscle tone to facilitate surgery, have the unpleasant consequence of obtunding an important sign of light CHAPTER 1. CLINICAL ANESTHESIA: TERMINOLOGY, CONCEPTS, AND ISSUES 6 anesthesia, i.e. movement. It is therefore possible for a patient to remain immobile while being aware of his/her surroundings, and experiencing the trauma caused by the surgery [10]. While not necessarily lethal, awareness with pain during surgery often results in severe psychological consequences. These cases are fortunately extremely rare (about 1 in 10,000 [11]) and result principally from faulty equipment or human error. However, limited intraoperative awareness without the presence of pain is more common, mostly when patients are maintained in a shallower depth of anesthesia or during emergencies. A number of such cases have been reported in the literature. Very recently, two extensive studies have been conducted to determine the incidence of intra-operative awareness: - S A F E 1 Trial: Sandin et al. [12] have undertaken an extensive survey of intra-operative awareness cases in Sweden, spanning 2 major hospitals and 2 years of clinical activity (from 1997 to 1999). 11,785 patients were interviewed on three separate occasions (immediately after the operation, 1 to 3 days after surgery, and then once again 7 to 14 days after). A total of 18 cases (0.15%) of awareness were reported, and 14 of these cases occurred in surgeries involving N M B s . 11 of these 14 paralyzed patients experienced either pain, anxiety, or delayed neurotic symptoms, while the 4 non-paralyzed patients did not suffer during their period of wakefulness. This study also revealed that patients might recall intra-operative events only after some time after the surgery. For instance, this study revealed that only 11 cases would have been identified if the patients had only been interviewed once immediately following the operation. - A I M Trial: A similar study was undertaken in the United States [13]. 19,576 patients over 7 different sites have been enrolled so far in this study (the study objective is 50,000). The incidence of intra-operative awareness has so far been found to be 1 to 2 cases per 1000 patients. Considering that there are 20 millions of anesthetic procedures carried out each year, it is estimated that about 100 patients suffer from intra-operative awareness every weekday worldwide. It is estimated that up to 50% of the patients who experience intra-operative awareness will develop Post-Traumatic Stress Disorders (PTSD) in the following year. 1.2 Modern Concepts In this section, we will discuss the practice of modern clinical anesthesia by introducing the concepts of the anesthesia triad and balanced anesthesia. These concepts play a key role in defining and understanding the proposed research project. 1.2.1 T h e Role of the Anesthesiologist In lightly anesthetized patients, surgical stimulation leads to movement, cardiorespiratory changes (i.e., rise in blood pressure, heart rate and respiration), and hormonal responses (release of adrenaline, Cortisol, CHAPTER 1. CLINICAL ANESTHESIA: TERMINOLOGY, CONCEPTS, AND ISSUES 7 etc.). During "normal" stress situations of course (i.e., no surgery, no anesthesia), these cardiorespiratory and hormonal responses act as the body's self-defense system in that they assist the subject to fight or flee. During surgery, however, "stress responses" can be dangerous, and even life-threatening, as they can provoke strokes and heart attacks in susceptible individuals. Intra-operative awareness is another potential risk in under-anesthetized patients. While not fatal, it can provoke postoperative post-traumatic stress disorders [14, 15]. In order to blunt the effect of surgical stimulation, anesthesiologists use a combination of drugs to block sensation. However, the very mechanism of action of these drugs make them particularly dangerous, as they deprive the central nervous system from the information necessary to control normal body functions (i.e., gag reflex, respiration, cardiac rhythm and blood pressure). A n overdose may then stop a patient's breathing and might even provoke a cardiovascular collapse. Overdoses are usually associated with a lack of balance between the anesthetic regimen and the patient's pharmacological needs. When there is no surgical stimulation, the patients' needs are low and a small amount of drug may be sufficient to make them comfortable. However, during noxious stimuli (i.e., stimuli associated with transmission of nerve pain signals), drug titration needs to be increased to limit the effect of surgery. As a result, a common side effect is the depression of the cardiorespiratory system when surgical stimulation suddenly disappears. Therefore, anesthesiologists try to keep a balance between the toxicity of anesthetic drugs and the noxious stimulation of surgery. 1.2.2 Functional Components of Clinical Anesthesia Although its scientific definition is still open to debate, anesthesia has been described as a state of "drug-induced unconsciousness, [where] the patient neither perceives nor recalls noxious stimuli" [16]. This func-tional definition proposed by Prys-Roberts in 1987 limits the term of anesthesia to an absence of both conscious awareness and memory formation (i.e., hypnosis and amnesia). However, the role of anesthesiologists goes beyond provoking a mere hypnotic state. They also ensure that autonomic reflexes involving the sympathetic and parasympathetic nervous system (to provide car-diorespiratory control) are not sensitive to surgical stress. This is achieved by inducing a state of analgesia. Opioid drugs are typically used to achieve this endpoint. When both the hypnotic and analgesic states are adequately reached, patients do not usually exhibit purposeful movement to surgical stress. However, intra-abdominal surgeries require the blockade of reflex muscle activity in the abdominal wall in order to permit surgical exposure. To attain such a state of paralysis, it is necessary to use neuro-muscular blocking agents (NMBs) . It is important to note that these drugs act peripherally at the level of the synaptic link between the nerve and the muscle, and not centrally in the brain or the spinal cord. To summarize, it is common in the literature to consider that the state of general clinical anesthesia results from the combination of three functional components, that is, hypnosis, analgesia and immobility. CHAPTER 1. CLINICAL ANESTHESIA: TERMINOLOGY, CONCEPTS, AND ISSUES 8 While this anesthesia triad [17] concept is somewhat simplistic, most engineering oriented authors contribut-ing in this field are setting their work within this functional framework since each individual component can be accomplished by separate classes of intravenous drugs [18]. Note, however, that inhaled anesthetics have minor analgesic and relaxant properties and large doses of opioids can depress cognitive function. Not surprisingly then, anesthetics and opioids enhance or potentiate each other's action in a synergistic fashion. Conversely, interactions between N M B s and opioids, or N M B s and intravenous hypnotics are minimal. 1.2.3 Modern Balanced Anesthesia As stated earlier, all anesthetics are by definition hypnotics. They first act at the level of cognition (cortex) by rendering patients unconscious. However, with increasing doses of anesthetic drugs, it is possible to go beyond hypnosis and blunt further the response to noxious stimuli. This is particularly true for inhalational anesthetics. T i l l the 1940s, it was common for anesthetists to use a single agent at high concentration to achieve adequate anesthesia. Unfortunately, using higher doses usually results in stronger side effects during surgery (ventilation depression, cardiac arrhythmia, etc.) and also during recovery (nausea, vomiting, etc.). To alleviate undesirable side effects George Crile advocated in 1911 the use of local analgesics as a complement to light general anesthesia. In 1926 the term balanced anesthesia was introduced by John Lundy to describe a combination of agents that would achieve adequate anesthesia. This concept can be well-understood when considering that analgesia and areflexia can be achieved separately using drugs such as opioids and N M B s . These drugs are not hypnotics in the sense that they do not provoke uncon-sciousness. Nowadays, the focus in clinical anesthesia is to achieve an adequate balanced anesthetic state using a combination of hypnotics (inhalational/intravenous anesthetics), opioids and N M B s . This tech-nique has the advantage that much lower concentrations of these drugs are then needed as they are used concurrently. Hence, side effects, recovery time, and post operative nausea and vomiting are considerably reduced. Balanced anesthesia is now the standard in the management and conduct of clinical anesthesia. 1.3 An Extensive Pharmacopoeia A combination of anesthetics and opioids, with or without N M B s , are administered together to create the state of general anesthesia. These drugs, even when taken within the same family, have different properties. Since they provide the actuators through which the patient's state can be regulated (i.e., allowing the control of the anesthetic state), it is necessary to provide control engineers with some knowledge of the mechanism of action of the most commonly used drugs. 1.3.1 Anesthetics (i.e. hypnotics) Inhaled Anesthetics W i t h the advent of fluorine technology in the 1940s, new inhaled anesthetics were developed. Compared to ether and chloroform, fluorine compounds have lower blood solubility (thus ensur-CHAPTER 1. CLINICAL ANESTHESIA: TERMINOLOGY, CONCEPTS, AND ISSUES 9 ing rapid induction and recovery), lower toxicity, are less irritating to the airway, and are not flammable. Nowadays three agents are commonly used with or without nitrous oxide: isoflurane, desflurane and sevoflu-rane. A l l these agents provoke a decrease in the mean arterial blood pressure when administered to healthy subjects. A major advantage of inhaled anesthetics is that the drug uptake in the arterial blood stream can be precisely titrated by measuring the difference between the inspired and expired concentrations. This measurement is done in real-time using a device such as a mass spectrometer. A t steady state, the expired concentration correlates with the brain concentration. As a result, inhaled gases are used extensively to maintain a desired depth of anesthesia. This discussion would not be complete without mentioning the concept of Minimum Alveolar Concen-tration ( M A C ) . A value of 1 M A C is the minimum alveolar concentration of an inhaled anesthetic that prevents purposeful movement in response to a noxious stimuli in 50% of patients. The M A C value is thus used to compare dosage between inhaled anesthetics. Intravenous Anesthetics Intravenous anesthetics can be classified into 5 families: Barbiturates (thiopen-tal), Benzodiazepines (midazolam, diazepam, lorazepam), Phencyclidines (ketamine), Carboxylated imida-zoles (etomidate), and Isopropylphenols (propofol). Compared with volatile agents, intravenous anesthetics (besides ketamine) do not provide analgesic effects - hence, they are defined as hypnotic rather than anesthetic drugs. However, opioids and intravenous anesthetics, when used in combination, are strongly synergistic, both in terms of hypnosis and analgesia. Propofol was introduced in the early 1990s and has become the intravenous drug of choice in anesthesia. Two particular characteristic of propofol are its fast redistribution and its metabolism. As a result it can be easily used in infusion schemes as it provides very fast emergence, without cumulative effect. Inhalational Anesthetics vs. Intravenous Agents Inhalational anesthetics are considered by many practitioners as near ideal anesthetics as they have both an hypnotic and an analgesic effect. This explains why past closed-loop anesthesia attempts used inhalational anesthetics as the sole actuator. Combined with the fact that the lung partial pressures of inhaled anesthetics is closely related to the vapour concentration in the brain, the control problem is significantly simplified since additional states are measurable. As opposed to inhaled anesthetics, the brain concentration of an intravenous drug cannot be easily measured. As a result, the titration of these drugs is more difficult as the anesthesiologist does not have any feedback on how much of the administered drug has been metabolized or stored in inactive tissues. In the majority of cases, intravenous agents are given as large boluses (large doses given over a very short amount of time, e.g., < 1 min) for the induction of anesthesia, while maintenance is accomplished with inhaled vapors (but more and more often with intravenous infusions). However, since intravenous agents are more specific than inhaled anesthetics, they give more flexibility in separately controlling the functional components of anesthesia. CHAPTER 1. CLINICAL ANESTHESIA: TERMINOLOGY, CONCEPTS, AND ISSUES 10 1.3.2 Opioids Opioids are unique in the sense that they provoke analgesia without loss of touch, temperature sensation and consciousness, when administered in small doses. They act as agonists at specific receptors within the Central Nervous System (CNS). Their principal effect may be the inhibition of neurotransmitter release, resulting in a significant analgesic effect. Contrary to most anesthetics, opioids do not directly depress the heart, and are thus particularly suitable for cardiac anesthesia. Opioids can produce unconsciousness when used in very large doses. This observation has led some authors to believe that opioids should be considered to be anesthetics, as they fit Prys-Roberts' definition [19]. However, the state of unconsciousness brought by opioids is not reliable. It has been shown for instance that they cannot fully replace inhaled vapours to provoke an adequate state of hypnosis. However, their use can reduce the requirements of inhaled anesthetics by up to 50% [20]. Also, the sedative effect of opioids is opposed by the presence of acute pain. Hence, even though patients in severe pain receive very large amount of opioids, they can remain aware. In current practice, therefore, opioids are almost always supplemented by other anesthetics. Five opioid compounds are used in clinical anesthesia: morphine, hydromorphone, fentanyl, sufentanil, and remifentanil. While they all have similar effects, their kinetic characteristics differ tremendously due to large differences in their lipid-solubility. Of particular interest is remifentanil, a new agent introduced in the mid 1990s. The potency of remifentanil is twice that of fentanyl and its effect-site equilibration time is slightly less than that of alfentanil (about 1.1 min). Remifentanil differs from all other opioids due to a unique structure where an ester linkage makes it susceptible to hydrolysis. This property results in a rapid degradation to inactive metabolites [21]. The main characteristics of remifentanil then are: rapid onset, brevity of action, noncumulative effects in inactive tissues and rapid recovery after termination of the infusion. Remifentanil is used mostly to provide the analgesic component of general anesthesia. Its brevity of action allows patients to recover rapidly from undesirable opioid-induced side-effects such as ventilatory depression. 1.3.3 Neuromuscular Blocking Agents (NMBs) N M B s act locally at the level of the neuromuscular junction by interrupting the transmission of nerve impulses. Their principal use is to produce skeletal muscle relaxation to facilitate intubation of the trachea and to provide optimal surgical conditions. N M B s do not have any analgesics or hypnotic properties. They also do not interact in a clinically significant way with opioids and anesthetics. Succinylcholine is used whenever a rapid onset and short duration of action is needed. Cisatracurium, vecuronium, mivacurium and rocuronium are used when a longer effect is desired. CHAPTER 1. CLINICAL ANESTHESIA: TERMINOLOGY, CONCEPTS, AND ISSUES 11 1.3.4 Drugs: Action, Effect and Interaction The development of safer and more potent intravenous agents with faster onset of effect and, in certain cases, shorter duration of action, has greatly impacted the practice of anesthesia. Nowadays, small drug quantities used in combination can produce a balanced state of anesthesia while minimizing side effects. In North America, inhalational gases are used to maintain the background of anesthesia. However, intra-venous agents are increasingly employed for maintenance. The use of intravenous agents is geared towards facilitating intubation, compensate for undesirable changes in patients' state and also to anticipate painful surgical stimuli. In this realm, the short acting characteristic of intravenous agents such as remifentanil and propofol indicates that these drugs should best be used in infusion regimens, since their administration as boluses often result in too strong effects for too short periods of time. The inability of measuring plasma drug concentration of intravenous agents raises questions. Currently this affects the ability of the anesthesiologist to set precise rates of infusion. Infusion regimens published in medical journals are prone to error due to model uncertainty. Hence, the resulting titration might not correspond to the patients' real needs. In the context of closed-loop control, and when using intravenous drug as "actuators", it is necessary to account for both patient variability and drug synergism. Patient variability results from differences in the way the drug distributes and is eliminated from the body which is further influenced by cardiac output, an individual's age, lean body mass, renal and liver functions, etc. Genetic differences and enzyme activity might also alter the sensitivity to the drug. For example, while some patients might be hyporeactive (e.g., tolerance due to addiction), others may be hypersensitive. When using different drugs in combination, significant interactions can be observed. A n additive effect signifies that a second drug taken concurrently with a first will produce an effect equal to the superposition of their effects (e.g. the anesthetic effects of two inhaled anesthetics are additive [22]), whereas a synergistic effect means that the resulting effect is greater than what could be expected from superposition. Drug synergism often appears when using hypnotics in combination with opioids. In some particular cases, drugs can also be antagonistic, in which case they counter-act each other when administered concurrently (e.g., an opioid antagonist such as naloxone provokes a complete reversal of the effect of opioids). From a control point of view, such drug interactions tend to generate significant cross-coupling in the multivariable system that models the patient. Only very few models of such couplings have been discussed in the literature; see [23, 24] and [25]. These models are mainly mathematical expressions that describe drug interaction in steady state. As such they do capture the static non-linearity of drug interactions, but fail to characterize the dynamics associated with state transitions during transitory events. CHAPTER 1. CLINICAL ANESTHESIA: TERMINOLOGY, CONCEPTS, AND ISSUES 12 1.4 Conduct of Anesthesia In most cases, we can divide the anesthesia procedure into 3 distinct phases: induction, maintenance and emergence. I n d u c t i o n Changing a patients consciousness from an alert to an anesthetized state involves a transition phase called induction. A n intravenous induction of anesthesia is typically used in adults since inhalational inductions are slower and can be associated with an intermediate "stage of excitation" with vomiting and/or spasm of the vocal cords. During intravenous induction with a hypnotic agent, the intermediate stage of excitation is generally not seen and airway, respiratory and cardiovascular reflexes are rapidly depressed with the sudden onset of unconsciousness. For example, the tongue relaxes to obstruct the airway and spontaneous breathing stops. As a result, the upper airway must be instrumented and ventilation must be, at least initially, manually or mechanically controlled. Airway instrumentation methods are invasive (i.e., like surgery) and generally determine the initial need for opioids and N M B s , which in turn may further depress airway and ventilatory control. There is no clinically useful brain monitor as yet that can assist in the initial titration of intravenous anesthetic drug dosage, so cardiorespiratory stability is used to infer an absence of consciousness. The induction of anesthesia is then a critical phase, as the first priority is to maintain the patient's airway. For example, major surgery generally requires the patient's airway to be secured by an endotracheal tube - a tube and air-filled cuff that is placed through the vocal cords - which then allows positive pressure and mechanical ventilation. A n additional advantage of the endotracheal tube is that its intra-tracheal cuff prevents the pulmonary aspiration of gastric contents in the case of either gastric-esophagal reflux or gastric vomiting with food ingestion prior to emergency surgery. As discussed earlier, the insertion of the endotracheal tube is particularly painful and stimulating (e.g., skin incision) as it passes through the vocal cords and thus can provoke strong gagging reflexes. To facilitate the insertion of the tube, it is common to first use an opioid such as fentanyl to blunt the pain, followed a few minutes later by a large bolus of propofol. The administration of such a bolus will rapidly provoke a deep level of hypnosis. To further blunt any motor reflex, a N M B is generally administered. A n alternative to the endotracheal tube is the Laryngeal Mask Airway - a tube and air-filled cuff that sits above the voice-box. Less invasive than the endotracheal tube, the laryngeal mask airway usually does not necessitate the use of N M B s for its insertion since it does not pass through the vocal cords. However, it is less secure in case of gastric reflux and does not therefore prevent pulmonary aspiration. As a result, laryngeal mask airways are used essentially on healthy patients undergoing minor elective surgery on empty stomachs. The induction of anesthesia is usually a very short event lasting only a few minutes. CHAPTER 1. CLINICAL ANESTHESIA: TERMINOLOGY, CONCEPTS, AND ISSUES 13 M a i n t e n a n c e The effect of the propofol induction dose being very short (due to re-distribution), it usually wears off once the tube or the laryngeal mask airway is inserted. Hence, it is necessary to administer more anesthetic to maintain an adequate depth of hypnosis. The maintenance of anesthesia is often obtained via an inhalational agent set around 1 M A C . Prior to surgical incision or any other stimulating surgical event, it is customary to reinforce the analgesic component of anesthesia by administering a small bolus of opioid. The inhalational gas being also an analgesic, another alternative is to increase the concentration of the inhaled agent to 1.5 M A C . If during the surgery there is a long period of low stimulation, the end tidal gas concentration may need to be reduced to lessen the level of cardiac depression. The widespread use of inhalational agents for the maintenance of anesthesia is mostly attributed to their ease of use. As mentioned earlier, the measurement of the end-expired vapour concentration provides the anesthesiologist with a reliable feedback quantity. In some cases, an infusion of propofol is given concurrently with an inhaled agent. Usually, the infusion is set at a fixed rate and the concentration of the inhaled anesthetic is tuned according to the patient's need. This method makes it possible for the anesthesiologist to reduce the concentration of the inhaled anesthetic to about 0.5 M A C , as the intravenous anesthetic will deepen the level of hypnosis. This usually results in less side-effects, especially post-operative nausea. In recent years, Total Intravenous Anesthesia (TIVA) has been strongly advocated by many anesthe-siologists [26], and has gained popularity with the introduction of faster agents with shorter duration of effect [27]. In T I V A schemes, no inhalational agent is used. Hypnosis and analgesia are achieved using only 2 agents: an anesthetic drug (usually propofol) and an opioid (alfentanil or remifentanil in recent years). The major difficulty for T IVA schemes is the lack of arterial blood measurement of the concentration of each drug. Without this feedback information, precise titration of the quantity of drug necessary to achieve a particular endpoint is not possible. Hence, the development of T I V A has been made easier through the development of pharmacokinetic model-driven infusion devices. These devices reach a desired plasma con-centration by using a computer-controlled infusion pump driven by the known pharmacokinetic parameters of the drug. The first Target Controlled Infusion (TCI) devices were made available in the late 1980s, but they have not been yet approved for use in North America. The main advantage of T I V A is that patients wake up faster and experience less postoperative nausea, thus shortening their stay in post operative care units. However, T I V A is still a marginal technique, as the anesthesiologist remains blinded to the true plasma concentration. Also, the incidence of intraoperative awareness and movement is higher. Finally, recent economic studies have concluded that T I V A usually results in higher costs (up to three times the cost of standard practice) [28], [29]. In today's cost sensitive environment, this is a major deterrent to this technique. Emergence The emergence from anesthesia is simply achieved by turning off the vaporizer as well as any infusion device used during the surgery. This is usually done during skin closure so that the patient wakes CHAPTER 1. CLINICAL ANESTHESIA: TERMINOLOGY, CONCEPTS, AND ISSUES 14 up faster at the end of surgery. A n additional bolus of a long acting opioid may be given for postoperative pain management. If the patient does not breathe spontaneously and still needs artificial ventilation, the anesthesiologist may use naloxone to antagonize the effect of opioids. However this practice is rare, as it may be associated with pain and nausea. 1.5 Drug Delivery and Monitoring in Anesthesia This introduction of clinical anesthesia would not be complete without introducing unfamiliar readers to the different devices that are currently used. Most of the information presented here is from [30] and [2]. 1.5.1 D r u g D e l i v e r y D e v i c e s W i t h the exception of nitrous oxide, all anesthetic agents can be found in liquid form. Their administration is therefore ensured by either a vaporizer, an infusion pump, or directly as a bolus (manual injection). T h e A n e s t h e s i a M a c h i n e The Anesthesia Machine delivers gaseous drugs and vapors (air, N 2 0) and has the ability to control ventilation. Its role is essential to ensure patient safety in terms of adequate respiration and blood oxygenation (even for procedures where no inhaled anesthetic is used, patients are intubated in order to support their breathing through mechanical ventilation). The anesthesia machine comprises four distinct devices: an oxygen and N 2 0 flow meter, a vaporizer, a breathing circuit, and a mechanical ventilator. The vaporizer mixes the anesthetic vapor with fresh gas. A knob permits the anesthesiologist to set the anesthetic vapour concentration at the outlet of the vaporizer. The breathing circuit provides the patient's lungs with the gas mixture via a face mask, the endotracheal tube or L M A . Pressure, volume, gas and vapour composition can be monitored and manually controlled by the anesthesiologist. Sophisticated breathing circuits can maintain temperature, water content, and usually incorporate a C 0 2 absorption cannister to remove the excess carbon dioxide in close circuit systems (systems in which the expired gas is recirculated for maximum efficiency in terms of anesthetic gas usage). The breathing cycle can be either controlled directly by the patient (spontaneous ventilation) or by a ventilator (mechanical ventilation) which maintains the respiration set by the anesthesiologist. Patients undergoing general anesthesia and unconscious sedation in the I C U are often connected to a mechanical ventilator. However, mechanical ventilation has truly an active role in two cases: - w h e n us ing opioids: large doses of opioids can significantly depress the respiration cycle. The brainstem mechanism regulating the blood C 0 2 concentration becomes less sensitive, thus inducing a slowing down of the respiration rate and a reduction of the inspired volume. In some cases, the administration of opioid boluses can be followed by a prolonged apneic state during which proper oxygen concentration in the brain could be seriously compromised. The mechanical ventilator then assumes the role of the autonomic function by ventilating the patient. CHAPTER 1. CLINICAL ANESTHESIA: TERMINOLOGY, CONCEPTS, AND ISSUES 15 - w h e n us ing N M B s : N M B drugs block the transmission of synaptic motor information. Hence, even though the autonomic respiratory function sends the appropriate activation signals to the lungs, their action is inhibited and no respiratory movement follows. As a result, mechanical ventilation is necessary. Using the ventilator, the anesthesiologist can decrease the blood CO2 concentration below the threshold where the autonomous respiratory function does not drive respiration. This can be helpful when the breathing pattern needs to be very regular and thereby avoid the use of N M B s or large opioid doses. In terms of anesthesia delivery, the anesthesia machine is mainly used to maintain the anesthetic ade-quacy during the maintenance phase. It is generally not used to induce the patient. Infusion P u m p s Most intravenous drugs are administered in the form of boluses to correct for inadequa-cies in the anesthesia regimen (e.g., during period of intense surgical stress). However, to limit hemodynamic depression, post-operative nausea and airway irritation, infusion pumps are used to supplement the anes-thetic gas with a continuous infusion of intravenous anesthetic such as Propofol. While opioids can also be delivered via an infusion pump, this practice is rare in the OR, and is limited to specific cases (e.g., neurosurgery). The infusion pump usually comprises a syringe plunger actuated by a high precision stepping motor. In newer versions, the anesthesiologist can input the type of drug infused, the syringe caliber, the patient's weight, and the desired infusion rate, usually in ug • m i n - 1 • k g 1 or in ng • m i n - 1 • k g 1 . Based on this information, the pump automatically controls the rate of descent of the plunger. A n alarm warns the practitioner when the plunger reaches the end position. These pumps are also usually fitted with a RS232 serial port. Using the proper communication protocol, it is possible to remotely control the pump from an external device such as a laptop. In some countries, T C I devices can be used. These devices target the blood plasma concentration by calculating the infusion profile according to the pharmacokinetic model of the infused drug. 1.5.2 A n e s t h e s i a M o n i t o r s Standard of care during anesthesia procedures includes the use of a number of sensors: - Pulse Oximeter: a light transducer clamped conveniently on one of the patient's fingers permits the calculation of heart rate and oxygen saturation. - Capnograph: a small tube connects the endotracheal tube or L M A to the capnograph. A sample of the gas expired by the patient is continuously analyzed. The capnograph permits the measurement of the end-tidal CO2. - Mass Spectrometer: this device is integrated to the breathing circuit. Samples of the inspired and expired gas mixture are analyzed by mass spectrometry to determine the inspired and expired con-CHAPTER 1. CLINICAL ANESTHESIA: TERMINOLOGY, CONCEPTS, AND ISSUES 16 centrations of the anesthetic gas. This is particularly convenient to calculate the quantity of the anesthetic agent that enters the patient's arterial blood stream. - Hemodynamics Monitors: E C G leads located on the patient's chest record the Electrocardiogram ( E C G ) . The E C G signal is mainly used to determine the heart rate. However, the very shape of the E C G (referred to as QRS complex) can be indicative of cardiovascular difficulties. The E C G is displayed continuously. A pressure cuff also allows the measurement of the brachial blood pressure. A rise or fall in systolic or diastolic blood pressure also indicates cardiovascular difficulties. This measurement is intermittent in order not to block systemic circulation in the instrumented limb. The information gathered by these sensors is displayed on a unique screen, usually referred to as the Anesthesia Monitor. The anesthesiologist's expertise in interpreting these readings permit the assessment of the patient's state and needs. 1.6 Summary Paradoxically, surgeons achieve healing by first inflicting injury. Anesthesiologists use anesthetics and opioids to prevent the awareness of pain and attenuate the body's stress response to injury. The uncon-sciousness produced by general anesthesia is accompanied by a depression of the respiratory, cardiovascular and endocrine responses to surgery. Since the degree of surgical stimulation changes during surgery, anes-thesiologists must constantly adjust the extent of anesthetic depression to avoid both under- and overdosing. Otherwise, excessive activation of the sympathetic nervous system or pharmacological depression could in turn lead to injury of critical organs, especially in patients with limited respiratory and cardiovascular reserves. Unti l the mid-20 t h century, only inhaled anesthetics were available to create the state of general anes-thesia. However, their onset of action was slow and often accompanied by vomiting and signs of respiratory irritation. In the 1940s, the introduction of intravenous agents revolutionized the medical specialty of anesthesia, and over the past 50 years, further refinements in anesthetic pharmacology, equipment and technology have occurred. Nowadays, anesthesiologists have access to agents that can act within a minute of their administration and can block specific mechanisms such as cognition, awareness, memory, stress response and muscle movement. These agents are further characterized by their fast metabolism and elim-ination. Hence, constant monitoring of their titration is necessary to provide patients with an adequate drug regimen during surgery. Chapter 2 Quantifying Depth of Anesthesia: a Review In today's anesthesia practice, patient monitoring is ensured by a set of transducers to measure cardio-vascular and respiratory parameters. These measurements enable experienced anesthesiologists to assess their patients' state and take appropriate actions. For instance, a rise in blood pressure and heart rate is usually associated with surgical stress. Conversely, a drop in blood pressure and respiratory rate is often associated with unnecessary overdosing. However, these physiological measurements must be interpreted only within the surgical context. Movements and other visual cues (lachrimation, grimacing, etc.) are often the ultimate warning signs of an inadequate anesthetic regimen. Immediate corrective actions must then be taken to avoid intra-operative awareness. W i t h the introduction of N M B drugs in the early 1940s, this safety net disappeared. A growing number of intra-operative cases associated with the use of N M B s were reported. The search for a Depth of Anesthesia (DOA) monitor to quantitatively assess the anesthetic state of the patient became an important focus in research. In terms of close loop anesthesia, the availability of a feedback measure is a sine que non condition to the success of such an endeavor. While we have developed our own solutions (see Chapter 5), this chapter presents a review and gives a historical perspective on this very active research field. Some notions about the use of surrogate measures are presented in Section 1. We will then take a more in-depth look into the techniques that have been developed to measure hypnosis. One of these technique, based on bispectral analysis, was first commercialized in 1996 and has received considerable attention from the clinical community. The experience gained by anesthesiologists in using this type of monitor sheds some light on the eventual advantages of automation in anesthesia. As such, the clinical relevance, interests and potentials of this type of monitoring is reviewed in Section 3. Finally, we will conclude this Chapter with a short section dedicated to the research in the field of analgesia measurement. This particular topic has not received much attention as of yet, which explains the rather limited amount of available information. 17 CHAPTER 2. QUANTIFYING DEPTH OF ANESTHESIA: A REVIEW 18 2.1 The Use of Surrogate Measures Hypnosis and analgesia are the result of different pharmacological mechanisms within the CNS. It is not possible to measure them directly. However, there are physiological signs that can be sufficiently correlated with these two anesthesia endpoints. For instance, it is known that cortical activity mirrors the state of hypnosis for a patient. The fact that the brain can or cannot process sensory information can be observed in the E E G . Similarly, it is a well known fact that opioids depress the respiratory autonomic response of the body in a dose dependent fashion. Monitoring the changes in respiratory rate and blood CO2 can thus provide an indirect measurement of the opioid effect. The use of surrogate measurements to assess hypnosis and/or analgesia is subject to a limited therapeutic window. When using the E E G to determine the effect of anesthetics and opioids, great care has to be taken into defining a proper therapeutic window. For instance, in their study of opioids, Bil lard and Shafer [31] observed that the E E G could represent adequately the opioid effect, but only for a limited range of doses, see Fig . 2.1. Small doses do not have a marked influence on the E E G , while larger doses used for cardiac surgery reach a maximum E E G effect. Therefore, the use of the E E G seems to be irrelevant in preoperative and postoperative I C U (where often only a shallow analgesic state is provided) and in surgeries where cardiac stability is ensured by large opioid amounts. E F F E C T SITE ALFENTANIL [ng/ml] Figure 2.1: Therapeutic window of the E E G for measuring the opioid effect (from [31]). Various authors ([32], [33], [34]) have been trying to define the characteristics of an ideal index of Depth Of Anesthesia (DOA) . These authors were referring to either hypnosis or analgesia. They all agree on the following characteristics that an ideal index of hypnosis must possess: i . vary in a predictable and consistent fashion as the drug concentration increases, i i . show similar changes for different agents of the same family (anesthetics, opioids), i i i . indicate shallow levels of hypnosis or analgesia, iv. have a stable baseline value that varies only minimally, CHAPTER 2. QUANTIFYING DEPTH OF ANESTHESIA: A REVIEW 19 v. characterize the maximal change induced by the drug, and, vi . recover to the baseline value on discontinuation of the drug, or upon termination of the effect. 2 .2 Q u a n t i f y i n g t h e D e p t h o f H y p n o s i s : a R e v i e w The Central Nervous System, and in particular the brain (cerebral cortex) are the target organs of anesthetic drugs. Since sensory perception and consciousness are processed at this level, it is reasonable to assume that electroencephalogram signals (EEG) - which reflect cortical activity - can be used to determine the depth of drug-induced unconsciousness. The effects of anesthetic drugs on the E E G have been known since the early 40s [35] where neurophys-iologists observed that the E E G of anesthetized patients contain slower waves with higher amplitudes, see Fig. 2.2. Later on, Faulconer and Bickford [36] published a review cataloguing the different patterns of the E E G associated with ranging depths of anesthesia for different anesthetic agents. Bickford also showed that the E E G measured from the scalp is similar and synchronous to the E E G that can be measured at depth in the human brain. A n extensive review published in 1973 by Clark and Rosner [37] concluded that all anesthetics have an effect on the E E G . But even though all studies have shown that there is a clear correlation between the E E G and the anesthetic depth, they have also shown that patterns generated through anesthesia were different according to the drug that was used. Some researchers concluded that neuroelectric recordings couldn't provide a simple and uniform measure of anesthetic depth. Stanski [33] concluded that in order to use parameters derived from the E E G as a measure of the anesthetic depth, one must first study, for each drug, the correlation between the plasma concentration, the E E G effect and the resulting clinical anesthetic state of the patient. He also stressed out that the pharmacokinetic and pharmacodynamic model of the drug must be associated with the parameter to accurately represent and predict the anesthetic depth. 2.2.1 R e v i e w o f T o o l s a n d T e c h n i q u e s There are a number of signal processing tools and techniques available to quantify the E E G in order to derive a surrogate measurement of hypnosis. Here is a brief summary of these tools: 2.2.1.1 Time Domain Methods Bickford, Faulconer and Soltero ([38], [39], [40]) quantified the energy of the E E G by rectifying and inte-grating the signal. When the resulting integrated signal reached a pre-specified value, the system generated a pulse and reset the integrator. The frequency at which pulses were generated was an indication of the hypnotic depth. Later on, Bellville [41] developed a method using cyclopropane anesthesia. He observed that increasing concentration of cyclopropane was also increasing the E E G amplitude. He thus used the raw amplified CHAPTER 2. QUANTIFYING DEPTH OF ANESTHESIA: A REVIEW 20 Left motor region Right motor region (a) (c) (e) Breathing cyclopropane Breathing air (d) (f) (g) (h) Figure 2.2: Simultaneous records from the left and right motor regions, illustrating the changes in frequency during cyclopropoane anesthesia and recovery, (a) awake (b) breathing cyclopropane for 2 min. (c) 1 min later (d) immedi-ately on substitution of room air for cyclopropane (e) 3 seconds later (f) 2 min. later (g) 5 min. later (h) complete recovery (from [35]). CHAPTER 2. QUANTIFYING DEPTH OF ANESTHESIA: A REVIEW 21 FREQUENCY [HZ] FREQUENCY [HZ] Figure 2.3: Awake and anesthetized E E G and their power spectrum, signal directly as a feedback quantity. 2.2.1.2 Power S p e c t r u m A n a l y s i s W i t h the development of microprocessors and signal processing tools, researchers have focused their atten-tion on Fourier analysis of the E E G . Power spectrum analysis is used to obtain a frequency distribution of the E E G . Hence, any change in the frequency content of the signal can be visualized. Pichlmayr et al. [42] published a thorough review of the effect of the different anesthetic agents on the E E G spectral distribution. The E E G spectral distribution is characterized by different modes. It is common practice to distinguish between 5 frequency bands: 5 band (0.25 Hz - 3.5 Hz), 0 band (3.5 Hz - 7.5 Hz), a band (7.5 Hz - 12.5 Hz), /? band (12.5 Hz - 32 Hz), and 7 band (32 Hz - 70 Hz). For a normal awake patient, the E E G activity is principally concentrated in the delta and alpha bands. Wi th increasing level of anesthetics, the activity of the alpha band tends to reduce, while the low frequency content of the delta band is increased, see Fig. 2.3. To quantify the effect of anesthetics on the E E G , researchers have tried to derive univariate indexes based on Fourier analysis that characterizes the spectral distribution. Among the parameters that have been thoroughly investigated, we can mainly distinguish between the following two: M e d i a n Frequency ( M E F ) The M E F is the frequency that splits the power spectrum distribution into two parts of equal power. This descriptor is advocated by Schwilden and co-workers ([43], [44], [45] and [46]) who used it for close loop anesthesia. S p e c t r a l Edge Frequency ( S E F ) The SEF is the frequency at which 95% of the E E G power is present. This descriptor was proposed in 1980 by Rampil et al. [47]. According to its authors, the S E F is highly repeatable but presents large inter-subject variations. One advantage of the S E F over the M E F is its high sensitivity to deepening hypnotic levels. As a result, this descriptor is still being used today, in research. CHAPTER 2. QUANTIFYING DEPTH OF ANESTHESIA: A REVIEW 22 Levy [48] argued that, since E E G activity exhibits multimodal patterns, no univariate descriptor based on the processed E E G can serve as a consistent and adequate representation of the power spectrum. Levy thus doubted the reliability of univariate descriptors of D O A based on the E E G power spectrum. 2.2.1.3 E E G modeling The idea of modeling the E E G using auto-regressive techniques dates back to the early 1970s. But its application to anesthesia is more recent ([49]). As compared to univariate descriptors, auto-regressive (AR) modeling generates a set of parameters that can further be correlated to anesthetic depth. In order to derive a single index from the A R parameters, a neural network is trained. Sharma et al. [50] showed that this technique can lead to accurate results for measuring analgesia in dogs using halothane. However, the size of the proposed network is rather large (43 inputs, 43 nodes in the primary hidden layer, 6 nodes in the secondary hidden layer and one output). One of these inputs is actually the M A C value. Results are seriously compromised when removing this information, thus indicating that the anesthetic depth measured by this technique depends heavily on the dose-response characteristic of the inhaled anesthetic. 2.2.1.4 Bispectral Analysis: The B I S ™ Monitor Rampil, a leading researcher in the field of neuro-anesthesia, recently argued that anesthetic agents tend to synchronize the generation of postsynaptic potentials [51], thus resulting in slower waves of higher amplitude in the E E G . W i t h increasing levels of anesthesia, it is expected that some of the frequency components of the E E G will shift in time (change of phase). This phenomenon should be mostly predominant in the lighter anesthetic states. This change in latency is not observable by spectral analysis, as phase information is usually discarded. To illustrate this point, Bowles et al. [52] proposed the following example: let us consider a signal composed of 3 sinusoid waveforms. The phases of two of the frequency components are fixed, while the phase of the third one is allowed to slowly drift (see Fig. 2.4). As a result, the shape of the signal itself changes according to the phase of the third component. Clearly, even though the shape of the signal is different, standard power spectral decomposition is not able to capture the phase drift. However, when calculating the bicoherence value, a large discrepancy can be noted when the phase of the third component is coupled with the phases of the first two components. Hence, bispectral analysis can differentiate between these two signals and capture phase drifts. Ning et al. [53] applied bispectral analysis to the E E G in order to characterize sleep patterns in rats. They realized that there was a strong coupling between the frequencies of 6 and 8 Hz during Rapid-Eye Movement ( R E M ) sleep. Sleep patterns being close to patterns obtained during anesthesia procedure, Ning assumed that this technique might lead to interesting results in monitoring the depth of anesthesia. Their CHAPTER 2. QUANTIFYING DEPTH OF ANESTHESIA: A REVIEW 23 • Random phase signal: Figure 2.4: Changes in latencies in a signal can be tracked by the bicoherence index (from [52]). assumption was validated in 1990, when Kearse et al. [54] reported that an index based on bispectral analysis was more accurate than the spectral edge frequency for opioid-induced anesthesia (alfentanil and sufentanil). These findings were confirmed by Sebel et al. [55] with isoflurane and later on with propofol [56]. A study by Muthuswamy and Roy in 1993 on anesthetized dogs (halothane) also led to similar conclusions [57]. However, they also noticed that the accuracy of the analysis could probably be improved by using bicoherence indexes in association with other parameters such as the M E F , S E F and hemodynamic parameters. Probably the most compelling result was obtained by a research team from Aspect Medical Systems Inc. who derived no less than 33 variables (11 bispectra, 11 bicoherence indexes, 11 power spectral values). These variables were combined into a dimensionless index in order to predict somatic responses to surgical incision [52] 1. This index, referred to as BIS (Bispectral Index Scale), is a weighted sum of each spectral and bispectral variable. In the earliest versions, the weights were tuned using discriminant analysis and based on data collected from 170 standard surgical procedures. In comparison with other methods (SEF, M E F , etc.), the accuracy of the BIS was significantly higher. A major advantage of bispectral analysis over more conventional techniques is its wider therapeutic window that allows the differentiation between different levels of sedation [58]. Based on these findings, Aspect Medical Systems has developed a monitor of consciousness, referred to as the BIS Monitor, which integrates this technology. A single E E G signal is recorded through contact 1Note once again the confusion between hypnosis and analgesia. CHAPTER 2. QUANTIFYING DEPTH OF ANESTHESIA: A REVIEW 24 BIS Awake • Memory Intact Light Hypnotic ^ State Moderate Hypnotic State Deep Hypnotic < State [20}-Light/Moderate Sedation r 1 Deep Sedation General Anesthesia • Deep Hypnosis Increasing Burst Suppression Low Probability of Explicit Recall Low Probability of Consciousness • Memory Function Lest Isoelectric EEG Figure 2.5: Bispectral Index Scale and its meaning (from [59]). electrodes placed on the patient's forehead. The BIS value and other relevant information are then displayed on the monitor screen. To simplify the interpretation, the BIS is scaled between 0 and 100. A value of 100 represents the awake state. Wi th increasing concentration of anesthetics, the index decreases. General anesthesia is obtained for an index between 60 and 40, see Fig. 2.5. Lower values represent deep hypnotic states, while values between 90 and 60 generally represent sedation levels. 2.2.1.5 Quantitative Evoked Potentials: the A - L i n e ™ Monitor Somatosensory information provoked by auditory, visual, or tactile stimulation generate transitory oscilla-tory signals within the E E G itself. It has been advocated that such transient signals, if properly analyzed, can reveal information concerning a patients' state of consciousness ([32, 60, 61, 62, 63]). For instance, midlatency auditory evoked potentials ( M L A E P ) have a very distinct shape whether the subject is awake or asleep, see Fig . 2.6. The most remarkable feature beside the change in amplitude of the signal is the change in latency of the Pa and P;, waves. A very interesting work by Huang et al. [64] has shown that it is possible to measure the hypnotic depth of a dog under anesthesia by using wavelet transform of the M L A E P signal and feeding the wavelet coefficients to a properly trained neural network. A significant work has been carried out in that particular field. One major disadvantage of using evoked potentials is their very low signal to noise ratio, which makes them particularly difficult to acquire, as they are embedded inside the E E G signal.. Averaging over a large number of samples is necessary in order to filter out the E E G noise from the signal. New techniques such as wavelet de-noising have been successfully applied to partially alleviate this problem [65]. The latest developments in A E P extraction using A R X modeling have yielded results that allow the CHAPTER 2. QUANTIFYING DEPTH OF ANESTHESIA: A REVIEW 25 4 2 + Awake: Cp—3 [ig/ml 2 4-Sleep: 0^ ,-5.5 [ig/ml -14- ms ms tr Auditory stimulus Auditory stimulus Figure 2.6: Midlatency Auditory Evoked Potentials obtained from a responsive (awake) and non-responsive (asleep) dog subject to tail clamping. The waveforms have been obtained after averaging of 1000 samples (adapted from [64]). estimation of the anesthetic state based on 6 second epochs [66, 67], thus necessitating the averaging of only 15 sweeps. This new technology has been embedded in the A-Line Monitor (Alaris Medical Systems, California) and is now commercially available. 2.2.1.6 Other Monitors Following the success of the BIS Monitor, medical companies have taken a keen interest in this new market segment. New monitors have recently become available: - The P S A 4000™ (Physiometrix, Massachusetts): this monitor displays an index of hypnosis calcu-lated on 1.25 s E E G epochs. Numerous quantitative measurements are derived (power gradients, absolute power in specific frequency bands, etc.). These measurements are further combined into a unique dimensionless index, referred to as the Patient State Index (PSI). The PSI is scaled between 100 (awake) and 0 (isoelectric E E G ) . Three extensive E E G databases were used for the tuning and calibration of the PSI. While the PSI is based on similar principles as the BIS Monitor (i.e. compos-ite index tuned from E E G databases and clinical cases), it differs in that it does not use bispectral parameters. Instead, the PSI focuses on the power shift in specific frequency components between the frontal cortex and the posterior lobes. This power shift is directly related to the hypnotic depth. Paradoxically, the commercial version of the P S A 4000 uses only frontal electrodes. - The Narcotrend® Monitor (MonitorTechnik, Bad Bramstedt, Germany): this monitor is based on a pattern recognition algorithm that classifies the E E G into 14 different classes using Kugler's clas-sification and denomination [68] (A, B0, B l , B2, CO, C I , C2, DO, D l , D2, E0, E l , F0, F l ) . The range D0-E1 indicates an adequate depth of anesthesia for surgery. The range F0-F1 indicates burst suppression and isoelectric E E G . The Narcotrend index exhibits step-wise changes during increasing or decreasing hypnotic depth. - The A S / 5 M-Entropy Module (Datex Ohmeda, Helsinski, Finland): Datex Ohmeda is a leading CHAPTER 2. QUANTIFYING DEPTH OF ANESTHESIA: A REVIEW 26 manufacturer of anesthesia machines and monitors. They recently developed an E E G measure based on entropy analysis that captures the complexity of the E E G signal [69]. During the awake state, the E E G is a highly complex signal, resembling that of noise (no apparent pattern). W i t h increasing depth of hypnosis, patterns emerge, thus reducing the complexity of the signal. 2.2.2 Clinical Relevance, Interest and Potentials of the BIS Monitor The BIS Monitor has received considerable attention from the research and clinical communities. Over 1,250 papers discussing its pros and cons and its clinical relevance have been published since its commercial in-troduction in 1996. The BIS Monitor was the only F D A approved monitor of drug-induced unconsciousness until 2002. Intra-operative Awareness The most publicized role of the BIS Monitor is the reduction of the in-cidence of intra-operative awareness. However, validating such an assertion is a difficult task due to the rather limited occurrence of such cases. It is only very recently that two large studies were able to assess the impact of BIS monitoring to prevent intra-operative awareness: - SAFE2 Trial: following the S A F E 1 trial (see Section 1.1.2.2), Sandin et al. conducted a prospective study where 5,057 patients undergoing surgery requiring general anesthesia including the use of muscle relaxants or intubation were monitored with the BIS monitor. Only 2 patients experienced recall of intra-operative events, as compared to 14 in the control group (from the S A F E 1 Trial). The authors concluded that the BIS monitor can reduce intra-operative cases by up to 80%. Even though the study population is still rather limited (only 5,000 cases so far), this result was deemed statistically relevant. - B-AWARE Trial: the B - A W A R E Trial targets a population of patients at risk undergoing relaxant general anaesthesia. A total of 2,503 patients were enrolled in this study. This population was divided into 2 groups: control, and BIS monitored. In the BIS group, the index was kept in the 40-60 range. The incidence of intra-operative recall was about 0.89% (11 cases) in the control group, and 0.16% (2 cases) in the BIS group. Clearly, the BIS monitor succeeded in reducing intra-operative awareness in a significant manner. Drug Usage, Post-operative Acute Care Unit Discharge Times and Outcome Another well-publicized advantage of a BIS guided titration (i.e., maintaining the BIS value in-between a pre-determined range during the maintenance phase) is the reduction of drug usage and discharge time. For instance, several studies [70, 71, 72] concluded that an average reduction of 20% to 30% in the total amount of anesthetic drugs (inhaled and/or intravenous) could be achieved during the maintenance phase by keeping the BIS value in the 40-60 range, while still ensuring the adequacy of the anesthetic state. These studies CHAPTER 2. QUANTIFYING DEPTH OF ANESTHESIA: A REVIEW 27 25 Figure 2.7: Time of arousal following discontinuation of propofol and remifentanil infusion (from [70]). imply that, on average, patients without BIS are overdosed. Processed electroencephalography variables such as the BIS can thus offer practitioners a guideline to titrate each patient according to his/her specific needs. W i t h the reduction of drug usage, it has been shown that patients awake sooner from their surgery and are discharged faster than in standard practice [73, 70, 74]. Discharge times in the Post Anesthetic Care Unit (PACU) are reduced by up to 50%, see Fig. 2.7. The incidence of post operative nausea and vomiting, and patients' comfort are also significantly improved. The reduction in drug usage, discharge times and adverse post-anesthetic events (nausea and vomiting) have the added advantage of reducing anesthesia-associated costs. In today's cost-sensitive society, this could be seen as a major advantage. However, some studies [75] have shown that the savings compensate for the cost of the BIS sensors and the maintenance of the monitor only for long surgeries (>4 hours). 2.2.2.1 Long Term Benefits of BIS Monitoring In a recent abstract, Weldon et al. [76] showed an interesting correlation between one year post-operative death rates and the average BIS value. The study included 907 patients scheduled for major surgical procedures. Logisitc regression modeling showed that increasing age and lower BIS levels were both inde-pendently associated with higher one year post-operative mortality rates. In the elderly population, the mortality rate was 16.7% for an average BIS value <40, as compared to 4.2% for an average BIS >60. This study implies that maintaining middle-aged.and older patients into very deep anesthetic levels might not be advisable. While this result was later confirmed by Lennmarken [77], details from these studies have still not been released. One could theorize that patients driven into deep hypnotic states received more inhaled anesthetic and perhaps less opioids, and as a result, suffered from cardiovascular instability. TIME [mini CHAPTER 2. QUANTIFYING DEPTH OF ANESTHESIA: A REVIEW 28 2.2.2.2 Bad Press There are still many controversies concerning the use and reliability of the BIS monitor. One major argument fuelling the debate is whether a single index derived from the E E G could reflect the anesthetic adequacy in all types of anesthetic procedures. Many argue that in order to validate such monitors one would need to study the effect of each anesthetic and opioid drugs (and combination thereof) on the index. A second argument is that the BIS is tuned to reflect the anesthesiologist's assessment. As such it is an empirical measurement which might not reflect the underlying neurophysiologic mechanism of anesthetics. Finally, while many studies have concluded that a BIS guided titration results in lower drug usage and improved outcome, some consider that such titration would increase the risk of intra-operative awareness [78]. However, this argument seems to have been refuted by the recent B - A W A R E trial. Particular cases where the BIS has shown paradoxical behavior were also reported: - Low BIS values during recovery: Sleigh and Donovan [79] report that 14 out of 37 patients could respond and follow verbal commands during emergence while the BIS value was <75. They observed that the BIS value regains its baseline value during emergence only 30 s after patients start reacting to verbal commands. They attributed this delay to the update delay of the BIS and its 30 s averaging window. A similar observation was later on reported by Lehmann [80]. - Effect of E E G amplitude on the BIS: Schnider et al. [81] reported that a consistent BIS value of 40-50 can be obtained if the amplitude of the E E G is low enough. Muncaster et al. [82] in 2003 also reported having a patient with a very low BIS value during emergence due to a low E E G amplitude. - N M B s and BIS: The effect of N M B s on the BIS value has also recently been investigated. In one study on volunteers, Messner et al. showed that a paralyzed fully conscious subject can have a BIS value as low as 9 [83]. If confirmed, this could be a very damaging result for a monitor that is supposed to detect episodes of consciousness in paralyzed patients. A similar finding was reported by Vivien et al. [84]. In their study of 45 patients in intensive car unit, they measured an average drop of 25% in the BIS value after the administration of N M B s . A similar decrease in electromyographic (EMG) activity was measured. A conclusion from these studies is that the BIS value is strongly correlated to the E M G and should not be used in sedated patients whenever N M B s are used concurrently. Note however that N M B s do not affect the BIS value in patients under general anesthesia [85]. This could be due to an already very low E M G activity, where the N M B impact remains very l imited 2 . One can therefore hypothesize that the BIS value cannot reach levels above 60-70 without some E M G activity. As a result, paralyzed patients should not be allowed to be maintained at levels close to 60 since an episode of intra-operative awareness might not be picked up by the monitor 3. 2The contribution of the EMG signal in the EEG remains unaffected by the NMB (however the drug still effectively blocks muscle movements). 3Personal interpretation. CHAPTER 2. QUANTIFYING DEPTH OF ANESTHESIA: A REVIEW 29 Finally, a cornerstone of the mistrust and recurring complaint by clinicians is that the BIS algorithm has never been made public, and hence cannot be rigourously evaluated by independent researchers [86]. 2.2.2.3 Comparison with Other Monitors Because of the relative novelty of the other monitors, only a few clinical studies have been published comparing the bispectral index with its competitors. So far, results from the P S A 4000 [87] and the Narcotrend Monitor [70] seem to show a good correlation between these indexes and the BIS. No clear advantages have been found in using one particular technology, and both the P S A 4000 and the NarcoTrend monitor reported similar findings in terms of drug usage and P A C U discharge times. 2.3 Quantifying the Depth of Analgesia The measurement of antinociception is part of the critical care given by anesthesiologists in their daily practice. Surgical trauma is usually accompanied with strong sympathetic and parasympathetic activity (e.g., heart rate and blood pressure changes, sweating, lachrimation, somatic movements, etc.). Deriving descriptors of antinociception that can be used as a feedback quantity is a challenge. For instance, blood pressure alone is not a reliable measure, as other parameters such as blood loss and the action of vasoactive drugs can affect the cardiovascular system. As a result, the analgesia functional component hasn't yet received enough attention. In recent years, only a few monitors have been introduced, without full disclosure of their scientific base. 2.3.1 Using Physiological Measures 2.3.1.1 The E E G as a Measure of Analgesia The effect of opioids on the E E G has been thoroughly investigated. In 1984, Smith et al. [88] concluded that the E E G probably reflects the depth of anesthesia with high-dose narcotics. Later studies by Scott et al. ([89] and [90]) resulted in similar conclusions. The same techniques ( M E F , S E F , BIS) described previously can be used for measuring the effect of opioids as well. However, the therapeutic window of these techniques is very limited. Also, many agree today that cortical activity, while appropriate for consciousness monitoring, cannot be used to assessed the effect of opioids since their site of action is principally localized at the level of the spinal cord [91]. 2.3.1.2 End-tidal C 0 2 A l l opioids depress ventilation in a dose-dependent fashion. As a result, the amount of carbon dioxide (CO2) in the arterial blood during spontaneous ventilation is an indication of the opioid effect. The advantage of using the CO2 arterial partial pressure is that it can be easily measured using the capnograph. However, the usefulness of this measurement is limited to cases where patients breathe spontaneously (a majority of CHAPTER 2. QUANTIFYING DEPTH OF ANESTHESIA: A R E V I E W 30 day care surgeries), thus limiting the therapeutic window of this parameter to doses of opioids that do not provoke apnea. 2.3.1.3 Heart Rate and Blood Pressure A rise in heart rate and blood pressure is usually indicative of increased sympathetic activity (i.e. stress). Monitoring these parameters can be useful to determine the level of stress. However, heart rate and blood pressure baselines are patient-dependent. Also, numerous drugs other than opioids can depress or alter these signs. Finally, changes in heart rate and blood pressure are only indicative of strong stimulation and do not provide a graded response to increasing stimulation. 2.3.2 Using Heart Rate Variability In healthy subjects, parasympathetic and sympathetic activity are balanced according to stress levels, circadian rhythms, position, etc. This balance is regulated by the CNS. In most subjects under stress, sympathetic activity dominates. Conversely, relaxation increases parasympathetic activity and depresses sympathetic tone. While the field of Heart Rate Variability (HRV) is quite prolific, most of the research done t i l l now was oriented towards the development of diagnosing tools for awake subjects. In standard analysis of the H R V , researchers differentiate between the H F band (0.2 - 0.5 Hz) which reflects parasympathetic activity, and the L F band (0.01 - 0.15 Hz) which is believed to reflect both sympathetic and parasympathetic activity (the exact mechanism remains unknown). High frequency activity is usually concentrated around the respiratory frequency. One explanation could be that the regulation of respiration is mediated by the same center that regulates the activation of the vagal nerve. Typically, the H F measurement usually quantifies the activity concentrated around ±0.06 Hz around the respiratory frequency. In order to normalize this measurement, standard practice usually involves the ratio L F / H F , which is believed to yield a more consistent measurement of the sympathetic/parasympathetic balance. The effect of anesthesia on the H R V signal has not yet been thoroughly investigated. A standard measurement practice for instance has not been adapted, since the respiratory rate in anesthetized patients can drop to a point where both the H F and L F activity share the same frequency band. Hence, it seems that normalizing one with respect to the other is not indicated. 2.3.2.1 Existing Technologies Two monitors have been advocated as analgesia monitors: the Anemon I, Medical Control System SA, and the Fathom, Amtec. There is, however, no published literature reviewing the merits of these technologies. The Fathom monitor is still at the prototype stage and is not yet commercially available. This particular device analyzes the respiratory pattern fluctuations. These fluctuations are difficult to measure. However, they strongly correlate with fluctuations of the cardiac rhythm called sinus arrythmia. Observing the CHAPTER 2. QUANTIFYING DEPTH OF ANESTHESIA: A REVIEW 31 changes of the interval between R-waves of the electrocardiogram and the respiration cycle forms the basis for this measurement ([92] and [93]). Clinical trials are currently being conducted in England to validate the efficacy of this technique. 2.4 Summary One hundred and fifty years after the discovery of anesthesia, the current standard of care in anesthesia management still does not include monitors to assess quantitatively the effect of anesthetic drugs on the target organ, i.e., the brain. The current practice still relies uniquely on secondary signs to warn the practitioner of either pharmacological toxicity or anesthetic inadequacy. It is only in the past 10 years that serious advances have been made. Since 1996, practitioners have access to monitors to measure the degree of drug-induced unconsciousness. As of today, monitors, such as the BIS, are still considered by many to be nothing more than gadgets. Also, the price of the BIS sensor deters many to use the device in their everyday practice. Further, the lack of predictability in the BIS behavior during transients (e.g., unconscious patients after induction can have a stable BIS >80 while a reacting patient after emergence can have a BIS <60) has brought many practitioners to question its reliability. However, with the new awareness studies, the use and interest of such monitors might be re-evaluated once the economic and social impact of intra-operative awareness are properly studied. In terms of analgesia, no real breakthrough has been achieved yet. This can be explained by the fact that anesthesiologists can assess this state better than hypnosis, since they monitor cardiorespiratory functions. Chapter 3 Modeling Anesthetic Drugs: the Traditional Approach In this chapter, we present the drug modeling approach favored in traditional pharmacological studies. In this approach, it is usual to first consider how the administered drug distributes within the body. This analysis leads to a pharmacokinetic model (PK) which can be used to predict the blood plasma concentration of the drug. The second step is to relate this concentration to the drug effect itself. This yields a second mathematical expression referred to as the pharmacodynamic model (PD). The organization of this chapter mirrors this approach. The pharmacokinetic concepts are covered in Section 3.1, while pharmacodynamic modeling is discussed in Section 3.2. Propofol and remifentanil being the fastest anesthetic and opioid agents available in the current practice, we believe that a close loop anesthesia delivery system will benefit the most from these agents. Hence, most of the discussion presented here is focused on these two drugs. 3.1 Pharmacokinetics Anesthesiologists act through the administration of drugs. A pharmacokinetic model of a drug is a mathe-matical expression relating the drug blood plasma concentration Cp(t) to the administered dose I(t). The aim of this section is thus to define the transfer function 1 PK(s): PK(s) = ^ - (3.1) While special emphasis will be given to propofol and remifentanil, most of the concepts presented here can be broadened to a wider variety of drugs. For a more in-depth look into these concepts, we invite 'In contrast with the usual approach in pharmacology, PK(s) is expressed here in the Laplace domain (or frequency domain). In linear system theory, the Laplace transform is used extensively to express time-domain differential equations into algebraic expressions, which are thus easier to solve. Continuous time domain system analysis and control design are often conducted in the Laplace domain. 32 CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 33 16 8 4 2 T I M E F R O M INJECTION [min] Figure 3.1: Simulation of the time course of the percentage of a. thiopental bolus in the central blood pool, viscera, lean tissue and fat as a function of time, assuming no elimination clearance (adapted from [97]). interested readers to refer to the textbooks from Stanski and Watkins [94], Prys-Roberts and Hug [95], and Fragen [96]. 3.1.1 Principles and Concepts Pharmacokinetics can best be understood when considering the time course of the concentration of any given drug within the plasma and other tissues of the human body (see Fig . 3.1). During the absorption phase following an intravenous bolus administration, a drug mixes rapidly within the central blood pool, resulting in a plasma peak concentration. While at first glance the peak plasma concentration seems to occur instantaneously, a small delay elapses between the actual injection of the drug and its mixing within the blood pool. Systemic circulation then distributes the drug to a variety of tissues within the body. A part of the drug enters these tissues through molecular diffusion, according to the affinity of the tissues for the drug, the rate of perfusion and the relative concentration of the drug in the blood and in the tissues. Highly perfused, relatively low volume tissues such as the brain equilibrate rapidly with the concentration of the central blood pool. Since the site of action of anesthetics and opioids is the central nervous system (CNS), this explains the rapid onset of effect of such drugs. Concurrently to this equilibrium, a quantity of the drug will also pass from the blood through other tissues such as muscles, viscera, fat and bones, hence slowly reducing the drug plasma concentration. As the concentration decreases, a part of the drug accumulated inside the highly perfused tissues will then be re-distributed inside the blood, hence decreasing the concentration of the drug in these tissues. The reversible transfer of the drug from one location to another is referred to as distribution, and involves molecular movement across lipid membranes and capillary walls. The distribution kinetics depend CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 34 then mostly on the lipid solubility characteristic of the drug. Another important aspect that affects the distribution phase of the drug is protein binding. While the bound fraction of a drug does not contribute to the pharmacologic action of the drug, it acts as a depot: as the free fraction of the drug is removed from the body due to metabolism and excretion, protein bindings are broken to reestablish the initial ratio between bound and free drug. The elimination of a drug is usually achieved through metabolism and excretion. Metabolism transforms lipophilic substances (such as most anesthetics and opioids) into hydrophilic substances, hence facilitating their excretion through renal mediation. This process, done at the hepatic level, usually leads to the inactivation of the drug. There can be a large individual variability in the activity of the enzymes responsible for metabolism. This variability can be as great as 40-times which might be the result of genetic differences between individuals. Renal and hepatic excretion are responsible for most of the drug elimination within the body. 3.1.2 Model ing The phenomenon of drug uptake, distribution and elimination can be expressed mathematically. To that end, a number of model architectures have been proposed. 3.1.2.1 Exponential Models For most drugs, the time course of their concentration within the blood plasma after rapid intravenous administration and uptake can be fitted to resemble a decaying function, with two distinct modes cor-responding to the distribution and elimination phase respectively, see Fig . 3.2. This behavior can be mathematically expressed as: Cp(t) = A • e~ at + B • e-P* [ng/ml] or [Mg/ml] (3.2) where Cp(t) is the time course of the drug concentration expressed in nanogram per milliliter (opioids) or microgram per milliliter (barbiturates, propofol), a is the rate constant of the distribution phase (slope of the tangent at the origin), and (3 is the rate constant of the elimination phase (slope of the curve once the distribution phase is completed). In many cases, a tri-exponential model will capture significantly better the kinetic of the drug 2: Cp(t) = P • e - ' - ' + A • e~ at + B • e - ^ ' [ng/ml] or [ug/m\] (3.3) A major advantage of exponential models is that they can be easily derived using graphical means. The identification can be carried out directly by using either bolus data and analyzing the decaying blood plasma characteristic (such as shown in Fig . 3.2), or by using infusion data and analyzing how the plasma concentration increases over time. 2 T h e use of /3 being usually reserved to model and represent the slow elimination phase, researchers have introduced the notation P and 7r to describe the fast dynamics corresponding to the distribution phase. CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 35 T T I M E A F T E R I N J E C T I O N Bolus administration Figure 3.2: Time course of plasma concentration with two distinct phases: distribution and elimination (adapted from [98]). Note that in pharmacology it is usual to replace the rate constants by half-lives which represent the amount of time that is necessary to observe a reduction of 50% of the initial dose during a given phase. Half-lives are calculated as: T f = ^ and r ? = — and i f = ^ [min] (3.4) 2 7T 2 Oi 2 p In terms of control and system engineering, the exponential model (3.3) can be directly expressed in the Laplace domain: PK(s) = ^ 4 = P - — + A- — + B--^— (3.5) v ; I(s) S + TT s + a s + /3 v / 3.1.2.2 Compartmental Models Pharmacologists usually consider exponential models to be counter-intuitive and prefer the use of mamillary compartmental models. Since the literature abounds in pharmacokinetic parameters given in the compart-mental framework, it is necessary to introduce these models here as well. Note that mamillary compart-mental models are used extensively in pharmacology, and are not restricted to modeling only anesthetic drugs. Pharmacokinetically, a compartment represents a tissue group that has similar kinetic characteristic. For instance, it is possible to consider that the body can be divided into three compartments: a small central compartment that contains the arterial blood and highly perfused tissues such as the brain and liver, a larger compartment which contains muscles and viscera, and a third compartment consisting mostly of fat and bones, see Fig . 3.3. Following the administration of a bolus, the amount of drug delivered into the central compartment is eliminated according to the rate constant fcio (usually expressed in m i n - 1 ) . This elimination corresponds to metabolism and hepatic and/or renal excretion. In parallel to the elimination CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 36 Rapidly equilibrating compartment Slowly equilibrating compartment Muscles Viscera Elimination Figure 3.3: 2-compartment pharmacokinetic model. A third compartment is often added for improved accuracy (adapted from [99]). process, the drug is distributed in the two peripheral compartments at a rate of k\2 and k\^. Following the administration of a drug, the concentration C\ of the first compartment decreases rapidly while the concentrations C2 and C3 of the second and third compartments rise. Once the concentrations in the central compartment and any of the peripheral compartments reach equilibrium, the distributive process reverses and the drug stored in the peripheral compartment returns to the central compartment at the rate of k2i or & 3 i . Since the blood of the central compartment acts as a carrier for the drug, we can assume that there is no direct drug exchange between the two peripheral compartments. For a similar reason, only the drug presents in the central compartment can be eliminated. The mathematical expressions governing this model can be obtained by writing the mass balance equa-tions (expressed here in a state space representation): = C 3 (*). 10 Cp(t) 1 0 0 -&12 + ki3 Ci{t) C3(t\ k2i k31 Ci(«) -&21 0 C2V) 0 -*3l. C 3 ( t ) . + l_ 0 0 i(t) (3.6) where V i is the volume of the central compartment. Note that, by definition, the plasma blood concentration equals the drug concentration of the central compartment, i.e., Cp(t) = C\(t). V o l u m e o f D i s t r i b u t i o n The volume of distribution Vj is a measurement of the extent of the distributive process. It is defined as the ratio between the total amount of the drug in the body versus the drug plasma concentration prior to the beginning of the elimination phase. The volume of distribution does not have any physiological significance, and it is also referred to as the apparent volume of distribution. CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 37 The volume of distribution is dependent on the size of the tissues into which the drug distributes, the partition coefficient of the drug, and the extent of protein binding. This volume can be extremely large, even exceeding the anatomic size of the body, meaning that most of the drug is distributed inside peripheral tissues. For a 3-compartment model, the volume of distribution is simply the sum of all three compartment volumes: Vd = V1 + V2 + V3 [1] (3.7) where V2 and V 3 are the volumes of the rapidly equilibrating and slowly equilibrating compartments. Clearance Another important concept used by pharmacologists is the intercompartmental and total drug clearance concept. It has been mentioned earlier that the elimination of a drug from the body is achieved through a variety of mechanisms, involving metabolism and excretion. Total drug clearance can thus be thought of as the ability of the body to remove a given drug from the blood or plasma. Clearance is usually expressed in [ml-h - 1 ] . Using a 2- or 3-compartment model, the total body clearance Cl\ can be calculated as: C h = Vi • kw [ml-h - 1] (3.8) The intercompartmental clearances Cli2, Cl2\, and C/13, CZ31 represent a similar concept. Instead of the elimination process, they quantify the ability of each compartment to exchange the drug with one another. Intercompartmental clearances are not very useful. However, considering that Cl\2 = Cl2\ = Cl2 and C Z 1 3 = CI31 = CI3, and that CUj = Vi • kij, we reach the following equalities: V2_ku V3_ki3 Vi " k21 V, ~ *si [ ' 3.1.2.3 Equiva lence Be tween E x p o n e n t i a l and C o m p a r t m e n t a l M o d e l s The models expressed in (3.3) and (3.6) are essentially equivalent as they describe the same input/output relationship. Because of this equivalence, a P K parameter set can take different forms: - {P, A, B, TJ, Ti, Ti}: standard exponential parameter set, 2 2 2 - {Vi,P*, A*,B*,TJ,T?,Tf }: exponential parameter set with normalized partial coefficients (i.e., 2 2 2 P* + A* + B* = 1), - { V i , kio, kn, fa, ^ 2 1 , fei}: standard compartmental parameter set, - { V i , V2, Vs,Cli,Cl2,CI3}: compartmental parameter set expressed as volumes and clearances. A l l of these sets are equivalent as they define the same P K relationship. Most of the P K parameter sets found in the literature are expressed in either of these forms. CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 38 In order to simplify the notations, we propose to express the P K model as a SISO transfer function PK{s) using both the exponential and compartmental parameters: r K [ S ) ~ I(s) " Vi " (s + TT) • (s + a) • (s + /?) (i-W> This hybrid notation has the advantage of clearly expressing each one of the system modes, hence allowing a trained observer to readily identify the frequency response of the system. 3.1.3 L i t e r a t u r e R e v i e w Medical literature is quite prolific in reviewing drugs and their P K properties. In this section, we present a summary of a survey of published P K parameter sets for propofol and remifentanil. The complete survey is presented in the Annex B . Each model is summarized under the hybrid form of (3.10) in Table B . l and B.3. P r o p o f o l Propofol has received considerable attention over the past 20 years. A large number of phar-macokinetic models have been published by various authors, and under various conditions and methods of identification. In the mid-1990s, researchers started to use the N O N M E M 3 software ( N O N M E M Project, University of California, San Francisco, C A ) to derive the P K models from data sampled across a large population of patients. This software package allows the user to identify the most significant covariates and include them in the final model. In terms of propofol, the age, weight and Lean Body Mass ( L B M ) of patients were found to alter the P K parameters. The sampling site (venous vs. arterial) and the method of administration (bolus vs. infusion) were also significant covariates. Schnider et al. [100], and later on Schiittler and Ihmsen [101], used this technique in their analysis. As of today, Schiittler and Ihmsen P K sets are the most up-to-date. In their study (spanning hundreds of patients and thousands of blood samples) the intra-patient variability was found to be less than 20%. The mean absolute weighted residual 4 ( M A W R ) was about 25% and the mean weighted residual ( M W R ) was -3.4%. The P K parameter set is presented in the Table 3.1 and Table 3.2. R e m i f e n t a n i l Remifentanil being a relatively new drug, the derivation of its P K model has directly benefitted from the N O N M E M analysis software. As such, the models published by Egan et al. in 1996 [102], and later on by Minto et al. [103], are consistent and accurate (the M A W R was less than 25% in both studies). As a result, no further attempts at modeling this drug has been carried out since then. Minto's P K parameter sets (see Table 3.3) have the advantage of accounting for the age and lean body mass of the patients. However, they can only be used to predict arterial blood concentration during N O N M E M : NON-linear Mixed Effect Models 4 T h e M A W R is a standard PK accuracy parameter. It is defined as M A W R = X!»=i N 52J=I ..M ™—?P?A> where ctj is Cpij the j t h measured plasma concentration of the i t h individual, and cpij denotes the corresponding predicted value. N is the total number of individuals and M is the number of samples. The M W R is a similar accuracy parameter without the absolute value. CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 39 P K P A R A M E T E R V A L U E UNITS CI i (BW/70)^ if age < 60 [l-min"1] 01 (BW/70)^ - (age - 60) • 0W if age > 60 P-min"1] Cl2 03 (BW/70) 9 s • (1 + ven • t?14) • (1 + bol • 0 1 6) [1-min-1] Ch 05 (BW/70) e " • (1 + bol • 618) [1-min-1] VL 02 (BW/70) e » • (age/30)e i3 • (1 +bol • 015) [1] v 2 04 (BW/70) e 9 • (1 + bol • 0i 7) W v 3 06 [1] Table 3.1: Propofol P K parameter sets from [101] where BW stands for body weight, ven= 1 for venous sampling, ven= 0 for arterial sampling, bol= 1 for bolus administration, and bol= 0 for infusion administration. P A R A M E T E R ESTI-MATES V A L U E S UNITS S E 0i 1.44 [1/min] 0.09 0i 9.3 W 0.9 03 2.25 [1/min] 0.31 04 44.2 M 6.1 05 0.92 P/min] 0.15 06 266 [1] 43 07 0.75 0.06 08 0.62 0.09 09 0.61 0.11 010 0.045 0.012 011 0.55 0.13 012 0.71 0.26 013 -0.39 0.15 014 -0.40 0.10 015 1.61 0.36 016 2.02 0.41 #17 0.73 0.23 #18 -0.48 0.12 Table 3.2: Parameter estimates from the N O N M E M analysis ([101]). CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 40 P K P A R A M E T E R V A L U E UNITS Cl2 Ch 2.2 - 0.0162 • age + 0.0191 • L B M 3.25 - 0.0301 • age 0.121 - 0.00113 • age [1-min-1] [1-min-1] [1-min"1] V-2 V 3 1.94 - 0.0201 • age + 0.072 • L B M 7.12 - 0.0811 • age + 0.108 - L B M 5.42 [1] [1] [1] Table 3.3: Pharmacokinetic parameters of remifentanil as a function of the age and lean body mass (from [103]). remifentanil infusion. This is hardly a limitation since infusion is the method of administration of choice for this drug. Also, due to the drug susceptibility to blood esterases, venous blood concentration significantly lags the arterial concentration [104]. As such, P K models based on venous blood sampling should not be used as a basis to predict drug effect. 3.1.4 Concluding Remarks 3.1.4.1 Drug Uptake In P K analysis and identification, it is usually assumed that the drug mixes instantaneously within the central compartment upon administration. However, this assumption is invalid when considering the kinetics of the drug during the uptake phase, i.e., right after administration. Since the uptake phase is extremely short, this phenomenon is difficult to measure properly as it would require numerous blood samples taken in the first minute after administration. However, some authors [105] consider that a time delay of a few seconds should be accounted for in the model to describe the initial drug uptake phase. 3.1.4.2 Linearity Drugs pharmacokinetics are not linear in the sense that the rate of injection of the drug affects the overall drug distribution. This aspect of pharmacokinetics is well known by clinicians and practitioners who do not titrate anesthetic drugs the same way whether they use boluses or infusions. In that respect, most of the recent P K studies make a clear distinction between P K models obtained from bolus data and P K models obtained from infusions. Note for example in F ig . 3.4 the two-times difference between the propofol impulse response of the infusion vs. bolus models (from Schiittler [101]). As expected, this difference occurs mostly within the first few minutes following the administration. One of the only study addressing this issue was the one published in 1998 by Schnider et al. [100] for propofol. The authors observed that propofol kinetics are linear for an infusion range of 25 to 200 /xg-kg" 1 min CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 41 T I M E [min] FREQUENCY [rad/s] Figure 3.4: Bolus vs. infusion pharmacokinetic of propofol (using [101]). (a) Impulse response (b) Frequency response. 3.1.4.3 Time Variance and Other Factors Affecting PK Models There are a number of physiological differences between individuals that can affect pharmacokinetic pa-rameters. Understanding and quantifying these differences can reduce the model uncertainty. For instance, most pharmacokinetic parameters are derived for a given age population. Pediatric, adult or elderly patients can have very different hepatic and renal extraction ratio, thus affecting the elimination of the drug. Also, since the distribution of a lipophilic drug is affected by the presence of adipose tissues, obesity in patients should be accounted for, at it increases the elimination half-life of the drug. Hence, weight, fat and lean body mass are often cited as covariates of pharmacokinetic parameters. Finally, in recent years, it has been shown that cardiac output also affects drugs kinetics ([106], [107]). For instance, it has been shown that the initial arterial concentration of propofol, following its rapid infusion in sheep, is inversely related to the cardiac output. This issue hasn't yet been widely documented. The question whether or not the inclusion of the cardiac output in the pharmacokinetic parameters would reduce the uncertainty in the model remains opened. A preliminary study by Adachi et al. [108] tends to indicate that pharmacokinetic variability can be reduced when accounting for the cardiac output. However, since cardiac output can only be measured by invasive means, this issue hasn't yet received much attention. While age, lean body mass, and eventually cardiac output are covariates of pharmacokinetic parameters in normal healthy subjects, factors such as renal dysfunction and hepatic diseases will also bear a tremendous impact on drugs kinetics, mostly during the elimination phase. We invite interested readers to refer to [98] for a more in-depth analysis on how these diseases can affect total drug clearance. 3.1.4.4 Pharmacokinetic Interactions Effect of opioids administration on propofol A limited number of studies have documented the pharmacokinetic interactions between propofol and opioids. A first study by Briggs et al. in 1985 [109] CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 42 concluded that fentanyl significantly increases the level of the propofol blood concentration profile. To reach this conclusion, they compared the parameters of a 3-compartment model of propofol pharmacokinetics obtained from two different groups of patients. In the first group, patients received a single bolus of propofol. In the second group, a bolus of fentanyl was administered as a premedicant before the injection of propofol. They observed that the apparent volume of distribution, as well as the central compartment volume, were reduced by 45% in the group that received fentanyl. The total drug clearance was also significantly reduced. A similar study by G i l l et al. in 1990 [110] yielded a very different conclusion. Their experimental protocol also involved a patient population divided into two groups (one receiving only propofol, and the other receiving both propofol and fentanyl). The authors reported no significant differences in the pharmacokinetic profile of propofol between the two groups. As of today, the effect of the administration of an opioid on propofol pharmacokinetics remain unclear and a subject of debate. Effect of propofol administration on opioids Conversely, studies by Gepts et al. [ I l l ] and, more recently, Mertens et al. [112] in 2001 reported that the administration of propofol significantly reduces the elimination and inter-compartmental clearances of alfentanil. The authors also mentioned that a reason for this interaction could be linked to the fact that propofol tends to reduce the mean arterial blood pressure, thus affecting opioid distribution and elimination. Further studies are necessary to evaluate the precise mechanisms that can cause these interactions [24]. 3.2 Pharmacodynamics The role of pharmacodynamic modeling is to mathematically express the observed effect of a drug as a function of its plasma concentration: ( 3 ' n ) where PD(s) is the pharmacodynamic model and E(s) is the drug effect. Note that for a given drug there can be a multiplicity of effects. Hence, pharmacodynamic models are not unique as they describe only one over the many different possible endpoints of a drug. 3.2.1 The Dose/Response Relationship: a Steady-State Model The characteristic of any given drug can be expressed by dose-response curves, see Fig . 3.5. As expected, at low doses, drugs have only a minimal effect. As the dose increases, the effect increases as well. Often, dose-response curves exhibit a linear characteristic (constant slope) for a given range of doses. A t higher doses, a saturation phenomenon appears, indicating that the full effect of the drug has been reached. The CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 43 Efficacy J Individual Variability D O S E OF D R U G Figure 3.5: Dose-response relationships differ according to the drug potency, efficacy, slope and subjects' variability, sigmoid characteristic of the dose/response is usually captured using the Hi l l equation: where E is the steady state effect, Cp is the steady state plasma concentration, and EC50 is the concentration which yields 50% of the maximal effect, or which yields the maximal effect in 50% of a given population of patients. The parameter 7 is a measure of the steepness of the dose-response curve. Finally, Eo and Emax are the minimum and maximum effects. When considering anesthetics and opioids, it can be difficult to quantify drug effects. Very often, only quantal (e.g., response vs. no response) observations can be made. For example, in the case of opioids, the absence of reaction (i.e., movement) to intubation, skin incision, and skin closure is often taken as a relevant endpoint (see Fig . 3.6). The results obtained from different patients are then pooled in order to obtain the sigmoid-like shape which represents the probability of response versus the administered dose or the drug plasma concentration. In this case, the dose/response curve can be used to determine the dose or plasma concentration which would yield the desired effect (e.g., no response to intubation or skin incision) in a given percentage of a population. This type of dose/response curve are the most commonly found in the literature preceding the mid-1990s. Conversely, it is sometime possible to target an effect which can be quantifiable. The resulting dose/response curve would then indicate the dose or plasma concentration which are necessary to obtain a partial effect (e.g., the plasma concentration of propofol to reach a certain BIS value). Most of the recent research done in pharmacodynamic modeling seems to follow this new trend. When a number of effects or observational endpoints can be selected to model a drug, it is important to carefully select the effect or endpoint which yields a useable characteristic with respect to the therapeutic window of the drug. For instance, effects leading to dose/response characteristics such as in F ig . 3.7.a-b are not suitable since they are either on/off or they fail to characterize the whole therapeutic range of the drug. A n ideal endpoint for modeling the dose/response characteristic would be an effect leading to the linear characteristic of Fig . 3.7.c. This characteristic simply implies that a variation of the effect can be E = E0 + E, 'max EClJ+Cl (3.12) C H A P T E R 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL A P P R O A C H 44 o r~ Intubation Skin Incision 200 1 No response Response 400 600 MIHM IIMl B ml 111 M l . , , W " | TTTn Skin Closure . " [ » * • [ « • g I I I to \ 0 O °1 § g 50, " I o 2 £ 25, No response Response No response Response P L A S M A A L F E N T A N I L |ng/ml] Figure 3.6: Drug dose-response curves can be obtained by observing whether a patient responds to a particular stimulus at a given concentration (from [113]). Figure 3.7: Different dose/response characteristics (a) On/Off effect (b) The observed endpoint cannot characterize the effect of small doses (c) Ideal characteristic. CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 45 EFFECT SITE ALFENTANIL [ng/ml] Figure 3.8: Dose/response model for Alfentanil and using the probability of maximal E E G response (from [31]). achieved by applying a proportional variation in the drug dosing regimen, independently of the operating point of the system. To illustrate this discussion, let us consider the use of the E E G in measuring the effect of opioid agents, see F ig . 3.8. In this example, the smallest doses used in cardiac surgeries already induce a maximal E E G response (pronounced delta waves). Hence, this particular endpoint is useless in cardiac surgeries. However, it might be appropriate in non-cardiac cases where the opioid doses are generally much lower. 3.2.2 The keo Parameter: a Measure of the Effect Dynamics Drug-response relationships are obtained based on steady state observations. As such, they do not capture the dynamic behavior that exists between the drug plasma concentration and the observed effect5. The effect of rapid intravenous administration is well known and well documented. When plotting the observed effect versus the plasma concentration, a hysteresis-like cycle can be observed (see for instance Fig. 3.9) indicating a disequilibrium between the plasma concentration and the effect site concentration of the drug. A direct consequence of this disequilibrium is that the effect lags the plasma concentration. A n explanation for this disequilibrium is that the drug concentration within the biophase6 lags the drug concentration in the blood. Hence, predicting the effect site concentration Ce is relevant in order to devise a proper infusion profile. In 1979, Sheiner et al. [114] proposed an addition to the conventional pharmacokinetic models in the form of an effect compartment connected to the central compartment. This hypothetical effect compartment models the temporal aspect of pharmacodynamics where a rate constant keo expresses the dynamics of transfer of the plasma concentration to the biophase. Therefore, it is possible to mathematically express the effect site drug concentration C e (s) as a function of the plasma concentration Cp(s) as: s + ke0 5Often in the clinical literature, pharmacodynamics models are limited to only the static Hill dose vs. response relationship. °The biophase is defined as the physiological milieu targeted by the drug. It usually cannot be accessed for sampling. CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 46 F E N T A N Y L SPECTRAL SPECTRAL [ng/ml] E D G E [HZ] E D G E [HZ] Figure 3.9: Effect and plasma concentration of fentanyl and alfentanil following rapid intravenous administration (adapted from Stanski [89]). This equation only captures the effect dynamic. It has to be used in conjunction with the non-linear gain of the steady-state model (3.12): E(s) = E 0 + Emax • F r 7 C f f r , 7 (3.14) This modeling approach to pharmacodynamics has become the mainstream approach followed by phar-macologists since the early 1980s. It will be clear by anyone well-versed in System Engineering that this approach suffers from a number of limitations. In particular, the model (3.13) was never rigorously validated and tested using proper identification data. For instance, the absence of a time delay in (3.13) may hide the inherent latency of the drug effect. In addition, the dynamics of the sensor used to measure the effect is not included in the model. Finally, the use of the Hi l l saturation characteristics as the primary modeling tool cannot be recommended, as it can falsely characterize linear dynamics as non-linear. These limitations are further addressed in Chapter 6, where a new modeling approach based on System Identification know-how is proposed. 3.2.3 Literature Review 3.2.3.1 Propofol Only 6 studies involving quantitative E E G have been published since the mid-1990s. These studies involve a limited number of volunteers and/or patients. In terms of the equilibration time constant keo, results vary from 0.1 to >0.6 m i n - 1 . This large variation can be attributed to differences in the type of identification methodology used in the studies. Also, the fact CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 47 Figure 3.10: Individual (dashed lines) and mean (solid lines) Hill non-linear characteristics of published pharmaco-dynamic models, (a) Propofol study from Kuizenga et al. [118] (b) Remifentanil study from Egan et al. [102], that each study was using a different P K set to calculate the blood plasma concentration, or used their own blood measurements, probably resulted in additional variability. The EC50 parameter is also significantly different between studies. Using auditory evoked potentials, White et al. [115] found that a concentration of 2.17 mg/ml was sufficient to provoke unconsciousness. A similar result (about 2 mg/ml) was reported by Schnider at al. [100] for their semilinear canonical correlation index. Much higher values (about 6 mg/ml) were reported by Kazama et al. [116] for the bispectral index 7. Conversely, Billard et al. [117] and Kuizenga et al. [118] reported intermediate values (3 to 4 mg/ml) for BIS, which are closer to the usual recommended concentration suitable for surgery. Probably the most remarkable limitation of the P D models presented in some of these studies is the rather large inter-individual variability which can be observed from the Hi l l non-linear characteristic. For instance, Kuizenga et al. published the Hi l l parameters for each one of their subjects, see F ig 3.10.a. According to the Hi l l relationship, some of their subjects experienced a complete depression of their cortical activity at concentrations that barely had any effects in other subjects. The averaged Hi l l relationship might thus be useless in the face of so much variability. The only study which reported a reasonable P D variability is the one by Kazama et al. [116]. One of the reason for the consistency of their results was that they compensated for the dynamics of the BIS sensor and removed the initial time delay between the drug administration and the start of the BIS descent. While this technique had the advantage of leading towards a more consistent P D model, the fact that part of the system dynamics was removed from the identification data limits the usability of the proposed models. However, one very interesting conclusion of the study was that age is a P D covariate. However, conversely to many other drugs, they found that the EC50 of propofol tends to increase with age, indicating that older 'Accordingly, a concentration of 6 mg/ml is necessary to reach a BIS of 50. CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 48 F E N T A N Y L C O N C E N T R A T I O N [ng/ml] F E N T A N Y L C O N C E N T R A T I O N [ng/ml] Figure 3.11: Drug interaction between propofol and fentanyl (adapted from [120]). patients are less sensitive to propofol than younger patients. This result, while surprising, was corroborated by another study on rats [119]. They also found that the average EC50 in the adult population was about 6 mg/ml, which is much higher than most of the other studies. R e m i f e n t a n i l The effect of remifentanil on the spectral edge frequency has been studied by both Egan et al. [102] and Minto et al. [103] in the late 1990s. Both studies also reported a considerable variability in the Hi l l parameters (see for instance Fig. 3.10.b). This variability could be due in part to the fact that the P D modeling was carried out using spectral edge frequency which is less sensitive to opioids than to anesthetics. In some cases, only a marginal decrease in the spectral edge could be observed. 3.2.4 C o n c l u d i n g R e m a r k s 3.2.4.1 P h a r m a c o d y n a m i c Interact ions Propofol and opioids are known to interact with each other in a synergistic fashion in their anesthetic/analgesic effect, but not in their toxicity. This observation constitutes the basic assumption on which the concept of balanced anesthesia is based. In an interesting set of studies, Smith et al. [120] have measured the reduction of propofol and fentanyl effective dose when both drugs were used concomitantly. Their results are displayed in F ig . 3.11. Quite clearly, there is an optimal titration that ensures adequate hypnosis and analgesia. For instance, the effect of the combined use of propofol and, in this case, fentanyl shows that the effective concentration EC50 needed to suppress movement in response to a noxious stimulus in 50% of patients can be drastically reduced. Where a concentration of 12 / i g - m l - 1 of propofol (or 6 n g - m l - 1 of fentanyl) was necessary, using both drugs simultaneously reduces this requirement to 4 /xg-ml - 1 of propofol combined with 1.5 n g - m l - 1 of fentanyl. CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 49 60 min infusion Figure 3.12: Simulation of propofol and remifentanil pharmacodynamic interactions (adapted from [23]). Synergism is also apparent in hypnosis, but to a lesser extent. Pharmacodynamic interactions between propofol and remifentanil have also been studied [23]. Surface plots were generated to estimate the optimal titration that minimizes the recovery time upon termination of the infusion, see Fig. 3.12. The curve in the lower plane represents the EC50 of both propofol and remifen-tanil. At the extremes, a concentration of either 12 yug-ml - 1 of propofol or 12.44 ng-ml""1 of remifentanil is required to provide adequate anesthesia (in this case: no movement to painful stimuli). In this simulation the infusion rate of both drugs is maintained constant. After 1 hour, the infusion is discontinued, resulting in a decrease in both drugs plasma concentrations. The thick line is the isobole that predicts emergence in 50% of the patients. A minima at EC^0 = 2.5 ^ g - m l - 1 (propofol) and ECl0 = 5 n g - m l - 1 (remifentanil) ex-its, indicating that this combination of drugs would optimize the emergence time, while providing adequate titration for the maintenance of anesthesia. Of interest, a 2000 study by Minto et al. [25] has extended this concept to 3 drugs (propofol, midazolam and alfentanil). While surface plots clearly stress out the synergism between drugs, their interest is limited to steady state. No dynamic model of pharmacokinetic/pharmacodynamic drug interactions have been published so far. This represents a very fertile area for medical research. 3.2.4.2 Time Variance Propofol In a recent study from Kuizenga [118], the time invariance of propofol pharmacodynamic models has been seriously questioned. Three propofol infusions were given to 10 patients. Each successive infusion CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 50 T I M E [min] Figure 3.13: Measured and predicted time courses of the BIS following a series of induction sequence administered to a volunteer (from [118]). A pharmacodynamic model obtained based on the first repeat was used to predict the BIS time course. Note the good agreement between measurement and prediction for the first repeat. Results worsen considerably for the two other repeats. was started immediately after the patient regained consciousness. A set of pharmacodynamic parameters were obtained from the first infusion data. This set was used to predict the effect of the subsequent infusions based on the plasma concentration. The results have shown an important discrepancy between the predicted effect and the observed effect, see Fig. 3.13. The authors concluded that pharmacodynamic models cannot be used accurately over time to predict the effect of propofol. In their study, they analyzed the effect of propofol on both the bispectral index and the amplitude of the patients' electroencephalogram in the 11 to 15 Hz band. A second set of pharmacodynamic parameters were derived based on the observation on the two successive infusions. Both sets are presented in the Table B.2. Note the large discrepancy between the first and the second set. R e m i f e n t a n i l The long term administration of opioids usually results in the development of tolerance to analgesia followed by a strong physical dependence to the drug. In the case of remifentanil, a study by Vinik and Kissin [121] involving 13 volunteers have shown that tolerance to analgesia is profound and develops very rapidly. A constant-rate infusion of 0.1 ^ g k g _ 1 - m i n _ 1 was administered over a period of 4 hours. The analgesic effect was evaluated by measuring pain tolerance to thermal (cold water) and mechanical (pressure) noxious stimuli. Tolerance to pain was measured as the length of time that the subject was able to tolerate the stimulation. A baseline value was recorded before the start of the infusion. Painful stimulation was applied every 30 minutes. The study has shown that tolerance to pain is maximum after a 60 to 90 min infusion. Once the peak is reached, analgesia begins to decline steadily, see Fig. 3.14. After CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 51 o H 70 60 50 40 30 20 10 0 mean +SD -30 _ l _ _1_ 60 90 120 150 TIME [min] 180 210 240 Figure 3.14: Development of acute tolerance to analgesia during short-term constant-rate infusion of remifentanil (adapted from [121]). 4 hours, the effect of remifentanil was no longer observable. The authors of the study concluded that calculations for target-controlled infusion of remifentanil must include corrections for tolerance. However, this would require models which have not yet been developed. 3.3 Summary Since its inception in the early 1950s by Dost, the study of pharmacokinetics has received considerable attention from the research community. Recent advances in non-linear regression analysis have allowed pharmacologists to define P K parameter sets which account for the effect of covariates such as age, weight, gender and ethnicity. It is now possible to predict with a fair accuracy the blood plasma concentration of most of the anesthetic drugs in use today. In terms of control engineering, there are three important aspects in pharmacokinetics, which require some attention: - P K parameters are dependent on the infusion rate of the drug. For instance, slow infusions (less than 200 ^ g - k g _ 1 - m i n _ 1 ) are linear with respect to the input amplitude, while bolus-type infusions (above 2500 /ng-kg - 1 -min - 1 ) have a different initial drug distribution, which results in a different P K parameters set. - A second aspect is the interaction between drugs, and in particular propofol and remifentanil. While the literature is very confused about this issue, most researchers tend to agree that, if any interac-tions occur at the level of pharmacokinetics, these interaction mostly affect the re-distribution and elimination of the drugs. - Finally, the blood plasma concentration of intravenous drugs cannot be observed nor controlled directly since there is no available real-time measure. The pharmacokinetics of intravenous drugs can thus only be used to predict the plasma concentration. Accounting for patient variability appears then necessary. CHAPTER 3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL APPROACH 52 Patients with hepatic or renal diseases may have a widely different P K characteristic. The P K models presented in this section should not be used in these instances. Since Schiittler and Ihmsen in 2000, no further studies were conducted to derive the P K model of propofol. While each one of the reported studies (see Section 3.1.3) have contributed in improving the goodness of the fit (mostly by adding the age, weight and lean body mass as covariates), it was argued that no further improvements can be made. The error between the measured plasma concentration and the predicted plasma concentration is now about the same when considering population-normed P K models or individualized P K models (i.e., P K models derived specifically for one single individual). This indicates that further improvements can only be achieved by considering more complex physiological models instead of a mamillary compartmental model. While dose-response curves have been used intensively in the past to derive dosing guidelines (one of the brightest example is the M A C concept for gaseous anesthetics), it is only recently that P D models based on quantifiable endpoints such as EEG-based parameters have been researched and published in the literature. Unfortunately, these P D models do not seem to be appropriate for use in control-oriented frameworks. For instance, and conversely to P K models now used extensively in Target Controlled Infusion (TCI) pumps in Europe, P D models have not yet been used successfully in any clinical application. It is likely that the large variability observed between patients and the discrepancies between the model parameters published in the literature explain the lack of enthusiasm in the anesthesia community for P D models. However, we believe that the poor performance of these models stems mostly from a poor choice of the model structure rather than an inherent limitation of the system. In Chapter 6, we revisit in details the problem of P D modeling and identification. In particular, we propose a new system-oriented approach that yields better prediction performance and a more accurate description of the frequency response of the system dynamics. Chapter 4 Closing the Loop in Anesthesia: a Review In terms of control engineering, the idea of the anesthesiologists striving to keep a balance between the pharmacological toxicity of anesthetic drugs and the effect of surgical noxious stimuli is particularly ap-pealing. The anesthetic and opioid titration needs to be constantly adjusted in order to avoid both under-and overdosing. The idea of an automated system that would regulate drug titration in order to maintain the adequacy of the anesthetic regimen is immediate. The concept of a closed-loop anesthesia system is thus very close to that of automated flight control. The role of the anesthesiologist is similar to that of a pilot: after take-off (induction), the pilot usually maintains an adequate flight trajectory (hypnosis, analgesia, paralysis). Nowadays such tasks are performed by flight controllers able to plan ahead, optimize fuel consumption and minimize the duration of the flight. Closed-loop anesthesia is somehow similar in the sense that, by changing the titration of intravenous drugs, the anesthesiologist can drive the patient into a deeper or lighter hypnotic ahd/or analgesic state, according to the requirements of the surgical procedure. Using proper feedback quantities and drug models the possibility to automate the drug titration and to allow the practitioner to concentrate on higher level tasks seems viable. In keeping the comparison, closed-loop anesthesia would not replace the anesthesiologist. On the contrary, the workload of the anesthesiologist will be reduced during the maintenance stage, leaving room for full attention in case of an emergency, or when the override of the controller is required. 4.1 Review The concept and implementation of closed-loop anesthesia has been investigated for the past half cen-tury via numerous attempts at controlling anesthetic drugs titration through feedback control. A selected chronological survey of this prior art is summarized in Table 4.1 and 4.2. For each mentioned attempt, we listed the feedback quantity, the drug used, and the control technique employed. In spite of a good number 53 CHAPTER 4. CLOSING THE LOOP IN ANESTHESIA: A REVIEW 54 of such attempts, no clinically satisfying results have been obtained so far. This survey is also limited to attempts which focused more precisely on the control of hypnosis. Other work focusing on the control of heart rate and blood pressure, as well as the control of muscle relaxation in the context of anesthesia has not been reported here. In recent years, researchers have been using either a simple PID or a lookup table of the drug pharma-codynamic model to set the target plasma concentration of a target controlled infusion devices in order to reach and maintain a given hypnotic reference. The successes reported by Struys et al. [122] and Absalom et al. [123] can be easily attributed to the fact that their system titrated the drug according to the index of consciousness provided by the bispectral monitor, rather than to the performance of the closed-loop controller itself. The results produced by controllers embedding advanced techniques, as shown in the work by Frei et al. [124] and Gentilini et al. [125], emphasize that the problem is far from being solved due to the aforementioned challenges posed by the intra- and inter-patient variability. There have also been attempts at closing the loop by Linkens et al. ([126], [127], [128], [129]) using a variety of intelligent control techniques such as expert systems and fuzzy logic. Linkens et al. were probably among the first to attempt the control of distinct anesthesia components simultaneously (analgesia and areflexia) using different agents (atracurium and isoflurane). A n in-depth analysis of such cases reveals the need for strong knowledge of the patient model. The intra- and inter-patient variability makes the establishment of a priori rules very difficult. Stability of fuzzy-based controllers may be impossible to prove. From the perspective of interaction between drugs, and of particular interest, is the attempt by Zhang et al. in 1998 [130] at controlling an intravenous anesthetic (propofol) together with an opioid (fentanyl). This approach was limited to the control of the plasma concentration of propofol and fentanyl in dogs, where the setpoints were chosen to minimize the wake up time. Probably the most interesting results obtained in recent years were those of Locher et al. and L i u et al, see Table 4.2. Not only have they shown that closed-loop control is particularly effective in maintaining a proper BIS setpoint during anesthesia, they have also proven that closed-loop control brings further clinical advantages in terms of drug consumption and wake up time, as compared to manual titration [142]. Closed-loop controllers also react faster to the abrupt changes brought by surgical stimulation [141]. 4.2 Summary We have seen in Chapter 2 that the use of the new monitors of hypnosis, such as the BIS monitor, results in faster wake up, fewer side effects, and a reduction in drug usage, while ensuring that patients remain unaware during their surgery. For many years, however, the advantage of closed-loop system based on these monitors remained unclear. In particular, while closed-loop systems inherit de facto the advantages brought by the feedback sensors, doubt still persists whether there is any further clinical advantage to having a fully closed-loop system as compared to manually adjusting the titration based on the sensor information [143]. STUDY FEEDBACK QUANTTTY(IES) CONTROLLED AGENT(S) CONTROL TECHNIQUE POPULATION COMMENTS AND LIMITATIONS Bickford et al, 1950-1960 [38] [39] [40] EEG energy in D . ., , ,, . , , ° ,, Pentothal or — , , the 4 to 12 Hz _, . On/Off control , , Ether ' band Rabbits, cats, monkeys. _ ... ,. , ± , , , , . ™. . , ^ . , -„ Osculations due to the control technique. Clinical trial on 50 w , ... . . . „ atients under oin Method sensitive to extraneous interferences, pa len s un ergomg ^ Qp^j^g n a V e been administered concurrently, thus . various abdominal . , .. . . . ,, . , , . . . , , surgeries seriously limiting this technique m today s practice. Bellville et al., 1955-1960 [41] EEG amplitude in a denned Cyclopropane Analog control (P or PI?) frequency band No results have been presented. The proposed technique is merely an improvement of Not Disclosed Bickford's servo anesthesia. The authors mention speed (infusion pump) and position (vaporizer) control, however no specific information are disclosed. Method limited to the control of a unique anesthetic agent. Schwilden et al, 1985-1995 [44] [45] [131] [132] Methohexital Model-Based Adaptive Control Median Edge Propofol . ,. , . c ,, p ° ^ Adaptation was done it the requency system output was diverging too far from its reference Alfentanil Isoflurane 13 volunteers (22-29 yr ; Results on volunteers have shown that a constant excitation 44-85 kg) is necessary the guarantee the reliability of the feedback quantity (otherwise the volunteers were drifting from a 11 volunteers (24-31 yr ; drug-induced unconsciousness into a natural sleep). 54-87 kg) This technique works also for opioids. The controlled drug was used as the only anesthetic agent 11 patients during the maintenance phase. 25 female patients (31-47 yr), ASA I or II Kenny et al., 1990-1995 [133] [134] [135] „ , , Outer PI controller setting the Auditory Evoked n e i c ms-n „ . , Propofol reference to an inner TCI Potentials , device. The authors recommended feedback control of anesthesia as a 27 patients research tool to better identify pharmacodynamic models and the interaction between drugs. Roy et al, 1995-2000 [136] [64] Halothane Fuzzy rule-based control Auditory Evoked system controlling either the Potentials Propofol vaporizer or giving a reference to a TCI device 10 sessions conducted on 6 mongrel dogs with tail These papers emphasize mostly the hypnosis index derived clamping stimulation from midlatency auditory potentials using wavelet analysis. . Due to the extensive averaging needed, a value quantifying 9 sessions conducted on 5 t h e [ e y e l o f h y p n o s i s w a s calculated every 3 minutes, mongrel dogs Gent.ilini et al, 1995-2000 [124] (125] [18] [137] [138] Mean Arterial , . „ ,. , . _ Model Predictive Control Pressure Cascade Internal Model , a Control. An inner loop Isonurane , - , , , T . controls the end-tidal Bispectral index . .. m, r c concentration. 1 he reference of the inner loop is set by an outer loop that regulates BIS variations 20 atients The control of inhalational gas such as isonurane has the advantage that the drug plasma concentration is closely related to the end-tidal expired gas concentration which is readily available. This is a clear advantage over intravenous agents for which measurement of drug concentration is impractical. The authors have shown that isonurane can both control the patients ( - yr) mean arterial blood pressure (MAP - sensitive to noxious stimuli) and the BIS (sensitive to the hypnotic level). The authors proposed to investigate a SIMO controller where isonurane would be used to control both the MAP and the BIS I o I O O PI to I 3 Table 4.1: Prior Art: a Literature Survey (Part I) STUDY FEEDBACK QUANTETY(IES) CONTROLLED AGENT(S) CONTROL TECHNIQUE POPULATION COMMENTS AND LIMITATIONS S t r u y s et al., 2001 [122] Bispectral Index Propofol A lookup table (Hill model) acquired during induction serves to calculate the required effect site concentration changes. A TCI device tracks these changes. 10 female patients (12-60 yr) A continuous infusion of remifentanil was started 2 min before induction, thus the necessity to acquire a Hill curve that models propofol effect on the BIS with remifentanil acting as a bias. The results clearly indicate that the closed-loop control of propofol significantly reduces recovery time as compared to the standard anesthesia practice. However, this benefit could only be the result of titrating the drug according to the BIS. Absalom [123] [139] (140] ei al., 2002-2003 Bispectral Index Propofol Outer PID controller setting the reference to a TCI device. There is an additional constraint limiting the maximum change of the infusion rate 10 patients (67 yr ±11 ; 79 kg ± 1 1 ) The authors present a comprehensive control algorithm based on a PID control structure. The way they calculated the gains and the time constant of the controller is not clear. 3 patients out of 10 presented severe oscillations where the BIS value was clearly leading the target plasma concentration, which is a clear sign of the instability of the outer PID controller. The authors also showed in a following publication that the same system can be effectively applied to sedation. Locher ei [141] al,, 2004 Bispectral Index Isoflurane Modified version of Gentilini's cascade controller. The outer controller is a standard PID controller. The inner controller uses the end-tidal expired gaz concentration as feedback. 10 patients (+13 control) (29-59 yrs old) This cascade control system makes use of a feedback measure of the end-tidal concentration in its inner loop. The authors report significant improvements as compared to the manual control of BIS. They also exercised caution when designing the outer PID controller to ensure robustness. However, no mathematical proof of stability was given. The authors reported that the controller reacted promptly during surgical disturbances. They also commented on the fact that the wake up time in the manually controlled group presented much more disparity as compared to the controlled group. Liu et al., [142] 2005 Bispectral Index Propofol Outer PID controller setting the reference to a TCI device. No further information were given as to the controller parameters. It is likely that the TCI pump was age and weight adjusted. 83 patients (closed-loop) ; 81 patients (manual TCI) This study is the largest study involving closed-loop control of propofol using BIS. The authors successfully demonstrated that the closed-loop system was able to maintain the patients within ± 1 0 of the BIS setpoint for 90% of the time, vs. 70% of the time during manual control of the TCI pump. In the closed-loop group, patients received in average 23% less drug and woke up in 30% less time. Note that the controller was also used to induce patients. 8 o o o CO I Table 4.2: Prior Art: a Literature Survey (Part II) CHAPTER 4. CLOSING THE LOOP IN ANESTHESIA: A REVIEW 57 In other words, is the rigorous maintenance of a BIS setpoint brought by a closed-loop system worth the risk of having a machine deciding on the titration profile? Only recent evidence in the literature suggest that there are significant clinical advantages to be gained when closing the loop. L iu et al, in particular, reported a 30% decrease in wake up time, and a 23% decrease in drug usage. Locher et al. have shown that the closed-loop system responds faster to changes in the patient's state and results in a more predictable recovery time. Closed-loop control in anesthesia remains essentially the domain of clinical researchers from European countries. Only the work of Gentilini et al. has been driven by a control engineering perspective. However, their work is mainly based on existing anesthesia sensors and drug effect models. To ensure the success of closed-loop control for anesthesia, it appears essential to re-investigate the anesthesia system from a control engineering oriented point of view. Before designing and testing controllers, it is first important to assess the performance and adequacy of today's anesthesia monitors for use as feedback sensors. For instance, it is interesting to note that none of the cortical monitors currently on the market has been properly identified (i.e., there is no published model of their dynamic response to changes in the patient's cortical state). This stems from the fact that their behavior is highly non-linear 1. The lack of a sensor model is a limiting issue from a control engineering perspective. Also, the review of the propofol pharmacodynamic model presented in the Annex B reveals that there is a large discrepancy in terms of model parameters between the different published models. While this is, in part, due to the large variability that exists between patients, this is also due to a poor model identification methodology. This issue must be revisited. Finally, the great majority of closed-loop attempts did not account for the uncertainty that stems from patient variability, even though some others do mention that their controllers were sufficiently derated to ensure robustness. As a result, some of these attempts resulted in instability characterized by an oscillatory behavior [123]. Analyzing patient variability in order to quantify the uncertainty that exists in anesthesia drug delivery systems is therefore a mandatory aspect of the control design. The goal of this thesis is therefore to assess the feasibility of obtaining a control design that can maintain the adequacy of anesthesia in a large population of patients, while meeting the minimum requirements set forth by anesthesiologists in terms of closed-loop performance during setpoint changes and disturbance rejection. One additional condition to a successful design is that the stability of the controller with regard to patient variability to drug effect must be mathematically proven. lThis issue has recently been brought up by Dr. G. Schneider during the 2005 Advanced Modeling and Control in Anesthesia conference (AMCA 2005), where he has shown that the inherent time delay of most cortical monitors depends on the direction of change of the patient's state (e.g., when patients loose consciousness, the BIS Monitor delay is in average 15 seconds, while the delay is about 30 seconds when patients regain consciousness. PART B: Sensing, Modeling and Control In the following 4 chapters, we review in details every aspect of the anesthesia system from a control engineering point of view. We bring new insights in terms of anesthesia sensors and models, and propose a control system design that is robustly stable and achieves the required clinical performance. In Chapter 5, the inadequacy of existing monitoring technologies for use in close loop systems has prompted us to develop our own solutions, both from the point of view of Hypnosis and Analgesia. Sig-nificant time and efforts were invested into developing new tools and validating them in clinical studies. We give a detailed presentation of the derivation of the W A V C N S (hypnosis). We also show how the same technique can be applied to the analysis of the Heart Rate Variability (HRV) signal for the quantification of analgesia. In Chapter 6, we propose a new approach to pharmacodynamic modeling. Conversely to traditional P D modeling, this new approach accounts directly for the sensor dynamics. We also thoroughly validate the model structure using residual analysis. Forty-four (44) propofol P D models were derived based on induction data obtained during the L M A study (see Annex D). This makes of this study the second largest published propofol P D study. The 44 P K P D patient models established in the previous chapter are then analyzed in order to quantify the expected system uncertainty originating from inter- and intra-patient variability. This analysis is carried out in Chapter 7. We found that uncertainty is cause for concern. We also investigate different method to reduce the uncertainty to a more manageable level. This analysis results in a number of nominal P K P D models and uncertainty weights that depend on the operating mode of the system, and the age of the patient. These models are fully disclosed in the Annex C. In Chapter 8, we design a number of controllers using either a classical P ID loop shaping design, or the more sophisticated Hoo weighted mixed-sensitivity design procedure. We show that controllers achiev-ing relevant clinical performance can be derived. We further show that the stability of these controllers is guaranteed, even under the very conservative assumptions made in Chapter 7. We conclude that pa-tient variability can be effectively dealt with when proper modeling approach and uncertainty analysis are implemented. 58 Chapter 5 Quantifying Cortical and Autonomic Activity Using Wavelets Closed-loop control performances rely directly on the availability and reliability of process output measure-ments. The selection of appropriate feedback sensors to perform these measurements is a critical aspect of any close loop design. In terms of anesthesia, this issue is particularly challenging. As mentioned in Chapter 4, the control of anesthesia can only be achieved through the control of both its hypnotic and analgesic components. Hence, developing appropriate feedback sensing strategies for both endpoints is a priority. Hypnosis and Cortical Activity In terms of hypnosis, most of the efforts pursued in recent years were based on the fact that hypnotic drugs exert their effect at the level of the Central Nervous System (CNS). Most researchers now agree that quantifying cortical activity provides a surrogate measure of hypnosis. One such measure is the BIS monitor (Aspect Medical Systems Inc., M A ) developed in the mid-1990s (see Section 2.2.1.4). This device has become the reference in consciousness monitoring, and, as a result, most research groups working in control of anesthesia are now using the BIS as a feedback sensor. However, its inadequacy for use in a closed-loop framework became apparent in the early stage of our work. The BIS algorithms were designed for stand alone monitoring. As a result, they suffer from many shortcomings that limit their applicability for use within a closed loop control framework. In particular, the BIS inherent non-linearity and time delay, as well as the lack of a closed form transfer function describing its dynamic behavior, imply serious closed loop performance limitations. Analgesia and Autonomic Activity There is currently no technology that can express the patient's analgesic state as a quantified metric. In recent years, however, a renewal of interest in the anesthesia research community for developing such technologies has been observed. One approach, investigated by a research group in Finland [144, 145], is to quantify the magnitude of the noxious stimulation perceived by the patient's body. When considering that surgical stimulation acts 59 CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 60 as an output disturbance, this technology would then provide a measure of such disturbance. However, since it cannot quantify analgesic-based pharmacological depression, this technology cannot be used as an analgesia feedback sensor for controlling analgesic drug titration. Another approach, which we favor, stems from the fact that the Autonomic Nervous System (ANS) is ultimately responsible for managing the response to noxious stimulation. It is also the target of opioids and other analgesic drugs. Furthermore, considering that the ability of anesthesiologists to maintain the autonomic balance dictates the success of the procedure, we believe that the development of a quantified metric of the autonomic activity can provide a surrogate feedback measure of the analgesic state. In 1985, it was proposed that Heart Rate Variability (HRV), i.e., the variability between consecutive R-peaks of the E C G , be used as a source signal to measure the autonomic balance. A number of patents have been issued, but only one intra-operative monitor using that signal (Anemon-I) was commercialized for that particular purpose. No clinical studies validating this technology were published. The Anemon-I has been withdrawn from the market since the early 2000's. However, recent evidences in the literature suggest that the H R V does carry all of the necessary information to derive a meaningful metric [146, 147, 148]. Quantifying both the C N S and A N S may therefore provide the necessary feedback sensors to close the loop and achieve meaningful clinical performance. In terms of the C N S quantification, the E E G has already been identified as a non-invasive signal containing information relative to the patient's hypnotic state. Some success has been achieved by commercial stand alone monitors. However, their utility as feedback sensors is questionable. In terms of A N S quantification, the H R V signal has emerged as a good candidate to quantify autonomic activity. However, it is a difficult signal to analyze and is sensitive to many different confounding factors. In this chapter, we present a signal processing technique based on wavelet analysis. This technology, referred to as Wavelet-based Anesthesia Value ( W A V ) , quantifies the state of a system with respect to its two most extreme states (e.g., awake (performing a mental task), and comatose ( E E G electrical quiescence) for quantifying cortical activity). As the signal evolves from one state to the next, the changes are quantified and expressed into a bounded scale. In this Chapter, we show that this technology can be used to analyze and characterize both the C N S and A N S activity. We therefore propose two new anesthesia indices: the W A V C N S , which quantifies cortical activity based on E E G analysis, and the W A V A N S , which quantifies autonomic activity based on H R V analysis. In Section 5.1, we present the wavelet transform, and, in particular, its ability for detecting rapid temporal changes in non-stationary signals. We present in details the W A V technology in Section 5.2 and its application to the quantification of cortical activity based on a single frontal E E G signal. The resulting index, referred to as W A V C N S , was implemented in real-time and coupled to an E E G signal data acquisition system for use in the operating room. The BIS monitor being a clinical reference for consciousness monitoring, this setup allowed us to carry out a clinical study aimed at comparing both CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 61 Figure 5.1: Dyadic frequency tiling with L = 3. sensors. The results of this study are summarized in this section. We then show in Section 3 how the W A V technology can be applied to the analysis of the H R V signal. For illustration purposes, we present results from 3 clinical cases, where we compared the W A V A N S to the W A V C N S and the anesthesiologist's assessment of the patient's analgesic state. We also compared the W A V A N S to the heart rate, blood pressure and drug titration recorded during the cases. These preliminary results confirmed the strong potential of the W A V A N S for assessing the autonomic balance. Finally, Linear Time Invariant (LTI) transfer functions modeling the dynamic behavior of both W A V sensors are presented in Section 4. 5.1 The Discrete Wavelet Transform This section is mainly an overview of wavelets and Wavelet Transform (WT) . It is an opportunity to introduce the notations that will be used throughout this Chapter. We also present the rationale behind the use of W T to analyze the E E G in order to estimate the hypnotic depth. We invite interested readers to refer in particular to [149] and [150] for a thorough review about W T and its applications. 5.1.1 Standard Wavelet Dyadic Decomposition A dyadic wavelet transform w is a transform that decomposes a given signal x into a number of coefficient sets. Each of these sets represents the activity of the signal in a particular frequency band: cl band: [0, /s /2 L ] d£ band: [fs/2L, fs/2L~1} w : x dj" band: [/s/2'" 1 , fs/21} (5.1) ^ df band: [/s/4,/s/2] For instance, a L-level dyadic wavelet decomposition yields L + 1 sets of coefficients corresponding to a dyadic tiling of the frequency domain, see Fig. 5.1. Furthermore, each set {ct} and {d;}i=1 . jr, is a time series of coefficients describing the time evolution of the signal in the corresponding frequency band. Hence, two signals having the same frequency content will yield different sets of coefficients if the signal patterns are different (i.e., if the frequency components have different phases). CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 62 Similarly to Fourier analysis, which was made practical by the discovery of the Fast Fourier Transform (FFT) algorithm, wavelet analysis remained the apanage of mathematicians untill the mid-1990s when the use of filter banks appeared to yield much faster and efficient ways of performing the necessary calculations. Before then, wavelet analysis existed only in its continuous form, involving a very large computational complexity ill-suited for real-time applications. In a Discrete Wavelet Transform (DWT) framework, the wavelet decomposition is obtained through the use of a filter bank and 2 F I R filters, HQ (low pass) and Hi (high pass). The characteristic of the coefficient sets {CL} and {di}i=\t...ti, therefore depends on the wavelet filters Ho and Hi associated to the transform w. 5.1.2 Change in Phase and Wavelets The coefficients obtained through the Wavelet Transform (WT) capture both the frequency and time information of a given signal. Wavelet analysis is therefore particularly well suited to track both time and frequency changes in a signal. The analysis of many biological signals have benefited from wavelet analysis. To illustrate this, let us consider a signal composed of 3 frequency components. The phase of one of the component is allowed to drift slowly in time from 0 to TT following the time course shown in Figure 5.2.d. Even though the change in the signal pattern is obvious (Figure 5.2.a), the power spectra (Figure 5.2.c) was unsuccessful in tracking the changes in the signal latency, while the wavelet decomposition (Figure 5.2.b) clearly shows an evolution of the signal over time. It is interesting here to make a parallel with the arguments presented by Bowles et al. [52] in favor of bispectral analysis (see Figure 2.4). Following Rampil's observation [51] that anesthetic drugs tend to synchronize the generation of postsynaptic potentials and thus affect the signal latency, it appears that the time-frequency localization property of the W T makes it well suited for capturing the evolution of the E E G with increasing anesthetic depth. This is illustrated in Figure 5.2. 5.2 Estimating the Anesthetic Drug Effect: the W A V C N S Note that contributions in the development of the W A V C N S are shared with Ms. T. Zikov [151]. 5.2.1 Concepts and Derivation The brain is the target organ of anesthetic drugs.. The electroencephalogram (EEG) signal being a non-invasive measurement of cortical activity, the effects of anesthetic drugs onto this signal has been thoroughly studied. Because of differences in the mechanism of action, the effect of different anesthetics usually results in different E E G patterns. However, for most drugs and with increasing blood concentration, the E E G signal evolves from a low-amplitude, large-bandwidth, noise-like signal, to a high-amplitude slow-waves signal. If large amounts of anesthetic drugs are used, the E E G activity eventually disappears, resulting in CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 63 Wavelet coefficients Frequency (Hz) Phase (rad) (b) (c) (d) Figure 5.2: Changes in latencies in a signal can also be tracked by the wavelet coefficients (from [151]). (a) Time series signal. Note how the signal pattern changes significantly when changing the phase of one of the frequency component, (b) 'Periodogram' of the wavelet coefficients (note how the change in phase is clearly noticeable in the wavelet domain), (c) Classical Fourier periodogram. (d) Time profile of the phase change. CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 64 0 X Figure 5.3: The function h quantifies the state of a system by attributing a unique value to each of the operating modes of the system, ff the system evolves in a monotonous fashion from state a to state b, it is required that h is also monotonous. an isoelectric signal (i.e., complete absence of cortical activity). Since this trend can be observed for most anesthetic agents, it is commonly assumed that the E E G can be used to estimate the hypnotic state of the patient, and thus provide a measurement of the drug effect. In order to avoid complex and time-consuming interpretations of raw E E G signals, a common approach is to extract a single univariate value that represents the patient's state. Since there is no gold standard for measuring depth of anesthesia, the various analysis techniques used to calculate this value are usually derived such that the resulting quantitative value is sufficiently correlated to the anesthesiologist's qualitative assessment of the patient's state. In order to illustrate our methodology for deriving an index of hypnosis, let us consider the patient to be a system which can evolve between two states a and b. Our goal is to derive an index that characterizes the operating mode of this system, i.e., finding a function h, which, given an observation x acquired while the system operates between a and b, yields a single value i comprised between two finite bounds, e.g., 1 and 0 (see Fig . 5.3): / i : x € [ a , 6 ] — > i e [ 0 , 1 ] . (5.2) In terms of hypnotic depth, the observation x can be a short E E G epoch, long enough to contain the necessary information, and short enough to allow frequent updates of i. The states a and b correspond to the two extremes of the signal observational window, i.e., the awake and isoelectric E E G . As far as anesthetic drugs are concerned, this window is sufficient since it is usually not desired to titrate patients beyond the point where cortical activity is totally suppressed. In this section, we will first give some very general guidelines on how to define the function h. We will then apply these guidelines to the estimation of the anesthetic depth, using data obtained from healthy subjects and anesthetized patients. 5.2.1.1 Main Concept and Approach In order to define the function h, we first need to assume that both states a and b correspond to well established operating modes from which observation data are available. In the case of E E G analysis, this CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 65 Figure 5.4: The characteristic of the function h depends on the selection of the feature function / . (a) Both hi and h2 are valid functions, whereas h$ and hA cannot be used to represent adequately intermediate states, (b) The linearity of a given function h can be assessed by introducing new data sets, (c) Variability is another important aspect. Por instance, even though hi is nearly linear, it fails in characterizing properly the intermediate states Ti and T-2- The function hi is preferred here since it allows for the proper discrimination between consecutive states. assumption is verified since a corresponds to the awake state (e.g., when able to perform a mental task) and b corresponds to an isoelectric signal. Since these states are well defined, we can then obtain two reference data sets corresponding to observations of the system in the states a and b: (5.3) Ra — (xa,fc, k = 1,2,... M} (state a), Rb = {xb,k, fc=l,2, . . . M } (state b), where each observation vector xQjfe and x ^ contains a finite number of samples and represent the kth epoch of the data sets Ra and Rf,. To characterize the data sets, a chosen feature f is calculated from each epoch. We denote the feature function / as: / : x —> /(x) = f (5.4) Each epoch of the reference data sets is then associated with a feature or f^j.. This feature can be either a scalar (such as average value, Root Mean Square (RMS) amplitude, maximum or minimum value, standard deviation, etc.) or a vector (such as histogram, probability density function, etc.) derived in the original signal domain or in a transformed signal domain. A particular state can then be characterized by averaging the feature sets {ia,k\k=\,...,M and {h,k}k=i,...,M- This results in two references which are the averaged features fa and ff,: > = • Eifcli fo,fe, and (5.5) < b = ~M ' ^ f c = i 6 ' f c ' Let us now assume that the system operates in the intermediate state c and that an observation x c is provided. Using the corresponding feature fc = /(x c), it is then possible to estimate how far the system has evolved from the state a towards the state b by comparing the value fc to the references fa and ff,. Two CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 66 indices ja and jp are defined such that: ja = < life - t i l l 3P = V life - f t l l l and (5.6) where the norm || • | | i for a vector v is defined as: N (5.7) The norm || • | | i quantifies the difference between fc and fa (or fj,) by integrating the distance between these two vectors. The indexes ja and jp thus measure a distance of the system from either state a or b. Note that higher degree norms can be used for this analysis, however, they would emphasize larger differences which might be due to outliers, hence leading to noisier indexes. In order to maximize the signal-to-noise ratio, the two indexes ja and jp can be combined to yield a single descriptor j = jp — ja. It is convenient here to define a function g that, given an observation x, associates the resulting index j: g : x — g(x) = j = ||/(x) - ffcHi - ||/(x) - fa\\i. (5.8) In order to scale the final index i such that i —> 1 (state a) and i —* 0 (state b), the function g is applied to the entire reference data sets. This yields two sets j a and j ; , : 9-Ra—>3a = {ja,k, k = l,2,...M} g:Rb—>jb = {jbik, k = l,2,...M} The function h is a scaled version of g, and is thus defined as: i Ti (5.9) h : x — > h(x) = i = g(x) Ja Jb Ja Jb (5.10) where: Ja and (5.11) T — X^M n Jb — ' l^k=\3b,k-5.2.1.2 Feature F u n c t i o n Selec t ion The performance of the index i as a substitute to the qualitative assessment of the patient's state clearly depends on the characteristic of the function h over the state continuum [a, b]. This is illustrated in Fig. 5.4.a where hi and h2 are adequate as they uniquely define each intermediate state. Conversely, a function such as / i 3 should be avoided as there exist different states which yield similar index values, thus leading to confusion in the interpretation of the index. As a rule of thumb, any function whose characteristic exhibits CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 67 a change of sign in its first derivative should be avoided. Similarly, a function such as h4 is not very useful since it does not discriminate properly between different states. The selection of / is critical in shaping the characteristic of h. For instance, a good discrimination between consecutive intermediate states is achieved if a feature function / is determined such that the evolution between the states a and b is linear. We define / (resp. h) such a function. While linearity is an important aspect, it is not sufficient to guarantee proper discrimination between states. Another important aspect is the variability of the index i while the system operates in a stationary mode. A large variability may cause a significant overlap between consecutive states, indicating that either the observation signal or the feature function are poorly suited for this analysis. In order to facilitate the selection of the feature function / , we then introduce two design parameters: L (linearity parameter) and V (variability parameter). Both parameters rely on the use of new data sets T i , T2, ..., T/v corresponding to intermediate states of operation {e.g., sedation, general anesthesia, deep hypnotic state, etc.). For each data set, R or T, we can associate a desired index value i: fi: {Ra, Tu..., TN, R b } —- {1, i\,..., iN, 0}. (5.12) Note that, by definition of h, ia = 1, ia = 1, if, = 0, and it = 0. L i n e a r i t y Linearity can be easily assessed by calculating the corresponding index for each epoch of each data set. We thus obtain a set of vectors { i n } of length equal to the number of epochs in the corresponding data set: h : {Ra,Ti,.. - ,TN, Rb} — • {ia , i i , • • -, i/v, ifc}- (5.13) We can assess the characteristic of h by simply comparing the average vector values in = mean(i„) with the set of desired values in in (5.12), see Fig. 5.4.b. The parameter L is the normalized root mean square error between the characteristics of h and h: N s l A ' • ( 5'1 4 ) L is bounded between 1 and 0. Note that to normalize L, we divided the root mean square error of (h — h) by the maximum possible root mean square error {e.g., the function /14 in Fig . 5.4.a. would typically have L = 1). To optimize the linearity of h, it is therefore necessary to minimize L. V a r i a b i l i t y The variability of h can be assessed by calculating for each vector i n of (5.13) the standard deviation cr(in). If a statistically significant discrimination of each intermediate state T i , T2, . . . , T/v is a critical feature of h, it is necessary to ensure that: tf(in) + °"(in+l) < l*n ~ (5.15) i.e., that the standard deviation intervals of two consecutive states does not overlap, see Fig . 5.4.c. It is, however, important to keep in mind that the standard deviation a{i) is strongly influenced by the eventual CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 68 non-stationarity of the data set T. This is mostly true when T is composed by data acquired from different systems or under different recording conditions. This is also true when the data sets are collected based on the evaluation of a human operator. We define the design parameter V as: where a(io) = cr(ia) and <r(ijv+i) = cr(h)- For the same reasons outlined above, it is advised to keep V < 1, or as small as possible if the stationarity of each set T cannot be properly guaranteed. 5.2.1.3 Data Sets In the search for a function h that would adequately estimate the anesthetic depth, 5 different observation sets corresponding to 5 distinct hypnotic states were recorded, see Fig. 5.5 for representative samples: - Ra (awake state): 15 minutes recorded from 5 healthy adult subjects (3 minutes each). Subjects were asked to keep their eyes shut and minimize muscle activity while concentrating on a mental task. - T i (light R E M sleep): 15 minutes recorded from 3 subjects (5 minutes each). Only epochs without ocular artifacts were included in this set. - T2 (general anesthesia): 18 minutes recorded from 6 patients (3 minutes each) undergoing minimally invasive arthroscopy surgery. - T3 (deep hypnotic state): 9 minutes recorded from 5 patients. Delta activity is usually transitory, which explains the limited amount of data available. - Rb (isoelectric E E G ) : f»5 minutes recorded from 2 patients exhibiting electrical quiescence. The data sets T2, T 3 and Rb were collected from surgical audit cases. While a variety of anesthetic regimens were used during these cases, only drugs provoking a concentration dependent depression of the cortical activity were used (vapour anesthetics (isofiurane, desflurane, sevoflurane, NO2), propofol and midazolam). While a simple observation of the E E G was sufficient to extract periods of isoelectric activity (Rb), the classification of E E G epochs between T2 and T3 was done following the anesthesiologist's assessment for each patient. The anesthesiologist was first asked to identify periods of adequate anesthetic depth, i.e., when patients were unconscious in the context of surgical anesthesia. The E E G corresponding to these periods were then reviewed post hoc to identify periods of deeper hypnosis (apparent delta waves) which were then classified as the T 3 set. The remaining E E G data were then included in the T2 set. Each E E G signal was recorded from a frontal differential channel ( F p i - A T i (the left outer malar bone) and F p z as ground), using the Crystal Monitor Model 1 6 ™ (Cleveland Medical Devices Inc., OH), a (5.16) CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 69 Figure 5.5: Normalized E E G signals at different anesthetic depths. Each 1-second epoch is first detrended and then normalized using the RMS amplitude, (a) Awake healthy subject (b) Light R E M sleep (the spikes due to R E M activity were manually removed) (c) General anesthesia (d) Deep anesthetic state (5-waves) (e) Isoelectric state (normalized by the RMS amplitude of the last non-isoelectric epoch). sampling frequency of 480 S/s and a resolution of 16 bits. The signals were then resampled at 256 S/s and notch-filtered at 50/60 Hz prior to analysis. The E E G signals were further divided into 1-second epochs (256 samples per epoch). Because the signals were measured through contact electrodes with a high input impedance amplifier, there are a number of factors (e.g., skin conductivity and thickness, electrode impedance and placement) which can affect the amplitude of the recorded signal. To minimize the influence of these factors, each epoch is first detrended and then normalized by its R M S value prior to analysis: (5.17) R M S ( x ) ' This also compensates for the bias and calibration errors that are expected between different acquisition devices. Note that for the isoelectric epochs, the normalization was performed using the R M S value of the last E E G epoch preceding the occurrence of isoelectricity. Finally, the desired values in terms of the mean index values were defined for each data set to reflect the anesthesiologist's assessment: ^ W A V C N S : {Ra,Ti,T2,T3, Rb} — • {1,0.75,0.5,0.25,0} (5.18) 5.2.1.4 Feature Function / W A V In order to extract information from the E E G signal pertinent to depth of anesthesia, the Probability Density Function (PDF) of the wavelet coefficient set df was chosen as the feature function, and we denote CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 70 it as / W A V C N S : / W A V C N S = * — • / W A V c n s ( X ) = PDF(dj"). (5.19) This choice is motivated by a number of facts. First, the information regarding the anesthetic depth of a patient lies in both the spectral distribution of the E E G and the pattern of the signal itself. The wavelet transform is therefore an appropriate tool to extract this time-frequency information. However, the resulting coefficient sets cannot be used right away. It is necessary to format the information they contain using statistical tools. The interest in using the P D F in this case lies in the fact that the area of a P D F curve is, by definition, always equal to 1. Therefore, the indexes ja and jp are naturally bounded between 0 and 2. Secondly, even though we have already chosen the particular structure of / , this choice offers design flexibility through the selection of the frequency band I and wavelet transform w. Thus, by iterating between different wavelet filters and frequency bands, a search algorithm can easily be written to determine the best / and w. 5.2.1.5 Se lec t ion o f the Frequency B a n d and Wavele t F i l t e r Let us denote with W A V C N S the index i defined in (5.10), and corresponding to the feature function / W A V C N S • To select the appropriate frequency band I and wavelet transform w, we assess the linearity and variability of the W A V C N S for different frequency bands and wavelet filters using the design parameters L and V defined in Section II (5.14) and (5.16). The linearity and variability of the W A V C N S were assessed for the three frequency bands df, d™ and dg, and for the first 16 Daubechies wavelet filters, see Fig. 5.6. In addition to the standard definition of L given in the previous section, we added a sign information to this otherwise positive value. A change in the sign of the linearity parameter displayed in Fig . 5.6 indicates a change in the dominant curvature of h (a positive value indicates concavity, while a negative value indicates convexity). This observation is of interest when considering that, depending on the concavity or convexity of the function, the resulting index will either be more sensitive to changes in the lighter hypnotic states and less sensitive to changes in the deeper states (convexity), or conversely, more sensitive to changes in the deeper hypnotic states and less sensitive to changes in the lighter states (concavity). The selection of the frequency band and filter allows for the modification of this characteristic. The linearity analysis indicates that, either the 7-band (32-64 Hz) in conjunction with any wavelet filter of order >4, or the E M G band (64-128 Hz) and the Daubechies filter #2, yield a nearly linear characteristic. These choices also agree with the variability criterion according to which the parameter V should remain well below 1. Note that these results are in agreement with the recent findings suggesting a fundamental relationship between consciousness and 7 activity. It has also been reported by John et al. [152] that loss of consciousness after anesthetic administration showed a marked drop in the 7-band activity, along with an increase in slower CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 71 ds band (16-32 \ di band -®—-© 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 DAUBECHIES W A V E L E T F I L T E R # Figure 5.6: Selection of the frequency band and wavelet filter for the sampling frequency of 256 S/s. (a) Linearity parameter L. A positive value indicates a concave characteristic, while a negative value indicates convexity. By appropriately choosing the wavelet filter, the characteristic of h can change from concave to convex, (b) Variability parameter. waveforms, which is in line with our findings. It is also interesting to note that the E M G band yields a similar result for the particular wavelet filter Db2. However, this is explained by the fact that this low order filter still preserves a significant amount of 7-band information. Higher order filters, which extract mostly E M G information, lead to a convex characteristic of the index which does not give a good discrimination between deeper states. Given these choices, it is advisable to use the 7-band instead of the E M G band. First, the use of the 7-band decreases computational complexity by allowing to lower the sampling rate to 128 S/s, and thus, use the d™ coefficient set instead of d™ via a single convolution operation. Also, the signal corresponding to the E M G band is of lower power than that of the 7-band, which makes it more susceptible to environmental fu: isoelectric state -0.5 0 0.5 W A V E L E T CO E F F I C I E N T S Figure 5.7: Awake and isoelectric reference PDFs based on / ^ A W A V C N S ' CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 72 0.8 0.6 & 0.4 U.2 'a i v V > S s's/l d: band iS. fom</ \ \ ^ \ (16-32H?) \ _ \ ; . ! di band ' sC*Ssv^- \ (64-128 H%) i— 1 1 1 1 1 7", r 2 T3 Rj Unfiltered index Filtered index (Average filter of 30 sec) Figure 5.8: Characteristic of the optimal function ^ W A V C N S ^ o r t ' l e sampling frequency of 256 S/s. (a) Linearity characteristic, (b) W A V C N S values for the data sets Ra, T±, T2, T3 and Rb. Note the larger variability in the intermediate data sets. noise. Furthermore, the use of the E M G band would also require to notch-filter both at 50/60 Hz and 100/120 Hz to reduce the electrical noise from the mains. Finally, the E M G band is more susceptible to neuromuscular blocking agents, which aim at suppressing muscle activity. A n index based on this band would therefore be particularly sensitive to this type of drugs. For the 7-band, an optimal result is obtained when using the Daubechies wavelet filter #6 (w[Dbe]). The optimal feature function / W A V C N S * s t n u s defined as: / w l v C N S = * — / w l v C N S ( x ) = PDF(dr[Dbel), (5.20) where the signal x is sampled at 128 S/s. Figure 5.7 illustrates the reference awake and isoelectric P D F s based on / ^ A V C N S a v e r a S e ( i o v e r the sets Ra and R/,. For more details on the derivation of the PDFs , see Section 5.2.2.2. The corresponding characteristic of ^ V V A V C N S ^ s P l ° t t e d m Fig. 5.8.a. Results show a nearly linear characteristic where: ^ W A V C N S : iR°> Tl> T2> T3> R b ) ~ * i 1 0 0 - 7 7> 5 4> 2 0> ° ) (5-21) Note that the values in (5.21) are average values obtained from the training data sets. For each data set, the corresponding W A V C N S values are plotted in Fig. 5.8.b. The variability in the awake and isoelectric states is small as compared to the variability observed from the intermediate states. This is due to the fact that the awake and isoelectric states are, by definition, very well identified. Conversely, the intermediate states are subject to interpretation. This results in a higher non-stationarity of the T data sets. CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 73 5.2.2 Pract ical Issues and Implementation This section addresses practical issues concerning the implementation of the W A V C N S algorithm. 5.2.2.1 Pre-Processing Notch Filtering The E E G signal is a low power signal particularly susceptible to environmental noise. For instance, electrical heaters used as warming units in the operating rooms draw large amount of current from the main electrical outlets. When these units are in proximity of the E E G acquisition device, strong 50/60 Hz and 100/120 Hz perturbations may be superimposed to the E E G signal. While appropriate shielding can alleviate this problem, it is necessary to use a 50/60 Hz notch filter to remove this noise. Resampling The E E G signal needs to be resampled to a rate suitable for the analysis. Since we have determined that the 7-band (32-64 Hz) is of interest, a sampling rate of 128 S/s is optimal. Note again that an advantage of using this sampling rate is that only one convolution of the filter w[Dbe] with the signal is needed to obtain the d™' D b 6' coefficients. When higher sampling rates are used, the wavelet coefficients in the band of interest are obtained through a series of filter banks, which adds to the computational complexity. Artifact Detection and De-noising Ocular Artifacts (OAs) and head movements can make the W A V C N S particularly unstable in lighter hypnotic levels (awake and light sedation), when the patient is still conscious. This difficulty arises mostly before and during induction. One strategy to deal with these artifacts is to systematically reject corrupted epochs from the analysis and simply set the index to its previous value. However, this can result in substantial data loss. A better approach is to remove the artifact from the raw E E G signal, without perturbing the high frequency content of the E E G . This can be easily done for two reasons. First, OAs and head movements typically perturb the E E G from 0 Hz up to 16 Hz. Thus, the high frequency information is not affected. However, they substantially affect the R M S value of the signal, and thus the entire analysis. Second, this type of artifact signal, which is superimposed with the true cortical E E G signal, is of high amplitude com-pared to the relatively small E E G amplitude. Using a redundant wavelet decomposition, which provides a good temporal resolution, this signal induces large coefficients that are easily distinguishable from the smaller background coefficients of the true E E G . A simple technique such as thresholding of wavelet coef-ficients in the lower frequency bands (up to 16 Hz) is sufficient to remove the perturbing artifacts from the source E E G , see Fig . 5.9.a [59]. This technique is particularly effective in removing most of the artifacts while keeping the exact amplitude and phase of the high frequency components, see Fig. 5.9.b. A difficulty arises during the state transition from awake to anesthetized. This transition is fairly rapid. Further, the initial induction bolus usually drives the patient into a deep hypnotic state. As a result, 5-waves are likely to occur within a few seconds from the loss of consciousness. The 5-waves being particularly CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 74 Figure 5.9: Effect of ocular artifact de-noising during an awake period, (a) Raw E E G and de-noised EEG. (b) Note how the high frequency information remains unaltered by the de-noising technique, both in terms of ampfitude and phase. similar to OAs (i.e., low frequency, high amplitude), it is necessary to switch off the de-noising technique as soon as the patient does not exhibit ocular activity. The moment of anesthesia induction, which corresponds to the bolused administration of an intravenous anesthetic, is suitable for switching off the de-noising technique as the patient rapidly looses consciousness. From the perspective of automation in clinical anesthesia, this solution is viable since the computer is allowed to directly command the infusion device. 5.2.2.2 The W A V C N S Algorithm Normalization For the reasons outlined in Section 5.2.1.3, the E E G epochs are first detrended and normalized prior to analysis. This normalization is carried out on E E G epochs with R M S amplitudes above 4 M V R M S -During electrical quiescence (i.e., when the R M S amplitude of the E E G is less than 4 UVRMS), the nor-malization using R M S values close to zero would lead to the amplification of measurement noise. Therefore, we directly assign the feature f*b (see Fig. 5.7) to any isoelectric E E G epoch: • ^ W A V C N S ' ^ ^ ^ * s o e ' e c * " c * ft>- (5.22) The Redundant Wavelet Transform For the derivation and implementation of the W A V C N S algorithm, we use a redundant D W T , so-called Stationary Wavelet Transform (SWT), which provides the same time resolution in all frequency bands of decomposition. There is a number of reasons for using a redundant transform. First, when choosing an optimal frequency band for the analysis, the S W T provides the same number of coefficients in each frequency band (i.e., it is equal to the sampling rate for a 1-second long epoch), and therefore enables easier comparison between them. (For a non-redundant transform, the number of coefficients decreases by a factor 2 for each level of decomposition.) In addition, a greater number of coefficients results in better temporal resolution in the decomposition bands, and, therefore, smother reference PDFs . Secondly, at the pre-processing stage, the better time resolution of the S W T provides an improved CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 75 characterization of artifacts, as well as a smoother estimation of the signal of interest (i.e., the true cortical E E G ) after thresholding in the wavelet domain [153, 154]. However, the utilization of the S W T somewhat increases the computational demand, since the complex-ity of the S W T is 0(N- log2N) for a signal of TV samples, in comparison to O(N) for the D W T . Finally, note that caution has to be taken when performing the wavelet decomposition on epochs of short length, since boundary effects can significantly perturb the result. Therefore, prior to the decomposition, each epoch is extended using a standard procedure such as symmetrization. After applying the S W T , only the coefficients that are not affected by boundary effects are kept for further analysis. Probability Density Function When calculating the P D F of the wavelet coefficients, special care has to be taken when selecting an appropriate number of bins. It has been shown [155] that the optimal histogram bin size, which provides the most efficient, unbiased estimation of the probability density function, is achieved when: Bw = 3.49 • a - . A T 1 / 3 , (5.23) where Bw is the width of the histogram bin, a is the standard deviation of the distribution and N is the number of available samples. However, when it comes to estimating the depth of consciousness, the statistics of the E E G signal can significantly vary between the different states. For example, the variance of the wavelet coefficients in the 7-band for the awake state is much greater than for the anesthetized state. Furthermore, for the isoelectric E E G , the variance is almost equal to zero (see Fig. 5.7). Hence, an optimal number of bins defined by (5.23) for the proper representation of the awake state might actually be too small to adequately represent deeper states, and might cause them to be mistaken for the isoelectric state. Similarly, if the number of bins is too large, the P D F waveform might exhibit a 'comb' effect in the awake and light sedation states, while precisely defining the deeper states. Therefore, a compromise needs to be reached when choosing the number of bins for the analysis so that all states are represented in a satisfactory manner. Based on the length of an epoch, a number of bins equal to the sampling frequency of the signal appears to provide an adequate compromise. In fact, the isoelectric reference P D F becomes zero everywhere except at the origin, where it is equal to one. 5.2.2.3 Post-Processing Trending Using an analysis epoch as short as 1 second can lead to a particularly noisy W A V C N S index. To overcome this problem, a filtering stage has to be added to the W A V C N S algorithm to extract the most significant information. There are a number of techniques which can be used. The simplest one is to use an averaging window, where the displayed index becomes the mean value of the previous index values over a fixed horizon. For example, the BIS monitor uses a window of 30 or 15 seconds, since it yields good results in terms of CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 76 o 2nd order low -20 0 180 0.001 0.01 0.1 F R E Q U E N C Y [ H Z ] Figure 5.10: Characteristic of a 30-second averaging filter and a 2 n d order low pass filter (uio = 0.02 Hz and C~l)-noise reduction. During the acquisition of our clinical data, the BIS monitor used a 30-second averaging period. Therefore, we have chosen the same averaging period for our analysis, in order to perform the direct comparison between the WAVCNS and the BIS. A property of this filter is that the filtered index only depends on the past 30 seconds of data, which makes its interpretation straightforward. However, while F IR filters are simple to implement, they introduce a non-minimum phase element which is detrimental in applications involving control and identification. In that respect, the use of an IIR filter is more suitable. For instance, a simple second order low pass filter characterized by its parameters UJ$ and C, such as described in Fig. 5.10, gives a better high frequency noise rejection profile, while being minimum phase. We found that an IIR filter with a cutoff frequency of 0.02 Hz and damping factor of 1 is an adequate alternative to the 30-second averaging window filter. Sca l ing To allow for a better comparison with the BIS Monitor (see Section VI), the W A V C N S was scaled between 100 (awake state) and 0 (isoelectric state) by multiplying the right hand side of (5.10) by 100. 5.2.3 Clinical Results To validate the WAVCNS, a clinical study was conducted in 2002 at the University of British Columbia Hospital ( U B C H ) . The aim of the study was to compare the real-time W A V C N S to the BIS v.3.4 (Aspect Medical Systems, M A ) within a clinical setting. Note that, to allow for a better comparison, the WAVCNS IIR trending filter was replaced by a 30 seconds averaging filter similar to that of BIS. The study protocol and demographics are summarized in the Annex D . l . We present in Fig. 5.11 two cases representative of the time course of both indices. Probably one of the most surprising results from this study is the striking overall similarity between the BIS and the WAVCNS • CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 77 Figure 5.11: Time course of the B I S and W A V C N S for Patient #22 and Patient #8. Yet, after closer inspection, differences during periods of large transients become apparent. This observation led us to separate steady state behavior from transitory behavior in the comparison of the two indices. Our results and conclusions are the subject of the two following subsections. We also observed that periods of burst suppression patterns are yielding different time courses. We felt that this result could carry strong clinical interest. We therefore present a short discussion concerning the behavior of both BIS and WAVCNS at the end of this section. 5.2.3.1 Steady State Behavior The BIS index is derived and tuned specifically to reflect the anesthesiologist assessment of the patients' anesthetic depth. The BIS value is nowadays a reference against which other indexes are compared [156, 157, 158]. However, this statement holds true only in steady state, i.e., when the patient's state of consciousness remains constant, since transitory states are affected by the time delay induced by the bispectral analysis [159]. In order to evaluate the accuracy of the WAVCNS, epochs of steady state behavior were extracted from each case using a decision algorithm based on the BIS value. A 1-minute sliding window was used to assess the stationarity of the BIS. If the BIS value remained constant within this window (with bounds set at ± 5 of the average BIS value) the epoch was classified as steady-state. Consecutive steady-state epochs were pooled together. The first (or last) 30 seconds corresponding to the beginning (or the end) of each steady-state period were removed under the assumption that they are transitory. A total of 13 hours of recorded data was identified as stationary (73% of the total recording time), representing more than 45,000 E E G epochs. The WAVCNS VS. BIS correlation density plot is presented in Fig . 5.12. There is a strong correlation CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 78 100, 80 L 60 L 40 20 L 2li t = 0.969, p = 0 95% conf. int. = [0.9685 0.9696] J I I 40 60 BIS v.3.4 00 80 100 100 HO 60 40 CO 2 20 > i 0 • -20 3 -40 -60 -80 •100 Total steady state time - 13 his 16 min = 72.5°, 'o of total time m e a n +2-SD = 1 2 . 5 / \ \ *• m e a n = 1. 4 m K i. m e a n - 2 - S D = - 9 . 8 1 I I I i 1 1 1 0 10 20 30 40 50 60 70 80 90 100 (BIS+WAVCNS)/2 (b) Figure 5.12: Correlation during steady-state operation between the BIS (v.3.4) and W A V C N S - A total of 13 hours of steady-state data were collected from the arthroscopy study, (a) Correlation density plot (b) Bland-Altman test. between the two values in the 30-60 range. The WAVCNS was slightly higher in the light anesthesia/awake range (60-100). However, most of the data points from this region were obtained from the emergence phase where the BIS value is known to be underestimated [79]. 5.2.3.2 T rans i to ry Behav io r The WAVCNS algorithm was designed to exhibit a smooth transitory behavior. In order to evaluate the performance of the WAVCNS during large transients, we observed the time courses of both indexes during induction and emergence. I n d u c t i o n For this analysis, the study population was separated into two groups depending on the re-action to airway manipulation and L M A insertion. Note that airway management and the insertion of the L M A typically occurred 40 to 80 seconds after L O C . Patients exhibiting no reaction or mild extremity movements were classified in Group 1 (n = 10). Patients exhibiting airway reaction (biting, coughing, swallowing, etc.) or purposeful motor movements were classified in Group 2 (n =9). One case (patient #15) was excluded due to the use of a neuro-muscular blocking agent during induction (the presence or absence of reaction could not be assessed). In Fig . 5.13, the time courses of the BIS and WAVCNS are plotted for both groups. Each individual time course is synchronized using the L O C event. In Group I, the induction is fast and profound. The WAVCNS typically starts decreasing with the L O C event. As compared to the WAVCNS, the BIS value exhibits a more erratic behavior during this transitory phase. Conversely to Group 1, the induction titration in Group 2 was not adequate to warrant the absence of CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 79 -Loss Of Consciousness (LOC) 100 50 I 0 100 Loss Of Consciousness (LOO 3 >. 50 t » - 1 \ " 5 1 1 -20 20 40 60 T I M E [S] (a) 80 100 120 Figure 5.13: Time courses of the B I S (v.3.4) and W A V C N S during induction for both patient groups (each case is synchronized at the L O C ) . (a) Population showing no reaction to L M A insertion (n=10). (b) Population showing some reaction to the L M A insertion (n=9). reaction to the insertion of the laryngeal mask. The post-insertion BIS and W A V C N S time courses exhibit more variability than those of Group 1. It is interesting to notice that there is no significant difference in the BIS time course prior to L M A insertion between the two groups. Conversely, the rate of descent and nadir value observed in the W A V C N S shows significant difference between the groups. Thus, it can be hypothesized that the W A V C N S time course during the awake-to-anesthetized transition might be capable of predicting the reaction to L M A insertion. This finding was further corroborated by Lundqvist et al. [160] in a study involving 50 patients. In F ig . 5.14.a, the average time courses of the BIS and W A V C N S are plotted together. We observe a marked delay between the L O C event and the BIS values indicating loss of consciousness (i.e., < 80). The W A V C N S precedes the BIS by an average of 15 seconds in both groups. Note again how the W A V C N S drop is well synchronized with the L O C event. Emergence Following the clinical protocol, 3 events were recorded during the anesthesia emergence: - Pnrt (Patient Not Reacting): the patient was not reacting to the anesthesiologist's touch and/or voice, - P r t (Patient Reacting): the patient was reacting (movement, vocalization) to the anesthesiologist's touch and/or voice, or to the surgical environment, - P r s (Patient Responding): the patient responded to a verbal command. The P r s event only indicates a time when the patient was able to follow verbal command. However, it is not the earliest time when the patient was able to process cognitive thoughts. This is due to the fact that the clinical protocol of this observational study did not make provisions for how frequently the anesthesiologist was to assess the patient's cognitive state during emergence. Therefore, the P r s event cannot be used as a marker of Return Of Consciousness (ROC). CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 80 Figure 5.14: Time delay between the BIS and W A V C N S - (a) Induction (both study groups), (b) Emergence. Conversely, the P r t event indicates the earliest time when a patient reaction was observed. This reaction was triggered either by the anesthesiologist's probing, or by environmental factors. Since the P r t corresponds to a well defined point in time, which is easily observable, it can be used to synchronize the WAVCNS and BIS time courses during emergence, see Fig. 5.14.b. Similarly to induction, we observed a delay in the BIS as compared to the WAVCNS (about 30 seconds in average). Another interesting difference between the BIS and WAVCNS is the difference of the index value during emergence. In particular, the WAVCNS was 90.6 (±10.2), while the BIS value remained in the range below 80 (79.0 ±18.8) at the time of the P r s event, i.e., when the patients were responsive. This result is comparable to the results reported by Sleigh et al. [79] who found BIS values under 80 while patients were responding to verbal command. 5.2.3.3 Behavior During Burst Suppression The WAVCNS was derived based on the methodology and data sets presented in Section 5.2.1. The frequency band and wavelet filter were selected to optimize the linearity of the index characteristic, while limiting its variability during periods of steady state. When considering the intermediate states T\, T2 and T 3 , the results presented in Fig. 5.8 show that the WAVCNS index has a nearly linear characteristic. However, this result does not necessarily hold true for other intermediate states. In particular, during burst suppression patterns, the EEG alternates between periods of fast, large-CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 81 S-waves burst-suppression with predominant isoelecticity burst-suppression with ! predominant bursts ! ->*^ >r*-gradual increase of 6 and a activity 3 4 T I M E [min] Figure 5.15: Time course of the BIS (v.3.4) and W A V C N S during an episode of burst suppression (patient #4). The corresponding E E G is also plotted for comparison purposes. amplitude waves, and periods of electrical quiescence. This non-stationarity of the E E G is a well known phenomenon, and has posed a challenge for assessing the patient's cortical state. Questions may therefore arise as to how the W A V C N S behaves during such episodes. Among the 20 patients of our study, 3 patients exhibited a short period of burst suppression. We present in Fig . 5.15 the time courses of the BIS and W A V C N S during such an episode, as well as the corresponding E E G . In this particular case, burst suppression with prominent isoelectric patterns occurred following a high-dose propofol induction (>5 mg/kg). This E E G pattern lasted 3 to 4 minutes during which the W A V C N S dropped to values of about 10. This level is consistent with the W A V C N S scale since this E E G signal corresponds to an intermediate state between T 3 and R^. As the periods of isoelectric E E G got shorter (possibly due to the elimination of the propofol from the blood plasma), the W A V C N S increased up to values of about 20. A gradual increase in 0 and a activity 5 minutes into this case provoked a further increase in the WAVCNS- Once steady state was reached (t > 6mm), both the W A V C N S and BIS time courses converged towards similar levels. During periods of burst suppression patterns, the W A V C N S m a Y exhibit an oscillatory behavior cor-responding to intermittent periods of bursts and electrical quiescence. In particular, we can observe an increase in W A V C N S during E E G bursts, and a decrease during periods of isoelectricity. This is consistent with the implementation of the index described in Section 5.2.2. The BIS time course does not seem to follow the changes in the E E G signal during periods of burst suppression in the similar way as the WAVCNS- This difference stems from the fact that the BIS monitor uses a different signal processing algorithm to deal with the non-stationarity of the E E G in this state. Conversely to the BIS, the W A V C N S algorithm is unique across the whole 100-0 scale, i.e., E E G epochs are processed similarly whether the patient is awake or experiences periods of bursts. We therefore expect the W A V C N S to differ from the BIS during burst suppression. However, note that further studies are needed to provide more clinical insight into the behavior of the W A V C N S during this particular state. CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 82 5.3 Measuring the Autonomic Nervous System (ANS) State: the W A V A N S Index Anesthesiology is commonly regarded as "the practice of autonomic medicine" [Lawson:2001], where the greater part of the anesthesiologist's expertise relies on his or her ability to control autonomic functions through the administration of drugs. Yet, there is currently no available monitor of the Autonomic Nervous System (ANS). Appropriate therapeutic actions are devised through the observation and interpretation of trends in patients' vital signs. The goal of this section is to close this gap by developing a feedback sensor for the direct assessment of the A N S state. This measure can then be used as a surrogate measure of the analgesia component of anesthesia. In particular, it is desired to obtain a bounded dimensionless index which would increase during au-tonomic activation (leading to heart rate and blood pressure increase) and decrease during autonomic depression (leading to heart rate and blood pressure decrease). In anesthetized patients, autonomic acti-vation can be obtained through noxious surgical stimulation and the use of anti-cholinergic drugs such as atropine. Conversely, the depression of the autonomic system can be achieved through the use of analgesic drugs, vapour anesthetics and adrenergic drugs such as /^-blockers. To derive this measure, we propose to use the H R V signal, and use the same signal processing technique as used previously for the derivation of the W A V C N S -We present results in the form of 3 case reports obtained from data gathered from clinical audits at the University of British Columbia Hospital. 5.3.1 The Use of Heart Rate Variability The H R V signal is derived directly from the E C G signal. The H R V is based on the R - R interval between two consecutive beats, see Figure 5.16.a. The last two decades have witnessed the recognition of the H R V as a promising quantitative markers of autonomic activity. The apparently easy derivation of this measure has popularized its use; from being an independent predictor of mortality following myocardial infarction, to detecting autonomic neuropathy in diabetic patients [161]. In terms of the quantification of the autonomic state, one particularly interesting aspect of the H R V signal is that its high frequency (HF) content is associated with vagal tone. When patients are under acute stress, such as surgical noxious stimulation, the vagal tone tends to be depressed, which is characterized by a decrease in H F activity. As a patient's body experiences increasing level of pain, beat-to-beat intervals tend to become more regular (see for instance Figure 5.16.b for illustration).'Conversely, during periods of relaxation, consecutive beat-to-beat intervals are more irregular. In addition to the effect of pain and noxious stimulation, there are a number of factors also known to affect the H R V high frequency content: - Respiration (or Respiratory Sinus Arrhythmia - RSA) : inspiration tends to increase the heart rate CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 83 R-R interval 1.2 0.8 'a\ 0£-Relaxed ai 0 100 200 300 400 V 1.2 r -0.4 0 1 2 Time [secj 0 100 200 #Beat 300 400 Figure 5.16: The HRV signal as a measure of autonomic activity, (a) 3 consecutive QRS complex (the HRV is obtained based on the detection of the R-wave and the measurement of the R-R interval between two consecutive R-waves). (b) Tachogram of a subject performing a relaxation exercise, and a patient reacting to surgical stimulation. while expiration tends to slow it down. As a result, the H R V exhibits an oscillatory behavior following the respiration pattern. Hence, changes in the respiratory pattern strongly influence the H F content of H R V and must be accounted for. - Opioids and inhalation agents: these agents are known analgesics. As such, they reduce the generation and transmission of pain signals and increase the H F of H R V . - /?-blocking drugs: these drugs are essentially used to depress the sympathetic system in order to prevent hypertension and tachycardia. However, evidence from the literature also support that the administration of (3-blocking drugs is associated with an in-crease in the H F content of the H R V [162]. - Anti-cholinergic drugs: drugs such as atropine antagonize the action of acetylcholine. This results in vagal blockade and the disappearing of H F in the H R V . The effect of pain, analgesic drugs, /?-adrenergic and anti-cholinergic agents all follow the same trend: pain and stress decreases the H F (associated with a rise in blood pressure and heart rate). Conversely, analgesics and /?-blockers increase the H F (associated with a decrease in blood pressure and heart rate). Therefore, the quantification of the H F activity of the H R V should provide a reliable metric of the effect of pain and analgesic/adrenergic drugs on the autonomic system. The fact that H R V analysis can yield a usable 'depth of anesthesia' index was already known in 1985. Unfortunately, the very nature of the H R V signal makes it particularly difficult to obtain a usable index while achieving real-time operation. There are a number of difficulties related to the acquisition and processing of the H R V signal: - E C G temporal resolution: the variability in terms of the R - R interval during periods of surgical stress can be as low as few milliseconds. As a result, the H R V should be derived based on an E C G signal CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 84 sampled at a high sampling rate (>500 samples per second) in order to gain good temporal resolution of the R-peak [161]. - Detection of the R-peak: the task of detecting the R-peak of the QRS complex is especially difficult during periods of artifactual activity. Electrosurgical Units (ESU) are the most significant source of interference in the operating room and completely obliterate the E C G [163]. This poses a severe obstacle for real-time intra-operative monitoring, since there is usually significant data loss during long periods of time. - Slow signal: probably the most damaging aspect of the H R V signal is the slow update rate linked to the heart beat. Any real-time analysis must be optimized in order to provide measurements which minimally lag relevant physiological events. In today's current state of the art, existing methods of 'short-time' H R V analysis usually require at least 5 minutes recordings to provide reliable information [161]. This type of performance is not suitable for a real-time application such as the one envisaged here. A part of the problem lies in the fact that neither time nor spectral analyses are suitable for non-stationary signals such as H R V . We believe that this difficulty can be overcome by using Wavelet Transform. 5.3.2 The W A V A N S Similarly to the W A V C N S index, the goal is to analyze the H R V signal in order to derive a metric, referred to as W A V A N S i representative of the patient autonomic state. This metric should represent the autonomic state of the patient, and be bounded between 2 extreme states: relaxed (i.e., no pain and comfortable) and extreme stress (i.e., high level of pain triggering a significant autonomous reaction). The W A V ANS index is obtained by applying to the H R V signal a methodology similar to that of the W A V C N S - A feature extracted from the H R V signal is compared to 2 reference features corresponding to the 2 extreme autonomic states. The W A V A N S is then the combination of the likelihood of the patient being relaxed, and the likelihood of the patient being under large surgical stress. To find the wavelet filter and frequency band which would carry information relative to the patient's autonomic state, we obtained 5 H R V data sets correspond-ing to different levels of stress, see Figure 5. The sets Ti and T 3 were obtained from 2 clinical patients during periods of surgical stress (T 2) and intense surgical stress leading to a significant autonomic response (T3). The sets Ra and T\ were obtained from a volunteer performing a breathing relaxation exercise (Ra, see Figure 5) and during a cold pressure test (i.e., moderate pain, T2). Finally, the set R^ is a simple constant value representing the absence of vagal activity. In all cases, patients and subjects were breathing at 8 breaths per minute (a metronome was used to guide the breathing pattern of the volunteer). A preliminary analysis based on these data revealed that the wavelet filter Daubechies #2 and the frequency band 0.2-0.4 Hz are optimal to discriminate between these 5 H R V signals (this correspond the d™ CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 85 R , T , T 2 T 3 R S (•) W Figure 5.17: Application of the WAV technology to the quantification of autonomic activity using the HRV signal, (a) Representative data sets of HRV signals for different stress levels: Ra: relaxed subject (y-scale: 0.2 s/div), T i : cold pressure test on volunteer (y-scale: 0.2 s/div), T^: clinical patient undergoing surgical stimulation (y-scale: 0.1 s/div), T 3 : clinical patient reacting to surgical stress (y-scale: 0.1 s/div). Rf. no vagal tone (i.e., no heart rate variability -theoretical state). (b) Result of the best wavelet and frequency band selection. Note that no post-processing filtering is required due to a long analysis epoch. high frequency band). The resulting W A V A N S metrics is expressed in the 0 - 1 0 0 % scale, where 0 % represents deep relaxation, and 1 0 0 % represents the total lack of parasympathetic activity. 5.3.3 Implementation Issues Even though the W A V A N S is based on similar principles than the W A V C N S J the nature of the H R V differs considerably to that of the E E G . As a result, the implementation of the W A V A N S needs to account for the various specificities of the H R V signal. HRV Signal Acquisition New H R V samples are only available at each heart beat. The R-peak of the E C G is used as a triggering event to measure the beat-to-beat interval. The temporal resolution of the R-peak is therefore particularly important to guarantee an adequate H R V resolution. In the current implementation, we use an E C G sampling rate of 960 S/s, which results in a temporal resolution close to 1 ms. A spike detection algorithm detects the peak of the R-wave. Artifacts are detected by analyzing the E C G signal between two consecutive R-peaks. If the beat-to-beat waveform differs considerably from a reference waveform, or if the beat-to-beat interval is much smaller or larger than previous measurements, an artifact flag indicating that the current R - R interval may be corrupted is activated. This flag is used in the processing algorithm to limit the effect of corrupted H R V samples on the W A V A N S -CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 86 D a t a flow In the current implementation, the H R V is sampled at the same sampling frequency than the E C G signal. A zero-order-hold is used to hold the R - R interval value between each heart beat. The signal is then resampled at a rate of 0.8 samples per second using a resampling filter bank. This choice of sampling frequency is motivated by the results found in Section 5.3.2. As compared to the W A V C N S implementation, the W A V A N S is updated at every new H R V sample. This results in an update rate of 1.25 second. Also, the S W T decomposition is calculated based on the last 2-minute of H R V data (96 samples) in order to have enough wavelet coefficients to derive a meaningful P D F . Note that, due to the large size of the analysis window (2 minutes), there is no need for further post-processing filtering. A r t i f a c t s Dealing with corrupting artifacts is also challenging. For E E G analysis, epochs containing samples corrupted by artifacts can be rejected without impairing significantly the availability and reliability of the W A V C N S - A similar strategy is not applicable for H R V analysis. It would indeed result in a large amount of data loss during, e.g., periods of electrocautery activity 1 , as it would be necessary to wait at least a complete H R V period (in this case 2 minutes) to obtain the next artifact-free H R V epoch. As a result, even H R V epochs containing corrupted samples are used in the W A V A N S calculation. It is therefore necessary to limit the effect of artifactual H R V samples on the W A V A N S calculation. In order to do so, wavelet coefficients obtained from corrupted samples must be removed from the P D F calculation. Because of the convolution operation, each corrupted H R V sample results in the removal of / wavelet coefficients, where I is the length of the high-pass wavelet filter. It is therefore particularly interesting to use low order wavelet filters to limit the effect of corrupted samples due to the convolution operation. In the current implementation, we use a 4 t h order F IR wavelet filter (Daubechies #2). Note that, if there are less than 32 valid wavelet coefficients left to carry out the P D F calculation, the H R V epoch is discarded. The W A V A N S value is then hold t i l l the next cycle. If more than 12 cycles (corresponding to 15 seconds) have elapsed since the last W A V A N S update, the monitor sends a warning message to the user and stop displaying the index. 5.3.4 Case Repor ts Three female patients undergoing hysterectomy procedures were induced with a combination of propofol and remifentanil (Case #1), or propofol and fentanyl (Case #2 and #3). After securing the airway, maintenance was insured by a background inhalational anesthetic supplemented by opioid boluses. In all cases, the patients were maintained in a hypocapneic state by mechanical ventilation (8 breaths/minute). Aside from the W A V A N S , all the data presented here were obtained on-line. The W A V C N S was obtained 'interferences due to electrosurgical units have been identified as the most perturbing operating room noise for ECG acquisition devices [163]. CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 87 directly from a 2-channel frontal electroencephalographic montage. The W A V A N S was calculated off-line based on a 3-lead E C G signal. Intra-operative events and the anesthesiologist's assessment of the patients' state were logged in by an independent observer. In Case #2 and #3, the observer also recorded manually the heart rate, blood pressure (each time this measurement became available), the M A C value and the anesthetic titration. Case #1 (see Figure 5.18) Following the initial propofol/remifentanil induction bolus, the patient's A N S state reached a low 'relaxed' level. Upon the start of surgery (multiple skin incisions), the patient's heart rate and blood pressure in-creased significantly. This increase was mirrored by a significant increase in the W A V A N S index. At 0h35, the patient was deemed 'light' by the anesthesiologist. Similarly, the W A V A N S reached a high value (>60). A t 0h45, a large bolus of fentanyl was administered, provoking a marked drop in the W A V A N S index. In 2 other instances, the anesthesiologist estimated that the patient was light (increase in heart rate and blood pressure). In these 2 instances, the W A V A N S also reached higher levels. From this case, it appears that a W A V A N S level of 70 (and above) is associated with a higher incidence of significant autonomic reaction to surgical stimulation. The large drop during emergence (2h35) can be due to a change in respiration pattern, as neuro-muscular blockade was reversed to allow extubation (the patient was spontaneously breathing at 2h30). Note that the W A V C N S reached a very low level (<30) during most of the procedure, thus indicating a rather deep hypnotic state. Following the indication of both C N S and A N S monitors, it appears that the sedation vs. analgesia balance was tilted towards the sedation endpoint. Ad hoc analysis of this case indicates that increasing the analgesia titration while decreasing hypnotic dosage would have been beneficial to this patient. The benefit of dual CNS and A N S monitoring is evident in this case. Case #2 (see Figure 5.19) The W A V A N S peaked 5 times during this case. The first 2 peaks (0h30 and 0h35) correspond to the most stimulating part of the surgery involving large skin incisions. Blood pressure and heart rate increased in response to the stimulation which prompted the anesthesiologist to give a fentanyl top-up dose. The next 3 W A V A N S peaks (0h55, lh lO, lh20) were also characterized by an increase in heart rate and were correlated to short episodes of electrocautery activity. In all instances, the duration of autonomic activation was very short (2-3 minutes) and may not have been captured by the blood pressure measurement (updated only every 5 minutes). Finally, at the time of 'adequate anesthesia' (as per the anesthesiologist's assessment), the W A V A N S was less than 50. Case #3 (see Figure 5.20) As compared to the previous 2 cases, the W A V A N S was, in average, less than 50 during most of the surgery. Similarly, there was no indication of a 'light' analgesic state during the whole procedure. Both the heart rate and blood pressure were stable. The absence of peak in the W A V A N S , or sudden increase in heart rate and blood pressure can be explained by the fact that this surgery is not as invasive as in the previous cases (limited skin incisions). Interestingly, the increase in W A V A N S towards CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 88 Figure 5.18: Case report #1 - Female patient, 55 years old, undergoing an hysterectomy procedure. The large gaps in W A V A N S index correspond to periods of time heavily corrupted by electrocautery artifacts (the H R V signal could not be obtained from the E C G signal). Electrocautery also affected the EEG signal and created data loss in the W A V C N S (indicated by vertical bars in the hypnosis index). The 'patient light' events correspond to the anesthesiologist's assessment of the patient during the surgery. Periods of intense electrocautery activity are indicated (it is assumed that they coincide with surgical noxious stimulation). the end of the case mirrors the decrease in M A C value. This is consistent with the fact that M A C is a measure of the analgesic potency of the administered vapour. Upon emergence, the patient reacted to the anesthesiologist's voice, but was not able to respond to verbal commands. These 3 cases illustrate how surgical stimulation triggers abrupt changes in heart rate and blood pressure. These abrupt changes are used by anesthesiologists to identify the presence of excessive autonomic acti-vation, and then to proceed to take appropriate pharmacological action. Conversely, the W A V A N S is a bounded and population-normed metric. As such, it can be used as a predictor of autonomic activation. For instance, in Case #1, the trend of the index shows that the patient's state is getting increasingly lighter (from 0h50 to lh40), up to the point where the patient did actually react to the surgical act. Also, the index was quite high during the second half of the surgery, indicating the need to deepen the analgesic state. Conversely, in Case #3, the index was quite low and steady throughout most of the surgery. This is consistent with the rather aggressive sevoflurane titration (MAC>1.0) in a surgery involving only minor skin incisions and moderate stimulation. This may also explain the prolonged emergence. Finally in Case #2, the number of W A V A N S 'spikes' and their rather short duration may indicate an overall shallow analgesic state that resulted in strong and acute autonomic activation whenever surgical stimulation was present. These preliminary results are very encouraging. It appears that the W A V A N S has the potential to C H A P T E R 5. QUANTIFYING CORTICAL AND AUTONOMIC A C T I V I T Y USING WAVELETS 89 SYS BP (mm Hg ) i 120 1 fit t - (b. - {(it 1IA* H O Induction Fentanyl bolus ^ . Fentanyl : „ bolus jSurgery B p Adequate Patient anesthesia Patient light responding Hydrotflorphone M A C 1.4 1.2 -1.0 - 100 0.8 _ 80 eg 0.6 ^60 0.4 < ^ 40 0.2 20 0.0 - 0 OOhOO OOhlO 00h20 00h30 00h40 00h50 OlhOO T I M E OlhlO 01h20 01h30 Figure 5.19: Case report #2 - Female patient, 39 years old, undergoing a hysterectomy procedure. The heart, rate (top graph) was calculated ad hoc based on the E K G . The blood pressure from the cuff sensor was recorded manually each time this measurement was updated by the anesthesia monitor. CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 90 OOhOO OOhlO 00h20 00h30 00h40 00h50 OlhOO OlhlO T I M E Figure 5.20: Case report #3 - Female patient, 38 years old, undergoing a laparotomy procedure. Due to the nature of the surgery, skin incision was kept to a minimum. provide a guide for analgesic titration and guarantee patients' safety from a hemodynamic standpoint. It is important to point out, however, that this study suffers from a number of limitations: - First, the W A V A N S was calculated off-line. Artifactual periods were visually reviewed in order to facilitate the detection of the R-peak and limit the amount of data loss. Unfortunately, the limitation of the front end amplifier in attenuating the ESU interference resulted in large gaps in the H R V signal, which precluded the calculation of the W A V A N S for long periods of time during the most important part of the surgery. - Also, the 3 patients were mechanically ventilated, i.e., the effect of a change in the respiration pattern cannot be inferred from this data. The analysis algorithm was specifically tuned for 8 breaths a minute, so a different respiratory rate may have resulted in different W A V A N S levels {e.g., the abrupt W A V A N S drop at the end of Case #1 can be attributed to an erratic and rapid breathing pattern). - Finally, this audit was simply an observational study, i.e., the administration of opioids and change in gas concentration often coincided with surgical stimulation. It is therefore difficult to identify the contribution of noxious stimuli from the effect of opioids from the analysis of the W A V A N S -It is evident that more clinical data obtained under precisely defined clinical protocols are needed in order to further develop the proposed W A V A N S index and assess its robustness in a clinical setting. In particular, it CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 91 will be important to clearly demonstrate the existence of a dose/response relationship between the W A V A N S and drug dosage. 5.4 Dynamic Behavior We conclude this Chapter by deriving deriving transfer functions for both indices in order to model their dynamic behavior in response to changes in the patient's state. W A V C N S Transfer F u n c t i o n The elements dictating the W A V C N S dynamic behavior are the second order IIR low pass filter used to extract the main index trend (characterized by a cutoff frequency wo and a damping factor C), and the 1-second update rate (zero-order-hold element). As a result, the sensor dynamics can be represented by the following LTI transfer function: Hcm(s) = , 0 / ° — j • — — (5.24) s z + 2 • C • wo • s + OJQ S W i t h the current W A V C N S implementation, £ = 1 and UJQ = 0.02 Hz: " C N S ( S ) = (8^W • ^  (5'25) If we assume no a priori knowledge about the W A V C N S sensor dynamics, an identification procedure can be carried out using a database of various E E G signals corresponding to known sensor values. To illustrate this procedure, we consider the 5 data sets used for the derivation of the W A V C N S in Section 5.2.1.3. Each of these data sets corresponds to a stationary hypnotic depth: W A V C N S : {Ra, Tu T2, T3, Rb} -+ {98,78, 56,19,0} [%] (5.26) A composite E E G signal is first obtained by arbitrarily selecting E E G epochs from each of the 5 data sets, see Fig . 5.21. A n identification input signal is then obtained by attributing to each of the E E G epoch the corresponding W A V C N S value according to (5.26). Finally, the composite E E G signal is processed by the sensor algorithm to yield an output identification signal. The sensor dynamic ffcNs(s) is further identified by a least square estimation algorithm. In this case, the identification procedure resulted in a second order transfer function plus delay: ^ w = a»- .+*-Ur ' .+ !)••"• ( 5 2 7 ) The difference between the identified (5.27) and the analytic (5.25) functions is negligible in the low fre-quency range, see Fig. 5.22. Differences are more apparent for higher frequencies. This discrepancy results from the limited input identification signal bandwidth (<1 r a d s - 1 in this case) and the identification method whose accuracy can not be guaranteed near the 1 Hz sampling frequency. Since it is unlikely that the controller cross-over frequency be placed in this frequency range, no further refinements are necessary. C H A P T E R 5. QUANTIFYING CORTICAL AND AUTONOMIC A C T I V I T Y USING WAVELETS 92 100 r 0 50 100 150 200 250 300 350 T I M E [S] Figure 5.21: Experiment design for identifying the sensor dynamics. The E E G signal is composed based on E E G epochs arbitrarily chosen from a database of signals. The input identification signal corresponds to the simulated patient's instantaneous state (direct correspondence with the source E E G signal), while the output identification signal is the W A V C N S -CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 9 3 > 50 o 0 « w -50 100 cn 80 Z U 60 > s 40 20 0 - W A V C N S i(measured) W A V C N S l(predictcd) / - f 3 4 Time [min] (a) 10 15 20 Time [min] (b) Figure 5.23: Accuracy of H C N S and -HANS for predicting the W A V C N S and W A V A N S time courses calculated for various step changes, (a) The -ffcNS is able to predict accurately the W A V C N S time course, (b) The -HANS is an approximation of the index dynamics. As such, some discrepancies between measured and predicted outputs exist. Also, there is a noticeable index variation in T 3 . This can be the result of the fact that the patient state was not as stable as in the other states.) W A V A N S Transfer Function Following the proposed implementation, the W A V A N S dynamic is equiv-alent to a 2 minutes moving average due to the selected length of the H R V period upon which the analysis is carried out. Considering that the index is updated every 1.25 seconds, the W A V A N S transfer function can be expressed as: 1 95 y» g-i-25-, # A N S ( S ) 1 96 .1 - e -1.25s (5.28) fc=0 The moving average filter can further be approximated as a first order transfer function whose time constant is one half the moving average window: HANS(S) 1 (60 • s + 1) 1 -e ' -1.25-s (5.29) The performance of H C N S and . H A N S was assessed by comparing the time course of the indexes vs. their predicted time courses, see Figure 5.23. This test was done based on source signals (EEG or R R ) obtained from a combination of the initial data sets. The HCNS was found to be remarkably accurate in predicting the W A V C N S time course. Conversely, the transfer function H A N S is only a coarse approximation of the W A V A N S dynamics. This result can be improved by reducing the W A V A N S analysis window and adding a post processing I IR trending filter. However, this would reduce the number of wavelet coefficients, which might in turn decrease the resolution of the PDF waveform and increase the measurement noise, while making the index more sensitive to artifacts. CHAPTER 5. QUANTIFYING CORTICAL AND AUTONOMIC ACTIVITY USING WAVELETS 9 4 5.5 Summary The search for anesthesia feedback sensors has motivated the development of a methodology to analyze and characterize non-stationary physiological signals. This methodology, referred to as W A V , uses wavelet decomposition to capture rapid signal changes in both time and frequency. The statistical representation of the wavelet coefficients using probability density functions allowed us to derive bounded metrics, thus quantifying the signal into a pre-defined scale. The flexibility of wavelets can be used as an advantage to tune the static characteristic of the index in order to reflect the anesthesiologist's assessment. We applied this methodology to the analysis of E E G signals in an attempt to quantify patients' cortical activity. The resulting index, the W A V C N S , has demonstrated superior dynamic performance as compared to today's leading monitoring technology. Also, to our knowledge, the W A V C N S is the only index of cortical activity whose dynamics are linear and can be expressed as an LTI transfer function. As such, it is an attractive solution for use as a feedback sensor for advisory and close loop systems, and to identify the pharmacodynamic models of anesthetic drugs. The W A V methodology has also demonstrated strong potential for use in the analysis and characteriza-tion of the H R V signal. The case reports presented here show that the W A V A N S was representative of the anesthesiologist's assessment of the patient's analgesic state. In particular, all intra-operative events leading to a significant change in heart rate and blood pressure were well correlated with an elevated W A V A N S value (for instance, a W A V A N S level of 7 0 was associated with a higher incidence of significant blood pressure and heart rate change). In addition, the administration of opioids (fentanyl, hydromorphone) was followed by a decrease in the W A V A N S index. Similarly, changes in M A C value seemed to be well correlated with changes in W A V A N S -It is likely that other applications can benefit from the W A V technology, e.g. automatic sleep scoring, depression and Alzheimer's decease diagnosing, etc... As of today, our Group has invested considerable time and efforts into the development and validation of the W A V C N S - The intellectual property is now protected [164], and a commercial device integrating this technology should become available in mid-2006 (NeuroSENSE™ Monitor, Cleveland Medical Devices Inc., OH). Chapter 6 A System Oriented Approach to Pharmacodynamic modeling We exposed in details in Chapter 3 the approach to pharmacokinetic (PK) and pharmacodynamic (PD) modeling favored by pharmacologists. In terms of pharmacokinetics, the models proposed in the literature are consistent. The use of the N O N M E M analysis software has enabled clinicians to improve the fit of P K models by including patients' biometric data (weight, lean body mass, age, etc.). The inclusion of P K models within commercial T C I systems are a testimony to their reliability and acceptance from the clinical community. In terms of pharmacodynamics, the situation is different. The latest P D studies carried out for intra-venous anesthetic drugs have yielded contradictory results and large inter-individual variability in some of the parameters. Furthermore, P D models in the literature were derived specifically for use with one particular sensor, and thus, cannot be readily used to predict the time course of an effect measured by different sensor. It is our belief that the traditional modeling approach deserves to be revisited from a system engineering perspective. In this chapter, we propose a system oriented approach to P D modeling. This approach is detailed in the first section following a discussion of the limitation of the traditional approach. The new approach is then applied to the modeling of propofol vs. W A V C N S using the induction data from 44 patients undergoing various ambulatory surgeries. The traditional approach is also carried out for comparison purposes. We show that the system oriented P D approach yields better fit and predictive performance. Further-more, the resulting P D model is independent of the sensing technology used to observe the effect. Finally, we also show that the traditional approach can lead to strongly colored residuals, which, in turn, affects the accuracy of the static dose vs. response characteristic of the drug. 6.1 A New System Oriented Approach to PD Modeling This section presents in detail a P D modeling approach based on standard system engineering know-how. 95 CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 96 PD(s) c, \ Quantified effect Pharmacodynamic Model Figure 6.1: Traditional PD model block diagram. To better grasp the fundamental differences between this system-oriented approach and the traditional pharmacological approach (also referred in the following to as Sheiner's approach), we first bring a critical look to the traditional approach previously discussed in Chapter 3. 6.1.1 Sheiner's Approach The approach described in Chapter 3 can be summarized under the form of the block diagram of Fig. 6.1 where the P D model expresses the relationship between a quantified effect Eq and the plasma concen-tration Cp of a given drug. This relationship is thus expressed as a Wiener type model, where a Linear Time Invariant (LTI) function is followed by a non-linear element to capture the static dose vs. response relationship. In this approach, PD(s) is a first order transfer function with a unity gain. It expresses the drug concentration Ce within the target organ as a function of the plasma concentration. The transfer function PD(s) is uniquely defined by the effect site rate constant keo, which is the rate at which the effect-site concentration equilibrates with the plasma concentration, see (3.13). The addition of this function to the classical dose/response curves was proposed by Sheiner et al. [114] in 1979 in an attempt to collapse the 'hysteresis' loop observed when plotting the effect vs. the plasma concentration during rapid administration of the drug, see Fig . 3.9.b. The non-linear function which relates the quantified effect to the effect site concentration is handled by a Hi l l equation in order to capture the sigmoid characteristic of the saturation of the dose vs. response curve. This saturation can be particularly severe when considering an effect that can only be measured over a limited portion of the drug therapeutic window. For slowly varying plasma concentrations, PD(s) can be neglected, in which case the P D model can be simplified to the classical dose vs. response curve expressed by the Hi l l equation. While Sheiner's approach has been the standard in P D studies since 1979, P D models are known to suffer from a limited accuracy, large patient variability and strong non-linearity. In particular, the traditional P D approach suffers from three major limitations, which might explain its overall poor performance: 1. use of a measured input signal even though this signal is not readily available, 2. use of a limited L T I model order that fails to capture all of the system dynamics, and, 3. inherent assumption that sensor dynamics are negligible. CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 97 These limitations are discussed in more detail below. For illustration purposes, we consider the 6 published P D studies reported in Table B.2 which involve the effect of propofol on electroencephalographic variables. These studies are recent (1997-2002) and originate from different research centers (USA, U K , Japan and The Netherlands). As such we believe that they are a good representation of the current state-of-the-art in P D modeling practices. P r e d i c t e d vs . M e a s u r e d P l a s m a Concen t r a t i on In most studies, the plasma concentration Cp is measured through blood sampling. There are different scenarios as for the use of these blood samples. Some authors (e.g., Schnider et ai, White et ai, and Ludbrook et al.) first derive their own P K model in order to predict the time course of the plasma concentration, which they further use with the corresponding measured effect to identify the P D parameter set. This method is referred to as the parameteric method. Others (Schnider et ai, White et ai, and Kuizenga et al.) use a connect-the-dots method to extract the time course of Cp and adjust the keo value to collapse the hysteresis loop, and further, identify the dose vs. response curve. In particular, Kuizenga uses linear interpolation between two consecutive measured samples during increases in blood concentration, as well as a logarithmic interpolation between measured concentrations during decreases in concentrations. This method is usually referred to as the non-parametric method. While both methods seem appropriate, the blood plasma concentration of the drug is usually not available in real time. Hence, the advantage of the parametric method lies in the fact that a P K model is given in order to predict the plasma concentration time course. The method can be used to implement on-line algorithms for drug effect prediction. However, questions can arise as to the validity of P K models derived on a very limited number of patients, as compared to dedicated P K models derived in studies involving hundreds of subjects and using the latest identification methods. As for the P D models developed by the non-parametric method (i.e., where the input identification data are measured instead of predicted), one can reasonably argue that these models are probably closer to expressing the true pharmacodynamics of the drug. However, their usefulness remains limited as long as no direct real-time plasma concentration measurement is available. In a minority of P D studies, authors use third party P K models to predict the plasma concentration time course. For instance, Kazama et al. have used the Gepts et al. propofol P K model (see Table B . l ) . The blood samples were then used to validate the P K model output. Cases which strayed too far from the predicted plasma concentrations were rejected from the analysis, thus yielding a more consistent P D parameter set. It appears that the parametric approach is the only approach yielding direct applications (providing the P K parameter set is disclosed by the authors). However, the limited amount of data used to derive the P K models can be a significant factor limiting the reliability of the overall P K P D model. In that respect, the effect of the P K model in P D studies has been discussed in a 2003 study by Minto et al. [165]. The authors concluded that substituting the P K set used to predict the plasma concentration in the P D study CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 98 with another P K set will yield significantly different results, mostly during transient phases. Inadequate L T I M o d e l S t ruc tu re The effect vs. plasma concentration 'hysteresis' phenomenon has received considerable attention since the original work of Sheiner. The idea of collapsing the hysteresis loop 1 using an additional fictitious P K effect-site compartment can be particularly attractive to clinicians and pharmacologists, as it provides a physical meaningful explanation for the effect lag, and a way to mathematically express it in the well-known framework of compartmental mamillary models. In Sheiner's approach, one would try to collapse the 'hysteresis' loop through a first order transfer function with unity gain (3.13). We can see in Figure 6.2.c-e how the selection of the model time constant keo affects the width of the loop. As keo decreases, the loop collapses and then inverts. When keo equals the system time constant, we obtain a perfect equilibration between the measured and predicted time courses. However, as with any identification procedure, the goodness of the fit between measured and predicted time courses depends both on the quality of the identification data, and the selection of the model order and structure. In the case presented in Figure 6.2, the model has a similar structure to that of the system, which explains the success of this approach. However, when the system dynamics are more involved than that of the model, the traditional approach can easily yield poor results. In particular, a mis-modelled time delay can result in a smaller time constant keo, since this constant now captures both the system first-order and time delay dynamics. This is illustrated in Figure 6.3. Another important implication is that the residuals (i.e., the error between the measured and predicted time courses) may be strongly biased. This can result in a stronger non-linear characteristic if a Hi l l non-linear saturation is derived to further minimize these residuals. For instance, in the example of Figure 6.3, a non-linear H i l l saturation was added to improve the goodness of the fit between the measured and predicted outputs. Even though the system is linear, the model now includes a non-linear element, which is a direct consequence of not including all of the linear dynamics in the LTI part of the model. A n important consequence of this result is that P D models derived based on Sheiner's approach may only provide an adequate fit for input signals similar to those used during the identification procedure. If a different input signal is presented to the system, the model will fail in predicting adequately the system output. For instance, it is clear in the example presented in Figure 6.3 that the system output will not be accurately predicted for a slow varying low amplitude input signal due to the Hi l l saturation. It appears therefore particularly important that a thorough analysis of input/output data be performed ^ote that the term 'hysteresis' is used here inadequately since the width of the loop depends on the rate of administration of the drug. For instance, the loop is much wider for bolus administration than for a constant drug infusion. It will then be clear to anyone well-versed in system modeling and identification that the 'hysteresis' phenomenon results mostly from unaccounted dynamics between the plasma concentration and the effect, rather than from a true non-linear hysteresis in the system. This can be easily illustrated in the example presented in Figure 6.2, where we consider a simple first order system. This system is subject to a given input signal. Using a measurement of the output, it is possible to derive an input vs. output plot which shows the disequilibrium generated by the system first order dynamic, see Figure 6.2.b. Even though this plot resembles to that of an hysteresis, there is no such element in the system, which is actually perfectly linear. CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 99 Measurement noise Model / / \ / T V Measured output signal • Predicted output signal Input signal (b) 3 O 1 =0.2 S" 1 r / / L m It*** Predicted output signal (c) 3 O -0 I 0=0.1 s-1 — - A Predicted output signal 3 O -O I rt U 0=0.0 5 s-1 .„ Predicted output signal («0 Figure 6.2: Example of the 'hysteresis' phenomenon observed between the drug plasma concentration and the corresponding effect, (a) In this example, an output signal is obtained from a first order system, (b) Plotting the output vs. input relationship yields the 'hysteresis' curve described by pharmacologists. By tuning appropriately the rate constant fceo of a first order model, one can collapse the loop, (c) Under-compensation: the rate constant is too large, (d) Perfect equilibration. Note that the measured vs. predicted characteristic is linear (i.e., in this case, a Hill saturation would not improve further the fit), (e) Over-compensation: the rate constant is too small, which results in an inverted loop. CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 100 Figure 6.3: Effect of a pure time delay, (a) A pure 10 seconds time delay was added to the true system of Figure 6.2. (b) Plotting the output us. input relationship yields a 'hysteresis' curve similar to that of the previous example (i.e., there is no discriminating feature which can stress out the presence of a delay), (c) Under-compensation. (d) Equilibration. Note, however, that the rate constant is now about 6 times smaller than the true time constant of the system. The first order model therefore mostly captures the time delay dynamics. The model is just a coarse approximation of the true system. Even though the 'hysteresis' loop is reduced, the error between the predicted output and the measured output becomes quite large, (e) Over-compensation. CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 101 PR's) \ PD(s) \ \ H(s) \ Scaling Junction —+  E l Quantified effect Pharmacokinetic Model Pharmacodynamic Model Sensor Figure 6.4: Block diagram of the proposed system oriented approach to PD modeling. to determine the optimal L T I model structure. The non-linear saturation should only be added once all of the system linear dynamics are accounted for by the LTI part of the model, that is, once the residuals are white and show no correlation to the input signal. While the example of Figure 6.3 illustrated the effect of a time delay unaccounted in the model structure, similar results are obtained when considering higher order systems. Note that Sheiner et al. in their original 1979 paper reported that a limitation of their technique is that "its kinetic simplicity may cause it to be inadequate to describe complex pharmacodynamics". Sensor D y n a m i c s A major part of the drug effect dynamics can be directly related to the sensor dynamic whose output filter extracts the main trend out of the raw data. Sensors typically introduce delays, low pass filtering, and sampling dynamics. They might also be responsible in part for the non-linearity of the dose vs. response relationship. In Sheiner's approach, these dynamics are implicitly included in the P D model. As a result, a P D model is limited to the feedback sensor used for its identification. Hence, each time a new sensor is introduced, the P D identification procedure must be repeated. Also, any change in the signal processing algorithm of the sensor will make the corresponding P D model obsolete. Out of the 6 P D studies, 4 used the BIS index, 1 used auditory evoked potentials and 1 used the semilinear correlation index. Differences between the P D parameter sets can be explained by differences in these feedback quantities. Hence, only the 4 studies involving BIS can be reasonably compared. This aspect of pharmacodynamics was already alluded to by Kazama et al. who tried to remove the influence of the BIS averaging filter by shifting in time the BIS data in order to synchronize them with the plasma concentration time course predicted by the T C I pump. Even though the resulting P D model does not exactly predict the BIS time course, it is more consistent than P D models obtained in other studies (lesser patient variability). 6.1.2 Proposed Methodo logy To alleviate these limitations, we propose to revisit pharmacodynamic modeling from a more systematic system engineering point of view. As a guide for this discussion, we consider again the example of the effect of an anesthetic drug onto the brain. While the exact mechanism of action of anesthetic drugs is still debated, it has been shown that drugs CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 102 such as barbiturates, benzodiazepines and propofol exert their pharmacological action at the level of the 7-aminobutyric acid ( G A B A ) receptors. The G A B A is one of the major inhibitory neurotransmitter in the brain. Its role is to counterbalance the influence of excitatory neurotransmitter. B y mimicking its structure, anesthetic drugs bind themselves to the cell receptors, thus sending false information and shifting the excitatory balance towards relaxation. Hence, the effect of anesthetic drugs can be well understood when considering the percentage of bound G A B A receptors of each individual cell in the brain. Depending on the proximity of blood vessels, and thus the bio-availability of the drug, this percentage is going to evolve in time. Unfortunately, the effect at the cellular level cannot be readily measured. However, the combination of each individual cellular effect will result in a number of macroscopic effects which can be observed and quantified. For example, as an increasing number of receptors is occupied by the drug, loss of consciousness and memory formation can be observed. The depression of the C N S will result in a drop in heart rate and blood pressure. Similarly, changes in the E E G patterns can be observed. These E E G changes can be further quantified into a unique value {e.g., the W A V C N S ) through the signal processing technique described in Chapter 5. Hence, when considering the pharmacodynamics of anesthetic drugs, one can consider the effect of a given drug concentration in the blood perfusing the brain onto the E E G signal itself. We can take here the assumption that this relationship can be adequately described by a LTI transfer function followed by a non-linear function. The LTI element captures the dynamics between the percentage of bound G A B A receptors and the plasma concentration. The non-linear function, in turn, models the saturation phenomenon proper to E E G changes: a low number of bound G A B A receptors will have no quantifiable effects on the E E G , while no increase in the number of bound receptors will provoke any further change in the E E G when the complete suppression of cortical activity is reached. The patient system can thus be modelled using the block diagram of Fig . 6.4. As compared to the traditional approach, the plasma concentration Cp(s) is obtained using a well defined P K model derived for arterial blood concentration. Furthermore, the pharmacodynamics are now limited to the effect of the drug onto an observable physiological signal (and not its quantified correspondent). In this scheme, the sensor dynamic is thus a separate entity all together. Also, no a priori decision is made as to the structure of PD(s), which can be any reorder plus delay transfer function: p D ( s ) = E 1 { s ) ^ e ^ . s b m - s m + b m . 1 - s ^ + . . . + b 1 - s + ao^ 1 w Cp(s) s n + o „ _ i • s""1 + . . . + ax • s + an 2 • E C 5 0 ' ~  y ' J Residual analysis should be performed to validate the selected model structure. In this approach, the pharmacodynamic model establishes the relationship between the observable effect E0bs and the plasma concentration. E0bs is bounded between 0 (no effect) and 1 (maximal effect which can be observed in the E E G signal). The term EC50 represents here the steady-state plasma concentration necessary to obtain 50% of the measurable effect. Conversely to the traditional approach, the non-linear element does not express the dose vs. response relationship. It only models the saturation of the effect CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 103 with respect to the dose. This saturation can be expressed using the Hi l l relationship, which is essentially a non-linear function that can be represented as a variable steady state gain in the frequency domain, see Section 6.3.4 and Section 7.2.1.2: El Eobs = o 5 7 + E i =* 0<Eobs<l (6.2) Note that the non-linear element is defined here entirely by the steepness coefficient 7. Finally, the quantified effect Eq is obtained through filtering and scaling. The sensor filter H(s) is essentially a unity gain element, while the scaling function is an affine function defined by two scaling parameters: Emax and EQ as follows: Eq = {Emax — Eo) • E0bs + EQ (6.3) It is interesting to note that, as compared to Sheiner's approach, the static dose vs. response curve is obtained by considering the non-linear function (6.2), the static gain of PD(s) and the scaling function of the sensor: C 7 Eq,ss = {Emax ~ EQ) • Y ' 7 b Eo (6.4) where CPiSS is the steady state blood plasma concentration of the drug. Similarities between the traditional approach and the new system oriented approach are obvious. The traditional approach can be viewed as a simplification of the proposed approach, where the non-linear function has become an output non-linearity, and where the traditional PD(s) now captures both the effect and the sensor dynamics. However, it is clear that this is an over-simplification, which results in significantly degraded modeling performances. 6.2 Application to Propofol and W A V C N S Index In this section, we derive a P D model for Propofol following both the traditional P D approach, and the new system-oriented methodology presented in Section 6.1.2. The models derived here aim at predicting the time course of the W A V C N S following intravenous propofol administration. 6.2.1 C l i n i c a l D a t a To identify the pharmacodynamic model of propofol, we are using here the induction data collected during the L M A study (see Annex D.2). In this study, the W A V time course during propofol induction was observed in 76 patients. The L M A insertion can provoke a strong reaction but is not as stimulating as the endotracheal intubation. The L M A study has revealed that patients that do react to the L M A insertion (jaw tone, coughing, swallowing, grimacing, etc.), had also a significantly different W A V time course. Furthermore, most of the reaction cases elicited an E M G activation, which resulted in an elevation of the W A V during the insertion itself. Hence, the L M A insertion can be considered to be a typical output CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 104 disturbance. Among the 76 patients of the study, 51 patients aged 18 to 60 years did not elicit any strong airway reaction to the insertion. To avoid having the L M A reaction acting as a disturbance in the output identification signal, only these 51 cases were considered for this analysis. In the L M A study, the start and end of the propofol injection were recorded as well as the L M A insertion time. Even though gaseous anesthetic drugs were given once the patient's airway was instrumented, it is usually assumed that these drug have only a marginal effect in the first 2 minutes following their administration. Hence, the identification window was chosen to start at the propofol injection, and to stop 90 seconds after the L M A insertion. Three cases had insufficient data to cover the whole identification window. These cases were withdrawn from the analysis. In order to keep consistent results, it was further decided to eliminate all cases whose W A V during induction reached values lower than 10. It is indeed possible that a saturation phenomenon can affect the W A V time course in the lower W A V values, which in turn would result in a significantly different dynamics. Two cases were excluded based on this criteria. Finally, two more cases were withdrawn from the study due to a large E M G activation observed after the L M A insertion. These cases were characterized by strong extremity movements, but no airway reaction. A total of 44 cases were thus compiled for the identification procedure. As compared to other published P D studies, this is the second largest study in terms of patient population. These cases were further subdivided into 4 age groups (18-29 yrs; 30-39 yrs; 40-49 yrs and 50-60 yrs) to assess if age represents a covariate factor in the P D models. Note that a potential limitation of this study is the co-administration of a fentanyl bolus as part of the induction sequence (usually about 1 minute prior to the propofol administration). In this analysis, we neglected the effect of this fentanyl dose. Even though some synergism might have occurred, it is a well documented fact that the propofol/opioid synergism is not very acute in terms of hypnosis. 6.2.2 P D Identification - Traditional Approach The identification procedure for the traditional P D approach can be divided into three stages. A n illustrative example is given in Fig . 6.5.a. Stage #1: Plasma Concentration Profile To identify the pharmacodynamic parameters, it is neces-sary to first calculate the time course of the propofol plasma concentration based on the infusion profile. To achieve this, a number of P K models can be used. We select here the parameter set from Schiittler and Ihmsen (see Section 3.1.3) since it is derived based on a large adult population (20 to 60 years) and incorporates both the age and weight as covariates. Also, one advantage of Schiittler's P K set is that the parameters can be adapted for a bolus-type administration, which is necessary when using induction data. Stage #2: Effect Dynamic Identification Once both the input Cp and output W A V C N S data are defined, the effect dynamic PD(s) is identified based on a least square optimization algorithm, where the CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 105 Figure 6.5: PD identification illustrative example (patient #52). (a) Traditional approach (b) System oriented approach square of the error between the measured and predicted effects is minimized. This yields the effect rate constant keo. This method is essentially identical to that of the 'hysteresis' collapse method. Once keo is identified, PD(s) is defined such as in (3.13) and used to estimate the effect-site concentration Ce. Stage # 3 : H i l l Pa ramete rs Ident i f ica t ion The Hi l l parameter set {EC50,7} is finally identified using Ce and the measured W A V C N S values. A search algorithm compares the predicted vs. measured W A V C N S and selects the Hi l l parameters to minimize the root mean square (RMS) of the residuals. During the identification analysis, the R M S of the residuals is used as a minimization criterion to optimize the fit between the predicted vs. measured effect. Other criteria, such as the performance indexes presented in Section 6.3.2, can also be used. 6.2.3 P D I d e n t i f i c a t i o n - N e w A p p r o a c h The identification procedure for the system oriented approach is more convoluted as it involves knowledge of the sensor dynamics and the identification of a time delay. A n illustrative example is given in Fig . 6.5.b. Stage # 1 : P l a s m a C o n c e n t r a t i o n Prof i le Similarly to the traditional approach, Schiittler's P K models are used to predict the time course of the propofol plasma concentration. CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 106 c,C<> c/(,> \ *• PK(s) • • J • E,.,W l/iura/W WAVas effect Knownfunctions Figure 6.6: Block diagram used in stage #3: the non-linear element is omitted and the sensor filter is used as an input filter. Stage # 2 : W A V C N S Dynamics Conversely to the traditional approach, the sensor dynamics are now a separate subsystem, which can be expressed as a transfer function i?cNs( s ) a n d a scaling function. In that respect, a major advantage of the W A V technology is that it allows us to express the W A V C N S dynamics as a linear time invariant transfer function. In particular, we found in Section 5.4 that: H c M = / ° • (6.5) where UJQ = 0.02 Hz and C — 1 (current implementation). Stage #3 : Effect Dynamic Identification Once the sensor filter i?cNs(s) is defined, the block diagram of Fig . 6.4 is rewritten as in Fig. 6.6, where the non-linear saturation function is omitted. The only unknown element is the transfer function PD(s) whose output Eq<u can be directly obtained from the W A V C N S values by applying the inverse of the scaling function. In this block diagram, we use the property of linear time invariant systems to bring the sensor filter as an input element. The transfer function i ? c N s ( s ) is then used to obtain a filtered version Cp of the plasma concentration. Using both C / and Eq%u as input and output identification signals respectively, the parameters of PD(s) can be estimated using a classical least square identification method. One difficulty arises in the determination of the order of PD(s). W i t h data such as the L M A data, it would be difficult to identify a high order function because of the limited input excitation. In this case, we choose PD(s) to be a first order plus delay function: PD(s)= E r ^ =e-TdS-— - (6 6) r i J [ S ) Cp(s) s + kd 2-EC50 [ b - b ) where kd is the effect dynamics and EC50 is the effect site concentration that yields 50% of the maximum observable effect. This choice is further discussed in Section 6.2.6. Stage #4 : Hi l l Parameters Identification Once the PD(s) function is identified, the system block diagram is returned to its original form of Fig. 6.4. At this point, only the non-linear element parameters are unknown. Similarly as in the traditional approach, a search algorithm determines the Hi l l parameters that minimize the R M S of the residuals. CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 107 6.2.4 Adequacy of the Identification Data When identifying the model parameters of a given system, it is usually recommended to design an experiment which will excite all the modes of the system in order to yield unbiased and rapidly converging parameter estimates. A large amplitude white noise is considered to be the ideal identification input signal, since it contains many different frequency components susceptible to stress all of the system dynamics. However, it is not always possible to subject a system to inputs such as Pseudo-Random Binary Sequences (PRBS) or filtered Gaussian signals. In the particular case of pharmacodynamics, the needs of the anesthesia procedure and surgery always supersede the needs of the identification itself. It is therefore important to assess a priori the adequacy of the input signal for estimating the LTI parameters of the models. In the induction data provided by the L M A study, the propofol is given as a single bolus over a given period of time. Using Sheiner's approach, the identification uses Cp(t) as an input signal to determine the time constant keo. This time course can usually be divided into two distinct parts: an initial large amplitude fast transient following the drug uptake, and a slower decreasing waveform corresponding to the initial distribution of the drug. The initial fast transient excites the system in a bandwidth roughly equivalent to 4 times the propofol administration time. This corresponds to a bandwidth of about 0.04 to 0.15 rad-s - 1 across the whole study population. In the proposed approach, we are using a filtered version of the plasma concentration to identify the gain, the time constant and the delay of PD(s) as defined in (6.6). This also roughly corresponds to an excitation bandwidth of 0.04 to 0.15 rad-s - 1 . Note that the second part of the input transient mostly reveals the steady state characteristic of the system. However, the identification window was too short for the system to settle. Therefore, a larger error in the dc gain is expected. The identification data provided by the L M A study appears therefore suitable to identify this system, providing that the P D dynamics are close or within the excitation bandwidth. Dynamics lower than the identification bandwidth may not be identified precisely due to the limited amount of excitation in that band. 6.2.5 Identification Results P D models were derived for each patient following the identification methods described previously. The results are summarized for each age group in Table 6.1. In the traditional P D approach, the models are fully defined by the three parameters {keo,ECso,-y} (for simplification, we assumed that EmaX=0 and Eo=100). As for the proposed P D approach, the models are defined by the parameter set {T</, kj,,EC^o,^} and the sensor dynamic H C N S ( S ) - In both cases, these models must be used in conjunction with Schuttler's pharmacokinetics. Note that, following the discussion in Section 6.2.1, the keo and kj, values are well contained within the identification bandwidth of the input signal. CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING T R A D I T I O N A L P D A P P R O A C H P R O P O S E D P D A P P R O A C H P A T I E N T # L T I fce0 E C 5 0 H i l l 7 Td L T I kd E C 5 0 H i l l 7 [ s -^ lO" 3 ] lM6/ml] M [ s - 1 -10- 3 ] G l : >18 - <30 years 007 10.0 1.8 2.9 22 133.5 3.2 4.7 008 13.8 2.7 3.2 4 44.4 3.1 2.5 010 0.8 0.4 3.3 44 25.0 2.4 1.9 015 1.9 1.1 1.8 45 51.5 3.8 1.2 016 4.2 1.7 1.9 39 85.7 3.8 2.3 023 10.7 3.1 2.8 18 82.5 3.9 2.1 030 5.6 1.8 3.1 32 44.4 2.9 2.8 035 8.4 1.9 3.8 12 26.7 1.9 2.3 038 11.8 3.0 2.6 7 35.2 3.4 1.9 046 9.9 2.4 3.9 9 32.8 2.8 2.8 048 8.6 2.2 3.3 17 46.4 2.8 2.3 053 9.9 2.1 3.4 4 26.2 2.4 2.5 058 ' 10.6 1.9 3.0 9 50.4 2.5 2.6 066 12.9 2.4 2.4 18 160.5 3.6 3.9 071 8.1 1.8 2.7 20 75.0 2.5 1.9 M e a n 8 . 5 2 . 0 2 . 9 2 0 . 0 6 1 . 4 3 .0 2 . 5 S D 3 .8 0 . 7 0 . 6 1 3 . 9 4 0 . 1 0 . 6 0 . 8 G2: >30 - <40 years 006 2.4 1.1 4.6 44 54.8 3.2 2.7 009 5.9 2.2 2.6 29 83.1 4.0 2.3 029 6.3 2.6 3.2 18 34.4 3.7 2.1 036 10.9 2.9 1.5 1 29.6 3.3 1.2 047 11.9 3.0 2.2 1 24.9 3.1 1.5 049 8.7 3.0 2.5 12 35.2 3.9 1.8 051 9.5 2.5 2.9 4 24.8 . 2.7 2.0 061 6.1 1.8 2.9 12 28.7 2.8 2.2 063 8.2 2.1 2.5 5 27.0 2.8 2.1 065 17.0 3.3 2.5 4 67.2 3.6 2.0 068 7.3 2.4 2.5 12 29.3 3.1 1.8 074 5.1 2.1 2.1 13 29.1 3.7 1.8 M e a n 8 . 3 2 . 4 2 . 7 1 2 . 9 3 9 . 0 3 . 3 2 . 0 S D 3 . 8 0 . 6 0 . 7 1 2 . 6 1 9 . 0 0 . 4 0 . 4 G3: >40 - <50 years 004 2.1 0.9 2.6 35 38.0 3.3 1.8 025 9.5 4.4 1.9 11 36.6 6.1 1.3 027 15.3 4.7 2.2 2 32.6 4.7 1.3 040 10.8 4.1 2.4 12 35.0 4.5 1.4 042 7.2 2.8 2.5 10 28.7 3.9 2.0 043 7.7 2.8 2.5 12 34.8 3.9 1.8 052 11.7 3.2 3.1 9 36.6 3.2 1.9 069 9.8 2.7 2.9 8 35.6 3.4 2.3 072 5.4 1.7 2.5 13 30.0 3.0 1.9 M e a n 8 . 8 3 . 0 2 . 5 1 2 . 4 3 4 . 2 4 . 0 1 .7 S D 3 .8 1.2 0 . 4 9 .1 3 .1 1.0 0 . 3 G4: >50 - <60 years 018 11.9 3.1 3.1 3 31.5 3.5 2.3 033 4.4 2.3 1.8 29 42.0 4.4 2.2 041 8.3 3.7 1.6 2 21.8 4.7 1.4 057 5.4 2.2 1.6 16 28.8 3.7 1.1 060 7.0 2.7 2.6 10 26.4 4.0 1.9 064 12.7 3.9 1.3 6 58.0 5.0 1.5 070 8.7 3.1 2.1 6 32.2 4.2 1.5 075 5.2 1.8 2.4 12 24.3 3.1 1.8 M e a n 8 . 0 2 . 9 2 .1 1 0 . 5 3 3 . 1 4 .1 1 .7 S D 3.1 0 . 7 0 . 6 8 . 8 1 1 . 8 0 . 7 0 . 4 Table 6.1: P D models obtained from the proposed approach with E m a x = 0 and EQ=100. CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 109 s + ke0 K E C 5 , .E, c- r , s 1 • 2- ECa, Figure 6.7: LTI part of the models, (a) In the traditional approach, the Hill equation is replaced by its static gain, (b) In the system-oriented approach, the sensor dynamics is to be explicitly included. 6.2.6 Model Validation Any identification procedure involving a priori model structures, such as the one proposed by Sheiner et al. or the one defined in (6.6), should be followed by a discussion on the adequacy of the model for representing the system. Otherwise, the model may not represents accurately the system when subjected to inputs different than the ones used for estimating the model parameters. Model validation can be carried out in a number of ways. It is usually recommended to test the model using a combination of tools. To give confidence in the models derived in this chapter, we present here validation results using some of these tools. D i r e c t R e s i d u a l A n a l y s i s In pharmacodynamic modeling, the analysis of the residuals should be carried out only on the linear part of the model. In the traditional approach, the LTI part of the model is limited to a first order element and the static gain K E C 5 0 °f the Hi l l dose vs. response relationship (see Section 6.3.4). In the proposed approach, the LTI part of the model includes both PD(s), as defined in (6.6), and the sensor dynamics, see Figure 6.7. To illustrate this discussion, we consider here the pharmacodynamic models obtained for Patient #65. In this particular case, no anesthetic gas was delivered to the patient after the insertion of the L M A due to a mechanical failure. As a result, the patient progressively regained consciousness as the propofol concentration decreased in the blood plasma. About 7 minutes into this case, arousal was observed and the situation was promptly corrected. In terms of identification data, this case provides us with both a fast and large amplitude transient following propofol administration, and a slower transient during the distribution and elimination of the propofol, see Figure 6.8. The model outputs were calculated following the parameters calculated in Table 6.1. The residuals are usually calculated as the difference between the measured and predicted outputs, see Figure 6.9.a. Using the proposed system-oriented model, it can be noted that the peak amplitude of the residuals during the fast and large input transient is similar to the amplitude obtained during the slower transient. Considering that the system noise characteristic remains identical over the whole time window, we can reasonably infer that most of the linear dynamics have been accounted for. Conversely, the residuals obtained from the traditional approach are substantially larger during the first input transient. Direct inspection of the residuals tend therefore to stress out the inadequacy of the traditional LTI model. C H A P T E R 6. A SYSTEM ORIENTED A P P R O A C H TO PHARMACODYNAMIC MODELING 110 Figure 6.8: Propofol plasma concentration (bottom plot) and W A V C N S time course (top plot) compared to the outputs obtained from the PD linear models (Patient #65). Fast and large ^ amplitude transient ^ Slow transient 30 20 3 1 0 3 o a I "10 -20 -30 3 Pt 2 1 1 a /-o 0 I -1 p 1 -2 -3 Systeir -orient* d apprc ach V • 7 A Tr aditiona. appro; ch («) / V / L / niinrif lf l i \ I 0 20 40 60 80 100 120 140 160 180 200 220 T I M E [sec] (b) Figure 6.9: Measured vs. predicted W A V C N S time courses for both models (Patient #65). (a) Standard residuals, (b) Prediction errors (assuming an autoregressive model). CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 111 P r e d i c t i o n E r r o r A n a l y s i s A better way of expressing this is to consider the prediction error ep. Considering that the system can be modelled by an autoregressive A R X model, we have: A ( q ) - W A V c m ( k - T s ) = B(q)-Cp(k-Ts-Td)+e(k-Ts), (6.7) where e is a white noise, A and B are polynomials of the delay operator q -1, and Ts is the sampling time. This model expresses the fact that the output is a weighted sum of past Cp, W A V C N S , and e values. The residuals e are calculated as: e(k • Ts) = W A V c m ( k • Ts) - • Cp(k • Ts - Td) - -^e(k • Ts), (6.8) In case the model perfectly describes the system, the residuals can be expressed as: e(k-Ts) = -^-ye(k-Ts), (6.9) Consequently, the prediction error e p defined as ep(k • Ts) = A(q) • e(k • Ts) should ideally have the char-acteristics of a white noise. If not, the model did n ° t capture all of the system dynamics. The prediction error obtained from the residuals calculated from both models are presented in Figure 6.9.b. As can be observed from this plot, both prediction errors behave like white noise during the slowest part of the input transient. However, during the initial large input transient, the prediction error obtained from the traditional model is clearly biased as it shows a large oscillation. W h i t e n e s s Test To further test the whiteness of the prediction error, it is usually recommended to calculate the autocorrelation of e p: ReJr) = E L i e P ( f c - r s ) - e p ( f e - T s - r ) £ L i e p ( * - r - ) 2 where N is the number of samples in ep. Note that (6.10) is here normalized (i.e., Rep(Q) = 1) by the autocorrelation coefficient at lag 0 (equivalent to the mean of e2). A n easy way of assessing whether ep is white is then to check whether the signal yields autocorrelation coefficients contained within the 99% confidence interval of a Gaussian distribution A/ r(0,1). The corresponding bounds for the normalized autocorrelation function can be calculated as: SD. = ^ 8 - (6.11) ep VN A large number of coefficients outside these bounds means that the residuals can be predicted by considering their past values. This is usually a sign that the model order is insufficient. We plotted the result of this analysis in Figure 6.10.a. It clearly appears that the model order in the traditional approach is not adapted to propofol pharmacodynamics. This result suggests that the order of the A polynomial in the traditional approach is insufficient. In that respect, the proposed system-oriented approach provides much improved results. CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 112 Independence Be tween P r e d i c t i o n E r r o r and Input Another useful test is to calculate the cross-correlation between the prediction error ep and the input signal Cp: Similarly to (6.10), the cross-correlation is normalized by considering both the autocorrelation of the pre-diction error at lag 0, and the autocorrelation of the input Cp at lag 0. Ideally, this should also have the characteristics of a white noise, otherwise not all of the input contributions to the output are accounted for by the model. Similarly, this test can be checked easily by using the 99% confidence intervals for the cross-correlation function: The results of this analysis are presented in Figure 6.10.b. It appears that both models have captured appropriately all the contributions of the input signal. This indicates that the order of the B polynomial in the traditional approach may already be sufficient. Note, however, that the delay in this particular case is rather small (4 seconds only). A traditional model derived on a case involving a larger delay (such as for patient #15) can fail this test. M o d e l V a l i d i t y Both the visual inspection of the residuals and the analysis of the prediction error tend to warn against the use of the traditional model for close loop control. A t best, traditional models can be used to predict the time course of the drug effect, but only for inputs similar to those used in the identification process. Conversely, the proposed system-oriented model structure, which incorporates both a first order plus delay and the sensor dynamics, yields adequate results, and accurately models the frequency response of the system. Note, however, that there is no guarantee that this model structure validated for the data collected on the patient #65 will also yield satisfactory results in other cases. Ideally, this analysis should be carried out for each individual cases. However, is is our experience that all cases presented in Table 6.1 are consistent with the example of patient #65. 6.2.7 N o m i n a l Propofo l P D Mode l s The models presented in Table 6.1 were obtained specifically for each patient. In order to expand on this result, we derive population-normed nominal models that can be used to predict the W A V C N S time course for any surgical patient. A usual approach would be to calculate the nominal model parameters by averaging the P D parameters over the whole study population. However, a particularly interesting result is the age dependency of some of the P D parameters. For instance, in Sheiner's approach, both the EC$o and the 7 parameters were shown to be statistically different between the 18-39 yrs and 40-60 yrs age groups (p < 0.05). Sfc=i  eP(k • Ts) • Cp(k - Ts — T) (6.12) (6.13) C H A P T E R 6. A SYSTEM ORIENTED A P P R O A C H TO PHARMACODYNAMIC MODELING 113 Fast and large amplitude transient < • Slow transient -0 § 0.5 S -0.5 y 1 System oriented model 99% covfidenci intervals • 1 1 i ; i 1 I 0.04 5 e> 2 U 0.02 a < 3 a o o z y H Fast and large amplitude transient < > Slow transient i 0.02 -0.04 Syst sm oi iente j mo del c« ink rvals V y 0 20 40 60 80 100 120 140 160 180 200 220 0 20 40 60 80 100 120 140 160 1 80 200 220 L A G [sec] L A G [sec] W (b) Figure 6.10: Standard residual analysis for model validation, (a) Whiteness test (note that the autocorrelation coefficient at lag 0 is not supposed to be contained within the confidence intervals). Note that the traditional PD model fails this test, indicating the need for a higher degree model, (b) Independence of the residuals with respect to the input signal Cp. CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 114 PD P A R A M E T E R S E T C O V A R I A T E S I G N I F I C A N C E Traditional PD feeO = 8.4 No covariate EC no ECX = 1.4 + 0.029 • age p < 0.05 between G1+G2 and G3+G4 7 7 = 3.5 - 0.023 • age p < 0.05 between G1+G2 and G3+G4 System oriented PD Td ka EC50 7 Td = 26 - 0.3 • age kd = 81 - 0.99 • age ECso = 1-94 + 0.043 • age 7 = 3.0 - 0.027 • age p < 0.05 between G l and G2+G3+G4 p < 0.05 between G1+G2 and G3+G4 p < 0.001 between G1+G2 and G3+G4 p < 0.005 between G1+G2 and G3+G4 Table 6.2: Population-normed PD models with age as a linear covariate. Similarly, all of the parameters obtained from the system-oriented approach have shown a similar de-pendence. In particular the EC$o parameter was shown to increase significantly with age, see F ig . 6.11 2. Conversely to age, neither weight nor gender could be found to be PD covariates. Using a basic linear fit, population-normed PD parameters were derived for both approaches with age as a covariate. These models are summarized in the Table 6.2. 6.3 Discussion In this section, we discuss the results of the results of the two identification procedures carried out to model the effect of propofol administration onto the W A V C N S -6.3.1 Clinical Adequacy of the Identification Data We already discussed in Section 6.2.4 the adequacy of the data in terms of their suitability for the identifi-cation procedure. In this section, we discuss their suitability for deriving relevant and meaningful clinical models for propofol pharmacodynamics. As compared to other published PD studies, it will be clear to anyone well versed in pharmacology that the identification data used here suffer from serious limitations. First, they are limited to a single bolus of propofol given within a very short period of time and in an uncontrolled manner (i.e., the rate of injection was assumed constant but may have varied). As a result, the Cp(t) time course may only be a rough estimate of the actual time course of the plasma propofol concentration (in a dedicated PD study, blood samples would have been obtained to check the adequacy of the PK model to capture the time course). Also, a small dose of a fast acting opioid was administered concomitantly, which might have interacted 2This result may be surprising. In Schnider et al. [166], the authors concluded that the propofol dose required to provoke unconsciousness decreases with age. As compared to Schnider's study, we are here calculating the average propofol dose which provokes a suitable anesthetic depth ( W A V C N S ) - We found that older patients actually require more propofol than younger patients, which is consistent with the results published by Kazama et al. CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 115 4* I*"1 ECS0 [ug/ml] -Ludbrooke/fl/.,2002 •Billard etal., 1997 + This study J (allgroups) le-3 2e3 3e-3 Kuizenga e/ a/., 2001 / (initial infusion) Kazama etal, 1999 — [allgroups) 4e3 -+-5e-5 H I 1 1 1 6e3 7e3 8e3 9c3 10e3 Kuizenga etal., 2001 J (sequential infusions) This study f (all groups) Kuizenga etal, 2001 (sequential infusions) Kuizenga et al., 2001 (initial infusion) i—Billard «/a/., 1997 — 1 H Kazama etal, 1999 (20-5?jro) J This study (a// groups) Kuizenga et al., 2001 / (initial infusion) ~ 1 Kazama etal., 1999 (allgroups) T J Billard 1997 Kuizenga etal., 2001 (sequential infusions) l i e 3 Figure 6.12: Comparison between the identified PD parameters from the traditional approach, and the parameters from literature. synergetically with the initial propofol bolus. It is therefore reasonable to question the reliability of the PD results presented in this chapter. In Fig. 6.12, we compare the PD parameters obtained from the traditional approach with other published parameter sets involving the bispectral index as an observed effect. The keo time constant is faster than most studies, but this could be due to the fact that W A V C N S has a faster dynamic response than BIS during induction. The E C 5 0 of our study was in the same range than Kuizenga et al. [118] and Billard et al. [117] studies. However, the values were significantly less than the values reported by Kazama et al. [116]. However, this can be explained by the difference between the PK parameter set used here, and the PK set used by Kazama. Finally, the steepness of the Hill characteristic is also reasonably similar to that of other studies. The stronger steepness reported by Kuizenga and Billard can be due to the switching characteristic of the BIS monitor as well as the inherent time delay of BIS (note that Kazama's study did not report such a large 7 value, probably because they compensated for the BIS and the PK delay). CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 116 In conclusion, even though the identification data were not obtained specifically for a P D analysis, the consistency of our results with other studies tends to indicate that they are appropriate. It is, however, reasonable to expect a better fit and less parameter variability when using data where propofol is the sole agent and where the rate of administration is controlled. 6.3.2 System Oriented vs. Traditional Approach To compare the proposed approach to the traditional approach, we define a number of performance indexes. These indexes can be classified into three categories: i . Performance Error (PE) based indexes For each patient i, PE{ is defined as the error between the measured effect and the predicted effect: PEi = {PEij = Wij - Wij}j=it2,...,Ni (6.14) where Wij is the jth measured W A V Q N S sample and w% is the corresponding predicted W A V C N S - Note that the total number of samples Ni is different for each patient. The mean value of the performance error represents the bias between measured and predicted effect: 1 N i Has(PEi) = — -J2PEij (6.15) 1 3=1 The accuracy of the model can be assessed by calculating the root mean square of the residuals: RMS(PEi) = K . 1 3=1 Note that the P D parameter sets in Table 6.1 have been derived by minimizing for each patient the R M S ( P £ i ) value. Some clinicians also use the Mean Absolute Residual ( M A R ) value as an indicator of model accuracy: 1 Ni MAR(PEi) = — -J2\PEij\ (6.17) 1 3=1 The M A R indicate the expected average error between the measured effect and the model output, i i . Correlation coefficient The correlation coefficient r? is defined for each individual i as: r- = 1 3 r~o (6.18) Ejii ( ^ • - ^ • £ , = 1 ^ 3 ) r? is a measure of the correlation between the measured and predicted W A V C N S time courses. Both time courses are identical when rf = 1. Note that this value can be negative. CHAPTERS. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 117 100 80 J 40 20 30 " 60 90 Time [s] Figure 6.13: Error in peak effect (time and value). In clinical studies, the median of rf over the study population is considered to be representative of the performance of a given model. This parameter is a meaningful performance index whenever the identification data cover a large range of plasma concentration and effect. This measure is thus particularly suitable in this analysis. i i i . Peak Effect The peak effect corresponds to the time tp and value wp of the maximum W A V C N S depth immediately following propofol administration and before airway management, see F ig . 6.13. Note that to reduce sensor noise, tp and wp correspond to the time and W A V C N S value when the index first comes within 2% of its maximum peak value. The time to peak effect has been recently used by Minto et al. [165] as a clinically relevant endpoint for assessing the performance of P K P D models. For each patients, we define two measurements: £i,t„ = (ti,P - U,P), and eiiWp = (witP - w"i;P) (6.19) Patient-specific models The performance of the P D models of Table 6.1 has been assessed using the performance indexes defined previously. The results of this analysis are summarized in the Table 6.3 for the traditional models and Table 6.4 for the system-oriented models. As compared to the traditional P D approach, the root mean square of the residuals are reduced by about 50%. Also, the overall accuracy of the fit is significantly improved with a median r 2 of 0.968 as compared to 0.879. Of particular interest, 4 cases in the traditional approach scored a correlation coefficient inferior to 0.7 which tends to indicate that the identification procedure failed in obtaining a proper model for these cases. Probably the most dramatic improvement is the reduction of the prediction error of the time to peak effect. The system oriented P D approach also seems to yield more consistent results as shown by the reduced standard deviation in all performance indexes. Population-normed models A similar performance evaluation was performed for the population-normed models presented in Table 6.5 and Table 6.6. C H A P T E R 6. A SYSTEM ORIENTED A P P R O A C H TO P H A R M A C O D Y N A M I C MODELING P E R F O R M A N C E E R R O R C O R R E L A T I O N P E A K E F F E C T bias R M S M A R CT(MAR) Median Best Worst N ei,tp W [%] °{F-i,u,p) [%] >18 - <30 yrs 0.9 10.6 8.7 6.0 0.859 0.964 0.689 -54.5 25.3 -4.1 5.3 >30 - <40 yrs 0.8 9.4 7.8 5.2 0.895 0.952 0.779 -46.7 16.2 -3.8 4.3 >40 - <50 yrs 0.5 8.8 7.3 4.8 0.889 0.934 0.838 -40.0 18.5 -5.5 4.6 >50 - <60 yrs 0.1 9.5 8.0 5.1 0.834 0.947 0.530 -43.8 15.7 -6.6 7.3 A l l groups 0.7 9.7 8.0 5.4 0.879 0.964 0.530 -47.4 20.3 -4.8 5.2 Table 6.3: Goodness of fit of patient-specific models obtained from the traditional PD approach. P E R F O R M A N C E E R R O R C O R R E L A T I O N P E A K E F F E C T bias R M S M A R a (MAR) Median Best Worst £i,tp M 'P M tp) [%] Wp CT(£i,»P) [%] >18 - <30 yrs -0.2 4.8 3.7 3.1 0.985 0.991 0.897 -10.1 8.3 -3.0 4.5 >30 - <40 yrs -0.3 5.1 4.2 3.0 0.974 0.995 0.833 -9.3 7.2 -3.8 2.5 >40 - <50 yrs -0.3 5 4.1 2.8 0.964 0.993 0.923 -7.8 7.8 -5.1 2.6 >50 - <60 yrs -0.7 6.1 5.1 3.4 0.910 0.993 0.799 -11.6 4.8 -4.5 4.2 A l l groups -0.4 5.2 4.2 3.0 0.968 0.995 0.799 -9.8 7.2 -3.9 3.6 Table 6.4: Goodness of fit of patient-specific models obtained from the system-oriented PD approach. C H A P T E R 6. A SYSTEM ORIENTED A P P R O A C H TO P H A R M A C O D Y N A M I C MODELING 1 P E R F O R M A N C E E R R O R C O R R E L A T I O N P E A K E F F E C T bias R M S M A R <r(MAR) Median Best Worst ei,tp W N Zi^p a - ( E I I U , P ) {%} [%] >18 - <30 yrs 3.2 14.8 11.5 9.2 0.727 0.951 -0.127 -47.9 13.0 -2.1 7.3 >30 - <40 yrs 3.7 12.9 10.2 7.7 0.797 0.911 0.511 -39.4 13.3 -1.9 5.8 >40 - <50 yrs 3.9 13.2 10.9 7.4 0.783 0.900 -0.188 -35.4 13.8 -3.2 10.7 >50 - <60 yrs 0.1 12.8 10.9 6.6 0.712 0.934 0.258 -37.5 9.9 -6.7 9.0 A l l groups 2.9 13.6 10.9 7.9 0.772 0.951 -0.188 -41.1 13.3 -3.1 8.0 Table 6.5: Goodness of fit of population-normed models obtained from the traditional PD approach. P E R F O R M A N C E E R R O R C O R R E L A T I O N P E A K E F F E C T bias R M S M A R CT(MAR) Median Best Worst £i,tp M tp *(£i,tp) Is] [%] Up CT(£i,l»p) [%1 >18 - <30 yrs 0.6 12.7 9.7 8.2 0.865 0.983 0.105 -6.0 13.9 -1.3 6.8 >30 - <40 yrs -0.5 10.4 8.2 6.3 0.891 0.952 0.538 -6.2 13.6 -2.6 5.1 >40 - <50 yrs -0.4 9.9 8.2 5.6 0.886 0.964 0.364 -7.2 13.7 -4.8 8.9 >50 - <60 yrs -4.0 11.1 9.3 5.9 0.756 0.989 0.310 -14.1 10.9 -9.2 8.8 A l l groups -0.7 11.2 8.9 6.7 0.867 0.989 0.105 -7.8 13.2 -3.8 7.6 Table 6.6: Goodness of fit of population-normed models obtained from the system-oriented PD approach. CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 120 The improvement of the new P D approach in terms of predictive performance is not as significant as with the patient-specific models. As compared to in the traditional approach, an overall improvement of just about 18% in the Performance Error indexes is obtained. And while the median correlation coefficient is also higher (0.867 instead of 0.772), 10 cases in the new P D modeling approach score a correlation coefficient inferior to 0.7. When reviewing these cases, it appeared than the main source for error is the under or over-estimation of the time delay. On a more positive note, the performances in terms of peak effect were maintained. 6.3.3 Dose vs. Response Characterist ics The static dose vs. response characteristics of propofol are presented for both approaches in F ig . 6.14. These characteristics should be equivalent since they were obtained based on the same identification data. However, they are clearly different. For instance, let us consider the propofol dose vs. response for a 20 years old patient. In the traditional approach, we find that a W A V C N S of 40 can be obtained with a plasma concentration of 2.2 /xg-ml - 1 . According to the system oriented dose vs. response, the same effect is obtained for a concentration of 3.3 /^g-ml - 1 . This is a significant difference considering that a C P j S S of 2.2 /xg-ml - 1 would yield a W A V C N S of 68 according to the new P D approach. This inconsistency between the two static models can be explained by considering that un-modelled dynamics in the traditional approach yield colored residuals, which, in turn, influence the derivation of the Hi l l parameters (see the example in Figure 6.3). This problem is particularly acute when identification data contain mostly transitory signals. Any study using the traditional approach and transitory identification data will suffer from similar limitations. A cautionary note relative to the accuracy of the dose vs. response relationship should be brought to the attention of the readers in this case. Wi th the traditional approach, a better identification method would be to use only steady state data, where the plasma concentration is fixed and the effect reaches steady state. This is, however, very difficult to obtain in a clinical setting. First, no other drug should be used concomitantly in order to limit potential synergism. Further, no surgical stress or excitation is allowed in order to maintain the effect in steady state. Finally, even with the use of a T C I pump to reach a steady state plasma concentration, each experiment would take considerable time in order to bypass the transient response. The only study achieving this is the one conducted by Kazama et al. in 1999. In this study, the authors induced patients using a T C I pump programmed to reach and maintain a fixed plasma concentration. Once steady state was achieved, an incremental concentration step was programmed. Even though the identification data did contain some transitory signals, the input and output signals contained large epochs of steady state data. When comparing Kazama's dose vs. response relationships to the ones presented in Fig . 6.14.b, a number of similarities become apparent. For instance, Kazama et al. were the only authors who established that propofol requirements increase with age, for effects compatible with an anesthetic depth suitable for surgeries. They also showed that this characteristic is inverted (i.e., older patients CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 121 100 90 SO 70 60 50 40 30 20 • 10 -0 J \ 18-29 years -f— —j 30-39 years •• i 40-49 years - • 50-60 years 0 1 2 3 4 5 6 7 8 9 10 STEADY S T A T E PROPOFOL PLASMA CONCENTRATION Cp [ug/ml] 100 p 90 • SO • 70 • CO g 60 • 1 50 • 40 -30 -20 -10 -0 -0 A . \ \ X . . 1 2 3 4 5 6 7 8 9 10 STEADY S T A T E PROPOFOL PLASMA CONCENTRATION Cp [(ig/ml] (b) Figure 6.14: Identified W A V C N S to, Propofol dose-response relationships, (a) Traditional approach (b) New ap-proach. Note that the thick lines represent averages over the corresponding age group. become more sensitive to propofol) when only a light effect is sought. Kazama et al. also reported low 7 values, from 2.3 to 2.0 , which is consistent with our findings (from 2.5 to 1.7). Conversely, all other studies have reported much higher values (from 4 .2 to 5.7) . Finally, they also reported larger EC50 values than most studies. However, this is mainly due to the pharmacokinetic set used to drive the T C I pump. This pharmacokinetic set (from Gepts et al. [Ill]) may have over-estimated the propofol plasma concentration since it was designed for infusions (it is possible that the high rate of infusion following a step change in the pump setpoint should have been considered a bolus instead of an infusion). 6 .3 .4 Frequency Response Characteristic This discussion would not be complete without presenting the frequency response of the models derived in Section 6.2 .5 and represented as block diagrams in Figure 6.15. First, we need to linearize the non-linear Hi l l static dose vs. response relationship. Let us denote by f(x) the Hi l l function. For both approaches, f{x) can be written as: x'< (6 .20) where X0 = 0 .5 (proposed approach) and X0 = EC^o (traditional approach). To linearize / ( x ) , we consider small input variations x around a given operating point x. Considering further that / is a function, and assuming small enough variations x, we can approximate f(x) as an afflne function: f(x) as f{x) + KT-x, where: —7-1 7 • x' (6 .21 ) (6 .22) CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 122 1 (s + k21)-(s + fc31) V, (s + it)-(s + a)-(s + p) • s + ktQ Pharmacokinetics y / Pharmacodynamics (a) W A V C N S Sensor 1 (s + k21)-(s + *31) c-r,s 1 V , (s + rt)-(s + a)-(s + P) 2-ECJO s + kd 1-e" ( 8 S + 1)2 - W A V C N S Pharmacokinetics Pharmacodynamics (b) Sensor Figure 6.15: P K P D block diagrams, (a) Traditional approach. Note that the Hill saturation is characterized by both the steepness coefficient y and the EC50 parameter, (b) System-oriented approach. Conversely to the traditional approach, the Hill saturation is defined uniquely by its steepness coefficient 7. 1 (s + fc21)-(s + fc31) kto V"., (s + it)-(s + a)-(s + P) • s + Ka 4-EC, W A V C N S Pharmacokinetics Pharmacodynamics (») 1 (s + fc21)-(s + fc31) C -T , s 1 fcJ V, (s + Jt)-(s + a)-(s + B) 2 - £ C 5 0 s + fcd Y 2 1 1-e"* • (8-s + 1)2 s • W A V C N S Pharmacokinetics Pharmacodynamics (b) Sensor Figure 6.16: Linearized P K P D block diagrams. Note that the scaling function was removed since it does not affect the dynamic response of the system, (a) Traditional approach, (b) System-oriented approach. The pharmacodynamic gain depends therefore directly on the operating mode x of the system. For very small (x —> 0) or very large (x —> 00) doses of propofol, the gain is close to 0, indicating that variations of the propofol dose will not translate into significant changes in the observed effect. Conversely, if the system already operates in a range suitable for general anesthesia (i.e., 4 0 < W A V C N S <60), the P D gain is close to its maximum. Considering a nominal operating point set at W A V C N S — 50, the system then operates at a steady state plasma concentration of EC50. As a result, the overall pharmacodynamic gain of the system is expressed as: K e c » = i~kr0 (6-23) Equation (6.23) holds true for both traditional and system oriented models. The linearized block diagrams can then be represented as in Fig. 6.16. CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 123 Figure 6.17: Frequency response of the PD models derived in Section 6.2.5. (a) Traditional approach (b) System-oriented approach (frequency response also models the sensor dynamics for consistency with the traditional approach). The frequency response of the pharmacodynamic models are presented in Figure 6.17. The sensor dynamics were added to the system-oriented models for consistency with the traditional models (these models inherently contain the sensor dynamics). In both cases, age only resulted in marked steady gain differences (this was already observed in the dose vs. response relationships). Otherwise, age did not result in any significant differences in terms of the frequency responses. To further compare the dynamics of the traditional vs. system-oriented models, we calculated the averaged frequency responses across all ages and plotted them in Figure 6.18. Beside providing a better temporal fit between measured and predicted time courses, the system oriented approach yields P D models with slightly larger bandwidth for control purposes. This comparison also reveals that the phase margin of the system-oriented approach model drops sharply due to the time delay. Using traditional P D models for designing the controller frequency response may therefore lead to unforeseen instability if the design is too aggressive. The phase characteristic of the system-oriented model clearly motivates control designers to exercise caution when tuning the cutoff frequency of the open loop system. 6.4 Summary Paradoxically, pharmacodynamics have been traditionally handled through static dose vs. response curves. However, with the advent of fast-onset short-acting drugs, a lag between effect and plasma concentration was observed and could not be accounted by a static characteristic. A model of the drug effect dynamics was required. CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 124 1000 | 100 o 10 o -55 -45 -8 a -90 ? £ -135 -180 0.001 0.01 0.1 1 F R E Q U E N C Y [rad • s1] Figure 6.18: Comparison of traditional vs. system-oriented frequency response models (all cases averaged). In an attempt to compensate for this lag, an addition to the mamillary pharmacokinetic compartmental model was proposed by Sheiner et al. It was shown that it was possible to collapse the effect vs. plasma concentration 'hysteresis' by introducing a new rate constant ke0. This actually resulted in the addition of a single order unity gain transfer function to the static dose vs. response relationship. This approach rapidly became the standard in P D studies, and has now been in use for the last 25 years. In this chapter, we have revisited pharmacodynamic modeling using an approach based on system identification know-how. In this approach, the plasma concentration is obtained using a well defined P K model and does not necessitate blood samples. Further, the pharmacodynamic LTI function - which captures the effect dynamic - is modified in order to account for the time delay that exists between the administration of the drug and the onset of effect. Finally, the sensor dynamics is now a distinct part of the model, which allows other sensing technologies to be used in conjunction with the P D models identified using the proposed approach. The system-oriented model yields an improved fit and more consistent performances than traditional P D models. Another advantage of the proposed approach is that the accuracy of the static dose vs. response relationship benefits significantly from the fact that all of the effect dynamics are now modelled through the LTI function PD(s) and the sensor dynamics. As a result, even short transitory identification data (e.g., induction data) can be used to identify the static dose vs. response relationship, thus simplifying the identification procedure in terms of clinical protocols and data acquisition (the traditional approach often requires volunteer-based studies). For instance, we show in Figure 6.19 a comparison between the static dose vs. response model derived in this study, and the static models derived for propofol and BIS (v.3.12) by Kazama et al. (note that we use here the fact that W A V C N S and BIS are roughly equivalent in steady state - this equivalence was already established in Figure 5.12). The similarity between the two relationships =. —1 - -11 4 f — P —! F-—-I ——1 | | | j Syste m-c nt a\ e d apprt jacl L F=i T 2 ditiona i i T I oach j i • rac Iti( in J /* ppr >a :h S Y t i orient i pproa :h CHAPTER 6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC MODELING 125 100 90 80 70 [ ^ 60 9 5/3 50 I-40 30 « 20 10 0 20-39 years Ok Kazama et al. 'This study (using Gepts etj al. PR parameter set) _j i i 0 1 2 3 4 5 6 7 8 9 10 STEADY STATE P R O P O F O L PLASMA CONCENTRATION C P [ng/ml] (a) 100 90 2 80 3 t> 70 5 ^ 60 S3 > 50 40 30 20 10 0 40-59 rears Kazama et al. i ! 7 ^ This study (using Gepts et 2X..PK'parameler.s'et) j... _i i i i i i ' i i • 0 1 2 3 4 5 6 7 8 9 10 STEADY STATE P R O P O F O L P L A S M A CONCENTRATION C p [ng/ml] (b) Figure 6 .19: Comparison of the static dose vs. response relationships obtained for propofol. In their study, Kazama et al. used Gepts PK parameter set (see Table B. l ) . This set was originally derived for infusion administration. As a result, the EC 5o found in their study are significantly different than the ones we originally found. However, this difference originates mostly from the choice of the P K set. In order to limit the influence of the P K set, we therefore re-processed all the cases using Gepts P K parameter set. We also added EQ as an identified parameter. We found that EQ is comprised between 0.02 (GI) and 0.06 (G4). When scaling the effect using the W A V C N S scale, the Eo values translate in a baseline W A V C N S of 98 (GI) to 94 (G4). These slight differences in the overall awake baseline values were found to be also a function of age. is remarkable, mostly when considering that Kazama's study involved for each patient an identification window of more than 9 0 minutes during which a TCI pump was used to target and maintain step-wise plasma concentrations (the surgery was not allowed to start during the recording of the identification data). In comparison, our method uses very limited data (2-4 minutes) which are readily available from any surgery involving limited amount of opioids and the placement of an LMA 3 . Another advantage of the proposed modeling approach concerns the reduction of system uncertainty. In particular, we will see in the next chapter (see Section 7 .2 .2) that the uncertainty related to the use of Sheiner's PKPD approach is higher than that of the proposed approach in the low frequency band. This implies limitation in terms of achievable performance when designing robust controllers based on quantified uncertainty. Probably the most compelling advantage of the proposed approach is that the model structure yields non-colored residuals, indicating that it effectively models all of the system's linear dynamics. Considering the relatively poor prediction performance of population-normed models to predict the effect time course, it is likely that PD models will best be used within a close loop control framework rather than an open loop framework. Having a precise description of the frequency response of the system will then be invaluable. 3 As compared to an endotracheal tube, the LMA limits the stimulation related to the airway management. Chapter 7 Managing PKPD Uncertainty Patient variability is probably the most challenging aspect in closing the loop in anesthesia. Quantifying this variability and expressing it as a system uncertainty is a necessity in order to prove stability when closing the loop. The goal of this chapter is to express patient variability as a quantified uncertainty structure suitable for use with control design tools. Quantifying uncertainty allows us to assess more precisely the importance of the different sources of uncertainty in P K P D models. This can ultimately help us decide which control strategy should be implemented to reduce uncertainty and provide improved performance. To carry out this task, a quantification and analysis method for system uncertainty is proposed and described in detail in Section 7.1. This method is based on the classical representation of uncertainty as a Nyquist disk whose radius encompasses the original uncertainty region defined by the frequency response plot of the system. This method, while conservative, appears to be the most suitable for this particular application. In Section 7.2, we apply the proposed method to the propofol P K P D models identified in the previous chapter. In particular, we show how serious patient variability can be in terms of close loop control. Finally, in Section 7.3, we propose a number of approaches to reduce the uncertainty to more manageable levels. Some of these approaches are rather simple and straightforward, while others require more efforts. Using the uncertainty quantification and analysis tool described in Section 7.1, we compare all methods with respect to their efficacy in reducing uncertainty vs. their cost in terms of design effort and practicality in a clinical setting. 7.1 Quantifying Uncertainty When considering system uncertainty, two approaches can be envisaged: - Parametric (structured) uncertainty: uncertainty can be defined by considering bounded real perturbations in some (or all) of the system parameters (gains and time constants). In this framework, it is usual to describe any given parameter T as T = TQ • (l + u- As), where ro is the nominal value of r , 126 CHAPTER 7. MANAGING PKPD UNCERTAINTY 127 and A s € l , | A s | < 1. The quantity u = (TMAX — T~MIN ) / (T~MAX + TMIN ) is usually less than 1, and expresses how much the parameter can drift from its nominal value. In this framework, it is assumed that each parameter can take any value in their uncertainty range. Linear fractional transformations are used to represent the system under a form more suitable for robustness analysis. - Unstructured uncertainty: in this framework, only the overall uncertainty in terms of the gain and phase of the system is defined. The unstructured framework is usually selected when the model structure itself is poorly denned, or when the uncertainty cannot be expressed as parametric uncer-tainty. Note that parametric uncertainty can also be lumped into an unstructured framework. This can lead to a more conservative design since the unstructured framework contains a large number of plants that do not belong to the original set defined by the parametric uncertainty. When considering the PK parameters defined in Table 3.1 and expressed as in (3.10), it is possible to define for each parameter of the hybrid set {Vi ,7r, a,P,k2\, k3i} the corresponding parametric uncertainty bounds. This can be done using the standard error values defined Table 3.2. For instance, the PK parameters of a 30 yrs old 70 kg individual can be defined as: 7T = 9.4 lO" 3 (1+0.40 A s ) Is"1] hi = 9.9 • 10" 4 • (1+0.51 Ag 2 1) [s"1 a = 4.8 10~4 (1+0.46 A?) [s"1] hi = 7.1 • 10" 5 • (1+0.59 Ag 3 1) [s-1 P = 4.0 10~5 (1 + 0.51 A?) [s-1] Vi = 9.3 • (1-r - 0.19 • A^1) [1] These values were calculated considering a two-times standard error. According to this parameter set, the maximum steady state gain1 of the PK model is 7.5 • 106 while the minimum gain is only 8.1 • 104. Even though we are considering a very limited and well defined adult patient population, note that this represents nearly a possible hundred-times difference in the steady state level of propofol plasma concentration. This result is clearly erroneous. Parametric uncertainty assumes, indeed, that each parameter of the model can take any value from the corresponding uncertainty range. This may lead to a combination of parameters that is not possible depending on how parametric uncertainty is originally defined. For instance, the hybrid PK parameters are all function of the same 6 coefficients defined in Table 3.2. Therefore, we cannot choose their values independently (i.e., it is not possible to have for the same patient the minimum k2i and ksi values, while having at the same time the maximum TT, a, (3 and V i values). This dependency between parameters makes here a significant difference in terms of uncertainty. For instance, when calculating all possible PK models according to the parametric uncertainty defined by Shiittler et al, we find that the ratio between the maximum and minimum steady state gains is only 4.6. In terms of the PD parameters presented in Table 6.1, parametric uncertainty may indeed be considered for kd and EC50 (a student t-test reveals that these two quantities are independent). However, unless the time delay dynamics are expressed as rational LTI transfer functions (e.g., using Pade approximation), the uncertainty related to the time delay cannot be directly expressed as parametric uncertainty. 'Steady state gain is (k2\ • fc3i)/(Vi • n • a • /?). CHAPTER 7. MANAGING PKPD UNCERTAINTY 128 G 0 G0(s) Figure 7.1: Unstructured uncertainty expressed in a multiplicative framework. When considering P K P D uncertainty, the unstructured framework appears therefore more suitable and straightforward. We present here the multiplicative unstructured system uncertainty approach, which we will be following throughout this chapter. This approach is discussed in details in Skogestad et al. [167]. 7.1.1 The Relative Uncertainty Framework Let us consider a system G(s) for which some uncertainty bounds exist. A common approach is to model the uncertainty using a structure such as the one of Fig . 7.1, and which can be expressed as: G(jw) = Go • G0[ju) • (1 + w(ju) • A„( jw)) , (7.2) where | | A u ( j w ) | | < 1, Go(jw) is the normalized nominal model of the system, and Go = \\G$(jw = 0)|| is the steady state gain of the system. The weight function w(jiv) quantifies the magnitude of the unstructured uncertainty, while Au(JUJ) expresses this uncertainty as a unity disk in the Nyquist plot. The multiplicative uncertainty gain w(ju) is also referred to as relative uncertainty. To better understand this concept, let us consider the example of Fig . 7.2.a-b. We assume that the uncertainty bounds of G(JCJ) are defined in the frequency domain and that a nominal model G$(jw) exists. If we map the frequency domain information into the complex Nyquist plane, the uncertainty at each frequency UJ can be represented as the section of a ring centered at the origin. The outer and inner radii of the ring are defined by the uncertainty magnitudes G M I N = G M I N / G O and G M A X = G M A X / G O , while the section is defined by the uncertainty angles <£>MIN and ipuAX-To quantify the uncertainty according to the framework of (7.2), we define the disk V centered on Go(jw) and of radius Z(w). This circle is the smallest circle centered on Gnfj'w) and which contains all of the original uncertainty surface. Following (7.2), the multiplicative weighting function W(JCJ) is defined as: 11 IIGoO)|| 7.1.2 Selection of a Near Optimal Nominal Model (7.3) It is important to note that the uncertainty quantified by (7.2) considers a larger uncertainty surface than the one originally defined from the frequency domain. A limitation of the method is therefore that its representation of the frequency domain uncertainty can be very conservative, mostly if the center of the disk Go(jw) is located close to the edges of the ring section. A control design based on the uncertainty defined by the disk T>(1(LJ),GO(JU)) may thus be unnecessarily conservative. CHAPTER 7. MANAGING PKPD UNCERTAINTY 1.29 Uncertainty bounds o f G(s) Unstructured uncertainty bounds <DfyH S0(jo>)) Uncertainty bounds of G(s) (b) Figure 7.2: Quantifying uncertainty, (a) Uncertainty expressed in the frequency domain via Bode plots, (b) Nyquist mapping of the uncertainty at the given frequency w. A simple way to reducing the uncertainty is then to minimize by selecting the best location for the center of the uncertainty disk. In other words, this leads to the optimization of the nominal model Go, see Figure 7.3. Note that this would still be a conservative characterization of the frequency domain uncertainty. A first step towards the optimization of the nominal model is to determine the optimal Nyquist path for Go(jw). We define the 4 following relevant coordinates: Ci G 3 = * i = GMAX • COS(<£>MIN) 2/I = GUAX • sin((^MiN) X3 = G M I N • C0S(<£MAX) 2/3 = G M I N • sin(<y5MAx) Co X 2 = GMAX • C0S(</?MAX) 2/2 = G M A X • s in^MAx) Xi = GMIN • COS(<^MIN) 2/4 = G M I N • sin(v?MiN) (7.4) The coordinates C\, C2, G3, and C4 mark the four corners of the ring. Depending on the phase uncertainty, we can distinguish between the three following cases: - (^MAX(W) — </?MIN(<*>)) < TT: in this case, it is possible to find a circle C that circumscribes the ring section (see Figure 7.3.a). The origin Ccir of C is located at the intersection of the perpendicular bisectors of the line segments {C\, C2} and {C2, G3} (or { G 4 , C i } ) . This can be analytically expressed as: Ccir b" - b' a — a" b" - V (7.5) CHAPTER 7. MANAGING PKPD UNCERTAINTY 130 Unstructured uncertainly bounds Optimized uncertainty bounds Figure 7.3: Effect of the nominal model Nyquist path with respect to the uncertainty radius: note that nominal models which are close to the uncertainty edges can significantly increase the unstructured uncertainty radius. where: x2 - x i 2/1 -2 /2 1 v\ ~vl-A 2/1 - 2/2 and X2 ~ X 3 J/3 " J / 2 1 y \ - y \ - x \ + x 2 J /3 -2 /2 (7.6) In most cases, the circumscribing circle minimizes the uncertainty radius l(ui). However, if the phase uncertainty covers a larger span than the gain uncertainty (i.e., the ring section is 'thin'), the center point of the line segment {C\, C 2 } may yield a smaller uncertainty circle (see Figure 7.3.b). In general, we find that the optimal center is defined as: O, opt c, opt ^opt — %cir Vopt — Vcir «4 * d r + Vcir < 2 ' \ / ( ^ l + X 2 ) 2 + ( j / l + J / 2 ) 2 Xopt = ( X l + X 2 ) / 2 Vopt = (j/1 + J/2)/2 (7.7) if: yf: Xlir + Vcir > o • V ( Z l + Z 2 ) 2 + (2/l+2/2) 2 ^ < (<PMAX(W) — (T ,MIN(W)) < 27r: when the phase uncertainty reaches 180 degrees, (7.5) yields the center of the circle that circumscribes the complimentary ring section (see Figure 7.3.c). The center of the minimizing circle is opposite to Ccir and located on the inner radius of the ring: opt Lopt Vopt \lxlir + yt 2 cir GMIN .2 'cir <jlr •2/c (7.8) CHAPTER 7. MANAGING PKPD UNCERTAINTY 131 Optimized unstructured Q uncertainty Circumscribing circle C w Figure 7.4: Relocation of the nominal model G 0 to the location G°PT to reduce the uncertainty radius, (a) (VMAX(W ) — V M I N ( ^ ) ) < T (large gain uncertainty relative to phase): in this case the circumscribing circle minimizes the uncertainty disc, (b) (<PMAX(^) — VMIN(W)) < 7r (small gain uncertainty relative to phase): in this case, the mini-mizing uncertainty disc center is located on the center-point of the {Ci , C2} segment, (c) (V>MAX(W) — <PMIN(V)) > TT: in this case, the minimizing uncertainty disc center is located on the inner radius of the ring section, directly opposite to the center of the circle which circumscribes the complementary ring section. - (lPuAx{i*>) — <^MIN(W)) > 2TT: in this trivial case, the ring section becomes a complete ring. As such, the center of the circle can be located anywhere on the inner edge of the ring. For convenience sake, we adopt the following convention: Copt — Xopt = GuiN • c Os ( (¥?MIN + VMIN)/2) yopt = G M I N • s in((<£MiN + VMIN)/2) (7.9) Note that this analysis is carried out for each frequency w of interest. This results in the definition of an optimal Nyquist path Copt(juj). Finding the Nyquist path describing the optimal nominal model is a first step. For control purposes, it is then necessary to express this path into a LTI transfer function which can be used in the design process. In most cases, it is not possible to find a rational LTI function that exactly describes the selected Nyquist path. Instead, a frequency domain identification can be carried out to determine a near optimal solution. The realization of the near optimal nominal model Gopt(ju) can be done in a number of ways: - Nominal parameters tuning: this tuning aims at minimizing the difference between the nominal Nyquist path and the optimal path by selecting the most appropriate model parameters. This yields a realization of the same order than the original nominal model. - Weight function: a LTI function wnom(ju)) is tuned in order for G°Qpt{JLo) = wncm(ju>)-Go(juj) to be as close as possible from the optimal Nyquist path. This approach yields a nominal model G°Gpt(juj) whose order is greater than that of the initial nominal model. CHAPTER 7. MANAGING PKPD UNCERTAINTY 132 - H y b r i d t un ing : a hybrid realization, where a combination of both methods are used, can also be investigated. However, this already requires much more efforts in terms of the identification procedure. Due to the large numbers of parameters in the initial P K P D nominal model, the tuning of a weight function wnom(jcj) appears to be a more attractive method to optimize the nominal Nyquist path. A general approach to identifying wnom(juj) is to minimize the difference between G°0pt(jui) and Copt(ju>) directly in the Nyquist complex plane by minimizing the cost function JNYQ '• J|cw*r»|| NYQ JO \wnom{juj)\\ -cosA</?(w) • du, (7.10) l | G o f » where A(p(u>) = arg(G 0(jw)) - &rg(Copt(juj)) + aig(wncm(ju)). The parameter A n o m ( w ) is a frequency dependent weight which can be used to improve the agreement between | |Go p t(jo;) | | and Copt(ju>) in the frequency band of interest (for instance, a good agreement in the controller bandwidth is more important than in the high frequency range). In practice, A n o m ( w ) = 1 in the lower frequency range and \nom(oj) = 0 for frequencies for which the system nominal phase reaches values below — 27r, or well beyond the expected controller cutoff frequency. In systems presenting a large time delay uncertainty, and whose nominal model is already of order >2, minimizing the Nyquist cost function JNYQ may be computationally demanding as the number of parameters of wnom(juj) may be large. A more effective method is to carry out the identification by first minimizing the gain difference between G°Qpt(juj) and C o p t ( j w ) by tuning only the poles and zeros of wnom(juj). This can be done by minimizing the cost function Jn0m'-2 nom f Jo Anom(^) *(^)U- • duj. (7.11) G 0 f » In order to obtain also a close agreement in terms of the phase, we add to wnom(juj) a time delay e ' ^ ' " calculated to minimize the phase difference. This simplified approach can only be done if the phase of e j T o p t u j s predominant to the phase of wnom(juj). 7.1.2.1 R e a l i z a t i o n o f the U n c e r t a i n t y Weigh t Once the near optimal nominal model is identified, the new uncertainty radius l(u>) can be easily calculated as: where: = \j(zoPtM - zi(w))2 + (j/opt(w) - yi(w))2 if: V M A X ( W ) - ^MIN(W) < 2?r l(ui) = G M I N ( W ) + G M A X ( W ) if: ^ M A X ( W ) - ^ M I N ( W ) > 2TT xrH = | | G f ( W ) | | - cx» (arg (GfH) ) (7.12) (7.13) [yfiw) = \\Gf{w)\\ .sin(arg(G?*H)) The resulting uncertainty weight magnitude || is directly calculated as the ratio between l(io) and the new nominal model GgP*(jw). CHAPTER 7. MANAGING PKPD UNCERTAINTY 133 Once Hw^ju;)!! is determined, an identification procedure is carried out to express w(juj) as a rational linear transfer function W(JUJ). This transfer function will be used in the design process to optimize the stability of the controller. The parameters of w(ju) are estimated by minimizing the criteria Jw: e(u>)2dw, where e(uj) — \ \w(jw) XW » 1, and: (7.14) (7.15) sign(e(w)) = 0 if e(w) < 0 [sign(e(w)) = 1 if e(w) > 0 Equation (7.14) represents essentially a least square minimization criterion. Note that the phase of w(jw) does not bear any relevance. The parameter Xw(UJ) simply adds a penalty weight on the realizations w(juj) whose magnitude at some frequency UJ is less than ||u>(ju;)||. Note that it is sometime necessary to have Xw(UJ) as a function of UJ when the uncertainty in the higher frequency range is large. In this case, we would choose Xw(UJ) and JW such that: V = JO 1 + A • sign(e(w)) e(u)2duj, where: < Xw(UJ) — 0 if UJ < 0.1 • uic (7.16) XW{UJ) » 1 if UJ G [0.1 • UJC; 10 • UJC] where UJC is the cutoff frequency of the controller. This choice essentially guarantees that the realization w(juj) captures adequately the system uncertainty within the controller bandwidth. Assuming a large phase margin in the low frequency region, the uncertainty is less of an issue for UJ < 0.1 -UJC. Finally, high frequency uncertainty is also not an issue since the cutoff frequency UJC of the controller guarantees a low gain in this frequency range. 7.1.2.2 Comparing Uncertainty Between Systems Note that any system presenting a relative uncertainty gain above 1 cannot be tightly controlled since the uncertainty disk contains the Nyquist origin. Therefore, there exists a particular plant for which the open loop gain can be 0, meaning that no control action can have any effect on the output. More generally, it can be shown that Robust Stability (RS) and Robust Performance (RP) are guaranteed if: (7.17) (RS)^\\w(juj)-T(juj)\\<l (RP) \\vjp(juj) • S(JUJ)\\ + \\w(juj) • T(jw)\\ < 1 where WP(JUJ) is the performance weight, and S and T are the sensitivity and complimentary sensitivity functions of the close loop system. Clearly, low ||«;(jw)| | values will benefit to the overall performance of • the system. Hence, we can compare the effect of system uncertainty between two systems SI and S2 by comparing directly their relative uncertainty weight: a system SI presents less uncertainty than a system S2 in a frequency bandwidth B W if ||wsi(iw)|| < ||tt*52(jw)||, for w G B W . CHAPTER 7. MANAGING PKPD UNCERTAINTY 134 7.2 Application to Propofol PKPD Models In this section, we calculate the relative uncertainty of the propofol P K P D models derived in the previous chapter. As a first step, it is necessary to clearly identify the sources of P K P D uncertainty to adequately account for them in our calculations and analysis. 7.2.1 Origins of P K P D Uncertainty In P K P D modeling, it is customary to distinguish between two different types of uncertainty: the uncertainty caused by inter-patient variability (i.e., the variability observed between different individuals), and the uncertainty originating from intra-patient variability (i.e., the variability observed within one particular individual). 7.2.1.1 Inter-patient Variability Inter-patient variability was already discussed in Chapter 3. It affects both the pharmacokinetics and pharmacodynamics of any given drug. For instance, it has been shown that age as well as weight, lean body mass, ethnicity, etc., are all factors of P K P D variability in humans. Co-existing illnesses involving either the liver and/or kidneys may also significantly alter the way drugs are metabolized and eliminated from the body. In general, 2 patients with similar physiological characteristics (age, weight, lean body mass, A S A ) may have largely different P K P D parameters. For instance, patient #15 in Table 6.1 (female, 21 yrs old, 53 kg, 157 cm, A S A I) and patient #53 (female, 21 yrs old, 67 kg, 163 cm, A S A I) have significantly different P K time delay (45 sec vs. 4 sec), EC50 parameter (3.8 ug/ml vs. 2.3 ug/m\), and saturation characteristics (Hill steepness of 1.2 vs. 2.5). Inter-patient variability can be easily characterized by considering the differences between P K P D models obtained over a large population of patients. In particular, Table 6.1 provides a good representative sample of an adult population with respect to the response to propofol administration. The parametric variability observed between the different P K P D models presented in this table is summarized in Table 7.1. Inter-patient variability clearly plays a prominent role in the overall system uncertainty. For instance, there is a significant difference in the P K time delay and P D time constant between patients. Also, while the EC50 variability is more limited, there is still a 6-times difference in terms of the overall P K P D steady state gain 2 . 7.2.1.2 Intra-patient Variability Intra-patient variability expresses the variability observed in the drug response within one particular subject. This variability originates from different factors. 2The overall PKPD gain regroups both PK and PD steady state gains, as well as the linearized Hill gain (operating point: W A V C N S = 5 0 ) . CHAPTER 7. MANAGING PKPD UNCERTAINTY 135 Table 7.1: Propofol PKPD Inter-patient Variability E C 5 0 7 M I N M A X Td o u M W [8]" [%] M I N M A X K D 0 u [s-1] Is"1] Is"'1] [%] M I N M A X E C 5 0 , o u [Mg/ml] [MS/ml] [/.g/ml] (%] M I N M A X 70 u [1] 11] [1] [%] A G E G R O U P G l G 2 G 3 G 4 4 45 24.5 83.7 1 44 22.5 95.6 2 35 18.5 89.1 2 29 15.5 87.1 25.0 160.5 92.7 73.0 24.8 83.1 53.9 54.0 28.7 38.0 33.3 13.9 21.8 58.0 39.9 45.4 1.9 3.8 2.8 33.3 2.8 4.0 3.4 17.6 3.0 6.0 4.5 33.3 3.1 5.1 4.1 24.4 1.9 4.7 3.3 42.4 1.2 2.7 1.9 38.5 1.3 2.3 1.8 27.8 1.1 2.3 1.7 35.3 P O P U L A T I O N 1 45 23 95.6 21.8 160.5 91.1 76.1 1.9 6.0 4.0 51.9 1.1 4.7 2.9 62.0 D r u g a d m i n i s t r a t i o n It is a well-documented fact that the pharmacokinetics of intravenous agents differ depending on the method of administration of the drug. Even though bolus and infusion PK models have the same steady state gain, the initial peak plasma concentration following a bolus administration is significantly over-predicted by the corresponding infusion model (see for instance Figure 3.4). During steady state (and for small setpoint changes and/or disturbances), it is likely that the controller will administer propofol at an infusion rate inferior to 0.5 mgmin _ 1 kg _ 1 . In this range, it is expected that the propofol pharmacokinetics will be accurately described by the infusion model. However, during large transients, the controller may have to output infusion rates above 1 mg-min_1-kg_1, in which case the propofol uptake and distribution may follow the behavior observed for bolus regimen. If the controller output is not constrained to infusion rates i