Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Monitoring nociception during general anesthesia with heart rate variability Brouse, Chris J. 2015

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

Item Metadata

Download

Media
24-ubc_2015_may_brouse_chris.pdf [ 3.32MB ]
Metadata
JSON: 24-1.0167194.json
JSON-LD: 24-1.0167194-ld.json
RDF/XML (Pretty): 24-1.0167194-rdf.xml
RDF/JSON: 24-1.0167194-rdf.json
Turtle: 24-1.0167194-turtle.txt
N-Triples: 24-1.0167194-rdf-ntriples.txt
Original Record: 24-1.0167194-source.json
Full Text
24-1.0167194-fulltext.txt
Citation
24-1.0167194.ris

Full Text

Monitoring Nociception DuringGeneral Anesthesia with Heart RateVariabilitybyChris J. BrouseB.A.Sc., The University of Ottawa, 2004M.A.Sc., The University of British Columbia, 2007A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Electrical & Computer Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)April 2015c Chris J. Brouse 2015AbstractWe have developed a novel real-time cardiorespiratory coherence (CRC)algorithm to monitor nociception during general anesthesia. We have madetwo novel and significant contributions to the field.We have developed a novel filter for measuring the respiratory sinusarrhythmia (RSA) in the heart rate variability (HRV) in real-time. Thefilter uses information from a secondary signal source (a respiration ratederived from a respiration signal) to track the RSA as it moves in timeand frequency. It then uses this information to dynamically tune the centrefrequency and bandwidth of a band-pass filter to isolate and measure theRSA. This filter is very e↵ective at tracking the RSA in time and frequency,and it may provide the most robust measure of RSA yet devised.We have integrated multiple signals and algorithms together into an end-to-end software system for robust continuous real-time nociception monitor-ing during general anesthesia. The software system incorporates not onlyour novel RSA filter to measure nociception, but also many peripheral algo-rithms for detecting and rejecting artifacts in the input signals. The inputsignals required for real-time nociception monitoring can be extremely noisy,and artifacts are a very significant challenge. To our knowledge, no othernociception monitor includes such robust artifact handling using redundantsignal channels.We estimated the sensitivity of our real-time CRC algorithm to noci-ception and antinociception, and compared it to traditional univariate HRVmeasures and standard clinical vital signs. Following ethics approval andinformed parental consent, data were collected from 48 children receivinggeneral anesthesia during dental surgery. A total of 42 dental dam insertion(nociception) and 57 anesthetic bolus (antinociception) events were noted.A nociception index was created for each HRV algorithm, ranging from 0(no nociception) to 100 (strong nociception). Dental dam insertion changedthe CRC nociception index by an average of 27 [95% CI from 21 to 33] (P< 0.000005), and a bolus dose of anesthetic changed it by an average of -19[-27 to -12] (P < 0.00003). Real-time CRC was more sensitive to nociceptionand antinociception than were the traditional measures.iiPrefaceChapter 2 is based on work conducted in UBC’s Electrical & ComputerEngineering for Medicine (ECEM) group by Chris J. Brouse, Dr. GuyA. Dumont, and Dr. J. Mark Ansermino. I was responsible for propos-ing cardiorespiratory coherence as a measure of nociception, experimentingwith wavelet transform coherence, and the development of the real-timecardiorespiratory coherence algorithm.Chapter 3 is based on work conducted in UBC’s Electrical & ComputerEngineering for Medicine (ECEM) group and the British Columbia Chil-dren’s Hospital’s Pediatric Anesthesia Research Team (PART), by Chris J.Brouse, Dorothy Myers, Erin Cooke, Jonathan Stinson, Joanne Lim, Dr.Guy A. Dumont, and Dr. J. Mark Ansermino. This research was approvedby the University of British Columbia and Children’s & Women’s HealthCentre of British Columbia Research Ethics Board. The project title is“Skin conductance fluctuations and heart rate variability as measures ofintraoperative nociception in children,” and the certificate number is H09-00905. I was responsible for part of the clinical study design, part of runningthe study and collecting data, and for the post hoc data analysis.Some aspects of Chapters 2 and 3 have been published:• Brouse CJ et al., “Real-Time Cardiorespiratory Coherence DetectsNociception During General Anesthesia,” IEEE TBME (in prepara-tion).• Brouse CJ et al., “Monitoring Nociception During General Anesthesiawith Cardiorespiratory Coherence,” JCMC 2013.• Brouse CJ et al., “Real-Time Cardiorespiratory Coherence is Blindto Changes in Respiration During General Anesthesia,” IEEE EMBC2013.• Brouse CJ et al., “Real-Time Cardiorespiratory Coherence DetectsAntinociception During General Anesthesia,” IEEE EMBC 2012.• Brouse CJ et al., “Measuring Adequacy of Analgesia with Cardiores-piratory Coherence,” STA 2012.iiiPreface• Brouse CJ et al., “Wavelet Transform Cardiorespiratory CoherenceDetects Patient Movement During General Anesthesia,” IEEE EMBC2011.• Brouse CJ et al., “Wavelet Transform Cardiorespiratory Coherence forMonitoring Nociception,” CinC 2010.I collected the clinical data along with Dorothy Myers, Erin Cooke, andJonathan Stinson. I performed all of the data analysis and wrote themanuscripts.Chapter 4 is based on work conducted in UBC’s Electrical & ComputerEngineering for Medicine (ECEM) group by Chris J. Brouse, Dr. Guy A.Dumont, and Dr. J. Mark Ansermino. I was responsible for proposing thenociception monitor architecture, experimenting with the individual dataanalysis algorithms, and the development of the real-time nociception mon-itoring system.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x1 Monitoring Nociception: A Review . . . . . . . . . . . . . . 11.1 Why Monitor Nociception? . . . . . . . . . . . . . . . . . . . 11.2 Heart Rate Variability . . . . . . . . . . . . . . . . . . . . . . 51.2.1 Time Domain Analysis . . . . . . . . . . . . . . . . . 71.2.2 Frequency Domain Analysis . . . . . . . . . . . . . . 71.2.3 Joint Time-Frequency Domain Analysis . . . . . . . . 91.2.4 Nonlinear Dynamics Analysis . . . . . . . . . . . . . 101.2.5 Problems with Time/Frequency HRV Theory . . . . . 101.3 Computational Eciency is Important . . . . . . . . . . . . 111.4 Existing Nociception Monitors . . . . . . . . . . . . . . . . . 131.5 Taylor’s Experiment & A New Understanding of RSA . . . . 141.6 Structure of the Dissertation . . . . . . . . . . . . . . . . . . 151.7 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . 162 Monitoring Nociception with Cardiorespiratory Coherence 192.1 A Review of Joint Time/Frequency Analysis . . . . . . . . . 192.1.1 The Heisenberg-Gabor Uncertainty Principle . . . . . 192.1.2 The Fourier Transform . . . . . . . . . . . . . . . . . 192.1.3 The Short Time Fourier Transform . . . . . . . . . . 212.1.4 The Gabor Transform . . . . . . . . . . . . . . . . . . 232.1.5 The Wigner-Ville Distribution . . . . . . . . . . . . . 232.1.6 The Discrete Wavelet Transform . . . . . . . . . . . . 24vTable of Contents2.1.7 The Continuous Wavelet Transform . . . . . . . . . . 252.1.8 The Stockwell Transform . . . . . . . . . . . . . . . . 262.2 What is Coherence? . . . . . . . . . . . . . . . . . . . . . . . 272.3 Wavelet Transform Cardiorespiratory Coherence . . . . . . . 302.3.1 Limitations of WTCRC . . . . . . . . . . . . . . . . . 322.4 Real-Time Cardiorespiratory Coherence . . . . . . . . . . . . 372.4.1 Analyzing Filter . . . . . . . . . . . . . . . . . . . . . 382.4.2 Smoothing Filter . . . . . . . . . . . . . . . . . . . . 412.4.3 Coherence Estimator . . . . . . . . . . . . . . . . . . 422.4.4 Real-Time Delay . . . . . . . . . . . . . . . . . . . . . 422.5 Deriving a Nociception Index . . . . . . . . . . . . . . . . . . 472.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 493 Clinical Validation . . . . . . . . . . . . . . . . . . . . . . . . . 513.1 Clinical Study of Nociception . . . . . . . . . . . . . . . . . . 513.1.1 Patients & Anesthesia . . . . . . . . . . . . . . . . . . 513.1.2 Pre-Induction Setup & Data Collection . . . . . . . . 523.1.3 Induction of Anesthesia & Sensor Attachment . . . . 543.1.4 Maintenance Phase of Anesthesia . . . . . . . . . . . 553.1.5 Emergence . . . . . . . . . . . . . . . . . . . . . . . . 563.2 CRC Analyzing Filter Tuning . . . . . . . . . . . . . . . . . 563.3 Nociception Index Tuning . . . . . . . . . . . . . . . . . . . . 563.3.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . 563.3.2 Derivation from a Bounded Measure . . . . . . . . . . 573.3.3 Derivation from an Unbounded Measure . . . . . . . 583.4 The Baseline-Stimulus-Response Experiment . . . . . . . . . 593.4.1 Data Segmentation . . . . . . . . . . . . . . . . . . . 593.4.2 Statistical Analysis . . . . . . . . . . . . . . . . . . . 593.5 Response to Dental Dam . . . . . . . . . . . . . . . . . . . . 603.5.1 Experimental Concept . . . . . . . . . . . . . . . . . 603.5.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . 613.6 Response to SNS . . . . . . . . . . . . . . . . . . . . . . . . . 623.6.1 Experimental Concept . . . . . . . . . . . . . . . . . 623.6.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . 673.7 Response to Bolus Dose of Anesthetic . . . . . . . . . . . . . 683.7.1 Experimental Concept . . . . . . . . . . . . . . . . . 683.7.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . 733.8 Response to Change in Respiration Rate . . . . . . . . . . . 743.8.1 Experimental Concept . . . . . . . . . . . . . . . . . 743.8.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . 79viTable of Contents3.9 Response to Change in Peak Airway Pressure . . . . . . . . 803.9.1 Experimental Concept . . . . . . . . . . . . . . . . . 803.9.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . 823.10 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 853.10.1 Major Findings . . . . . . . . . . . . . . . . . . . . . 853.10.2 Clinical Relevance . . . . . . . . . . . . . . . . . . . . 903.10.3 Limitations . . . . . . . . . . . . . . . . . . . . . . . . 913.11 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 924 A Real-Time Nociception Monitor . . . . . . . . . . . . . . . 944.1 Real-Time System Architecture . . . . . . . . . . . . . . . . 944.2 Real-Time Heart Rate Measurement . . . . . . . . . . . . . . 954.2.1 Real-Time HR from the ECG . . . . . . . . . . . . . 954.2.2 Real-Time PR from the PPG . . . . . . . . . . . . . . 984.3 Real-Time Resampling . . . . . . . . . . . . . . . . . . . . . 1004.3.1 HR Resampling . . . . . . . . . . . . . . . . . . . . . 1004.3.2 Respiration Resampling . . . . . . . . . . . . . . . . . 1024.3.3 RR Resampling . . . . . . . . . . . . . . . . . . . . . 1034.4 Real-Time Artifact Detection . . . . . . . . . . . . . . . . . . 1034.4.1 HR Artifact Detection . . . . . . . . . . . . . . . . . 1034.4.2 CO2 Artifact Detection . . . . . . . . . . . . . . . . . 1054.4.3 RR Artifact Detection . . . . . . . . . . . . . . . . . 1054.5 Real-Time Artifact Correction . . . . . . . . . . . . . . . . . 1074.5.1 HR Artifact Correction . . . . . . . . . . . . . . . . . 1074.5.2 CO2 Artifact Correction . . . . . . . . . . . . . . . . 1124.5.3 RR Artifact Correction . . . . . . . . . . . . . . . . . 1204.6 Real-Time Signal Source Selection . . . . . . . . . . . . . . . 1204.6.1 We Need Secondary Redundant Signal Sources . . . . 1204.6.2 Real-Time Signal Quality Estimation . . . . . . . . . 1224.6.3 Real-Time Signal Source Selection . . . . . . . . . . . 1244.6.4 The E↵ect of Sampling Frequency on HRV . . . . . . 1244.7 Real-Time Delay . . . . . . . . . . . . . . . . . . . . . . . . . 1314.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1325 Conclusion & Future Work . . . . . . . . . . . . . . . . . . . . 1355.1 Feedback Control of Antinociceptive Drugs . . . . . . . . . . 1365.2 Enhanced Nociception Monitoring with Sensor Fusion . . . . 1365.2.1 Combining HRV Measures . . . . . . . . . . . . . . . 1385.2.2 Cardiovascular Respiratory Coherence . . . . . . . . . 1385.2.3 End-Tidal CO2 . . . . . . . . . . . . . . . . . . . . . 141viiTable of Contents5.2.4 Respiration Rate . . . . . . . . . . . . . . . . . . . . . 1435.3 Further Validation . . . . . . . . . . . . . . . . . . . . . . . . 1445.3.1 Spontaneous Ventilation . . . . . . . . . . . . . . . . 1445.3.2 In Adults and in Disease States . . . . . . . . . . . . 1495.3.3 Sensitivity to Nociception at Di↵erent Levels of Antinoci-ception . . . . . . . . . . . . . . . . . . . . . . . . . . 1495.3.4 Respiration Rate Ramp Down . . . . . . . . . . . . . 1505.3.5 Outcomes Studies . . . . . . . . . . . . . . . . . . . . 151Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153viiiList of Tables3.1 Response to a dental dam insertion. Asterisks (*) denote sta-tistically significant results, factoring in the Bonferroni cor-rection. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 623.2 Response to an SNS event. Asterisks (*) denote statisticallysignificant results, factoring in the Bonferroni correction. . . . 683.3 Response to a bolus dose of anesthetic drugs. Asterisks (*)denote statistically significant results, factoring in the Bon-ferroni correction. . . . . . . . . . . . . . . . . . . . . . . . . . 743.4 Response to a change in RR from 8 to 16 breaths/min. As-terisks (*) denote statistically significant results, factoring inthe Bonferroni correction. . . . . . . . . . . . . . . . . . . . . 803.5 Response to a change in PPaw from 15 to 12 cmH2O. As-terisks (*) denote statistically significant results, factoring inthe Bonferroni correction. . . . . . . . . . . . . . . . . . . . . 854.1 Theoretical error in the HR and the HRV-based nociceptionindices arising from the sampling frequency of the source ECGsignal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1274.2 Actual error in the HR/PR and the HRV-based nociceptionindices arising from the source signal and its sampling frequency.128ixList of Figures1.1 Anatomy of the autonomic nervous system (ANS). The ANSa↵ects the behaviour of nearly every system in the body. . . . 31.2 Example series of instantaneous HR (thin blue line), used forHRV analysis, as compared to 10 s averaged HR (thick blackline). The rapid fluctuations around t = 50 s are missing en-tirely from the mean HR series. Even the slower fluctuationsaround t = 100 are not fully captured in the mean HR series. 61.3 HR series power spectrum under normal conditions (a) andunder autonomic blockade (b) (From [2]). . . . . . . . . . . . 81.4 The real-time CRC algorithm. The RR tunes the analyzingfilters to isolate the RSA while rejecting other HRV com-ponents. This dynamic RSA tracking algorithm is a novelcontribution to the field. . . . . . . . . . . . . . . . . . . . . . 171.5 The nociception monitor system architecture. Each box rep-resents an algorithm required for real-time monitoring of noci-ception. Many are peripheral algorithms used for pre-processingof the signals or artifact detection/rejection. The algorithmsare described in the sections referenced. This highly redun-dant and robust architecture is a novel contribution to thefield. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.1 Heisenbox representing the spread in time and frequency of aquantum of information in the time/frequency plane. . . . . . 202.2 Heisenboxes representing the Fourier transform’s localizationin time and frequency. . . . . . . . . . . . . . . . . . . . . . . 212.3 Heisenboxes representing the short time Fourier transform’slocalization in time and frequency. . . . . . . . . . . . . . . . 222.4 Filter bank implementation of the DWT. . . . . . . . . . . . 242.5 Heisenboxes representing the discrete wavelet transform’s lo-calization in time and frequency. . . . . . . . . . . . . . . . . 25xList of Figures2.6 Heisenboxes representing the continuous wavelet transform’slocalization in time and frequency. . . . . . . . . . . . . . . . 262.7 Example of coherent HR and respiration signals from the clin-ical study of nociception (see Section 3.1). Top plot: HR.Bottom plot: respiration. The light grey region denotes aninspiration phase, and the dark grey region denotes an ex-piration phase. Note that the sinusoidal pattern from therespiration signal is evident in the HR signal. The patient’srespiration is having a strong e↵ect on the HR. . . . . . . . . 282.8 Example of incoherent HR and respiration signals from theclinical study of nociception (see Section 3.1). Top plot: HR.Bottom plot: respiration. The light grey region denotes aninspiration phase, and the dark grey region denotes an expi-ration phase. Note that the sinusoidal pattern from the res-piration signal is not present in the HR signal. The patient’srespiration is having little e↵ect on the HR. . . . . . . . . . . 292.9 CWT decomposition scales used in the WTCRC analysis,and their corresponding Fourier frequencies. Only the circledpoints are analyzed. . . . . . . . . . . . . . . . . . . . . . . . 312.10 Example HR wavelet power spectrum, as an intermediate stepin WTCRC. Data are from the clinical study of nociception(see Section 3.1). Top plot: Heart rate. Bottom plot: HRpower spectrum, calculated with the CWT. Bright areas in-dicate strong power. Vertical dashed lines (- -) denote clinicalevents and the dot-dash line (- - .) denotes a rezeroing ar-tifact in the capnogram. The patient movement event is asign of strong nociception. Note the change in HR frequencycontent (a shift to a lower frequency) during the period ofnociception from t = 75 - 175. . . . . . . . . . . . . . . . . . . 33xiList of Figures2.11 Example respiration wavelet power spectrum, as an interme-diate step in WTCRC. Data are from the clinical study ofnociception (see Section 3.1). Top plot: Capnogram. Bot-tom plot: Capnogram power spectrum, calculated with theCWT. Bright areas indicate strong power. Vertical dashedlines (- -) denote clinical events and the dot-dash line (- - .)denotes a rezeroing artifact in the capnogram. The patientmovement event is a sign of strong nociception. Note the dis-turbances in respiration frequency content during the periodof nociception from t = 75 - 175. During this period, thepatient is resisting the mechanical ventilator. The dominantrespiration frequency is 0.2 Hz (12 breaths/min), and thereare harmonics at higher frequencies. . . . . . . . . . . . . . . 342.12 Example cross wavelet power spectrum, as an intermediatestep in WTCRC. Data are from the clinical study of noci-ception (see Section 3.1). Top plot: Heart rate. Middle plot:Capnogram. Bottom plot: Cross power spectrum, calculatedwith the CWT. Bright areas indicate strong power. Verticaldashed lines (- -) denote clinical events and the dot-dash line(- - .) denotes a rezeroing artifact in the capnogram. Thepatient movement event is a sign of strong nociception. Notethe disturbances in frequency content during the period ofnociception from t = 75 - 175. . . . . . . . . . . . . . . . . . . 352.13 Example anesthetic bolus event analysis with WTCRC. Dataare from the clinical study of nociception (see Section 3.1).Plots from top to bottom: a) Heart rate; b) Capnogram; c)Cardiorespiratory coherence map. Bright areas indicate highcoherence. Horizontal continuous line (green) indicates therespiratory frequency obtained from b); d) WTCRC nocicep-tion index. Vertical dashed lines (- -) denote clinical eventsand the dot-dash line (- - .) denotes a rezeroing artifact inthe capnogram. The patient movement event is a sign ofstrong nociception. The nociception index is high in the pe-riod preceding the bolus dose of anesthetic drugs (baselineperiod, highlighted in red), and low following it (response pe-riod, highlighted in green). The capnometer rezeroing artifactleads to a transient false increase in the nociception index. . . 362.14 The real-time CRC algorithm. The RR tunes the analyzingfilters to isolate the RSA while rejecting other HRV components. 38xiiList of Figures2.15 The complex Morlet analyzing filter in the time domain. Sam-pling frequency fs is 4 Hz, centre frequency fc is 0.25 Hz, andcuto↵ is at 2a. . . . . . . . . . . . . . . . . . . . . . . . . . . 392.16 The complex Morlet analyzing filter in the frequency domainat di↵erent centre freqencies. Cuto↵ is at 2a. . . . . . . . . . 402.17 Example analysis of the HR signal using the real-time CRCanalyzing filter. Data are from the clinical study of nocicep-tion (see Section 3.1). Top plot is the raw heart rate. Bottomplot is the filtered heart rate. The solid blue line is the realcomponent, and the dashed green line is the imaginary com-ponent. The missing data at the start of the filtered signal iscaused by the width of the analyzing filter. . . . . . . . . . . 412.18 Example analysis of the respiration signal using the real-timeCRC analyzing filter. Data are from the clinical study ofnociception (see Section 3.1). Top plot is the raw CO2 signal.Bottom plot is the filtered CO2. The solid blue line is thereal component, and the dashed green line is the imaginarycomponent. The missing data at the start of the filtered signalis caused by the width of the analyzing filter. . . . . . . . . . 422.19 Example signal powers from the real-time CRC analyzing fil-ter examples in Figs. 2.17 and 2.18. Top plot is the HRpower. Middle plot is the CO2 power. Bottom plot is thecross power. The solid blue line is the real component, andthe dashed green line is the imaginary component. The miss-ing data at the start of the signals is caused by the width ofthe analyzing filter. . . . . . . . . . . . . . . . . . . . . . . . . 432.20 Frequency response of the causal Gaussian smoothing filter atvarious smoothing levels. The full (noncausal) filter responseis also shown for comparison. Note that the causal filter haszero time delay, but is less ecient at rejecting high freqencies. 442.21 Characteristics of the real-time CRC analyzing filter, withcentre frequency of 0.15 Hz and sampling frequency of 2.5Hz. Top plot: analyzing filter in the time domain. Bottomplot: group delay of the analyzing filter. . . . . . . . . . . . . 452.22 Group delay of the real-time CRC analyzing filter across arange of respiration frequencies from 0.05 - 0.50 Hz (3 - 30breaths/min). . . . . . . . . . . . . . . . . . . . . . . . . . . . 46xiiiList of Figures2.23 Characteristics of the real-time CRC smoothing filter, withthe bandwidth at  = 15, the cuto↵ at 3, and samplingfrequency of 2.5 Hz. Top plot: smoothing filter in the timedomain. Bottom plot: group delay of the smoothing filter. . . 482.24 Example CRC and nociception index following from the ex-amples in Figs. 2.17, 2.18, and 2.19. Top plot is the CRC.Bottom plot is the CRC nociception index. The missing dataat the start of the signals is caused by the combined widthsof the analyzing and smoothing filters. Unlike in the previousexample figures, this width is not directly related to real-timedelay. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 493.1 The nociception study timeline. . . . . . . . . . . . . . . . . . 523.2 Two synchronization points are used to synchronize data fromthree monitors in the nociception study. Timeline is notdrawn to scale. . . . . . . . . . . . . . . . . . . . . . . . . . . 543.3 The overall MSE of CRC results, using a range of analyzingfilter quality values. The error is calculated relative to ananalyzing filter with quality of 5 a. . . . . . . . . . . . . . . 573.4 Data segment structure for a general BSR experiment. . . . . 593.5 Data segment structure for a dental dam event. . . . . . . . . 613.6 Example dental dam event, analyzed with real-time CRC.This analysis corresponds to event number 39 in Fig. 3.7.Vertical lines denote clinical events. The green (left) andred (right) vertical bands represent the baseline and responseperiods, respectively. The nociception index is low in theperiod preceding the dental dam, and high during it. Themissing CRC data at the beginning of the window correspondsto the combined length of the analyzing and smoothing filters. 633.7 CRC nociception index response for each dental dam event,sorted by response strength. Arrows indicate the direction ofthe response from the baseline period to the response period.Upward arrows (blue) indicate a change in the expected direc-tion (increase), and downward arrows (red) indicate a changein the opposite direction (decrease). . . . . . . . . . . . . . . 64xivList of Figures3.8 Boxplots of responses to dental dam insertion for all measures.The central red bar represents the median (second quartile).The box edges are the first and third quartiles. The whiskersextend 1.5x the interquartile range (IQR) beyond the boxedges. Points beyond the whiskers are drawn as red plussymbols (+), and may be considered outliers. Such possibleoutliers were not excluded from analysis. . . . . . . . . . . . . 653.9 Comparison of mean nociception index responses to a dentaldam insertion. Longer bars are better. . . . . . . . . . . . . . 663.10 Data segment structure for an SNS event. . . . . . . . . . . . 673.11 Example SNS event, analyzed with real-time CRC. This anal-ysis corresponds to event number 107 in Fig. 3.12. Verticallines denote clinical events. The green (left) and red (right)vertical bands represent the baseline and response periods,respectively. The nociception index is low in the period pre-ceding the SNS, and high during it. The missing CRC dataat the beginning of the window corresponds to the combinedlength of the analyzing and smoothing filters. . . . . . . . . . 693.12 CRC nociception index response for each SNS event, sortedby response strength. Arrows indicate the direction of the re-sponse from the baseline period to the response period. Up-ward arrows (blue) indicate a change in the expected direction(increase), and downward arrows (red) indicate a change inthe opposite direction (decrease). . . . . . . . . . . . . . . . . 703.13 Boxplots of responses to the SNS for all measures. The cen-tral red bar represents the median (second quartile). Thebox edges are the first and third quartiles. The whiskers ex-tend 1.5x the interquartile range (IQR) beyond the box edges.Points beyond the whiskers are drawn as red plus symbols(+), and may be considered outliers. Such possible outlierswere not excluded from analysis. . . . . . . . . . . . . . . . . 713.14 Comparison of mean nociception index responses to the SNS.Longer bars are better. . . . . . . . . . . . . . . . . . . . . . . 723.15 Data segment structure for an anesthetic bolus event. . . . . 73xvList of Figures3.16 Example anesthetic bolus event, analyzed with real-time CRC.This analysis corresponds to event number 60 in Fig. 3.17.Vertical dashed lines (- -) denote clinical events and the dot-dash line (- .) denotes a rezeroing artifact in the capnogram.The patient movement event is a sign of strong nociception.The nociception index is high in the period preceding thebolus dose of anesthetic drugs, and low following it. The cap-nometer rezeroing artifact leads to a transient false increasein the nociception index. The missing CRC data at the be-ginning of the window corresponds to the combined length ofthe analyzing and smoothing filters. . . . . . . . . . . . . . . 753.17 CRC nociception index response for each anesthetic bolusevent, sorted by response strength. Arrows indicate the direc-tion of the response from the baseline period to the responseperiod. Downward arrows (blue) indicate a change in the ex-pected direction (decrease), and upward arrows (red) indicatea change in the opposite direction (increase). . . . . . . . . . 763.18 Boxplots of responses to a bolus dose of anesthetic drugs forall measures. The central red bar represents the median (sec-ond quartile). The box edges are the first and third quartiles.The whiskers extend 1.5x the interquartile range (IQR) be-yond the box edges. Points beyond the whiskers are drawn asred plus symbols (+), and may be considered outliers. Suchpossible outliers were not excluded from analysis. . . . . . . . 773.19 Comparison of mean nociception index responses to a bolusdose of anesthetic drugs. Longer bars are better. . . . . . . . 783.20 Data segment structure for a change in respiration rate event. 793.21 Example change in RR event, analyzed with real-time CRC.This analysis corresponds to event number 19 in Fig. 3.22.Vertical lines denote clinical events. The yellow vertical bands(both left and right) represent the baseline and response pe-riods, respectively. The nociception index is approximatelythe same in the period preceding the change in RR as it isfollowing it. The transient increase in the nociception in-dex immediately following the change in RR is caused by adelay in the RR estimation. For this brief period, the real-time CRC algorithm is looking for the RSA at the wrongfrequency. The missing CRC data at the beginning of thewindow corresponds to the combined length of the analyzingand smoothing filters. . . . . . . . . . . . . . . . . . . . . . . 81xviList of Figures3.22 CRC nociception index response for each change in RR event,sorted by response strength. Arrows indicate the direction ofthe response from the baseline period to the response period. 823.23 Boxplots of responses to a change in RR for all measures.The central red bar represents the median (second quartile).The box edges are the first and third quartiles. The whiskersextend 1.5x the interquartile range (IQR) beyond the boxedges. Points beyond the whiskers are drawn as red plussymbols (+), and may be considered outliers. Such possibleoutliers were not excluded from analysis. . . . . . . . . . . . . 833.24 Comparison of mean nociception index responses to a changein RR. Shorter bars are better. . . . . . . . . . . . . . . . . . 843.25 Data segment structure for a change in peak airway pressureevent. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 843.26 Example change in PPaw event, analyzed with real-time CRC.This analysis corresponds to event number 16 in Fig. 3.27.Vertical lines denote clinical events. The yellow vertical bands(both left and right) represent the baseline and response pe-riods, respectively. The nociception index is approximatelythe same in the period preceding the change in PPaw as it isfollowing it. The missing CRC data at the beginning of thewindow corresponds to the combined length of the analyzingand smoothing filters. . . . . . . . . . . . . . . . . . . . . . . 863.27 CRC nociception index response for each change in PPawevent, sorted by response strength. Arrows indicate the direc-tion of the response from the baseline period to the responseperiod. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 873.28 Boxplots of responses to a change in PPaw for all measures.The central red bar represents the median (second quartile).The box edges are the first and third quartiles. The whiskersextend 1.5x the interquartile range (IQR) beyond the boxedges. Points beyond the whiskers are drawn as red plussymbols (+), and may be considered outliers. Such possibleoutliers were not excluded from analysis. . . . . . . . . . . . . 883.29 Comparison of mean nociception index responses to a changein PPaw. Shorter bars are better. . . . . . . . . . . . . . . . . 89xviiList of Figures4.1 The nociception monitor system architecture. Each algorithmis represented as a box, and is implemented as a plugin insidethe iAssist software framework. The algorithms are describedin the sections referenced. . . . . . . . . . . . . . . . . . . . . 964.2 Example analysis with the SimpleDi↵ algorithm. The topplot shows a somewhat noisy ECG, with significant baselinewander. The middle plot shows how the noise has been re-jected, leaving only the relevant high frequency content. Thebottom plot shows the squared di↵erence and R-peak detec-tion threshold (the horizontal red dotted line). . . . . . . . . 974.3 Example analysis with the SimpleZeroCrossings algorithm.The PPG signal has been heavily filtered by the Masimooximeter prior to collection, allowing us to detect pulses assimple zero crossings (at the horizontal red dotted line). . . . 994.4 Example HR resampling condition, where the local windowfalls entirely within a single heart period. . . . . . . . . . . . 1014.5 Example HR resampling condition, where the local windowspans two heart periods. This condition can be extended tohandle the condition where the local window spans multipleheart periods. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1024.6 Example of CO2 downsampling. The top plot is the originalCO2 signal at 300 Hz (note that the source was recorded at 25Hz, but the S/5 software upsamples to 300 Hz). The CO2 sig-nal is somewhat nonstationary, as the patient is fighting theventilator. The middle plot is the CO2 signal resampled to2.5 Hz using lowpass filtering and decimation (calculated withthe resample function in Matlab). The bottom plot is theCO2 signal resampled to 2.5 Hz using decimation only (cal-culated with the downsample function in Matlab. Note thatthe lowpass filtering introduces small distortions, unlike thesimple decimation-only approach. Simple decimation outper-forms the standard approach, even under the nonstationaryconditions shown here. Regardless, the di↵erence between thetwo results is small. . . . . . . . . . . . . . . . . . . . . . . . . 104xviiiList of Figures4.7 Example analysis with the DetectHRArtifacts SimpleDiffalgorithm. The top plot shows a HR with three typical arti-facts. The first artifact (left) results from a false beat, second(middle) results from a misplaced beat, and third (right) re-sults from a missed beat. Since the artifacts are so large inmagnitude, it is dicult to discern other details of the HR.The same HR is shown later in Fig. 4.9, with full details vis-ible. The middle plot shows how the artifacts have been iso-lated, leaving only the relevant high frequency content. Thebottom plot shows the squared di↵erence and HR artifact de-tection threshold (the horizontal red dotted line). Verticalred bars denote the detected artifact regions. . . . . . . . . . 1064.8 Example of an HR artifact spreading out in time during aCRC analysis. The top plot shows an HR with an artifactwith width 2 s. The middle plot shows the HR after be-ing processed with the CRC analyzing filter. The solid blueline is the real component, and the dashed green line is theimaginary component. The bottom plot shows the CRC out-put after the HR has been processed with the analyzing andsmoothing filters. The red regions approximate the artifactwidth. The artifact width increases with each stage of pro-cessing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1084.9 Example cubic interpolation of HR artifacts. The top plotshows a gold standard HR without artifacts. The second plotshows the same HR, with simulated artifacts removed. Theseare the same artifacts shown earlier in Fig. 4.7. The thirdplot shows the HR with artifacts corrected using cubic inter-polation. The bottom plot shows the HR di↵erence betweenthe gold standard and corrected HRs. Vertical red bars de-note the artifact regions. . . . . . . . . . . . . . . . . . . . . . 1114.10 HR interpolation peak error results. Each thin line representsthe peak error at a range of artifact lengths for a single HRsegment. The thick red line represents the mean peak error,averaged across all segments. As the artifact width increases,the peak error also increases. . . . . . . . . . . . . . . . . . . 113xixList of Figures4.11 CRC peak error results from HR interpolation. Each thinline represents the peak error at a range of artifact lengthsfor a single CRC segment. The thick red line represents themean peak error, averaged across all segments. As the artifactwidth increases, the peak error also increases. The verticalscaling encompasses the entire possible range of the CRC no-ciception index. . . . . . . . . . . . . . . . . . . . . . . . . . . 1144.12 CRC error results in time from HR interpolation. Each thinline represents the mean error across all CRC segments at asingle artifact length, over time. The HR artifact begins at t= 0 (far left of the plot). As the artifact width increases, theerror lines climb higher (greater error). The vertical scalingencompasses the entire possible range of the CRC nociceptionindex. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1154.13 Example repeater interpolation of a CO2 rezeroing artifact.The top plot shows a gold standard CO2 signal without ar-tifacts. The second plot shows the same CO2 signal, with asimulated rezeroing artifact lasting 10 s. The third plot showsthe CO2 with artifact corrected using a signal repeater. Thebottom plot shows the CO2 di↵erence between the gold stan-dard and corrected CO2 signals. The vertical red bar denotesthe artifact region. . . . . . . . . . . . . . . . . . . . . . . . . 1174.14 CO2 interpolation peak error results. Each blue circle repre-sents the peak error for a single CO2 segment. The horizon-tal dotted red line represents the mean peak error, averagedacross all segments. . . . . . . . . . . . . . . . . . . . . . . . . 1184.15 CRC error results in time from CO2 interpolation. The linerepresents the mean error across all CRC segments, over time.The CO2 artifact begins at t = 0 (far left of the plot). Thevertical scaling encompasses the entire possible range of theCRC nociception index. . . . . . . . . . . . . . . . . . . . . . 1194.16 Example of electrocautery noise corrupting an ECG signal.Electrocautery noise is present for the first 3 s of the signal. 1224.17 Example of rezeroing in the capnometer (top plot, from 35- 50 s). The CO2 and O2 signals are temporarily unavail-able. The spirometry signals (bottom plot) are still availablethroughout. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1234.18 Example of HR series derived from a NeuroSense ECG at 900Hz, a Datex/Ohmeda ECG at 300 Hz, and a Datex/OhmedaPPG at 100 Hz. . . . . . . . . . . . . . . . . . . . . . . . . . . 128xxList of Figures4.19 Example of sampling frequency and source signal e↵ect onthe HRV-based nociception indices. Top plot: HR; Secondplot: SDNN nociception index; Third plot: LF/HF nocicep-tion index; Bottom plot: CRC nociception index. In all plots,the blue line is generated from the NeuroSense ECG at 900Hz, the green line from the Datex/Ohmeda ECG at 300 Hz,and the red line from the Datex/Ohmeda PPG at 100 Hz.Note that in all plots the blue and green lines are very closetogether (small error) but the red line is farther away (largererror). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1294.20 The real-time nociception monitor, pictured during use inthe operating room at British Columbia’s Children’s Hospital.The nociception monitor is the blue outlined screen next tothe anesthesiologist at right. . . . . . . . . . . . . . . . . . . . 1345.1 Boxplots of CRC during the expected “no nociception” pe-riods of each experiment from Chapter 3. The central redbar represents the median (second quartile). The box edgesare the first and third quartiles. The whiskers extend 1.5xthe interquartile range (IQR) beyond the box edges. Pointsbeyond the whiskers are drawn as red plus symbols (+), andmay be considered outliers. Such possible outliers were notexcluded from analysis. . . . . . . . . . . . . . . . . . . . . . 1375.2 Example PPG signal for pulse oximetry. Two light wave-lengths (red light at 660 nm and infrared light at 880 nm) areused. Rapid oscillations represent arterial blood pulsationscaused by heart beats. Slow oscillations (baseline wander)are caused by respiration. . . . . . . . . . . . . . . . . . . . . 1405.3 Example response to a dental dam insertion event, analyzedwith real-time CRC, and showing the time-varying EtCO2.The plot of EtCO2 is second from the bottom. Vertical linesdenote clinical events. The level of nociception increases fol-lowing dental dam insertion, as reflected in the increase inboth the CRC nociception index and in the EtCO2. . . . . . 142xxiList of Figures5.4 Example response to cold pressor event, analyzed with real-time CRC. Vertical lines denote clinical events. The green andred vertical bands represent the baseline and response peri-ods, respectively. The nociception index increases strongly(from 50.7 to 78.3) after the subject places his/her hand inthe ice water. The missing CRC data at the beginning of thewindow corresponds to the combined length of the analyzingand smoothing filters. Missing data in the middle of the plotcorresponds to a temporary cessation of respiration, when theRR could not be accurately calculated. . . . . . . . . . . . . . 1455.5 Example response to cold pressor event, analyzed with real-time CRC. Vertical lines denote clinical events. The greenand red vertical bands represent the baseline and responseperiods, respectively. The nociception index unexpectedlydecreases (from 78.3 to 52.8) after the subject places his/herhand in the ice water. The missing CRC data at the be-ginning of the window corresponds to the combined lengthof the analyzing and smoothing filters. Missing data in themiddle of the plot corresponds to a temporary cessation ofrespiration, when the RR could not be accurately calculated. 1465.6 Example response to cold pressor event, analyzed with real-time CRC. Vertical lines denote clinical events. The greenand red vertical bands represent the baseline and responseperiods, respectively. The nociception index increases slightly(from 83.9 to 87.3) after the subject places his/her hand inthe ice water; however, the index is very high throughout.The missing CRC data at the beginning of the window corre-sponds to the combined length of the analyzing and smooth-ing filters. Missing data in the middle of the plot correspondsto a temporary cessation of respiration, when the RR couldnot be accurately calculated. . . . . . . . . . . . . . . . . . . 147xxiiList of AbbreviationsAbbreviation DescriptionANS Autonomic Nervous SystemASA American Society of AnesthesiologistsATP Adenosine TriphosphateBSR Experiment Baseline-Stimulus-Response ExperimentCI Confidence IntervalCRC Cardiorespiratory CoherenceCVRC Cardiovascular Respiratory CoherenceCWT Continuous Wavelet TransformDFA Detrended Fluctuation AnalysisECG ElectrocardiogramED Exponential DistributionEEG ElectroencephalogramESU Electrosurgical UnitET Endotracheal TubeEtCO2 End-Tidal CO2FFT Fast Fourier TransformFT Fourier TransformFTLE Finite Time Lyapunov ExponentsHF Band High Frequency BandHR Heart RateHRV Heart Rate VariabilityIQR Inter-Quartile RangeJTF Joint Time/FrequencyLF Band Low Frequency BandLTI Linear Time InvariantMSE Mean Squared ErrorNI Nociception IndexNIBP Noninvasive Blood PressureNMB Neuromuscular BlockadeNN Interval Normal-to-Normal IntervalxxiiiList of AbbreviationsOR Operating RoomPC Personal ComputerPONV Post-Operative Nausea and VomitingPPG PhotoplethysmogramPPM Physiological Patient MonitorPR Pulse RateRMSE Root Mean Squared ErrorRR Respiration RateRSA Respiratory Sinus ArrhythmiaSA Node Sino-Atrial NodeSDNN Standard Deviation of NN IntervalsSNS Standardized Nociceptive StimulusSQI Signal Quality IndexTCI Target Controlled InfusionWT Wavelet TransformWTC Wavelet Transform CoherenceWVD Wigner-Ville DistributionxxivChapter 1Monitoring Nociception: AReview1.1 Why Monitor Nociception?Nociception is the physiological perception of noxious stimuli. Noxious stim-uli may be thermal (extreme heat or cold), chemical (e.g. environmental irri-tants), or mechanical (extreme pressure or deformation). Nociception is de-tected by specialized nerve endings called nociceptors, which exist through-out most of the body. Nociceptors transmit a signal to the brain when theydetect a noxious stimulus. The nociception signal can alter the state of theautonomic nervous system (ANS) by acting on lower brain regions, mostnotably the brainstem. In conscious subjects, the signal is also relayed viathe thalamus to the neocortex, and is experienced as pain. Pain is entirelya conscious experience. However, nociception can still a↵ect the body evenin the absence of consciousness.Nociception evokes a physiological stress response. The ANS increasessympathetic tone and decreases parasympathetic tone. Parasympathetictone is conveyed by the vagus nerve, which acts on the body via acetyl-choline. Sympathetic tone is conveyed by sympathetic nerves, which acton the body via norepinephrine. Sympathetic tone is also conveyed by theendocrine system, which acts on the body via epinephrine in the blood-stream. The sympathetic nervous system invokes the body’s stress response(so-called “fight-or-flight” response), while the parasympathetic nervous sys-tem invokes the body’s relaxation response (so-called “rest-and-digest” re-sponse). The two ANS branches operate in complementary opposition toeach other. Their relative action constitutes the body’s autonomic state,also called the autonomic tone or autonomic balance. The ANS state dy-namically adjusts to maintain homeostatic balance. As external or internalconditions change, the ANS adapts to compensate. Nociception shifts thebalance towards a more sympathetic stress state.The sympathetic response a↵ects nearly every system in the body (see11.1. Why Monitor Nociception?Fig. 1.1). The heart beats faster, becomes more contractile, and the cardiacoutput increases. The vasculature constricts, pulling the blood away fromthe periphery and towards the core, increasing blood pressure. Digestion isslowed as blood is redistributed towards the skeletal muscles. The pupilsdilate and skin perspiration increases. Conscious attention and awareness isnarrowed and focused. In this state, the body is primed for “fight-or-flight”actions. However, during surgery this state may be maladaptive.The anesthesiologist administers drugs to prevent awareness duringsurgery. Awareness may be conscious of unconscious (physiological). Gen-eral anesthesia includes immobility (muscle relaxation), hypnosis (uncon-sciousness) and antinociception. Hypnotic drugs, such as propofol, suppressconscious awareness and memory of the surgery. Research suggests thatthe hypnotic agents are very e↵ective, with only 0.007 - 0.15% (between1 in 14,560 and 1 in 655) of patients reporting any conscious awarenessduring surgery [48, 57, 61, 63]. Nevertheless, devices that monitor depthof hypnosis are becoming part of routine practice to further reduce intra-operative conscious awareness. Such devices include the Bispectral Index(BIS) (Aspect Medical Systems Inc., Newton, MA, USA), Entropy (GEHealthcare, Helsinki, Finland), and NeuroSense (NeuroWave Systems Inc.,Cleveland Heights, OH, USA) monitors. The anesthesiologist also adminis-ters antinociceptive drugs to suppress the body’s physiological awareness ofnociception, which can lead to a sympathetic stress response. Antinocicep-tive drugs suppress both the nociceptive a↵erent signals (signals sent fromthe nociceptors to the brain) and the sympathetic e↵erent signals (signalssent from the brain to the body). In so doing, antinociceptive drugs helpmaintain hemodynamic stability. Morphine, fentanyl, and remifentanil areexamples of popular modern antinociceptive drugs. It is common to usepropofol and remifentanil together, as their e↵ects have been found to besynergistic. That is, when used together, less is required of each drug thanwould be required separately. Furthermore, an increase in propofol willhave antinociceptive properties, and an increase in remifentanil will alsohave hypnotic properties [51].Delivering the right amount of antinociceptive drugs is more art thanscience. Drug dosages are usually calculated based on patient weight, butthere is a huge variation in individual response to drugs. Sex, genetics, andenvironment all contribute to the variation in response [66]. Some groupshave genetic variations that predispose them to require more or less of certaindrugs. For example, redheads may be both more sensitive to certain typesof nociception, and less sensitive to certain antinociceptive drugs [40, 41].Environmental di↵erences can also a↵ect the required drug dosage. For21.1. Why Monitor Nociception?HeartKidneysStomachLiverIntestinesLungsArteriesSweat GlandsSympathetic Nervous SystemParasympathetic Nervous SystemPupilsAutonomicNervous SystemFigure 1.1: Anatomy of the autonomic nervous system (ANS). The ANSa↵ects the behaviour of nearly every system in the body.31.1. Why Monitor Nociception?example, chronic opioid drug users acquire a tolerance to the drugs’ e↵ects,leading them to require higher doses of opioid antinociceptive drugs. Furthercomplicating matters, the level of surgical stress changes over the course ofsurgery, and the level of drugs should change in step with the patient’sresponse to the stress. It is dicult to predict both the level of nociceptionat any point in time, as well as how the individual patient will respond to agiven dose of drugs, and thus how much anesthetic is required at any timeduring surgery. Patients may be underdosed or overdosed during surgery.Underdosing and overdosing of antinociceptive drugs are both danger-ous. An underdose leads to inadequately suppressed nociception, whichresults in an excessive and/or prolonged stress response, and possible phys-iological damage. Some patient populations are particularly vulnerable tosuch damage, including the very elderly and those with compromised car-diovascular systems. Studies have linked surgically-induced stress hormoneswith perioperative myocardial infarction in this population [21, 39]. It isbelieved that stress-induced coronary spasm may increase shear stress oncoronary plaques, leading to their disruption and thus to myocardial infarc-tion. In contrast, an overdose of antinociceptive drugs can induce apneaand can lead to post-operative nausea and vomiting (PONV). PONV is oneof the most common complications of surgery, leading to delayed recovery,increased hospital stays, and additional nursing time [30]. PONV can be asignificant contributor to the cost of health care delivery. Reducing the doseof anesthetic agents can reduce the risk of PONV; however, optimal anes-thesia must still be the highest priority [28]. Optimal anesthesia may alsoreduce the risk of postoperative delirium (POD) and postoperative cognitivedysfunction (POCD). POD has been associated with delayed recovery andincreased hospital stays, while POCD has been associated with increasedmortality and long-term cognitive disability [22]. Both the surgical stressresponse and anesthetic drugs have been found to contribute to POD [32].A recent study has shown that BIS-guided control of hypnosis may decreasePOD and POCD [17]. Similar guidance for antinociception may decreasethem further.A direct, real-time nociception monitor would be highly beneficial tothe field of anesthesia. It would eliminate any uncertainty in the drug doselevels. A nociception monitor would ensure the patient receives the precisedrug dose required for his/her physiology and for the current level of surgicalstress. It would reduce the risk of overdosing or underdosing. A nociceptionmonitor would thus increase patient safety during general anesthesia. Heartrate variability (HRV) shows promise for monitoring nociception.41.2. Heart Rate Variability1.2 Heart Rate VariabilityThere are many di↵erent possible measures of nociception, but HRV hasbeen directly associated with poor health and negative clinical outcomes.Over the last 30 years, various features of HRV have been associated withcardiovascular disease, sepsis, diabetes, Alzheimer’s disease, emotional dis-turbances, arthritis, and even periodontal disease [73]. Research has shownthat HRV is a significant independent predictor of mortality in post-myocardial infarction patients [35]. There is growing evidence linking HRVmeasures to sudden cardiac death [37], and to prolonged postoperative my-ocardial ischemia [38].HRV refers to fluctuations in the instantaneous heart rate (HR) overtime. Most of these fluctuations are too fast to be observed in the standard5 - 10 s HR trends reported by physiological monitors. The instantaneousHR (for HRV) contains important information that is missing from the meanHR. Figure 1.2 illustrates the di↵erence between the instantaneous HR andthe mean HR.HRV reflects the state of the ANS. The sino-atrial (SA) node acts as theheart’s natural pacemaker. Cells in the SA node have an unstable restingmembrane potential, causing them to gradually depolarize in the absenceof any external stimulus. When the membrane potential reaches a thresh-old, the cells fire an action potential, which triggers a heart beat. The SAnode cells then rapidly repolarize, and restart the process of gradual depo-larization. The sympathetic nervous system releases norepinephrine nearthe SA node cells, which acts to accelerate the depolarization. In contrast,the parasympathetic nervous system releases acetylcholine near the SA nodecells, which acts to decelerate the depolarization. The time between heartbeats thus reflects the SA node’s intrinsic firing rate, as well as the inputsfrom the sympathetic and parasympathetic nervous systems.All HRV analyses operate on normal-to-normal (NN) heartbeat intervals.NN intervals are those between adjacent heartbeats that result from SAnode depolarization. Beats that result from other causes (e.g. prematureventricular complexes) are not used in HRV analysis because they do notreflect the autonomic state.Dozens of mathematical techniques have been applied to HRV analysis[71]. Most are minor variants on four major method themes: analyses in thetime domain, frequency domain, joint time-frequency domain, and nonlineardynamics. This section explores these themes and some of their comprisingmethods.51.2. Heart Rate Variability0 50 100 150 200 250 300707580859095100105110115Time (s)Heart Rate (bpm)Figure 1.2: Example series of instantaneous HR (thin blue line), used forHRV analysis, as compared to 10 s averaged HR (thick black line). Therapid fluctuations around t = 50 s are missing entirely from the mean HRseries. Even the slower fluctuations around t = 100 are not fully capturedin the mean HR series.61.2. Heart Rate Variability1.2.1 Time Domain AnalysisThe earliest techniques analyze the HR solely in the time domain. Thesetechniques are mainly statistical, analyzing the NN intervals or the di↵erencebetween NN intervals over a sliding window.The standard deviation of NN intervals (SDNN) is one of the most widelyused measures of HRV [71]. SDNN estimates the overall variability of theHR. As the autonomic balance shifts toward a more sympathetic state, theSDNN decreases. The opposite is true for more parasympathetic states.SDNN is typically used in analysis windows 5 mins long, however, it shouldcontinue to function when analyzed in shorter windows.Another common measure is the root mean square of successive di↵er-ences between NNs (RMSSD) [71]. By calculating the beat-to-beat changesin HR, RMSSD estimates the short term high frequency variation in theHR.Time domain techniques such as SDNN and RMSSD have been criticizedfor being too global; they fail to capture the underlying structure of theHR [45]. More advanced techniques have been developed to address thislimitation.A relatively new point process modeling approach is particularly inter-esting [4]. This method models the HR samples as random variables drawnfrom a statistical point process. The random variables are modeled usingthe Inverse Gaussian distribution, whose parameters are estimated from therecent history of beat times. The parameters are continually updated withnew estimations each time a new beat occurs. This method is interestingfor two reasons. First, it converts a discrete, unevenly sampled signal (theHR series) into a continuous time signal with arbitrary resolution. That is,it can estimate HRV parameters between beats. It may thus be useful forultra short term HRV analysis. Second, it can be used to predict the timeof the next heart beat. Future beat prediction could be used to help rejectHR artifacts, and to detect abrupt changes in autonomic state. Unfortu-nately this point process modeling approach is somewhat computationallycomplex, requiring a convex optimization procedure (e.g. Newton-Raphson)to be run after each new heart beat arrives. It may be dicult to implementthis method in real-time.1.2.2 Frequency Domain AnalysisFrequency domain techniques reveal an underlying structure to the HR se-ries. The Fourier transform (FT) transforms the HR series into the frequency71.2. Heart Rate Variability0.020.010.1 0.3 0.50.2 0.40Frequency (Hz)Power Spectrum of HR0.020.010.1 0.3 0.50.2 0.40Frequency (Hz)Autonomic BlockadeHigh frequencypeakLow frequencypeakVery low frequencypeakFigure 1.3: HR series power spectrum under normal conditions (a) andunder autonomic blockade (b) (From [2]).domain, and the resulting power spectrum can be divided into di↵erent fre-quency bands. The power in these bands has been shown to reflect theautonomic state [2] (Fig. 1.3).Traditional theory states that each ANS branch drives power in a cor-responding frequency band. The high frequency band (HF, 0.15-0.4 Hz)is caused by respiratory sinus arrhythmia (RSA) and is believed to reflectparasympathetic tone. The low frequency band (LF, 0.04-0.15 Hz) is mainlycaused by the baroreflex, and is believed to reflect sympathetic tone. A shiftin balance to a more sympathetic state would thus cause a decrease in HFpower and an increase in LF power. A shift to a more parasympathetic statewould cause the opposite e↵ect. The frequency bands can be measured inraw power (in ms2) or in normalized units (LFnu and HFnu), in which theLF or HF band power is divided by the total power from 0.04 Hz and up.The LF/HF power ratio is the most common measure of autonomicbalance [71]. It is believed to calculate the ratio of sympathetic tone toparasympathetic tone. The power ratio is often normalized to make it blindto changes in overall power [45]. Larger power ratios reflect sympatheticdominant states, while smaller power ratios reflect parasympathetic domi-nant states.The FT-based power ratio has some limitations, however. The FT re-quires evenly sampled data, but the HR series is unevenly sampled. A simpleand fast algorithm exists for resampling the HR series onto an even grid [6].81.2. Heart Rate VariabilityWe refer to this algorithm herein as “Berger’s algorithm,” named after itsauthor. The HR is typically resampled to 4 Hz. Alternatively, the powerspectrum of an unevenly sampled signal can be estimated directly using astatistical least squares technique called Lomb-Scargle [44, 62]. However,Lomb-Scargle can be very computationally intensive.Furthermore, frequency domain analysis requires long data windows andassumes stationary time series. The FT basis function, the complex expo-nential, is highly localized in frequency but poorly localized in time. As aresult, the FT works best with long input time series. HR windows of at least60 s are recommended to produce accurate power spectra. If the windowsare reduced much further, the frequency domain measure will have insuf-ficient resolution to accurately separate the LF and HF frequency bands.The FT assumes that the input signal is stationary over the analysis win-dow, but in reality the HR is highly nonstationary. These shortcomingslimit the usefulness of frequency analysis in nociception monitoring.1.2.3 Joint Time-Frequency Domain AnalysisJoint time-frequency (JTF) domain techniques reduce the stationarity as-sumption, and can operate on much shorter analysis windows. They achievethis by providing tunable localization in both time and frequency. JTFtechniques use variable basis functions that can shift their energy betweenthe time and frequency domains. In so doing, they produce power spectranearly as precise as the FT, but with a vastly improved time localization.JTF techniques can respond much faster to changes in autonomic balance.Several JTF techniques have been successfully applied to HRV analysis.The wavelet transform (WT) has been shown to produce results compara-ble to the FT, yet can do so over much shorter analysis windows [5]. TheWigner-Ville and Exponential distributions (WVD & ED) have also shownpromising results. These treat the HR series as a stochastic process, andcalculate its changing power spectral density over time [3, 78]. WVD andED exhibit some problems with cross-term interference, and are more com-putationally intensive than the WT [25].JTF techniques are typically applied in much the same way as FT anal-ysis. The HR series must first be resampled onto an even grid, using e.g.Berger’s algorithm [6]. The power in the LF and HF bands is then calcu-lated, and the LF/HF power ratio is used to quantify autonomic balance.91.2. Heart Rate Variability1.2.4 Nonlinear Dynamics AnalysisNonlinear dynamics analyses detect chaos in the HR. It has been proposedthat a healthy cardiovascular system is chaotic, and that a decrease in non-linear variability is associated with cardiovascular diseases [23]. A chaoticHR produces a fractal (self-similar) HR series.Many nonlinear dynamics methods have been applied to HRV analysis,often achieving greater sensitivity than FT and JTF methods. Approximateentropy (ApEn) estimates the HR’s entropy [55], which is a measure of thesignal’s randomness, unpredictability, or information content. It has beenshown to decrease significantly after fentanyl-based induction of anesthesia[68] and in children experiencing painful stimulation during propofol anes-thesia [75]. Finite time Lyapunov exponents (FTLE) quantifies chaos fromthe HR’s phase portrait (a graphical map of the changing system statesover time). It has been found to predict sudden cardiac death in implanteddefibrillator patients with greater accuracy than frequency analysis [79]. De-trended fluctuation analysis (DFA) detects fractal scaling exponents in theHR. In sleep studies, DFA has been found to better detect changes in HRVthan frequency analysis [54]. Similar multifractal analyses can be performedwith the WT. Nonlinear methods are often considered to give deeper insightinto the ANS than existing linear methods.Nonlinear dynamics analyses require long data sets for analysis, and canbe very computationally intensive [19]. These requirements may precludetheir use for nociception monitoring, particularly in short duration surgicalprocedures.1.2.5 Problems with Time/Frequency HRV TheoryThe underlying theory of the LF/HF power ratio has questionable validity.Some believe that the LF band power is a↵ected solely by sympathetictone, while others believe it is a↵ected by a combination of sympatheticand parasympathetic tone [7]. Most researchers believe that the HF bandpower is a↵ected solely by parasympathetic tone [7]. Some research has castdoubt on this simple view of the HF band, however. A study of respiration’se↵ect on HRV under sympathetic blockade suggests that RSA is driven byparasympathetic tone and restrained by sympathetic tone. The observed HFband power reflects the net balance of the two ANS branches [72]. Indeed,it may be impossible to separate the contributions of the two ANS branchesin the HR. The traditional LF/HF power ratio may thus be a poor measureof autonomic balance.101.3. Computational Eciency is ImportantProblems with HRV analysis extend beyond just the theory, however.Rigid HF band limits at 0.15 and 0.4 Hz are not guaranteed to captureRSA. RSA will only fall within the HF band if the respiration rate (RR) isbetween 9 and 24 breaths/min. If the RR drops below 9 breaths per minute,the RSA power will fall in the LF band, significantly distorting the standardLF/HF power ratio metric. If the RR rises above 24 breaths/min, the HFmeasure of parasympathetic tone will be decreased. Worse, if the RR dropsbelow 9 breaths/min, not only will the HF measure be falsely decreased, butthe LF measure will be falsely increased as well. The parasympathetic tonewill be measured as sympathetic tone.Some work has been done to calculate dynamic HF band limits, basedon the HR and RR [3]. The HR defines a hard upper limit on the HRsampling frequency; frequencies above half the HR cannot be used, as theyfall outside the Nyquist limit. The RR, measured from a di↵erent source,defines the centre frequency of the HF band. The band location and widthare adjusted over time as the RR changes.We can develop a better measure of RSA, using precise a priori knowl-edge of the RR. Most HRV studies use conscious, spontaneously ventilatedsubjects, leading to a highly variable RR. In our studies, subjects are un-conscious and mostly mechanically ventilated, and thus have an RR that isconstant over long periods of time. We can eliminate the concept of the HFband completely, since the RSA power will exist entirely at the known RRwith limited variability in frequency.RSA is often neglected in HRV analyses, but we believe it may hold thekey to quantifying autonomic balance. Respiration is an external perturba-tion to the cardiovascular system and the ANS. By analyzing the perturba-tions caused by respiration, revealed in RSA, we might be able to measurethe autonomic balance as it changes in time.1.3 Computational Eciency is ImportantComputational eciency is an important consideration in a nociceptionmonitor. The algorithms developed and tested in this dissertation wereexecuted on standard Macintosh and Windows PCs, which were runningno other concurrent software and have no significant resource constraints.These devices could easily support computationally complex algorithms,such as the O(n2) Lomb-Scargle periodogram [44, 62] or nonlinear mea-sures such as approximate entropy [55]. However, eventually nociceptionmonitoring would need to be performed on a physiological patient monitor111.3. Computational Eciency is Important(PPM), which has very di↵erent resource constraints.In a PPM, the computational resources will be shared with many otheralgorithms performing physiological measurement in real-time. A PPMmay include pre-processing algorithms for raw 12-lead ECG signals, multi-channel invasive arterial blood pressure signals, capnography, spirometry,and the photoplethysmogram (PPG). It may also include algorithms formeasuring heart rate, blood pressure, blood oxygen saturation, core tem-perature, respiration rate, exhaled carbon dioxide, and minute ventilation,derived from the aforementioned pre-processed signals. It may further in-clude algorithms for detecting dangerous physiological conditions, includingcardiac arrhythmias (e.g. ventricular fibrillation, atrial fibrillation, prema-ture ventricular contractions, asystole, and many more), hemodynamic in-stability, respiratory apnea, and blood oxygen desaturation. If the monitoris an ultra-portable device with highly constrained resources, it may also in-clude algorithms for compressing data. Finally, the monitor would of courseneed to include the algorithms for monitoring nociception. Most algorithmson a PPM must run in so-called “hard” real-time, with strict deadlines fortiming. There are standards and regulatory bodies (e.g. the US Food &Drug Administration or FDA) that dictate the maximum delay allowed fordetecting life-threatening events such as asystole or ventricular fibrillation.Medical device manufacturers must therefore leave significant headroom intheir PPMs’ internal resources, to guarantee that they will always meet thedeadlines under all loading conditions. This means that in a PPM in thereal world, all the algorithms running concurrently in real-time cannot useall the available resources.To further complicate matters, there is a continual push to developsmaller PPMs. Some might argue that riding the curve of Moore’s Lawwill render today’s computationally inecient algorithms tractable in thenot-too-distant future. Moore’s Law is an observation by Intel co-founderGordon Moore in 1965 that the number of transistors on a computer chip,and thus performance, will double approximately every 18-24 months [52].Moore’s Law has held up very consistently until the present day, and willlikely continue for some decades in the future. However, the advances in com-puting performance provide engineers with two diverging options: increasecomputing power within existing form factors, or maintain similar comput-ing power in smaller and/or more energy ecient form factors. Over thelast 20 years, medical device manufacturers have been designing PPMs thatexploit the second option (smaller, lower power). Today’s state-of-the-artPPMs are handheld devices based on embedded system-on-package designs.Such embedded systems have significant resource constraints as compared121.4. Existing Nociception Monitorsto standard desktop PCs. This miniaturization trend will likely continuein the future, with body-worn sensors and PPMs next on the horizon, andeven implantable devices.By designing a nociception monitoring system to require minimal re-sources, we can help ensure it meets the constraints imposed in current andfuture PPM designs.1.4 Existing Nociception MonitorsPreviously published research has described nociception monitors, but nonehas seen widespread acceptance in general anesthesia.The PhysioDoloris Analgesia Nociception Index (ANI) (MetroDoloris,Loos, France) measures changes in the magnitude of RSA as a reflection ofnociception. The algorithm applies a wavelet bandpass filter to the HR toisolate the RSA, then calculates the area under the RSA curve (AUC). TheANI is derived from changes in the short-term and long-term AUC [43]. TheANI assumes fixed frequency bands, as with most traditional time/frequencydomain HRV measures. The wavelet filters assume the RSA will exist in thefrequency range from 0.167 - 0.667 Hz.The Surgical Stress Index (SSI) (GE Healthcare, Helsinki, Finland) mea-sures changes in both the HR and the PPG amplitude as a reflection ofnociception [70]. The parameters are normalized and fused into the SSI,which produces a bounded index from 0 (no nociception) to 100 (strong no-ciception). Changes in HR are considered a reflection of both sympatheticand parasympathetic tones. Changes in PPG amplitude are considered areflection of the heart’s stroke volume and the vascular resistance, both ofwhich are primarily controlled by the sympathetic tone.The MedStorm Stress Detector (MedStorm Innovation, Oslo, Norway)measures the number of fluctuations of skin conductance per second (NFSC)as a reflection of nociception. Sympathetic nerves release acetylcholine,which acts on skin muscarinic receptors to release sweat. When the sweatreaches the skin surface, it increases skin conductivity. MedStorm useschanges in skin conductance as a reflection of sympathetic tone [69]. Earlystudies showed promise for this device, but a later study has shown thatNFSC correlates poorly with post-operative pain measures in children [18].The Brain Anaesthesia Response (BAR) monitor (Cortical DynamicsLtd., North Perth, Australia) measures cortical input as a reflection of no-ciception. The system performs autoregressive modeling of the electroen-cephalogram (EEG) signal, allowing it to estimate the state of the cortex131.5. Taylor’s Experiment & A New Understanding of RSAand the cortical input from the thalamus. The thalamus relays inputs fromlower brain regions (e.g. the brainstem), which themselves relay signalsfrom the body (e.g. nociceptors). Cortical input is then considered a directmeasure of nociception [42].1.5 Taylor’s Experiment & A NewUnderstanding of RSAThere is a strong link between RSA and vagal (parasympathetic) tone. Therelationship between the two phenomena has been shown in spontaneouslybreathing anesthetized dogs [31]. Furthermore, a progressive decrease inRSA power has been shown during progressive cholinergic blockade in hu-mans [26, 59]. Taylor et al. conducted an experiment to further explore therelationship between respiration, RSA, and the sympathetic and parasym-pathetic tones [72].Ten healthy subjects (8 male, 2 female; 20-34 years of age) participated inthe study. Subjects stepped their breathing down from 15 to 3 breaths/min(0.25 - 0.05 Hz) in sequential decrements of 1 breath/min for 96 s at eachstep. A computer display provided visual cues to control subjects’ respira-tion. In one experiment, subjects maintained constant alveolar ventilationby increasing tidal volume as the respiration slowed. In a second experiment,subjects maintained a constant high tidal volume as the respiration slowed.Subjects performed the ramped breathing experiments under three condi-tions in fixed order: saline injection (control), atenolol injection (0.2 mg/kg,-adrenergic blockade), and atropine sulfate injection (0.04 mg/kg, mus-carinic cholinergic blockade). The atenolol was chosen to selectively blockthe e↵ects of the sympathetic nervous system on the heart. The atropine sul-fate was chosen to block the e↵ects of the parasympathetic nervous systemson the heart. Together, they blocked both the sympathetic and parasympa-thetic e↵ects. Subjects lay supine for 10 minutes prior to injections, and thedrugs were given 7 minutes to take e↵ect, after which subjects were tiltedhead-up 40. All measurements were taken with subjects tilted, starting 3minutes after the tilt. Measurements were taken during constant and in-creasing tidal volume, performed in random order. Subjects rested for 10minutes supine between breathing ramps.RSA power increased during -adrenergic blockade, and decreased dur-ing double autonomic blockade, as compared with the control. The RSApower was consistently stronger during -adrenergic blockade at each res-piration frequency. The increase in RSA power was evident during both141.6. Structure of the Dissertationconstant and increasing tidal volume. The RSA power was nearly elimi-nated during double autonomic blockade.RSA power was found to increase with decreasing frequency. This ef-fect was observed during both constant and increasing tidal volume experi-ments. The trend of increasing RSA power held for respiratory frequenciesfrom approximately 0.25 Hz down to 0.08 Hz (15 breaths/min down to 5breaths/min). The trend was reversed at frequencies below this point, withthe RSA power gradually declining from respiratory frequencies of 0.08 Hzdown to 0.05 Hz (5 breaths/min down to 3 breaths/min). The experi-ment did not measure data outside this range, but presumably the trendshold. The same trends of changes in RSA power were also evident during-adrenergic blockade.The experiments revealed that RSA power is driven by parasympathetictone, and restrained by sympathetic tone. Administering a drug that se-lectively blocks the action of the sympathetic nervous system caused anincrease in RSA power. This contrasts with traditional beliefs that sym-pathetic tone has no e↵ect on RSA power. RSA thus reflects the balancebetween the sympathetic and parasympathetic tones, and thus the overallautonomic tone.The experiments also revealed that respiration frequency a↵ects the RSApower. This suggests that a nociception monitor based on the strength ofthe RSA power (such as the PhysioDoloris ANI [43]) may report biasedresults at di↵erent respiration frequencies. Furthermore, the experimentdemonstrated that the RSA can exist at frequencies outside the standardHF band from 0.15 - 0.4 Hz.Taylor’s experiment suggests it should be possible to design a nocicep-tion monitor based on RSA alone, provided it does not simply measure theRSA power. A technique called cardiorespiratory coherence may meet theserequirements.1.6 Structure of the DissertationThis dissertation is organized in five chapters.Chapter 1 describes the problem of nociception monitoring, andpresents some HRV analysis approaches that may be used to solve the prob-lem.Chapter 2 describes how cardiorespiratory coherence can be used tomonitor nociception. It describes the problem of joint time/frequency anal-ysis and presents several alternative methods to solve the problem. It de-151.7. Contributionsscribes a novel real-time filter developed to track the RSA as it moves intime and frequency, so that coherence can be measured between the HR andthe respiration in real-time.Chapter 3 describes a clinical study of nociception and antinocicep-tion during general anesthesia that was conducted to test the nociceptionmonitor. Several di↵erent experiments are described, and their results arepresented. The major findings of the study are detailed and put into clinicalcontext.Chapter 4 describes the architecture of the real-time nociception moni-tor. It describes many di↵erent algorithms that were developed to detect andreject artifacts in the physiological signals, so that cardiorespiratory coher-ence can be measured accurately. The architecture is presented graphicallyin Fig. 4.1.Chapter 5 concludes the dissertation, and presents suggestions for fu-ture work in nociception monitoring.1.7 ContributionsWe have made two novel and significant contributions to the field. The firstis a novel filter for measuring RSA in real-time. The second is a novel end-to-end software system that integrates multiple signals and algorithms forrobust continuous real-time nociception monitoring during general anesthe-sia.We have developed a novel filter for measuring the RSA in the HRVin real-time. The filter uses information from a secondary signal source (arespiration rate derived from a respiration signal) to track the RSA as itmoves in time and frequency. It then uses this information to dynamicallytune the centre frequency and bandwidth of a band-pass filter to isolateand measure the RSA while excluding other components of HRV. A blockdiagram of the filter is shown in Fig. 1.4. The details are described inSection 2.4. Previous work in the field of HRV assumes the RSA will fallwithin a strict frequency range; however, we show in this dissertation thatthe standard assumption is invalid (see Sections 3.8, 3.9, and 3.10). Thisfilter is very e↵ective at tracking the RSA in time and frequency, and it mayprovide the most precise measure of RSA yet devised.We have integrated multiple signals and algorithms together into an end-to-end software system for robust continuous real-time nociception monitor-ing during general anesthesia. The software system incorporates not onlyour dynamically-tuned RSA filter to measure nociception, but also many161.7. ContributionsInstantaneous HRCoherence &SmoothingRespirationRespirationRateAnalyzingFilterAnalyzingFilterFigure 1.4: The real-time CRC algorithm. The RR tunes the analyzingfilters to isolate the RSA while rejecting other HRV components. This dy-namic RSA tracking algorithm is a novel contribution to the field.peripheral algorithms for detecting and rejecting artifacts in the input sig-nals. The input signals required for real-time nociception monitoring can beextremely noisy, and artifacts are a very significant challenge. We have de-veloped methods to address the artifacts and we have deployed them withinthis software system. A block diagram of the system architecture is shownin Fig. 1.5. The details are described in Chapter 4. To our knowledge,no other nociception monitor includes such robust artifact handling usingredundant signal channels.171.7.ContributionsDetectZeroCrossingsSection 4.2.2 ResampleHRSection 4.3.1ResampleRespSection 4.3.2ResampleRRSection 4.3.3 DetectHRArtifactsSection 4.4.1DetectRespArtifactsSection 4.4.2DetectRRArtifactsSection 4.4.3 CorrectHRArtifactsSection 4.5.1CorrectRespArtifactsSection 4.5.2CorrectRRArtifactsSection 4.5.3 RealTimeCRCSection 2.4 NociceptionIndexSection 2.5PhysiologicalMonitor(Raw Data) PPG CO2 or O2Paw or FlowDetectRPeaksSection 4.2.1 ResampleHRSection 4.3.1 DetectHRArtifactsSection 4.4.1 CorrectHRArtifactsSection 4.5.1ECG EstimateSQISection 4.6.2EstimateSQISection 4.6.2 SelectHRSourceSection 4.6.3ResampleRespSection 4.3.2 DetectRespArtifactsSection 4.4.2 CorrectRespArtifactsSection 4.5.2 EstimateSQISection 4.6.2EstimateSQISection 4.6.2 SelectRespSourceSection 4.6.3RRFigure 1.5: The nociception monitor system architecture. Each box represents an algorithm required for real-time monitoring of nociception. Many are peripheral algorithms used for pre-processing of the signals or artifactdetection/rejection. The algorithms are described in the sections referenced. This highly redundant and robustarchitecture is a novel contribution to the field.18Chapter 2Monitoring Nociception withCardiorespiratory Coherence2.1 A Review of Joint Time/Frequency Analysis2.1.1 The Heisenberg-Gabor Uncertainty PrincipleIn 1946, Gabor showed that the Heisenberg uncertainty principle also ap-plied in time/frequency signal analysis [27]. In signal processing, a rectangleof localization over time and frequency can be interpreted as a “quantum ofinformation” in the time/frequency plane [46]. Fig. 2.1 shows an example ofa “Heisenbox” of joint localization in time and frequency. The parameterst and f denote the standard deviation of its spread in time and frequency,respectively.Gabor’s extension of the Heisenberg uncertainty principle, commonlycalled the Heisenberg-Gabor uncertainty principle, places a limit on thejoint time/frequency localization. The limit is defined by [46]:tf  14⇡ . (2.1)Much of the modern work in joint time/frequency analysis involves de-signing analysis kernels or basis functions to trade o↵ localization in onedomain for another. For example, a particular transform might shrink thewidth of the Heisenbox (i.e. smaller t) at the expense of an increasedheight (i.e. larger f ). These basis functions have a distribution in the jointtime/frequency plane, as defined by the shape of their Heisenbox [20].In the next subsections, we will explore a number of di↵erent jointtime/frequency transforms in context of their localization in the two do-mains.2.1.2 The Fourier TransformThe Fourier transform (FT) is the original time/frequency analysis method.While the FT in its abstract sense is technically not a joint time/frequency192.1. A Review of Joint Time/Frequency AnalysisTimeFrequency tfFigure 2.1: Heisenbox representing the spread in time and frequency of aquantum of information in the time/frequency plane.method, it is instructive to review its strengths and weaknesses. The FTconvolves the signal of interest with a series of sinusoids of di↵erent frequen-cies: F{x(t)} = Z +11x(t)e2⇡itkdt, (2.2)where x(t) is the signal being analyzed, the complex exponential is theFourier analysis kernel, and k is a real number defining the analysis fre-quency.The abstract FT provides exquisite localization in the frequency do-main, but zero localization in the time domain. The energy of the complexexponential is spread across all time (zero localization), at a single point infrequency (perfect localization). Fig. 2.2 illustrates the FT’s localization intime and frequency.The abstract FT is generally inappropriate for real-time analysis ofbiomedical signals. Real-time analysis requires good localization in the timedomain. Furthermore, biomedical signals can be highly nonstationary. In-deed, many features of interest are measured as transient phenomena oc-curring at specific points in time. The abstract FT, which spans across alltime, cannot accurately measure such features. Practical implementationsof the FT use windowing techniques to limit the FT in time, and will bediscussed in the next section.202.1. A Review of Joint Time/Frequency AnalysisTimeFrequencyFigure 2.2: Heisenboxes representing the Fourier transform’s localization intime and frequency.Another (related) weakness of the FT is that it assumes the signal isstationary and extends indefinitely for all time. Real signals are limitedin time, so when the FT analyzes a finite signal window it assumes thatthe signal repeats. That is, it assumes that the last point in the windowwraps around and connects to the first point in the window. Typically thesepoints will not align perfectly, leading to a perceived sharp discontinuity inthe signal and thus so-called edge artifacts in the frequency domain repre-sentation. Windowing techniques can be used to mitigate this problem, andwill be discussed in the next subsection.2.1.3 The Short Time Fourier TransformThe short time Fourier transform (STFT) was developed to address theshortcomings of the traditional FT. The STFT segments the time domainsignal into short windows that are narrow enough to be assumed as station-ary, and then computes the FT on them. A windowing function must bechosen to suppress edge artifacts from the analysis. The STFT is definedas: STFT{x(t)}(t, f) = Z +11x(⌧)w(⌧  t)e2⇡i⌧kd⌧, (2.3)212.1. A Review of Joint Time/Frequency AnalysisTimeFrequencyFigure 2.3: Heisenboxes representing the short time Fourier transform’s lo-calization in time and frequency.where x(t) is the signal being analyzed, the complex exponential is theFourier analysis kernel, k is a real number defining the analysis frequency,and w(t) is a windowing function.Typical windowing functions include (but are not limited to) the Han-ning, Gaussian, Hamming, and Blackman [46]. The windowing function ischosen to compromise between time/frequency localization and suppressionof spectral leakage from edge artifacts. The windowing functions cannoteliminate edge artifacts entirely; they can only suppress them.The STFT increases the localization in the time domain, at the expenseof the localization in the frequency domain. Fig. 2.3 illustrates the STFT’slocalization in time and frequency. The localization in the time domain canbe improved by shortening the window, but this will necessarily come at theexpense of decreased localization in the frequency domain. The engineerchooses the window length to suit the needs of the specific signal beinganalyzed.The STFT’s fixed analysis windows cannot provide both good time reso-lution at high frequencies and good frequency resolution at low frequencies.Since the STFT uses analysis windows of a fixed length, its localization isfixed across time and frequency. Notice in Fig. 2.3 that both boxes haveexactly the same size and shape at both points in time and frequency.222.1. A Review of Joint Time/Frequency Analysis2.1.4 The Gabor TransformThe Gabor transform is a special case of the STFT, where the windowingfunction is a Gaussian [47]. Gabor proposed the transform as a wave de-composition in terms of Gaussian wave packets. These wave packets are theresult of a sinusoid multiplied by a Gaussian. The Gaussian windowing func-tion localizes the signals in time. The Gaussian has an interesting propertyin that it is a Gaussian in both the time domain and the frequency domain.As a result, the Gabor transform allows for similarly localized measurementsin both time and frequency.The Gabor transform is given by the following equation:G{x(t)}(t, f) = Z +11x(⌧)e⇡(⌧t)2e2⇡i⌧kd⌧, (2.4)where x(t) is the signal being analyzed, the complex exponential is theFourier analysis kernel, k is a real number defining the analysis frequency,and the real exponential is the Gaussian windowing function.2.1.5 The Wigner-Ville DistributionThe Wigner-Ville Distribution (WVD) is similar to a spectrogram. It de-scribes the signal’s energy distribution in the joint time/frequency domain.The WVD provides greater time and frequency resolution than the STFT[58]. The WVD is computed as:WVD{x(t)}(t, f) = Z +11x⇣t+ ⌧2⌘x⇤ ⇣t ⌧2⌘ e2⇡i⌧kd⌧, (2.5)where x(t) is the signal being analyzed, the asterisk (*) denotes the complexconjugate, the complex exponential is the Fourier analysis kernel, and k isa real number defining the analysis frequency.The improved resolution of the WVD comes at the expense of additionalartifacts. Multiple spectral peaks or harmonics will produce artifacts alongthe frequency axis. Analogously, multiple events in time or multiple phasesarriving from the same event will produce artifacts along the time axis. Thespectral energy computed from the WVD can be negative, which has nophysical meaning. Finally, the most important limitation of the WVD iscross-term interference. Cross-term interference results in a product of thesignal and additive noise together, which makes it dicult to later separatethe noise from the spectrum. Components that exist at di↵erent times232.1. A Review of Joint Time/Frequency AnalysisHhigh(t)x(t)Hlow(t)22Level 1 CoeffsHhigh(t)Hlow(t)22Level 2 CoeffsHhigh(t)Hlow(t)22Level 3 Coeffs. . .Figure 2.4: Filter bank implementation of the DWT.and frequencies can become mixed together, and appear as artifacts in thetime/frequency representation. If there is noise for a brief time in the signal,the WVD will cause that noise to be manifested at other times, and possiblyfor all time in the case of signals of infinite duration [20].These issues with WVD are well known, and limit its application tosolving real problems [58]. Indeed, real-time nociception monitoring requiresthat transient signal features (e.g. a change in RSA) be localized in time.The WVD may cause these transients to appear at multiple times, obscuringthe interpretation of the results.2.1.6 The Discrete Wavelet TransformThe discrete wavelet transform (DWT) divides the frequency domain intodyadic slices [46]. Each slice exists at a separate “scale”, loosely analogousto a frequency band.WT{x(t)}(t, s) = 1|s| Z +11 x(⌧) ✓⌧  ts ◆ d⌧, (2.6)where x(t) is the signal being analyzed,  is the wavelet analysis kernel, ands is the wavelet scale.In the DWT, the scale parameter always takes on a value of a power oftwo (e.g. 1, 2, 4, 8, ...). This power progression decomposes the frequencyaxis into dyads, or regions whose bandwidth is related by a power of two.The DWT is typically implemented as a filter bank, where the signal isrecursively split into high frequency and low frequency halves and down-sampled. The low frequency halves continue to be split in this fashion. Thefilter bank is illustrated in Fig. 2.4.Many unique mother wavelets have been designed for use with the DWT.These include the Daubechies, Haar, biorthogonal, and many, many more.Each wavelet family is designed to meet di↵erent mathematical criteria, andcan be used to detect di↵erent kinds of features in signals. Most of these242.1. A Review of Joint Time/Frequency AnalysisTimeFrequencyFigure 2.5: Heisenboxes representing the discrete wavelet transform’s local-ization in time and frequency.mother wavelets do not have direct links to Fourier frequency. While scalecan be loosely related to frequency, it is often dicult to link the two directly.Fig. 2.5 illustrates the DWT’s localization in time and frequency. No-tice that it has greater time localization at high frequencies than at lowfrequencies, and greater frequency localization at low frequencies than athigh frequencies. Notice also how the DWT decomposes the time/frequencyplane diadically (i.e. in powers of two). Finally, note that the term “fre-quency” is not entirely precise when speaking generally about the DWT,but it is a fair approximation in context of the other transforms describedin this section.The DWT has important limitations. First, its dyadic decompositionof the time/frequency plane can be somewhat restrictive, particularly if wewish to analyze features that span the border between two dyads. Secondly,its “scales” can be dicult to interpret in terms of traditional Fourier fre-quencies typically used in signal analysis.2.1.7 The Continuous Wavelet TransformThe continuous wavelet transform (CWT) is an extension of the DWT thatallows for non-dyadic scales [46]. Its transform equation is the same as thatof the DWT (see Eqn. 2.6), but the scale parameter s can take on anypositive real number value. Many of the discrete mother wavelets cannot be252.1. A Review of Joint Time/Frequency AnalysisTimeFrequencyFigure 2.6: Heisenboxes representing the continuous wavelet transform’slocalization in time and frequency.extended to the continuous case, including the popular Daubechies wavelet.The CWT thus uses a di↵erent set of wavelets, including the Morlet andMexican hat.The CWT uses a scaling function to determine how the wavelet coversthe frequency spectrum. The scaling function is often denoted as (t). TheCWT requires a scaling function because unlike the DWT, its coverage ofthe time/frequency plane is not pre-determined.Fig. 2.6 illustrates the CWT’s localization in time and frequency. No-tice that it has greater time localization at high frequencies than at lowfrequencies, and greater frequency localization at low frequencies than athigh frequencies.2.1.8 The Stockwell TransformThe Stockwell transform (S-transform) is an extension of the CWT. It can bethought of as a CWT with a specific mother wavelet (the Morlet), multipliedby a phase correction factor [67]. The S-transform is also related to theGabor transform, through its use of a Gaussian windowing function.The S-transform maintains a direct link to the Fourier frequency spec-trum. Wavelets analyze a signal at di↵erent scales, which may not havedirect correspondence to frequencies. However, the Morlet basis functionused in the S-transform does have a direct correspondence. The S-transform262.2. What is Coherence?features amplitude and modulation adjustment in phase space, and a linearfrequency scale. The S-transform is given by the following equation:S{x(t)}(t, f) = Z +11x(⌧) |f |p2⇡ e (⌧t)2f22 e2⇡i⌧fd⌧, (2.7)where x(t) is the signal being analyzed, the complex exponential is theFourier analysis kernel, and the real exponential is the Gaussian windowingfunction.The S-transform has been growing in popularity. In one application, itwas used to discriminate between normal and seizure conditions in neonatalEEGs [56].Historically, the S-transform was very computationally intensive. A re-cent advancement has decreased its complexity by several orders of magni-tude [16]. Soon it may become competitive with other transforms that areimplementable with the fast Fourier transform algorithm.2.2 What is Coherence?Coherence is a measure of the coupling between two signals. It estimateshow much of one signal is reflected in the other. Cardiorespiratory coherence(CRC) estimates how much of the HR can be explained by respiration.Examples of coherent and incoherent HR and respiration signals are shownin Figs. 2.7 and 2.8, respectively.Coherence has traditionally been measured using Fourier analysis, butmore recently the wavelet transform has begun to find application in thefield. The continuous wavelet transform (CWT) is a widely used tool for vi-sualizing and analyzing signals across time and frequency. Much of the workon applying wavelets to coherence analysis is based on the work of Torrence& Compo [74]. The wavelet transform (either continuous or discrete) is amathematical technique that converts signals between the time and scale(somewhat analogous to frequency) domains, with varying degrees of local-ization in each domain. With the CWT, coherence can be measured in thejoint time/frequency plane. Recently, coherence has also been measured us-ing the WVD [53]. However, in light of the problems and limitations of theWVD described in the previous section, we will not explore its use further.Initially, a simple formulation of the CWT was developed by Torrence& Compo, with some early discussion of applications in coherence analy-sis [74]. The authors distributed a free software library for performing theCWT. Another group (Grinsted et al.) then extended the software library,272.2. What is Coherence?0 5 10 15 20 25 3076788082848688HR (bpm)0 5 10 15 20 25 30024Respiration (% CO2)Time (s)Figure 2.7: Example of coherent HR and respiration signals from the clinicalstudy of nociception (see Section 3.1). Top plot: HR. Bottom plot: respira-tion. The light grey region denotes an inspiration phase, and the dark greyregion denotes an expiration phase. Note that the sinusoidal pattern fromthe respiration signal is evident in the HR signal. The patient’s respirationis having a strong e↵ect on the HR.282.2. What is Coherence?0 5 10 15 20 25 3095100105HR (bpm)0 5 10 15 20 25 300246Respiration (% CO2)Time (s)Figure 2.8: Example of incoherent HR and respiration signals from theclinical study of nociception (see Section 3.1). Top plot: HR. Bottom plot:respiration. The light grey region denotes an inspiration phase, and the darkgrey region denotes an expiration phase. Note that the sinusoidal patternfrom the respiration signal is not present in the HR signal. The patient’srespiration is having little e↵ect on the HR.292.3. Wavelet Transform Cardiorespiratory Coherenceand used wavelet transform coherence (WTC) to analyze geophysical timeseries for correlation and causation [29]. Keissar et al. then used the WTCsoftware library to analyze cardiorespiratory coherence. They first comparedWTC during exercise and at rest [33]. Later, they performed a more de-tailed mathematical analysis of WTC, and applied it to detecting changes inposture (supine vs. standing) [34]. The work of Keissar et al. demonstratedthat WTC can reflect changes in autonomic tone.We propose that it may be possible to use cardiorespiratory coherencefor monitoring nociception. In our work, we have demonstrated that WTCresponds to nociception and antinociception. This work will be described inthe coming chapters of this dissertation. We began by modifying the soft-ware library from Grinsted et al. to create our own enhanced wavelet-basedalgorithm called wavelet transform cardiorespiratory coherence (WTCRC).We found that WTCRC is moderately correlated with the LF/HF ratio (R= 0.55) [11], that it can detect patient movement during general anesthe-sia (a sign of strong nociception) [15], and that it responds to nociceptionand antinociception during general anesthesia [14]. We then adapted thealgorithm for application in real-time, and showed that it responds to no-ciception and antinociception more strongly than traditional HRV-basedmeasures [13], [12].The next sections describe the WTCRC and real-time CRC algorithms,and the next chapter describes the clinical experiments performed to char-acterize their performance as nociception monitors.2.3 Wavelet Transform CardiorespiratoryCoherenceWTCRC is a measure of CRC using the wavelet transform to analyze thesource signals in the joint time/frequency domain. The WTCRC algorithmdecomposes the HR and the respiration signals using the CWT. At any givenscale, the wavelet transform (in discrete time) is given by:Wn(s) = N1Xn0=0xn0 ⇤ (n0  n)ts  , (2.8)where xn is the input signal, n is the time index, s is the scale, t is thesampling time, and the asterisk (⇤) is the complex conjugate operator. Acomplex Morlet wavelet is used as  , as its scales are directly related toFourier frequencies. The HR and respiration signals are decomposed to 96302.3. Wavelet Transform Cardiorespiratory Coherence0 10 20 30 40 50 60 70 80 90 10000.20.40.60.811.21.41.61.82Wavelet ScaleFourier FrequencyFigure 2.9: CWT decomposition scales used in the WTCRC analysis, andtheir corresponding Fourier frequencies. Only the circled points are ana-lyzed.scales, at 12 scales per octave. The CWT scales and corresponding fre-quencies are shown in Fig. 2.9. These parameters provide good localizationacross the relevant time/frequency spectrum of HRV. The result is a 2D ma-trix of wavelet coecients at di↵erent times and frequencies. The waveletcoecients for the HR and respiration are denoted as W Tn and WRn , respec-tively.From the wavelet coecients, the algorithm calculates the wavelet powerspectrum for each signal, as well as the cross power spectrum:W TTn (s) = W Tn (s)W T⇤n (s),WRRn (s) = WRn (s)WR⇤n (s),W TRn (s) = W Tn (s)WR⇤n (s). (2.9)The HR power, respiration power, and cross power are shown in Figs. 2.10,2.11, and 2.12, respectively. Power spectra are then smoothed in time with a312.3. Wavelet Transform Cardiorespiratory CoherenceGaussian window (en2/2s2) and in scale with a rectangular window (length0.6 x scale).The algorithm uses the smoothed power spectra to calculate the coher-ence estimator: Cˆ2n(s) = ⌦W TRn (s) · s1↵2h|W TTn (s)| · s1i h|WRRn (s)| · s1i , (2.10)where the angled brackets (hi) denote the smoothing operator. The coher-ence estimator is a 2D matrix of coherence values at di↵erent times andfrequencies. It was calculated using crosswavelet and wavelet coherencesoftware provided by A. Grinsted.Finally, the algorithm extracts the coherence values at the known respira-tory frequency at each point in time, using respiratory rate values calculateda priori from the respiration signal. WTCRC outputs a 1D vector of time-varying coherence values at the respiratory frequency, denoted Cˆ2. WTCRCis strictly bounded between 0 (no coherence) to 1 (strong coherence).2.3.1 Limitations of WTCRCWTCRC has some limitations that hinder its usefulness as a nociceptionmonitor.WTCRC is computationally inecient. The underlying CWT operatorsimultaneously analyzes signals across a wide range of times and frequencies.This makes it an ideal tool for visualizing signal power and coherence in thejoint time/frequency plane, and exploring the relationship between HR andrespiration. However, the CWT performs many redundant calculations thatare not required to produce an index of nociception. For example, if thepatient is breathing at 0.25 Hz, we only require the CRC at this frequency.Nevertheless, WTCRC will still analyze the signals across the full spectrumfrom 0 - 2 Hz.WTCRC is imprecise. The CWT analyzes signals at frequencies on afixed grid, as shown in Fig. 2.9. The respiration frequency is unlikely tocoincide perfectly with a point on the grid, so the CWT will rarely capturethe full RSA power in the HR and respiration signals. The RSA poweris typically spread across adjacent scales in the grid. WTCRC smoothsthe power across scales in an attempt to compensate for the inaccuracy.However, this solution is not ideal.WTCRC su↵ers significant real-time delay. The smoothing is overly ag-gressive, in part to compensate for the imprecision described above. WTCRC322.3. Wavelet Transform Cardiorespiratory Coherence50 100 150 200 250 300 350 4008090100110Heart Rate (bpm)Time (s)Frequency (Hz)50 100 150 200 250 300 350 400  1.94  1.29 0.862 0.576 0.384 0.256 0.171 0.1140.07620.0509patient moving  124 s 180 s  bolus dose of propofol & remifentanil 330 s  capnometer rezeroFigure 2.10: Example HR wavelet power spectrum, as an intermediate stepin WTCRC. Data are from the clinical study of nociception (see Section3.1). Top plot: Heart rate. Bottom plot: HR power spectrum, calculatedwith the CWT. Bright areas indicate strong power. Vertical dashed lines(- -) denote clinical events and the dot-dash line (- - .) denotes a rezeroingartifact in the capnogram. The patient movement event is a sign of strongnociception. Note the change in HR frequency content (a shift to a lowerfrequency) during the period of nociception from t = 75 - 175.332.3. Wavelet Transform Cardiorespiratory Coherence50 100 150 200 250 300 350 4000246Respiration (% CO2)Time (s)Frequency (Hz)50 100 150 200 250 300 350 400  1.94  1.29 0.862 0.576 0.384 0.256 0.171 0.1140.07620.0509patient moving  124 s 180 s  bolus dose of propofol & remifentanil 330 s  capnometer rezeroFigure 2.11: Example respiration wavelet power spectrum, as an interme-diate step in WTCRC. Data are from the clinical study of nociception (seeSection 3.1). Top plot: Capnogram. Bottom plot: Capnogram power spec-trum, calculated with the CWT. Bright areas indicate strong power. Verticaldashed lines (- -) denote clinical events and the dot-dash line (- - .) denotesa rezeroing artifact in the capnogram. The patient movement event is a signof strong nociception. Note the disturbances in respiration frequency con-tent during the period of nociception from t = 75 - 175. During this period,the patient is resisting the mechanical ventilator. The dominant respirationfrequency is 0.2 Hz (12 breaths/min), and there are harmonics at higherfrequencies.342.3. Wavelet Transform Cardiorespiratory Coherence50 100 150 200 250 300 350 4008090100110Heart Rate (bpm)50 100 150 200 250 300 350 4000246Respiration (% CO2)Time (s)Frequency (Hz)50 100 150 200 250 300 350 400  1.94  1.29 0.862 0.576 0.384 0.256 0.171 0.1140.07620.0509patient moving  124 s 180 s  bolus dose of propofol & remifentanil 330 s  capnometer rezeroFigure 2.12: Example cross wavelet power spectrum, as an intermediate stepin WTCRC. Data are from the clinical study of nociception (see Section 3.1).Top plot: Heart rate. Middle plot: Capnogram. Bottom plot: Cross powerspectrum, calculated with the CWT. Bright areas indicate strong power.Vertical dashed lines (- -) denote clinical events and the dot-dash line (- - .)denotes a rezeroing artifact in the capnogram. The patient movement eventis a sign of strong nociception. Note the disturbances in frequency contentduring the period of nociception from t = 75 - 175.352.3. Wavelet Transform Cardiorespiratory Coherence50 100 150 200 250 300 350 4008090100110Heart Rate (bpm)50 100 150 200 250 300 350 4000246Respiration (% CO 2)Respiration Rate(breaths/min)50 100 150 200 250 300 350 400 11677.551.734.5  2315.410.36.854.573.0550 100 150 200 250 300 350 400020406080100Time (s)WTCRCNociception Indexpatient moving  124 s 180 s  bolus dose of propofol & remifentanil 330 s  capnometer rezeroa)b)c)d)Figure 2.13: Example anesthetic bolus event analysis with WTCRC. Dataare from the clinical study of nociception (see Section 3.1). Plots from top tobottom: a) Heart rate; b) Capnogram; c) Cardiorespiratory coherence map.Bright areas indicate high coherence. Horizontal continuous line (green) in-dicates the respiratory frequency obtained from b); d) WTCRC nociceptionindex. Vertical dashed lines (- -) denote clinical events and the dot-dash line(- - .) denotes a rezeroing artifact in the capnogram. The patient movementevent is a sign of strong nociception. The nociception index is high in theperiod preceding the bolus dose of anesthetic drugs (baseline period, high-lighted in red), and low following it (response period, highlighted in green).The capnometer rezeroing artifact leads to a transient false increase in thenociception index.362.4. Real-Time Cardiorespiratory Coherenceis somewhat slow to respond to changes in coherence, and brief changes maynot elicit much response at all.We can address these limitations by improving the WTCRC signal ana-lyzing kernels. The improved system is described in the next section.2.4 Real-Time Cardiorespiratory CoherenceThe joint time/frequency analysis methods described in Section 2.1 are use-ful for exploring signals in both time and frequency, but this is not necessaryfor nociception monitoring. These methods decompose a given signal at dif-ferent points in time and frequency, and measure the energy localized at eachpoint in a 2-dimensional plane. The time/frequency plane is divided intodi↵erent tiles, with the tile size and shape determined by the Heisenboxescharacteristic of the analysis method used. This approach was used withthe WTCRC algorithm, and was certainly e↵ective. Overall, this decompo-sition approach is very useful for a general class of problems exploring timeand frequency. However, we have information available that renders thisapproach unnecessary and inecient. Recall that RSA reflects the balancebetween the sympathetic and parasympathetic nervous systems, and thusnociception, as discussed in Section 1.5. It is unnecessary to analyze the in-stantaneous HR series across the entire frequency axis when we only requirethe RSA signal energy, which exists at the respiratory frequency. Using asecondary signal source that estimates the respiratory frequency, we can an-alyze the instantaneous HR only at this specific frequency to measure theRSA and thus the level of nociception.The real-time CRC algorithm is a novel adaptation of the existingWTCRC algorithm. We have developed a new kind of filter for measur-ing the RSA in the HRV in real-time. The filter uses information from asecondary signal source (a respiration rate derived from a respiration sig-nal) to track the RSA as it moves in time and frequency. It then uses thisinformation to dynamically tune the centre frequency and bandwidth of aband-pass filter to isolate and measure the RSA while excluding other com-ponents of HRV. Previous work in the field of HRV assumes the RSA willfall within a strict frequency range; however, we will show later in this dis-sertation that the standard assumption is invalid (see Sections 3.8, 3.9, and3.10). This filter is very e↵ective at tracking and isolating the RSA in timeand frequency. In so doing, it improves the precision of RSA measurementas compared to the previous CWT-based coherence estimator. Moreover,it may be the most robust RSA measurement yet devised. In estimating372.4. Real-Time Cardiorespiratory CoherenceInstantaneous HRCoherence &SmoothingRespirationRespirationRateAnalyzingFilterAnalyzingFilterFigure 2.14: The real-time CRC algorithm. The RR tunes the analyzingfilters to isolate the RSA while rejecting other HRV components.the coherence, we also use a custom-designed smoothing filter to reducereal-time delay.The real-time CRC algorithm is depicted visually in Fig. 2.14. The RRtunes the centre frequency and bandwidth of a bandpass analyzing filter.The analyzing filter isolates the RSA while rejecting other HRV components.The filtered signals are then used for coherence estimation and smoothing.2.4.1 Analyzing FilterReal-time CRC begins by analyzing the HR and respiration signals using anovel filter. The filter is an extension of both the Stockwell transform andthe CWT. Like the S-transform, the CRC analyzing filter maintains a directrelationship with Fourier frequencies. Like the CWT, the filter employs anonlinear frequency bandwidth. As the centre frequency decreases, so toodoes the bandwidth. This has the advantage of suppressing low frequencyHRV components. As the RSA frequency decreases, it may approach thespectral region of low frequency and very low frequency fluctuations in theHR. By tightening the frequency bandwidth in these regions, the analyzingfilter is better able to reject undesired HR fluctutations that are adjacent tothe RSA power.The filter analyzes the signals using a customized complex Morlet basisfunction. This filter is a complex exponential modulated by a Gaussian: (t) = e2⇡ifctefbt2 , (2.11)382.4. Real-Time Cardiorespiratory Coherence−4 −3 −2 −1 0 1 2 3 4−0.2−0.100.10.2AmplitudeTime (s)  Real ψ(t)Imaginary ψ(t)ModulatorFigure 2.15: The complex Morlet analyzing filter in the time domain. Sam-pling frequency fs is 4 Hz, centre frequency fc is 0.25 Hz, and cuto↵ is at2a.where fc is the centre frequency and fb is the freqency domain bandwidth.The bandwidth term is defined as:fb = 2fc/fs, (2.12)where fs is the sampling frequency. As fb increases, so does the filter’sfrequency domain bandwidth. The bandwidth term is complementary tothe variance of the Gaussian modulator:2a = 1/2fb, (2.13)where 2a denotes the variance, and a is the standard deviation of theanalyzing filter. As a increases, so does the filter’s time domain bandwidth.The variance specifies the time domain bandwidth, which is inversely relatedto the frequency domain bandwidth (by the Heisenberg-Gabor uncertaintyprinciple). The analyzing filter’s time domain representation is shown inFig. 2.15.As fc decreases, fb decreases and 2a increases (Eqns. 2.12 and 2.13).This means that as the centre frequency decreases, the filter grows widerin the time domain and narrower in the frequency domain. This tradeo↵392.4. Real-Time Cardiorespiratory Coherence0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.805101520Frequency (Hz)Filter Power  fc = 0.5 Hzfc = 0.25 Hzfc = 0.15 HzFigure 2.16: The complex Morlet analyzing filter in the frequency domainat di↵erent centre freqencies. Cuto↵ is at 2a.is modeled after the principles of wavelet analysis. The analyzing filter’sfrequency domain representation is shown in Fig. 2.16.This custom analyzing filter tracks the RSA as it moves across thetime/frequency plane. The respiratory frequency (fr) is used as the fil-ter’s centre freqency (fc = fr). This is where the RSA power resides in thefrequency domain. An example analysis of an HR and respiration signal isshown in Figs. 2.17 and 2.18.The analyzing filter is applied to the HR and respiration signals to cal-culate their individual spectral powers and cross powers. The filtered HRand respiration signals are denoted as FHRt (fr) and FRespt (fr), respectively.Their individual and cross powers are then given by:PHRt (fr) = FHRt (fr)FHR⇤t (fr),PRespt (fr) = FRespt (fr)FResp⇤t (fr),PXt (fr) = FHRt (fr)FResp⇤t (fr), (2.14)where the superscript X denotes the cross power. The parameter fr andsubscript t denote that the signal powers exist at the respiratory frequencyat di↵erent points in time, both of which are variable. The ⇤ denotes thecomplex conjugate operator. An example of signal and cross powers is shown402.4. Real-Time Cardiorespiratory Coherence50 100 150 200 250 300 350 4008090100110Heart Rate(bpm)50 100 150 200 250 300 350 400−10−50510Heart Rate Filtered(bpm)Time (s)Figure 2.17: Example analysis of the HR signal using the real-time CRCanalyzing filter. Data are from the clinical study of nociception (see Section3.1). Top plot is the raw heart rate. Bottom plot is the filtered heart rate.The solid blue line is the real component, and the dashed green line is theimaginary component. The missing data at the start of the filtered signal iscaused by the width of the analyzing filter.in Fig. 2.19.2.4.2 Smoothing FilterThe spectral powers in Eqn. 2.14 are then smoothed in time using the left(causal) half of a Gaussian filter:St = et2/22s , for t  0, (2.15)where s is the filter’s standard deviation, and defines the level of smoothing.As s increases, the filter’s frequency bandwidth decreases, and the level ofsmoothing increases. As with the analyzing filter, the smoothing filter mustbe truncated. A cuto↵ point of 3s is used, producing an error of approxi-mately 1 %. Since the smoothing filter is causal, it produces no additionaltime delay. The loss of the non-causal (right) half of the filter results inlower frequency domain precision, however, allowing some higher frequencyenergy to bleed through. Fig. 2.20 illustrates the frequency response of thesmoothing filter. This design trades o↵ real-time delay for filter quality.412.4. Real-Time Cardiorespiratory Coherence50 100 150 200 250 300 350 4000246Respiration(% CO2)50 100 150 200 250 300 350 400−505Respiration Filtered(% CO2)Time (s)Figure 2.18: Example analysis of the respiration signal using the real-timeCRC analyzing filter. Data are from the clinical study of nociception (seeSection 3.1). Top plot is the raw CO2 signal. Bottom plot is the filteredCO2. The solid blue line is the real component, and the dashed green lineis the imaginary component. The missing data at the start of the filteredsignal is caused by the width of the analyzing filter.2.4.3 Coherence EstimatorFinally, the algorithm calculates the coherence estimator:Cˆ2t (fr) = ⌦PXt (fr)↵2⌦PHRt (fr)↵ DPRespt (fr)E , (2.16)where the angled brackets denote the smoothing operator.The result is a series of real-time cardiorespiratory coherence values atthe time-varying respiratory frequency. Coherence can range from 0 (nocoherence, strong nociception) to 1 (perfect coherence, no nociception). Anexample of the CRC resulting from the real-time algorithm is shown in Fig.2.24.2.4.4 Real-Time DelayThe real-time CRC analyzing filter can be interpreted in terms of its groupdelay. We can calculate the group delay using Matlab’s grpdelay() function.The group delay for the analyzing filter is shown in Fig. 2.21. In this422.4. Real-Time Cardiorespiratory Coherence50 100 150 200 250 300 350 40050100150200Heart Rate Power(bpm2)50 100 150 200 250 300 350 400020406080Respiration Power(% CO22 )50 100 150 200 250 300 350 400−100−50050100Cross Power(bpm x % CO2)Time (s)Figure 2.19: Example signal powers from the real-time CRC analyzing filterexamples in Figs. 2.17 and 2.18. Top plot is the HR power. Middle plotis the CO2 power. Bottom plot is the cross power. The solid blue line isthe real component, and the dashed green line is the imaginary component.The missing data at the start of the signals is caused by the width of theanalyzing filter.432.4. Real-Time Cardiorespiratory Coherence0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1010203040506070Frequency (Hz)Filter Power  Full filter, σ = 5Causal filter, σ = 5Causal filter, σ = 10Causal filter, σ = 15Figure 2.20: Frequency response of the causal Gaussian smoothing filter atvarious smoothing levels. The full (noncausal) filter response is also shownfor comparison. Note that the causal filter has zero time delay, but is lessecient at rejecting high freqencies.442.4. Real-Time Cardiorespiratory Coherence0 2 4 6 8 10 12 14 16 18−0.2−0.100.10.20.3Sample NumberSample Weight  RealImaginary0 0.2 0.4 0.6 0.8 1 1.288.599.5Frequency (Hz)Group Delay (samples)Figure 2.21: Characteristics of the real-time CRC analyzing filter, with cen-tre frequency of 0.15 Hz and sampling frequency of 2.5 Hz. Top plot: ana-lyzing filter in the time domain. Bottom plot: group delay of the analyzingfilter.example, the analyzing filter cuto↵ is at 1.44a, the centre frequency is at0.15 Hz, and the sampling frequency is at 2.5 Hz. These are the settingswe will later derive in Chapters 3 and 4. The top plot shows the filtercoecients in the time domain. The bottom plot shows the group delay ofthe analyzing filter. The relationship between the filter’s centre frequencyand the group delay is shown in Fig. 2.22.From an alternative perspective, the analyzing filter’s group delay de-pends on the filter length. Longer filters produce longer group delay. Theanalyzing filter’s centre point is sampled, and the remaining points are sam-pled symmetrically away from the centre point at a rate of fs. The filterhas odd length, with the centre tap corresponding to the current point intime. Half of the filter analyzes the signal’s past, and half analyzes its fu-ture. The delay is caused by the right half of the filter, which operates onfuture (non-causal) signal values. The analyzing filter’s group delay is given452.4. Real-Time Cardiorespiratory Coherence0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.50246810121416Respiration Frequency (Hz)Group Delay (samples)Figure 2.22: Group delay of the real-time CRC analyzing filter across arange of respiration frequencies from 0.05 - 0.50 Hz (3 - 30 breaths/min).462.5. Deriving a Nociception Indexby: td = (N  1)/2fs, (2.17)where N is the total number of filter taps.The group delay can be mitigated by truncating the filter in the timedomain. The analyzing filter does not have compact support, because bothof its components (the complex exponential and Gaussian) extend to infinity.Truncating the filter beyond a cuto↵ point limits it to a finite number oftaps (finite N in Eqn. 2.17). This gives it a finite time delay, but introducessome error into the analysis. The cuto↵ point can be tuned to compromisebetween measurement error and real-time delay.The real-time CRC smoothing filter is best interpreted in terms of itsgroup delay. The group delay for the smoothing filter was calculated usingMatlab’s grpdelay() function, and is shown in Fig. 2.23. In this example,the smoothing filter bandwidth is at  = 15 and the cuto↵ is at 3. Theseare the settings we will later use in Chapters 3 and 4. The top plot showsthe filter coecients in the time domain. The bottom plot shows the groupdelay of the smoothing filter. Note that although the filter is entirely causalin the time domain, it still produces some group delay from its smoothingaction. Note also that the smoothing filter is not linear phase; the groupdelay decreases as frequency increases.Real-time delay is a function of the group delays of the CRC analyz-ing and smoothing filters. The analyzing filter’s delay ranges from 5 - 16samples, depending on the patient’s respiration rate. The smoothing filter’sdelay is fixed with regard to respiration rate, but it does vary with frequency.The worst case smoothing filter delay is 29 samples. Therefore, the totaldelay of the real-time CRC filters ranges from 34 to 45 samples, or 13.6 to18 seconds at the chosen sampling frequency of 2.5 Hz.2.5 Deriving a Nociception IndexA nociception index was defined as a dimensionless linear scale ranging from0 (no nociception) to 100 (strong nociception). This definition is intuitive toanesthesiologists, who are already familiar with existing indices of conscious-ness that operate on a similar scale. Consciousness indices are reported bydevices such as the Bispectral Index (BIS) (Aspect Medical Systems Inc.,Newton, MA, USA), Entropy (GE Healthcare, Helsinki, Finland), and Neu-roSense (NeuroWave Systems Inc., Cleveland Heights, OH, USA). The no-ciception index would be a natural complement to these devices.472.5. Deriving a Nociception Index0 20 40 60 80 1000.050.10.15Sample NumberSample Weight0 0.2 0.4 0.6 0.8 1 1.20510152025Frequency (Hz)Group Delay (samples)Figure 2.23: Characteristics of the real-time CRC smoothing filter, with thebandwidth at  = 15, the cuto↵ at 3, and sampling frequency of 2.5 Hz.Top plot: smoothing filter in the time domain. Bottom plot: group delay ofthe smoothing filter.482.6. Conclusion50 100 150 200 250 300 350 40000.20.40.60.81CRC50 100 150 200 250 300 350 400020406080100CRC Nociception IndexTime (s)Figure 2.24: Example CRC and nociception index following from the ex-amples in Figs. 2.17, 2.18, and 2.19. Top plot is the CRC. Bottom plot isthe CRC nociception index. The missing data at the start of the signalsis caused by the combined widths of the analyzing and smoothing filters.Unlike in the previous example figures, this width is not directly related toreal-time delay.Given an algorithm that measures nociception, its nociception index (NI)is derived as the fractional distance between the algorithm’s minimum andmaximum levels of nociception:NI = 100 xCurrentValueMinLevelMaxLevelMinLevel(2.18)CRC ranges from 0 to 1, with a MinLevel (lowest nociception) at 1 andMaxLevel (highest nociception) at 0, so its nociception index is calculatedas:CRC = 100 x (1 Cˆ2), (2.19)where Cˆ2 is the coherence estimator from Eqn. 2.16.An example of the di↵erence between CRC and its nociception index isshown in Fig. 2.24.2.6 ConclusionIn this chapter, we charted the history of joint time/frequency signal anal-ysis and the concept of cardiorespiratory coherence, and we described our492.6. Conclusioncontributions to its evolution. We derived and detailed a novel real-timealgorithm for measuring CRC and proposed its application as a nociceptionindex.Our real-time CRC analyzing filter is a novel and significant contributionto the field. We developed this novel filter for measuring the RSA in theHRV in real-time. The filter uses information from a secondary signal source(a respiration rate derived from a respiration signal) to track the RSA as itmoves in time and frequency. It then uses this information to dynamicallytune the centre frequency and bandwidth of a band-pass filter to isolate andmeasure the RSA while excluding other components of HRV (see Section2.4). Previous work in the field of HRV assumes the RSA will fall withina strict frequency range; however, we show in this dissertation that thestandard assumption is invalid (see Sections 3.8, 3.9, and 3.10). This filteris very e↵ective at tracking the RSA in time and frequency, and it mayprovide the most robust measure of RSA yet devised.We used the new analyzing filter as part of a real-time measure of car-diorespiratory coherence. We propose that real-time CRC can be used formonitoring nociception during general anesthesia. In the next chapter, wewill test real-time CRC in a series of clinical experiments on real patientsreceiving general anesthesia during surgery.50Chapter 3Clinical Validation3.1 Clinical Study of NociceptionThe nociception study timeline is illustrated in Fig. 3.13.1.1 Patients & AnesthesiaThis study was approved by the University of British Columbia ClinicalResearch Ethics Board. Parents of all subjects provided written informedconsent to participate in the study. Subjects were receiving general anes-thesia during restorative dental surgery. All subjects were between 3 - 6years of age, and had ASA physical status I or II. Exclusion criteria in-cluded cardiorespiratory disease, developmental delay, ANS dysfunction,neuromuscular disease, cutaneous disease, chronic pain, a history of headinjury, contraindication to propofol or remifentanil, anticipated dicult air-way management, or use of any medications that may alter ANS function(e.g. anticholinergics, alpha-agonistors, anticonvulsants).Subjects received standard preoperative oral analgesic medications ofacetaminophen (15 mg/kg, up to a maximum of 375 mg) and ibuprofen (10mg/kg, up to a maximum of 250 mg).Propofol and remifentanil were administered for general anesthesia. Anes-thesia was induced with a bolus dose of 1% propofol based on a manuallyadjusted pediatric target controlled infusion (TCI) algorithm using the Paed-fusor pharmacokinetic model [1]. A target plasma level of 4 µg/mL of propo-fol was initially selected. The patient’s trachea was intubated with a nasalendotracheal tube (ET). Following insertion of the ET, the target propofoldrug concentration was reduced to 3.5 µg/mL (maintenance hypnosis) andan infusion of remifentanil was initiated and manually adjusted to achievea target blood concentration of 2 ng/mL (maintenance analgesia). Propo-fol and remifentanil were infused separately using Alaris (CareFusion, SanDiego, CA, USA) automatic infusion pumps. Infusion rates were manuallyadjusted every 5 minutes in accordance with the pharmacokinetic model.The anesthesiologist could deliver rescue bolus doses of propofol513.1. Clinical Study of Nociception3 minInductionIncrease RRDecrease PPawEmergenceDental DamIntubationSNS 13 min 5 min SNS 25 min 5 min SNS 35 min 5 min SNS 4Propofol InfusionRemifentanil InfusionRemiBolus RR 8 brpmFigure 3.1: The nociception study timeline.(1 mg/kg) and/or remifentanil (0.5 µg/kg) at their discretion. These res-cue boluses were administered to suppress strong conscious or autonomicarousal in the patient. The anesthesiologist could also deliver bolus doses ofmorphine and/or fentanyl towards the end of the case to assist with post-operative pain management.3.1.2 Pre-Induction Setup & Data CollectionPhysiological data were recorded and annotated throughout each case, us-ing three separate monitors. MedStorm Stress Detector software, runningon a standard Windows notebook PC, recorded skin conductance from threepalmar electrodes applied to the patient’s right hand. The MedStorm soft-ware was used to create annotations in real-time during the procedures. ADatex/Ohmeda S/5 physiological monitor (GE Healthcare, Helsinki, Fin-land) was used to record a 5-lead (3-channel) electrocardiogram (ECG) at300 Hz, capnometry (CO2, O2) and spirometry (flow, airway pressure) at25 Hz, PPG at 100 Hz, and noninvasive blood pressure (NIBP) at 1/180Hz (the NIBP cu↵ was cycled every three minutes). The S/5 monitor cal-culated and recorded 10 s averaged trends from the recorded signals, andupsampled all waves to match the maximum sampling frequency of 300 Hz.A NeuroSense monitor (NeuroWave Systems Inc., Cleveland Heights, OH,USA), consisting of two data collection devices and specialized software run-ning on a touch screen Windows PC, recorded two EEG channels from fourelectrodes applied to the patient’s forehead and a 3-lead (1-channel) ECG,both at 900 Hz.The three monitors were synchronized for post hoc data analysis. Near523.1. Clinical Study of Nociceptionthe beginning of each case, after data recording was initiated on all threedevices and before the patient entered the operating room (OR), two setsof markers were created to provide synchronization points. The process re-quired two research assistants: one at the MedStorm monitor and the otherat the Datex monitor. The MedStorm and Datex monitors were synchro-nized first. The research assistants would synchronize themselves verbally(e.g. “1, 2, 3, Press!”), then each press a button to create a new markerat the same time. The MedStorm marker was uniquely labeled to indicateits purpose (e.g. “Sync with Datex”), and the Datex marker number wasnoted (usually it was the first marker). The time di↵erence between thetwo markers corresponds to the system clock o↵set between the MedStormand Datex monitors. The Datex and NeuroSense monitors were synchro-nized next. A research assistant would synchronize him/herself with theNeuroSense monitor by watching its clock seconds counting, then choose anupcoming synchronization time (e.g. the next 10 s rollover). When the cho-sen time arrived, the research assistant would press a button on the Datexmonitor to create a new marker. The research assistant would read out theNeuroSense time to the research assistant at the MedStorm monitor, whowould then create a new MedStorm marker. The time of this MedStormmarker was not important. The marker was uniquely labeled to indicateits purpose (e.g. “Sync with NeuroSense”), along with the synchronizationtime read from the NeuroSense monitor. The Datex marker number wasnoted (usually it was the second marker). The time di↵erence between theDatex marker and the recorded NeuroSense time corresponds to the systemclock o↵set between the Datex and NeuroSense monitors. With these twoo↵set times, the three monitors can be synchronized in time, and their dataanalyed together (Fig. 3.2).A research assistant annotated the data in real-time with markers iden-tifying significant clinical events, including dental dam insertions, surgicalstimulation, rescue bolus doses of anesthetics, and changes in ventilation.Annotation times were accurate to within 10 s of the recorded time, limitedby the Datex synchronization resolution described above. Some events wereannotated very precisely, approaching this accuracy limit. For example thestart/end of the standardized nociceptive stimulus (SNS) (see Section 3.1.4)were very clear, precise, and unambiguous events. Such precision was im-possible with other events, when their occurrence itself did not exist at asingle point in time. For example, bolus doses of anesthetics were markedat the end of the injection, though the injection typically began severalseconds earlier. Other events, such as cavity drillings, often spanned manyminutes but could only be annotated periodically to provide an approximate533.1. Clinical Study of NociceptionSynchronization Point 1Synchronization Point 2MedStorm DataDatex DataNeuroSense DataTimeCase Begins Case EndsFigure 3.2: Two synchronization points are used to synchronize data fromthree monitors in the nociception study. Timeline is not drawn to scale.understanding of the current level of surgical stimulation. Precision of theannotations thus depended on the type of annotation.3.1.3 Induction of Anesthesia & Sensor AttachmentAfter induction of anesthesia, recording devices were attached to the patient.Three RAs and a nurse were typically required for this process. One researchassistant would apply the three MedStorm palmar electrodes, typically onthe right hand. Next, the research assistant would apply the two tetanicnerve stimulator electrodes, typically on the left forearm. The MedStormand nerve stimulator were always applied to opposite sides of the patient,to reduce the e↵ect of electromagnetic interference. At the same time, asecond research assistant would apply the ECG electrodes. Five electrodeswere attached for the Datex/Ohmeda monitor, at the left arm (LA), leftleg (LL), right arm (RA), right leg (RL), and chest (C) positions. Threeadditional electrodes were attached for the NeuroSense high frequency ECGmonitor, at the LA, LL, and RL positions. At the same time, a third researchassistant would apply the NeuroSense EEG electrodes. Four EEG electrodeswere attached for the NeuroSense monitor, on the forehead. The skin wasabraded with sandpaper, cleaned with an alcohol swab, dried with a drycotton swab, and finally electrodes were applied. At the same time, a nursewould attach a PPG clip, typically to the left index finger. The nurse wouldalso attach a noninvasive blood pressure cu↵, typically to the right upperarm.Following electrode attachment, a standardized nociceptive stimulus(SNS) was administered (see Section 3.6). This SNS occurred when the543.1. Clinical Study of Nociceptionpatient was receiving propofol anesthesia only, and theoretically providedthe most reliable nociceptive stimulus. Unfortunately it also occurred priorto intubation, when the patient was receiving manual bag mask ventilation,and CRC measurement was not possible during this period. CRC performsbest when the respiration is stationary, and the bag mask ventilation ishighly nonstationary.Following a bu↵er period after the SNS, a bolus dose of remifentanil wasadministered. The trachea was then intubated with a nasal ET. Completionof intubation marked the start of the maintenance phase of anesthesia.3.1.4 Maintenance Phase of AnesthesiaIn the maintenance phase of anesthesia, the patient is anesthetized andmechanically ventilated. Ventilation is consistent and stationary for longperiods, which is important for reliable CRC measurement.Multiple periods of nociception and antinociception occurred during themaintenance phase of anesthesia. A dental dam was inserted at least oncein each case, providing a strong nociceptive stimulus. An SNS was admin-istered up to three times during the maintenance phase of each case (seeSection 3.6). Other nociceptive stimuli included tooth extractions, cavitydrillings, and cap/crown insertions. Antinociception interventions, in theform of rescue bolus doses of anesthetics, often followed nociceptive periods.The study protocol included three ventilatory maneuvers to test di↵erentconditions on RSA.In the first ventilatory maneuver, called V1, the RR was set to8 breaths/min. V1 was set shortly after completion of tracheal intuba-tion, shortly after the start of the maintenance phase of anesthesia. Aftera five minute stabilization period, a second SNS was administered for 60 s.Another five minute bu↵er period was provided after the end of the SNS.In the second ventilatory maneuver, called V2, the RR was increased to16 breaths/min. V2 was set after completion of V1. After a five minutebu↵er period, a third SNS was administered for 60 s. Another five minutebu↵er period was provided after the end of the SNS.In the third ventilatory maneuver, called V3, the PPaw was decreasedto 12 cmH2O. V3 was set after completion of V2. After a five minute bu↵erperiod, a fourth SNS was administered for 60 s. Another five minute bu↵erperiod was provided after the end of the SNS.The protocol allowed the anesthesiologist to override the ventilatory ma-neuvers if the patient’s end-tidal CO2 (EtCO2) deviated from an acceptablerange of 35 - 45 mmHg. To minimize the impact of overrides on the study,553.2. CRC Analyzing Filter Tuningthe research assistant would ask the anesthesiologist to adjust the PPawduring V1 and V2, and to adjust the RR during V3.3.1.5 EmergenceWhen the dentist had completed surgery, the infusions of propofol andremifentanil were discontinued to begin emergence. Upon spontaneous respi-ration, the anesthesiologist would extubate the trachea of the patient. Afterextubation, all electrodes and sensors were removed and the patient wasmoved to a recovery room. A research assistant would stop data collectionon the MedStorm, Datex, and NeuroSense monitors, and copy the data filesto a portable disk for post hoc analysis.3.2 CRC Analyzing Filter TuningThe CRC analyzing filter length was tuned to compromise between real-time delay and filter quality. A cuto↵ point of 5a was used as the goldstandard, since cuto↵ at this point produces nearly zero error. Cuto↵ valuesranging from 1a to 5a were tested, in increments of 0.01a (401 values ofa tested in total). CRC was calculated for each time point in the dentaldam and anesthetic bolus data segments described below (Sections 3.5 and3.7), using each of these di↵erent cuto↵ values. CRC was averaged in 60s windows and compared to the gold standard. The mean squared error(MSE) was calculated across all windows for each cuto↵ value.An analyzing filter quality of 1.44a was chosen as a good tradeo↵ be-tween CRC estimation accuracy and real-time delay (Fig. 3.3). This ana-lyzing filter produces a theoretical real-time delay ranging from 2.2 - 3.4 s inthe standard HF band (0.15 - 0.4 Hz, or 9 - 24 breaths/min). All subsequentCRC analyses were performed using this value.3.3 Nociception Index Tuning3.3.1 DefinitionA nociception index was defined as a dimensionless linear scale boundedfrom 0 (no nociception) to 100 (strong nociception). This definition is in-tuitive to anesthesiologists, who are already familiar with existing indicesof consciousness that operate on a similar scale. Consciousness indices arereported by devices such as the BIS (Aspect Medical Systems Inc., Newton,MA, USA), Entropy (GE Healthcare, Helsinki, Finland), and NeuroSense563.3. Nociception Index Tuning1 1.5 2 2.5 3 3.5 4 4.5 5020406080100120Mean Squared ErrorAnalyzing Filter Quality (SD)Figure 3.3: The overall MSE of CRC results, using a range of analyzingfilter quality values. The error is calculated relative to an analyzing filterwith quality of 5 a.(NeuroWave Systems Inc., Cleveland Heights, OH, USA). The nociceptionindex would be a natural complement to these devices.Given an algorithm that measures nociception, its nociception index isderived as the fractional distance between the algorithm’s minimum andmaximum levels of nociception:NI = 100 xCurrentValueMinLevelMaxLevelMinLevel(3.1)Some measures are bounded and others are unbounded. The nociceptionindex will be derived di↵erently in the next sections depending on theirboundedness.3.3.2 Derivation from a Bounded MeasureCRC ranges from 0 to 1, with a MinLevel (lowest nociception) at 1 andMaxLevel (highest nociception) at 0, so its nociception index is calculatedas:CRC = 100 x (1 Cˆ2). (3.2)573.3. Nociception Index TuningHFnu and LFnu (the normalized powers in the HF and LF bands, re-spectively) similarly range from 0 to 1, but have opposite MinLevel andMaxLevel. HFnu should decrease in response to nociception (as with CRC),while LFnu should increase.3.3.3 Derivation from an Unbounded MeasureThe remaining traditional HRV measures are not strictly bounded, so limitswere estimated empirically. The LF power, HF power, LF/HF ratio, SDNN,and RMSSD all have a lower bound of zero, but have no upper bound. Eventhe lower bound may never be reached in practice. The algorithm resultswere calculated for each time point across all cases in the dataset. Any timepoints a↵ected by HR artifact (missed or false detected beats) were omitted.Any time points with respiration frequency outside the standard HRV HFband (0.15 - 0.4 Hz) were excluded from consideration in the frequencydomain measures (LF, HF, LFnu, HFnu, and LF/HF ratio), since this isknown to corrupt the results [11]. The largest value observed for each of thealgorithms was used as its upper bound, and the smallest as its lower bound.The HF power should decrease in response to nociception, so the upperbound was used as the minimum level of nociception. All other measuresshould increase in response to nociception, so their upper bound was usedas the maximum level.The LF, HF, and LF/HF ratio were measured on a logarithmic scale.Prior informal testing has revealed that these measures have a very highdynamic range, and produce a strongly skewed nociception index. A loga-rithmic transformation produces a better nociception index. The 60 s slidingwindow was selected as a compromise between frequency localization andreal-time delay. Shorter windows impose shorter delays, but provide lessability to discriminate between frequency bands.Nociception indices for the unbounded traditional HRV measures werederived empirically across all 48 nociception study cases (see Section 3.1).The frequency domain measures were run against 32 hours of clean data.The time domain measures were run against 52 hours of clean data. The[minimum, maximum] observed values were [85.06, 170.1] log(ms2) for LFpower, [83.18, 166.2] log(ms2) for HF power, [-49.01, 45.08] for the LF/HFratio, [0.000, 157.3] ms for SDNN, and [0.000, 108.1] ms for RMSSD. All sub-sequent traditional HRV analyses were performed using nociception indicesderived from these bounds.583.4. The Baseline-Stimulus-Response ExperimentTT - TBB T + TRBT - TBB - TBBaseline PeriodResponse PeriodStimulus EventSegment Start Segment EndT + TRB + TRBaseline Buffer PeriodResponse Buffer PeriodEdge ArtifactBuffer PeriodTBB TRB TRTBT + TRB + TR + TEBT - TBB - TB - TEBEdge ArtifactBuffer PeriodTEB TEBFigure 3.4: Data segment structure for a general BSR experiment.3.4 The Baseline-Stimulus-Response ExperimentThe Baseline-Stimulus-Response (BSR) experiment is designed to estimatethe change in a nociception index in response to a given stimulus. A no-ciceptive stimulus, for example, should elicit an increase in a nociceptionindex, and an antinociceptive stimulus should elicit a decrease. The BSRexperiment measures the nociception index value before a known stimulus(in a baseline period) and after it (in a response period). The di↵erencebetween the two values is attributed to the stimulus.3.4.1 Data SegmentationA data segment is analyzed around each stimulus event (see Fig. 3.4).Bu↵er periods may be applied before and/or after the stimulus, to ensure thestimulus has time to a↵ect the ANS or to ensure the baseline and responseperiods don’t cross contaminate. The widths of the bu↵ers can be adjusteddepending on the nature of the stimulus. The widths of the baseline andresponse periods can likewise be adjusted. Bu↵er periods are also applied atthe start and end of each data segment, to ensure the baseline and responseperiods are not corrupted by edge artifacts.3.4.2 Statistical AnalysisThe response to any one event is unpredictable, but the mean response toa large number of events should estimate the true response. Like flippinga slightly weighted coin, a single trial does not provide much information,593.5. Response to Dental Dambut as the number of trials grows larger the coin’s bias should be revealed.There are many external factors that may influence the outcome of anyone trial. It is known that stimuli unrelated to the experimental processare being applied to the patient at various times throughout each case.These independent stimuli may occur at the same time as the experimentalstimulus, and may influence the results. The BSR experiment is conducted“blind” to such independent stimuli, however. Consider, for example, theSNS response experiment described in Section 3.6. The SNS is a nociceptivestimulus, so it is expected that a nociception index should increase from thebaseline period to the response period. However, if the dentist happensto extract a tooth during the baseline period of one event, this event mayactually show a decrease in the nociception index. The BSR experimenttreats such external stimuli as independent noise with zero mean. That is,the external stimuli may occur during either the baseline or response periodswith equal probability. Then, as the number of events grows large, thereshould be just as many unexpected stimuli in both the baseline and responseperiods, and their net e↵ect should cancel out. Individual trials may showan unexpected result, but the overall mean change in response to a stimulusconverges on the true value.The mean change in response to a stimulus is measured. The algorithmis used to analyze each data segment, and its results are averaged over thebaseline and response periods. The change from the baseline to the responseperiod is calculated. Changes are averaged over all events, to arrive atan overall mean change. The percent change is calculated from the meanbaseline and change results.A 95% confidence interval is estimated about the mean change, using acorrected percentile bootstrapping method (the bootci function in Matlab).A Wilcoxon rank-sum test is applied to estimate the statistical significanceof the mean change. A response below the 5% level is considered statisticallysignificant. A Bonferroni correction is applied to account for the multiplesignificance tests being performed. Nine tests will be performed in eachof the experiments that follow, leading to a corrected significance level of0.05/9 = 0.0056.3.5 Response to Dental Dam3.5.1 Experimental ConceptA dental dam inserted into the mouth is a strongly nociceptive stimulus.Dental dams are inserted prior to dental surgery in anesthetized children to603.5. Response to Dental DamTT - 60 s T + 60 sT - 120 sT - 240 s T + 180 sBaseline PeriodResponse PeriodDental Dam EventSegment Start Segment EndFigure 3.5: Data segment structure for a dental dam event.stop debris from entering the airway and to keep the teeth dry. The jaw iswidely distracted with a retractor and a rubber dam is clamped to the teethand gums with steel clips. The dental dam insertion is expected to followa period of low nociception (such as performing dental X-rays), and to leadto a period of increased nociception. Nociception indices should increaseduring the dental dam insertion compared to baseline before insertion.Case annotations were searched to find the first dental dam event in eachcase. The first dental dam always occurred during the maintenance phase ofanesthesia, and the patient was always mechanically ventilated. Events wereexcluded if the ECG or respiration signals contained significant artifacts.Data segments were extracted around each dental dam event (see Fig.3.5). Each segment comprised a baseline period and a response period. Theresponse period consisted of a 60 s analysis window immediately followingthe start of the dental dam event marker time. The baseline period consistedof a 60 s analysis window ending 60 s before the start of the dental dam event.This allowed for a 60 s bu↵er between the periods, to ensure the analysis wasnot corrupted by cross contamination. The analysis windows were paddedwith 120 s of data on each end, to ensure the analysis was not corrupted byedge artifacts.3.5.2 ResultsForty-two dental dam events were analyzed. Five were excluded due tosignificant ECG artifacts, where the HR was impossible to discern, and one613.6. Response to SNSMeasure Mean Mean 95 % CI P-value % "Baseline ChangeCRC 36 27 20 to 32 < 0.000005* 90WTCRC 17 14 7.4 to 19 < 0.0002 * 79LF 58 4.8 1.7 to 7.6 < 0.05 62HF 56 0.97 -2.5 to 3.5 > 0.4 64LFnu 74 8.3 4.3 to 13 < 0.03 83HFnu 81 6.4 3.6 to 10 < 0.02 83LF/HF 67 5.2 2.8 to 7.4 < 0.03 83SDNN 16 5.4 2.6 to 7.9 < 0.002* 74RMSSD 17 -1.1 -3.8 to 2.0 > 0.4 45Table 3.1: Response to a dental dam insertion. Asterisks (*) denote statis-tically significant results, factoring in the Bonferroni correction.was excluded due to a lost CO2 waveform.The CRC nociception index inceased significantly in response to dentaldam insertion. Real-time CRC was more sensitive to the stimulus thanwere the traditional HRV measures, and it was even more sensitive thanthe o✏ine WTCRC algorithm. Dental dam insertion changed the CRCnociception index by an average of 27 [95% CI from 20 to 32], WTCRC by14 [7.4 to 19], LF by 4.8 [1.7 to 7.6], HF by 0.97 [-2.5 to 3.5], LFnu by 8.3[4.3 to 13], HFnu by 6.4 [3.6 to 10], LF/HF by 5.2 [2.8 to 7.4], SDNN by 5.4[2.6 to 7.9], and RMSSD by -1.1 [-3.8 to 2.0] (Table 3.1, Figs. 3.6, 3.7, 3.8,3.9).3.6 Response to SNS3.6.1 Experimental ConceptAn SNS is a mildly nociceptive stimulus. A percutaneous tetanic ulnarnerve stimulator (DualStim Plus, Life-Tech Inc., Houseton, Texas, USA)delivers an electrical current of 50 mA with a frequency of 100 Hz throughthe patient’s ulnar nerve for 60 s. This current intensity and duration isroutinely used to monitor depth of neuromuscular blockade (NMB) duringgeneral anesthesia when NMB agents are used. The stimulus induces tetaniccontractions of the muscles in the forearm, producing strong discomfort andmild nociception. This stimulus has been shown to provide a reproducibleexperimental pain model for surgical pain during propofol/remifentanil anes-thesia [60]. The nociceptive stimulus is standardized, theoretically produc-623.6. Response to SNS50 100 150 200 250 300 350 40080859095100Heart Rate(bpm)50 100 150 200 250 300 350 400012345Respiration(% CO2)50 100 150 200 250 300 350 400020406080100CRCNociception IndexTime (s)start dental dam  240 s end dental dam  347 sFigure 3.6: Example dental dam event, analyzed with real-time CRC. Thisanalysis corresponds to event number 39 in Fig. 3.7. Vertical lines denoteclinical events. The green (left) and red (right) vertical bands represent thebaseline and response periods, respectively. The nociception index is low inthe period preceding the dental dam, and high during it. The missing CRCdata at the beginning of the window corresponds to the combined length ofthe analyzing and smoothing filters.633.6. Response to SNS5 10 15 20 25 30 35 400102030405060708090100Dental Dam Event NumberReal−Time CRC ResponseFigure 3.7: CRC nociception index response for each dental dam event,sorted by response strength. Arrows indicate the direction of the responsefrom the baseline period to the response period. Upward arrows (blue)indicate a change in the expected direction (increase), and downward arrows(red) indicate a change in the opposite direction (decrease).643.6.ResponsetoSNS0102030405060708090100Real−Time CRC0102030405060708090100WTCRC0102030405060708090100LF0102030405060708090100HF0102030405060708090100LFnu0102030405060708090100HFnu0102030405060708090100LF/HF0102030405060708090100SDNN0102030405060708090100RMSSDFigure 3.8: Boxplots of responses to dental dam insertion for all measures. The central red bar representsthe median (second quartile). The box edges are the first and third quartiles. The whiskers extend 1.5x theinterquartile range (IQR) beyond the box edges. Points beyond the whiskers are drawn as red plus symbols (+),and may be considered outliers. Such possible outliers were not excluded from analysis.653.6. Response to SNS0510152025Response to Dental Dam (Index Change)RMSSDSDNNLF/HF HFnu LFnu HF LFWTCRCReal−Time CRC−1.1+5.4 +5.2+6.4+8.3+0.97+4.8+14+27Figure 3.9: Comparison of mean nociception index responses to a dentaldam insertion. Longer bars are better.663.6. Response to SNSTT - 60 s T + 60 sT - 120 sT - 240 s T + 180 sBaseline PeriodResponse PeriodStart SNS EventSegment Start Segment EndFigure 3.10: Data segment structure for an SNS event.ing the same level and duration of nociception each time administered ineach patient. The SNS was administered three times during the maintenancephase of anesthesia in each case. The SNS follows a period of unknown no-ciception (since it is administered during surgery), and is expected to leadto a period of increased nociception.Case annotations were searched to find all three SNS events during themaintenance phase of anesthesia in each case. The patient was always me-chanically ventilated during these events. Events were excluded if the ECGor respiration signals contained significant artifacts.Data segments were extracted around each SNS event (see Fig. 3.10).Each segment comprised a baseline period and a response period. The re-sponse period consisted of a 60 s analysis window immediately following thestart of the SNS marker time. The baseline period consisted of a 60 s analy-sis window ending 60 s before the start of the SNS event. This allowed for a60 s bu↵er between the periods, to ensure the analysis was not corrupted bycross contamination. The analysis windows were padded with 120 s of dataon each end, to ensure the analysis was not contaminated by edge artifacts.3.6.2 ResultsA total of 116 SNS events were analyzed. Eleven were excluded due tosignificant ECG artifacts, where the HR was impossible to discern, two wereexcluded because they occurred in conjunction with a change in respirationfrequency (causing unreliable CRC measurements), two were excluded due to673.7. Response to Bolus Dose of AnestheticMeasure Mean Mean 95 % CI P-value % "Baseline ChangeCRC 34 9.8 6.2 to 13 < 0.0005* 71WTCRC 14 5.2 2.5 to 8.0 < 0.008 68LF 55 2.3 0.26 to 4.6 > 0.1 55HF 49 2.3 0.62 to 4.3 > 0.2 63LFnu 57 7.2 3.9 to 11 < 0.03 60HFnu 67 6.4 3.8 to 9.1 < 0.02 66LF/HF 58 4.2 2.3 to 5.9 < 0.03 62SDNN 18 2.3 0.27 to 4.1 < 0.01 69RMSSD 22 -1.5 -4.1 to 1.2 > 0.2 41Table 3.2: Response to an SNS event. Asterisks (*) denote statisticallysignificant results, factoring in the Bonferroni correction.a lost CO2 waveform, eight were missing because the research assistant wasunable to administer them (e.g. the case ended early), three were missingdue to a failure of the data recording equipment, and two were missing dueto a failure of the tetanic nerve stimulator.The CRC nociception index inceased significantly in response to the SNS.Real-time CRC was more sensitive to the stimulus than were the traditionalHRV measures, and it was even more sensitive than the o✏ine WTCRCalgorithm. The SNS changed the CRC nociception index by an average of9.8 [95% CI from 6.2 to 13], WTCRC by 5.2 [2.5 to 8.0], LF by 2.3 [0.26 to4.6], HF by 2.3 [0.62 to 4.3], LFnu by 7.2 [3.9 to 11], HFnu by 6.4 [3.8 to9.1], LF/HF by 4.2 [2.3 to 5.9], SDNN by 2.3 [0.27 to 4.1], and RMSSD by-1.5 [-4.1 to 1.2] (Table 3.2, Figs. 3.11, 3.12, 3.13, 3.14).3.7 Response to Bolus Dose of Anesthetic3.7.1 Experimental ConceptA bolus dose of anesthetic drug (e.g. propofol, remifentanil, fentanyl, mor-phine) is a strongly antinociceptive stimulus. Bolus doses of anesthetics aretypically given to rescue the patient from a strong sympathetic response tonociception. As such, rescue boluses are expected to follow a period of strongnociception, and to lead to a period of decreased nociception. Nociceptionindices should decrease after a bolus dose of anesthetic drugs compared tobaseline before the bolus.683.7. Response to Bolus Dose of Anesthetic50 100 150 200 250 300 350 4007075808590Heart Rate(bpm)50 100 150 200 250 300 350 400012345Respiration(% CO2)50 100 150 200 250 300 350 400020406080100CRCNociception IndexTime (s)start SNS  240 s 301 s  end SNSdrilling  107 sFigure 3.11: Example SNS event, analyzed with real-time CRC. This anal-ysis corresponds to event number 107 in Fig. 3.12. Vertical lines denoteclinical events. The green (left) and red (right) vertical bands represent thebaseline and response periods, respectively. The nociception index is low inthe period preceding the SNS, and high during it. The missing CRC dataat the beginning of the window corresponds to the combined length of theanalyzing and smoothing filters.693.7. Response to Bolus Dose of Anesthetic10 20 30 40 50 60 70 80 90 100 1100102030405060708090100SNS Event NumberReal−Time CRC ResponseFigure 3.12: CRC nociception index response for each SNS event, sortedby response strength. Arrows indicate the direction of the response fromthe baseline period to the response period. Upward arrows (blue) indicatea change in the expected direction (increase), and downward arrows (red)indicate a change in the opposite direction (decrease).703.7.ResponsetoBolusDoseofAnesthetic0102030405060708090100Real−Time CRC0102030405060708090100WTCRC0102030405060708090100LF0102030405060708090100HF0102030405060708090100LFnu0102030405060708090100HFnu0102030405060708090100LF/HF0102030405060708090100SDNN0102030405060708090100RMSSDFigure 3.13: Boxplots of responses to the SNS for all measures. The central red bar represents the median (secondquartile). The box edges are the first and third quartiles. The whiskers extend 1.5x the interquartile range (IQR)beyond the box edges. Points beyond the whiskers are drawn as red plus symbols (+), and may be consideredoutliers. Such possible outliers were not excluded from analysis.713.7. Response to Bolus Dose of Anesthetic012345678910Response to SNS (Index Change)RMSSDSDNNLF/HF HFnu LFnu HF LFWTCRCReal−Time CRC−1.5+2+4.2+6.5+7.3+2.5 +2.2+5.4+9.7Figure 3.14: Comparison of mean nociception index responses to the SNS.Longer bars are better.723.7. Response to Bolus Dose of AnestheticTT - 60 s T + 60 s T + 120 sT - 180 s T + 240 sBaseline PeriodResponse PeriodAnesthetic Bolus EventSegment Start Segment EndFigure 3.15: Data segment structure for an anesthetic bolus event.Case annotations were searched to find all anesthetic bolus events. Bo-lus doses of propofol, remifentanil, fentanyl, and morphine were consideredto have antinociceptive properties. Multiple bolus doses of anesthetics de-livered in quick succession (<120 s apart) were considered to constitute asingle bolus event. Events were only retained for analysis if they occurredduring the maintenance phase of anesthesia, when the patient was mechan-ically ventilated, and when the ECG and respiration signals were free ofsignificant artifacts.Data segments were extracted around each anesthetic bolus event (seeFig. 3.15). Each segment comprised a baseline period and a response period.The baseline period consisted of a 60 s analysis window immediately pre-ceding the anesthetic bolus event marker time. The bolus dose of anestheticdrug was given a bu↵er period of 60 s to take e↵ect. The response periodconsisted of a 60 s analysis window immediately following the bu↵er period.The analysis windows were padded with 120 s of data on each end, to ensurethe analysis was not corrupted by edge artifacts.3.7.2 ResultsFifty-seven anesthetic bolus events were analyzed. One was excluded dueto significant ECG artifacts, where the HR was impossible to discern, threewere excluded due to significant CO2 artifacts, two were excluded due toa baseline or response period that extended outside the maintenance phaseof anesthesia, two were excluded because they occurred in conjunction with733.8. Response to Change in Respiration RateMeasure Mean Mean 95 % CI P-value % #Baseline ChangeCRC 57 -19 -26 to -12 < 0.00003* 70WTCRC 31 -15 -21 to -9.3 < 0.00002* 75LF 54 -6.5 -10 to -3.0 < 0.02 63HF 60 -1.5 -4.5 to 2.0 > 0.2 53LFnu 73 -12 -19 to -6.9 < 0.002* 75HFnu 80 -9.8 -15 to -5.2 < 0.002* 72LF/HF 68 -7.2 -11 to -4.0 < 0.002* 72SDNN 18 -4.2 -6.5 to -2.0 < 0.005* 68RMSSD 11 0.49 -1.8 to 3.3 > 0.3 61Table 3.3: Response to a bolus dose of anesthetic drugs. Asterisks (*) denotestatistically significant results, factoring in the Bonferroni correction.a change in respiration frequency (causing unreliable CRC measurements),and two were missing due to a failure of the data recording equipment.The CRC nociception index decreased significantly in response to a bolusdose of anesthetic drugs. Real-time CRC was more sensitive to the stimulusthan were the traditional HRV measures, and it may have been even moresensitive than the o✏ine WTCRC algorithm. A bolus dose of anestheticschanged CRC by an average of -19 [95% CI from -26 to -12], WTCRC by -15[-21 to -9.3], LF by -6.5 [-10 to -3.0], HF by -1.5 [-4.5 to 2.0], LFnu by -12[-19 to -6.9], HFnu by -9.8 [-15 to -5.2], LF/HF by -7.2 [-11 to -4.0], SDNNby -4.2 [-6.5 to -2.0], and RMSSD by 0.49 [-1.8 to 3.3] (Table 3.3, Figs. 3.16,3.17, 3.18, 3.19).3.8 Response to Change in Respiration Rate3.8.1 Experimental ConceptA change in RR is a stimulus that is neither nociceptive nor antinociceptive.The anesthesiologist will often change the RR to compensate for an EtCO2that has exceeded acceptable limits. Increasing the RR will decrease theEtCO2, and vice versa. When the RR changes, so too does the location of theRSA in the frequency domain. Nociception indices should be independentof the RR, remaining the same from before to after a change in RR.Case annotations were searched to find all RR change events. Each casecontained a change in RR from 8 to 16 breaths/min, as specified in thestudy protocol. The marked event times were not always exact. The time of743.8. Response to Change in Respiration Rate50 100 150 200 250 300 350 4007580859095100105110Heart Rate(bpm)50 100 150 200 250 300 350 4000123456Respiration(% CO2)50 100 150 200 250 300 350 400020406080100CRCNociception IndexTime (s)patient moving  124 s 180 s  bolus dose of propofol & remifentanil 330 s  capnometer rezeroFigure 3.16: Example anesthetic bolus event, analyzed with real-time CRC.This analysis corresponds to event number 60 in Fig. 3.17. Vertical dashedlines (- -) denote clinical events and the dot-dash line (- .) denotes a reze-roing artifact in the capnogram. The patient movement event is a sign ofstrong nociception. The nociception index is high in the period precedingthe bolus dose of anesthetic drugs, and low following it. The capnometerrezeroing artifact leads to a transient false increase in the nociception index.The missing CRC data at the beginning of the window corresponds to thecombined length of the analyzing and smoothing filters.753.8. Response to Change in Respiration Rate5 10 15 20 25 30 35 40 45 50 550102030405060708090100Anesthetic Bolus Event NumberReal−Time CRC ResponseFigure 3.17: CRC nociception index response for each anesthetic bolus event,sorted by response strength. Arrows indicate the direction of the responsefrom the baseline period to the response period. Downward arrows (blue)indicate a change in the expected direction (decrease), and upward arrows(red) indicate a change in the opposite direction (increase).763.8.ResponsetoChangeinRespirationRate0102030405060708090100Real−Time CRC0102030405060708090100WTCRC0102030405060708090100LF0102030405060708090100HF0102030405060708090100LFnu0102030405060708090100HFnu0102030405060708090100LF/HF0102030405060708090100SDNN0102030405060708090100RMSSDFigure 3.18: Boxplots of responses to a bolus dose of anesthetic drugs for all measures. The central red barrepresents the median (second quartile). The box edges are the first and third quartiles. The whiskers extend 1.5xthe interquartile range (IQR) beyond the box edges. Points beyond the whiskers are drawn as red plus symbols(+), and may be considered outliers. Such possible outliers were not excluded from analysis.773.8. Response to Change in Respiration Rate02468101214161820Response to Anesthetic Bolus (Index Change)RMSSDSDNNLF/HF HFnu LFnu HF LFWTCRCReal−Time CRC+0.49−4.2−7.2−9.8−12−1.5−6.5−15−19Figure 3.19: Comparison of mean nociception index responses to a bolusdose of anesthetic drugs. Longer bars are better.783.8. Response to Change in Respiration RateTT - 60 s T + 60 sT - 120 sT - 240 s T + 240 sBaseline PeriodResponse PeriodChange RR EventSegment Start Segment EndT + 120 sFigure 3.20: Data segment structure for a change in respiration rate event.change in RR was refined post hoc by manually inspecting the CO2 signal,and finding the first breath at the new RR. The refined time is accurateto within 1/2 of a respiration cycle, or approximately 2 s. The change inRR always occurred during the maintenance phase of anesthesia, when thepatient was mechanically ventilated. Events were only retained for analysisif the ECG and respiration signals were free of significant artifacts.Data segments were extracted around each change in RR event (see Fig.3.20). Each segment comprised a baseline period and a response period. Thebaseline period consisted of a 60 s analysis window ending 60 s before thechange in RR event. The response period consisted of a 60 s analysis windowstarting 60 s after the change in RR event. This allowed for a 60 s bu↵erbefore and after the change in RR, to ensure the analysis was not corruptedby cross contamination. The analysis windows were padded with 120 s ofdata on each end, to ensure the analysis was not corrupted by edge artifacts.3.8.2 ResultsA total of 43 change in RR events were analyzed. One was excluded due tosignificant ECG artifacts, where the HR was impossible to discern, two wereexcluded due to significant artifacts in the CO2 signal, one was missing dueto a failure of the data recording equipment, and one was missing becausethe case ended early.The CRC nociception index did not change significantly in response to achange in RR. A change in RR from 8 to 16 breaths/min changed CRC byan average of -2.2 [95% CI from -10 to 4.7], WTCRC by 0.28 [-5.2 to 5.3],LF by -9.2 [-13 to -5.6], HF by -8.2 [-11 to -5.5], LFnu by -29 [-37 to -22],793.9. Response to Change in Peak Airway PressureMeasure Mean Mean 95 % CI P-value % " / #Baseline ChangeCRC 34 -2.2 -10 to 4.7 > 0.3 51 / 49WTCRC 14 0.28 -5.2 to 5.3 > 0.1 65 / 35LF 59 -9.2 -13 to -5.6 < 0.008 26 / 74HF 55 -8.2 -11 to -5.5 > 0.05 16 / 81LFnu 74 -29 -37 to -22 < 0.000004* 19 / 81HFnu 82 -26 -33 to -19 < 0.000002* 16 / 84LF/HF 67 -16 -20 to -12 < 0.000002* 16 / 84SDNN 18 1.4 -0.74 to 3.8 > 0.4 53 / 47RMSSD 19 5.1 1.8 to 8.7 > 0.3 72 / 26Table 3.4: Response to a change in RR from 8 to 16 breaths/min. Aster-isks (*) denote statistically significant results, factoring in the Bonferronicorrection.HFnu by -26 [-33 to -19], LF/HF by -16 [-20 to -12], SDNN by 1.4 [-0.74 to3.8], and RMSSD by 5.1 [1.8 to 8.7] (Table 3.4, Figs. 3.21, 3.22, 3.23, 3.24).3.9 Response to Change in Peak Airway Pressure3.9.1 Experimental ConceptA change in PPaw is a stimulus that is neither nociceptive nor antinoci-ceptive. The anesthesiologist will often change the PPaw to compensatefor an EtCO2 that has exceeded acceptable limits. Increasing the PPawwill decrease the EtCO2, and decreasing the PPaw will increase the EtCO2.When the PPaw changes, so too does the magnitude of the RSA power, evenwhen the ANS state remains unchanged. Nociception indices should be in-dependent of the PPaw, remaining the same from before to after a changein PPaw.Case annotations were searched to find all PPaw change events. Eachcase contained a change in PPaw from 15 to 12 cmH2O, as specified in thestudy protocol. The marked event times were not always exactly recorded.The time of change in PPaw was refined post hoc by manually inspectingthe Paw signal, and finding the first breath at the new PPaw. The refinedtime is accurate to within 1/2 of a respiration cycle, or approximately 2s. The change in PPaw always occurred during the maintenance phaseof anesthesia, when the patient was mechanically ventilated. Events wereonly retained for analysis if the ECG and respiration signals were free of803.9. Response to Change in Peak Airway Pressure50 100 150 200 250 300 350 400 4506065707580Heart Rate(bpm)50 100 150 200 250 300 350 400 450012345Respiration(% CO2)50 100 150 200 250 300 350 400 450020406080100CRCNociception IndexTime (s)240 s  change RR from 8 to 16 breaths/minFigure 3.21: Example change in RR event, analyzed with real-time CRC.This analysis corresponds to event number 19 in Fig. 3.22. Vertical linesdenote clinical events. The yellow vertical bands (both left and right) repre-sent the baseline and response periods, respectively. The nociception indexis approximately the same in the period preceding the change in RR as itis following it. The transient increase in the nociception index immediatelyfollowing the change in RR is caused by a delay in the RR estimation. Forthis brief period, the real-time CRC algorithm is looking for the RSA at thewrong frequency. The missing CRC data at the beginning of the windowcorresponds to the combined length of the analyzing and smoothing filters.813.9. Response to Change in Peak Airway Pressure5 10 15 20 25 30 35 400102030405060708090100Change in RR Event NumberReal−Time CRC ResponseFigure 3.22: CRC nociception index response for each change in RR event,sorted by response strength. Arrows indicate the direction of the responsefrom the baseline period to the response period.significant artifacts.Data segments were extracted around each change in PPaw event (seeFig. 3.25). Each segment comprised a baseline period and a response period.The baseline period consisted of a 60 s analysis window ending 60 s beforethe change in PPaw event. The response period consisted of a 60 s analysiswindow starting 60 s after the change in PPaw event. This allowed for a 60s bu↵er before and after the change in PPaw, to ensure the analysis was notcorrupted by cross contamination. The analysis windows were padded with120 s of data on each end, to ensure the analysis was not corrupted by edgeartifacts.3.9.2 ResultsA total of 35 change in PPaw events were analyzed. Two were excludeddue to significant ECG artifacts, where the HR was impossible to discern,one was missing due to a failure of the data recording equipment, four weremissing because either the start or end PPaw was unachievable, and six weremissing because the case ended early.The CRC nociception index did not change significantly in response to a823.9.ResponsetoChangeinPeakAirwayPressure0102030405060708090100Real−Time CRC0102030405060708090100WTCRC0102030405060708090100LF0102030405060708090100HF0102030405060708090100LFnu0102030405060708090100HFnu0102030405060708090100LF/HF0102030405060708090100SDNN0102030405060708090100RMSSDFigure 3.23: Boxplots of responses to a change in RR for all measures. The central red bar represents the median(second quartile). The box edges are the first and third quartiles. The whiskers extend 1.5x the interquartilerange (IQR) beyond the box edges. Points beyond the whiskers are drawn as red plus symbols (+), and may beconsidered outliers. Such possible outliers were not excluded from analysis.833.9. Response to Change in Peak Airway Pressure051015202530Response to Change in RR (Index Change)RMSSDSDNNLF/HF HFnu LFnu HF LFWTCRCReal−Time CRC+5.1+1.4−16−26−29−8.2 −9.2+0.28−2.2Figure 3.24: Comparison of mean nociception index responses to a changein RR. Shorter bars are better.TT - 60 s T + 60 sT - 120 sT - 240 s T + 240 sBaseline PeriodResponse PeriodChange PPaw EventSegment Start Segment EndT + 120 sFigure 3.25: Data segment structure for a change in peak airway pressureevent.843.10. DiscussionMeasure Mean Mean 95 % CI P-value % " / #Baseline ChangeCRC 32 5.4 -1.0 to 11 > 0.1 74 / 26WTCRC 16 4.1 -0.95 to 9.1 > 0.05 63 / 37LF 59 -2.8 -7.7 to 1.4 > 0.3 43 / 57HF 37 5.3 2.8 to 8.4 > 0.1 71 / 26LFnu 43 4.1 -3.6 to 12 > 0.2 54 / 46HFnu 54 5.2 -1.2 to 12 > 0.1 60 / 40LF/HF 51 2.6 -1.5 to 6.4 > 0.1 60 / 40SDNN 27 -5.9 -11 to -2.1 > 0.2 43 / 54RMSSD 33 -7.2 -14 to -2.7 > 0.2 34 / 60Table 3.5: Response to a change in PPaw from 15 to 12 cmH2O. Asterisks (*)denote statistically significant results, factoring in the Bonferroni correction.change in PPaw. A change in PPaw from 15 to 12 cmH2O changed CRC byan average of 5.4 [95% CI from -1.0 to 11], WTCRC by 4.1 [-0.95 to 9.1], LFby -2.8 [-7.7 to 1.4], HF by 5.3 [2.8 to 8.4], LFnu by 4.1 [-3.6 to 12], HFnuby 5.2 [-1.2 to 12], LF/HF by 2.6 [-1.5 to 6.4], SDNN by -5.9 [-11 to -2.1],and RMSSD by -7.2 [-14 to -2.7] (Table 3.5, Figs. 3.26, 3.27, 3.28, 3.29).3.10 Discussion3.10.1 Major FindingsThe results of these experiments demonstrate that a nociception index basedon real-time CRC is more sensitive to nociception and antinociception thanare traditional time and frequency domain HRV measures. Furthermore,the real-time CRC algorithm is even more sensitive than its o✏ine forebear,WTCRC. The real-time CRC filters impose only a very short delay, allowingit to respond quickly to changes in autonomic state.The CRC nociception index increased significantly in response to noci-ception. Dental dam insertions provided a strong nociceptive stimulus, andthe SNSs provided a smaller nociceptive stimulus. Real-time CRC exhib-ited a statistically significant increase to both stimuli. WTCRC and mostof the traditional HRV measures also increased significantly, though theirincreases were of a smaller magnitude. We can thus conclude that CRC ismore sensitive to nociception than are the alternative indices.The CRC nociception index decreased significantly in response to antinoci-ception. Bolus doses of anesthetic drugs provided a strong antinociceptive853.10. Discussion50 100 150 200 250 300 350 400 45060657075Heart Rate(bpm)50 100 150 200 250 300 350 400 45001234Respiration(% CO2)50 100 150 200 250 300 350 400 450020406080100CRCNociception IndexTime (s)240 s  change PPaw from 15 to 12 cmH2OFigure 3.26: Example change in PPaw event, analyzed with real-time CRC.This analysis corresponds to event number 16 in Fig. 3.27. Vertical linesdenote clinical events. The yellow vertical bands (both left and right) repre-sent the baseline and response periods, respectively. The nociception indexis approximately the same in the period preceding the change in PPaw asit is following it. The missing CRC data at the beginning of the windowcorresponds to the combined length of the analyzing and smoothing filters.863.10. Discussion5 10 15 20 25 30 350102030405060708090100Change in PPaw Event NumberReal−Time CRC ResponseFigure 3.27: CRC nociception index response for each change in PPaw event,sorted by response strength. Arrows indicate the direction of the responsefrom the baseline period to the response period.stimulus. Real-time CRC exhibited a statistically significant decrease tothis stimulus. WTCRC and most of the traditional HRV measures also de-creased significantly, though their decreases were of a smaller magnitude.We can thus conclude that CRC is more sensitive to antinociception thanare the alternative indices.The CRC nociception index did not change significantly in response tonon-nociceptive stimuli. Changes in RR and in PPaw provided strong stim-uli that were neither nociceptive nor antinociceptive. Real-time CRC wasapproximately the same before and after the stimuli. WTCRC and the tra-ditional statistical HRV measures were also unchanged. Most of the tradi-tional frequency domain HRV measures did change significantly in responseto a change in RR, however.Reliance on fixed frequency band limits caused the traditional frequencydomain HRV measures to perform poorly as indices of nociception. Manyevents in this study occurred during periods of low respiration frequency,including 36 of the dental dam events, 39 of the SNS events, and 16 of theanesthetic bolus events. In such events, the RSA is not found in the expectedHF band, but is instead mostly in the LF band. This condition can leadto a muted response, or even a response in the direction opposite of that873.10.Discussion0102030405060708090100Real−Time CRC0102030405060708090100WTCRC0102030405060708090100LF0102030405060708090100HF0102030405060708090100LFnu0102030405060708090100HFnu0102030405060708090100LF/HF0102030405060708090100SDNN0102030405060708090100RMSSDFigure 3.28: Boxplots of responses to a change in PPaw for all measures. The central red bar represents the median(second quartile). The box edges are the first and third quartiles. The whiskers extend 1.5x the interquartile range(IQR) beyond the box edges. Points beyond the whiskers are drawn as red plus symbols (+), and may be consideredoutliers. Such possible outliers were not excluded from analysis.883.10. Discussion01234567Response to Change in PPaw (Index Change)RMSSDSDNNLF/HF HFnu LFnu HF LFWTCRCReal−Time CRC−7.2−5.9+2.6+5.2+4.1+5.3−2.8+4.1+5.4Figure 3.29: Comparison of mean nociception index responses to a changein PPaw. Shorter bars are better.893.10. Discussionexpected. The LF power, HF power, LFnu, HFnu, and LF/HF ratio were alla↵ected by this shortcoming. The existing HRV-based nociception monitors(ANI and SSI, see Section 1.4) are also likely a↵ected by this shortcoming,though they were not tested directly.CRC is unique in combining information from the HR and respirationsignals, operating as a form of sensor fusion. CRC makes no assumptionon the location of the RSA, completely eliminating the concept of fixedfrequency bands. The algorithm tracks the RSA as it moves in time andfrequency using the respiration frequency calculated from the respiration sig-nal (e.g. CO2). CRC has information on the exact location of the RSA. Wehave previously shown that CRC continues to function when the respirationfrequency is low (< 0.15 Hz), while the LF/HF ratio fails [11]. CRC obtainsfurther benefits from analyzing the shape of the respiration signal. Duringperiods of strong nociception the patient may struggle against the mechani-cal ventilator, leading to irregularities in the shape and frequency content ofthe respiration waveform. This feature is evident in Fig. 3.16 in the periodfrom approximately t = 80 - 170 s. CRC measures these changes, leading toimproved sensitivity in detecting nociception/antinociception. Indeed, wehave previously shown that CRC is sensitive to patient movement events,which are a sign of strong nociception [15]. Sensor fusion between HR andrespiration has a drawback, however. Unlike other methods, CRC requiresa clean respiration signal. Dependence on multiple signals increases the riskof artifacts a↵ecting the measurements. A good example of a common arti-fact in the CO2 waveform is shown in Fig. 3.16 at around t = 330 s. Thecapnometer periodically rezeros itself to ensure accurate measurements, andthis feature can lead to false detections of nociception. Appropriate artifacthandling will be essential for clinical implementation.3.10.2 Clinical RelevanceOur experiments were performed only on healthy pediatric subjects in anarrow age range. Children have a more active RSA, and the RSA responseto changing autonomic balance has been shown to decrease with age [36].Further studies will be required to confirm these results in adult subjects,and in subjects with diseases a↵ecting autonomic function (e.g. diabetesmellitus) or medications (e.g. atropine). We have recently shown that CRCresponds to circadian changes in autonomic balance in nine healthy adultsubjects aged 26.3 ± 4.6 years [8]. The response of CRC to these autonomicchanges provides evidence that it may be used to monitor nociception inadult subjects.903.10. Discussion3.10.3 LimitationsThe scaling of the nociception indices used in these experiments overesti-mates the sensitivity of most traditional HRV measures. The lower andupper bounds of the algorithms were set empirically based on the smallestand largest values observed in the dataset. The bounds correspond to thestrongest possible responses, at nociception indices of 0 and 100. It is likelythat smaller and larger values exist for the algorithms, but such conditionswere not observed in our dataset. Wider bounds would compress the algo-rithm’s nociception index scale, decreasing its apparent sensitivity. The LFpower, HF power, and LF/HF ratio are particularly problematic, as theyhave a very large dynamic range. Since CRC naturally produces a boundedscale, no empirical derivation was necessary and the measured sensitivity islikely more accurate.The experiments underestimate the true sensitivity of all measures. Theyassume that insertion of a dental dam or administration of an SNS willalways increase the patient’s level of nociception, and that a bolus dose ofanesthetic drugs will always decrease it. However, these assumptions maynot always hold. For example, in some cases the surgical stimulation canbe stronger before the dental dam than during its insertion. This appearsto be the situation in dental dam event number 1 (see Fig. 3.7). In thisevent, nearly all measures reported a decrease in nociception. Likewise, abolus dose of anesthetic may be insucient to counteract increasingly strongnociceptive stimuli. The patient may experience stronger nociception afterthe anesthetic bolus event than before it. This appears to be the situationin anesthetic bolus event numbers 1 and 2 (see Fig. 3.17). In these events,nearly all measures reported an increase in nociception. These unexpectedresponses are likely to be correct reflections of the change in autonomic state,but since they act in the wrong direction they serve to decrease the reportedoverall mean sensitivity.While our results show that CRC is not a↵ected by a change in RR, theexperiment was only performed on the range of 8 - 16 breaths/min. It ispossible that there is a small e↵ect that may only become measurable atlarger changes in RR. Furthermore, we have shown that CRC can operateacross a wider range than the traditional frequency domain HRV measures,but this range is not infinite. The respiration and HR can only couple acrossa finite frequency range. Indeed, no two oscillators can couple at any arbi-trary frequency. If the patient breathes too slowly or two quickly, the HRwill be unable to couple to it and coherence will drop to zero (CRC nocicep-tion index of 100). The exact lower and upper RR bounds are unknown, and913.11. Conclusionprobably vary across the population. Age, health, and individual physiol-ogy all likely contribute to variation in RSA coupling bounds. These boundscould be tested in the future, by ramping the RR under otherwise steadystate conditions and observing the resulting change in CRC.CRC may have actually responded to a change in PPaw, though theresponse was too small to be statistically significant. Although the meanchange was small (5.4), the majority of events showed an increase (26/35 or74%) in the nociception index. There is a possible theoretical mechanismto explain this e↵ect. The airway pressure in the lungs provides an exter-nal perturbation to the ANS. The perturbation acts both directly (throughstretch receptors in the lungs) and indirectly (through changes in intratho-racic pressure and thus baroreflex). The cyclic perturbation causes corre-sponding oscillations in HR (i.e. RSA, see Section 1.5). When the PPawis decreased, the perturbation e↵ect is likewise decreased. In the extreme,when the PPaw drops to zero, there is no perturbation, no RSA, and noCRC. Decreasing the PPaw may actually decrease CRC, and thus increasethe nociception index. We may test this hypothesis in the future, by de-creasing the PPaw under otherwise steady state conditions and observingthe resulting change in CRC.Our theoretical real-time delay underestimates the true delay. We havecalculated the response delay of the CRC filters to changes in HR and res-piration, but there is a further delay in the response of these physiologicalprocesses as well. When the ANS changes state (e.g. in response to noci-ception), it changes levels of neurotransmitters that a↵ect the functioningof the heart. These neurotransmitters (primarily norepinephrine and acetyl-choline) can take up to 10 s to elicit a change in HR. The total system delayfrom state change (e.g. nociception) to measured response includes both thephysiological and CRC filter delays.Our experiments could not measure the specificity of the nociceptionindices. We have used dental dam insertions, SNSs, and bolus doses of anes-thetics as indicators of true positive responses to nociception and antinoci-ception. No corresponding indicator of true negative responses exists, how-ever, so measuring specificity is beyond the scope of this study.3.11 ConclusionIn this chapter, we described the clinical validation of real-time CRC as anociception index. We described a clinical study of nociception during gen-eral anesthesia, as well as several experiments for testing its performance in923.11. Conclusioncomparison to traditional clinical and HRV-based measures of nociception.CRC shows promise as a real-time monitor of nociception during generalanesthesia. We have shown that CRC responds to both nociception andantinociception, and that CRC is more sensitive than traditional HRV-basedmeasures. In the future, CRC could provide anesthesiologists with feedbackabout the adequacy of analgesia in real time, increasing patient safety duringsurgery.93Chapter 4A Real-Time NociceptionMonitor4.1 Real-Time System ArchitectureA practical nociception monitor requires far more than just an algorithm formeasuring nociception. The analyses performed in Chapter 3 were executedin pseudo-real-time, on manually pre-processed data. The instantaneous HRseries were derived only semi-automatically; that is, the derivation startedwith automatic processing and the results were then cleaned by manual in-spection. False or missed beats in the ECG were manually corrected sothat gold standard instantaneous HR series were used for analysis. Therespiration signals were visually inspected for artifact (e.g. zeroing in thecapnometer). Data were excluded from analysis where respiration or HRartifacts could not be perfectly corrected. This is acceptable practice whenassessing algorithm performance, but it is not acceptable for a real-timesystem. A fully functioning nociception monitor requires many peripheralalgorithms for signal pre-processing. It requires algorithms for detectingbeats and calculating the instantaneous HR series (see Section 4.2), resam-pling the signals onto a common grid (see Section 4.3), detecting artifacts inthe signals (see Section 4.4), correcting any detected artifacts (see Section4.5), and selecting the cleanest signals for analysis (see Section 4.6). Onlyafter all of these pre-processing algorithms have run can we actually measurenociception using the real-time CRC algorithm. Creating and testing thesealgorithms is a significant e↵ort. This chapter describes these peripheralalgorithms.The real-time nociception monitor was built using the iAssist softwareframework [9]. iAssist is designed to be extensible, flexible, and scalable.Small software modules (called plugins) can be written and connected tothe framework. Plugins can be interconnected in a directed acyclic graphstructure, where the outputs of one plugin feed into the inputs of others.Each plugin executes as a separate thread, and the threads run concurrently.944.2. Real-Time Heart Rate MeasurementThe nociception monitor system architecture is illustrated in Fig. 4.1.Each algorithm was developed as a plugin in iAssist, and connected in thegraph structure shown. The system input is nearly raw data from a physio-logical monitor, and the system output is a real-time measure of nociception.Each of the algorithms used in the real-time nociception monitor aredescribed in the sections below. The individual algorithms may be updatedand changed over time as improved methods are developed.4.2 Real-Time Heart Rate Measurement4.2.1 Real-Time HR from the ECGA custom designed R-peak detection algorithm was developed for operationin real-time. It runs in O(n) (linear) time and is simple programmatically.The algorithm is called SimpleDiff.The SimpleDi↵ algorithm detects ECG R-peaks by calculating a simpledi↵erentiation (hence the name). This is equivalent to performing a firstlevel wavelet decomposition with a Haar wavelet. The R-wave is a highfrequency component of the ECG. No other normal ECG feature approachesthis frequency. The algorithm iterates over each ECG sample, and calculatesthe di↵erence from the previous sample and its square:ECGDi↵ = ECG( n ) - ECG( n - 1 ), (4.1)ECGDi↵Squared = ECGDi↵2, (4.2)where n is the current index. Squaring the di↵erence makes all values posi-tive, and increases the contrast between small and large changes. By keepingonly high frequency components of the ECG, the algorithm retains the R-waves while rejecting other unnecessary features such as T-waves, P-waves,and baseline wander (see Fig. 4.2).When the squared di↵erence exceeds a threshold, calledECGDiffThreshold, the algorithm detects that it is in the region of an R-peak. In this region, there may be multiple points that exceed the threshold,as the R-wave climbs sharply, then falls back to baseline. There may alsobe points in this region that do not exceed the threshold. The algorithmthus allows for threshold gaps in the region of an R-peak, with maximumwidth specified by a parameter called MaxGapWidth. When the maximumgap width is exceeded, the algorithm detects that it has left the region ofan R-peak. It retains the start and end indices of the region.The algorithm chooses a peak in the given region as the R-peak. Thealgorithm first detects candidate R-peaks as slope changes in the ECG.954.2.Real-TimeHeartRateMeasurementDetectZeroCrossingsSection 4.2.2 ResampleHRSection 4.3.1ResampleRespSection 4.3.2ResampleRRSection 4.3.3 DetectHRArtifactsSection 4.4.1DetectRespArtifactsSection 4.4.2DetectRRArtifactsSection 4.4.3 CorrectHRArtifactsSection 4.5.1CorrectRespArtifactsSection 4.5.2CorrectRRArtifactsSection 4.5.3 RealTimeCRCSection 2.4 NociceptionIndexSection 2.5PhysiologicalMonitor(Raw Data) PPG CO2 or O2Paw or FlowDetectRPeaksSection 4.2.1ResampleHRSection 4.3.1DetectHRArtifactsSection 4.4.1CorrectHRArtifactsSection 4.5.1ECGEstimateSQISection 4.6.2EstimateSQISection 4.6.2SelectHRSourceSection 4.6.3ResampleRespSection 4.3.2DetectRespArtifactsSection 4.4.2CorrectRespArtifactsSection 4.5.2EstimateSQISection 4.6.2EstimateSQISection 4.6.2SelectRespSourceSection 4.6.3RRFigure 4.1: The nociception monitor system architecture. Each algorithm is represented as a box, and is imple-mented as a plugin inside the iAssist software framework. The algorithms are described in the sections referenced.964.2. Real-Time Heart Rate Measurement1 2 3 4 5−50510Amplitude (ms)ECG1 2 3 4 5−2−1012Difference (ms)ECG Difference1 2 3 4 50123456Squared Difference (ms2 )ECG Squared DifferenceTime (s)Figure 4.2: Example analysis with the SimpleDi↵ algorithm. The top plotshows a somewhat noisy ECG, with significant baseline wander. The middleplot shows how the noise has been rejected, leaving only the relevant highfrequency content. The bottom plot shows the squared di↵erence and R-peak detection threshold (the horizontal red dotted line).974.2. Real-Time Heart Rate MeasurementWhen the ECGDiff changes sign within an R-peak region, the index is flaggedas a candidate R-peak. Small peaks are then rejected. The ECG amplitudeat the candidate peak is measured, as is the minimum ECG amplitude withina region surrounding the candidate peak. If the di↵erence between theminimum and the peak value is too small, the peak is rejected. In some cases,such as small ECG artifacts, there may be no remaining candidates. Usuallythere remains at least one candidate, however. The algorithm chooses theR-peak as the candidate peak with the largest amplitude.Artifacts are inevitable. The algorithm may miss some R-peaks, or de-tect false R-peaks, in the presence of ECG noise. Some may be corrected byfurther logical extensions to the algorithm. For example, candidate R-peakscould be better rejected by modeling the expected location of the next peakbased on the locations of prior peaks. In some cases, however, the ECG isso noisy that R-peaks are not present, and no additional logic can improvethe performance. Even if an algorithm works correctly 99% of the time, itwill still produce an artifact every 60 - 90 s. Artifacts are inevitable, soproper detection and correction are essential components of any clinicallyuseful nociception monitor.4.2.2 Real-Time PR from the PPGA custom designed algorithm was developed to detect pulses in the PPGsignal in real-time. It runs in O(n) (linear) time and is very simple pro-grammatically. The algorithm is called SimpleZeroCrossings.The SimpleZeroCrossings algorithm detects a pulse each time the PPGsignal amplitude crosses from a positive value to a zero or negative value.This simple algorithm is only possible because Masimo pulse oximeters wereused exclusively in the study. Masimo oximeters employ advanced adaptivefiltering algorithms to clean the PPG signal. Raw PPG signals can be verynoisy. The baseline wanders as the temperature changes, as tissue becomesmore or less dense, as the local perfusion increases or decreases, and whenthe sensor or patient’s appendage is moved. Masimo’s filtering removes allof these e↵ects, resulting in a very smooth and regular PPG signal centredaround zero (see Fig. 4.3).The time di↵erence between pulses represents the pulse period, and itsreciprocal represents the pulse rate (PR). The instantaneous PR and HRseries are both representative of autonomic state. The PR tends to be anoisier signal, due mainly to its lower sampling frequency. See the upcomingsection on the e↵ect of sampling frequency on HRV for more details; Section4.6.4. For the purposes of this work, we will consider the PR and HR984.2. Real-Time Heart Rate Measurement1 2 3 4 5 6 7 8 9 10−4−3−2−101234Time (s)Amplitude (Arbitrary Units)Figure 4.3: Example analysis with the SimpleZeroCrossings algorithm. ThePPG signal has been heavily filtered by the Masimo oximeter prior to collec-tion, allowing us to detect pulses as simple zero crossings (at the horizontalred dotted line).994.3. Real-Time Resamplingequivalent.As with the HR from the ECG, artifacts in the PR are inevitable. Thealgorithm may miss some pulses, or detect false pulses, in the presence ofPPG noise that Masimo was unable to filter out. Even if an algorithm workscorrectly 99% of the time, it will still produce an artifact every 60 - 90 s.Artifacts are inevitable, so proper detection and correction are essentialcomponents of any clinically useful nociception monitor.4.3 Real-Time ResamplingThe input signals are resampled for internal processing in iAssist. The HRsignal is resampled onto an even grid for time/frequency analysis (Section4.3.1). The CO2 signal is downsampled to match the new HR samplingfrequency (Section 4.3.2). The RR signal is upsampled to match the newHR sampling frequency (Section 4.3.3).An internal resampling frequency of 2.5 Hz is used. This frequency isdi↵erent than that used previously (see Chapter 3), because it has benefitsin a real-time environment. It is an integer fraction of the CO2 signal’ssampling frequency of 25 Hz, allowing for very simple downsampling. It alsoproduces a slight benefit in computational complexity, with fewer data pointsto analyze in real-time. There is no drawback to this reduced frequency.The maximum frequency that can be analyzed (1.25 Hz) is still significantlyhigher than any reasonable RSA frequency that is analyzed within the CRCsystem. It allows for a maximum RR of 75 breaths/min, which lies beyonda reasonable physiological range.4.3.1 HR ResamplingThe HR resampling algorithm is a real-time adaptation of Berger’s resam-pling algorithm [6]. It resamples the unevenly-spaced heart beats onto anevenly-sampled grid of arbitrary frequency.Each point in the resampled timeline exists at the centre of a local win-dow. The local window width encompasses two resampling periods:LocalWindowWidth = 2/ResamplingFreq (4.3)The algorithm calculates the number of heart periods within the localwindow, NumHPsInLocalWindow. For each point in the resampled timeline,the algorithm finds the last beat time at or before the local window start, andthe first beat time at or after the local window end. It then calculates the1004.3. Real-Time ResamplingHeart PeriodLocal WindowWidthLocal Window Start Local Window EndCurrent SampleFigure 4.4: Example HR resampling condition, where the local window fallsentirely within a single heart period.heart periods from all beats between them (including the end beats them-selves). The heart periods are stored in an array called LocalWindowHPs.There are two possible analysis conditions.If the local window falls entirely within a single heart period (see Fig.4.4), then it is calculated as:NumHPsInLocalWindow =LocalWindowWidthLocalWindowHPs(1), (4.4)where the array index “1” retrieves the first (and only) local window heartperiod.Otherwise, the local window spans two or more heart periods (see Fig.4.5), and it is calculated as:NumHPsInLocalWindow =aLocalWindowHPs(1)+ bLocalWindowHPs(end) , (4.5)where a and b are shown in Fig. 4.5. The array index “1” retrieves the firstlocal window heart period, and index “end” retrieves the last. In the rareevent that the local window also spans one or more complete heart periods,these heart periods are simply counted and added to NumHPsInLocalWindow.Finally, the algorithm calculates the resampled HR as:HRResampled = 60 x fr x NumHPsInLocalWindow2 , (4.6)1014.3. Real-Time ResamplingLocal Window Start Local Window EndCurrent SampleLocal WindowWidthHeart Period 1 Heart Period 2a bFigure 4.5: Example HR resampling condition, where the local window spanstwo heart periods. This condition can be extended to handle the conditionwhere the local window spans multiple heart periods.where fr is the resampling frequency. The factor of 60 is used to convertthe HR into units of beats/min (bpm).The Berger resampler imposes a small real-time delay. To estimate theHR at any one point in time, the algorithm requires one future heart beat.The real-time delay is thus variable, ranging anywhere from zero to one fullheart period.4.3.2 Respiration ResamplingThe CO2, O2, Paw, or Flow signal is downsampled from 25 Hz to the internalanalysis frequency of 2.5 Hz.Standard downsampling practice involves lowpass filtering and decima-tion, but this approach has some drawbacks. The lowpass filter, used forantialiasing, su↵ers edge artifacts. Windowing (e.g. applying a Hammingwindow) can reduce the edge artifacts, but causes distortions in the passbandfrequency content. Fig. 4.6 shows an example of a downsampled CO2 signal,with distortions from the lowpass filter. Filtering has a computational costthat depends on the implementation. Filtering can have complexity O(nlog n) when using fast Fourier transforms (FFTs) or complexity O(n) whenusing real-time direct form. Finally, filtering also imposes a real-time delay,1024.4. Real-Time Artifact Detectionregardless of windowing or implementation.The lowpass filter may not be necessary. When the patient is mechani-cally ventilated during general anesthesia, the respiration pattern is smoothand mostly stationary. The high frequency components in the CO2 signalare typically very small. Any aliasing will have a correspondingly small e↵ecton the resampled signal. Simple decimation can be used without incurringsignificant error (see Fig. 4.6). The complexity is significantly reduced, asa simple decimator operates in O(n) (linear) time. The algorithm does notimpose any real-time delay.4.3.3 RR ResamplingThe RR trend is upsampled using a repeater. It is recorded from the monitorat 0.1 Hz. The trend value is simply repeated at the internal samplingfrequency until a new value arrives from the monitor. This algorithm isvery simple, operating in O(n) (linear) time. It does not impose any real-time delay.4.4 Real-Time Artifact Detection4.4.1 HR Artifact DetectionA custom designed HR artifact detection algorithm was developed for op-eration in real-time. It runs in O(n) (linear) time and is simple program-matically. Its concept is similar to the SimpleDiff algorithm described inSection 4.2.1.The DetectHRArtifacts SimpleDiff algorithm detects HR artifacts bycalculating a simple di↵erentiation. This is equivalent to performing a firstlevel wavelet decomposition with a Haar wavelet. HR artifacts have verystrong high frequency components. These components are much more pow-erful than anything else in the HR signal. The algorithm iterates over eachHR sample, and calculates the di↵erence from the previous sample and itssquare:HRDi↵ = HR( n ) - HR( n - 1 ), (4.7)HRDi↵Squared = HRDi↵2, (4.8)where n is the current index. Squaring the di↵erences makes all valuespositive, and increases the contrast between small and large changes. Bykeeping only high frequency components of the HR, the algorithm retains1034.4. Real-Time Artifact Detection10 20 30 40 50 600123456CO 2 (%)Original CO210 20 30 40 50 600123456CO 2 (%)CO2 Downsampled with Lowpass Filtering & Decimation10 20 30 40 50 600123456CO 2 (%)CO2 Downsampled with Decimation OnlyTime (s)Figure 4.6: Example of CO2 downsampling. The top plot is the originalCO2 signal at 300 Hz (note that the source was recorded at 25 Hz, butthe S/5 software upsamples to 300 Hz). The CO2 signal is somewhat non-stationary, as the patient is fighting the ventilator. The middle plot is theCO2 signal resampled to 2.5 Hz using lowpass filtering and decimation (cal-culated with the resample function in Matlab). The bottom plot is theCO2 signal resampled to 2.5 Hz using decimation only (calculated with thedownsample function in Matlab. Note that the lowpass filtering introducessmall distortions, unlike the simple decimation-only approach. Simple dec-imation outperforms the standard approach, even under the nonstationaryconditions shown here. Regardless, the di↵erence between the two results issmall.1044.4. Real-Time Artifact Detectionthe artifacts while rejecting other normal HR features such as RSA andMeyer waves (see Fig. 4.7).When the squared di↵erence exceeds a threshold, calledHRDiffThreshold, the algorithm detects that it is in the region of an arti-fact. In this region, there may be multiple points that exceed the threshold,as the HR may rise or fall for several samples before returning to baseline.There may also be points in this region that do not exceed the threshold.The algorithm thus allows for threshold gaps in the region of an artifact,with maximum width specified by a parameter called MaxGapWidth. Whenthe maximum gap width is exceeded, the algorithm detects that it has leftthe region of an artifact. It retains the start and end indices of the region.This algorithm can be applied to both the HR (derived from the ECG)or the PR (derived from the PPG).4.4.2 CO2 Artifact DetectionA custom designed CO2 artifact detection algorithm was developed for op-eration in real-time. It runs in O(n) (linear) time and is very simple pro-grammatically.The DetectCO2Artifacts SimpleZeros algorithm detects artifacts suchas CO2 rezeroing by finding zero points in the signal. Zeros rarely occur ina clean CO2 signal. A small number of zero values may occur during theinspiration phase when the RR is very low (< 8 breaths/min). Thus, asingle zero does not necessarily indicate an artifact. The algorithm waits fora predetermined period of successive zeros before flagging the CO2 regionas an artifact. This tuning parameter, called MaxCleanZeroLengthSecs, ismeasured in seconds and is user-tunable.The algorithm has a variable real-time delay. In normal clean conditions,the algorithm imposes zero delay. In conditions when the RR is very low,any clean zero points will be delayed until the next non-zero point arrives.At the start of a CO2 rezeroing artifact, the algorithm will impose a delay ofup to the value specified in MaxCleanZeroLengthSecs. Later in the artifact(once it has been identified as an artifact) the delay returns to zero.4.4.3 RR Artifact DetectionWhen the Datex/Ohmeda S/5 monitor is unable to determine the RR, itreports a negative error value. Di↵erent negative values indicate di↵erenterror conditions, but for the purpose of real-time CRC analysis the specificcause of the error is irrelevant.1054.4. Real-Time Artifact Detection10 20 30 40 50406080100120140160Heart Rate(bpm)Heart Rate10 20 30 40 50−40−2002040Difference (bpm)Heart Rate Difference10 20 30 40 50050010001500Squared Difference(bpm2)Heart Rate Squared DifferenceTime (s)Figure 4.7: Example analysis with the DetectHRArtifacts SimpleDiff al-gorithm. The top plot shows a HR with three typical artifacts. The firstartifact (left) results from a false beat, second (middle) results from a mis-placed beat, and third (right) results from a missed beat. Since the artifactsare so large in magnitude, it is dicult to discern other details of the HR.The same HR is shown later in Fig. 4.9, with full details visible. The middleplot shows how the artifacts have been isolated, leaving only the relevanthigh frequency content. The bottom plot shows the squared di↵erence andHR artifact detection threshold (the horizontal red dotted line). Verticalred bars denote the detected artifact regions.1064.5. Real-Time Artifact CorrectionThe DetectRRArtifacts algorithm detects negative RR values, and re-ports them as artifacts. It runs in O(n) (linear) time and is very simpleprogrammatically. It imposes no real-time delay.4.5 Real-Time Artifact CorrectionThe resampled HR and respiration signals are the simplest target for artifactinterpolation. The CRC algorithm’s analyzing and smoothing filters act tospread out the artifact in time. As a result, the artifact’s e↵ect is greater onthe filtered HR and respiration, and greater still on the final CRC output(see Fig. 4.8). Interpolating artifacts in the original resampled HR andrespiration signals, before they are processed by the CRC algorithm, thusgives the best chance of reliably correcting artifacts.4.5.1 HR Artifact CorrectionCubic InterpolationA custom designed HR artifact correction algorithm was developed usingcubic interpolation. Cubic interpolation is performed with a third degreepolynomial and its derivative:f(x) = ax3 + bx2 + cx+ d,f 0(x) = 3ax2 + 2bx+ c. (4.9)Solving for the polynomial’s parameters, we obtain:a = 2f(0) 2f(1) + f 0(0) + f 0(1),b = 3f(0) + 3f(1) 2f 0(0) f 0(1),c = f 0(0),d = f(0). (4.10)Since the derivative is unknown for the sampled HR, it must be estimatedfrom the data: f 0(0) = (f(1) f(1)) /2,f 0(1) = (f(2) f(0)) /2. (4.11)The polynomial that results from estimated derivatives is called a Catmull-Rom spline.1074.5. Real-Time Artifact Correction30 60 90 120 150 1807274767880828486Heart Rate(bpm)Heart Rate30 60 90 120 150 180−10−50510Heart Rate(bpm)Heart Rate Filtered30 60 90 120 150 180020406080100CRCCRCTime (s)Figure 4.8: Example of an HR artifact spreading out in time during a CRCanalysis. The top plot shows an HR with an artifact with width 2 s. Themiddle plot shows the HR after being processed with the CRC analyzingfilter. The solid blue line is the real component, and the dashed green lineis the imaginary component. The bottom plot shows the CRC output afterthe HR has been processed with the analyzing and smoothing filters. Thered regions approximate the artifact width. The artifact width increaseswith each stage of processing.1084.5. Real-Time Artifact CorrectionThe cubic function parameters are thus estimated from four standard-ized points x1 = -1, x2 = 0, x3 = 1, and x4 = 2. In discrete time, thesepoints are labeled as n1, n2, n3, and n4, respectively. The cubic function isconfigured such that the artifact exists on the domain (x2, x3). This domainmust be scaled to fit individual artifacts, which can have di↵erent widthsdepending on artifact severity and sampling frequency. If the artifact startsat ArtifactStart and ends at ArtifactEnd, then the artifact width andcorresponding four standardized points (in discrete time) are given by:ArtifactWidth = ArtifactEndArtifactStart + 1,n1 = ArtifactStartArtifactWidth 2,n2 = ArtifactStart 1,n3 = ArtifactEnd + 1,n4 = ArtifactEnd + ArtifactWidth + 2. (4.12)This arrangement results in a real-time delay of 2 x ArtifactWidth + 1samples. This is the time from ArtifactStart to n4.We can reduce the time delay by warping the timeline. The standardizedpoints become x1 = 0-, x2 = 0, x3 = 1, and x4 = 1+, where “0-” denotesthe time of the last sample before 0 and “1+” denotes the time of the nextsample after 1. The four standardized points (in discrete time) are thengiven by: n1 = ArtifactStart 2,n2 = ArtifactStart 1,n3 = ArtifactEnd + 1,n4 = ArtifactEnd + 2. (4.13)This warped timeline reduces the real-time delay to just ArtifactWidth +2 samples. Despite the smaller delay, it produces nearly identical results tothe standard Catmull-Rom spline formulation in Eqn. 4.12.Substituting Eqns. 4.13 and 4.11 into Eqn. 4.10, the polynomial’s pa-rameters then become:a = 0.5 HR(n1) + 1.5 HR(n2) 1.5 HR(n3) + 0.5 HR(n4),b = HR(n1) 2.5 HR(n2) + 2 HR(n3) 0.5 HR(n4),c = 0.5 HR(n1) + 0.5 HR(n3),d = HR(n2). (4.14)1094.5. Real-Time Artifact CorrectionFinally, by substituting the parameters from Eqn. 4.14 into f(x) inEqn. 4.9 and evaluating from x = ArtifactStart to x = ArtifactEnd, weestimate the missing samples in the artifact region.Examples of cubic interpolation are shown in Fig. 4.9.ExperimentThe HR artifact correction algorithm was tested using simulated artifacts.The anesthetic bolus data segments were used as gold standard clean HRsignals (see Section 3.7). Each of the anesthetic bolus data segments containsa baseline and response period known to be free of HR artifacts (see Fig.3.16). Each anesthetic bolus event thus provides two gold standard HRsegments, each 60 s wide. The segments should provide a good randomizedsample of possible HR conditions, including di↵erent means and variances,di↵erent levels of nociception, di↵erent respiratory phases, etc. Artifactsranging from 1 - 10 samples in length (0.25 to 2.5 s) were simulated atthe beginning of each segment, and corrected using the cubic interpolationalgorithm. Their e↵ect on the HR and CRC nociception index were thenmeasured.The peak HR artifact correction error was measured. Artifact correctionerror was measured as the di↵erence between the interpolated and goldstandard HRs for each segment at each artifact length. The peak error wascalculated for each segment at each artifact length, then averaged across allsegments.The CRC nociception index error was measured. CRC was calculatedover each of the anesthetic bolus data segments, both with the original goldstandard HR and with the simulated HR artifacts corrected. Artifact cor-rection error was measured as the di↵erence between the CRC calculatedfrom the interpolated and gold standard HRs for each segment at each arti-fact length. The peak error was calculated for each segment at each artifactlength, then averaged across all segments.ResultsA total of 57 anesthetic bolus events were analyzed, providing 111 goldstandard HR segments. Three segments had unusually large peak HR error,and were excluded as outliers.The mean peak HR error (averaged across all 111 segments) ranged from0.15 to 1.5 bpm, resulting from artifact lengths of 0.25 to 2.5 s, respectively.The error gradually increased as the artifact length increased (see Fig. 4.10).1104.5. Real-Time Artifact Correction10 20 30 40 50758085Heart Rate(bpm)Gold Standard Heart Rate10 20 30 40 50758085Heart Rate(bpm)Heart Rate with Artifact Gaps10 20 30 40 50758085Heart Rate(bpm)Interpolated Heart Rate10 20 30 40 50−4−202Heart Rate Difference(bpm)Interpolation ErrorTime (s)Figure 4.9: Example cubic interpolation of HR artifacts. The top plot showsa gold standard HR without artifacts. The second plot shows the same HR,with simulated artifacts removed. These are the same artifacts shown earlierin Fig. 4.7. The third plot shows the HR with artifacts corrected using cubicinterpolation. The bottom plot shows the HR di↵erence between the goldstandard and corrected HRs. Vertical red bars denote the artifact regions.1114.5. Real-Time Artifact CorrectionThe mean peak CRC error (averaged across all 111 segments) rangedfrom 0.12 to 6.0, resulting from artifact lengths of 0.25 to 2.5 s, respectively.The error gradually increased as the artifact length increased (see Fig. 4.11).The CRC error increased sharply in response to a corrected artifact, thengradually returned to zero in less than 60 s (see Fig. 4.12).DiscussionA custom designed HR artifact correction algorithm was developed for oper-ation in real-time. It uses cubic interpolation to correct small HR artifacts.The algorithm imposes a real-time delay of just ArtifactWidth + 2 samples,runs in O(n) (linear) time, and is very simple programmatically.Cubic interpolation can only correct short artifacts. The cubic functioncreates a smooth curve. When interpolating a sinusoid, it cannot accuratelymodel more than half a cycle. Since the RSA is the only important HRcomponent for measuring CRC, the algorithm works best when interpolatingno more than half a respiration cycle. As the respiration frequency decreases,the respiration cycle length increases and the algorithm can better handlelonger artifacts.Most artifacts tend to be short, lasting just two cardiac cycles. Suchartifacts result from a single false beat, misplaced beat, or missed beat. Theartifact duration is thus typically less than the 2.5 s tested herein (equivalentto a HR of 48 bpm). Longer artifacts can occur, however, and may not beaccurately corrected with the cubic interpolation algorithm.This algorithm can be applied to both the HR (derived from the ECG)or the PR (derived from the PPG).A more advanced artifact correction algorithm may be developed basedon predictive filtering using maximum entropy [76, 77]. This predictivefiltering scheme may prove very e↵ective in cases where the HR exhibits aclean sinusoidal pattern, such as during periods of strong RSA. We leave itto future work to implement this algorithm.4.5.2 CO2 Artifact CorrectionSignal RepeaterA custom designed CO2 artifact correction algorithm was developed usinga signal repeater. The algorithm converts the RR trend into a respirationperiod, which indicates the number of samples in a respiration cycle. It thenextracts the last full respiration cycle from the CO2 signal, immediatelypreceding the artifact region. The algorithm then repeats this respiration1124.5. Real-Time Artifact Correction0.25 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5123456Artifact Length (s)HR Error (bpm)Figure 4.10: HR interpolation peak error results. Each thin line representsthe peak error at a range of artifact lengths for a single HR segment. Thethick red line represents the mean peak error, averaged across all segments.As the artifact width increases, the peak error also increases.1134.5. Real-Time Artifact Correction0.25 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.50102030405060708090100Artifact Length (s)CRC ErrorFigure 4.11: CRC peak error results from HR interpolation. Each thinline represents the peak error at a range of artifact lengths for a singleCRC segment. The thick red line represents the mean peak error, averagedacross all segments. As the artifact width increases, the peak error alsoincreases. The vertical scaling encompasses the entire possible range of theCRC nociception index.1144.5. Real-Time Artifact Correction5 10 15 20 25 30 35 40 45 50 55 600102030405060708090100CRC Mean ErrorTime (s)Figure 4.12: CRC error results in time from HR interpolation. Each thinline represents the mean error across all CRC segments at a single artifactlength, over time. The HR artifact begins at t = 0 (far left of the plot). Asthe artifact width increases, the error lines climb higher (greater error). Thevertical scaling encompasses the entire possible range of the CRC nociceptionindex.1154.5. Real-Time Artifact Correctioncycle as many times as is necessary to fill the artifact region, truncating anyoverlap at the end.The signal repeater imposes zero real-time delay. When an artifact isencountered, it can be immediately corrected with CO2 values from thepreceding respiration cycle.An example of signal repeater interpolation is shown in Fig. 4.13.ExperimentThe CO2 artifact correction algorithm was tested using simulated artifacts.The anesthetic bolus data segments were used as gold standard clean CO2signals (see Section 3.7). Each of the anesthetic bolus data segments containsa baseline and response period known to be free of CO2 artifacts (see Fig.3.16). Each anesthetic bolus event thus provides two gold standard CO2segments, each 60 s wide. The segments should provide a good randomizedsample of possible CO2 conditions, including di↵erent exhaled CO2 levels,di↵erent RRs and PPaws, di↵erent respiratory phases, etc. Artifacts 10 s inlength were simulated at the beginning of each segment, and corrected usingthe signal repeater algorithm. Their e↵ect on the CO2 and CRC nociceptionindex were then measured.The peak CO2 artifact correction error was measured. Artifact correc-tion error was measured as the di↵erence between the interpolated and goldstandard CO2 signals for each segment. The peak error was calculated foreach segment, then averaged across all segments.The CRC nociception index error was measured. CRC was calculatedover each of the anesthetic bolus data segments, both with the original goldstandard CO2 and with the simulated CO2 artifacts corrected. Artifact cor-rection error was measured as the di↵erence between the CRC calculatedfrom the interpolated and gold standard CO2 signals for each segment. Thepeak error was calculated for each segment, then averaged across all seg-ments.ResultsA total of 57 anesthetic bolus events were analyzed, providing 114 goldstandard CO2 segments.The mean peak CO2 error (averaged across all 114 segments) was 1.5%.The error distribution was bimodal, with most errors very small (< 0.5%)and some very large (> 5.0%) (see Fig. 4.14).1164.5. Real-Time Artifact Correction10 20 30 40 50 60123456CO 2 (%)Gold Standard CO210 20 30 40 50 600246CO 2 (%)CO2 with Simulated Rezeroing Artifact10 20 30 40 50 60123456CO 2 (%)Interpolated CO210 20 30 40 50 60−0.0500.050.1CO 2 Difference (%)Interpolation ErrorTime (s)Figure 4.13: Example repeater interpolation of a CO2 rezeroing artifact.The top plot shows a gold standard CO2 signal without artifacts. Thesecond plot shows the same CO2 signal, with a simulated rezeroing artifactlasting 10 s. The third plot shows the CO2 with artifact corrected usinga signal repeater. The bottom plot shows the CO2 di↵erence between thegold standard and corrected CO2 signals. The vertical red bar denotes theartifact region.1174.5. Real-Time Artifact Correction10 20 30 40 50 60 70 80 90 100 1101234567Artifact NumberCO 2 Error (%)Figure 4.14: CO2 interpolation peak error results. Each blue circle repre-sents the peak error for a single CO2 segment. The horizontal dotted redline represents the mean peak error, averaged across all segments.The mean peak CRC error (averaged across all 114 segments) was 3.8.The CRC error increased in response to a corrected artifact, then graduallyreturned to zero in less than 60 s (see Fig. 4.15).DiscussionA custom designed CO2 artifact correction algorithm was developed for oper-ation in real-time. It uses a periodic signal repeater to interpolate rezeroingartifacts. The algorithm imposes zero real-time delay, runs in O(n) (linear)time, and is very simple programmatically.The signal repeater algorithm assumes that the CO2 signal is stationarybefore, during, and immediately following the artifact. It can only functionproperly when the patient is mechanically ventilated. Spontaneous ventila-tion is too erratic, making it impossible to extract a clean respiratory cyclefor repetition. Even during mechanical ventilation, nonstationary featurescan occur and will lead to incorrect interpolation. Possible nonstationaritiesinclude changes in RR or warped breaths, caused by the patient fightingthe ventilator (see e.g. the baseline period in Fig. 3.16, highlighted in red).Indeed, the bimodal distribution of errors shown in Fig. 4.14 results from1184.5. Real-Time Artifact Correction5 10 15 20 25 30 35 40 45 50 55 600102030405060708090100CRC Mean ErrorTime (s)Figure 4.15: CRC error results in time from CO2 interpolation. The linerepresents the mean error across all CRC segments, over time. The CO2artifact begins at t = 0 (far left of the plot). The vertical scaling encompassesthe entire possible range of the CRC nociception index.1194.6. Real-Time Signal Source Selectionthe di↵erence between stationary and nonstationary CO2 signals. Station-ary signals produce the points with small error, while nonstationary signalsproduce the points with large error. When the CO2 signal is nonstationary,the interpolation may produce discontinuities as the respiratory cycles donot align perfectly. These may lead to significant errors in the correctedCO2 signal and resulting CRC nociception index. Furthermore, changes inthe exhaled CO2 concentration (e.g. changing EtCO2) produce a nonsta-tionary component in the signal amplitude. While typically a much smallere↵ect than changes in RR, it still produces some small discontinuities in theinterpolated signal. The points with small error in Fig. 4.14 are caused bythis e↵ect.The algorithm can correct CO2 artifacts of any length, but longer arti-facts increase the probability of encountering nonstationary features. Reze-roing artifacts tend to be short, on the order of 10 - 12 s.A more advanced artifact correction algorithm may be developed basedon predictive filtering using maximum entropy [76, 77]. This predictivefiltering scheme may prove very e↵ective for the CO2 and O2 signals, whichtypically have a clean sinusoidal pattern during mechanical ventilation. Weleave it to future work to implement this algorithm.4.5.3 RR Artifact CorrectionA custom designed RR artifact correction algorithm was developed for op-eration in real-time. It simply repeats the last clean value to replace anyartifact value. The algorithm runs in O(n) (linear) time, is very simpleprogrammatically, and imposes zero real-time delay.4.6 Real-Time Signal Source Selection4.6.1 We Need Secondary Redundant Signal SourcesThe instantaneous HR is a very noisy signal. Noise in the source signal (e.g.an ECG) often leads to missed or false beats, which manifest as artifactsin the instantaneous HR measurements. It can be impossible to accuratelydetect every R-peak in the ECG, while simultaneously never detecting anyfalse peaks.Movement artifacts can lead to incorrect HR measures. The patient isunlikely to move much autonomously during surgery when properly anes-thetized. However, the surgical team may occasionally lean on the patient’schest, pull on the ECG electrode cables, or roll the patient over onto his/her1204.6. Real-Time Signal Source Selectionside. During these movement events, the ECG baseline may swing verystrongly upward or downward. These swings may occur at relatively lowfrequencies, which usually do not interfere with R-peak detection. Theymay also occur at high frequencies, which may cause false R-peaks to bedetected and thus cause artifacts in the instantaneous HR series.Electrocautery is an example of a very challenging noise source in theECG. Surgeons typically no longer use traditional scalpels for making in-cisions. Instead, they use a tool called an electrosurgical unit (ESU) thatsimultaneously cuts the tissue and cauterizes the blood vessels. The ESUincludes a large grounding plate that has a low impedance contact withthe patient’s skin. The surgeon uses a “knife” (a pen-shaped device with afine point) to make the incisions. The ESU passes high frequency currentthrough the knife tip, through the patient’s tissue, and into the groundingplate. Since the grounding plate is so large, the current density is low andthe contacting tissue is undamaged. Since the knife is so small, the currentdensity is high and the contacting tissue heats up significantly. In so doing,the ESU ablates the tissue and coagulates the blood. While the ESU is avery useful tool for the surgeon, the high frequency current injected intothe patient’s body wreaks havoc with electrical sensors such as the ECG.Fig. 4.16 shows an example of electrocautery noise corrupting an ECG. Theelectrocautery noise is typically at a very high frequency, and can mimic theshape of the ECG R-peak. Beat detection algorithms can often mistake theelectrocautery noise for valid beats, resulting in significant artifacts in theinstantaneous HR series. These artifacts then cause significant errors in theCRC nociception index measures. While it is relatively easy to detect theelectrocautery noise in the ECG [10], removing it is much more challenging.The ECG and the PPG are unlikely to experience the same noise at thesame time. While the ECG sensors are axed to the patient’s thorax, thePPG sensor is typically axed to a distal appendage such as a finger, toe,or earlobe. Movement artifacts in the ECG are often not manifested in thePPG. Furthermore, the PPG is measured as an optical signal, so it is lessa↵ected by electrocautery noise than the ECG. The PPG can be a↵ectedby other noise sources, including movement of the limbs (if it is attached toe.g. a finger) or ambient light. However, these noise sources typically do nota↵ect the ECG. Thus, at any given point in time it is likely that at least oneof the ECG or PPG signals should be clean enough to measure an accurateinstantaneous HR series.Likewise, respiration signals can be a↵ected by artifact. Particularlyproblematic is re-zeroing in the capnometer, which occasionally leads to abrief period with no CO2 and O2 signals. Spirometry signals (e.g. Paw or1214.6. Real-Time Signal Source Selection1 2 3 4 5−12−10−8−6−4−202468Time (s)AmplitudeFigure 4.16: Example of electrocautery noise corrupting an ECG signal.Electrocautery noise is present for the first 3 s of the signal.Flow) can be used in place of capnometry during these periods of artifact.Fig. 4.17 shows an example of the spirometry signal availability duringcapnometry rezeroing.Redundant signal sources allow for clean and uninterrupted nociceptionmonitoring during challenging noisy conditions.4.6.2 Real-Time Signal Quality EstimationA signal quality index (SQI) can be derived for the HR and respiration sig-nals. During periods of artifact, it is impossible to know exactly how wellthe signal was reconstructed by the artifact correction algorithms (see Sec-tion 4.5). However, due to the ease and high reliability of artifact detection,it is possible to confidently know how much of a given signal window is po-tentially corrupted. We define the SQI as the percentage of a given signalwindow that is not a↵ected by artifacts.We have developed an SQI estimation algorithm. It uses a variable calledTotalNumSamples which stores the total number of samples in the mostrecent signal window required for CRC analysis. It uses a second variablecalled NumCleanSamples which stores the number of clean (i.e. not a↵ectedby artifact) samples within this same window. The SQI is then calculated1224.6. Real-Time Signal Source Selection10 20 30 40 50 60 70 80 90010203040506070Time (s)Concentration (%)  CO2O210 20 30 40 50 60 70 80 90−15−10−5051015Time (s)Pressure and ∆ Volume  PawFlowFigure 4.17: Example of rezeroing in the capnometer (top plot, from 35 -50 s). The CO2 and O2 signals are temporarily unavailable. The spirometrysignals (bottom plot) are still available throughout.1234.6. Real-Time Signal Source Selectionas:SQI = 100 xNumCleanSamplesTotalNumSamples, (4.15)4.6.3 Real-Time Signal Source SelectionThe simplest source selector chooses the redundant signal with the highestSQI. Redundant HR signals include the HR from the ECG and the PRfrom the PPG. Redundant respiration signals include the CO2 or O2 fromcapnometry and the Paw or Flow from spirometry.A more nuanced selection algorithm may also consider other signal qual-ity metrics. For example, in the case of the HR signal the source samplingfrequency has a strong e↵ect on the accuracy of the HR measures. Thesource signal sampling frequency thus a↵ects the accuracy of the CRC no-ciception index, but this is not reflected in the SQI as we have defined ithere. See the upcoming section on the e↵ect of sampling frequency on HRVfor more details; Section 4.6.4. A more nuanced selection algorithm mayconsider the tradeo↵ between SQI and these other metrics, though we leavethis to future work.Theoretically, individual HR/PR samples can be individually swappedin and out, such that the signal from the highest frequency source (e.g.the ECG) is used except for those samples that are a↵ected by artifact.However, in practice there are synchronization challenges that make thissubstitution dicult. The pulses in the PPG signal always lag the beatsin the ECG signal, because the blood pressure wave takes time to travelfrom the heart to the periphery where the PPG sensor detects them. Thisphenomenon is called pulse transit time [65], and it varies with di↵erentconditions in the patient’s body. We cannot simply assume the pulses andbeats occur at the same time, and we cannot assume they will have a fixedo↵set. It may be possible to dynamically align the beats and pulses with acorrelation method. However, we leave it to future work to implement thismore advanced selection algorithm.4.6.4 The E↵ect of Sampling Frequency on HRVIntroductionThe relation between ECG sampling frequency and HR signal quality isimportant, yet it remains poorly understood even after decades of work onHRV analysis. Investigation of this relationship began in the early 1990s.When the Task Force of the European Society of Cardiology released their1244.6. Real-Time Signal Source Selectionseminal 1996 report on HRV [71], they included only a brief mention of ECGsampling frequency, reproduced here verbatim: “The sampling rate must beproperly chosen. A low sampling rate may produce jitter in the estimationof the R-wave fiducial point, which alters the spectrum considerably. Theoptimal range is 250 to 500 Hz or perhaps even higher, while a lower samplingrate (in any case 100 Hz) may behave satisfactorily only if an algorithm ofinterpolation (parabolic) is used to refine the R-wave fiducial point.” Thereport o↵ers no concrete justification for its choice of bounds. Its wording isvague, and falls short of making a firm recommendation for an ECG samplingfrequency. Furthermore, the report considers only spectral analysis of HRV.The report’s conclusions are based primarily on theoretical work from 1990[50].We wish to test the Task Force’s statement, to measure the e↵ect ofsource signal and sampling frequency on the HRV-based nociception indicesdevised in Chapters 2 and 3. Is a 300 Hz ECG and/or a 100 Hz PPG signaladequate for monitoring nociception with CRC?MethodSome cases from the Nociception Study (see Chapter 3) dataset included a900 Hz ECG recorded from the NeuroSense monitor, a 300 Hz ECG recordedfrom the Datex/Ohmeda monitor, and a 100 Hz PPG recorded from theDatex/Ohmeda monitor. These data can provide estimates of the e↵ect ofsource signal and sampling frequency on the nociception indices.Data Preparation Cases were inspected to find clean segments of mul-tiresolution data for analysis. For each case, an HR/PR series was derivedfor the NeuroSense and Datex/Ohmeda ECG and PPG signals. ECG R-peaks and PPG zero crossings were detected and manually inspected andcorrected to arrive at the best possible HR/PR series. In certain circum-stances beats were impossible to discern even visually (e.g. during occasionalzeroing in the PPG signal). Clean segments with duration at least 5 min-utes were visually identified. Segments were only used if all of the followingsignals were clean: the HR from NeuroSense and Datex/Ohmeda ECG, thePR from PPG, and the CO2 respiration.Theoretical multiresolution HR series were calculated to estimate theexpected HR error. The 900 Hz ECG signal segments were downsampledusing Matlab’s resample function, which applies an anti-aliasing low-passfilter then decimates. The ECG was downsampled by a factor of 3 (to 300Hz), 9 (to 100 Hz), and 27 (to 33 Hz). Beat times were likewise decimated,1254.6. Real-Time Signal Source Selectionthen corrected to the nearest local maximum in the downsampled ECGs toaccount for any misplacement.Data segments were synchronized across the actual 900 Hz ECG, 300 HzECG, and 100 Hz PPG. For each segment, the HR/PR series were visuallyinspected to identify a synchronization point. A typical synchronizationpoint consisted of a sharp change in HR, often a transient deceleration.The synchronization point was meant to be easily and uniquely identifiedin all three HR/PR series. Synchronization allowed the HR/PR series to beanalyzed and directly compared. Synchronizing locally within each segment(rather than in each case, which often comprised multiple segments) helpedreduce the e↵ect of signal drift. The sampling frequencies of the di↵erentsignals are not known with perfect precision. If signals were synchronizednear the start of a case, and a segment was analyzed near the end of thecase (with hours of data between), the signals may no longer be perfectlyaligned. For example, a sampling frequency imprecision of just 0.1% betweentwo signals would lead to a drift of nearly one second every fifteen minutes.Data Analysis The HR/PR error was calculated from the theoretical andactual signals. The HR was derived from the 900 Hz ECG and treated as areference gold standard. The HR was derived from the theoretical ECGs at300 Hz, 100 Hz, and 33 Hz and were compared to the gold standard. TheHR/PR was derived from the actual 300 Hz ECG and 100 Hz PPG andwere also compared to the gold standard. The di↵erences were calculatedfor each beat of each segment. The overall root mean squared error (RMSE)was calculated for each theoretical and actual signal.The HRV-based nociception index error was calculated from the theoret-ical and actual signals. The nociception indices from CRC, SDNN, RMSSD,LF, LFnu, HF, HFnu, and LF/HF were calculated as described in Chapters2 and 3. The nociception indices were derived from the 900 Hz ECG andtreated as reference gold standards. The nociception indices were derivedfrom the theoretical ECGs at 300 Hz, 100 Hz, and 33 Hz and were com-pared to the gold standard. The nociception indices were derived from theactual 300 Hz ECG and 100 Hz PPG and were also compared to the goldstandard. The di↵erences were calculated for each sampling point of eachsegment. The overall root mean squared error (RMSE) was calculated foreach nociception index from each theoretical and actual signal.1264.6. Real-Time Signal Source SelectionMeasure RMSE at RMSE at RMSE at300 Hz 100 Hz 33 HzHR 0.20 bpm 0.66 bpm 2.5 bpmCRC 2.5 8.2 24LF 0.25 0.77 4.0HF 0.89 4.2 15LFnu 2.1 12 34HFnu 1.5 4.6 13LF/HF 0.82 3.7 12SDNN 0.075 0.71 4.7RMSSD 0.40 3.6 17Table 4.1: Theoretical error in the HR and the HRV-based nociceptionindices arising from the sampling frequency of the source ECG signal.ResultsA total of 26 clean multiresolution data segments from eight subjects wereidentified in the data. The theoretical error increased for all measures asthe sampling frequency of the source ECG decreased. Likewise, the actualerror was greater when using the 100 Hz PPG signal than when using the300 Hz ECG signal. The errors were slightly greater in the actual 300 HzECG as compared to the theoretical one. The errors were much greater inthe actual 100 Hz PPG than in the theoretical 100 Hz ECG. An example ofthe e↵ect of the source signal on the HR series is shown in Fig. 4.18 and onthe HRV-based nociception indices in Fig. 4.19. Full results are shown inTables 4.1 and 4.2.DiscussionThe results show that the error in the HR/PR and in the HRV-based no-ciception indices increases as the sampling frequency of the source signaldecreases. The results also show that the actual error was greater than thetheoretical error.The results suggest that an ultra-high-frequency ECG signal is not nec-essary for reasonable accuracy in nociception monitoring. They also suggestthat the PPG may not be an adequate substitute for the ECG. These find-ings seem to confirm the statement from the Task Force HRV document[71].The CRC nociception index su↵ered greater degradation from reduced1274.6. Real-Time Signal Source SelectionMeasure RMSE from RMSE from300 Hz ECG 100 Hz PPGHR 0.34 bpm 1.9 bpmCRC 3.4 16LF 0.29 3.0HF 1.1 14LFnu 2.2 32HFnu 1.7 11LF/HF 0.96 11SDNN 0.097 5.9RMSSD 0.47 15Table 4.2: Actual error in the HR/PR and the HRV-based nociception in-dices arising from the source signal and its sampling frequency.5 10 15 20 25 30 35 40 45 50889092949698Beat NumberHR/PR (bpm)  HR from ECG at 900 HzHR from ECG at 300 HzPR from PPG at 100 HzFigure 4.18: Example of HR series derived from a NeuroSense ECG at 900Hz, a Datex/Ohmeda ECG at 300 Hz, and a Datex/Ohmeda PPG at 100Hz.1284.6. Real-Time Signal Source Selection125 250 375 500 62580859095HR / PR (bpm)125 250 375 500 625020406080100SDNNNociception Index125 250 375 500 625020406080100LF/HFNociception Index125 250 375 500 625020406080100CRCNociception IndexTime (s)Figure 4.19: Example of sampling frequency and source signal e↵ect on theHRV-based nociception indices. Top plot: HR; Second plot: SDNN noci-ception index; Third plot: LF/HF nociception index; Bottom plot: CRCnociception index. In all plots, the blue line is generated from the Neu-roSense ECG at 900 Hz, the green line from the Datex/Ohmeda ECG at300 Hz, and the red line from the Datex/Ohmeda PPG at 100 Hz. Notethat in all plots the blue and green lines are very close together (small error)but the red line is farther away (larger error).1294.6. Real-Time Signal Source Selectionsampling frequency than did most other indices. The degradation was ev-ident in both the theoretical signals and actual signals. While the exactcause is unknown, it is likely that the reduced time resolution of low fre-quency signals reduces the accuracy of the high frequency HR content, suchas RSA. That is, the HRV feature required for CRC is most susceptible toerror when the source signal sampling frequency is reduced.Some of the increased error in the actual 300 Hz ECG as compared to thetheoretical 300 Hz ECG may have been caused by signal drift. The two ECGsources may have exhibited some small drift over the data segments, whichtypically had length of 5 - 10 minutes. The drift would have manifested itselfas a small (though not true) error between the two HR series. Nevertheless,the errors diverged by at most 50%, suggesting that the theoretical modelis a reasonable prediction of the actual results.The greatly increased error in the actual 100 Hz PPG as compared tothe theoretical 100 Hz ECG may be explained by physiological and/or tech-nological factors. It is unlikely that signal drift from sampling frequencyimprecision alone would explain the error. There may have been an addi-tional source of pulse time variability in the PPG as compared to the ECG.The ECG is a nearly instantaneous measure of the heart’s beating. ThePPG, in contrast, was measured at an extremity (typically a finger) andwill thus su↵er a delay between the beat time and the pulse time. Thisdelay, called pulse transit time, is a↵ected by the patient’s blood pressure[65]. Changes in blood pressure (blood pressure variability) may have ledto changes in pulse transit time. These changes can be thought of as anadditional source of error in the PR that is not present in the HR and wouldnot have been accounted for in the theoretical results. There may also be atechnological explanation for the increased error. The PPG signal is heavilyfiltered and digitally manipulated before being recorded for our experiments.The Datex/Ohmeda monitor does not record the raw PPG signal. Instead,it uses an oximeter module from Masimo which employs advanced filters toremove noise, baseline drift, motion artifacts, and more. Oximeters oftenuse adaptive filters that strongly distort the waveform. The internal signalmay have a low sampling frequency, as most of the relevant informationlies below 5 Hz. The actual internal sampling frequency is unknown, butthe PPG signal is later upsampled to 100 Hz for recording within the Da-tex/Ohmeda monitor. Finally, the S/5 Collect software upsamples the PPGsignal to match the 300 Hz ECG, which may introduce further noise.It is possible that the HR series derived from the PPG may be improvedwith more advanced signal processing techniques. Interpolation may im-prove the HR estimates, bringing them closer to the real values. Statistical1304.7. Real-Time Delayfiltering, such as with a Kalman filter, may remove some of the high fre-quency noise in the HR series, again reducing some of the error. It appearsas though errors in the HR series are magnified during HRV analysis. Forexample, an error of less than 2 bpm RMS in the HR leads to an error ofup to 32 RMS on the LFnu nociception index. Note that this is nearly onethird of the full scale of the index. A small reduction in the HR error cantherefore lead to a large reduction in the nociception index error.The experiment was limited by the signals available for processing. Afuture study could be performed to record a raw unfiltered PPG signal at100 Hz for more direct comparison with an ECG.ConclusionIn this section, we have investigated the e↵ects of source signal and samplingfrequency on the accuracy of the HRV-based nociception indices. We foundthat any error is reasonably small when using a 300 Hz ECG, but that a100 Hz PPG has much greater error. We also found that the source signaland sampling frequency have a much stronger e↵ect on the error of thenociception index from CRC than from most other methods. A 300 HzECG should be adequate for CRC, and a 100 Hz PPG should probablyonly be used when the HR signal from the ECG is corrupted by artifacts.These findings seem to confirm the statement made in the Task Force HRVdocument [71].4.7 Real-Time DelayWhen monitoring nociception, the real-time delay includes contributionsfrom many di↵erent components. While we previously looked at the delayof the CRC analyzing and smoothing filters in Section 2.4.4, these are notthe only sources of delay.The body’s nociceptors entail a delay in transmitting information fromthe damaged tissue to the brainstem. Nociceptive A fibres conduct theirsignals very quickly, in much less than a second. For the purposes of no-ciception monitoring, the A conduction time is negligible. In contrast,nociceptive C fibres conduct their signals much more slowly, taking as longas 2 s to reach the brainstem. The total real-time delay from the nociceptorsthus ranges from 0 - 2 s.The autonomic nervous system entails a delay in its e↵ect on the heart.The brainstem responds to a nociceptive stimulus by adjusting its output ofacetylcholine and norepinephrine. Acetylcholine acts rather quickly on the1314.8. Conclusionheart, taking only about 2 s to a↵ect the HRV. In contrast, norepinephrineacts more slowly on the heart, taking about 10 s to a↵ect the HRV. Thisslower action is due largely to the fact that norepinephrine is synthesizedby the adrenal gland and delivered to the heart via the blood. In contrast,acetylcholine is delivered directly via nerve fibres to the heart’s SA node.The total real-time delay from the neurotransmitter action on the heart thusranges from 2 - 10 s.The real-time CRC analyzing and smoothing filters entail a delay inmeasuring the change in HRV. As calculated in Section 2.4.4, the analyzingfilter group delay ranges from 2 - 6.4 s, depending on the respiration rate.The smoothing filter’s maximum group delay is fixed at 11.6 s. The totaldelay from measuring CRC thus ranges from 13.6 to 18 s.The total delay from nociceptive stimulus to detection in the nociceptionmonitor thus ranges from 15.6 to 30 s, depending on the type of nociception,the response of the neurotransmitters, and the patient’s respiration rate.The total delay for real-time CRC can be dramatically shorter than thedelay of a human anesthesiologist. The anesthesiologist observes changesin parameters like mean heart rate and blood pressure to detect a stressresponse to nociception. We demonstrated in earlier work that these pa-rameters have a very low sensitivity to nociception [14]. By the time theychange enough to be noticed, far more time has passed than that required forreal-time CRC. Alternatively, the anesthesiologist may observe behaviouralend points like patient movement or the patient resisting the mechanicalventilator. However, by the time these behaviours become noticeable, thepatient had likely already been experiencing some level of nociception forsome time. Real-time CRC can detect these events much sooner.4.8 ConclusionIn this chapter, we described a novel architecture for a real-time nociceptionmonitor. This is a very challenging problem and a significant contributionto the field. There are many practical challenges to implementing real-timeCRC in a working nociception monitor. Artifacts pose great diculties forreal-time monitoring. Our novel architecture includes signal pre-processing,artifact detection/rejection, and source selection, and we described examplealgorithm implementations for each of these areas. Using this architectureand algorithms, we can mitigate the e↵ects of artifacts to make for a morereliable real-time nociception monitor. Solutions such as these are requiredto bring the real-time CRC algorithm from theory into reality.1324.8. ConclusionThe working nociception monitor is shown in Fig. 4.20.1334.8.ConclusionFigure 4.20: The real-time nociception monitor, pictured during use in the operating room at British Columbia’sChildren’s Hospital. The nociception monitor is the blue outlined screen next to the anesthesiologist at right.134Chapter 5Conclusion & Future WorkWe have described in this dissertation the design, development, and clinicalvalidation of a real-time nociception monitor for use during general anesthe-sia.We have made two novel and significant contributions to the field. Thefirst is a novel filter for measuring RSA in real-time. The second is a novelend-to-end software system that integrates multiple signals and algorithmsfor robust continuous real-time nociception monitoring during general anes-thesia.We have developed a novel filter for measuring the RSA in the HRVin real-time. The filter uses information from a secondary signal source (arespiration rate derived from a respiration signal) to track the RSA as itmoves in time and frequency. It then uses this information to dynamicallytune the centre frequency and bandwidth of a band-pass filter to isolate andmeasure the RSA while excluding other components of HRV (see Section2.4). Previous work in the field of HRV assumes the RSA will fall withina strict frequency range; however, we show in this dissertation that thestandard assumption is invalid (see Sections 3.8, 3.9, and 3.10). This filteris very e↵ective at tracking the RSA in time and frequency, and it mayprovide the most robust measure of RSA yet devised.We have integrated multiple signals and algorithms together into an end-to-end software system for robust continuous real-time nociception monitor-ing during general anesthesia. The software system incorporates not onlyour dynamically-tuned RSA filter to measure nociception, but also manyperipheral algorithms for detecting and rejecting artifacts in the input sig-nals. The input signals required for real-time nociception monitoring can beextremely noisy, and artifacts are a very significant challenge. We have de-veloped methods to address the artifacts and we have deployed them withinthis software system. A block diagram of the system architecture is shownin Fig. 1.5. The details are described in Chapter 4. To our knowledge,no other nociception monitor includes such robust artifact handling usingredundant signal channels.In the remainder of this final chapter, we propose future work that may1355.1. Feedback Control of Antinociceptive Drugsbe performed to further improve nociception monitoring during general anes-thesia. We explore di↵erent experiments that may be performed to furthervalidate the nociception monitor under di↵erent clinical conditions and indi↵erent patient populations.5.1 Feedback Control of Antinociceptive DrugsThe ultimate future goal of the research is to develop a feedback controller forautomatically delivering antinociceptive drugs to the patient. CRC wouldserve as the feedback loop for the controller. In the simplest implementa-tion, the controller would compare the current CRC value to an internalsetpoint, and adjust the delivery of antinociceptive drugs accordingly. If theCRC rises above or falls below the setpoint, the controller would increase ordecrease the delivery of antinociceptive drugs, respectively. We must there-fore derive a setpoint for the controller, based on our understanding of theCRC nociception index.We can derive a CRC setpoint by estimating its value during non-nociceptive periods. The dental nociception dataset can be analyzed toderive this setpoint estimate. Data from the experiments in Chapter 3 areused to calculate a non-nociceptive baseline value. The Dental Dam and SNSBaseline periods, the Anesthetic Bolus Response period, and both Changein RR and PPaw Baseline and Response periods are considered as generallynon-nociceptive. The CRC nociception index distributions for each experi-ment, as well as the overall distribution, are shown in Figure 5.1. The overallmedian and mean CRC values were found to be 29.9 and 34.6, respectively.A CRC nociception index between 30 - 35 appears to be a reasonable firstestimate for a feedback controller setpoint. However, while we assume allperiods analyzed are non-nociceptive, there is surely some nociception inmany individual cases, and some cases have very strong nociception (seeChapter 3 for a detailed explanation). We should therefore consider themetrics derived herein as an upper bound for the setpoint.5.2 Enhanced Nociception Monitoring withSensor FusionOther parameters may be used in conjunction with HRV to better monitornociception. Sensor fusion between multiple parameters may provide greatersensitivity, accuracy, or reliability.1365.2.EnhancedNociceptionMonitoringwithSensorFusion0102030405060708090100Dental DamBaseline0102030405060708090100Anesthetic BolusResponse0102030405060708090100SNSBaseline0102030405060708090100Change RRBaseline0102030405060708090100Change RRResponse0102030405060708090100Change PPawBaseline0102030405060708090100Change PPawResponse0102030405060708090100OverallResponseFigure 5.1: Boxplots of CRC during the expected “no nociception” periods of each experiment from Chapter 3.The central red bar represents the median (second quartile). The box edges are the first and third quartiles. Thewhiskers extend 1.5x the interquartile range (IQR) beyond the box edges. Points beyond the whiskers are drawnas red plus symbols (+), and may be considered outliers. Such possible outliers were not excluded from analysis.1375.2. Enhanced Nociception Monitoring with Sensor Fusion5.2.1 Combining HRV MeasuresA combination of multiple HRV measurement techniques may be employedto improve nociception monitoring. Combining measures may provide fur-ther insight into nociception and the ANS. It may also be an e↵ective so-lution for monitoring nociception where CRC is dicult or impossible tomeasure. For example, it is dicult to reliably measure CRC during spon-taneous ventilation (see Section 5.3.1). Worse, it is impossible to measureCRC when there is no respiration signal, such as at the beginning of surgeryprior to intubation. Of all the traditional HRV analysis techniques, SDNNappears to be the most reliable and useful measure. However, nonlinear dy-namics measures hold promise for providing a deeper understanding of HRVand the autonomic balance. While many of these techniques are too com-putationally complex for real-time operation (e.g. FTLE or ApEn), futureadvances in processing power may make them more practical. In time, theyshould be re-examined. They may allow for monitoring nociception whenthe link between HR and respiration is nonlinear, such as during sponta-neous ventilation (see Section 5.3.1).5.2.2 Cardiovascular Respiratory CoherenceThe principle of CRC may be extended to incorporate signals from thevasculature.The respiration a↵ects the distribution of blood volume within the body.During spontaneous inhalation, the diaphragm tightens and drops, expand-ing the chest cavity and thus inducing negative pressure within the thorax.The lungs expand and air rushes in down the pressure gradient. Duringspontaneous exhalation, the diaphragm relaxes and the chest cavity col-lapses, inducing positive pressure within the thorax. The lungs contractand air rushes out down the pressure gradient. The respiratory cycle thusinduces a pressure oscillation within the thorax. In mechanical ventilation,the thoracic pressure oscillation is somewhat di↵erent. During inspiration,the ventilator forces air into the lungs with positive pressure. During expi-ration, the ventilator returns to normal atmospheric pressure, allowing thebody to force out the air as usual. Mechanical ventilation thus induces en-tirely positive pressure within the thorax (strong positive pressure duringinhalation and weaker positive pressure during exhalation). However, bothtypes of respiration still induce a pressure oscillation within the thorax. Thisoscillation a↵ects the blood distribution. Changes in pressure expand andcollapse the vasculature in the thorax, causing blood to rush in and out.1385.2. Enhanced Nociception Monitoring with Sensor FusionBlood volume oscillations can be measured using photoplethysmography.In photoplethysmography, light is passed through a perfused (vascularized)tissue, and the transmitted or reflected light is detected. The resulting signalis called a PPG. Typically one or more light emitting diodes (LEDs) are usedas light sources, and a photodiode is used as the detector. Some of the lightis absorbed by the tissue, and some is transmitted or reflected. The tissueis mostly static, absorbing a fixed percentage of the incident light (e.g. skin,bone, fat, and muscle). In contrast, the blood volume is dynamic, producingchanges in the relative absorption and transmission/reflection of light. Theblood volume oscillates with the pulse, increasing at systole (high arterialpressure) and decreasing at diastole (low arterial pressure). When there ismore arterial blood within the tissue (at systole), more light is absorbed andless light reaches the detector, and the PPG signal decreases. The oppositee↵ect occurs at diastole. For example, see the rapid oscillations (havingperiod ⇠ 1 s) in Fig. 5.2. The venous blood is considered more or lessstatic, since it is not pulsatile like the arterial blood. However, the venousblood does oscillate with respiration, as it is pulled into and pushed outof the thorax. The arterial blood is less a↵ected by respiratory oscillations,since it is held at a higher pressure overall. The result is an oscillation in therelative absorption and transmission/reflection of light with the respirationcycle. For example, see the slow oscillations (having period ⇠ 10 - 15 s) inFig. 5.2.Oscillations in the PPG reflect the sympathetic tone. The sympathetictone controls vasoconstriction and vasodilation, changing the perfusion ofthe tissues. When the sympathetic tone is high, the vasculature constricts,peripheral resistance increases, and perfusion decreases. With less blood inthe tissue, more light is transmitted and the PPG signal increases. Thisprinciple is used in the Surgical Stress Index [70]. Furthermore, an in-crease in peripheral resistance decreases the blood oscillations observed inthe PPG. The pulse amplitude (from systole to diastole) decreases, as doesthe respiration cycle amplitude. The e↵ect of respiration on the PPG signalis thus weaker during more sympathetic states, and stronger during moreparasympathetic states.The coherence function may be extended to incorporate the PPG signal.This extended function may be termed cardiovascular respiratory coherence(CVRC). CVRC would be a direct sensor fusion between the HR, respiration,and PPG signals. Further research would be required to determine whetheror not CVRC would o↵er greater insight into the autonomic nervous systemor serve as a better monitor of nociception. However, it may be worth inves-tigating in the future. One important caveat: CVRC would require a raw,1395.2. Enhanced Nociception Monitoring with Sensor Fusion0 5 10 15 20 25 30202530354045Time (s)Signal Amplitude (arbitrary units)  Red signal (~660 nm)Infrared signal (~880 nm)Figure 5.2: Example PPG signal for pulse oximetry. Two light wavelengths(red light at 660 nm and infrared light at 880 nm) are used. Rapid os-cillations represent arterial blood pulsations caused by heart beats. Slowoscillations (baseline wander) are caused by respiration.1405.2. Enhanced Nociception Monitoring with Sensor Fusionunfiltered PPG signal for processing, which may be dicult to obtain. Manypulse oximeter systems perform heavy filtering on the PPG. The signals aretypically band-pass filtered in an attempt to isolate the cardiac pulsations.Other oximeter manufacturers, most notably Masimo, filter the PPG sig-nals with adaptive filters for motion artifact tolerance. Unfortunately, thisfiltering eliminates the respiratory oscillations required by CVRC.5.2.3 End-Tidal CO2End-tidal CO2 (EtCO2) reflects autonomic tone. EtCO2 is the concentra-tion of CO2 measured at the end of exhalation. It is typically the highestmeasured concentration of CO2 in the respiration cycle. CO2 is a wasteproduct of human metabolism. The body’s cells consume O2 to produceadenosine triphosphate (ATP, a source of cellular energy), expelling CO2as a waste product. As a cell’s metabolism increases, it di↵uses more CO2into the bloodstream to be exhaled. When the body’s overall metabolismincreases, a measurable increase in EtCO2 results. The autonomic tone con-tributes to the body’s overall metabolism. In a more sympathetic state, theadrenal glands secrete higher concentrations of epinephrine, norepinephrine,and cortisol, which increase overall cellular metabolism [24, 49]. The oppo-site e↵ect occurs when the body is in a more parasympathetic state. Thus,a change in EtCO2 may reflect a change in autonomic tone. Indeed, changesin EtCO2 appear to have accompanied periods of strong nociception in ourexperimental data. An example of this e↵ect is shown in Fig. 5.3.However, it would be dicult to use EtCO2 alone as a nociception mon-itor. There are many confounding factors that may contribute to changesin EtCO2. Decreasing the respiration rate, airway pressure, or minute ven-tilation will all cause an increase in EtCO2, because of decreased alveolarventilation. That is, the rate of clearance of EtCO2 from the plasma is sim-ply reduced, while the rate of EtCO2 production is unchanged. Inceasingthese parameters (and with them increasing alveolar ventilation) will like-wise cause a decrease in EtCO2. Even changes in the ventilator’s I:E ratiomay a↵ect the EtCO2, since this ratio can a↵ect the alveolar ventilation.Notice, for example, the rapid decrease in EtCO2 that accompanied an in-crease in RR in Fig. 3.21. In contrast, notice the rapid increase in EtCO2that accompanied a decrease in PPaw in Fig. 3.26. These examples are adirect consequence of increased and decreased alveolar ventilation, respec-tively, and are unrelated to any change in nociception or autonomic balance.Furthermore, increased muscular activity (shivering), malignant hyperther-mia, bicarbonate infusion, or tourniquet release will all cause an increase1415.2. Enhanced Nociception Monitoring with Sensor Fusion50 100 150 200 250 300 350 400 4508090100110Heart Rate(bpm)50 100 150 200 250 300 350 400 4500246Respiration(% CO2)50 100 150 200 250 300 350 400 45055.25.45.65.86EtCO2 (%)50 100 150 200 250 300 350 400 450020406080100CRCNociception IndexTime (s)240 s  start of strong nociception (dental dam insertion)Figure 5.3: Example response to a dental dam insertion event, analyzed withreal-time CRC, and showing the time-varying EtCO2. The plot of EtCO2is second from the bottom. Vertical lines denote clinical events. The levelof nociception increases following dental dam insertion, as reflected in theincrease in both the CRC nociception index and in the EtCO2.1425.2. Enhanced Nociception Monitoring with Sensor Fusionin EtCO2, yet these are unrelated to nociception. In contrast, decreasedmuscular activity (e.g. administration of muscle relaxants), hypothermia,pulmonary embolism, or bronchospasm will all cause a decrease in EtCO2.Furthermore, changes in EtCO2 observed in our nociception dataset wereoften subtle, and may be dicult to detect in general. EtCO2 may still serveas a useful adjunct to HRV.HRV and EtCO2 may be indirectly fused to arrive at a better measure ofnociception. In one possible application, EtCO2 may serve as a mechanismfor confirming or refuting a nociception response observed in the HRV. Ifan HRV analysis reports a response to nociception, and if the ventilatoryparameters (e.g. RR, PPaw, MVexp, I:E) have all been kept constant inrecent history, then the EtCO2 should increase shortly afterward. If anincrease in the EtCO2 is indeed observed, then it would serve as confirmationthat the patient is responding to a nociceptive stimulus. In contrast, a failureto detect an increase in EtCO2 would suggest that the HRV analysis may bein error, perhaps due to artifact in the HR source signal. In another possibleapplication, the HRV is used to rapidly detect a response to nociception, andthe subsequent change in EtCO2 is used to measure the magnitude of thechange in autonomic balance. For example, a 25% increase in EtCO2 mayreflect a 25% increase in metabolism, which may correlate to the magnitudeof the change in autonomic balance. However, further research would berequired to understand the exact relationship between EtCO2, metabolismlevel, and autonomic balance.5.2.4 Respiration RateThe RR can reflect the autonomic balance in spontaneously ventilated pa-tients. An increase in sympathetic tone will increase overall metabolism, asdiscussed in the previous section. This metabolic shift will decrease plasmaconcentrations of O2 and increase plasma concentrations of CO2 and H+,thus decreasing the plasma pH (increasing acidity). The body will detectthe changing concentrations of O2, CO2, and H+, and respond by increasingthe RR. However, this e↵ect does not occur during mechanical ventilation.Mechanically ventilated patients are typically receiving drugs that suppressspontaneous ventilation. The ventilator controls the inspiration and expira-tion, and is driven independently of the patient’s metabolism.As with EtCO2, the RR is likely best fused indirectly with other param-eters in a nociception monitor. RR shares many of the same confoundingfactors with EtCO2; indeed, these two measures are intimately linked. TheRR could serve as a valuable confirmation parameter, as with EtCO2. If1435.3. Further Validationthe patient is breathing spontaneously and an HRV analysis reports a re-sponse to nociception, then the RR should increase shortly afterward. If anincrease in RR is indeed observed, then it would serve as confirmation thatthe patient is responding to a nociceptive stimulus.5.3 Further ValidationThe experiments presented in this dissertation serve as a good starting point,but further validation is required.5.3.1 Spontaneous VentilationCRC must be validated on spontaneously ventilated patients. Initial exper-imentation suggests that CRC may not perform as well under these condi-tions. In one experiment, cardiorespiratory data were collected from healthyadult volunteers undergoing a cold pressor test. Subjects’ non-dominanthand was placed in ice water to invoke a sympathetic response to nocicep-tion. CRC was calculated before and during the cold pressor test, followingthe standard Baseline-Stimulus-Response experiment template described inSection 3.4. Sometimes CRC increased as expected (see Fig. 5.4). Othertimes it actually decreased (see Fig. 5.5). Often the CRC nociception indexwas erratic and high throughout the experiment (see Fig. 5.6). Overall, CRCdid not show a statistically significant change to the stimulus, respondingmore or less at random. There are several potential explanations for CRC’sfailure in this experiment.CRC may have diculty tracking the respiration rate. It is dicult toaccurately estimate the instantaneous RR during spontaneous ventilation.The RR trend recorded from the Datex/Ohmeda monitor adapts much tooslowly to changing respiration conditions to accurately track the RR. Themonitor appears to use a 60 s moving average filter, while the RR oftenchanges much more quickly than this during spontaneous ventilation. Forexample, notice the erratic respiration in the CO2 signal in Figs. 5.4 and 5.6.In an attempt to circumvent this problem, a custom-written RR estimationalgorithm was developed. The algorithm calculates the FFT over a short(15 s) window, and estimates the RR as the frequency with the peak power.While not very robust, the algorithm does work in most cases and is asimple first attempt at estimating very short term RR values. However, dueto the Heisenberg-Gabor uncertainty principle, reducing the window lengthin the time domain necessarily reduces precision in the frequency domain.The short-term RR estimation algorithm may adapt quickly to changes in1445.3. Further Validation50 100 150 200 250 300 350 400 4505060708090100110Heart Rate(bpm)50 100 150 200 250 300 350 400 450051015202530Respiration (% CO 2)RR (breaths/min)  RespirationRR50 100 150 200 250 300 350 400 450020406080100CRCNociception IndexTime (s)240 s  start cold pressor (hand in ice water)Figure 5.4: Example response to cold pressor event, analyzed with real-timeCRC. Vertical lines denote clinical events. The green and red vertical bandsrepresent the baseline and response periods, respectively. The nociceptionindex increases strongly (from 50.7 to 78.3) after the subject places his/herhand in the ice water. The missing CRC data at the beginning of the windowcorresponds to the combined length of the analyzing and smoothing filters.Missing data in the middle of the plot corresponds to a temporary cessationof respiration, when the RR could not be accurately calculated.1455.3. Further Validation50 100 150 200 250 300 350 400 4505060708090Heart Rate(bpm)50 100 150 200 250 300 350 400 4500510152025Respiration (% CO 2)RR (breaths/min)  RespirationRR50 100 150 200 250 300 350 400 450020406080100CRCNociception IndexTime (s)240 s  start cold pressor (hand in ice water)Figure 5.5: Example response to cold pressor event, analyzed with real-timeCRC. Vertical lines denote clinical events. The green and red vertical bandsrepresent the baseline and response periods, respectively. The nociceptionindex unexpectedly decreases (from 78.3 to 52.8) after the subject placeshis/her hand in the ice water. The missing CRC data at the beginningof the window corresponds to the combined length of the analyzing andsmoothing filters. Missing data in the middle of the plot corresponds to atemporary cessation of respiration, when the RR could not be accuratelycalculated.1465.3. Further Validation50 100 150 200 250 300 350 400 450304050607080Heart Rate(bpm)50 100 150 200 250 300 350 400 45005101520Respiration (% CO 2)RR (breaths/min)  RespirationRR50 100 150 200 250 300 350 400 450020406080100CRCNociception IndexTime (s)240 s  start cold pressor (hand in ice water)Figure 5.6: Example response to cold pressor event, analyzed with real-timeCRC. Vertical lines denote clinical events. The green and red vertical bandsrepresent the baseline and response periods, respectively. The nociceptionindex increases slightly (from 83.9 to 87.3) after the subject places his/herhand in the ice water; however, the index is very high throughout. The miss-ing CRC data at the beginning of the window corresponds to the combinedlength of the analyzing and smoothing filters. Missing data in the middle ofthe plot corresponds to a temporary cessation of respiration, when the RRcould not be accurately calculated.1475.3. Further Validationrespiration, but it can provide only approximate RR values. Unfortunately,this reduced RR precision translates to CRC filters that may be looking forRSA power at the wrong frequency.Even short-term RR estimates may not be valid over the full durationof the CRC analysis window. The respiration may be nonstationary evenover very short CRC analysis windows. Recall that the Fourier transformassumes the analyzed signal is stationary for all time (infinite time domainwindow). Joint time/frequency analysis techniques like CRC relieve thisassumption by shrinking the analysis window in the time domain, tradingfor reduced localization in the frequency domain. However, these techniquesstill assume the analyzed signal is stationary over their support width in thetime domain. While this assumption is generally valid during mechanicalventilation, spontaneous ventilation can be so erratic as to be highly nonsta-tionary even over very short analysis windows. Therefore, as with the RRtracking problems described above, the CRC filters may again be lookingfor RSA power at the wrong frequency.Furthermore, rapid changes in respiration may cause nonlinear responsesin HRV, which CRC cannot measure. Changes in RR may be too fastfor the HRV can adapt to them. The sympathetic and parasympatheticfeedback loops and the response of the heart to changing autonomic tone allcontribute a delay in the system, limiting its ability to track rapid changesin respiration. The output of the system is nonlinear while it tries to adaptto the rapid changes in input. However, CRC is a linear method. Coherenceestimation is based on a linear time-invariant (LTI) system model, withnonlinear e↵ects modeled as uncorrelated noise. When the HRV response tochanging respiration becomes nonlinear, the apparent coherence decreasesand uncorrelated noise increases. The result is a false increase in the CRCnociception index. Indeed, Keissar et al. discovered that a rapid change inposture from supine to standing produced a strong nonlinear response in theANS, which manifested as a loss of coherence. Likewise, a deep sigh (onetype of rapid change in respiration) also caused a loss of coherence [34].Finally, in certain cases the nociceptive stimulus appears to have eliciteda parasympathetic response, as opposed to the expected sympathetic re-sponse. In Fig. 5.5, for example, the mean HR briefly increased after thehand was inserted into the ice water, but quickly returned to baseline. TheHRV and respiration both showed signs of a shift towards a more parasym-pathetic state. The RSA power increased, and the RR decreased. The CRCnociception index reflected this change by decreasing, and it appears to havebeen the correct response.1485.3. Further Validation5.3.2 In Adults and in Disease StatesCRC has been shown to work on healthy children, but further studies willbe required to confirm these results in adult subjects, and in subjects withdiseases a↵ecting autonomic function (e.g. diabetes mellitus) or medications(e.g. atropine). We have recently shown that CRC responds to circadianchanges in autonomic balance in nine healthy adult subjects aged 26.3 ± 4.6years [8]. The response of CRC to these autonomic changes provides evi-dence that it may be used to monitor nociception in healthy adult subjects.Indeed, it also suggests that CRC can work during spontaneous ventilation,though perhaps only averaged over longer time periods. Controlled nocicep-tion studies in adult subjects should be conducted to validate the results inthis population.The baseline nociception index (i.e. during periods where there is nonociception) can be measured and compared across age stratifications. Pa-tients of di↵erent age groups undergoing routine surgery can be enrolled ina study. For example, the study might include four strata: children, youngadults, middle-aged, and elderly. Following induction of general anesthesiaand intubation, but before the beginning of surgery, the nociception indexcan be measured and recorded. The nociception index can be comparedacross the population strata to determine if there is a change in baselineindex with age.Similarly, the baseline nociception index can be measured and comparedacross patients with and without disease in the same age range. Patients inthe same age group undergoing routine surgery can be enrolled in a study.A control group will have no conditions a↵ecting autonomic or cardiores-piratory function. A second group will have a specific disease or conditionthat is known to a↵ect autonomic or cardiorespiratory function. Followinginduction of general anesthesia and intubation, but before the beginningof surgery, the nociception index can be measured and recorded. The no-ciception index can be compared between the test and control groups todetermine if there is a change in baseline index with disease.5.3.3 Sensitivity to Nociception at Di↵erent Levels ofAntinociceptionA study of patient movement during abdominal incision could be a good testof a nociception monitor. This type of study has been used in the past foridentifying the variables useful for measuring nociception [64]. The studywould require a healthy patient population undergoing abdominal surgery.1495.3. Further ValidationSubjects should not have any diseases a↵ecting autonomic or cardiorespi-ratory function. Patients would receive general anesthesia prior to surgery.Experimental equipment would be set up during the period between induc-tion of anesthesia and the initial incision in the abdomen. The patient wouldbe connected to the nociception monitor and would be visually observed byresearch sta↵ during initial abdominal incision. The abdominal incision isexpected to be a very strong nociceptive stimulus. The nociception indexwould be compared during a period immediately before, during, and imme-diately after the incision. Research sta↵ would identify patient movementduring the incision. The population can be classified post hoc into “movers”and “non-movers,” and the nociception index values compared between thetwo groups. We expect that the “movers” group should have a higher noci-ception index on average.A more advanced study could test the relationship between antinocicep-tive drugs and measured nociception. The study would involve a furtherpopulation sub-division into high- and low-dose antinociception. The pa-tient population would be randomized into two groups. One group receivesa standard infusion of remifentanil, and the other receives a reduced infu-sion. The study is then performed the same as described above. Duringpost hoc analysis, the nociception index values are also compared betweenthe high- and low-dose antinociception groups. We expect to see fewer pa-tients moving in the high-dose group as compared to the low-dose group. Weexpect that the nociception index should be consistently lower during andafter abdominal incision in the high-dose group as compared to the low-dosegroup. Finally, we expect no significant di↵erence in the nociception indexbetween the groups before abdominal incision. The index might be slightlylower in the high-dose group, as it would suppress the baseline sympathetictone. However, in the absence of nociception the antinociceptive dose shouldnot have a significant e↵ect on the index.5.3.4 Respiration Rate Ramp DownWe showed in this dissertation that the CRC nociception index is not signif-icantly a↵ected by a change in RR (see Section 3.8). However, this experi-ment was limited to testing a change from 8 breaths/min to 16 breaths/min.It is unknown how cardiorespiratory coupling behaves outside this range.A respiration rate ramp down study would resolve this uncertainty.The study would require a healthy patient population undergoing routinesurgery. Subjects should not have any diseases a↵ecting autonomic or car-diorespiratory function. Patients would receive general anesthesia prior to1505.3. Further Validationsurgery. The experiment would be conducted during the period betweeninduction of anesthesia and the start of surgery. The ventilator would beconfigured with a set peak airway pressure and inspiration/expiration ra-tio. It would then begin with a fast respiration rate, on the order of 30breaths/min. The respiration rate would be ramped down slowly, in incre-ments of e.g. 2 breaths/min. The ventilator would pause for approximatelyone minute at each step, to ensure the patient’s cardiorespiratory systemhas time to adapt to the change and reach steady state. The study wouldend after the respiration rate reaches a very low value (e.g. 2 breaths/min).As the respiration rate decreases, the minute ventilation will also decrease.Reduced ventilation can be harmful for the patient. The patient’s end-tidalCO2 must be observed to ensure it is in a safe range. The anesthesiologistmay adjust the inspired O2 to compensate for the reduced minute venti-lation. Unfortunately, s/he must not adjust other parameters that woulda↵ect the CO2 (e.g. peak airway pressure or tidal volume) as these maya↵ect a nociception measure based on respiration. The study may need tobe stopped if the CO2 exceeds safe limits. The complete ramp down shouldtake no more than 15 minutes, so it should be fairly unobtrusive to theclinical sta↵.5.3.5 Outcomes StudiesThe ultimate test of a nociception monitor is its e↵ect on end-point clinicaloutcomes. Recall from earlier discussion in Section 1.1 that underdosing andoverdosing of antinociceptive drugs are both dangerous. An underdose leadsto inadequately suppressed nociception, which results in an excessive and/orprolonged stress response, and possible physiological damage. Some patientpopulations are particularly vulnerable to such damage, including the veryelderly and those with compromised cardiovascular systems. Studies havelinked surgically-induced stress hormones with perioperative myocardial in-farction in this population [21, 39]. In contrast, an overdose of antinocicep-tive drugs can lead to post-operative nausea and vomiting (PONV). PONVis one of the most common complications of surgery, leading to delayed re-covery, increased hospital stays, and additional nursing time [30]. PONVcan be a significant contributor to the cost of health care delivery. Reducingthe dose of anesthetic agents can reduce the risk of PONV; however, opti-mal anesthesia must still be the highest priority [28]. Optimal anesthesiamay also reduce the risk of postoperative delirium (POD) and postoperativecognitive dysfunction (POCD). POD has been associated with delayed re-covery and increased hospital stays, while POCD has been associated with1515.3. Further Validationincreased mortality and long-term cognitive disability [22]. Both the surgi-cal stress response and anesthetic drugs have been found to contribute toPOD [32]. A recent study has shown that BIS-guided control of hypnosismay decrease POD and POCD [17]. Similar guidance for antinociceptionmay decrease them further.A study may be designed to track the di↵erence in clinical outcomesresulting from the use of a nociception monitor during general anesthesia. Apatient population can be recruited and divided into two groups. One groupwould be monitored by the nociception monitor during general anesthesia.The other group would serve as a control, receiving standard (traditional)monitoring only. The type of surgery should be kept consistent among allsubjects. Clinical outcomes can then be compared between the two groups.We expect to see a reduction in PONV, POD, and POCD in the nociceptionmonitoring group as compared to the control group. We also expect to see areduction in recovery time, post-operative hospital stay. Finally, we expectto see a reduction in perioperative myocardial infarction. This study mayrequire a large population to show a statistically significant di↵erence inclinical outcomes.152Bibliography[1] A Absalom, D Amutike, A Lal, M White, and G N C Kenny. Accuracyof the ’Paedfusor’ in children undergoing cardiac surgery or catheteri-zation. British Journal of Anaesthesia, 91(4):507–513, October 2003.[2] S Akselrod, D Gordon, F Ubel, D Shannon, A Berger, and R Cohen.Power spectrum analysis of heart rate fluctuation: a quantitative probeof beat-to-beat cardiovascular control. Science, 213(4504):220–222, July1981.[3] R Bailon, P Laguna, L Mainardi, and L Sornmo. Analysis of heartrate variability using time-varying frequency bands based on respiratoryfrequency. Conf Proc IEEE Eng Med Biol Soc, pages 6675–8, January2007.[4] R Barbieri, E C Matten, A A Alabi, and E N Brown. A point-processmodel of human heartbeat intervals: new definitions of heart rate andheart rate variability. Am J Physiol Heart Circ Physiol, 288:424–435,2005.[5] N Y Belova, S V Mihaylov, and B G Piryova. Wavelet transform: Abetter approach for the evaluation of instantaneous changes in heartrate variability. Autonomic Neuroscience: Basic & Clinical, 131(1-2):107–22, January 2007.[6] R D Berger, S Akselrod, D Gordon, and R J Cohen. An EcientAlgorithm for Spectral Analysis of Heart Rate Variability. IEEE TransBiomed Eng, 33(9):900–4, 1986.[7] G G Berntson, J T Bigger, D L Eckberg, P Grossman, P G Kaufmann,M Malik, H N Nagaraja, S W Porges, J P Saul, P H Stone, and M W vanDer Molen. Heart rate variability: origins, methods, and interpretivecaveats. Psychophysiology, 34(6):623–48, November 1997.153Bibliography[8] P Boudreau, C J Brouse, G A Dumont, and D B Boivin. Sleep-stateand circadian variation of cardiorespiratory coherence. In Conf ProcIEEE Eng Med Biol Soc, page in press, 2012.[9] C J Brouse, G Dumont, P Yang, J Lim, and J M Ansermino. iAssist:a software framework for intelligent patient monitoring. In Conferenceproceedings: IEEE Engineering in Medicine and Biology Society, vol-ume 2007, pages 3790–3, January 2007.[10] C J Brouse, G A Dumont, F J Herrmann, and J M Ansermino. AWavelet Approach to Detecting Electrocautery Noise in the ECG. IEEEEngineering in Medicine and Biology Magazine, pages 76–82, 2006.[11] C J Brouse, G A Dumont, D Myers, E Cooke, and J M Ansermino.Wavelet Transform Cardiorespiratory Coherence for Monitoring Noci-ception. In Computing in Cardiology, pages 713–6, 2010.[12] C J Brouse, G A Dumont, D Myers, E Cooke, J Stinson, J Lim, andJ M Ansermino. Real-Time Cardiorespiratory Coherence Detects No-ciception During General Anesthesia. IEEE Trans Biomed Eng, page(submitted), 2013.[13] C J Brouse, W Karlen, G A Dumont, D Myers, E Cooke, J Stin-son, J Lim, and J M Ansermino. Real-Time Cardiorespiratory Co-herence Detects Antinociception During General Anesthesia. In ConfProc IEEE Eng Med Biol Soc, pages 3813–3816, 2012.[14] C J Brouse, W Karlen, G A Dumont, D Myers, E Cooke, J Stinson,J Lim, and J M Ansermino. Monitoring Nociception During GeneralAnesthesia with Cardiorespiratory Coherence. J Clin Monit Comput,page (accepted for publication), 2013.[15] C J Brouse, W Karlen, D Myers, E Cooke, J Stinson, J Lim, G ADumont, and J M Ansermino. Wavelet transform cardiorespiratorycoherence detects patient movement during general anesthesia. In ConfProc IEEE Eng Med Biol Soc, volume 33, pages 6114–7, 2011.[16] R A Brown and R Frayne. A Fast Discrete S-Transform for BiomedicalSignal Processing. Conf Proc IEEE Eng Med Biol Soc, 2008:2586–9,January 2008.[17] M T V Chan, B C P Cheng, T M C Lee, and T Gin. BIS-guided anes-thesia decreases postoperative delirium and cognitive decline. Journalof neurosurgical anesthesiology, 25(1):33–42, January 2013.154Bibliography[18] E K Choo, W Magruder, C J Montgomery, J Lim, R Brant, and J MAnsermino. Skin conductance fluctuations correlate poorly with post-operative self-report pain measures in school-aged children. Anesthesi-ology, 113(1):175–82, July 2010.[19] F Christiansen and H H Rugh. Computing Lyapunov spectra withcontinuous GramSchmidt orthonormalization. Nonlinearity, 10:1063–1072, 1997.[20] L Cohen. Time-frequency distributions: a review. Proceedings of theIEEE, 77(7):941–81, 1989.[21] M M Dawood, D K Gutpa, J Southern, A Walia, J B Atkinson, and K AEagle. Pathology of fatal perioperative myocardial infarction: implica-tions regarding pathophysiology and prevention. International Journalof Cardiology, 57:37–44, November 1996.[22] S Deiner and J H Silverstein. Postoperative delirium and cognitive dys-function. British journal of anaesthesia, 103 Suppl:i41–46, December2009.[23] T A Denton, G A Diamond, R H Helfant, S Khan, and H Karagueuzian.Fascinating rhythm: a primer on chaos theory and its application tocardiology. Am Heart J, 120:1419–1440, 1990.[24] J P Desborough. The stress response to trauma and surgery. Brit JAnaesth, 85(1):109–17, July 2000.[25] S Elsenbruch, Z Wang, W C Orr, and J D Z Chen. Time-frequencyanalysis of heart rate variability using short-time Fourier analysis. Sig-nal Processing, 21:229–240, 2000.[26] F M Fouad, R C Tarazi, C M Ferrario, S Fighaly, and C Alicandri.Assessment of parasympathetic control of heart rate by a noninvasivemethod. American Journal of Physiology - Heart and Circulatory Phys-iology, 246:838–42, 1984.[27] D Gabor. Theory of Communication. Journal of the Institution ofElectrical Engineers, 93(26):429–441, 1946.[28] T J Gan. Postoperative Nausea and Vomiting - Can It Be Eliminated?Jama, 287(10):1233–6, March 2002.155Bibliography[29] A Grinsted, J C Moore, and S Jevrejeva. Application of the crosswavelet transform and wavelet coherence to geophysical time series.Nonlinear Proc Geoph, 11:561–566, 2004.[30] J Hirsch. Impact of postoperative nausea and vomiting in the surgicalsetting. Anaesthesia, 49 Suppl:30–3, January 1994.[31] P G Katona and F Jih. Respiratory sinus arrhythmia: noninvasive mea-sure of parasympathetic cardiac control. Journal of Applied Physiology,39(5):801–805, 1975.[32] H Kehlet. Multimodal approach to control postoperative pathophysi-ology and rehabilitation. Brit J Anaesth, 78:606–17, 1997.[33] K Keissar, L R Davrath, and S Akselrod. Time Frequency WaveletTransform Coherence of Cardio-Respiratory Signals during Exercise. InComputers in Cardiology, pages 733–6, 2006.[34] K Keissar, L R Davrath, and S Akselrod. Coherence analysis betweenrespiration and heart rate variability using continuous wavelet trans-form. Phil Trans R Soc A, 367:1393–406, April 2009.[35] R E Kleiger, J P Miller, J T Bigger, and A J Moss. Decreased heartrate variability and its association with increased mortality after acutemyocardial infarction. The American Journal of Cardiology, 59:256–62,February 1987.[36] T Laitinen, L Niskanen, G Geelen, and E La. Age dependency ofcardiovascular autonomic responses to head-up tilt in healthy subjects.J Appl Physiol, 96:2333–2340, 2004.[37] T Laitio, J Jalonen, T Kuusela, and H Scheinin. The role of heartrate variability in risk stratification for adverse postoperative cardiacevents. Anesthesia & Analgesia, 105(6):1548–60, December 2007.[38] T T Laitio, H V Huikuri, T H Makikallio, J Jalonen, E S H Kentala,H Helenius, O Pullisaar, J Hartiala, and H Scheinin. The Breakdownof Fractal Heart Rate Dynamics Predicts Prolonged Postoperative My-ocardial Ischemia. Anesthesia & Analgesia, 98:1239–44, 2004.[39] G Landesberg, W S Beattie, M Mosseri, A S Ja↵e, and J S Alpert.Perioperative myocardial infarction. Circulation, 119:2936–44, June2009.156Bibliography[40] E B Liem, T V Joiner, K Tsueda, and D I Sessler. Increased Sensitiv-ity to Thermal Pain and Reduced Subcutaneous Lidocaine Ecacy inRedheads. Anesthesiology, 102(3):509–514, 2005.[41] E B Liem, C-M Lin, M-I Suleman, A G Doufas, R G Gregg, J M Veau-thier, G Loyd, and D I Sessler. Anesthetic Requirement is Increased inRedheads. Anesthesiology, 101(2):279–283, 2004.[42] D T J Liley, N C Sinclair, T Lipping, B Heyse, H E M Vereecke, andM M R F Struys. Propofol and remifentanil di↵erentially modulatefrontal electroencephalographic activity. Anesthesiology, 113(2):292–304, August 2010.[43] R Logier, M Jeanne, A Dassonneville, M Delecroix, and B Tavernier.PhysioDoloris: a monitoring device for Analgesia / Nociception balanceevaluation using Heart Rate Variability analysis. In Conf Proc IEEEEng Med Biol Soc, pages 1194–7, 2010.[44] N R Lomb. Least-Squares Frequency Analysis of Unequally-SpacedData. Physics, 39(1964), 1976.[45] F Lombardi, A Malliani, M Pagani, and S Cerutti. Heart rate variabilityand its sympatho-vagal modulation. Cardiovascular Research, 32:208–216, 1996.[46] S Mallat. A Wavelet Tour of Signal Processing. Elsevier, London, 2ndedition, 1999.[47] G F Margrave, P C Gibson, J P Grossman, D C Henley, V Iliescu,and M P Lamoureux. The Gabor transform, pseudodi↵erential opera-tors, and seismic deconvolution. Integrated Computer-Aided Engineer-ing, 12:43–55, 2005.[48] G A Mashour, L Y-J Wang, C R Turner, J C Vandervest, A Shanks, andK K Tremper. A retrospective study of intraoperative awareness withmethodological implications. Anesthesia and Analgesia, 108(2):521–6,February 2009.[49] D E Matthews, G Pesola, and R G Campbell. E↵ect of epinephrine onamino acid and energy metabolism in humans. Am J Physiol-Endoc M,258:E948–56, June 1990.157Bibliography[50] M Merri, D C Farden, J G Mottley, and E L Titlebaum. Samplingfrequency of the electrocardiogram for spectral analysis of the heart ratevariability. IEEE Transactions on Biomedical Engineering, 37(1):99–106, January 1990.[51] M J Mertens, E Olofsen, F H M Engbers, A G L Burm, J G Bovill, andJ Vuyk. Propofol Reduces Perioperative Remifentanil Requirements ina Synergistic Manner. Anesthesiology, 99(2):347–359, August 2003.[52] G Moore. Cramming More Components onto Integrated Circuits. Elec-tronics Magazine, 38(8), 1965.[53] M Orini, R Bailo´n, L T Mainardi, P Laguna, and P Flandrin. Charac-terization of Dynamic Interactions Between Cardiovascular Signals byTime-Frequency Coherence. IEEE Transactions on Biomedical Engi-neering, 59(3):663–73, March 2012.[54] T Penzel, J W Kantelhardt, L Grote, J H Peter, and A Bunde. Com-parison of Detrended Fluctuation Analysis and Spectral Analysis forHeart Rate Variability in Sleep and Sleep Apnea. IEEE Transactionson Biomedical Engineering, 50(10):1143–1151, 2003.[55] S M Pincus. Approximate entropy as a measure of system complexity.Mathematische Annalen, 88:2297–2301, 1991.[56] K-K Poh and P Marziliano. Analysis of Neonatal EEG Signals usingStockwell Transform. Conf Proc IEEE Eng Med Biol Soc, 2007(1):594–7, January 2007.[57] R J Pollard, J P Coyle, R L Gilbert, and J E Beck. IntraoperativeAwareness in a Regional Medical System - A Review of 3 Years’ Data.Anesthesiology, 106(2):269–74, December 2007.[58] S Qian and D Chen. Joint Time-Frequency Analysis. IEEE SignalProcessing Magazine, pages 52–67, 1999.[59] M Raczkowska, D L Eckberg, and T J Ebert. Muscarinic cholinergicreceptors modulate vagal cardiac responses in man. Journal of theAutonomic Nervous System, 7:271–8, 1983.[60] M Rantanen, H Yppa¨rila¨-Wolters, M van Gils, A Yli-Hankala, M Huiku,M Kyma¨la¨inen, and I Korhonen. Tetanic stimulus of ulnar nerve as apredictor of heart rate response to skin incision in propofol remifentanilanaesthesia. Br J Anaesth, 99(4):509–13, October 2007.158Bibliography[61] R H Sandin, G Enlund, P Samuelsson, and C Lennmarken. Aware-ness during anaesthesia: a prospective case study. Lancet, 355:707–11,February 2000.[62] J D Scargle. Statistical Aspects of Spectral Analysis of Unevenly SpacedData. The Astrophysical Journal, 263:835–853, 1982.[63] P S Sebel, T A Bowdle, M M Ghoneim, I J Rampil, R E Padilla, T JGan, and K B Domino. The incidence of awareness during anesthesia: amulticenter United States study. Anesthesia and analgesia, 99(3):833–9, September 2004.[64] E R J Seitsonen, I K J Korhonen, M J van Gils, M Huiku, J M PLo¨tjo¨nen, K T Korttila, and A M Yli-Hankala. EEG spectral entropy,heart rate, photoplethysmography and motor responses to skin incisionduring sevoflurane anaesthesia. Acta Anaesth Scand, 49:284–92, 2005.[65] R P Smith, J Argod, J L Pe´pin, and P A Le´vy. Pulse transit time: anappraisal of potential clinical applications. Thorax, 54(5):452–7, May1999.[66] U M Stamer and F Stu¨ber. The pharmacogenetics of analgesia. Expertopinion on pharmacotherapy, 8(14):2235–45, October 2007.[67] R G Stockwell, L Mansinha, and R P Lowe. Localization of the ComplexSpectrum: The S Transform. IEEE Transactions on Signal Processing,44(4):998–1001, 1996.[68] R J Storella, R B Kandell, J C Horrow, T S Ackerman, M Polansky, andS Zietz. Nonlinear Measures of Heart Rate Variability Fentanyl-BasedInduction of Anesthesia. Methods, 81:1292–1294, 1995.[69] H Storm. Changes in skin conductance as a tool to monitor nociceptivestimulation and pain. Current Opinion in Anaesthesiology, 21(6):796–804, 2008.[70] M M R F Struys, C Vanpeteghem, M Huiku, K Uutela, N B K Blyaert,and E P Mortier. Changes in a surgical stress index in response tostandardized pain stimuli during propofolremifentanil infusion. Brit JAnaesth, 99(3):359–67, September 2007.[71] Task Force of the European Society of Cardiology and the North Amer-ican Society of Pacing Electrophysiology. Heart rate variability: stan-159Bibliographydards of measurement, physiological interpretation, and clinical use.European Heart Journal, 17:354–381, 1996.[72] J A Taylor, C W Myers, J R Halliwill, H Seidel, and D L Eckberg.Sympathetic restraint of respiratory sinus arrhythmia: implications forvagal-cardiac tone assessment in humans. Am J Physiol Heart CircPhysiol, 280(6):2804–14, June 2001.[73] J F Thayer, S S Yamamoto, and J F Brosschot. The relationship ofautonomic imbalance, heart rate variability and cardiovascular diseaserisk factors. International Journal of Cardiology, 141:122–31, May 2010.[74] C Torrence and G P Compo. A Practical Guide to Wavelet Analysis.Bulletin of the American Meteorological Society, 79(1):61–78, 1998.[75] D L Toweill, W D Kovarik, R Carr, D Kaplan, S Lai, S Bratton, andB Goldstein. Linear and nonlinear analysis of heart rate variabilityduring propofol anesthesia for short-duration procedures in children.Pediatric Critical Care Medicine, 4(3):308–14, July 2003.[76] T J Ulrych. Maximum entropy power spectrum of truncated sinusoids.Journal of Geophysical Research, 77(8):1396–1400, 1972.[77] T J Ulrych, D E Smylie, O G Jensen, and G K C Clarke. Predictivefiltering and smoothing of short records by using maximum entropy.Journal of Geophysical Research, 78(23):4959–4964, 1973.[78] H G Van Steenis and J H M Tulen. The Exponential Distribution Ap-plied to Nonequidistantly Sampled Cardiovascular Time Series. Com-puters and Biomedical Research, 193(29):174–193, 1996.[79] N Wessel, U Meyerfeldt, A Schirdewan, J Kurths, and A Voss. Short-term Forecasting of Life-threatening Arrhythmias with Finite TimeLyapunov Exponents. Conference proceedings: IEEE Engineering inMedicine and Biology Society, 20(1):326–329, 1998.160

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.24.1-0167194/manifest

Comment

Related Items