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

Download

Media
831-ubc_2004-0452.pdf [ 4.39MB ]
Metadata
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
831-1.0065520-fulltext.txt
Citation
831-1.0065520.ris

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 FACULTY OF 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 FACULTY 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^rV funtx 34 / f l G / loo H Name of Author (please print)J Date (dd/mm/yyyy) Title of Thesis: Mv^ory Sy <?fen* -por A d^^'yf ration o-f fhe^i fr.pkr!ca e D e 9 r e e : M r V S c _ _ _ _ _ Year: 2qq <j-D e p a r t m e n t o f E l e c t f l c ^ l P ^ Q M T W h^l^r.'A. The University of British Columbia ' J J Vancouver, BC Canada grad.ubc.ca/forms/?formlD=THS page 1 of 1 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 Sec-tion. Hypotension can cause serious fetus hypoxia, therefore maternal blood pressure must be kept above a minimum level. Phenylephrine dosage is mainly determined in a heuris-tic 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 qual-ity 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 phenyle-phrine'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 T a b l e o f C o n t e n t s Table of Contents iii 1 Introduction and Background 1 1.1 Background 1 1.2 Objectives ' 3 1.3 Thesis Overview and Contributions 5 2 Noninvasive BP Measurement by Pulse Transit Time 7 2.1 Literature Review on Pulse Transit Time 8 2.2 Mathematical Model 9 2.3 Pulse Transit Time Computation 15 2.3.1 R peak Detection 15 2.3.2 PPG Maximum Slope Detection 17 2.3.3 PTT Computation Algorithm 18 2.4 Limitations of the PTT-BP Model 18 2.5 PTT-BP Model Validation 19 2.6 Conclusions 25 3 Patient Modeling N 27 3.1 MISO System Identification 28 3.1.1 MISO Method 28 3.1.2 MISO Model and Model Uncertainty 30 3.1.3 Results of the MISO System Identification 34 3.2 Single Variable Approach Using Least Squares Estimation on the SAS Re-sponse 35 3.2.1 SISO System Identification Method 38 3.2.1.1 SISO Model and Model Uncertainty 40 3.2.1.2 Results of the SISO System Identification 43 3.2.2 Positive SISO Models 44 3.2.2.1 Positive SISO Model and Model Uncertainty 47 3.2.2.2 Results of the Positive SISO System Identification 48 3.3 Conclusions 51 iii CONTENTS iv 4 Advisory System Design 53 4.1 Advisory System User Interface 55 4.2 Controller Design and Tuning 56 4.3 Safety Concerns 60 4.4 PE Gain Adaptation and SAS Time Delay Estimation 60 4.5 Advisory System Simulation 65 4.6 Conclusions 70 5 Conclusions 71 5.1 Main Contributions 71 5.2 Limitations of the Advisory System 72 5.3 Direction for Future Research 73 A Data Acquisition 74 A . l Procedures 74 A. 2 Conclusions 75 B Advisory System Structure Overview 77 B. l SBP Regulator Interface 77 B.2 SBP Regulator 80 C Variable Regression Estimation 82 D Matlab Code for Positive System Identification 83 D.l processdata.m 83 D. 2 removespinal.m 85 E Matlab Code for Noninvasive BP Measurement by PTT 88 E. l PTT2BP.m 88 E.2 rdetect.m 88 E.3 plethdetect.m 94 E.4 findPTT.m 99 E.5 findSBP.m 101 E.6 recalibrate.m 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 Sphygmomanom-etry during Cesarean Section (Case 10) 25 3.1 MISO Patient Model 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 Signif-icant 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 Fig-ure 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 PE Model with Tuning Q=1000, R=l , Hu=l, Hp=6 59 4.3 Combining BP measured by cuff with BP 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 a2 = lOmmHg Advisory System Simulation. . . . 66 4.9 Doubled Gain Model Advisory System Simulation 66 4.10 Oscillating Model in Section 3.2.1.1 Advisory System Simulation 67 4.11 Perfect Model Closed Loop Simulation 67 4.12 Perfect Model with Noise a2 = lOrnmHg Closed Loop Simulation 68 4.13 Doubled Gain Model Closed Loop Simulation 68 4.14 Oscillating Model in Section 3.2.1.1 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 List of Tables 4.1 Controls on the Advisory Page of the System 55 4.2 Indicators on the Advisory Page of the System . 56 4.3 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, computer-aided 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 en-gineering 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 Sec-tion 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 hy-potension [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 un-der 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 hu-man 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 advi-sory 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, con-straints, 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 demon-strates 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 Anesthesia -1 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 PE 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 Overv iew and Contr ibut ions 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 sphyg-momanometry 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 Ce-sarean 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. More-over, 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 measure-ment. Therefore, a noninvasive method for measuring BP is desirable for Cesarean Section patients without significant organ dysfunction during short duration surgery but with ex-pected 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 conjunc-tion 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 (2.1.1) P blood density a Aa/AP arterial compliance 2.2 Mathematical Model 9 As the compliance Ca is denned as A a / A P , the above equation can be used to relate pulse wave speed to BP. Although the model is capable of inferring the BP if the parameters are known (a, p, Ca and Aa) , those parameters are not available unless we use technology 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 PTT 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: 2.2 Mathematical Model 10 F • d = ^mv2 + mgh where F = force exerted on blood d = distance from heart to fingertip m = mass of blood (2.2.1) v = pulse wave velocity g = 9Mm/s2 h = height difference between two sites 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 (2.2.3) la - a a • a The term = 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 l ^ + pgh (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/m3 [5]. 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) (2-2.5) A i p — PTT'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)2 x £ (2.2.6) From the above calculations, BP can be estimated from PTT and a few average em-pirical 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 = pj?r-> + B, as shown in Figure 2.2. 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 _ A , p f(L) = E i l i l K I I (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: B = BP PTT BP BPl BP2 BPN . PTTX PTT2 PTT2 1...0 : -.0 0...1 where (2.2.8) PTTN 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 0...1 = UY,VT (2.2.9) B = VY,-WT BP-PTT2! 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 A P = ^ A a (2.2.10) a 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 = ^Aa p y y 3 + Diastolic BP Then one may be able to estimate the systolic BP by measuring arterial cross section area. However, noninvasive measurement of a and Aa 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 15 2.3 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 PPG 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 differen-tiated 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 compo-nents. 2.3.1 R peak D e t e c t i o n 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 imple-mentation 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 PTT from E C G and P P G 100 1550.2 1550.4 - - P P G O P P G max Slope E C G * E C G R peak 1550.8 1551 Seconds Figure 2.3: Definition of PTT. performed. The highest frequency decomposition captures much of the noise, and the low-est 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 17 2.3.2 P P G M a x i m u m S lope D e t e c t i o n 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; If dPPG(i) > temp_max & dPPG(i) < MAX & dPPG(i) > MIN Yes Store temporary maximum temp_max = dPPG(i); temp_max_index = i; up_slope_counter++; flat_counter - -; look_forward_counter=0; _NoJ If up_slope_counter > UPSLOPEMIN & down_slope_counter > DOWNSLOPEMIN & flat_counter < FLATMAX & look forward counter > HORIZON Yes Nc4 Store temp_max_index)-*| as Peak If temp_max==null No I look forward counter++ TJ Reset Set up_slope_counter =0; down_slope_counter =0; flat_counter =0; look_forward_counter =0; temp_max = null; temp_max_index= null; Yes _Yes. Yes lfdPPG(i)==0 1 No~ flat_counter++; W If flat_counter >= FLATMAX No lfdPPG(i)>0 No I . . I Yes down_slope_counter++; up_slope_counter++; Figure 2.4: Flow Chart of Maximum Slope Detection in PPG. 2.4 Limitations of the P T T - B P Model x 18 2.3.3 P T T Computat ion Algor i thm In our current set of data, ECG and PPG were collected at 300Hz though the PPG 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 = 30 • RRI, DOWNSLOPEMIN = 5 • RRI, FLATMAX = 15 • RRI and HORIZON = 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 x PTT Right finger x x xx xx * x * x x xSxv&x*," x v X X Ear < x * x><x x x Left finger x x v v X>X X y < „ X . . X Toe 100 200 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 significant fluctuations in 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 140 i > i i i — SBP by PPT SBP by Cuff — Invasive SBP 120 - -100 - -O) I E E 80 I i r i ! i V s i 1 i i iff ; " ' V f | ^ A ^ . ! 1 60 i i j i i ( I i i i i i i 1 i i i 40 -i i i i i -I I I I I 1 — 1000 2000 3000 4000 5000 6000 Seconds 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-2 0 0 0 -1 5 0 0 h 1 0 0 0 r -5 0 0 r-- 6 0 - 4 0 - 2 0 0 2 0 4 0 6 0 8 0 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 PTT 60 80 100 120 140 160 180 Avg of S B P by cuff and S B P by PTT (mmHg) 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 140 120h ixl Q_ O O) — 100 o a) £ at a) § 80 ca c 'CL <n "o * 60 O) X E E 40 20 500 1000 Seconds 1500 — SBP by PTT - - SBP by Cuff O Spinal Anesthetic Phenylephrine 2000 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 den t i f i ca t i on 28 In Section 3.1, the multivariable approach is discussed and essential results are pre-sented. 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 M I S O M e t h o d Initially, as shown in [7], it was believed that BP is affected by the injection of PE 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 MISO 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. C D -PE z2+b21.z+b31 z 3-alz 2-a2z-a3 PE model SAS Dead time Discrete Integrator Z2+b22.Z+b32 z 3-alz 2-a2z-a3 L SAS model [HR Baseline] [SBP Baseline] 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 sys-tem identification because controller design is performed in state space. Since the sub-space 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 ARMAX 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 31 ' al 1 0 " " 1.2456 1 0 " A = a2 0 1 = 0.0919 0 1 (3.1.1) a3 0 0 . -0.3523 0 0 . B Ki K2 K3 K\b2i K2b22 K3b23 Kihi K2b32 K3b33 0.3893 -0.6213 0.0091 0.2842 -1.0247 0.0062 -0.0436 0.0901 -0.0016 (3.1.2) C=[ 1 o 0 ] (3.1.3) Ks = 0.0153 (3.1.4) PE SAS HR u BjKK> dB dAH C SBP 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 ARMAX identification and it represents an average of results from the subspace 3.1 M I S O S y s t e m I den t i f i ca t i on 32 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. Distribution of a1 2r-1.5 -1 -0 .5 --2.5 -2 -1.5 -1 -0.5 0 Distribution of a2 -1 -0.5 0 0.5 1 1.5 2 Distribution of a3 -0.4 -0.2 0 0.2 0.4 0.6 0.8 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 vari-ances and they are related to the linear fractional uncertainty model as the 95% confidence interval presented in Figure 3.2. da\ 0 o • " 1.288 0 0 " dA = ± da2 0 0 = ± 2.173 0 0 . da3 0 0 . . 0.935 0 0 . 3.1 MISO System Identification 33 Figure 3.4: Distribution of First Column of B. Distribution of Spinal Anesthesia's K2 -10 - 5 0 5 10 15 Distribution of Spinal Anesthesia's K2b22 -10 -8 - 6 -4 - 2 0 2 4 6 Distribution of Spinal Anesthesia's K2b32 -10 - 5 0 5 10 15 20 Figure 3.5: Distribution of Second Column of B. 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 d613 db2i dbii dbn dbn dbs2 dhs 6.485 4.129 5.923 7.299 5.510 8.270 0.554 0.766 0.555 3.1.3 R e s u l t s of the M I S O S y s t e m Ident i f i ca t ion (3.1.6) 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 {K2). 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 W 100 80 model actual " -\ r t /i 2500 1000 1500 Time (sec) Figure 3.7: Case 14: Individual Model Validation. 3.2 Single Variable Approach Using Least Squares Estima-tion 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 K2 during Model Validation. 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 S I S O S y s t e m Ident i f ica t ion 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 re-sponse 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 least-squares 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 non-hypotensive patient seems to have a BP drop greater than 20% of the her baseline. How-ever, 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 80 a. o o i 60 oo » 4 0 o * 20 S A S + P E 0 200 400 600 800 1000 1200 1400 1600 S e c o n d s 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 SAS Response 0 -0.05 -CD c >^ -0 1-co u-' CO Xi a. CQ CO C 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 devel-oped (see Section 3.2.2). 3.2.1.1 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 1 1 1 i 1 1 1 1 1 600 800 1000 1200 1400 1600 1800 2000 2200 2400 Seconds Figure 3.14: Subtraction of the SAS Response From the raw data (same case as in Fig-ure A. l ) . Pole-Zero Map 11 1 1 1 r - r ^ 1 • - •. . i . . 1 1 1 1 0.8- : 0.6-0.4 / . --04 - '•. y --0.6- : -0.8- '• . .,1 1 1 1—"••••I — I — I- ••• "—I 1 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Real Axis Figure 3.15: Distribution of the Two Poles of PE. 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™ ' Z (3.2.1) (z-p)(z-p) KPE = 2.795e~4, 0 < KPE < 1.6917e-3 (3.2.2) p = 0.9726 -MO. 1478, 0.9414 < Re(p) < 0.9923, (3.2.3) 0.0381 < Im(p) < 0.2083 The SAS ARX model is: SAS model(z) = ., K s A S ' z _ (3.2.4) v ; zk(z- 1)(Z-PSAS) Therefore, each patient has an individual gain and time delay, PSAS = 0.863 (3.2.5) (3.2.6) (3.2.7) KSAS = -2.4195e-4, -7.2951e"4 < KSAS < -2.7158e-5 k = 20.56, 0 < k < 56 3.2.1.2 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 PE 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 Model Output SBP data from PTT + PE 0 SAS £ 80 ^ + 4 - 4 - + 4 o oi zs I 6 0 -E <D 40 - + 20 - + 1 1 I 1 1 I l I 800 1000 1200 1400 1600 1800 2000 2200 Seconds Figure 3.18: Case 10 PE Model Validation. 3.2.2 P o s i t i v e S I S O M o d e l s While the model presented in Section 3.2.1.1 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 ARX 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 1500 1600 1700 1800 Seconds Figure 3.19: Case 12 PE Model Validation. Model Validation Figure 3.20: Case 14 PE Model Validation. 3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 46 Model Validation "n r~ co •s £ 60 e 40 Model Output — SBP data from PTT + PE O SAS 1000 1200 1400 1600 1800 2000 2200 2400 2600 Seconds Figure 3.21: Case 19 PE Model Validation. Model Validation n r~ 120 E 60 E Model Output SBP data from PTT + PE O SAS 800 900 1000 1100 1200 1300 1400 1 500 1600 1700 1800 Seconds Figure 3.22: Case 12 PE Model Validation with Customized Gain. 3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 47 positive system is also expected from compartmented models in pharmacokinetics. PE model(z) = , K p S z , = £PEZ (3.2.8) 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 a2 that would 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. Pole-Zero Map 1 1 i i 1 i i i~~ 4- First pole x Second Pole 0.6 -0.4 -0.2 -0 - x x x x x x x a c x + x f x x * x-tf xf++ -j* -t*f + +-tHH--0.2 --0.4 --0.6 -i i i i i i 1 1— 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 Figure 3.23: Distribution of the Two Real Poles of Positive PE Model. 3.2.2.1 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 3.2.1.1. 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. KPE = 0.0022, (3.2.9) 0 < KPE < 0.1334 pi = 0.8149, pi = 0.7156, (3.2.10) 0.6940 <pi< 0.8782, 0.5764 <p2< 0.8461 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. There-fore, the average model is kept as the nominal PE patient model. 3.2.2.2 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 Gain 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 - - SBP data from PTT + PE o SAS 801 20h 800 1000 1200 1400 1600 1800 2000 2200 Seconds Figure 3.25: Case 10 PE Model Validation. 3.2 Single Variable Approach Using Least Squares Estimation on the SAS Response 50 Model Validation 120 E 60 E 20 Model Output SBP data from PTT + PE O SAS 800 900 1000 1100 1200 1300 1400 1500 1 600 1700 1800 Seconds Figure 3.26: Case 12 PE Model Validation. Model Validation 120 E 60 E 20 — Model Output - - SBP data from PTT + PE O SAS + + 1000 1200 1400 1600 1800 2000 2200 2400 2600 Seconds Figure 3.27: Case 14 PE Model Validation. 3.3 Conclusions 51 140 ! 100 I 60 £ Model Validation Model Output SBP data from PTT + PE o SAS 1000 1200 1400 1600 1800 2000 2200 2400 2600 Seconds 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 52 for real-time advisory control application of the patient model. 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 user-friendly 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 P l a n t S e c o n d s | C o n t r o l s - S e c t n o s j A d a c t a a m A d v i s o r y j s u n u a c w i j Information aeset rh(? ^yitem after modifying 'Ir^ronftohon' F* m the patients bese'ine Log Path C B D C ^ f a l i n a «e: \parryWC*_=ro1ett \ J«j*J D D r D d i c l l M " i sBP.-eaJatof'/eco-i.txl Controls SBP Set Point Total Phenyteplrre Dose _ 1 AdmnistBfed(LQ);|So3r Toggle svrftch Last Dose:|ioo,oo Automatic -i Advi iory Dose Suggest Phenylephrine Dose(ug): : 0 5 Over r ide Phenylephr ine DoseCug): g, 0 ( 1 | S6Pf.Ptwakdars wstilrifeOffln CumtttTrw.. ObsawSBP Prnfinir-i. with K* Infaion 80.0 •:;3:10PW 05:!8:00PH OS-.2JJ0OPM 0504:00PM 05:27:00PM C5:30:Q0W Time Phenylephrine Infuson Rate Phenylephrine Plot second 10 0-7 5.-Jf s.o-2.5-1 o.o-1 OSlUliCPM Phenylephrine Range o w a o o m c m i i o o p M m&mm vw*om <B:».mm Time ~~"jn)ercMi Ij Constraints M a x Bolus M a x Infusion Rate r l 100.00 ug l-i 0.83 ug/see Max Total Phenylephrine ii 1000.00 ug Adaptive IMIwI^ I Disturbance ^ A ^ s m e M . j o S ''• DeadTime. 16.00 learnt Mbce! Ihftiaf Hadet E n a b t e G a i n A d a p t a t i o n : # o f New G a i n A d a p t e d i Pebad Irfltai Model A d a p t A f t e r j I S O ( S e c o n d s ) Figure 4.1: Screen Shot of the User Interface of the Advisory System. 4.1 Advisory System User Interface 55 4.1 Advisory System User Interface There are five pages in the user interface: 'Plant Settings', 'Controller Settings', 'Adapta-tion', '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 on the Advisory Page of the System. Controls SBP Setpoint The user chooses the desirable BP setpoint. Max Bolus The user selects the input constraint on maximum bolus. Max Infusion Rate The user selects the input constraint on maximum infusion rate. Max Infusion Rate The user selects the input constraint on maximum total amount of PE. Toggle Switch The user changes the system between advisory and automatic mode. Go! Dose The user can administer the dose displayed in Suggest Phenylephrine Dose(ug) by pressing the button. Override Phenylephrine Dose(ug) The user can input the override dose. Override The user can administer the dose displayed in Override Phenylephrine Dose(ug) by pressing the button. Enable Gain Adaptation: 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. Reload Initial Model The user can load the original model and discard the learnt model by pressing the button. Adapt After The system observes the response after each individual PE 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 Dose(ug) The optimal PE dose suggested by the system. 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 Adapted The indicator displays the number of adaptations the system has 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) = Ax(k) + Bu(k) y(k) = Cx(k)+d(k) SBP(k) = y(k) • SBP baseline + SBP baseline sampling interval = 10 seconds 0 1 1 (4.2.2) -0.5832 1.5305 0 0.0022 1 0 A B C In this SISO system, u(k) is the PE infusion rate in jUg/second and y(k) is the (SBP-SBP 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 (4.2.3) y(k) = A 0 B CA I CB Au(k) 0 I S(k) + d(k) 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 10 - 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 BP 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 -* 2 0 r -10 --20 --30 -0 F = -45 --90 -3 "135 -~i ~m' a -225 --270 --315 --360 L 10~3 10"2 10"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 anesthe-siologist 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 ei-ther 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 t^g of PE, and the anesthesiologist's direct decision. When the system predicts hypotension, a continu-ous 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 PE 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 61 If PTT SBP available If previous PTT SBP available If both previous and current cuff SBP available Use PTT SBP If cuff SBP available Use cuff SBP SBP not available Stop Advising T •Use.-. • . ::;•:] • previous PTT SBP + (cuff SBP - previous cuff SBP) Figure 4.3: Combining BP 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 BP reading arrives from the monitor, the computa-tion has to converge rapidly. The V R E algorithm calculates the cross correlation function between the SAS input and the BP 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 Adapta t ion and S A S T ime Delay Est imat ion 62 Table 4.3: Comparison of SAS Time Delay Estimation by V R E and Expert Knowledge. Case number Delay by Expert Knowledge (10 Sec) Delay by V R E Estimation (10 Sec) 10 4 4 12 18 27 14 18 9 17 25 26 • 19 15 8 20 36 36 23 10 6 26 38 50 36 16 14 37 4 5 45 24 26 46 29 36 47 56 4 49 15 4 50 37 12 52 18 7 53 16 17 54 48 4 56 18 9 60 10 14 63 10 14 64 0 5 It can be observed that in three particular cases, 47 and 50 and 54, the V R E estima-tion 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 dis-played 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 BP 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 PE Gain Adaptation and SAS Time Delay Estimation 6 3 given the lack of system excitation. 140 120 Cesarean Section Case 1500 Seconds SBP by PTT case47 O Spinal Anesthetic + Phenylephrine 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 PE administration is closely followed by another one. 4.4 PE 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. BP Adapt Data 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 65 4.5 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 per-formance 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 ~ ~ \ + + + + + + + + + + + + + + + + + + + + + 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 simu-lation. 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 120 UJ 80 40 20 4.5 Advisory System Simulation 66 Mode! Simulation Model Output — Observed Output + PE o Spinal Set Point Figure 4.8: Perfect Model with Noise a2 — 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 3.2.1.1 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 Model Output Observed Output + PE o Spinal — Set Point ^ ^ ^ ^ 20 ++ + + X -it- + ++ +-+ + + + + + + -* fH-H- # + + + • + + + + + ^ 1500 2000 Figure 4.12: Perfect Model with Noise a2 = lOmmHg Closed Loop Simulation. Model Simulation Figure 4.13: Doubled Gain Model Closed Loop Simulation. 4.5 Advisory System Simulation 69 Model Simulation 140 ! 100 ra 60 Model Output - - Observed Output + PE O Spinal — Set Point + 4-+ S + + A t + + + + -h + + + + 1000 1600 Seconds Figure 4.14: Oscillating Model in Section 3.2.1.1 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 3.2.1.1, 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 70 4.6 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 intro-duces 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 indepen-dently 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 phenyle-phrine. 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 cor-related 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 hy-potension 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 prop-erties. 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 ran-domly selected patients undergoing Cesarean Section spinal anesthesia at the British Columbia Women's Hospital. Patients with significant medical disease, pregnancy com-plications 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 admin-istration 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 7 5 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 per-formed 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 A.2 Conclusions 76 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 -60 -w < 40 -20 -500 SAS -l- PE 1000 1000 1500 Seconds 2500 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 RRI and differ-entiated PPG peaks calculated by Find Pleth Peaks. The resulting PTT is used by the 77 B . l S B P Regulator Interface 78 Datex-Ohmeda S5 Collect Time, Physiological Signals SBP Controlled Variable SBP Regulator Interface ECG, PPG E C G to RRI Wavelet | , SWTdec FGfi S-levfll SWT ECG Detect R Peaks RRI Sample Intervals R Peaks, RRI N Find PTT Find Pleth Peaks Differentiated p p G ^ Differentiated; PPG Peaks Flip Conv 3 SWT l_» Peakdetect [Approximation PTT Compute SBP by PTT SBP by cuff Differentiated PPG Peaks Recalibrate H Calibrated Offset SBP for Control SBP Controlled Variable SBP Regulator Phenylephrine Grasby 4300 Infusion Pump Phenylephrine Patient Infusion Rate Figure B . l : Advisory System Data Flow Diagram 1. B . l SBP Regulator Interface 79 Advisory Advisory Input Plant Dynamics Setpoint, Internal States Advisory Dose Prediction Plant Dynamics, Internal States, Advisory Dose Dose Prediction Display Plant Dynamics, Internal States, Feedback Gain Model Display Continuous .Prediction S B P Regulator C O R E : Get User Input and BP Index. Compute Phenyle-phrine Dose Prediction Adaptation / Initialization Plant Dynamics Kalman_design Observer Plant Dynamics 1 mpc_design Feedback Gain Plant Dynamics Plant Dynamics, Impulse Response! Peia^ Adapt_SBP Detrended S B P New Plant Dynamics File Path Simulation Phenylephrine Dose ReadDose From File Disturbance ReadSpinal FromFile Figure B.2: Advisory System Data Flow Diagram 2. B .2 SBP Regulator 80 SBP Regulator Interface to derive the noninvasive BP index discussed in Chapter 2. E C G to RRI - 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 - SBP 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) = $(t)9 = [-y{t-l),...,-y(t-n),u(t-k-l),...,u(t-k-m)} (C.0.1) 9 = [a1,...,an,bi,...,bm]T The cross correlation between input and the output increments is given by: Ei{tM) = iE5=o«(3-*i)(y(j' + i ) - » ( j ) ) ( C 0 2 ) 1 < ki < time of first PE dose — time of SAS A different cross-correlation E\ is computed from each guess of fc;. 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 D . l p r o c e s s d a t a . m 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 85 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 D .2 removespinal.m 86 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 87 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); Appendix E Matlab Code for Noninvasive BP Measurement by PTT 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, PTTJC, 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 elec-trocardiogram signal plotenable = 0; upsamplerate = 5; interpolateenable = 1; Rpeak = []; errors = []; avgs = []; 88 E.2 rdetect .m 89 % 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); 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 rdetect .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, so find the % 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 92 % 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); E.2 rdetect .m 93 % 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 ); E.3 plethdetect .m 9 4 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 E .3 plethdetect.m 95 % 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 = []; E.3 plethdetect .m 96 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; 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 98 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 E.4 f i n d P T T . m 99 E.4 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, SBPJreq, height) %Derive SBP from P T T and SBP measured by cuff using the formula SBP = A/PTT2+ 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 DBP %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; 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); E .5 findSBP.m 103 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 = P T T - S t a r t + 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 104 E.6 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 106 end temp_SBP_PTT=0; temp_SBP=[]; else A = [A, A (end)]; B = [B, B(end)]; end end Appendix F Acronyms ARX - 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 VRE - Variable Regression Estimation 107 Bibliography [1] J.M. Bland and D.G. Altman. Measuring agreement in method comparison studies. Statistical Methods in Medical Research, 8:135-160, 1999. [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] Li Cuiwei, Zheng Chongxun, and Tai Changfeng. Detection of ECG characteris-tic 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 2tfh Conference on Decision and Control, 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 B I B L I O G R A P H Y 109 Blood pressure evaluation based on arterial pulse wave velocity. Computers in Car-diology, pages 397-400, 1996. [9] P. Fung, G. Dumont, M. Ansermino, M. Huzmezan, and A. Kamani. Toward an ad-visory 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&h Annual 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&h Annual International Conference of the IEEE Engineering in Medicine and 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. Good-man 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 iden-tification 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. B I B L I O G R A P H Y 1 1 0 [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 detec-tion. 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):385-393, 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 En-gineering, 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. Psychophysiol-ogy, 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., MA. Matlab user's guide. Wavelet toolbox, 1997. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065520/manifest

Comment

Related Items