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 BRITISH 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 i n 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 closedloop 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 transitions, 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  iii  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. thorough introduction to clinical anesthesia.  Part A contains the first 4 chapters and presents a  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  ii  Table of Contents  iv  List of Tables  ix  List of Figures  xi  List of Acronyms  xx  Acknowledgements  xxii  Dedication  xxiii  P A R T A : The Anesthesia Framework  1  1  Clinical Anesthesia: Terminology, Concepts, and Issues  2  1.1  Clinical Anesthesia: A K e y Specialty in Critical Care  2  1.1.1  The Early Days  3  1.1.2  Risks and Outcome in Anesthesia  3  1.2  1.3  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  A n Extensive Pharmacopoeia  8  1.3.1  Anesthetics (i.e. hypnotics)  8  1.3.2  Opioids  1.3.3  Neuromuscular Blocking Agents (NMBs)  10  1.3.4  Drugs: Action, Effect and Interaction  11  •  10  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 2  AND THESIS ORGANIZATION  v  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  2.4 3  29  2.3.1  Using Physiological Measures  29  2.3.2  Using Heart Rate Variability  30  Summary  31  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  3.3 4  Quantifying the Depth of Analgesia  Pharmacodynamics  42  3.2.1  The Dose/Response Relationship: a Steady-State Model  42  3.2.2  The k o Parameter: a Measure of the Effect Dynamics  45  3.2.3  Literature Review  46  3.2.4  Concluding Remarks  48  e  Summary  51  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 i n Phase and Wavelets  62  5.2  Estimating the Anesthetic Drug Effect: the  WAVCNS  62  5.2.1  Concepts and Derivation  62  5.2.2  Practical Issues and Implementation  73  5.2.3  Clinical Results  76  ABSTRACT 5.3  AND THESIS  ORGANIZATION  vi  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  5.3.3  Implementation Issues  85  5.3.4  Case Reports  86  W A V A N S  84  5.4  Dynamic Behavior  91  5.5  Summary  94  6 A System Oriented Approach to Pharmacodynamic modeling 6.1  6.2  6.3  6.4  A New System Oriented Approach to P D Modeling  95  6.1.1  Sheiner's Approach  96  6.1.2  Proposed Methodology  Application to Propofol and  W A V C N S  :  101  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  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  Summary  123  7 Managing P K P D Uncertainty 7.1  7.2  7.3  7.4  95  126  Quantifying Uncertainty  126  7.1.1  The Relative Uncertainty Framework  128  7.1.2  Selection of a Near Optimal Nominal Model  128  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  Reducing Uncertainty  141  7.3.1  Passive Methods  141  7.3.2  Total and Partial Self Tuning  Summary  ,  144 149  ABSTRACT 8  ORGANIZATION  vii  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  8.3  8.4  8.5  8.6 9  AND THESIS  P I D Loop Shaping Design  164  8.2.1  First Design  164  8.2.2  Second Design  166  8.2.3  T h i r d Design  166  Hoo Control Design  168  8.3.1  Hoo Design Principle  168  8.3.2  Implementation  170  8.3.3  Application to Anesthesia Control  170  Reducing the Uncertainty  172  8.4.1  Accounting for age  174  8.4.2  Identifying the P K time delay  174  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  Summary  184  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  Bibliography  192  197  ABSTRACT  A N D THESIS  ORGANIZATION  A Pharmacopoeia  B  viii 209  A.l  Vapour Anesthetics  209  A.2  Intravenous Anesthetics  209  A . 3 Opioids  210  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  D.2  D.3 E  Arthroscopy Observational Study  229  D.l.l  Protocol  229  D.l.2  Results  230  D.l.3  Comments and Unexpected Results  230  Laryngeal Mask Airway Study  230  D.2.1 Protocol  232  D.2.2  235  LMA  Results  Electro-Convulsive Shock Therapy Study  Survey  236 242  List of Tables 1.1  Contemporary anesthetic mortality rates (adapted from Brown [5])  3.1  Propofol P K parameter sets from [101] where B W stands for body weight, ven= 1 for venous sampling,  5  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  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 P D 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 P D approach  119  6.6  Goodness of fit of population-normed models obtained from the system-oriented P D 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  ix  m a 3 :  = 0 and Eo=100  108  180  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 P K 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 D. l  228  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  D.5  E C T Study Population  238  D.6  Thiopental pharmacokinetic models  240  D.7  Thiopental identified pharmacodynamic models  241  235  List of Figures 2.1  Therapeutic window of the E E G for measuring the opioid effect (from [31])  2.2  Simultaneous records from the left and right motor regions, illustrating the changes in frequency  18  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]). . .  3.2  Time course of plasma concentration with two distinct phases: distribution and elimination (adapted from [98])  3.3  35  2-compartment pharmacokinetic model. A third compartment is often added for improved accuracy (adapted from [99])  3.4  33  36  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])  3.7  44  Different dose/response characteristics (a) On/Off effect (b) The observed endpoint cannot characterize 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 pharmacodynamic 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  5.3  63  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  5.4  64  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  5.5  65  Normalized E E G signals at different anesthetic depths. Each 1-second epoch is first detrended and then normalized using the R M S 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).  5.6  69  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  5.7  Awake and isoelectric reference PDFs based on / ( y  71 AV  71  LIST OF 5.8  xiii  FIGURES  Characteristic of the optimal function characteristic, (b)  WAVQNS  ^YVAV NS C  ^  o r t  ^ sampling frequency of 256 S/s. (a) Linearity e  values for the data sets R ,Ti,T^, a  T3  and R^. Note the larger variability  in the intermediate data sets 5.9  72  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 5.11 Time course of the BIS and  WAV NS C  n d  order low pass filter (u>o = 0.02 Hz and C=l).  for Patient #22 and Patient #8  5.12 Correlation during steady-state operation between the BIS (v.3.4) and of steady-state data were collected from the arthroscopy study,  77 WAVCNS-  A total of 13 hours  (a) Correlation density plot (b)  Bland-Altman test  78  5.13 Time courses of the BIS (v.3.4) and  WAVCNS  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) 5.14 Time delay between the BIS and  76  WAVCNS-  5.15 Time course of the BIS (v.3.4) and  79  (a) Induction (both study groups), (b) Emergence. . . .  WAVCNS  80  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: R : relaxed subject (y-scale: a  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  WAVANS  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  WAVCNS  (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  89  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  90  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  92  WAVCNS  5.22 Bode plots of the sensor dynamic (analytic and identified models) 5.23  Accuracy of  HQNS  and  .HANS  various step changes, (a) The HANS  for predicting the HQNS  WAVCNS  and  92  WAV NS A  is able to predict accurately the  time courses calculated for  WAVCNS  time course, (b) The  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.)  93  6.1  Traditional P D model block diagram  96  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 k  e0  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. 6.3  .  99  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  6.6  105  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 6.7  xv  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.8  Propofol plasma concentration (bottom plot) and  WAVCNS  109  time course (top plot) compared to the  outputs obtained from the PD linear models (Patient #65) 6.9  Measured vs. predicted  WAVCNS  110  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 P D model fails this test, indicating the need for a higher degree model, (b) Independence of the residuals with respect to the input signal C  113  p  6.11 Effect of age on the  ECBO  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) 6.14 Identified  WAVCNS  VS.  117  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 P D models derived in Section 6.2.5. (a) Traditional approach (b) Systemoriented 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 WAVCNS  7.1  WAVCNS  scale, the Eo values translate in a baseline  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  Unstructured uncertainty expressed in a multiplicative framework  128  LIST OF FIGURES 7.2  xvi  Quantifying uncertainty, (a) Uncertainty expressed in the frequency domain via Bode plots, (b) Nyquist mapping of the uncertainty at the given frequency u>  7.3  129  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  7.4  130  Relocation of the nominal model Go to the location (<^MAX(W)  — <£>MIN(W)) <  7r  to reduce the uncertainty radius,  (a)  (large gain uncertainty relative to phase): in this case the circumscribing  circle minimizes the uncertainty disc, (b)  (^MAX(W)  —  ^MIN(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)  (^MAX(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)  WAVCNS  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 7.7  137  Normalized Propofol P K P D uncertainty bounds expressed in the frequency domain. The P K part of the models was derived based on the published P K parameter sets of Schiittler et al. (a) Based on Traditional P D 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 7.9  P  '(JU;)  of the nominal model  140  Nominal model selection in the frequency domain. Note that the rapid gain drop in the high frequency range could not be captured adequately by w (ju). nom  Also, some small phase differences between  the optimal model and the realization G Q ' at about w = 1 • 10~ rad-s P  3  -1  results in 'bumps' which  significantly increase the uncertainty gain 7.10 Effect of age consideration on P K P D relative uncertainty  141 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 7.12 Effect of operating  WAVCNS  range reduction on uncertainty  143 144  LIST OF  xvii  FIGURES  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 ; 2-10 ] -4  rad-s  -1  -3  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  WAVCNS  range)  152  7.17 The frequency uncertainty bounds are represented as ring sections in the Nyquist plot. In the classical 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  WAVCNS  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 8.11 Measured vs. predicted  WAVCNS  161  time courses (LMA patient #065). The noise filter T^  oiee  switched from 90 s to 170 s. The noise filter T ^  oise  was  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 WAVCNS  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  WAVCNS  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  WAVCNS  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 GQ (jui). The second design is pt  unstable since there exists a frequency range where the complementary sensitivity is greater than the inverse of the uncertainty weight 8.17 Infusion rates and  WAVCNS  167  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  WAVCNS  unit  8.21 Close loop sensitivity analysis  171 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  WAVCNS  time courses in the 3 test patients. Accounting for age does result in a  significant increase in performance 8.25 Infusion rates and  WAVCNS  175  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 8.28 Infusion rates and  WAVCNS  181  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 of the  D.3  WAVCNS  WAVCNS  and  AWAVCNS  for both groups.  AWAVCNS  is the rate of change  index. Note the significant difference between the two groups before the L M A insertion. 236  Correlation between  WAVCNS  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 C is estip  mated through the 3-compartment P K model proposed by Stanski et al. This model is weight and age dependent, (a) Successful identification, (b) Successful identification despite limited fasciculation. (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 ACRONYM  DESCRIPTION  ACL  Anterior Cruciate Ligament  ASA  American Society of Anesthesiologists  ASSR  Auditory Steady State Response  BIS  Bispectral Index Scale  BSR  Burst Suppression Ratio  CNS  Central Nervous System  DOA  Depth of Anesthesia  DWT  Discrete Wavelet Transform  ECG  Electrocardiogram  ECT  Electro-Convulsive Therapy  EEG  Electroencephalogram  EMG  Electromyography  EOG  Electrooculogram  ESU  Electro-surgical unit  ERP  Evoked Response Potential  FDA  Food and Drug Administration  FIR  Finite Impulse Response  HRV  Heart Rate Variability  ICU  Intensive Care Unit  IIR  Infinite Impulse Response  LMA  Laryngeal Mask Airway  LOC  Loss Of Count  MAC  Minimum Alveolar Concentration  X X  LIST  OF  ACRONYMS ACRONYM  DESCRIPTION  MEF  Median Frequency  MLAER  Mid-Latency Auditory Evoked Response  NMB  Neuromuscular Blocking Agent  OA  Ocular Artifact  OAs  Ocular Artifacts  OR  Operating Room  PD  Pharmacodynamics  PDF  Probability Density Function  PET  Positron Emission Tomography  PID  Proportional-Integrative-Derivative control law  PK  Pharmacokinetics  PKPD  Pharmacokinetic - Pharmacodynamic  PS  Patient Simulator  REM  Rapid Eye Movement  RMS  Root Mean Square  RSA  Respiratory Sinus Arrhythmia  RWT  Redundant Wavelet Transform  sec  Semilinear Canonical Correlation  SEF  Spectral Edge Frequency  SEP  Somatosensory Evoked Potential  SP  Smith Predictor  STFT  Short Time Fourier Transform  SWT  Stationary Wavelet Transform  UBC  The University of British Columbia  UBCH  The University of British Columbia Hospital  VGH  Vancouver General Hospital  WAV  Wavelet-based Anesthesia Value Wavelet-based  WAVCNS  W A V  WT  Anesthesia Value for Central  Nervous System monitoring Wavelet-based Anesthesia Value for Autonomic  C  N  S  Nervous System monitoring Wavelet Transform  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 M s . Tatjana Zikov, who shared with me her scientific expertise i n 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 D r . Kris Lundqvist, D r . Henrik Huttunen, Dr. Jason Waechter, and M s . 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 particular D r . Mohammad Modarreszadeh and M r . Robert Schmidt, for their interest i n carrying forward the development of the N e u r o S E N S E ™ 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  a truly  righteous  to the memory and altruistic  xxiii  of my grandfather, man,  Albert  who will always  Vande  live in our  Kerckhove, memories.  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 i n 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 i n 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 i n 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 1 9  th  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 i n 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.1.1  1. CLINICAL  ANESTHESIA:  TERMINOLOGY,  CONCEPTS,  AND  ISSUES  3  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 William 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 i n 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 1 9  th  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. B u t 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)  1.1.2  was founded.  Risks and Outcome in Anesthesia  1.1.2.1  Mortality rate  Development of clinical anesthesia since the 1 9  th  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. B u t 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 1 9  th  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. ii. 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. iii. 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 t c . . . ) .  CHAPTER 1. CLINICAL ANESTHESIA: TERMINOLOGY, CONCEPTS, AND ISSUES STUDY  DATE  TOTAL  MORTALITY  CASES  RATE  1:1,560  Beecher Clifton  1948-52 1952-62  599,548 205,640  Harrison  240,483  Hatton  1967-76 1977  1:3,955 1:4,537  Lunn  1979  190,389 1,147,362  1:2,885 1:6,789  Eichhorn Eichhorn CEPOD CEPOD  1976-85 1985-88 1986 1986  757,000 244,000 486,000  1:151,400 0 1:185,056  486,000  1:1,185  5  a  a  b  a  e  d  d  e  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 ( A S A Physical Status I)  or i n 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  I n t r a o p e r a t i v e awareness  While patient's safety is still an important issue, new concerns arose from the use of neuromuscular blocking agents (NMBs) i n 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 intraoperative 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 i n defining and understanding the proposed research project. 1.2.1  T h e R o l e 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 i n 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 i n 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 "druginduced unconsciousness, [where] the patient neither perceives nor recalls noxious stimuli" [16]. This functional 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 cardiorespiratory 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 i n 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 contributing i n 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 i n the sense that they do not provoke unconsciousness. 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 technique 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 i n the management and conduct of clinical anesthesia.  1.3  A n 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 ing  1. CLINICAL  ANESTHESIA:  TERMINOLOGY,  CONCEPTS,  AND  ISSUES  9  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 sevoflurane. All 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 Concentration ( 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 i n 50% of patients. The M A C  value is thus  used to compare dosage between inhaled anesthetics. Intravenous Anesthetics tal),  Intravenous anesthetics can be classified into 5 families: Barbiturates (thiopen-  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 i n 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 i n 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. A s 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 i n 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.3.2  1. CLINICAL  ANESTHESIA:  TERMINOLOGY,  CONCEPTS,  AND  ISSUES  10  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 ( N M B s )  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 i n 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.3.4 The  1. CLINICAL  ANESTHESIA:  TERMINOLOGY,  CONCEPTS,  AND  ISSUES  11  Drugs: Action, Effect and Interaction 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, intravenous 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 i n 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. A s 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.4  1. CLINICAL  ANESTHESIA:  TERMINOLOGY,  CONCEPTS,  AND  ISSUES  12  Conduct of Anesthesia  In most cases, we can divide the anesthesia procedure into 3 distinct phases: induction, maintenance and emergence. Induction  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. A s 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. An  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. A s 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,  Maintenance  CONCEPTS, AND ISSUES  13  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. A s 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 anesthesiologists [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 I V A 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 concentration 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 Delivery Devices  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). The Anesthesia Machine  The Anesthesia Machine delivers gaseous drugs and vapors (air, N 0 ) and 2  has the ability to control ventilation. Its role is essential to ensure patient safety i n 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 0 flow meter, a vaporizer, a 2  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 . by  Pressure, volume, gas and vapour composition can be monitored and manually controlled  the anesthesiologist.  Sophisticated breathing circuits can maintain temperature, water content, and  usually incorporate a C 0 absorption cannister to remove the excess carbon dioxide in close circuit systems 2  (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 u s i n g o p i o i d s : large doses of opioids can significantly depress the respiration cycle. The brainstem mechanism regulating the blood C 0 concentration becomes less sensitive, thus inducing 2  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  - when using N M B s :  ANESTHESIA:  TERMINOLOGY,  CONCEPTS,  AND  ISSUES  15  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 adequacy 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 i n 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 anesthetic 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 O R , 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 or in ng • m i n 1  - 1  information, the pump automatically controls the rate of descent of the plunger.  • k g . Based on this 1  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  Anesthesia Monitors  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 Q R S 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 i n 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 unconsciousness 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, anesthesiologists 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. Until the m i d - 2 0 century, only inhaled anesthetics were available to create the state of general anesth  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 elimination. 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 cardiovascular 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. A s 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.1  2. QUANTIFYING  DEPTH  OF ANESTHESIA:  A  REVIEW  18  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 i n 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, Billard and Shafer [31] observed that the E E G could represent adequately the opioid effect, but only for a limited range of doses, see F i g . 2.1. Small doses do not have a marked influence on the E E G , surgery reach a maximum E E G  effect. Therefore, the use of the E E G  while larger doses used for cardiac seems to be irrelevant in preoperative  and postoperative I C U (where often only a shallow analgesic state is provided) and i n surgeries where cardiac stability is ensured by large opioid amounts.  E F F E C T SITE ALFENTANIL  Figure 2.1: Therapeutic window of the EEG  [ng/ml]  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 ( D O A ) . 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,  ii.  show similar changes for different agents of the same family (anesthetics, opioids),  iii.  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  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  19  Quantifying the D e p t h of Hypnosis: 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 ( E E G ) - 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 neurophysiologists 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  EEG  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 But  EEG.  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 of Tools a n d Techniques  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 integrating 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  Left motor region  DEPTH OF ANESTHESIA:  A REVIEW  20  Breathing cyclopropane  Right motor region  (a) Breathing air  (c)  (d)  (e)  (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) 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]).  CHAPTER  2. QUANTIFYING  DEPTH  OF ANESTHESIA:  FREQUENCY [HZ]  A REVIEW  21  FREQUENCY [HZ]  Figure 2.3: Awake and anesthetized EEG  and their power spectrum,  signal directly as a feedback quantity. 2.2.1.2  Power Spectrum Analysis  W i t h the development of microprocessors and signal processing tools, researchers have focused their attention on Fourier analysis of the E E G . the E E G .  Power spectrum analysis is used to obtain a frequency distribution of  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. W i t h 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. To  quantify the effect of anesthetics on the E E G ,  2.3.  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 F r e q u e n c y ( 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. Spectral Edge Frequency ( S E F )  The S E F 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. A s 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 DOA 2.2.1.3 The  based on the E E G power spectrum.  E E G modeling  idea of modeling the E E G using auto-regressive techniques dates back to the early 1970s. B u t its  application to anesthesia is more recent ([49]). As compared to univariate descriptors, auto-regressive (AR) can  modeling generates a set of parameters that  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 i n dogs using halothane. However, the size of the proposed network is rather large  (43 inputs, 43 nodes i n the primary hidden layer, 6 nodes i n 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 i n 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 F i g . 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 i n 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 (REM)  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. [56].  [55] with isoflurane and later on with propofol  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] . This index, referred to as BIS (Bispectral Index Scale), is a weighted sum of each spectral 1  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 1  Note once again the confusion between hypnosis and analgesia.  CHAPTER  2.  QUANTIFYING  DEPTH  OF ANESTHESIA:  A  REVIEW  Awake  • Memory Intact  24  BIS  Light Hypnotic ^ State Moderate Hypnotic State  Light/Moderate Sedation r 1  Deep Sedation  Low Probability of Explicit Recall  General Anesthesia •  Low Probability of Consciousness  Deep Hypnosis Deep Hypnotic < [20}State  • Memory Function Lest  Increasing Burst Suppression  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. W i t h increasing concentration of anesthetics, the index decreases. General anesthesia is obtained for an index between 60 and 40, see F i g . 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 oscillatory signals within the E E G can  itself. It has been advocated that such transient signals, if properly analyzed,  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 F i g . 2.6. The most remarkable feature beside the change i n amplitude of the signal is the change i n latency of the P and P;, waves. A very interesting work by Huang et al. [64] has shown that a  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 i n 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  2 +  DEPTH  OF ANESTHESIA:  A 4  Awake: Cp—3 [ig/ml  2  4-  REVIEW Sleep: 0^,-5.5 [ig/ml  ms  -14-  tr  25  ms  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 i n this new market segment. New monitors have recently become available: - The P S A 4000™ (Physiometrix, Massachusetts): this monitor displays an index of hypnosis calculated on 1.25 s E E G epochs. Numerous quantitative measurements are derived (power gradients, absolute power i n 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 EEG).  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 i n 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, B a d Bramstedt, Germany): this monitor is based on a pattern recognition algorithm that classifies the E E G into 14 different classes using Kugler's classification and denomination [68] (A, B0, B l , B 2 , CO, C I , C2, DO, D l , D2, E0, E l , F 0 , 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 introduction 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 - A W A R E 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) i n 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 wellpublicized 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% i n the total amount of anesthetic drugs (inhaled and/or intravenous) could be achieved during the maintenance phase by keeping the B I S value i n 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  TIME [mini  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 ( P A C U ) 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 independently 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.  CHAPTER 2.2.2.2  2. QUANTIFYING  DEPTH  OF ANESTHESIA:  A  REVIEW  28  B a d 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. A s 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 ( E M G ) 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 limited . One can 2  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  2  The contribution of the EMG  muscle movements). Personal interpretation. 3  signal in the EEG  remains unaffected by the NMB (however the drug still effectively blocks  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 PSA  4000 and the NarcoTrend  monitor reported similar findings i n terms of drug usage and P A C U discharge times.  2.3 The  Quantifying the Depth of Analgesia measurement of antinociception is part of the critical care given by anesthesiologists i n 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 2.3.1.1 The  U s i n g Physiological Measures T h e E E G as a Measure of Analgesia  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 (MEF,  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  All opioids depress ventilation in a dose-dependent fashion. A s 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 REVIEW  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 i n 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 C N S . 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 till 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 i n 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 S A , 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. A s 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 ( P K ) 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 ( P D ) . 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 mathematical expression relating the drug blood plasma concentration C (t) to the administered dose I(t). The p  aim of this section is thus to define the transfer function  PK(s)  1  = ^-  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  16  DRUGS:  8  4  THE TRADITIONAL  APPROACH  33  2 T I M E FROM 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 F i g . 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. A s 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  Modeling  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 For  Exponential Models  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 corresponding to the distribution and elimination phase respectively, see F i g . 3.2.  This behavior can be  mathematically expressed as:  C (t) p  = A • e~  at  + B • e-P*  [ng/ml] or [Mg/ml]  (3.2)  where C (t) is the time course of the drug concentration expressed in nanogram per milliliter (opioids) or p  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  C (t) p  = 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 F i g . 3.2), or by using infusion data and analyzing how the plasma concentration increases over time. 2  T h e use of /3 b e i n g usually reserved to m o d e l a n d represent the slow e l i m i n a t i o n phase, researchers have i n t r o d u c e d the  n o t a t i o n P a n d 7r to describe the fast d y n a m i c s corresponding to the d i s t r i b u t i o n phase.  CHAPTER  3. MODELING  ANESTHETIC  DRUGS: THE  T  TRADITIONAL  APPROACH  35  T I M E AFTER INJECTION  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:  Tf = ^  and  r? =  7T  2  —  and  if = ^  Oi  2  (3.4)  [min] p  2  In terms of control and system engineering, the exponential model (3.3) can be directly expressed in the Laplace domain: PK(s)  3.1.2.2  v  ;  =  ^I(s)4 = P -  —  S +  TT  + A- — s +  + B--^— a  s +  (3.5)  /3  v  /  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 compartmental framework, it is necessary to introduce these models here as well. Note that mamillary compartmental 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 ) . This - 1  elimination corresponds to metabolism and hepatic and/or renal excretion. In parallel to the elimination  CHAPTER  3. MODELING ANESTHETIC DRUGS: THE TRADITIONAL  Rapidly equilibrating compartment  APPROACH  36  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 C and C3 of the second and third compartments rise. Once the concentrations in the central 2  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 k i 2  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 equations (expressed here in a state space representation): 10  -&12 +  = C (*).  k  Ci(«)  -&21  0  C2V)  -*3l.  C (t).  2  0  ki3  3  ki  31  l_  +  3  0  i(t)  0  (3.6)  Ci{t) C (t) p  1 0 0 C (t\ 3  where Vi 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., C (t) = C\(t). p  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:  V = V +V +V d  1  2  3  [1]  (3.7)  where V and V3 are the volumes of the rapidly equilibrating and slowly equilibrating compartments. 2  Clearance  clearance  Another important concept used by pharmacologists is the intercompartmental and total drug  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 [ m l - h ] . Using a 2- or 3-compartment model, the total body clearance Cl\ can be calculated -1  as:  C h = Vi • k The intercompartmental clearances Cli , 2  [ml-h ]  (3.8)  -1  w  Cl \, and C/13, CZ31 represent a similar concept.  Instead  2  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\ = Cl \ = Cl 2  2  2  and C Z 1 3 = CI31 = CI3, and that CUj = Vi • kij, we reach the following equalities:  V _ku  V3_ki3  2  Vi " k21 3.1.2.3  V, ~ *si  [  '  Equivalence Between Exponential and Compartmental 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}: 2  - {Vi,P*,  2  standard exponential parameter set,  2  }: exponential parameter set with normalized partial coefficients (i.e.,  A*,B*,TJ,T?,Tf 2  2  2  P* + A* + B* = 1), - { V i , kio, kn, fa, ^ 2 1 , fei}: standard compartmental parameter set, - { V i , V , Vs,Cli,Cl ,CI3}: 2  2  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  38  APPROACH  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)  " V i " (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  Literature Review  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. Propofol  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 B o d y 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. A s 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 ( M A W R ) was about 25% and the mean weighted residual ( M W R ) 4  was -3.4%. The P K parameter set is presented in the Table 3.1 and Table 3.2.  Remifentanil  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. A s 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). A s 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 P K accuracy parameter. It is defined as M A W R = X!»=i  N 52J=I ..M ™—?P?A  >  where ctj is  Cpij  the j  th  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  PK  ANESTHETIC  DRUGS:  PARAMETER  THE TRADITIONAL  APPROACH  VALUE  UNITS  Ch  (BW/70)^ if age < 60 01 (BW/70)^ - (age - 60) • 0 if age > 60 03 (BW/70) • (1 + ven • t? ) • (1 + bol • 0 ) 05 ( B W / 7 0 ) " • (1 + bol • 6 )  [l-min" ] P-min" ] [1-min- ] [1-min- ]  VL  02 (BW/70) » • (age/30)  v v  2  04 (BW/70)  3  06  [1] W [1]  CI i  W  Cl  9s  2  14  16  e  18  e  e9  39  • (1 +bol • 0 ) • (1 + bol • 0i ) ei3  15  7  1  1  1  1  Table 3.1: Propofol P K parameter sets from [101] where B W 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.  PARAMETER MATES  0i 0i 03 04 05 06 07 08 09 010 011 012 013 014 015 016 #17 #18  ESTI-  VALUES  1.44 9.3 2.25 44.2 0.92 266 0.75 0.62 0.61 0.045 0.55 0.71 -0.39 -0.40 1.61 2.02 0.73 -0.48  UNITS  SE  [1/min]  0.09 0.9 0.31 6.1 0.15 43 0.06 0.09 0.11 0.012 0.13 0.26 0.15 0.10 0.36 0.41 0.23 0.12  W [1/min] M P/min] [1]  Table 3.2: Parameter estimates from the N O N M E M analysis ([101]).  CHAPTER  3. MODELING  ANESTHETIC  PK  PARAMETER  Cl Ch 2  V-2 V  3  DRUGS:  THE TRADITIONAL VALUE  2.2 - 0.0162 • age + 0.0191 • LBM 3.25 - 0.0301 • age 0.121 - 0.00113 • age 1.94 - 0.0201 • age + 0.072 • L B M 7.12 - 0.0811 • age + 0.108 - L B M 5.42  APPROACH  40  UNITS  [1-min- ] [1-min- ] [1-min" ] 1  1  1  [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 3.1.4.1  Concluding Remarks 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 i n 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 i n F i g . 3.4 the two-times difference between the propofol impulse response of the infusion vs. bolus models (from Schiittler [101]). A s 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" min 1  CHAPTER  3. MODELING  ANESTHETIC  DRUGS:  THE TRADITIONAL  T I M E [min]  41  FREQUENCY [rad/s]  Figure 3.4: Bolus vs. infusion pharmacokinetic of propofol (using [101]). response.  3.1.4.3  APPROACH  (a) Impulse response (b) Frequency  Time Variance and Other Factors Affecting P K Models  There are a number of physiological differences between individuals that can affect pharmacokinetic parameters. 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 i n 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, i n 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. i n 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 F i g . 3.5. A s 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 Hill equation: E = E + E,'max 0  EClJ+Cl  (3.12)  where E is the steady state effect, C is the steady state plasma concentration, and EC50 is the concentration p  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  E  max  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 F i g . 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 i g . 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  CHAPTER  3. MODELING  ANESTHETIC  o r~  Intubation  Skin Incision  Skin Closure  to  O  DRUGS:  200  THE TRADITIONAL  400  1  MIHM IIMl B ml 111  Response  ,W"|  ."[»*•[«•  TTTn  g  44  600  No response  Ml.,  APPROACH  I  No response Response I  No response  I  Response  \0  °1  § g  50,  "o 2I £  25,  PLASMA ALFENTANIL  |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 EEG  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 i g . 3.8. In this example, the smallest doses used i n 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 k o Parameter: a Measure of the Effect Dynamics e  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 effect . 5  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. An  explanation for this disequilibrium is that the drug concentration within the biophase lags the drug 6  concentration in the blood.  Hence, predicting the effect site concentration C  e  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  k o expresses the dynamics of transfer of the plasma concentration to the biophase. Therefore, it is possible e  to mathematically express the effect site drug concentration C ( s ) as a function of the plasma concentration e  C (s) p  as:  s+ k  e0  Often 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. 5  CHAPTER  3. MODELING  ANESTHETIC  DRUGS:  THE TRADITIONAL  APPROACH  FENTANYL  SPECTRAL  SPECTRAL  [ng/ml]  E D G E [HZ]  E D G E [HZ]  46  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 +  E  0  max  •  C F  r  7  f f  r  ,  (3.14)  7  This modeling approach to pharmacodynamics has become the mainstream approach followed by pharmacologists 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 Hill 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 i n Chapter 6, where a new modeling approach based on System Identification know-how is proposed.  3.2.3 3.2.3.1  Literature Review 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 k o, results vary from 0.1 to >0.6 m i n . This large variation - 1  e  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 pharmacodynamic 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 . Conversely, Billard et al. [117] and Kuizenga et al. [118] reported intermediate values (3 7  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 H i l l non-linear characteristic. For instance, Kuizenga et al. published the Hill parameters for each one of their subjects, see F i g 3.10.a. According to the H i 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 Hill 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  FENTANYL CONCENTRATION  DRUGS:  THE TRADITIONAL  [ng/ml]  APPROACH  48  FENTANYL CONCENTRATION  [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 m g / m l , which is much higher than most of the other studies. Remifentanil  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 H i l l parameters (see for instance F i g . 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 3.2.4.1  Concluding Remarks P h a r m a c o d y n a m i c Interactions  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 i n F i g . 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 of fentanyl) was necessary, using both drugs  simultaneously reduces this requirement to 4 / x g - m l  - 1  -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. A t the extremes, a concentration of either 12 yug-ml of propofol or 12.44 ng-ml"" of remifentanil is -1  1  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% its,  of the patients. A minima at EC^ = 2.5 ^ g - m l 0  - 1  (propofol) and ECl  0  = 5 ng-ml  -1  (remifentanil) ex-  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 Propofol  Time Variance 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  TIME  APPROACH  50  [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.  Remifentanil  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 - m i n _ 1  The  _ 1  was administered over a period of 4 hours.  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  70  DRUGS:  THE TRADITIONAL  APPROACH  51  mean +SD  60 50 40 30 20  o H  10 0 -30  60  90 TIME  _l_ 120  _1_  150  180  210  240  [min]  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 i n 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 i n 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 - m i n ) are linear with respect to the input amplitude, while bolus-type infusions _1  _1  (above 2500 / n g - k g - m i n ) have a different initial drug distribution, which results in a different P K -1  -1  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 interactions 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 i n 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 i n 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 appealing. The anesthetic and opioid titration needs to be constantly adjusted in order to avoid both underand 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. Closedloop 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, closedloop anesthesia would not replace the anesthesiologist. O n 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 century 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 i n the context of anesthesia has not been reported here. In recent years, researchers have been using either a simple P I D or a lookup table of the drug pharmacodynamic 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 We  Summary have seen i n 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)  Bickford et al, 1950-1960 [38] [39] [40]  EEG energy in ,, . , , ° ,, the 4 to 12 Hz , , band  Bellville et al., 1955-1960 [41]  EEG amplitude in a denned frequency band  Schwilden et al, 1985-1995 [44] [45]  CONTROLLED AGENT(S)  CONTROL TECHNIQUE  . ., , Pentothal or _, . Ether  — , , On/Off control '  D  Cyclopropane  POPULATION  Rabbits, cats, monkeys. ™. . , ^ . , -„ Clinical trial on 50 atients under oin pa len s un ergomg . various abdominal surgeries  Propofol ^  Not Disclosed  Model-Based Adaptive Control  13 volunteers (22-29 yr ; 44-85 kg)  [131]  . ,. , . ,, Adaptation was done it the system output was diverging too far from its reference c  Alfentanil  [132]  _ ... ,. , , , , , . Osculations due to the control technique. , ... . . . „ Method sensitive to extraneous interferences, ^ Qp^j^g e been administered concurrently, thus . , .. ... ,, . , , . . . , , seriously limiting this technique m today s practice. ±  I  w  n a V  No results have been presented. The proposed technique is merely an improvement of 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.  Analog control (P or PI?)  Methohexital Median Edge p ° requency  COMMENTS AND LIMITATIONS  Results on volunteers have shown that a constant excitation is necessary the guarantee the reliability of the feedback quantity (otherwise the volunteers were drifting from a drug-induced unconsciousness into a natural sleep). This technique works also for opioids. The controlled drug was used as the only anesthetic agent during the maintenance phase.  11 volunteers (24-31 yr ; 54-87 kg) 11 patients  o  I O O  25 female patients (31-47 yr), ASA I or II  Isoflurane  PI to  Kenny et al., 1990-1995 [133] [134] [135]  „ , , Auditory Evoked „ . , Potentials  Roy et al, 1995-2000 [136] [64]  Gent.ilini et al, 1995-2000 [124]  (125] [18] [137] [138]  n ei Propofol  Halothane Auditory Evoked Potentials  Propofol  Mean Arterial _ Pressure  , . Bispectral index T  c  , Isonurane a  Outer PI controller setting the c ms-n reference to an inner TCI , device.  27 patients  The authors recommended feedback control of anesthesia as a research tool to better identify pharmacodynamic models and the interaction between drugs.  Fuzzy rule-based control system controlling either the vaporizer or giving a reference to a TCI device  10 sessions conducted on 6 mongrel dogs with tail clamping stimulation . 9 sessions conducted on 5 mongrel dogs  These papers emphasize mostly the hypnosis index derived from midlatency auditory potentials using wavelet analysis. Due to the extensive averaging needed, a value quantifying calculated every 3 minutes,  , . „ ,. , . Model Predictive Control  20  Cascade Internal Model Control. An inner loop , -, , controls the end-tidal . .. m, concentration. 1 he reference of the inner loop is set by an outer loop that regulates BIS variations r  atients  patients (  -  Table 4.1: Prior Art: a Literature Survey (Part I)  yr)  t  h  e  [ e y e l  o  f  h  y  p  n  o  s  i  s  w  a  I  s  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 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  3  FEEDBACK QUANTETY(IES)  STUDY  S t r u y s et al., [122]  2001  Bispectral Index  CONTROLLED AGENT(S)  CONTROL TECHNIQUE  Propofol  A lookup table (Hill model) acquired during induction serves to calculate the required effect site concentration changes. A TCI device tracks these changes.  Propofol  Outer PID controller setting the reference to a TCI device. There is an additional constraint limiting the maximum change of the infusion rate  Absalom ei al., 2002-2003 [123] [139]  Bispectral Index  (140]  POPULATION  COMMENTS AND LIMITATIONS  10 female patients (12-60 yr)  10 patients (67 yr ± 1 1 ; kg ± 1 1 )  79  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. 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.  8 o  o o  CO  Locher ei al,, 2004 [141]  Bispectral Index  Isoflurane  Liu et al., 2005 [142]  Bispectral Index  Propofol  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.  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.  10 patients (+13 control) (29-59 yrs  old)  83 patients (closed-loop) ; 81 patients (manual TCI)  Table 4.2: Prior Art: a Literature Survey (Part II)  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. 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.  I  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 i u et al, in particular, reported a 30% decrease in wake up time, and a 23% decrease i n drug usage. Locher et al. have shown that the closed-loop system responds faster to changes in the patient's state and results i n 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 . The lack of a sensor model is a limiting issue from a control engineering 1  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. A s 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. l  This 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. Significant 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 i n 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 i n the Annex C. In Chapter 8, we design a number of controllers using either a classical P I D loop shaping design, or the more sophisticated Hoo weighted mixed-sensitivity design procedure. We show that controllers achieving 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 patient 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 measurements. 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. A s 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 CNS 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 i n 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 i n this section. We then show i n 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 ( W T ) . 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:  w : x  cl  band: [ 0 , / s / 2 ]  d£  band: [fs/2 , fs/2 ~ }  L  L  L  1  dj" band: [/s/2'" , fs/2 } 1  ^ df  (5.1)  1  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 { c t } and {d;}  i=1  . jr, is a time series of coefficients describing the time evolution  of the signal i n 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 i n its continuous form, involving a very large computational complexity ill-suited for real-time applications. In a Discrete Wavelet Transform ( D W T ) 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 \ ... i, =  t  t  therefore depends on the wavelet  filters Ho and Hi associated to the transform w.  5.1.2 The  Change in Phase and Wavelets coefficients obtained through the Wavelet Transform ( W T ) 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  EEG  with increasing anesthetic depth. This is illustrated in Figure 5.2.  5.2  Estimating the Anesthetic Drug Effect: the  Note that contributions in the development of the  5.2.1 The  WAVCNS  WAVCNS  are shared with Ms.  T. Zikov [151].  Concepts and Derivation brain is the target organ of anesthetic drugs.. The electroencephalogram ( E E G ) 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]—>ie[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 h 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. 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. 2  A  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: R — (x ,fc, a  a  k = 1,2,... M}  Rb = {xb,k, fc=l,2, . . . M }  (state a),  (5.3)  (state b),  where each observation vector x fe and x ^ contains a finite number of samples and represent the k  th  Qj  epoch  of the data sets R and Rf,. a  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 {i ,k\k=\,...,M  and {h,k}k=i,...,M- This results in two references which are the  a  averaged features f and ff,: a  >=  • Eifcli o,fe, and f  (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 is c  provided. Using the corresponding feature f = /(x ), it is then possible to estimate how far the system c  c  has evolved from the state a towards the state b by comparing the value f to the references f and ff,. Two c  a  CHAPTER indices j  5. QUANTIFYING a  CORTICAL  AND AUTONOMIC  ACTIVITY  USING  WAVELETS  66  and jp are defined such that: and  ja = life - t i l l  (5.6)  <  V  3P = life - f t l l l  where the norm || • | | i for a vector v is defined as: N  (5.7) The  norm || • | | i quantifies the difference between f and f (or fj,) by integrating the distance between these c  two vectors. The indexes j  a  and jp thus measure a distance of the system from either state a or b. Note  a  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 j  and jp can be combined to yield a  a  single descriptor j = jp — j .  It is convenient here to define a function g that, given an observation x,  a  associates the resulting index j:  g : x — g(x) = j = ||/(x) - ffcHi - ||/(x) - f \\i. a  (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 and j;,: a  9-Ra—>3a  = {ja,k,  g:R —>j  = {j ,  b  The  b  bik  k=  l,2,...M}  k=  l,2,...M}  (5.9)  function h is a scaled version of g, and is thus defined as:  h : x — > h(x) = i = g(x)  i Ja  Ti Jb  Ja  Jb  (5.10)  where: and  Ja  (5.11) T —  Jb —  5.2.1.2 The  X^  M  n  ' l^k=\3b,k-  Feature Function Selection  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 h are adequate as they uniquely define each intermediate state. Conversely, a function such 2  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 h is not very useful 4  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 , T , ..., T/v corresponding to intermediate states of operation {e.g., sedation, general anesthesia, deep 2  hypnotic state, etc.).  For each data set, R or T , we can associate a desired index value i:  fi: {R , T ..., T , R } —- {1, i\,..., i , 0}. a  u  N  b  (5.12)  N  Note that, by definition of h, i = 1, i = 1, if, = 0, and it = 0. a  Linearity  a  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 } of length equal to the number of epochs i n the corresponding n  data set:  h : {R ,Ti,.. a  - ,T , Rb} — • {ia,ii, • • -, i/v, ifc}-  (5.13)  N  We can assess the characteristic of h by simply comparing the average vector values i  n  the set of desired values i  n  = mean(i„) with  i n (5.12), see F i g . 5.4.b. The parameter L is the normalized root mean square  error between the characteristics of h and h:  N slA' •  (5  '  14)  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 F i g . 5.4.a. would typically have  L = 1). To optimize the linearity of h, it is therefore necessary to minimize L. The variability of h can be assessed by calculating for each vector i  Variability  n  of (5.13) the standard  deviation cr(i ). If a statistically significant discrimination of each intermediate state T i , T2, . . . , T/v is a n  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 F i g . 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:  (5.16) where a(io) = cr(i ) and <r(ijv+i) = cr(h)- For the same reasons outlined above, it is advised to keep V < 1, a  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: - R  a  (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 T , T 3 and Rb were collected from surgical audit cases. While a variety of anesthetic 2  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 T and T3 was done following the anesthesiologist's assessment for 2  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 T set. 2  Each E E G signal was recorded from a frontal differential channel ( F i - A T i (the left outer malar bone) p  and F  p z  as ground), using the Crystal Monitor Model 1 6 ™ (Cleveland Medical Devices Inc., O H ) , a  CHAPTER  5. QUANTIFYING  CORTICAL  AND AUTONOMIC  ACTIVITY  USING  WAVELETS  69  Figure 5.5: Normalized EEG 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 REM 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  EEG  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 RMS  value prior to analysis:  RMS(x)'  (5.17)  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: ^WAVCNS  5.2.1.4  Feature Function  :  {Ra,Ti,T ,T , 2  3  Rb} — • {1,0.75,0.5,0.25,0}  (5.18)  /WAV  In order to extract information from the E E G signal pertinent to depth of anesthesia, the Probability Density Function ( P D F ) of the wavelet coefficient set df was chosen as the feature function, and we denote  CHAPTER  5. QUANTIFYING  it as / W A V C N S  CORTICAL  AND AUTONOMIC  ACTIVITY  USING  WAVELETS  70  :  /WAVCNS = * — • / 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 i n the fact that the area of a P D F curve is, by definition, always equal to 1. Therefore, the indexes j  a  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  Selection of the Frequency B a n d and Wavelet Filter  Let us denote with W A V C N S the index i defined in (5.10), and corresponding to the feature function /WAVCNS • 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 F i g . 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  ds  USING  WAVELETS  71  band  (16-32  \  1  2  di  3  -®—-©  band  4  5  6  7  9 10 11 12 13 14 15 16  8  DAUBECHIES W A V E L E T FILTER #  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 COEFFICIENTS  Figure 5.7: Awake and isoelectric reference PDFs based on /W^ AAV C N S  '  CHAPTER  5. QUANTIFYING  CORTICAL  AND AUTONOMIC  ACTIVITY  USING  i v  V>Ss  's/l  \  d: band  \  Filtered index (Average filter of 30 sec)  iS.  ^  fom</  \  \  0.6  &  72  Unfiltered index  'a  0.8  WAVELETS  (16-32H?)  _\ ; .  0.4  !  U.2  di band (64-128  1  '  C* v^-  s  i—  1  1  r  7",  \  Ss  H%)  1  T  2  1 Rj  3  Figure 5.8: Characteristic of the optimal function ^ W A V C N S ^ ' sampling frequency of 256 S/s. (a) Linearity characteristic, (b) W A V C N S values for the data sets Ra , T±, T2 , T 3 and Rb . Note the larger variability in the intermediate data sets. o rt  l e  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  *  /WAV NS C  /wlv  s t  n  u  defined as:  s  /wlv  =* —  C N S  C N S  PDF(dr ), [Dbel  (x) =  (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  /^AVCNS  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 P D F s , see  Section 5.2.2.2. The corresponding characteristic of  ^VVAV NS C  ^ Pl° s  t t e  d  F i g . 5.8.a.  m  Results show a nearly linear  characteristic where: ^WAVCNS  :  i °> R  T  l> 2> T  T  3>  R  b )  ~  *  i  1  0  0  - > > > °) 77  54  20  (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 F i g . 5.8.b.  The variability i n  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 i n a higher non-stationarity of the T data sets.  CHAPTER 5.2.2  5. QUANTIFYING  CORTICAL  AND AUTONOMIC  ACTIVITY  USING  WAVELETS  73  P r a c t i c a l Issues a n d 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 i n proximity of the E E G acquisition device, strong 50/60 H z and 100/120 H z 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™' ' coefficients. When higher sampling rates are used, the wavelet coefficients Db6  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 i n 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, O A s and head movements  typically perturb the E E G from 0 H z 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 compared 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 i n the lower frequency bands (up to 16 Hz) is sufficient to remove the perturbing artifacts from the source E E G ,  see F i g . 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. A s a result, 5waves 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 EEG 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 O A s (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 N S  Normalization  C  Algorithm  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 RMSV  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 F i g . 5.7) to any isoelectric E E G epoch: •^WAVCNS  The Redundant Wavelet Transform  ' ^^^*  s o e  ' *" e c  c  *  ft>  -  (5.22)  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 i n better temporal resolution in the decomposition bands, and, therefore, smother reference P D F s . 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 EEG)  after thresholding i n the wavelet domain [153, 154].  However, the utilization of the S W T somewhat increases the computational demand, since the complexity of the S W T is 0(N- log N) 2  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:  B  w  where B  w  = 3.49 • a - . A T / , 1  3  (5.23)  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 EEG,  the variance is almost equal to zero (see F i g . 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 Trending  Post-Processing 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  2 order low nd  o -20  0  180 0.001  0.1  0.01 FREQUENCY  [HZ]  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 I R 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. Scaling  To allow for a better comparison with the BIS Monitor (see Section V I ) , 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 N S to the BIS v.3.4 (Aspect C  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 F i g . 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 F i g . 5.12. There is a strong correlation  CHAPTER  5. QUANTIFYING  CORTICAL  AND AUTONOMIC  100,  ACTIVITY  Total steady state time  L  - 13 his 16 min  60  CO  L  2  i  mean+2-SD  /  20  >  40  3 L  *•  2li  J 40  =12.5  \ mean =  mean-2-SD =  -40  -9.8  I I I  1  -80 •100  BIS v.3.4  00  0  1. 4 m  K i.  -20  -60 t = 0.969, p = 0 95% conf. int. = [0.9685 0.9696] I I 60 100 80  \  0  •  20  78  = 72.5°, 'o of total time  40 60  WAVELETS  100 HO  80  USING  10  20  30  40  50  60  i  70  1  80  1  90  1  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  Transitory Behavior  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.  Induction  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. A s 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  -Loss Of Consciousness (LOC)  ACTIVITY  USING  WAVELETS  79  Loss Of Consciousness (LOO  100 50 I  0  100  »-  1  3 >. 50  \"  t  -20  20  5  1  40  60  80  100  1  120  T I M E [S]  (a)  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] i n a study involving 50 patients. In F i g . 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 Emergence  event.  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 t (Patient Reacting): the patient was reacting (movement, vocalization) to the anesthesiologist's r  touch and/or voice, or to the surgical environment, - P The  P  r s  r s  (Patient Responding): the patient responded to a verbal command.  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 marker of Return O f Consciousness ( R O C ) .  r s  event cannot be used as a  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 event indicates the earliest time when a patient reaction was observed. This reaction rt  was triggered either by the anesthesiologist's probing, or by environmental factors. Since the P corresponds rt  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 event, i.e., when the patients were responsive. This result is comparable rs  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\, T and T3, the 2  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  S-waves  CORTICAL  AND AUTONOMIC  burst-suppression with predominant isoelecticity  !  ACTIVITY  burst-suppression with predominant bursts !  ->*^  USING  WAVELETS  81  gradual increase of 6 and a activity  >r*-  4  3 TIME  [min]  Figure 5.15: Time course of the BIS (v.3.4) and W A V during an episode of burst suppression (patient #4). The corresponding E E G is also plotted for comparison purposes. C N S  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 F i g . 5.15 the time courses of the BIS and W A V C N S during such an episode, as well as the corresponding EEG.  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.3  5. QUANTIFYING  CORTICAL  AND AUTONOMIC  ACTIVITY  USING  WAVELETS  82  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 i n 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 autonomic 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 activation 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 HRV  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 - R S A ) : 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\ -Relaxed 0£  ai  0 V  1.2  100  200  300  400  300  400  r  -0.4 0  1  2  Time [secj  0  100  200 #Beat  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 HRV  exhibits an oscillatory behavior following  the respiration pattern. Hence, changes in the respiratory pattern strongly influence the H F content of HRV  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  HRV.  - /?-blocking drugs: these drugs are essentially used to depress the sympathetic system i n 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 HRV  [162].  - Anti-cholinergic drugs: drugs such as atropine antagonize the action of acetylcholine. This results i n vagal blockade and the disappearing of H F i n the The  HRV.  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 i n blood pressure and heart rate). Therefore, the quantification of the H F activity of the HRV  should provide a reliable metric of the effect of  pain and analgesic/adrenergic drugs on the autonomic system. The  fact that HRV  analysis can yield a usable 'depth of anesthesia' index was already known i n 1985.  Unfortunately, the very nature of the HRV  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 - ECG can  signal: temporal resolution: the variability in terms of the R - R interval during periods of surgical stress be as low as few milliseconds. A s a result, the HRV  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 Q R S complex is especially difficult during periods of artifactual activity. Electrosurgical Units (ESU) are the most significant source of interference i n 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. A n y 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 ) and intense 2  surgical stress leading to a significant autonomic response (T3). The sets R  a  and T\ were obtained from  a volunteer performing a breathing relaxation exercise (R , see Figure 5) and during a cold pressure test a  (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  R,  T  ,  USING  T  2  WAVELETS  T  3  85  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: R : relaxed subject (y-scale: 0.2 s/div), Ti: 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. a  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 . A s 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. H R V 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  ECG  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 i n the processing algorithm to limit the effect of corrupted H R V samples on the W A V A N S -  CHAPTER Data  5. QUANTIFYING flow  CORTICAL  AND AUTONOMIC  ACTIVITY  USING  WAVELETS  86  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 i n 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 PDF. Note that, due to the large size of the analysis window (2 minutes), there is no need for further postprocessing filtering. Artifacts  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 , as it would be necessary to wait at least a complete H R V period (in this case 2 minutes) to obtain 1  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 I R 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 till 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 R e p o r t s  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 i n 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. A t 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 i n 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 C N S and A N S monitoring is evident in this case.  Case #2 and  (see Figure 5.19)  The W A V A N S peaked 5 times during this case. The first 2 peaks (0h30  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, l h l O , lh20) were also characterized by an increase i n 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, i n 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 i n 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 WAVANS index correspond to periods of time heavily corrupted by electrocautery artifacts (the H R V signal could not be obtained from the ECG signal). Electrocautery also affected the EEG signal and created data loss in the WAVCNS (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 activation, 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  CHAPTER  5. QUANTIFYING  CORTICAL  Fentanyl bolus  SYS BP (mm H g )  Induction  i  AND AUTONOMIC  ^  USING  WAVELETS  Adequate anesthesia Patient light  .  Fentanyl  :„ bolus jSurgery  120 - (b. {(it 1 fit t-  ACTIVITY  89  Patient responding  Hydrotflorphone B  p  1IA* HO  MAC 1.4 1.2  -  1.0  -  0.8  _  100 eg  80  0.6  ^60  0.4  <  0.2  ^ 40 20  0.0  -  0  OOhOO  OOhlO  00h20  00h30  00h40 00h50 TIME  OlhOO  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 EKG. The blood pressure from the cuff sensor was recorded manually each time this measurement was updated by the anesthesia monitor.  CHAPTER  5. QUANTIFYING  OOhOO  CORTICAL  OOhlO  00h20  AND AUTONOMIC  00h30 TIME  Figure 5.20: Case report #3 - Female patient, 38 years old, of the surgery, skin incision was kept to a minimum.  ACTIVITY  00h40  USING  00h50  WAVELETS  OlhOO  90  OlhlO  undergoing a laparotomy procedure. Due to the nature  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  WAVANS  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 WAVANS  WAVANS  levels {e.g., the abrupt  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  WAVANS-  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). A s a result, the sensor dynamics can be represented by the following L T I transfer function:  Hcm(s)  ,  =  / °  0  — j •— —  s + 2 • C • wo • s + z  (5.24)  S  OJQ  W i t h the current W A V C N S implementation, £ = 1 and UJQ = 0.02 Hz:  " If we assume no a priori  C  N  S  (  S  )  =  (8^W • ^  '  (5 25)  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, T T , T , R } -+ {98,78, 56,19,0} u  2  3  b  [%]  (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 F i g . 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 i n a second order transfer function plus delay:  ^w  =  a » - . + * - U r ' . + !)••"•  ( 5 2 7 )  The difference between the identified (5.27) and the analytic (5.25) functions is negligible in the low frequency 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.  CHAPTER  100  5.  QUANTIFYING  CORTICAL  AND AUTONOMIC  ACTIVITY  USING  WAVELETS  92  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  93  50  o 0 «w -50 100  -  80 cn Z U  >  s  WAVCNS  i(measured)  W A V C N S l(predictcd)  60 40  f  -  /  20 0  3 4 Time [min]  10  15  20  Time [min]  (a)  (b)  Figure 5 . 2 3 : 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.) WAVANS  Transfer Function  Following the proposed implementation, the  dynamic is equiv-  WAVANS  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  WAVANS  transfer function  can be expressed as: 1 #ANS(S)  1 96  95 y» g-i-25-,  .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  the transfer function  HANS  is only a coarse approximation of the  be improved by reducing the  WAVANS  WAVCNS  WAVANS  time course. Conversely,  dynamics. This result can  analysis window and adding a post processing  IIR  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.5  5. QUANTIFYING  CORTICAL  AND AUTONOMIC  ACTIVITY  USING  WAVELETS  94  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 WAV, 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 L T I transfer function. A s 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 i n the analysis and characterization 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... A s 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 ( N e u r o S E N S E ™ Monitor, Cleveland Medical Devices Inc., O H ) .  Chapter 6  A System Oriented Approach to Pharmacodynamic modeling We exposed in details in Chapter 3 the approach to pharmacokinetic ( P K ) 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 intravenous anesthetic drugs have yielded contradictory results and large inter-individual variability i n 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. Furthermore, 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, i n 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 i n detail a P D modeling approach based on standard system engineering know-how. 95  CHAPTER  6. A SYSTEM ORIENTED APPROACH TO PHARMACODYNAMIC  PD(s)  MODELING  96  c, \ Quantified effect  Pharmacodynamic M o d e l  Figure 6.1: Traditional P D 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 A p p r o a c h  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 E  q  tration C  p  and the plasma concen-  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 C within the target organ as a function of the plasma concentration. The transfer function e  PD(s)  is uniquely defined by the effect site rate constant k o, which is the rate at which the effect-site e  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 F i g . 3.9.b. The non-linear function which relates the quantified effect to the effect site concentration is handled by a H i 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 H i 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. Measured P l a s m a Concentration  In most studies, the plasma concentration C is p  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 C  p  and adjust the k o value to collapse the hysteresis loop, and further, identify the dose e  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 i n P D studies has been discussed i n 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 Structure  The effect vs. plasma concentration 'hysteresis' phenomenon has  received considerable attention since the original work of Sheiner. The idea of collapsing the hysteresis loop using an additional fictitious P K effect-site compartment can be particularly attractive to clinicians 1  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 k o affects the width of the loop. A s k o decreases, the loop collapses and then inverts. When k o equals the e  e  e  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 k o, since this constant now captures both the system first-order e  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 H i l l nonlinear 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 L T I 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 i n predicting adequately the system output. For instance, it is clear in the example presented i n Figure 6.3 that the system output will not be accurately predicted for a slow varying low amplitude input signal due to the H i 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  /  /  / V  \  T  Measured output signal •  Predicted output signal  Input signal  Model (b) 1  =0.2  3 O  S"  0=0.1 s-  0=0.05 s-  1  1  r  / /L  m  3 O -O  3 O -0  I  1  —  -A  .„  I rt U  It***  Predicted output signal (c)  Predicted output signal  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 fc o 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. e  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  PR's)  ORIENTED  \  Pharmacokinetic Model  PD(s)  APPROACH  \  Pharmacodynamic Model  TO PHARMACODYNAMIC  \  H(s)  MODELING  101  \ —+ l Quantified effect ScalingJunction E  Sensor  Figure 6.4: Block diagram of the proposed system oriented approach to P D 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 L T I 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. i n 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 i n part for the non-linearity of the dose vs. response relationship. In Sheiner's approach, these dynamics are implicitly included i n the P D model. A s 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 To  Proposed Methodology  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 i n 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 i n 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 nonlinear function. The L T I 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 i n 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. A s compared to the traditional approach, the plasma concentration C (s) is obtained using a well defined P K model derived for p  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  PD(s),  decision is made as to the structure of  which can be any reorder plus delay transfer function: p D ( s ) w  =  E {s)^ ^. b -s 1  e  s  m  C (s)  m  + b . - s ^ + . . . + b - s + ao^ m  1  p  1  1  s + o „ _ i • s"" + ... + a • s + an n  1  x  2 •EC ' 5 0  ~  y  '  Residual analysis should be performed to validate the selected model structure. In this approach, the pharmacodynamic model establishes the relationship between the observable effect  E bs 0  and the plasma concentration. E b 0  s  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  J  CHAPTER  6. A SYSTEM ORIENTED APPROACH  TO PHARMACODYNAMIC  MODELING  103  with respect to the dose. This saturation can be expressed using the Hill relationship, which is essentially a non-linear function that can be represented as a variable steady state gain i n the frequency domain, see Section 6.3.4 and Section 7.2.1.2:  El Eobs = o  5 7  +  E  i  =*  0<E <l  (6.2)  obs  Note that the non-linear element is defined here entirely by the steepness coefficient 7. Finally, the quantified effect E  q  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: E  max  and EQ as follows:  E = {E q  max  — Eo) • E b + EQ 0  (6.3)  s  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  Eq,ss = {Emax ~ EQ) • where C  PiSS  7  Y  '  7  b Eo  (6.4)  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  WAVCNS  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  Clinical Data  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  LMA  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 WAV  during  induction reached values lower than 10. It is indeed possible that a saturation phenomenon can affect the WAV  time course in the lower WAV  Two  values, which in turn would result in a significantly different dynamics.  cases were excluded based on this criteria. Finally, two more cases were withdrawn from the study due to a large E M G  the L M A  activation observed after  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 PD  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 i n terms of hypnosis. 6.2.2 The  P D Identification - T r a d i t i o n a l A p p r o a c h  identification procedure for the traditional P D approach can be divided into three stages. A n illustrative  example is given in Fig. Stage #1:  6.5.a.  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 C  p  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  Figure 6.5: PD identification illustrative example (patient #52). approach  MODELING  105  (a) Traditional approach (b) System oriented  square of the error between the measured and predicted effects is minimized. This yields the effect rate constant k o. This method is essentially identical to that of the 'hysteresis' collapse method. e  Once k o is identified, PD(s) is defined such as in (3.13) and used to estimate the effect-site concentration e  C. e  Stage # 3 : H i l l P a r a m e t e r s I d e n t i f i c a t i o n  C  e  The Hill parameter set {EC ,7} 50  is finally identified using  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 Hill 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 Identification - 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 Concentration Profile  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 c,C<>  MODELING  106  c/(,>  *• PK(s)  \  •  •  •  E,.,W  l/iura/W WAVas effect  J 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 i n Fig. 6.6, where the non-linear saturation function is omitted. The only unknown element is the transfer function PD(s) whose output E  q<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 ( ) is then used s  to obtain a filtered version Cp of the plasma concentration. Using both C / and E  q%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)= r  i  J  [  S  )  ^ =C (s)  E  r  -—  TdS  e  p  s+ k  d  2-EC  (6 6) [ b  50  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 :  H i 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. A t this point, only the non-linear element parameters are unknown. Similarly as i n the traditional approach, a search algorithm determines the H i l l parameters that minimize the R M S of the residuals.  CHAPTER 6.2.4  6. A SYSTEM  ORIENTED  APPROACH  TO PHARMACODYNAMIC  MODELING  107  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 ( P R B S ) 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 L T I 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 C (t) as an input signal to determine p  the time constant k o. This time course can usually be divided into two distinct parts: an initial large e  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 r a d - 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 r a d - 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 i n that band.  6.2.5 PD  Identification Results  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 (for  simplification, we assumed that Ema =0 X  {k o,ECso,-y} e  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 k o and kj, values are well contained within the e  identification bandwidth of the input signal.  CHAPTER  6. A SYSTEM  ORIENTED  APPROACH  TO PHARMACODYNAMIC  TRADITIONAL P D APPROACH PATIENT  PROPOSED P D APPROACH  Hill  LTI  #  fc 0  EC  e  [s-^lO" ] 3  5  0  MODELING  Hill  LTI T  7  d  M  lM6/ml]  kd [s -10- ] - 1  E C  5  0  7  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  10.6  1.9  3.0  9  50.4  2.5  2.6  12.9  2.4  2.4  18  160.5  3.6  3.9  058 066  '  071  8.1  1.8  2.7  20  75.0  2.5  1.9  Mean  8.5  2.0  2.9  20.0  61.4  3.0  2.5  SD  3.8  0.7  0.6  13.9  40.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  Mean  8.3  2.4  2.7  12.9  39.0  3.3  2.0  SD  3.8  0.6  0.7  12.6  19.0  0.4  0.4  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  1.8  G3: >40 - <50 years  072  5.4  1.7  2.5  13  30.0  3.0  1.9  Mean  8.8  3.0  2.5  12.4  34.2  4.0  1.7  SD  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  Mean  8.0  2.9  2.1  10.5  33.1  4.1  1.7  SD  3.1  0.7  0.6  8.8  11.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  KE C ,  s+ k  5  APPROACH  .E,  TO PHARMACODYNAMIC  •  e0  c  -r,  s  MODELING  109  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  M o d e l 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. Direct Residual Analysis  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 Hill dose vs. response relationship (see Section 6.3.4).  In the proposed approach, the L T I 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. A s 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. inspection of the residuals tend therefore to stress out the inadequacy of the traditional L T I model.  Direct  CHAPTER  6. A SYSTEM  ORIENTED  APPROACH  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 P D linear models (Patient #65).  ^  Fast and large amplitude transient  Slow transient  ^  30 20  3 3 a  Systeir -orient* d apprc ach  1 0  V  o •7  I "  10  A  -20  Tr aditiona. appro; ch  -30  («)  3 Pt  1  a /o I p  1  / V  2 1 L  0  /  / niinriflfli  -1 \  -2  I  -3 0  20  40  60  80  100  120  TIME  [sec]  140  160  180  200  220  (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  Prediction Error Analysis  APPROACH  TO PHARMACODYNAMIC  MODELING  111  A better way of expressing this is to consider the prediction error  e. p  Considering that the system can be modelled by an autoregressive A R X model, we have:  A(q)-WAV (k-T ) c m  = B(q)-C (k-T -T )+e(k-T ),  s  p  s  (6.7)  s  d  where e is a white noise, A and B are polynomials of the delay operator q , and T is the sampling time. -1  s  This model expresses the fact that the output is a weighted sum of past C , W A V C N S , and e values. The p  residuals e are calculated as:  e(k • T ) = W A V s  In case the model  c m  (k  • T) -  • C (k • T - T ) - -^e(k  s  p  s  d  • T ),  (6.8)  s  perfectly describes the system, the residuals can be expressed as:  e(k-T )  = -^-ye(k-T ),  s  (6.9)  s  Consequently, the prediction error e defined as e (k • T ) = A(q) • e(k • T ) should ideally have the charp  p  s  acteristics of a white noise. If not, the model  s  did ° t capture all of the system dynamics. The n  prediction error obtained from the residuals calculated from both models are presented in Figure 6.9.b. A s 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  R Jr) e  = ELie (fc-r )-e (fe-T -r) P  s  p  £Liep(*-r-) where N is the number of samples in e . p  s  2  Note that (6.10) is here normalized (i.e., R (Q) ep  autocorrelation coefficient at lag 0 (equivalent to the mean of e ). 2  e  p  = 1) by the  A n easy way of assessing whether  is white is then to check whether the signal yields autocorrelation coefficients contained within the  99% confidence interval of a Gaussian distribution A/ (0,1). The corresponding bounds for the normalized r  autocorrelation function can be calculated as: SD.  ep  = ^8-  VN  (6.11)  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  Independence Between P r e d i c t i o n E r r o r and Input  MODELING  112  Another useful test is to calculate the cross-  correlation between the prediction error e and the input signal C : p  p  Sfc=i (k • T ) • C (k - T — T) e  P  s  p  (6.12)  s  Similarly to (6.10), the cross-correlation is normalized by considering both the autocorrelation of the prediction error at lag 0, and the autocorrelation of the input C at lag 0. Ideally, this should also have the p  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: (6.13) 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 i n 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. Model Validity  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 The  Nominal Propofol P D Models  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 N S time course C  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).  CHAPTER  6. A SYSTEM  ORIENTED  Fast and large amplitude transient < •  APPROACH  TO PHARMACODYNAMIC  Fast and large amplitude transient  Slow transient  y  1  -  <  0.04  System oriented model  Slow transient  >  Syst sm oi iente j mo del  2 U 0.02  a< 3 a o z  § 0.5  c« ink rvals  o  y H  S  V  i •  -0.5  0  1  20  40  60  1  i  80  100 LAG  ;  120 [sec]  W  i  140  1  160  0.02  y  -0.04  I  180  113  5 e>  99% covfidenci intervals  0  MODELING  200  220  0  20  40  60  80  100  120  LAG  [sec]  140  160  1 80  200  220  (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 P D model fails this test, indicating the need for a higher degree model, (b) Independence of the residuals with respect to the input signal C . p  CHAPTER  6. A SYSTEM  ORIENTED PD  APPROACH PARAMETER  TO PHARMACODYNAMIC  SET  COVARIATE  MODELING  114  SIGNIFICANCE  Traditional P D feeO =  8.4  EC = 1.4 + 0.029 • age 7 = 3.5 - 0.023 • age  EC no  X  7  No covariate p < 0.05 between G1+G2 and G3+G4 p < 0.05 between G1+G2 and G3+G4  System oriented P D = 26 - 0.3 • age  T  T  ka  k = 81 - 0.99 • age  EC50  ECso = 1-94 + 0.043 • age 7 = 3.0 - 0.027 • age  d  d  d  7  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 P D models with age as a linear covariate. Similarly, all of the parameters obtained from the system-oriented approach have shown a similar dependence. In particular the EC$o  parameter was shown to increase significantly with age, see F i g . 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  6.3.1 We  WAVCNS-  Clinical Adequacy of the Identification Data  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). A s a result, the C (t) time course may only be a rough p  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 2  This 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  (WAVCNS)-  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  -Ludbrooke/fl/.,2002 •Billard etal., 1997  +  4* I*"  1  le-  3  2e  3e-  3  3  4e  3  -+5e-5  H  6e  This study (allgroups)  J  I  7e  3  1  8e  3  Kuizenga e/ a/., 2001 / (initial infusion)  3  115  1  9c  3  1  10e  3  lie  3  Kuizenga etal., 2001 J (sequential infusions)  Kazama etal, 1999 — [allgroups) Kuizenga et al., 2001 (initial infusion) i—Billard «/a/., 1997 EC  S0  [ug/ml]  — This study (all groups)  f  1  H Kazama etal, 1999J (20-5?jro)  Kuizenga etal, 2001 (sequential infusions) This study (a// groups) /  Kuizenga et al., 2001 (initial infusion)  ~  1  T  Kazama etal., 1999 J (allgroups)  Billard  1997  Kuizenga etal., 2001 (sequential infusions)  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 k o time constant is faster than e  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  EC50  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  S y s t e m Oriented vs.  Traditional A p p r o a c h  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 where Wij is the j  th  = {PEij =  Wij  -  (6.14)  Wij}j=i 2,...,Ni t  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(PE ) = — -J2PE i  The  (6.15)  ij  3=1  1  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 (MAR) Ni  1  MAR(PE ) i  = — -J2\PE \  MAR  (6.17)  ij  1  The  value as an indicator of model accuracy:  3=1  indicate the expected average error between the measured effect and the model output,  ii. Correlation coefficient The correlation coefficient r? is defined for each individual i as: r- = 1  3  Ejii  r~o  (6.18)  ( ^ • - ^ • £ , = 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 Time [s]  90  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. iii. Peak  Effect  The peak effect corresponds to the time t and value w of the maximum W A V N S depth p  p  C  immediately following propofol administration and before airway management, see F i g . 6.13. Note that to reduce sensor noise, t and w correspond to the time and W A V C N S value when the index first p  p  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„ =  Patient-specific models  (ti, - U, ), P  P  and  e  iiWp  = (w  itP  - w" )  (6.19)  i;P  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.  CHAPTER  6. A SYSTEM  ORIENTED  PERFORMANCE  APPROACH  ERROR  TO P H A R M A C O D Y N A M I C  CORRELATION  MODELING  PEAK  EFFECT  i,tp  e  bias  RMS  M A R  CT(MAR)  Median  Best  Worst  °{ -i,u,p) F  N  W  [%]  [%]  >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 offitof patient-specific models obtained from the traditional P D approach. PERFORMANCE  bias  RMS  ERROR  MAR  CORRELATION  PEAK  EFFECT Wp  'P  a  (MAR)  Median  Best  Worst  i,t  p  M  (£i,» )  tp)  £  M  CT  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 offitof patient-specific models obtained from the system-oriented PD approach.  CHAPTER  6. A SYSTEM  ORIENTED  PERFORMANCE  APPROACH  ERROR  TO P H A R M A C O D Y N A M I C  CORRELATION  PEAK  MODELING  1  EFFECT  i,tp  e  bias  RMS  M A R <r(MAR)  Median  Best  Worst  Zi^p  W  N  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 P D approach. PERFORMANCE  ERROR  CORRELATION  PEAK  EFFECT  t  Up  p  bias  RMS  M A R  CT(MAR)  Median  Best  Worst  i,t  *( i,tp)  M  Is]  £  ( i,l»p)  £  p  CT  [%]  £  % [ 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 P D approach.  CHAPTER The  6. A SYSTEM  ORIENTED  APPROACH  TO PHARMACODYNAMIC  MODELING  120  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. A n d 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. O n a more positive note, the performances in terms of peak effect were maintained.  6.3.3 The  D o s e vs. Response Characteristics static dose vs. response characteristics of propofol are presented for both approaches in F i g . 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  WAVCNS  of 40 can be obtained with a  plasma concentration of 2.2 / x g - m l . According to the system oriented dose vs. response, the same effect -1  is obtained for a concentration of 3.3 / ^ g - m l . This is a significant difference considering that a C -1  2.2 / x g - m l  -1  would yield a  WAVCNS  P j S S  of  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 H i l l parameters (see the example in Figure 6.3). This problem is particularly acute when identification data contain mostly transitory signals. A n y 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. W i t h 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 i n 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  100 90 SO 70 60 50 40 30 20 • 10 0 0 1  ORIENTED J  APPROACH  TO PHARMACODYNAMIC 100 p 90 • SO •  \ 18-29 years  -f— —j 30-39 years ••  i 40-49 years - • 50-60 years  2  3  4  5  6  7  8  9  10  70 • 60 • g 50 • 40 30 20 10 00  MODELING  121  A.  1 CO  \  1  2  3  \  4  X..  5  6  7  8  9  10  S T E A D Y S T A T E P R O P O F O L PLASMA C O N C E N T R A T I O N  STEADY STATE PROPOFOL PLASMA CONCENTRATION  Cp [ug/ml]  Cp [(ig/ml]  (b) Figure 6.14: Identified W A V C N S to, Propofol dose-response relationships, (a) Traditional approach (b) New approach. 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 Hill static dose vs. response relationship. Let us denote by  f(x)  the H i l l function. For both approaches, f{x) can be written as:  x'<  (6.20)  where X = 0 . 5 (proposed approach) and X = EC^o (traditional approach). To linearize / ( x ) , we consider 0  0  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) + K -x, T  (6.21)  where: —7-1  7 • x'  (6.22)  CHAPTER  6. A SYSTEM  ORIENTED  1 V,  APPROACH  TO PHARMACODYNAMIC  (s + k )-(s + fc ) 21  y  31  •  (s + it)-(s + a)-(s + p)  MODELING  /  122  WAVCNS  s+ k  tQ  Pharmacokinetics  Pharmacodynamics  Sensor  (a)  1  (s + k )-(s + * )  V,  (s + rt)-(s + a)-(s + P)  21  31  c  -r,s  1-e"  1 2-ECJO  Pharmacokinetics  s+ k  ( 8 S + 1)  d  -WAVCNS  2  Pharmacodynamics  Sensor  (b)  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 + fc )-(s + fc )  V".,  (s + it)-(s + a)-(s + P)  21  ko  31  •  WAVCNS  t  4-EC,  s + Ka  Pharmacokinetics  Pharmacodynamics  (») 1  (s + fc )-(s + fc )  V,  (s + Jt)-(s + a)-(s + B)  21  31  C  1 2-£C  -T,s  Pharmacokinetics  J s + fc  Y  fc  5 0  2  d  Pharmacodynamics  1  •  1-e"*  (8-s + 1)  2  s  • WAVCNS  Sensor  (b)  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 i n 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 plasma concentration of EC50.  WAVCNS  — 50, the system then operates at a steady state  As a result, the overall pharmacodynamic gain of the system is expressed  as: K  e  c  »  =  i~kr  -  (6 23) 0  Equation (6.23) holds true for both traditional and system oriented models. The linearized block diagrams can then be represented as in F i g . 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) Systemoriented 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 i n 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 systemoriented 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 1000  APPROACH  =.  —1  TO PHARMACODYNAMIC  - -114 f — P  —!  |  F-—-I |  |  100  ——1  F=i  T  10  124  j System-c nt a\ e d L apprtjacl  |  o  MODELING  ii  ditiona I oach  T  2  j i  o • rac Iti( in J  -55 -45  -8 a ?  -90  £  -135  /*  ppr >a :h  SY t i orient i pproa :h  -180 0.001  0.01  0.1  F R E Q U E N C Y [rad  •  1  s] 1  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 k .  This actually resulted in the addition  of a single order unity gain transfer function to the static dose vs.  response relationship. This approach  e0  rapidly became the standard i n 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 L T I 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 i n 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  WAVCNS  and BIS  are roughly equivalent in steady state  - this equivalence was already established in Figure 5.12). The similarity between the two relationships  CHAPTER  6. A SYSTEM  ORIENTED  APPROACH  TO PHARMACODYNAMIC  100  20-39 years  90  80  2  3  70 [  t> 5 ^  60 50 I-  9  Ok  40  5/3  30  «  20  125  100  90  ^  MODELING  Kazama et al.  S3 >  'This study (using Gepts etj al. PR parameter set)  40-59 rears  80  Kazama et al.  70 60 50  i  40  7^  !  This study (using Gepts et 2X..PK'parameler.s'et) j...  30 20 10  10  _j  0 0  1  2  3  4  5  6  7  8  i  9  i 10  STEADY STATE P R O P O F O L P L A S M A C O N C E N T R A T I O N  C [ng/ml]  0 0  _i  i  1  2  i 3  i  i  4  5  i 6  ' 7  i 8  i  9  •  10  STEADY STATE P R O P O F O L P L A S M A C O N C E N T R A T I O N  C [ng/ml]  P  p  (b)  (a)  Figure 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 EC o 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. 5  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 P K P D Uncertainty Patient variability is probably the most challenging aspect i n 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 i n 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 i n 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 i n 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 i n Section 7.1, we compare all methods with respect to their efficacy in reducing uncertainty vs. their cost i n 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- A ), where ro is the nominal value of r , s  126  CHAPTER 7. MANAGING PKPD UNCERTAINTY and A € l , s  | A | < 1. The quantity u = s  (TMAX  —  127  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 uncertainty. 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,k \, k i} the corresponding parametric uncertainty 2  3  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: 9.4  lO"  a = 4.8  10~  P = 4.0  10~  7T  =  3  4  5  s)  Is" ]  hi  = 9.9 • 10" • (1+0.51 Ag )  [s"  (1+0.46 A?)  [s" ]  hi  = 7.1 • 10" • (1+0.59 Ag )  [s-  (1 + 0.51 A?)  [s- ]  Vi = 9.3 • (1-r - 0.19 • A^ )  (1+0.40  A  1  1  1  4  21  5  31  1  1  1  [1]  These values were calculated considering a two-times standard error. According to this parameter set, the maximum steady state gain of the PK model is 7.5 • 10 while the minimum gain is only 8.1 • 10 . Even 1  6  4  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 k i and ksi values, while having at the same time the maximum TT, a, (3 and Vi values). This 2  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  G  0  UNCERTAINTY  128  G (s) 0  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 • G [ju) • (1 + w(ju) • A „ ( j w ) ) ,  (7.2)  0  where | | A ( j w ) | | < 1, Go(jw) is the normalized nominal model of the system, and Go = \\G$(jw = 0)|| is the u  steady state gain of the system. The weight function w(jiv) quantifies the magnitude of the unstructured uncertainty, while A (JUJ) expresses this uncertainty as a unity disk in the Nyquist plot. The multiplicative u  uncertainty gain w(ju) is also referred to as relative uncertainty. To better understand this concept, let us consider the example of Fig. 7