UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Advisory system for administration of Phenylephrine following spinal anesthesia for cesarean section Fung, Parry 2004

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

Item Metadata


831-ubc_2004-0452.pdf [ 4.39MB ]
JSON: 831-1.0065520.json
JSON-LD: 831-1.0065520-ld.json
RDF/XML (Pretty): 831-1.0065520-rdf.xml
RDF/JSON: 831-1.0065520-rdf.json
Turtle: 831-1.0065520-turtle.txt
N-Triples: 831-1.0065520-rdf-ntriples.txt
Original Record: 831-1.0065520-source.json
Full Text

Full Text

Advisory System for Administration of Phenylephrine Following Spinal Anesthesia for Cesarean Section by .  Parry Fung  B.A.Sc, University of British Columbia, 2002  A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF T H E REQUIREMENTS FOR T H E D E G R E E OF Master of Applied Science in T H E F A C U L T Y O F G R A D U A T E STUDIES (Department of Electrical and Computer Engineering) We accept this thesis as conforming to the required standard  The University of British Columbia July 2004 © Parry Fung, 2004  THE UNIVERSITY OF BRITISH  COLUMBIA  F A C U L T Y OF G R A D U A T E STUDIES  Library Authorization  In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission.  PcL^ V r  34 /flG / loo H  funtx  Name of Author (please print)J  Title of Thesis:  D e  9  r e e :  D e p a r t m e n t o f  Mv^ory  Date (dd/mm/yyyy)  Sy <?fen*  -por  A d^^'yf  MrVSc Electflc^l  The University of British Columbia Vancouver, BC Canada  grad.ubc.ca/forms/?formlD=THS  ration  P ^ Q M T W  page 1 of 1  fhe^i  fr.pkr!ca  e  2qq <j-  _ _ _ _ _ Year:  '  o-f  h^l^r.'A. J  J  last updated: 24-Aug-04  Abstract Phenylephrine is a drug used at British Columbia Women's Hospital, Vancouver British Columbia to treat maternal hypotension induced by spinal anesthesia for Cesarean Section. Hypotension can cause serious fetus hypoxia, therefore maternal blood pressure must be kept above a minimum level. Phenylephrine dosage is mainly determined in a heuristic manner by the anesthesiologist's experience and observation. Since an overdose of phenylephrine can result in maternal bradycardia and hypertension, an advisory system is developed to recommend the optimal dosage of phenylephrine that ensures an appropriate blood pressure response. Data was collected from patients undergoing Cesarean Section at the British Columbia Women's Hospital for drug response modeling. Preliminary results indicated that the quality of noninvasive blood pressure measurement by the existing cuff sphygmomanometry was a prominent source of model uncertainty. Therefore an algorithm that improves the resolution of the blood pressure measurement using pulse transit time, the travelling time of a pulse wave between two sites, was developed and is presented in this thesis. The refined blood pressure reading was used for patient modeling. Separating phenylephrine's blood pressure response from the spinal anesthesia's is the main challenge for this system identification. Various techniques are discussed and validated in Chapter 3. The advisory system based on the results of Chapter 3 was then developed. When hypotension occurs, the advisory model predictive controller recommends an adequate phenylephrine dose according to the identified internal patient model. The design, tuning and online adaptation of the system are illustrated in Chapter 4. Chapter 5 concludes this thesis and points towards future research in this field.  ii  Table of Contents  Table of Contents 1  2  iii  Introduction and Background  1  1.1 1.2 1.3  1 3 5  Background Objectives Thesis Overview and Contributions  '  Noninvasive B P Measurement by Pulse Transit Time  7  2.1 Literature Review on Pulse Transit Time 2.2 Mathematical Model 2.3 Pulse Transit Time Computation 2.3.1 R peak Detection 2.3.2 PPG Maximum Slope Detection 2.3.3 PTT Computation Algorithm 2.4 Limitations of the PTT-BP Model 2.5 PTT-BP Model Validation 2.6 Conclusions 3  Patient Modeling  8 9 15 15 17 18 18 19 25 N  3.1  MISO System Identification 3.1.1 MISO Method 3.1.2 MISO Model and Model Uncertainty 3.1.3 Results of the MISO System Identification 3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 3.2.1 SISO System Identification Method SISO Model and Model Uncertainty Results of the SISO System Identification 3.2.2 Positive SISO Models Positive SISO Model and Model Uncertainty Results of the Positive SISO System Identification 3.3 Conclusions  iii  27  28 28 30 34 35 38 40 43 44 47 48 51  CONTENTS  iv  4 Advisory System Design 4.1 Advisory System User Interface 4.2 Controller Design and Tuning 4.3 Safety Concerns 4.4 PE Gain Adaptation and SAS Time Delay Estimation 4.5 Advisory System Simulation 4.6 Conclusions  53 55 56 60 60 65 70  5 Conclusions 5.1 Main Contributions 5.2 Limitations of the Advisory System 5.3 Direction for Future Research  71 71 72 73  A Data Acquisition A . l Procedures A. 2 Conclusions  74 74 75  B Advisory System Structure Overview B. l SBP Regulator Interface B.2 SBP Regulator  77 77 80  C Variable Regression Estimation  82  D Matlab Code for Positive System Identification D.l processdata.m D. 2 removespinal.m  83 83 85  E Matlab Code for Noninvasive BP Measurement by P T T E. l PTT2BP.m E.2 rdetect.m E.3 plethdetect.m E.4 findPTT.m E.5 findSBP.m E.6 recalibrate.m  88 88 88 94 99 101 104  F Acronyms  107  Bibliography  108  List of Figures 1.1 Block Diagram of the Advisory System  4  2.1 Illustration of an Arterial Pulse Wave  10  2.2 Total Least Squares between the BP Measured by Cuff and PTT  12  2.3 Definition of PTT  16  2.4 Flow Chart of Maximum Slope Detection in PPG  17  2.5 PTT and its BP Measured at Different Sites from a Single Subject  20  2.6 Comparison between Systolic BP by PTT and Invasive Systolic BP from a Single Patient  22  2.7 Histogram on the Difference between Systolic BP by Cuff and by PTT from 22 Patients Undergoing Cesarean Section  23  2.8 Limits of Agreement between Systolic BP by Cuff and by PTT from 22 Patients Undergoing Cesarean Section  24  2.9 Comparison between Systolic BP by PTT and BP by Cuff Sphygmomanometry during Cesarean Section (Case 10) 3.1 MISO Patient Model  25 29  3.2 Linear Fractional Uncertainty Model with the Time Delay and Integrator Separated  31  3.3 Distribution of First Column of A  32 v  LIST OF FIGURES  vi  3.4  Distribution of First Column of B  33  3.5  Distribution of Second Column of B  33  3.6  Distribution of Third Column of B  34  3.7  Case 14: Individual Model Validation  35  3.8  Case 14: Nominal Model Validation  36  3.9  Case 17: Offset Mismatch during Model Validation  36  3.10 Case 2: Integrating Action Mismatch in PE during Model Validation. . . .  37  3.11 Case 12: Underestimate of K2 during Model Validation  37  3.12 Non-hypotensive Case (Case 8) with No Administration of PE, but Significant Decrease in BP  39  3.13 Least Squares Fit of the SAS Response  40  3.14 Subtraction of the SAS Response From the raw data (same case as in Figure A.l)  :  41  3.15 Distribution of the Two Poles of PE  41  3.16 Distribution of the DC Gain of PE  42  3.17 Patient Model in the Matlab Simulink Representation  42  3.18 Case 10 PE Model Validation  44  3.19 Case 12 PE Model Validation  45  3.20 Case 14 PE Model Validation  45  3.21 Case 19 PE Model Validation  46  3.22 Case 12 PE Model Validation with Customized Gain  46  3.23 Distribution of the Two Real Poles of Positive PE Model  47  3.24 Distribution of the DC Gain of PE  49  3.25 Case 10 PE Model Validation  49  3.26 Case 12 PE Model Validation  50  3.27 Case 14 PE Model Validation  50  LIST OF FIGURES  vii  3.28 Case 19 P E Model Validation  51  4.1  Screen Shot of the User Interface of the Advisory System  54  4.2  Frequency Response of the Open Loop and Closed Loop Nominal P E Model with Tuning Q=1000, R = l , H u = l , Hp=6  59  4.3  Combining B P measured by cuff with B P measured by P T T  61  4.4  Case 47 where V R E Estimation Does Not Agree with Expert Knowledge.  . 63  4.5  Case 54 where V R E Estimation Does Not Agree with Expert Knowledge.  . 64  4.6  Illustration on the Advisory System Adaptation Scheme  64  4.7  Perfect Model Advisory System Simulation  65  4.8  Perfect Model with Noise a = lOmmHg Advisory System Simulation. . . .  66  4.9  Doubled Gain Model Advisory System Simulation  66  2  4.10 Oscillating Model in Section Advisory System Simulation  67  4.11 Perfect Model Closed Loop Simulation  67  4.12 Perfect Model with Noise a = lOrnmHg Closed Loop Simulation  68  4.13 Doubled Gain Model Closed Loop Simulation  68  4.14 Oscillating Model in Section Closed Loop Simulation  69  A. l  Case 10. A Typical Case.  76  B. l  Advisory System Data Flow Diagram 1  78  B.2  Advisory System Data Flow Diagram 2  79  2  List of Tables 4.1  Controls on the Advisory Page of the System  4.2 Indicators on the Advisory Page of the System 4.3  55 . 56  Comparison of SAS Time Delay Estimation by VRE and Expert Knowledge. 62  viii  C h a p t e r  1  Introduction and Background 1.1  Background  During the past decades, engineering knowledge has been widely applied to solve many problems in medicine and biology. Signal processing of electrophysiological signals, computeraided diagnosis tools, digital image processing, control of biological parameters are just a few examples. There are many areas waiting to be explored related to the biomedical engineering field. In this thesis, a particular problem of patient blood pressure (BP) control is undertaken. Maternal hypotension is one of the most common complications during Cesarean Section performed under spinal anesthesia.  Before performing Cesarean Section surgery,  patients are anesthetized by the administration of a spinal anesthetic. The sympathetic blockade due to the spinal anesthetic causes peripheral vasodilatation, which induces hypotension [13]. The low BP causes reduced blood flow and less oxygen delivery to the fetus, which could suffer from significant hypoxia. To avoid this undesirable outcome, hypotension must be closely monitored and treated without delay. Phenylephrine (PE), acting directly on alpha-adrenergic receptors, can be used to treat hypotension [13]. PE raises BP by constricting resistance and capacitance of blood vessels and increasing total peripheral resistance. Whenever the anesthesiologist observes or predicts hypotension, he  1  1.1 Background  2  or she administers a dose of PE to the patient. This causes a transient increase of BP and maintenance doses are needed from time to time. While PE can treat hypotension effectively, the current practice lacks a well-defined guideline on administration dosage. Although it is widely accepted that 20-100/ug of PE is adequate in this application, the large range lacks exactitude and cannot avoid unfavourable outcome following over or under dose. Insufficient PE fails to treat hypotension but excess causes reflex bradycardia and hypertension. Once injected, it is impossible to withdraw the excess PE from the body, therefore it is important to find the optimum dose for the targeted BP response. Previous studies have documented success with closed loop control of mean arterial BP using continuous drug infusion (see  [18]  and  [20]).  Clinical results proved the effectiveness  of the methodical engineering approach. However, in spite of those encouraging outcomes, the Ethical Committee of the BC Women's Hospital were acutely aware of issues in human research. The Committee recommended that patients were not put at risk because of the intended research. Therefore, the protocol of this study was designed such that patients recruited do not experience a higher risk than normal Cesarean Section. Medical procedures remained largely unchanged and researchers did not alter the PE dosage and injection time for the sake of research. Rather, they cope with the limited number of PE administration in the study. That leads to a problem of poor system excitation for patient modeling, which was discovered to be the main challenge of the project. On the other hand, accurate but invasive BP measured by an arterial catheter cannot be used due to its risk. The insertion of the catheter may cause skin infection and arterial injury. Therefore, the group investigated a continuous noninvasive algorithm to infer BP based on pulse transit time (PTT). As closed-loop feedback to BP control cannot be directly applied until the performance of the system is verified by a preliminary study, an advisory system was developed to guide the current open loop dosage decision. The advisory system passively suggests a PE  1.2 Objectives  3  dose for hypotension treatment to the attending anesthesiologist. The advisory dose is filtered by the anesthesiologist's expert knowledge so that the system can be evaluated in a safe environment with the express expert endorsement. After clinical evaluation of the advisory system, the system can be tested in closed-loop operation. As the half-life of PE is short [13] and its administration is frequent (roughly a dose every 5 minutes), a previous study suggests that continuous PE infusion is more appropriate for hypotension treatment [22]. Ultimately, the patient's BP can be controlled automatically in a closed-loop manner with minimal supervision.  1.2  Objectives  The advisory system developed in this thesis aims at solving the PE dosage problem systematically by applying modern control theory. Since the criterion for the system is analogous to the objective of controller design - minimizing tracking error and amount of input, tackling the problem from the control perspective is rational. The goal of the advisory system is to promote safe drug delivery tailored to each patient, without modifying the clinical procedures or introducing new monitor equipment. The setup time has to be relatively short and the interface has to be user-friendly. As the Datex AS5 monitor is used in the operating room at the British Columbia Women's Hospital, the group decided to employ the plugin feature of the Datex S5 Collect data collection software (Datex, Finland). S5 Collect allows real-time data access from the existing Datex monitors. As a result, with an electronic infusion device and a laptop, the advisory system can be implemented without additional costs. There are three design criterions for this project: First, the controller has to utilize a patient model that clinicians can examine and validate during adaptation hence the controller has to be model-based. Second, dosage constraints, crucial safety reasons during drug administration, have to be imposed. Last but not the least, the prediction of BP  1.2 Objectives  4  trend has to be displayed in order to make the system more user-friendly and trustworthy to the anesthesiologists. An adaptive model predictive controller (MPC) was chosen as the backbone since it meets those design requirements. The first criterion rules out the use of a simplistic PID (Proportional Integral Derivative) controller. As constraints can be easily implemented as part of the MPC framework, the second criterion is satisfied. Lastly, the prediction nature of the MPC makes the forecast of the BP trend readily available. With a nominal patient model, the MPC of the system is designed such that the inter-patient variability, constraints, controllability and the frequency response of the system are carefully considered. In real-time, the robustness of the system is improved by a Kalman filter. Other than noise filtering, the Kalman filter also acts as a state estimator which is able to predict the BP trend in conjunction with the MPC prediction. The block diagram in Figure 1.1 demonstrates how the advisory system aids in deciding the optimal dosage of PE administration assessing the real-time BP of the patient.  PC Based Model Predictive Controller  Suggested Feedback Control Input  Anesthesiologists' Approval  Phenylephrine Infusion t  Patient  -, Sensors  -1  Anesthesia Monitor  Blood Pressure Measurement, Electrocardiogram, Plethysmogram. etc  Figure 1.1: Block Diagram of the Advisory System. As a patient model is essential for the MPC, input (PE), output (BP) and disturbance (spinal anesthetic) signals were recorded and analyzed using system identification tools. PTT was employed to improve the sampling rate and resolution of the BP measurement. An algorithm was developed to separate the response caused by input and disturbance.  1.3 Thesis Overview and Contributions  5  Afterwards, conventional linear time invariant system identification is performed on the data. After data collection and patient modeling, an advisory system meeting the three design criterions was developed. With the initial nominal patient model, the system was able to compute the optimal dose once it had determined the accurate patient model through online adaptation. It is also able to regulate BP in closed loop. The system is expected to minimize the undesirable effects of under and over dose and calculate the correct dose of PE for treating maternal hypotension induced by spinal anesthesia for Cesarean Section.  1.3  Thesis Overview and Contributions  The thesis describes in full details the development of the advisory system for Cesarean Section spinal anesthesia. Here are the three main contributions of this thesis: Continuous noninvasive BP measurement - Following an introduction, Chapter 2 reveals how the need for improvement in the sampling rate of the noninvasive BP initiated the development of the continuous noninvasive BP measurement by PTT. Although the PTT has long been known to vary with BP, their relationship was not clear. The proposed simplified physiologic model captures the critical and fundamental properties of the arterial blood flow. The model was proven to be capable of improving the resolution of current noninvasive BP measurement [10]. Patient model for P E and spinal anesthetic - In Chapter 3, different system identification techniques and patient models are presented along with the corresponding uncertainty and validation. The models not only benefit the advisory system design, but also provide a well-defined guideline to anesthesiologists on drug response. Advisory system for Cesarean Section hypotension treatment - The advisory system developed based on a MPC is able to treat hypotension timely and accurately.  1.3 Thesis O v e r v i e w a n d C o n t r i b u t i o n s  6  Utilizing control engineering, the system is expected to deliver better patient care and improve anesthesiologists' efficiency. In Chapter 4, the design and tuning of the predictive controller in the advisory system is discussed. The technical details of the system such as the MPC and Kalman filter structure and the gain adaptation are disclosed. Simulations of the system are presented to verify its performance. At the end, conclusions and remarks summarize the contributions of the system, state the limitations of the system and suggest directions for future research.  Chapter 2  Noninvasive B P Measurement by Pulse Transit Time A continuous measurement of BP is desirable for consistent patient care and monitoring, especially when treating hypotension with PE. The data collected, described in Appendix A, suggests that PE is administered roughly every 300 seconds by the anesthesiologist (see Figure A . l in Appendix A). Conventional noninvasive BP measured by cuff sphygmomanometry updates a BP reading at the most once every minute (usually 3 minutes), which is not frequent enough for hypotension treatment following spinal anesthesia for Cesarean Section. Although BP can be measured continuously by an intra-arterial catheter, this costly and invasive method introduces risk to the patient and workload for physicians. Risks of arterial injury and skin infection do not justify its use in Cesarean Section. Moreover, the rather complex setup procedure for arterial catheter insertion may take up to 30 minutes. Expensive disposable equipment adds to the cost of direct arterial measurement. Therefore, a noninvasive method for measuring BP is desirable for Cesarean Section patients without significant organ dysfunction during short duration surgery but with expected fluctuation in BP. An algorithm for inferring BP from PTT was developed for the advisory system which aims to treat maternal hypotension induced by spinal anesthesia performed for Cesarean Section.  7  2.1 Literature Review on Pulse Transit Time  8  In this chapter, the literature on PTT is briefly discussed. Then a human model of BP and PTT, published in [10]), is presented. The model simplifies the body structure and relates BP to PTT by fundamental physics (i.e. the conservation of energy). To reduce inter-patient variability, the model customizes parameters by using an accessible physical property of the patient (i.e. height). The practicality of the model is retained by ignoring or approximating unavailable parameters such as the compliance of arterial wall, heart pre-ejection period (PEP) and blood density. This chapter also includes a method for computing PTT from the electrocardiogram (ECG) and plethysmogram (PPG). The Matlab implementation of the algorithm is included in Appendix E. Because ECG and PPG are two standard signals in most anesthesia monitors, the PTT-BP model in conjunction with the PTT computation algorithm provides a low cost alternative to continuous invasive BP monitoring.  2.1  Literature Review on Pulse Transit Time  It is commonly accepted that PTT, the time the pulse wave takes to travel between two sites, varies inversely with BP [24]. A number of commercial BP monitors use PTT to infer BP [16], [21]. However, the details of the actual algorithms have not been disclosed by the manufacturers. The first pulse wave model was published back in 1922 by Bramwell and Hill [3] using the following equation to relate pulse wave velocity with blood density, arterial cross section area and compliance:  v where  v  wave speed  a  arterial cross section area  P  blood density  a  Aa/AP  arterial compliance  (2.1.1)  2.2 Mathematical Model  9  As the compliance C is denned as A a / A P , the above equation can be used to relate a  pulse wave speed to BP. Although the model is capable of inferring the BP if the parameters are known (a, p, C and Aa), those parameters are not available unless we use technology a  such as magnetic resonance imaging (MRI) and ultrasound. Since the proposal of the Bramwell and Hill model there have been numerous studies investigating the relationship between PTT and BP but most of them presented only the correlation between P T T and BP in their results [8], [23], [26]. There is no doubt that changes in BP are reflected in PTT. However, for monitoring purposes, a BP estimate that closely resembles the existing noninvasive or invasive measurement is demanded. A practical mathematical model published that can convert PTT into a easy-to-use BP index is still to be found.  2.2  Mathematical Model  The model assumes laminar blood flow from the heart chamber to the fingertip through a rigid pipe, the artery. While it is well known that the artery wall expands and contracts, its small compliance of 0.0018 liter/mmHg on average [15] justifies the assumption that the energy lost during expansion and contraction should be neglected. The model, using the pulse wave velocity, estimates the pressure difference between two sites: the heart and the fingertip. A pulse wave travels from the heart to the fingertip, along the artery and its average velocity can be estimated from the distance travelled divided by the PTT. The relationship between the PTT and BP is demonstrated in the following postulate. The work done by the pulse wave can be expressed in terms of the kinetic energy of the wave and the gravitational potential energy:  10  2.2 Mathematical Model  F •d where F d m v g h  = = = = = = =  ^mv + mgh force exerted on blood distance from heart to fingertip mass of blood pulse wave velocity 9Mm/s height difference between two sites 2  (2.2.1)  2  Figure 2.1: Illustration of an Arterial Pulse Wave. Figure 2.1 illustrates how a packet of blood travels along the artery as a pulse wave. In reality, the blood mass m does not move all the way from the heart to the fingertip. Instead, the pulse is transmitted as a mechanical pressure wave. Still, the transmitted energy is essentially the same as the kinetic energy of the blood mass m, with the assumption of non-viscous flow, hence Equation (2.2.1) holds. Like the Bramwell and Hill model [3], Equation 2.2.1 ignores the speed of the actual blood flow, which is negligible compared to the pulse wave speed. The pumping force from the heart can also be written in terms of BP difference:  2.2 Mathematical Model  11  F = ABP • a  (2.2.2)  where a is defined as the cross section area of the artery. By substituting Equation (2.2.2) into (2.2.1) and after rearranging we find: ABP = \—y + ~gh la - a a• a The term  (2.2.3)  = p represents the density of blood. The pulse wave velocity v can be  expressed as -p^f, so the fictitious mass m disappears from the equation: ABP = ±  P  ^ + pgh  l  (2.2 A)  The distance d can be approximated from the patient's height. PTT is the PTT in seconds. The average blood density for human p is 1035 kg/m [5]. 3  As illustrated in Figure 2.1, the pressure drop in the arterial side of circulation accounts for roughly 70% of the total pressure drop in the body [15]. The venous BP, about 2% of the overall BP [15], can be ignored. Therefore the patient's overall BP is approximately:  BP  = ABP/0.7 = MhPprn+Pah)  i  A —  PTT'  (2-2.5)  p  2  In summary, the BP can be written in terms of PTT, with two parameters namely A and B, where A can be estimated from the subject height as: A = (0.6 x height) x £ 2  (2.2.6)  From the above calculations, BP can be estimated from PTT and a few average empirical values. Recursive calibration between the estimated BP and the BP measured by a cuff is necessary in order to obtain the absolute BP of the patient using estimated values.  2.2 Mathematical Model  12  Calibration is performed using total least squares as it is an optimum method for reducing the uncertainty of two noisy signals [2]. Total least squares minimizes the sum of perpendicular distance from a data point to the line L = j? -> + B, as shown in Figure 2.2. p  r  Figure 2.2: Total Least Squares between the BP Measured by Cuff and PTT. The line L and its corresponding cost function f(L) are expressed as: r  _  f(L)  = EililKII =  A  , p  (--) 2  2  7  J2iLi II (BPi - (pff? + B)) x ^,^2+1 II  In order to minimize f(L) for size N vector BP and PTT, the following least squares  2.2 Mathematical Model  13  problem has to be solved: 1...0 B =  BP  : -.0  PTT  2  where  0...1 BP  l  BP  BP  2  (2.2.8) BP  .  N  PTT  X  PTT  PTT  2  PTT  N  Since the above least squares problem is relatively simple, singular value decomposition (SVD) [12], though computationally expensive, can be used to find the exact solution for parameter B: 1...0 : •••0  =  UY,V  T  (2.2.9)  0...1 B =  VY,-W BPT  PTT ! 2  The parameter B is updated as the weighted average of the current and new value, whenever a new BP measurement is available. The calibration uses the current 1000 seconds of BP readings in a sliding window to update the parameter B. As the sampling interval of BP is 10 seconds, total least squares fit utilizes the 100 data points in the sliding window such that the fit is statistically significant. Since A does not vary significantly between subjects, performing an adaptation on A may lead to unstable estimation. The calibration will adapt B only. A test was performed to confirm that adapting A is not feasible. The adaptation of B may still absorb some of the mismatch in A.  2.2 Mathematical Model  14  The BP in Equation (2.2.5) describes the mean BP assuming laminar flow. The cuff measurement of diastolic BP, and hence the estimation of mean BP, may be impossible in pregnancy due to the changes induced in the arterial wall by gestational hormones. When the BP cuff operates in continuous mode during the treatment of hypotension induced by spinal anesthesia for Cesarean Section, the constraint of measuring mean BP would result in a significant reduction in data points. Systolic BP measurement is also used by anesthesiologists to detect hypotension in current clinical practice. Therefore, based on the assumption that systolic BP is highly correlated with mean BP, PTT is used to infer systolic BP. An equation similar to Equation (2.2.5) can be also derived directly from the Bramwell and Hill equation, therefore one can rearrange Equation (2.1.1) and get:  2  AP=^Aa a  (2.2.10)  As in Figure 2.1, A P is the pressure difference between the systolic and diastolic BP. A P influences the arterial expansion by a factor of the compliance Ca- The systolic BP can be written in terms of A P and diastolic BP:  Systolic BP  = A P + Diastolic BP =  i ^ A a + Diastolic BP  (2.2.11)  a  =  ^ A a y y 3 + Diastolic BP p  Then one may be able to estimate the systolic BP by measuring arterial cross section area. However, noninvasive measurement of a and A a were not possible until the recent technological advances in MRI and ultrasound, which are still not commonly available in most operating room due to their high cost, ft is worthwhile to mention the structural similarity between the two Equations, (2.2.4) and (2.2.11). Both of them state that BP is inversely proportional to the square of PTT.  2.3 Pulse Transit Time Computation  2.3  15  Pulse Transit Time Computation  PTT in this thesis is defined as the time between the ECG R peak and the corresponding maximum inclination in the PPG, as illustrated in Figure 2.3. The R peak of ECG marks the contraction of the heart and the maximum inclination of PPG stamps the arrival time of the pulse wave at the fingertip. It is believed that the effect of pulse wave reflection is negligible at the rising edge of PPG. Therefore the PTT detected is the PTT between the heart and the fingertip. Theoretically the foot of PPG, i.e. the beginning of the rising edge, could be an alternative marker for pulse arrival. However, since the magnitude of PPG is calibrated continuously by the sensor, the administration of PE often causes saturation in the PPG signal due to a rapid increase in BP and in P P G magnitude. To avoid the effect of this saturation, the maximum slope of PPG was chosen as marker. The detection of PTT involves peak detection in both ECG and the discretely differentiated PPG. A wavelet-based approach was taken to decompose the signals into different frequency bands and a rule-based expert system was used to analyze the wavelet components. 2.3.1  R peak Detection  The detection of the R peak has been extensively covered by previous authors. A wavelet implementation based upon the techniques described in [4] was implemented by Chris Mott, a graduate student at the Computer and Electrical Engineering Department of University of British Columbia [10]. Numerous approaches have been proven successful and a good review of algorithms is provided in [17]. For this application an algorithm providing efficient real-time implementation was desired. A wavelet implementation based upon the techniques described in [4] was chosen. The ECG signal was decomposed into 512 sample segments, and a five level stationary wavelet transform (SWT) of Daubechies-1 and biorthogonal-1.5 [27] was  2.3 Pulse Transit Time Computation  16  P T T from E C G and P P G -  - PPG P P G max Slope ECG * E C G R peak  O  100  1550.2  1550.4  1550.8 Seconds  1551  Figure 2.3: Definition of PTT. performed. The highest frequency decomposition captures much of the noise, and the lowest frequency decomposition captures baseline trends. The intermediate frequency bands thus contain most of the energy of the QRS waveform and a recursive, rule-based pattern recognition algorithm identifies the locations of the R-peak patterns. These are primarily characterized by a local minimum followed by a local maximum in the stationary wavelet decomposition. One of the significant modifications made from the approach suggested in [4] is that cubic spline interpolation is used to upsample the RR-intervals and improve the accuracy of the peak location measurement. After an R-peak has been identified by its pattern in the wavelet decompositions, a segment of the original E C G ' s samples surrounding that location is upsampled by a factor of ten, and the point of the local maximum is recorded. With an original sampling rate of 300Hz the resolution after upsampling is 0.33 ms. The sequence of R-peak locations is provided continuously to the PTT algorithm.  2.3 Pulse Transit Time Computation  2.3.2  17  P P G M a x i m u m Slope Detection  The detection of the location of the maximum slope of the PPG signal is also accomplished using a wavelet decomposition. Since the PPG is a relatively simple signal, a 2-level SWT with Daubechies-3 wavelet was used to smoothen out and filter out the noise in the 512 sample segments of the differentiated PPG. The rule-based system, as displayed by a flow chart in Figure 2.4, processes the low frequency wavelet components and detects the maximum slope for each cycle of the PPG. The algorithm detects a pattern of upward sloping followed by downward sloping and identifies the peak between the two slopes. Inactive signal (no measurement) and artifacts are also identified and discarded by the system. The flow chart presented in Figure 2.4 shows the infinite loop offering the real-time application of the proposed algorithm. Initialize dPPG by filtering differentiated PPG by SWT  up_slope_counter =0; down_slope_counter =0; flat_counter =0; look_forward_counter =0; temp_max = null; temp_max_index= null;  Store temporary maximum  If dPPG(i) > temp_max & dPPG(i) < MAX & dPPG(i) > MIN  temp_max = dPPG(i); temp_max_index = i; up_slope_counter++; flat_counter - -; look_forward_counter=0;  Yes  _NoJ  If up_slope_counter > UPSLOPEMIN & down_slope_counter > DOWNSLOPEMIN & flat_counter < FLATMAX & look forward counter > HORIZON  Nc4  If temp_max==null  No  I  look forward counter++  TJ  Yes  Reset Store Set up_slope_counter =0; temp_max_index)-*| down_slope_counter =0; as Peak flat_counter =0; look_forward_counter =0; temp_max = null; temp_max_index= null;  Yes _Yes.  flat_counter++; lfdPPG(i)==0 1 No~ lfdPPG(i)>0 No I . . I Yes up_slope_counter++; down_slope_counter++;  Yes W If flat_counter >= FLATMAX  Figure 2.4: Flow Chart of Maximum Slope Detection in PPG.  No  x 2.4 Limitations of the P T T - B P Model  2.3.3  18  P T T Computation Algorithm  In our current set of data, ECG and PPG were collected at 300Hz though the P P G was interpolated by the monitor from 100Hz data. The threshold constants described in the flow chart are a set of values in terms of the patient's current RR interval (RRI), the inverse of heart rate, where UPSLOPEMIN FLATMAX  = 30 • RRI, DOWNSLOPEMIN  = 15 • RRI and HORIZON  = 5 • RRI,  = 100 • RRI. MAX and MIN are 25% and  5% of the saturation of the PPG signal. The dPPG is the filtered discretely differentiated PPG and i is a dPPG array counter. After R peaks of ECG and maximum slopes of PPG are detected, the corresponding pairs are mapped together to compute PTT. For each R peak, the algorithm checks the next maximum slope and validates whether the PTT lies within 20% to 85% of the RR interval. The confirmed PTT is then stored for BP computation in Equation (2.2.5). Theoretically, the pulse to pulse BP can be computed. In practice, a sampling interval of 5-10 seconds for BP is clinically adequate. Extreme values (noise) are filtered out by a probability density filter. A median filter takes the median of several PTT readings in each sampling interval into one BP index. The parameter B in Equation (2.2.5) has to be calibrated consistently whenever a new cuff measurement is available. As described in the previous section, recursive total least squares calibration is performed to converge the BP measured by PTT to the BP measured by cuff.  2.4  Limitations of the P T T - B P Model  There are several limitations to the model. On one hand, the three signals, ECG, PPG and cuff BP have to be relatively noise free. Since the PPG is sensitive to subject movement, the BP index from PTT can be inaccurate or missing at times of vigorous muscle activity. Body movement also affects the height differential between the two sites hence the potential  2.5 P T T - B P Model Validation  19  energy in Equation (2.2.1). The recursive calibration can only adapt slowly to the change in height. In addition, it is mandatory that all three signals, especially the ECG and PPG, be synchronized correctly. It was observed that some PPG sensor units have a built-in filter. The filter generates an artificial phase lag between the ECG and PPG, which inadvertently lengthens the PTT and distorts the BP index. On the other hand, as the R peak marks the electric excitation of the heart contraction, there is a small delay before mechanical contraction. This delay is the PEP. As we cannot noninvasively measure PEP, the PTT detected includes the PEP, which is negligibly short compared to the actual PTT. It is speculated that for subjects with fast heart rate, the PEP could become more significant. Since PTT is probably more sensitive to heart rate than PEP, PEP may become relatively longer for these subjects. Therefore the model may have to estimate the PEP based on the RR interval for these subjects.  2.5  P T T - B P Model Validation  In this section, the PTT model is validated by measuring PPG at different parts of the body. Then the PTT technique is compared with existing invasive and noninvasive BP measurements. First, to validate the relationship described in Equation (2.2.6), a test case with PPG measured at different sites was recorded, as shown in Figure 2.5. On the same subject, the PPG was measured at the right finger, ear, left finger and toe. According to the model, the parameter A of those different sites should be directly proportional to the square of its distance from heart. Assuming the ratio of the distance from the heart to the fingertip is 1, the ear 0.6 and the toe 1.5, we substitute the ratios into Equation (2.2.6) to calculate A. The BP estimated by PTT is plotted at the bottom of Figure 2.5. The consistent results reinforce the value of the model and the assumptions made. Next, the algorithm estimates BP from real patient data. During surgery, both invasive  2.5 PTT-BP Model Validation  20  PTT of right finger, ear, left finger and toe  Right finger  x  x  xx xx xSxv&x*,"  Left finger  x  * * xv X X x  x  PTT  x  x  x  x  v X>X X y< „ X  v  Toe  ..X  Ear < x  100  200  *  x><x x  300 Seconds  400  500  600  Pulse Transit Time to SBP SBP by PPT SBP by Cuff  100  200  300 Seconds  400  500  600  Figure 2.5: PTT and its BP Measured at Different Sites from a Single Subject.  2.5 P T T - B P Model Validation  21  and cuff-noninvasive systolic BP's were recorded every 5 seconds and the ECG and PPG were sampled at 300Hz from a Datex AS5 monitor. The parameter B in Equation (2.2.5) was kept constant for the entire case after it was estimated off-line once, as we intended to demonstrate the ability and robustness of the algorithm in capturing the changes in BP. The parameter B in Figure 2.6 is a constant offset and the variations in BP by PTT are not affected by the calibration. The cuff inflation produces the dips in the invasive BP due to the fact that the cuff was on the same arm as the arterial catheter. The case is over 100 minutes long and the cuff was cycled every 10 minutes. The results displayed in Figure 2.6 illustrate how effectively the algorithm can trace changes in BP. If the arterial BP signal is unavailable, the sudden increase or decrease of BP can be detected much earlier than using the cuff alone and appropriate procedures can be carried out to correct BP. One must note the constraint on the PTT-BP algorithm in that it can only be as accurate as its calibration standard, the cuff measurement. There is a -5mmHg bias between the cuff and the direct arterial measurement. As a result, the same bias can be observed between the PTT and the arterial measurement, though the mean difference between the two is less than O.OOlmmHg. Then 22 cases out of the 59 hypotensive cases described in Appendix A were processed. The ECG and PPG were sampled at 300Hz and the cuff sphygmomanometry systolic BP were recorded every 10 seconds. Even though the cuff BP was sampled every 10 seconds, the cuff inflates either every minute in normal mode or as frequently as it is able to in continuous inflation mode, fn normal mode, BP data remains unchanged until the next inflation. In continuous inflation mode, BP data recorded between inflations are interpolated by the monitor. During the operation of the study cases, spinal anesthesia and PE were administered, causing sudden but significantfluctuationsin BP. The 22 cases formed a total of 4657 data points where both the cuff systolic BP and the new noninvasive BP by PTT are available.  2.5 P T T - B P Model Validation  22  Pulse Transit Time to SBP i  >  i  i  i  —  140  —  SBP by PPT SBP by Cuff Invasive SBP  120 -  -  100 O) I E E  I i iff ; " ' V f | ^ A ^ .  80 r  i j 60 i  i i  i  i  I  i  V (i  i i i i  40 1000  !  1  iI 2000  i  iI I 3000 4000 Seconds  !  s iI  1  i i  i i  i  I  5000  1  i  i1 — 6000  Figure 2.6: Comparison between Systolic BP by PTT and Invasive Systolic BP from a Single Patient.  2.5 P T T - B P Model Validation  23  Figures 2.7 and 2.8 display the histogram of the difference and the limits of agreement [1] on the two BP measurements. Limits of agreement are often used by clinicians to compare the agreement between two measurements. It is a plot of the difference between the two measurements against the average of the two. Assuming a normal distribution, limits of agreement estimates the confidence interval from the standard deviation of the differences.  The mean difference is -0.0873 mmHg with a standard deviation of 11.49  mmHg. The offset parameter B was recursively calibrated by the total least squares, instead of being kept constant. Histogram of S B P by cuff minus S B P from PTT  2 5 0 0 r-  2000 1500h  1000r-  5 0 0 r-  -60  -40  -20  0  20  40  60  80  mmHg  Figure 2.7: Histogram on the Difference between Systolic BP by Cuff and by PTT from 22 Patients Undergoing Cesarean Section. The limits of agreement in Figure 2.7 can be misleading because, operating in the normal mode, the cuff only updates the BP reading immediately after an inflation. In continuous inflation mode, the monitor interpolates the systolic BP between cuff inflation, at each sampling interval. Figure 2.9 best describes the phenomenon. Before the spinal  2.5 PTT-BP Model Validation  24  Limits of Agreement on S B P by cuff and S B P by P T T  60  80  100 120 140 Avg of S B P by cuff and S B P by PTT (mmHg)  160  180  Figure 2.8: Limits of Agreement between Systolic BP by Cuff and by PTT from 22 Patients Undergoing Cesarean Section.  2.6 Conclusions  25  anesthetic is injected, the monitor was in normal mode, updating the BP every minute, so there is a delay between the two measurements. Cesarean Section Case — SBP by PTT - - SBP by Cuff O Spinal Anesthetic Phenylephrine  140  120h ixl  Q_ O O)  aat) £ a) ca  — 100  o  §  80  c 'CL  <n "o *  60  O)  X E E  40 500  1000  1500  2000  Seconds 20  Figure 2.9: Comparison between Systolic BP by PTT and BP by Cuff Sphygmomanometry during Cesarean Section (Case 10).  2.6  Conclusions  It is clear that the BP estimated by PTT is capable of early capturing the hypotension and hypertension caused by spinal anesthesia and PE, respectively. Continuous monitoring of BP will result in timely treatment of induced hypotension and avoidance of hypertension by careful titrated administration of PE. A simplistic, yet practical, model was developed to relate PTT and BP. With a more sophisticated calibration and noise filtering method, the real-time algorithm can easily be integrated into anesthesia monitors. Preliminary results suggest that the algorithm should  2.6 Conclusions  26  be an integral part of all anesthesia monitors as it requires no extra sensors to improve the quality and sampling rate of noninvasive BP measurement. The algorithm can also be utilized to trigger the recycling of the noninvasive cuff measurement, when a major BP fluctuation is detected.  Chapter 3  Patient Modeling The main challenge for identification of the patient model system is the lack of proper excitation. The measured disturbance, spinal anesthesia solution (SAS), is only injected once in each case. Since PE is used to reject the SAS disturbance, the patient is always in transient state when PE administration begins. Therefore no direct method can be used to separate the response of SAS and PE. Data collection procedures employed can be found in Appendix A. The patient model is assumed to be linear because of PE's dose-related increase in peripheral vascular resistance [13]. The model can be considered time invariant because of the short duration of the surgery (usually 30-60 minutes). No clear correlation was observed between the drug response and the physical properties of patients, such as actual weight, pre-pregnant weight, height, weight/height ratio and age. The problem was first handled as a multivariable system identification problem where PE, SAS and heart rate (HR) are inputs. The systolic blood pressure (SBP) is the output. This resulted in a model with a large degree of uncertainty. Therefore an attempt was made to estimate the response of SAS from the data of a non-hypotensive case, where no PE was administrated but the drop in BP was significant. The main assumption used is that each patient reacts the same way except for a different gain and time delay.  27  3.1 M I S O S y s t e m I d e n t i f i c a t i o n  28  In Section 3.1, the multivariable approach is discussed and essential results are presented. Then it is shown in Section 3.2 how a SISO system identification coupled with the SAS response estimation can improve the patient model by reducing uncertainty and determine that inter-patient variability primarily comes from the gain. These results were just published in [11] and [9].  3.1  MISO System Identification  This is the first attempt to analyze the effect on BP from PE, using black box modeling based on the collected data. The model performance was not satisfactory. Nonetheless this was essentially the first step to learn and understand the dynamics of the hypotension treatment. 3.1.1  MISO Method  Initially, as shown in [7], it was believed that BP is affected by the injection of P E and SAS and HR. Therefore the patient model was designed to have three inputs (PE, SAS, HR) and one output (BP), fntuitively three poles were chosen for the model to describe its dominant dynamics. A casual three-pole model can only have two or less zeros. The model was represented as a discrete transfer function matrix, see Figure 3.1. Two zeros and three poles were observed to be sufficient to approximate dynamics of PE and heart rate dynamics. For the SAS, an integrator and a delay of 18 samples were used to augment to the 2-zero-3-pole system. The integrator models the data trend and the delay accounts for the 180 seconds dead time in the SAS response.The sampling rate was set at 10 seconds. The 180 seconds dead time was an estimate drawn from observation and clinical data. Each input has specific zero location and gains but common dynamics (i.e. shares the same poles). The gain Ks exists only in the SAS input corresponding to the integrator gain. The HR input and SBP output are measured as the net value above patient's HR  3.1 M I S O System Identification  29  and SBP baseline; the SBP baseline is the SBP measurement of the patient prior to the surgery. If the delay in SAS input is considered external to the system model, the three discrete transfer functions can be converted into a 4-state MISO system. z +b21.z+b31 z -alz -a2z-a3 PE model 2  CDPE  3  2  Z +b22.Z+b32 z -alz -a2z-a3 SAS model 2  SAS  3  Dead time  Discrete Integrator  2  L  [HR B a s e l i n e ] [SBP B a s e l i n e ]  Figure 3.1: MISO Patient Model. Each case, see Figure A . l for an example, was identified simultaneously as a MISO system. The subspace method was chosen among others available in Matlab for the system identification because controller design is performed in state space. Since the subspace method, unlike other methods (e.g. autoregressive exogeneous input (ARX) and autoregressive moving average exogeneous input (ARMAX)) is non-iterative, a reasonable solution could be found even when the input excitation was limited to only an isolated excitation from the SAS and a few from the PE administration [14]. Even though this multivariable method was superior to conventional SISO identification, a number of cases were still discovered to have a negative gain for PE and positive gain for SAS, a physical impossibility. After further investigation, it was concluded that the oscillations coupled with the complex components in the identified poles and zeros were the cause of the incorrect modelling. In order to further refine the model, a structured ARMAX method was then employed  3.1 MISO System Identification  30  to parameterize the gains and the zeros but not the poles. The locations of the poles were first computed in average using the cases presenting real poles. This method was based on the results of the subspace method previously employed. The gain and the zeros location were estimated again using the ARMAX method. The gains and zeros were identified for each case. Their mean was derived based on the nominal model. The structured A R M A X approach helped identifying much fewer cases with complex zeros. The refinement technique successfully eliminated the number of outliers to three and ten out of sixty, in the PE channel and SAS channel, respectively. Outliers were discarded in the mean model estimation because of their complex components. Due to this fact the iterative search of the structured A R M A X method involved a certain degree of randomness. The results were slightly different every time but the variations were small compared to the analyzed uncertainties as discussed in the next section. 3.1.2  M I S O M o d e l and M o d e l Uncertainty  The model presented in this section was published in [11] derived from the first 40 cases of the 70 mentioned in Appendix A. At that time, the BP measurement relying on PTT was not yet developed and hence the BP represents the measurement taken by cuff sphygmomanometry. The nominal discrete transfer function model can be rearranged into a canonical state space model for uncertainty analysis, see Figure 3.2. The relationship between the state space and transfer function model parameters is demonstrated in the Equations (3.1.1), (3.1.2) and (3.1.3). For simplicity, the integrator augmenting the SAS channel is not included in the uncertainty analysis. Nonetheless the variance for the integrator gain, Ks, is 1000 times less than the mean. Therefore, neglecting this uncertainty for Ks was considered acceptable.  3.1 MISO System Identification  A=  31  " 1.2456 ' al 1 0 " 0.0919 a2 0 1 = a3 0 0 .  Ki K\b i Kihi  B  2  C=[  1  o  K Kb Kb  2  2  22  2  32  1 0 " 0 1 -0.3523 0 0 .  (3.1.1)  K Kb Kb  (3.1.2)  0.3893 -0.6213 0.0091 0.2842 -1.0247 0.0062 -0.0436 0.0901 -0.0016  3  3  23  3  33  (3.1.3)  0 ]  (3.1.4)  Ks = 0.0153  PE SAS HR  u  C  BjKK>  SBP  dB dAH  Figure 3.2: Linear Fractional Uncertainty Model with the Time Delay and Integrator Separated. As illustrated in Figure 3.2, the dA term contains uncertainty on the system poles and dB contains uncertainty on system zeros and gains. Since all three inputs share the same poles, only the first column of dA is nonzero. The distribution histograms for the twelve coefficients inside the A and B matrices are plotted in Figure 3.3, 3.4, 3.5 and 3.6, respectively. It is worthwhile to mention that the A matrix was not refined by the structured A R M A X identification and it represents an average of results from the subspace  32  3.1 M I S O S y s t e m I d e n t i f i c a t i o n  identification. Therefore, variables within the B matrix are more normally distributed and hence considered more accurate than the ones in the A matrix. The experiment was carried out to refine the A matrix with structured ARMAX when variables in the B matrix were fixed. The poor results during model validation procedure suggested that further refining the A matrix can degrade the model.  2r-  Distribution of a1  1.5 1 0.5-2.5  -2  -1  -1.5  -0.5  -0.4  0  -0.2  0  -1 Distribution of a2  0.5 Distribution of a3  0.2  0.4  -0.5  1  0  1.5  0.6  0.8  2  1  Figure 3.3: Distribution of First Column of A. According to common clinical practice, physicians often rely on the 95% confidence interval on dose effectiveness. If parameters are assumed to be normally distributed, the 95% confidence interval is between the mean plus or minus twice the standard deviation. The uncertainty representations dA and dB are determined from their corresponding variances and they are related to the linear fractional uncertainty model as the 95% confidence interval presented in Figure 3.2.  dA = ±  da\ da . da 2  3  0 0 0  o • 0 0 .  =±  " 1.288 0 2.173 0 . 0.935 0  0 " 0 0 .  3.1 MISO System Identification  33  Figure 3.4: Distribution of First Column of B.  Distribution of Spinal Anesthesia's K2  -10  -10  -10  -5  -8  0 5 Distribution of Spinal Anesthesia's K2b22  -6  -5  -4 -2 0 Distribution of Spinal Anesthesia's K2b32  0  5  10  10  15  2  4  15  Figure 3.5: Distribution of Second Column of B.  6  20  3.1 MISO System Identification  34  Distribution of HR's K3  -0.8  -0.6  -0.4  -0.2  Figure 3.6: Distribution of Third Column of B.  dB = ±  dbn db\i d6 db i dbii dbn 13  2  dbn  3.1.3  dbs2  dhs  6.485 4.129 5.923 7.299 5.510 8.270 0.554 0.766 0.555  (3.1.6)  R e s u l t s of the M I S O S y s t e m Identification  Simulations were run to validate the patient model. Most clinical data match the BP trend predicted by the MISO patient model, as in the left figure of Figure 3.7. The solid line is the simulated systolic BP and the dash line represents the clinical, data. In some cases, the model failed to reflect: 1. a constant offset; 2. an integrating factor in the PE; 3. estimate of the SAS gain {K ). 2  Some of these characteristics mentioned can be observed from the model validations results presented in Figures 3.9 to 3.11. An individual model was derived for the case displayed in the Figure 3.7 and its model validation is shown in Figure 3.8. The nominal  3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 35  model was comparable to the individual model in terms of capturing features of the overall response. The low fit percentage for the case-specific model led us to recast the system identification solution problem in a totally different manner. 160 140 m 120 100 80  model actual "  W  -  \ r  t /i  2500  1000  1500 Time (sec)  Figure 3.7: Case 14: Individual Model Validation.  3.2  Single Variable Approach Using Least Squares Estimation on the SAS Response  After the first MISO trial in Section 3.1, it was determined that the response of the PE must be separated from the SAS in a more effective way. However, the excitation of the SAS is followed closely by PE, due to the nature of this application. As there is always one and only one SAS injection in each Cesarean Section, the identification of SAS response is limited to at most one or two parameters. To compromise, an assumption was made that all patients have the same reaction to SAS except a different gain and time delay. Although the assumption seems coarse at first glance, the validation results prove its legitimacy. In the new SISO approach, the HR was dropped as one of the system inputs since its  3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 36  Figure 3.9: Case 17: Offset Mismatch during Model Validation.  3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 37  Figure 3.11: Case 12: Underestimate of K during Model Validation. 2  3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 38  effect was thought to be insignificant. Its components (i.e third column in the B matrix described in Equation (3.1.2)) are less than 10 times those representing PE and SAS. 3.2.1  SISO System Identification M e t h o d  To further reduce patient uncertainty, all data are treated as the fraction above or below the patient's baseline. Without subtracting the full response of SAS from data, the response of PE will be under-estimated if raw data is used for system identification. The first and most challenging step is to estimate the SAS response. The full SAS response can only be observed in non-hypotensive cases, where the gain is usually small and hence the low signal to noise ratio. However there is one case presented in Figure 3.12 where the patient's starting BP is high enough that even with a BP drop of roughly 30mmHg, PE treatment was not required. Using the data from that particular case, a discrete first order ARX model, augmented with an integrator and time delay is estimated for the SAS. The assumption made is that the SAS models for all patients share the same pole but use a different gain and time delay. With this assumption, the time delay is estimated experimentally based on a simple correlation analysis and the gain is estimated by a leastsquares fit between the impulse response of the SAS model and the data, as illustrated in Figure 3.13. The inapparent time delay often requires some clinical knowledge to be clearly identified. Although in a real-time application one needs a systematic estimating method as discussed in Section 4.4, tedious manual estimation is still preferred in the modeling stage. As a side note on potential protocol violation, according to Figure 3.12, the nonhypotensive patient seems to have a BP drop greater than 20% of the her baseline. However, her systolic BP baseline recorded outside of the operating room was in fact 106mmHg. Therefore all procedures of that case followed the designed and approved protocol - only giving PE when systolic BP is less than 20% of patient's baseline or lOOmHg, which ever  3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 39  the lowest. A Non-hypotensive C a s e  200  400  600  800  1000  1200  1400  1600  100 uj a. o  +  80  SAS PE  oi 60 oo » o  4 0  *  20  0  200  400  600  800 Seconds  1000  1200  1400  1600  Figure 3.12: Non-hypotensive Case (Case 8) with No Administration of PE, but Significant Decrease in BP. The PE response can then be extracted from data by subtracting the individual SAS response. Figure 3.14 displays the difference between raw data and data being used for PE system identification. The simplest model to capture essential properties of the PE response has two poles and a gain. All 22 cases with ECG and PPG properly recorded and PTT available were processed with the 22 ARX PE model. As shown in the pole map in Figure 3.15, the conjugate poles are localized around 0.97 ± i0.15 and most of them are stable. This justifies the assumption that all patients share the same pole in the SAS model estimation, although it may also simply indicate too fast a sampling rate. The imaginary poles in the PE model lead to unexplainable oscillatory PE impulse response. Since clinical knowledge suggests that PE does not cause any BP drop, even  3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 40  Least Square Fit of S A S Response  0 -0.05 -  C ccD o C O -' C Q C O C ^> -0 u1Xi  a.  o  •3 -0.15 -  u_  -0.2 -  -0.25 L  700  750  800  850  900  950  1000  Seconds  Figure 3.13: Least Squares Fit of the SAS Response. temporarily, a system identification forcing a positive system (real poles only) was developed (see Section 3.2.2).  SISO Model and Model Uncertainty  Although the inter-patient variability is small at the poles, it is not so for the gain. The DC gain varies between -0.0049 and 1.962, see Figure 3.16. The uncertainty in gain is significant even though it can be assumed that the gain for PE can never be negative. No clear correlation is observed between the PE gain and the physical properties of the patient, such as actual weight, pre-pregnant weight, height, weight/height ratio and age. It is concluded that a gain adaptation is necessary in the advisory system. Figure 3.17 summarizes the modeling of PE and SAS. A nominal model of PE, PE model(z), can be derived from the average poles and the gain:  3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 41  Separation of the PE and SAS Response 0.6 Fraction of SBP baseline data with SAS response subtracted Fraction of SBP baseline data 0.5h  _0.11 600  1  1  1  800  1000  1200  i  1  1400 1600 Seconds  1  1  1  1  1800  2000  2200  2400  Figure 3.14: Subtraction of the SAS Response From the raw data (same case as in Figure A.l). Pole-Zero Map  11  1  1  1  r-r^  0.8-  1  • - •. .i..  1  1  1  1  :  0.6-  0.4  /  -04 -  '•.  .  y -  -0.6-  :  -0.8.,1 -1  -  '• . 1 -0.8  1 -0.6  1—"••••I -0.4 -0.2  — I — 0  I- ••• "—I 0.2 0.4  1 0.6  1 0.8  Real Axis  Figure 3.15: Distribution of the Two Poles of PE.  1  3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 42  Figure 3.17: Patient Model in the Matlab Simulink Representation.  3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 43  PE model(z) = . ™ ' K  (3.2.1)  Z  (z-p)(z-p)  K  =  PE  0  2.795e~ , 4  <K  <  PE  p  1.6917e-  =  0.9414  (3.2.2)  3  0.9726 -MO. 1478,  < Re(p) < 0.9923,  (3.2.3)  0.0381 < Im(p) < 0.2083 The SAS ARX model is: SAS model(z) = v  .,  K  z (z-  ;  s  A  S  '  _  z  (3.2.4)  1)(Z- SAS)  k  P  Therefore, each patient has an individual gain and time delay, PSAS =  0.863  (3.2.5)  KSAS  -7.2951e" k  =  0 < k<  -2.4195e- , 4  = 4  <  KSAS  <  -2.7158e-  (3.2.6)  5  20.56, 56  (3.2.7)  Results of the SISO System Identification  The model is validated by comparing an open loop simulation with the actual data of a few cases, shown in Figures 3.18, 3.19, 3.20 and 3.21. Some of them are the cases which were validated by MISO modeling in the previous section. The manually estimated time delay was used in the model. Data used in system identification were employed in this model validation because extra data were not available. The model mismatch lies primarily in the SAS model. For example, the P E model would produce a much better fit in cases 14 and 19 of Figure 3.20 and 3.21 if the SAS model  3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 44  were more exact. From a control perspective, the disturbance model is less important and the mismatch can be compensated by a Kalman estimator [25]. The PE model fits the data well, given that a gain is estimated for each particular case. Using the nominal PE gain in the validation of case 12 of Figure 3.19 certainly produces a poor fit of data. The fit, shown in Figure 3.22, is improved if a customized gain is used. It is observed that a gain adaptation in the advisory system can further improve the performance of BP control. Model Validation  + 0  £ 80 ^ o oi zs  I  E  +  4  -  4  -  Model Output S B P data from PTT PE SAS  +  4  60<D  40 -  +  20 -  + 1  800  1  1000  I  1200  1  1  1400 1600 Seconds  I  l  I  1800  2000  2200  Figure 3.18: Case 10 PE Model Validation.  3.2.2  Positive SISO Models  While the model presented in Section demonstrates a good fit of data, the two imaginary poles of the system are introducing an oscillatory impulse response. Since the anesthesiologist does not expect that PE causes such an oscillation, one has to constrain the system's poles to be real. The constrained approach is similar to the conventional A R X least squares approach but it performs a linear search until the poles in Equation (3.2.1) become real. A similar  3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 45  Model Validation  40 k  800  900  1000  1100  1200  1300 1400 Seconds  1500  1600  1700  Figure 3.19: Case 12 PE Model Validation.  Model Validation  Figure 3.20: Case 14 PE Model Validation.  1800  3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 46  Model Validation  r~  "n  — + O  Model Output SBP data from PTT PE SAS  co  •s  £  e  60  40  1000  1200  1400  1600  1800 Seconds  2000  2200  2400  2600  Figure 3.21: Case 19 PE Model Validation.  Model Validation  n  r~ + O  Model Output S B P data from PTT PE SAS  120  E E  60  800  900  1000  1100  1200  1300 Seconds  1400  1 500  1600  1700  1800  Figure 3.22: Case 12 PE Model Validation with Customized Gain.  3.2 Single Variable Approach Using Least Squares Estimation on the S A S Response 47  positive system is also expected from compartmented models in pharmacokinetics.  PE model(z) = ,  K p  S  z  ,  =  (3.2.8)  £ Z PE  The solution of the quadratic solution determined that the system poles are real if a\ > 4 • 0,2- Therefore the modeling software searches for the maximum a that would 2  still result in two real poles. Figure 3.23 presents the two poles identified from the 22 cases. The Matlab implementation of this constrained system identification is included in Appendix D. 1  i  1  Pole-Zero Map  i  1  i  i 4x  i~~  First pole Second Pole  0.6 -  0.4 -  0.2 -  0-  x  x  xxxx  xacx+xf  x x * x-tf  x f + + -j* -t*f +  +-tHH-  -0.2 -  -0.4 -  -0.6 -  i  0.55  i 0.6  i  0.65  i  0.7  i  0.75  i  0.8  1  0.85  1—  0.9  Figure 3.23: Distribution of the Two Real Poles of Positive PE Model.  Positive SISO Model and Model Uncertainty  Model uncertainty is slightly increased using this approach but the model is more realistic. The uncertainties are laid out in a similar fashion as described in Section However, forcing the system to be positive sometimes results in an extremely small gain, which  3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 48  occurs when the identification uses a first order approximation (an almost horizontal line) instead of a higher order curve to fit the data. For that reason, gain distribution in Figure 3.24 is different from the one from the previous section and summarized in Figure 3.16. Therefore the average model uses the mean of the gain from only 6 cases with a DC gain greater than 0.0001.  K  =  PE  <K  0  PE  0.0022, (3.2.9)  <  0.1334  pi  =  0.8149,  pi  =  0.7156, (3.2.10)  0.6940 <pi<  0.8782,  0.5764 <p <  0.8461  2  An attempt was made to use the median poles and gain instead of the average values. The resulting poles and gain are: K~P~E = 0.0073,  (3.2.11)  pT = 0.8112, (3.2.12) pi = 0.7139, However, the two models show no noticeable difference during model validation. Therefore, the average model is kept as the nominal PE patient model.  Results of the Positive SISO System Identification  The same validation in Section 3.2 is repeated here in Figures 3.25, 3.26, 3.27 and 3.28. The performance of the positive models is slightly degraded but was more appreciated by our colleagues and anesthesiologists as more clinically relevant.  3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 49  Distribution of P E G a i n  0.02  0.04  0.06  0.08  0.1  0.12  0.14  Figure 3.24: Distribution of the DC Gain of PE.  Model Validation — Model Output - - S B P data from PTT + PE o SAS  801  20h 800  1000  1200  1400 1600 Seconds  1800  2000  Figure 3.25: Case 10 PE Model Validation.  2200  3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 50  Model Validation  + O  Model Output S B P data from PTT PE SAS  120  E E  60  20 800  900  1000  1100  1200  1300  1400  1500  1 600  1700  1800  Seconds  Figure 3.26: Case 12 PE Model Validation.  Model Validation — - + O  Model Output S B P data from PTT PE SAS  120  E E  60  +  +  20 1000  1200  1400  1600  1800  2000  2200  2400  Seconds  Figure 3.27: Case 14 PE Model Validation.  2600  51  3.3 Conclusions  Model Validation  140  +  o  Model Output S B P data from PTT PE SAS  ! 100  I  £  60  1000  1200  1400  1600  1800 Seconds  2000  2200  2400  2600  Figure 3.28: Case 19 PE Model Validation.  3.3  Conclusions  The main obstacle of the system identification problem is the separation of the input (PE) and disturbance (SAS) response since the excitations the patient can be subjected to are highly limited. The strict clinical protocol prohibits any additional artificial input for research purpose. Two system identification methods were proposed, the multivariable and SAS least squares estimation method. The model derived from the latter method performs a better fit on patient data. Three different patient models were identified. The oscillating SISO model has the best fit but it is believed that the positive SISO model is the most appropriate choice for BP control. Since there is no clinical explanation for the oscillation, a positive model is physically more plausible and thus would yield a more robust closed-loop. All three models have the gain as the main common source of model uncertainty. We believe that more data collection can further reduce model uncertainty, but a gain adaptation is essential  3.3 Conclusions  for real-time advisory control application of the patient model.  52  Chapter 4  Advisory System Design The advisory system in this Section determines the optimal PE dosage by applying modern control theory. The goal of the advisory system for use in British Columbia Women's Hospital is to reduce anesthesiologists' workload and improve patient care with minimum purchase of new equipment. Therefore, the system is designed to use two common pieces of hardware, a laptop and an electronic infusion pump. The laptop displays a userfriendly interface and drives the infusion pump for PE injection. The system is based on state-space model-predictive control framework, which was robustified to cope with large variation in the patient model, noise and disturbances. The system also has the ability to customize the PE dose for individual patient through online adaptation. The details of the programming logic can be found in Appendix B. Other technical features such as controller tuning, the Kalman filter design and online adaptation are to be discussed in this chapter. At the end of this chapter, the advisory system is tested with Matlab Simulink simulation environment to demonstrate its benefits and limitations. The system, installed on the laptop, is developed to run as a plugin under the Datex S5 Collect program. The anesthesiologists are familiar with S5 Collect software as it has been used for data collection in the hospital. To access the data from an existing Datex monitor in real-time, the system requires a Labview interface to connect to S5 Collect. The advisory system was coded in Labview 6.1, the programming language supported by the S5 Collect 53  54  software. The user interface (Figure 4.1) presents the amount of PE administered and the current systolic BP trend, using both the cuff and PTT measurement. For initialization, the system requires the systolic BP baseline of the patient, the target setpoint of BP, which can be changed at any time, and the input constraints (e.g. maximum infusion rate, bolus and total PE administered). As removing PE from the patient is obviously unfeasible, the minimum bolus or infusion rate is set to zero. Consequently, the system can only treat hypotension but not hypertension.  Pt SBP_regulator_firal_Vcrsi<m.vi  jvj  _* Plant Seconds  | Controls- Sectnos  j Adactaam  j sunuacwi  Advisory  j  S6Pf.Ptwakdars wstilrifeOffln  Information  CumtttTrw..  aeset rh(? ^yitem after modifying 'Ir^ronftohon' F* m the patients bese'ine C B D DDr  ObsawSBP  Log Path  C ^ f a l i n a D d i c l l M "  «e:\parryWC*_=ro1ett\ i sBP.-eaJatof'/eco-i.txl  Prnfinir-i. with K * Infaion  J«j*J  Controls SBP Set Point  Total Phenyteplrre Dose _  1  AdmnistBfed(LQ);|So3r 80.0 •:;3:10PW  Toggle svrftch  Phenylephrine Plot  -i A d v i i o r y Dose  Suggest  Phenylephrine  Dose(ug):  :  0  05:!8:00PH OS-.2JJ0OPM 0504:00PM 05:27:00PM C5:30:Q0W Tm i e  P h e n y e lp h n r ie Infuson R a e t  second  Automatic  Last Dose:|ioo,oo  Jf  5  10 07 5.s.o2.5-  Override  Phenylephrine  DoseCug):  g,  1  1  0(1|  o.oOSlUliCPM  o w a o o m c m i i o o p M m&mm T m i e  Phenylephrine Range  Adaptive IMIwI^I  Constraints Max Bolus rl 100.00  M a x Infusion Rate  ug  l-i 0.83  Max Total  vw*om  <B:».mm  ~~"jn)ercMi  Ij  learnt Mbce! Ihftiaf Hadet  ug/see  Phenylephrine  Enabte Gain Adaptation:  ii 1000.00 ug # o f New G a i n A d a p t e d  Disturbance ^A^smeM.joS  ''•  DeadTime. 16.00  i Adapt After  Pebad Irfltai Model  j ISO  Figure 4.1: Screen Shot of the User Interface of the Advisory System.  (Seconds)  4.1 Advisory System User Interface  4.1  55  Advisory System User Interface  There are five pages in the user interface: 'Plant Settings', 'Controller Settings', 'Adaptation', 'Advisory' and 'Simulation'. For real-time use, the user should go to the 'Advisory' page as shown in Figure 4.1: The user should ignore other pages once the controller is properly tuned. Since those pages contain information on patient dynamics and controller tuning, they should remain unchanged without expert knowledge. Plant Settings - The user chooses the number of input, output, disturbance and state. Controller Settings - The user specifies the patient state space model. Adaptation - The page displays the current internal state-space patient model. Advisory - The 'Advisory' page groups all essential controls and indicators compactly. Tables 4.1 and 4.1 summarize functions of the controls and indicators located on this page: Table 4.1: Controls SBP Setpoint Max Bolus Max Infusion Rate Max Infusion Rate Toggle Switch Go! Dose Override Phenylephrine Dose(ug) Override Enable Gain Adaptation: Reload Initial Model Adapt After  Controls on the Advisory Page of the System. The user chooses the desirable BP setpoint. The user selects the input constraint on maximum bolus. The user selects the input constraint on maximum infusion rate. The user selects the input constraint on maximum total amount of PE. The user changes the system between advisory and automatic mode. The user can administer the dose displayed in Suggest Phenylephrine Dose(ug) by pressing the button. The user can input the override dose. The user can administer the dose displayed in Override Phenylephrine Dose(ug) by pressing the button. The user can enable online adaptation by pressing the button. The system updates the gain of the patient model according to its observation. Adaptation is only carried in advisory mode. The user can load the original model and discard the learnt model by pressing the button. The system observes the response after each individual P E administration and updates the model after Adapt After seconds.  4.2 Controller Design and Tuning  56  Table 4.2: Indicators on the Advisory Page of the System. Indicators Suggest Phenylephrine The optimal PE dose suggested by the system. Dose(ug) SBP Plot The user sees the BP history as well as the future prediction. The red current time marker marks the current time and the prediction is plotted to its right. There are two different predictions, one assuming the suggested dose in Suggest Phenylephrine Dose(ug) is administered and one assuming no treatment. Phenylephrine Plot The user can review the PE administration history. In automatic (closed loop) mode, the plot also displays the PE dose prediction according to the most recent information and the internal patient model. # of New Gain The indicator displays the number of adaptations the system has Adapted performed so far. Simulation - Simulation can be enabled on the 'Simulation' page. Two text data files are needed to process the simulation. Each file contains the PE administration or spinal anesthetic administration information. The graphical user interface is easy to follow and the anesthesiologist is required to connect the electronic infusion pump, the laptop and the monitor together prior to use. The PE solution is prepared as usual but is injected through the infusion pump. The overall anesthesiologists' workload is dramatically reduced.  4.2  Controller Design and Tuning  There are three design criterions for the advisory system: 1) the system is model-based so that clinicians can verify the patient model employed; 2) PE dosage constraints has to be imposed; 3) In order to assist decision making, the system has to forecast future BP response. As the system has to be model-based and can handle constraints and make predictions, the MPC is a suitable choice. As MPC calculates its control moves by minimizing error over a prediction horizon, its structure satisfies those design criteria. A state space model predictive controller published by Maciejowski in [19] was implemented,  4.2 Controller Design and Tuning  57  following the confirmation that the patient model identified in Section 3.2.2 is controllable. The internal states are estimated by a Kalman filter and the control input is computed by minimizing the receding horizon cost function: Hp  Hu  (4.2.1) i=l The r) is the estimated BP observed by the Kalman filter described in Equation (4.2.4) and u is the input. The controller can be tuned by the cost function matrices Q and R as well as the predicting horizon Hp. Better performance can be achieved by increasing the weight in the Q matrix, at the cost of robustness. Increasing the weight in R detunes the controller, robustifying it. Increasing Hp enhances the stability and robustness at the expense of computational complexity. First, the patient model is arranged into time invariant discrete state space form: x(k + l)  =  y(k)  =  SBP(k)  =  y(k) • SBP baseline + SBP baseline  sampling interval  =  10 seconds  A  B C  Ax(k)  +  Bu(k)  Cx(k)+d(k)  0  1  1  (4.2.2)  -0.5832 1.5305 0 0.0022 1 0  In this SISO system, u(k) is the PE infusion rate in jUg/second and y(k) is the (SBPSBP baseline)/SBP  baseline at time k. d(k) is the output disturbance that the observer  aims to detect. Then a new augmented state is defined and the observer and controller are designed  4.2 Controller Design and Tuning  58  using an incremented input:  Ax(k)  =  x(k)-x(k-l) Ax(k)  C(k)  r,(k)  where r](k) =  y(k) — setpoint(k) is the observed output  y(k) =  A  0  B  CA  I  CB  (4.2.3)  Au(k)  S(k) + d(k)  0 I  A Kalman filter is then designed for state estimation:  A  0  CA  I  -L  0  5  /  Au(k) + Ly(k) (4.2.4)  The tuning strategy is discussed extensively in in [19]. The tuning is based on the frequency response (Figure 4.2) of the nominal model, using incremented input. With the open-loop cut-off frequency at about 6 x 1 0  -2  rad/s, the time constant is about 5 to 6  seconds in open loop. Thus the prediction horizon of the cost function (4.2.1) is designed to be about one minute (Hp=6), ten times as much as the open loop time constant. In other words, the controller minimizes the tracking error between chosen setpoint and actual B P output for the upcoming 60 seconds. The Q and R matrices in Equation (4.2.1) are set to 1000 and 1 for adequate performance of the controller, which is extremely important to this application. There is no side effect associated with a large control input except severe overdose, avoided by the hard input constraint. With the open-loop gain margin of the system increased from -20dB to the closed-loop 5dB, the phase margin from - 6 5 ° to 80°, the stability of the closed loop response is enhanced. The systolic BP output is conditioned by a conventional Kalman filter [25], before being used for feedback control. Its noise and disturbance covariance matrices are set to 1 and 1 respectively. The tuning is based on the occurrence of disturbance and noise are  4.2 Controller Design and Tuning  59  Frequency Response of Open and Closed Loop  — - Open Loop —  Closed Loop  50 40 -  *  20  r  -10  -  -20  -  -30  0F= -45  3  -  -90  -  "135  -  ~i ~ ' m  a  -225  -  -270  -  -315  -  L  -360 10~  3  10"  10"  2  1  Frequency (rad/sec)  Figure 4.2: Frequency Response of the Open Loop and Closed Loop Nominal PE Model with Tuning Q=1000, R = l , Hu=l, Hp=6. of similar magnitude. The main source of uncorrelated disturbance and noise is expected to be body movements. Since this is an advisory system, the control input has to be approved by the anesthesiologist based on clinical knowledge. However, the anesthesiologist cannot examine each dose at each sampling interval. Therefore, at all times, the system displays the suggested control input. The system predicts the future according to the current state and calculates an open loop dose that brings the BP back to the setpoint. The anesthesiologist can either follow the suggested dosage, or administer his/her own discretionary dosage through the advisory system, as discussed in Section 4.1. All PE administration is logged for the anesthesiologist's record. To help justifying the decision, the system predicts the future systolic BP output for two scenarios: 1) when the suggested dose is injected; 2) when nothing is injected. The system updates its prediction whenever a new BP measurement arrives. The gain of the internal PE model is recursively updated, which will be discussed  4.3 Safety Concerns  60  in details in Section 4.4.  4.3  Safety Concerns  There are two safety issues related to the functionality of this system: PE overdose and hypotension. Overdosing of PE is avoided by a hard input constraint, lOO^tg of PE, and the anesthesiologist's direct decision. When the system predicts hypotension, a continuous flashing of the 'drug injection' button on the user interface serves as a silent alarm to the anesthesiologist. In cases of system failure, the anesthesiologist in-charged could disconnect the infusion pump and administer PE in the usual way. There is a possibility that the BP index derived in Chapter 2 is unavailable due to ECG or PPG artifacts. In that case, the system tries to combine the previous and current cuff BP measurement with the previous BP estimation. The details on estimating the BP for control are laid out in Figure 4.3. The most undesirable case is when both BP signals are unavailable. Then due to incomplete information the advisory system will stop advising and will display a warning message on the screen. At that point no control action is suggested by the system but the anesthesiologist can still decide a dose based on patient history and other physiological data.  4.4  P E Gain Adaptation and SAS Time Delay Estimation  According to the results on patient modeling in Section 3, model uncertainty exists mostly in the PE gain. Consequently online adaptation must be performed. Since changes in the input are in an impulse form, conventional least squares techniques utilizing a covariance matrix will not deliver good results. Alternatively, the system only examines the peak response for each PE dose and compares it with the peak of the current internal model. The gain adaptation of PE faces the same challenge as during the system identification phase - separating the PE and SAS response. While the same assumption as in Section 3.2  4.4 P E Gain Adaptation and SAS Time Delay Estimation  If PTT S B P available  If previous PTT S B P available  If both previous and current cuff S B P available  61  Use PTT S B P  If cuff SBP available  Use cuff S B P  S B P not available Stop Advising  T •Use.-. • . ::;•:] • previous PTT S B P + (cuff S B P - previous cuff SBP)  Figure 4.3: Combining B P measured by cuff with BP measured by P T T . can be made and the same technique can be employed, the time delay of the SAS can no longer be determined by expert knowledge in a real-time application. Therefore, the Variable Regression Estimation (VRE) [6] is used to compute the time delay. The V R E method was selected because of its simplicity. As the system must complete the time delay estimation before next B P reading arrives from the monitor, the computation has to converge rapidly. The V R E algorithm calculates the cross correlation function between the SAS input and the B P output increment at each sampling time when the first dose of P E is administered. The details of V R E are included in Appendix C. The SAS time delay corresponding to the maximum cross correlation function is chosen. Then the gain for SAS is again estimated by least squares and the estimated SAS response is subtracted from the adaptation data. Table 4.4 lists the time delay estimation of the 22 cases used for system identification in Section 3.2 by expert knowledge and V R E .  4.4 P E G a i n A d a p t a t i o n a n d S A S T i m e Delay E s t i m a t i o n  62  Table 4.3: Comparison of SAS Time Delay Estimation by V R E and Expert Knowledge. Case number 10 12 14 17 19 20 23 26 36 37 45 46 47 49 50 52 53 54 56 60 63 64  Delay by Expert Knowledge (10 Sec) 4 18 18 25 15 36 10 38 16 4 24 29 56 15 37 18 16 48 18 10 10 0  Delay by V R E Estimation (10 Sec) 4 27 9 26 • 8 36 6 50 14 5 26 36 4 4 12 7 17 4 9 14 14 5  It can be observed that in three particular cases, 47 and 50 and 54, the V R E estimation is significantly different from that suggested by expert knowledge. One would find that sometimes the time delay cannot be clearly identified, especially in the case 47 displayed in Figure 4.4. The body movement following the SAS injection produces a lot of disturbance in the system, so the time delay estimation can sometimes be difficult, even for the anesthesiologist.  For the case 54 (Figure 4.5), the problem can be eliminated if  the Kalman filtered output is used for V R E instead of the actual output. The sharp dip of B P is caused by a temporary disconnection of sensor when the patient laid down after SAS injection. In case 50 (Figure ??) the sensor was unplugged for an even longer period causing a loss of BP data. This problem is solved for the advisory system by the scheme that fills in missing BP measurement described in Figure 4.3. Other than the three cases, Table 4.4 suggests that V R E is a simple yet effective solution for time delay estimation  4.4 P E Gain Adaptation and SAS Time Delay Estimation  63  given the lack of system excitation. Cesarean Section Case O +  140  SBP by PTT case47 Spinal Anesthetic Phenylephrine  120  1500 Seconds  2000  2500  Figure 4.4: Case 47 where V R E Estimation Does Not Agree with Expert Knowledge.  case50.eps.7Case 50 where V R E Estimation Does Not Agree with Expert Knowledge. The response of each P E administration is used for adaptation, after being filtered by the Kalman filter designed based on Equation (4.2.4). The simple adaptation scheme in advisory mode is presented in Figure 4.6. The peak response is compared with the peak response of the internal model. The gain for the adaptation is calculated based on the fixed model structure. The new gain of the internal model is then computed as the average of the adaptation data gain and the previously stored gain. Computing the new gain from both the current model and new observation would make the system less vulnerable to noise and artifacts. The adaptation updates the internal model of the controller, which consequently updates its control gain. However, to avoid interference, the adaptation will not take place when a P E administration is closely followed by another one.  4.4 P E Gain Adaptation and SAS Time Delay Estimation  64  Cesarean Section Case  Figure 4.5: Case 54 where V R E Estimation Does Not Agree with Expert Knowledge.  Data BP  Adapt  New Model Impulse Response  \ o i d Model Impulse Response Time 300 Seconds  Figure 4.6: Illustration on the Advisory System Adaptation Scheme.  4.5 Advisory System Simulation  4.5  65  Advisory System Simulation  The advisory system was tested in a Matlab Simulink framework before the Labview embedded system implementation was pursued. Figures 4.7 to 4.14 demonstrate the performance of the advisory system, both in advisory and closed loop mode, under different circumstances, such as the existence of noise and model mismatch. In advisory mode, the simulator arbitrary assumes that the anesthesiologist selects a setpoint of lOOmmHg and administers the suggested dose every 100 seconds. Model Simulation —  Model Output Observed Output + PE o Spinal Set Point I  120  ~ ~ \  UJ  80  40  +  + + + + + + + + + + + + + + + + + + + +  20  0  500  1000  1500  2000  2500  Seconds  Figure 4.7: Perfect Model Advisory System Simulation. When the system has the perfect model, it performs flawlessly both in guided open loop and closed loop (Figures 4.7 and 4.11). Continuous PE infusion in the closed-loop simulation eliminates the undesirable BP oscillation shown in the advisory mode simulation. The system can also deal with measurement white noise (Figures 4.8 and 4.12). As shown in Figure 4.9, when there is a gain mismatch, the gain adaptor quickly learns the actual dynamics of the patient and revises the internal model accordingly. The gain  4.5 Advisory System Simulation  66  Mode! Simulation  —  Model Output Observed Output  + PE o  Figure 4.8: Perfect Model with Noise a  2  —  Spinal Set Point  lOmmHg Advisory System Simulation.  Model Simulation  1000  1500 Seconds  Figure 4.9: Doubled Gain Model Advisory System Simulation.  4.5 Advisory System Simulation  67  Figure 4.10: Oscillating Model in Section Advisory System Simulation.  Model Simulation  1000  1500  Seconds  Figure 4.11: Perfect Model Closed Loop Simulation.  4.5 Advisory System Simulation  68  Mode! Simulation  140 + o —  Model Output Observed Output PE Spinal Set Point  ^ ^ ^ ^  20  ++  + X  + +  + -it-  +  ++ + + +  ++ +  +  + -*fH-  H- #  +  1500  • + +  ++ + ^ 2000  Figure 4.12: Perfect Model with Noise a = lOmmHg Closed Loop Simulation. 2  Model Simulation  Figure 4.13: Doubled Gain Model Closed Loop Simulation.  4.5 Advisory System Simulation  69  Model Simulation - + O —  140  Model Output Observed Output PE Spinal Set Point  ! 100  + ra 60  4-  + +  S +  A t  + +  +  +  +  -h  +  + +  1000  1600  Seconds  Figure 4.14: Oscillating Model in Section Closed Loop Simulation. adaptor has not been developed for the closed loop mode and hence no improvement can be observed at the end of the closed loop simulation, presented in Figure 4.13. When the patient's BP oscillates as described in Section, the system can become unstable if operated in closed loop, see Figure 4.14. This phenomenon is not surprising. Another simulation was carried out to show that the M P C with a perfect model can control the oscillatory patient model. However, the M P C controller with a doubled gain mismatch already drives the oscillatory patient model to become unstable. Closing the loop in a oscillatory plant with slight model mismatch can easily lead to instability; this oscillatory patient model is no exception. Furthermore, according to the anesthesiologist's advice, oscillation does not exist and that is the reason for the positive model derivation in Section 3.2.2. Investigation on the oscillation is needed if the anesthesiologist decides to close the loop. As an advisory system, it is ready for clinical verification.  4.6 Conclusions  4.6  70  Conclusions  Although at its current stage the system is literally a guided open loop application, it has the potential for closed-loop predictive control. Once anesthesiologists test and gain confidence in the system, they can monitor the process in a closed loop manner. According to the patient model, PE is a fast response drug. A dose of PE is required for almost every minute if BP has to be maintained during hypotension. The ability to run the system in fully automatic mode is potentially an improvement in patient care.  Chapter 5  Conclusions 5.1  Main Contributions  Advisory system for Cesarean Section hypotension treatment - An advisory sys-  tem for hypotension treatment in Cesarean Section was developed. The advisory system is expected to be an elegant solution for hypotension treatment. As the advisory system has knowledge of the drug's transient response, it is able at all times to calculate the necessary PE dosage for different circumstances. Furthermore, the gain adaptation can ensure that the dose is tailored to each individual patient. When the drug delivery is handled by reliable electronics, consistency and gradual change in dosage are enforced. The safety aspects related to the human subjects are dealt by constraints implemented within the advisory system. The advisory system should reduce the workload of the anesthesiologist by introducing a guideline for PE administration, while improving patient care. Continuous noninvasive blood pressure measurement - The poor quality of the  noninvasive blood pressure measurement using the noninvasive blood pressure cuff introduces an infrequent sampling problem. This problem was solved by data conditioning using pulse transit time. The continuous blood pressure reading makes the locations of the pole of the patient model obvious. The model gain is the main source of inter-patient  71  5.2 Limitations of the Advisory System  72  variability. The algorithm itself can be employed in general patient monitoring independently from the advisory system. As the cuff sphygmomanometry, electrocardiogram and plethysmogram are commonly available, blood pressure can be measured continuously at nominal extra cost. Patient model for phenylephrine and spinal anesthetic - The blood pressure  response to phenylephrine following spinal anesthesia was analyzed systematically. This achievement marks a new page in pharmacokinetics and pharmacodynamics for phenylephrine. As a result, the drug description is no longer limited to plasma half life and ED95 (Effective Dose at 95% confidence interval). Useful information contained in the proposed state space patient model can assist researchers in future studies using phenylephrine following spinal anesthesia.  5.2  Limitations of the Advisory System  Although a closed-loop predictive controller was designed, the patient model was based on bolus injection of phenylephrine. The actual behaviour of continuous phenylephrine infusion is still unknown. A new infusion phenylephrine model may be necessary if the system is to be operated in closed loop with continuous drug infusion. The gain for the phenylephrine patient model varies significantly between patients. Currently, there is no systematic method to predict the gain based on patient physical properties. It was observed from the data obtained through the approved clinical study that actual weight, pre-pregnant weight, height, weight/height ratio and age are not correlated with the gain. Therefore in real-time applications, the gain is customized to individual patients through online adaptation. As a result, the first few doses suggested by the system largely rely on the nominal model, which could be quite different from the actual unknown model. Nevertheless, the proposed system does not introduce this inherent problem of guessing the phenylephrine response before the administration, a task  5.3 Direction for Future Research  73  anesthesiologists have tackled in every case.  5.3  Direction for Future Research  Patient variability is a problem in biomedical engineering in general, especially in the hypotension treatment problem. The dynamics of the proposed patient model is not complex yet uncertainty is substantial. Future studies may refine the patient model by additional data and investigation of the relationship between model parameters and the patient properties. Once the actual patient model is identified, the control of blood pressure becomes relatively straightforward based on the advisory system proposed. Measuring the blood pressure noninvasively and continuously is also a relevant future research topic. It is believed that improvements can be made to both the mathematical model and the peak detection algorithm. For example, parameters like arterial cross section area and pre-ejection period can be assessed through MRI or ultrasound technology. Moreover, there could be other techniques of noninvasive blood pressure measurement such as the pulse pressure evaluation. Better blood pressure monitor can easily be integrated into the advisory system with potential for enhanced performance. Lastly, the advisory system has to be tested for its performance, safety and consistency. Upcoming clinical trials of the advisory system should confirm the value of this project in the near future.  Appendix A  Data Acquisition Clinical data collection was the first essential step for patient modeling. For ethical reasons, researchers had to avoid alternation to current medical procedures and observe the case passively. No drug was administered for research purpose. With the help of computers and electronics, subtle changes can be detected and analyzed by reviewing physiological data saved in digital format. Compared to traditional data recording by hand writing, the automated data collection improves accuracy and reduces human errors.  A.l  Procedures  Following institutional approval and informed consent, data was collected from 70 randomly selected patients undergoing Cesarean Section spinal anesthesia at the British Columbia Women's Hospital. Patients with significant medical disease, pregnancy complications or extreme weight and height were excluded. All the patients were preloaded with normal saline, 30 ml/kg (pre-pregnant weight) within 30 minutes prior to the administration of spinal anesthesia. Administration of spinal anesthetic was standardized for all patients for dose, position and intervertebral space. The standardized SAS consisted of hyperbaric bupivacaine 0.75%, preservative free morphine 200/xg and fentanyl 15/xg. Total volume of the solution was 2.0ml. Each episode of hypotension, defined as systolic BP < lOOmmHg or 80% of the patient's 74  A.2 Conclusions  75  baseline BP recorded prior to the surgery, which ever was lowest, was treated with PE during the time period between the injection of spinal anesthetic and the birth of infant. Fifty nine cases out of the 70 (84%) were hypotensive. The dosage of PE depended on a random grouping of the patients. There are 4 groups, 20-40-60-80/ig. However, if it was observed that the dose was not strong enough to treat hypotension, a doubled dose of up to 100/ig was injected. The study follows a strict protocol approved by the ethical review board of the British Columbia Women's Hospital. Noninvasive BP and HR measurements along with other available physiological signals were recorded every 10 seconds from a Datex AS5 monitor (Datex, Finland), connected via a serial port to a Dell laptop running Datex S5 Collect software. Even though the cuff BP was sampled every 10 seconds, the cuff inflates either every minute in normal mode or as frequent as it is able to in continuous inflation mode. In normal mode, BP data remains unchanged until the next inflation. In continuous inflation mode, BP data recorded between inflations come from interpolation performed by the monitor internal software. The monitor operated in continuous mode after the injection of spinal anesthetic. The ECG and PPG are recorded by the same software with a sampling rate of 300Hz. The administrations of drug or special events were marked with the snapshot function of the monitor and hence are synchronized with the data. Physical properties such as pre-pregnant weight, pregnant weight, height, SBP baseline and HR baseline were also recorded before the surgery in a patient report.  A.2  Conclusions  Despite the 10-second sampling rate, the minimum repeat inflation time for the cuff was roughly 35 seconds. Data points between cuff inflations come from an interpolation performed by the monitor. Such data are observed to be less reliable and poorly conditioned, as indicated in Figure A . l . A more frequent BP measurement was essential for modeling  76  A.2 Conclusions  purpose and achieved by the PTT method, developed and described in Chapter 2 A Typical Case 140 120 X  E E  100 80 60 0  100 r of "5 80 -  SAS -l-  PE  60 w  < 40 20 500  1000  1000  1500  2500  Seconds  Figure A . l : Case 10. A Typical Case. The database of 70 patients was established. While the data was employed for the patient modeling, it can also be used for other studies.  Appendix B  Advisory System Structure Overview The software of the system is composed of two main parts, the Datex S5 Collect interface and the SBP regulator. The interface is responsible to communicate with S5 Collect, retrieve physiological data, carry out the PTT algorithm described in Chapter 2 and compute a BP index for control purpose. After the user specifies the patient dynamics, constraints, patient baseline and Kalman filter / controller tuning in the initialization, the SBP regulator interacts with the user and begins BP control. The details of the programming logic are illustrated in Figures B.l and B.2.  B.l  SBP Regulator Interface  SBP Regulator Interface is an interface between the SBP Regulator and the Datex  S 5 Collect. The interface is built to follow Datex protocol and acquire real-time data from Datex S5 Collect. With the help of other subroutines, it produces a BP index for BP control in SBP Regulator. Datex S 5 Collect - is the data collection software developed by Datex-Ohmeda. It connects with Datex anesthesia monitor via a serial port and retrieves data in real-time. Find P T T - detects the PTT from R peaks calculated by E C G to R R I and differentiated PPG peaks calculated by Find Pleth Peaks. The resulting PTT is used by the 77  B . l S B P R e g u l a t o r Interface  78  SBP Regulator Interface E C G to RRI  Datex-Ohmeda S5 Collect  ECG,  Wavelet SWTdec  | ,  FGfi S-levfll SWT  Peaks  ECG Detect R Time, Physiological Signals  R Peaks, RRI N Find PTT  RRI  PTT  Compute SBP by PTT  Sample SBP by cuff Recalibrate  Intervals  Find Pleth Peaks Differentiated G ^ p p  Flip Conv PPG  3  Calibrated Offset  Differentiated; PPG Peaks Differentiated PPG Peaks  SBP for Control  SWT l_» Peakdetect [Approximation  SBP Controlled Variable  SBP Controlled Variable  SBP Regulator  Phenylephrine Infusion Rate  H  Grasby 4300 Infusion Pump  Phenylephrine  Patient  Figure B . l : Advisory System Data Flow Diagram 1.  B . l SBP Regulator Interface  79  S B P Regulator Adaptation / Initialization CORE:  Kalman_design  Advisory  Advisory Input  Plant Dynamics  Get User  Plant Dynamics Setpoint, Internal States  Observer  1  Input and B P Index. Plant Dynamics  mpc_design  Compute Phenylephrine Dose  Advisory Dose  Feedback Gain  Plant Dynamics Impulse Response! Plant Dynamics,  Prediction Plant Dynamics, Internal States, Advisory Dose  Detrended S B P  Dose Prediction Display Plant Dynamics, Internal States, Feedback Gain  Peia^  Adapt_SBP  New Plant Dynamics  Prediction  Model .Prediction Display Continuous  File Path Phenylephrine Dose  Simulation ReadDose From File Disturbance  Figure B.2: Advisory System Data Flow Diagram 2.  ReadSpinal FromFile  B.2 SBP Regulator  80  SBP Regulator Interface to derive the noninvasive BP index discussed in Chapter 2. E C G to R R I - Given the ECG, it decomposes the signal into 5 level SWT using the WaveletSWTdec (see Section 2.3). Then a rule-based system in E C G Detect R carries out the R peak detection in ECG. With the computed R peaks, Sample Intervals interpolates a R-R interval time series for PTT computation. Find Pleth Peaks - discretely differentiates the PPG and filters noise by 1 level SWT (completed by Flip Conv). Then a rule-based system in Peakdetect detects the differentiated PPG peaks for later PTT computation. Recalibrate - calibrates the offset parameter B described in Section 2.2 using total least squares. The algorithm minimizes the discrepancy between the BP measured by PTT and by cuff. SBP for Control - combines the BP measured by PTT and by cuff to produce a BP index for control purpose. The details can be found in Section 4.3.  B.2  S B P Regulator  - S B P Regulator operates in two modes, the advisory and automatic (closed loop) mode, fn advisory mode, the system displays a PE dose for hypotension treatment as a suggestion, which the user can follow or override. In automatic mode, the system regulates patient's BP continuously by PE infusion. The BP and dose history along with future prediction are displayed graphically for reference. The user first specifies the dynamics of the plant in state space form (A, B and C matrices), the input constraints, tuning of controller (Q, R, Hu and Hp) and Kalman filter (Q and R) and the patient's BP baseline (see Section 4.2). The system then computes a feedback gain for the patient and is ready to run. Throughout the process, the system is able to adapt to the actual patient dynamics and update the internal model in the MPC if adaptation is enabled. Kalman_design - designs an observer described in Section 4.2 according to the tuning  B.2 SBP Regulator  81  of noise and disturbance covariance matrix specified by the user. mpc_design - designs a feedback gain described in 4.2 according to user's input of the cost function matrix (Q and R) and the prediction and control horizon (Hp and Hu). Model Display - computes the impulse response of a the patient's PE model. It is used to display the current adapted model as well as the nominal model. Adapt_SBP - observes the PE drug response whenever a dose is administered. Then it updates the gain of the MPC internal patient model, which in turn updates the feedback gain. The adaptation algorithm is disclosed in Section 4.4. Advisory Input - computes the advisory PE dose from user defined SBP setpoint. It forecasts future SBP trend from the estimated states and derive a necessary dose from the MPC internal patient model that is just enough to drive BP back to its setpoint. Dose Prediction Display - In advisory mode, it forecasts the future BP trend for two situation, 1) when the advisory PE dose is administered 2) when no action is taken. The information is plotted in a graph for easy interpretation. Model Display Continuous - In automatic mode, it predicts the amount of PE that is expected to be administered and the resulting BP response according to the most recent information and the internal patient model. The information is displayed graphically for user's reference. ReadDoseFromFile and ReadSpinalFromFile - They are for simulation use only. They read the PE dose and spinal injection from a file and notifies the system about the administration. The simulation mode should be disabled for real-time use.  Appendix C  Variable Regression Estimation V R E calculates the cross correlation function between the input and the output at each sampling time. The time delay corresponding to the maximum cross correlation function is chosen. No priori knowledge is required for the estimation. If the system is described by the A R M A X model with a delay k:  y(t)  9  =  $(t)9  =  [-y{t-l),...,-y(t-n),u(t-k-l),...,u(t-k-m)}  =  (C.0.1)  [a ,...,a ,bi,...,b ]  T  1  n  m  The cross correlation between input and the output increments is given by:  Ei{tM) 1<  = ki <  iE5=o«(3-*i)(y(j' + i ) - » ( j ) )  (  C  0  2  )  time of first PE dose — time of SAS  A different cross-correlation E\ is computed from each guess offc;.The estimated time delay is the ki that corresponds to the maximum of E\. In this particular problem of estimating time delay of spinal anesthetic, the delay k is expected to range from 1 to the time of first phenylephrine administration minus time of spinal anesthetic injection.  82  Appendix D  Matlab Code for Positive System Identification p r o c e s s d a t a . m  D.l  num=64; %case number removespinal; upsampleX = 1 ; PEmodel = []; for i=l:length(PEs)/upsampleX tempPE = 0; for j=l:upsampleX if PEs(j+(i-l)*upsampleX)i,0 tempPE = PEs(j+(i-l)*upsampleX); end end PEmodel = [PEmodel; tempPE]; end BPmodel = detrend(resample(BPs, 1, upsampleX)); BPmodel = BPmodel(l:length(PEmodel)); Ts = 10*upsampleX;  83  D . l processdata.m  84  orderA = 2; %orderA^,=orderB orderB = 1; [filtB, filtA] = butter(l, .007/.05, 'low'); Hq = sos(qfilt('referencecoerncients', filtB, filtA)); BPfilt = filter(Hq, BPmodel); [ufiltB, ufiltA] = butter(l, .025/.05, 'low'); uHq = sos(qfilt('referencecoefficients', ufiltB, ufiltA)); PEfilt = PEmodel; Y = BPfilt(orderA+l:end); % leave orderA amount of history for estimation N = length(Y); X = zeros(N, orderA+orderB); for i=l:N X(i,:) = [-fliplr(BPfilt(i:i+orderA-l)'), fliplr(PEfilt(i+orderA-orderB:i+orderA-l)')/Ts] end A = D ;  b=D;  i=0; while ij50 lb = [-inf; -.5; 0; 0]; ub = [inf; (l-.02*i); inf; inf];  % A = [0 1; % 1 0];  % % b = [0;0];  theta = lsqlin(X,Y, A, b, [], [], lb, ub); PEsys = tf([theta(orderA+l:end)'], [1 theta(l:orderA)'], Ts);  D.2 removespinal.m  poles = pole(PEsys); if imag(poles(l))==0 i=50; else i=i+l; end end simBP = 80/Ts*impulse(PEsys, (length(BPmodel)-l)*Ts); D.2  r e m o v e s p i n a l . m  dtsp = dead_time(num); cmd = sprintf('BP = transpose(SBP_PTT%i);',num); eval(cmd); cmd = sprintf('HR = data%i(:,l);',num); eval(cmd); cmd = sprintf('Spinal = data%i(:,3);',num); eval(cmd); cmd = sprintf('Phenylephrine = data%i(:,4);',num); eval(cmd); cmd = sprintf('DBP = data%i(:,5);',num); eval(cmd); i=l;  SpinalStart=l; PhenyStart=size(Phenylephrine,l); dBP = [0; diff(BP)]; while 1  85  86  D.2 removespinal.m  if Spinal(i,l) =0 SpinalStart = i; break; end i=i+l; end if Hypo(num) =0 %Hypo occurs i=l; while 1 if Phenylephrine(i,l) =0 PhenyStart = i; break; end i=i+l; end else PhenyStart = size(Spinal,l); end %  Spinal Estimation  BPsp = BP(SpinalStart-3+dtsp:end)-baseline(num)/baseline(num); Spinalsp = Spinal(SpinalStart-3:end-dtsp); %  End of Spinal Estimation  y = 100*impulse(SPmodel(l), (length(BP)-SpinalStart-dtsp)*10); %from Spinal to end % fit y with BP(SpinalStart+dtsp:PhenyStart) % K*y = BPs  D.2 removespinal.m  temp = (BP(SpinalStart+dtsp:PhenyStart-l)-baseline(num))/baseline(num); K=y(l:PhenyStart-SpinalStart-dtsp)\(temp-temp(l)) if KiO BPs = (BP(SpinalStart+dtsp:end)-baseline(num))/baseline(num) - K*y; else BPs — (BP(SpinalStart-|-dtsp:end)-baseline(num))/baselme(num); end BPso = (BP(SpinalStart+dtsp:end)-baseline(num))/baseline(num); PEs = Phenylephrine(SpinalStart+dtsp:end); HRs = HR(SpinalStart+dtsp:end)-HRbaseline(num);  87  Appendix E  Matlab Code for Noninvasive BP Measurement by P T T E.l PTT2BP.m SBP_freq=10; %SBP sampling frequency freq = 300; %ECG and Pleth data collected at 300Hz [Rpeak, errors, avgs] = rdetect (ecg, freq); [dpleth.peaks, dpleth_peak_amp, dpleth] = plethdetect(pleth, freq); [PTT, PTT_x] = findPTT(Rpeak, dpleth_peaks, HRbaseline(casenum)); [SBP_PTT_raw, SBP.X, B, A] = findSBP(PTT, P T T J C , filt_BP, DBP, SBPJreq, height(casenum))  SBP.PTT = SBP_PTT_raw + B; E.2  r d e t e c t . m  function [Rpeak, errors, avgs] = rdetect (ecg, fs) % detect R-peaks of an ECG signal, fs is the sampling frequency and ecg is the electrocardiogram signal plotenable = 0; upsamplerate = 5; interpolateenable = 1; Rpeak = []; errors = []; avgs = [];  88  E.2 r d e t e c t . m  % loop through wave form in 512 sample chunks, this should % be big enough to always capture at least one beat len = 512; % this defines the overlap at the end of the sample, no peaks are detected within % this range to prevent strange things happening at the edges. edge_buffer = 16; % blank out range, the heart cannot beat within 200msec of the previous beat blankout = 0.200 * fs; % range to look for modulus minimum prior to R-peak. detection of this will % differentiate from downward spikes % it should occur within 0.03 of a second before R-peak mod_min_range = round(0.03 * fs); % initial value for standard dev and peak amplitude for running average filtering coeLstd = 200; peak_ampl = 1000; wavename = 'dbl'; levelqstart = 3; levelrpeak = 3; leveltstart = 5; maxlevel = max([leveltstart, levelrpeak, levelqstart]) ; % start at the first non-zero value segmentstart —find(ecg);segmentstart = segmentstart(l); % initialize variables peakblockindex = []; peakfound = 0; catchtime = 188 * 300 * inf; h = waitbar(0, 'Locating R peaks...', 'Position', [100 100 270 50]); while segmentstart + len - 1 j length(ecg) % take a sample of the larger signal s = ecg(segmentstart:segmentstart-|-len-l);  89  E.2 rdetect.m  90  if segmentstart i catchtime pic-tenable = 1; end if plotenable figure(12)  plot(s) end % decompose it using stationary wavelet transform [swa,swd] = swt(s,maxlevel,wavename); % look for R peak % extract the coefficients at the desired level coef = swd(levelrpeak,:); % calculate stdev of wavelet coefficients, weight with past history to % filter large devatiations from spikes noise_warning = 0; noise_energy = 0; newstd = std(coef); % if potential noise problem do extra checking to cancel operation if it is a bad segment if newstd ^ 1.5 * coeLstd — newstd j 0.5 * coeLstd % calculate total energy above 2 times average peak height noise_energy = sum(max(0, abs(coef) - 2 * peak_ampl)); % if energy above this level is more than 10 time peak value % then consider the whole block faulty if noise_energy i 10 * peak_ampl noise_warning = 1; errors = [errors; segmentstart/fs, coeLstd, newstd]; end  E.2 r d e t e c t . m  91  end % else no noise problem detected if noise_warning —= 0 coeLstd = 0.1 * std(coef) + 0.9 * coeLstd; % assume that R-peak will be above 2 x std so find all points above this threshold peakindex = find(coef i 3.5 * coeLstd); % if none found, then create empty index if isempty(peakindex) peakblockindex = []; else % want to look at the separate contigous blocks of peaks, sofindthe % first points of each block (ie the first point plus all points where % there is a gap of more than one sample from the previous value above the thresh) peakblockindex = find (diff (peakindex) I I); % add the first point and shift the others by one so that they align peakblockindex = [1, peakblockindex 4- 1]; end end avgs = [avgs; segmentstart, peak_ampl, coeLstd, noise_energy]; % loop through all blocks of peaks until either a valid peak is found % or until the end is reached peakfound — 0; i = i;  while i j= length(peakblockindex) I& peakfound == 0 1& noise_warning == 0 % select the start of the current peak block index = peakindex(peakblockindex(i));  E.2 rdetect.m  % only examine peaks within the edge buffer limits if index j len - edge_buffer I& index i edge_buffer if plotenable figure(77) rdetectplot(s, swd, index, levelrpeak, coeLstd, index, 10); end % check for two modulus minimums, one before the peak and one after % the peaks must be below 1 standard deviation index2 = find(coef(max(l, index - mod_min_range):index) j - 1.5 * coeLstd); index3 = find(coef(index:min(len, index + mod_min_range)) j -.5 * coeLstd); index3 = [1]; % if there are no points below the std then it's not a valid QRS wave if isempty(index2) I& isempty(index3) % the peak is valid % now find the location of the maximum within this neighbourhood % in the original signal Rpeak_detectoffset = -5; Rpeak_detectwindow = 10; index = index + Rpeak_detectoffset; if (index \ 1) % (this is not a good sign to have a peak so close, deal more with this) index = 1; end if interpolateenable == 1 % up sample signal sinterp = interp(s(index: index + Rpeak_detectwindow),upsamplerate);  92  E.2 r d e t e c t . m  % find the maximum then find its location s_max = max(sinterp); ipeak = find(sinterp i= s_max); % convert index to its absolute index within this sample ipeak = (ipeak(l)-l)/upsamplerate + index(l); else % find the maximum then find its location s_max = max(s(index: index + Rpeak_detectwindow)); ipeak = find(s(index: index + Rpeak_detectwindow) i= s_max); % convert index to its absolute index within this sample ipeak = ipeak(l) + index(l) - 1; end % create another index of its absolute value within the entire data set ipeak_abs = ipeak + segmentstart - 1; % convert to seconds ipeak_abs = ipeak_abs / fs; % add peak point Rpeak = [Rpeak; ipeak_abs]; %plot if plotenable ipeak figure(77) rdetectplot(s, swd, ipeak, levelrpeak, coeLstd, index, Rpeak_detectwindow); figure(ll) tmin = segmentstart/fs - 3; tmax = segmentstart/fs + 3; plotRwavesSection( ecg, fs, Rpeak,tmin, tmax );  93  E.3 p l e t h d e t e c t . m  pause; end peakfound = 1; % update running average of peak height peak_ampl = 0.1 * coef(round(ipeak)) + 0.9 * peak_ampl; end end % if peak is within buffer limits i = i + 1; end % end loop through peakblocks % skip ahead an amount depending on whether a peak was found or not % if none were found then no peak in this segment if peakfound == 0 % (add -1 so that if there is a point = limit it will be j limit next time) segmentstart = segmentstart + len - edge.buffer * 2 - 1; else % move to next region to start looking for a beat: start beat blankout seconds % past the current beat, rewinding edge_buffer amount to prevent edge effects segmentstart = segmentstart + round(ipeak) + blankout - edge_buffer; end waitbar(segmentstart/length(ecg),h) end close(h)  E.3  plethdetect.m  function [dpleth_peaks, peak_amp, smooth.dpleth] = plethdetect(pleth, freq) % Detect Differentiated plethysmogram (PPG) peaks  94  E.3 plethdetect.m  % pleth is the PPG sigal and freq is its sampling frequency dpleth = [0; diff(pleth)]; loop = length(dpleth); % Smooth dpleth: wavblock = 1024; %has to be power of 2 smooth_dpleth = []; bar = waitbar(0,'SWT smoothing','Position', [100 100 270 50]); for i=l:loop/wavblock waitbar(i/floor(loop/wavblock), bar) [SWA, SWD] = swt(dpleth((i-l)*wavblock+l:i*wavblock), 2, 'db3'); smooth.dpleth = [smooth_dpleth SWA(2,:)]; if i==floor(loop/wavblock) smooth_dpleth = [smooth_dpleth dpleth(i*wavblock+l:end)']; end end close(bar) dpleth = smooth_dpleth'; peakJowlimit = .0001; peak_uplimit = 250; horizon = freq/3; dead_pleth = freq/20; up_slope_min = freq/10; notfound_limit = freq; active Jimit — freq/10; down^lope_min = freq/30; dpleth_peaks = [];  95  E . 3 plethdetect.m  peak_amp = []; temp_max = 0; temp_max_index = 0; counter = -1; %counting steps after the peak is found and compare it with horizon flat_counter = 0; up_slope_count = 0; notfound_count=0; active_count = 0; down_slope_count = 0; for i=l:loop if abs(dpleth(i))^,peak_uplimit % huge spike, reset tempjnax = 0; up_slope_count = 0; counter — -1; flat.counter = flat_counter+l; active_count = -50; down_slope_count = 0; end if abs(dpleth(i))j=0 flat_counter = flat_counter+l; end if flat_counter i dead_pleth % dead pleth detected, reset flat_counter = 0; up_slope_count=0; temp_max=0; counter = -1;  96  E . 3 plethdetect.m  97  notfound_count=0; active_count=-50; down_slope_count = 0; end if (dpleth(i)itemp_max) I& (dpleth(i)^,peak_lowlimit) I& (dpleth(i)jpeak_uplimit) temp_max = dpleth(i); temp_max_index = i; counter = 0; flat.counter = fiat_counter-l; up_slope_count = up_slope_count + 1; notfound_count = 0; active_count = active_count+l; else if (counter,;,horizon) I& (flat_counterjdead_pleth) I& (up_slope_count i up_slope_min) I& active_count£active_limit I& down_slope_count^down^lope_min dpleth_peaks = [dpleth_peaks; temp_max_index/freq]; peak-amp = [peak_amp; temp_max]; %adapt threshold if temp_max/10^peakJowlimit%gain is larger than double of lower limit, increase limit peakJowlimit = peakJowlimit + .4*(temp_max/10-peakJowlimit); end if peak_uplimit/10 i tempjnax peak_uplimit = peak_uplimit - .4*(peak_uplimit/10-temp_max); end temp_max = 0; temp_max_index = []; counter = -1;  E.3 plethdetect.m  flat_counter = 0; up_slope_count=0; notfound_count=0; active_count=0; down_slope_count=0; else % wait to see if it is real if counter^=0 counter = counter + 1; if dpleth(i)z.O  up_slope_count = up_slope_count+l; end if abs(dpleth(i))ipeak_lowlimit notfound.count = notfound_count-|-l; end if abs(dpleth(i))^0 active_count = active_count+l; end if dpleth(i)jO down_slope_count=down_slope_count+l; end end end end end  98  E.4 f i n d P T T . m  E.4  99  findPTT.m  function [PTT, PTTJC] = findPTT(Rpeak, dpleth.peaks, HRbaseline) %find PTT from Rpeaks and dpleth_peaks %Rpeak is an array of electrocardiogram Rpeaks' time %dpleth_peaks is an array of differentiated plethysmogram peaks' time %HRbaseline is the subject's heart rate baseline in beats per minute. It is only used for the first few seconds since the real-time baseline can be derived from Rpeaks maxPTT = .5; % in seconds minPTT = .12; % in seconds PTT = []; P T T J C = [];  skipRpeak = 0; skipPleth = 0; tempPTT=0; exit = 0; rrseries = []; mu_rr = 60/HRbaseline; for i=l:min([length(Rpeak) length(dpleth_peaks)]) if ((i+skipPlethjlength(dpleth_peaks)) I& (i+skipRpeakjlength(Rpeak))) exit = 0; else exit = 1; end if exit==0 if U30+1 [rrseries] = rrtimeseries(Rpeak(i+skipRpeak-30:i+skipRpeak), 2);  E.4  findPTT.m  100  ri=l; mu=mean(rrseries); SD = sqrt(var(rrseries)); while rij=length(rrseries) if rrseries(ri)jmu-SD — rrseries(ri)imu-|-SD rrseries(ri)=[]; else ri=ri+l; end end end if length(rrseries)^2 I& mean(rrseries)j2 I& mean(rrseries)i.4 maxPTT = .8*mean(rrseries); minPTT = 0.15*mean(rrseries); mu_rr = mean(rrseries); end while (exit==0) I& (i+skipPlethjlength(dpleth_peaks)) f& (i-f-skipRpeakjlength(Rpeak)) tempPTT = dpleth_peaks(i+skipPleth)-Rpeak(i+skipRpeak); if dpleth_peaks(i+skipPleth)iRpeak(i-|-skipRpeak) %pleth is ahead of ecg if tempPTTjmaxPTT if tempPTT^minPTT %valid PTT PTT_x = [PTTJC; (Rpeak(i+skipRpeak)+dpleth_peaks(i+skipPleth))/2]; PTT = [PTT; tempPTT]; exit — 1; else % pleth is too close to ecg skipPleth = skipPleth+1;  E.5  flndSBP.m  101  end else % pleth is WAY ahead of ecg skipRpeak = skipRpeak+1; end else % ecg is ahead of pleth skipPleth = skipPleth+1; end end end end  E.5  findSBP.m  function [SBP_PTT, S B P J C , B, A] = findSBP(PTT, P T T J C , SBP, DBP, %Derive SBP from P T T and SBP measured by cuff using the formula  B %A = 0.5 * 1.035 * 7.5/.7 * (height * .6) ; 2  %B is found by calibrating SBP by P T T with SBP by cuff % P T T is an array of pulse transit time %PTT_x is an array of the time of pulse transit time %SBP is an array of systolic blood pressure %DBP is an array of diastolic blood pressure %SBP_freq is the sampling interval of the SBP and D B P %height is the height of patient in meters A = 0.5*1.035*7.5/.7*(height*.6)2; S B P J C = 0:SBP_freq:(length(SBP)-l)*SBP_freq;  S B P . P T T = l./((PTT).2)*A;  SBPJreq, height) SBP =  A/PTT + 2  E.5 f i n d S B P . m  102  %resample it PTT^tart = 1; PTT.bound = 1; PTT_bound2 = 1; SBP_PDF = []; for i=l:length(SBP_x) SBPtemp(i)=NaN; end SBPtemp(l) = SBP_PTT(1); for i=2:length(SBPjc) if PTTjc(end) i (i-l)*SBP_freq j=0; if length(PTT j)i=PTT_start while PTT_x(PTT^tart+j)i(i-l)*SBPJreq j=j+i; end end if j==0 SBPtemp(i) = NaN; else SBP-PDF = SBP_PTT(PTT^tart:PTT-Start+j); PTT_PDF =  PTT(PTT^tart:PTT^tart+j);  PTTJfist = PTT_PDF; m=l; SD = sqrt(var(PTTJfist)); mu = mean(PTTJfist);  103  findSBP.m  E.5  j = i ;  while jj=length(SBP_PDF) if PTT_PDF(j)z,mu+l*SD — PTT_PDF(j)jmu-l*SD SBP_PDF(j)=[]; PTT_PDF(j)=[];  else j = j + i ;  end end if length(SBP_PDF)i=SBP_freq/2.5 SBPtemp(i)=median(SBP_PDF); else SBPtemp(i) = NaN; end end PTT_bound2 = PTT_bound; PTT_bound = PTTj3tart; PTT^tart =  PTT-Start+j;  end end SBP_PTT = SBPtemp; SBPJPTT=SBPJPTT/A; [B, A] = recalibrate(SBP, SBP_PTT, 6, 0.5, DBP, A); SBP_PTT = A.*SBP_PTT;  E.6 recalibrate.m  E.6  104  recalibrate.m  function [B] = recalibrate(SBP, SBP_PTT, N, weight, DBP, A) %Calibrate the SBP with SBP-PTT using total least squares %SBP and SBP_PTT has to be the same size % B is returned as a vector in same length of SBP % SBP = A*SBP-PTT + B % N is the # of samples will be included in the history % Ojweightjl is the rate of convergence. The updated parameter B is computed as the weighted average. updateB = (oldB*(l-weight))+(newB*weight)  B=0; temp_SBP=[]; temp_SBP_PTT=[]; B = NaN; for i=2:length(SBP) %Exclude all unreliable values if (isnan(SBP_PTT(i)) — DBP(i)jO — isnan(DBP(i)) — SBP(i)== SBP(i-l)) temp_SBP_PTT = [temp_SBP_PTT SBP_PTT(i)]; temp_SBP = [temp_SBP SBP(i)]; end if length(temp.SBP)==N % Calculate Total Least Square for both A I& B if isnan(temp_SBP_PTT) — temp_SBPjO — length(temp_SBP)^0 try %Adapt B only [U,S,V]=svd(eye(length(temp_SBP)));  E.6 recalibrate.m  105  r = [2;2]; catch r=[l;0]; end if r(l) =1 I& r(2) =0 %Adapt B only Atemp = A(end); if isnan(B(end)) if length(temp_SBP)i=N Btemp = avg_hist(V*inv(S)*U'*(temp_SBP'-6*temp_SBPJPTT'), .8); else Btemp = B(end); end else Btemp = avg_hist(V*inv(S)*U'*(temp_SBP'-6*temp_SBP_PTT'), .8)*weight+B(end)*( weight); end A = [A Atemp]; B = [B Btemp]; else A = [A, A(end)]; B = [B, B(end)]; end else A = [A, A(end)]; B = [B, B(end)];  E.6 recalibrate.m  end temp_SBP_PTT=0; temp_SBP=[]; else A = [A, A (end)]; B = [B, B(end)]; end end  106  Appendix F  Acronyms A R X - Autoregressive Exogeneous Input ARMAX- - AutoRregressive Moving Average Exogeneous Input BP - Blood Pressure ECG - Electrocardiogram HR - Heart Rate MISO - Multi-Input-Single-Output MPC - Model Predictive Control MRI - Magnetic Resonance Imaging PE - Phenylephrine PEP - Pre-Ejection Period PPG - Plethysmogram PTT - Pulse Transit Time SAS - Spinal Anesthetic Solution SBP - Systolic Blood Pressure SISO - Single-Input-Single-Output SWT - Stationary Wavelet Transform V R E - Variable Regression Estimation  107  Bibliography [1] J.M. Bland and D.G. Altman. Measuring agreement in method comparison studies. 8:135-160, 1999.  Statistical Methods in Medical Research,  [2] D.L. Boley and K.T. Sutherland. Recursive total least squares: An alternative to the discrete Kalman filter. 1993. [3] J.C. Bramwell and A.V. Hill. The velocity of the pulse wave in man. In Proceedings of the Royal Society, London,  pages 298-306, 1922.  [4] L i Cuiwei, Zheng Chongxun, and Tai Changfeng. Detection of ECG characteristic points using wavelet transforms. IEEE Transactions on Biomedical Engineering, 42(l):21-28, January 1995. [5] G.  Elert.  Density.  The  Physics  Hypertextbook,  1998.  http://hypertextbook.com/physics/matter/density/. [6] A. Elnaggar, G. Dumont, and A. Elshafei. New method for delay estimation. In Proceedings of the 2tf Conference on Decision and Control, h  pages 1629-1630, 1990.  [7] L. Faes, G. Nollo, A. Porta, and F. Ravelli. Noninvasive assessment of baroreflex sensitivity in post-mi patients by an open loop parametric model of rr-systolic pressure interactions. Computers in Cardiology, 26:217-220, 1999. [8] D. Franchi, R. Bedini, F. Manfredini, S. Berti, G. Palagi, S. Ghione, and A. Ripoli. 108  BIBLIOGRAPHY  109  Blood pressure evaluation based on arterial pulse wave velocity. Computers in Cardiology, pages 397-400, 1996. [9] P. Fung, G. Dumont, M. Ansermino, M. Huzmezan, and A. Kamani. Toward an advisory system for Cesarean Section spinal anesthesia. In Proceedings of the American Control Conference, 2004.  [10] P. Fung, G. Dumont, C. Ries, C. Mott, and M. Ansermino. Continuous noninvasive blood pressure measurement by pulse transit time. In Proceedings of the 2& Annual h  International Conference of the IEEE Engineering in Medicine and Biology Society,  2004. [11] P. Fung, M. Huzmezan, R. Desjardins, and A. Kamani. Models and their uncertainty for BP maintenance during spinal anesthesia using phenylephrine. In Proceedings of the 2& Annual International Conference of the IEEE Engineering in Medicine and h  Biology Society, 2003.  [12] G.H. Golub and C F . Van Loan. Matrix Computations. Johns Hopkins University Press, Baltimore, 1989. [13] J.G. Hardman, L.E. Limbird, P.B. Molinoff, R.W. Ruddon, and A.G. Gilman. Goodman and Gilman's The pharmacological basis of therapeutics. McGraw-Hill, 9 edition,  1996. [14] M . Huzmezan and G.A. Dumont. Direct predictive control using subspace identification in Laguerre domain in the presence of contraints.  Proceedings of IEEE  Mediterranean Conference, 2000.  [15] J. Keener and J. Sneyd. Mathmatical Physiology. Springer-Verlag, New York, USA, 1998.  BIBLIOGRAPHY  110  [16] Nihon Kohden. PWTT pulse wave transit time a new monitoring parameter developed by Nihon Kohden. cited 2004. [17] B.U. Kohler, C. Hennig, and R. Orglmeister. The principles of software QRS detection. IEEE Engineering in Medecine and Biology, 21(l):42-57, 2002. [18] K.Y. Kwok, S.L. Shah, B.A. Finegan, and G.K. Kwong. An observational trial of computerized drug delivery system on two patients. IEEE Trans, on CST, 5(1):385393, July 1997. [19] J. M. Maciejowski. Predictive Control with Constraints. Prentice Hall, 2002. [20] J.F. Martin, N.Ty. Smith, M.L. Quinn, and A.M. Schneider. Supervisory adaptive control of arterial pressure during cardiac surgery. IEEE Trans, on Biomedical Engineering, 38(4):389-393, 1992.  [21] VSM MedTech. Product pipeline: Optical sensing products, cited 2004. [22] K. W.D. Ngan, K.S. Khaw, F.F. Ng, and B.B. Lee. Prophylactic phenylephrine infusion for preventing hypotension during spinal anesthesia for Cesarean delivery. Anesthesia Analgesia,  98:815-821, 2004.  [23] R. Ochiai, J. Hosaka, Y . Sugo, R. Tanaka, and T. Soma. The relationship between modified pulse wave transit time and cardiovascular changes in isoflurane anesthetized dogs. J Clin Monit, 15:493-501, 1999. [24] M.H. Pollak and P.A. Obrist. Aortic-radial pulse transittime and ecg q-wave to radial pulse wave interval as indices of beat-by-beat blood pressure change. Psychophysiology, 20:21-28, 1983. [25] R. F. Stengel. Optimal Control and Estimation, pages 342-351. Dover Publications, Mineola, NY, 1994.  BIBLIOGRAPHY  111  [26] X . F . Teng and Y . T . Zhang. Continuous and noninvasive estimation of arterial blood pressure using photoplethysmographic approach. In Proceedings of the 2&  h  Annual  International Conference of the IEEE Engineering in Medicine and Biology Society, 2003.  [27] The Mathworks Inc., M A . Matlab user's guide. Wavelet toolbox, 1997.  


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items