Citation
Analysis of neonatal heart rate variability and cardiac orienting responses

Material Information

Title:
Analysis of neonatal heart rate variability and cardiac orienting responses
Creator:
Lebrun, Devin N. ( Dissertant )
van Oostrom, Johannes H. ( Thesis advisor )
Place of Publication:
Gainesville, Fla.
Publisher:
University of Florida
Publication Date:
Copyright Date:
2003
Language:
English

Subjects

Subjects / Keywords:
Datasets ( jstor )
Electrocardiography ( jstor )
Heart rate ( jstor )
Infants ( jstor )
Neonates ( jstor )
Outliers ( jstor )
Quadrants ( jstor )
Signals ( jstor )
Spectroscopy ( jstor )
Standard deviation ( jstor )
Heart rate monitoring ( lcsh )
Dissertations, Academic -- UF -- Electrical and Computer Engineering
Electrical and Computer Engineering thesis, M.E
Electrocardiography -- Data processing ( lcsh )
Signal processing -- Digital techniques ( lcsh )
Genre:
bibliography ( marcgt )
non-fiction ( marcgt )
theses ( marcgt )

Notes

Abstract:
We have applied signal-processing methods to evaluate the clinical significance of heart rate variability (HRV) features, assess developmental changes in HRV, evaluate vagal tone, and demonstrate a cardiac orienting response (COR). The study is based on a sample of 28 low-risk premature newborns randomly assigned to one of two groups. Group 1 is exposed to an auditory stimulation beginning at 28 weeks GA and Group 2 begins exposure at 32 weeks GA. An eight-minute ECG recording is taken once a week from 28-34 weeks GA. The recordings were analyzed in the time and frequency domain using MATLAB. Extracted features are used to classify developmental changes in HRV across individuals and subject groups. In addition, we have used a staggered grouping approach to separate cardiac orienting responses from reflexive responses. At the conclusion of this thesis, the sample size (n=5) was not large enough to elicit any developmental trends in GA. However, we were able to extract useful features to describe developmental changes in HRV and we have acknowledged Empirical Mode Decomposition as an effective measure of a cardiac orienting response. Also, we have identified numerous exogenous interactions and perturbations which can interfere with the baroreceptor and chemoreceptor responses controlling HRV. These interactions were not considered in the original design of the study but will now be considered during final analysis.
Subject:
density, ecg, heart, neonate, power, premature, spectral, stationarity, variability
General Note:
Title from title page of source document.
General Note:
Includes vita.
Thesis:
Thesis (M.E.)--University of Florida, 2003.
Bibliography:
Includes bibliographical references.
General Note:
Text (Electronic thesis) in PDF format.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright Lebrun, Devin N.. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Embargo Date:
9/9/1999
Resource Identifier:
002950135 ( aleph )
53208737 ( OCLC )

Downloads

This item has the following downloads:

PDF ( 3 MBs ) ( .pdf )

lebrun_d_Page_082.txt

lebrun_d_Page_072.txt

lebrun_d_Page_004.txt

lebrun_d_Page_033.txt

lebrun_d_Page_057.txt

lebrun_d_Page_060.txt

lebrun_d_Page_084.txt

lebrun_d_Page_089.txt

lebrun_d_Page_020.txt

lebrun_d_Page_093.txt

lebrun_d_Page_002.txt

lebrun_d_Page_048.txt

lebrun_d_Page_010.txt

lebrun_d_Page_098.txt

lebrun_d_Page_091.txt

lebrun_d_Page_036.txt

lebrun_d_Page_024.txt

lebrun_d_Page_037.txt

lebrun_d_Page_046.txt

lebrun_d_Page_073.txt

lebrun_d_Page_064.txt

lebrun_d_Page_049.txt

lebrun_d_Page_099.txt

lebrun_d_Page_069.txt

lebrun_d_Page_066.txt

lebrun_d_Page_062.txt

lebrun_d_Page_051.txt

lebrun_d_Page_012.txt

lebrun_d_Page_075.txt

lebrun_d_Page_041.txt

lebrun_d_Page_035.txt

lebrun_d_Page_044.txt

lebrun_d_Page_065.txt

lebrun_d_Page_025.txt

lebrun_d_Page_067.txt

lebrun_d_Page_014.txt

lebrun_d_Page_056.txt

lebrun_d_Page_031.txt

lebrun_d_Page_088.txt

lebrun_d_Page_034.txt

lebrun_d_Page_086.txt

lebrun_d_Page_103.txt

lebrun_d_Page_008.txt

lebrun_d_Page_097.txt

lebrun_d_Page_030.txt

lebrun_d_Page_102.txt

lebrun_d_Page_047.txt

lebrun_d_Page_003.txt

lebrun_d_Page_096.txt

lebrun_d_Page_087.txt

lebrun_d_Page_009.txt

lebrun_d_Page_054.txt

lebrun_d_Page_050.txt

lebrun_d_Page_016.txt

lebrun_d_Page_026.txt

lebrun_d_Page_043.txt

lebrun_d_Page_045.txt

lebrun_d_Page_019.txt

lebrun_d_Page_100.txt

lebrun_d_Page_092.txt

lebrun_d_Page_022.txt

lebrun_d_Page_042.txt

lebrun_d_Page_068.txt

lebrun_d_Page_074.txt

lebrun_d_Page_013.txt

lebrun_d_Page_018.txt

lebrun_d_Page_055.txt

lebrun_d_Page_001.txt

lebrun_d_Page_085.txt

lebrun_d_Page_059.txt

lebrun_d_Page_063.txt

lebrun_d_Page_052.txt

lebrun_d_Page_080.txt

lebrun_d_Page_005.txt

lebrun_d_Page_094.txt

lebrun_d_Page_029.txt

lebrun_d_Page_017.txt

lebrun_d_Page_079.txt

lebrun_d_Page_028.txt

lebrun_d_Page_023.txt

lebrun_d_Page_095.txt

lebrun_d_Page_071.txt

lebrun_d_Page_040.txt

lebrun_d_Page_007.txt

lebrun_d_Page_061.txt

lebrun_d_Page_076.txt

lebrun_d_Page_083.txt

lebrun_d_Page_011.txt

lebrun_d_Page_081.txt

lebrun_d_Page_038.txt

lebrun_d_Page_053.txt

lebrun_d_Page_078.txt

lebrun_d_pdf.txt

lebrun_d_Page_104.txt

lebrun_d_Page_021.txt

lebrun_d_Page_006.txt

lebrun_d_Page_039.txt

lebrun_d_Page_032.txt

lebrun_d_Page_015.txt

lebrun_d_Page_070.txt

lebrun_d_Page_090.txt

lebrun_d_Page_058.txt

lebrun_d_Page_101.txt

lebrun_d_Page_027.txt

lebrun_d_Page_077.txt


Full Text












ANALYSIS OF NEONATAL HEART RATE VARIABILITY AND CARDIAC
ORIENTING RESPONSES















By

DEVIN N. LEBRUN


A THESIS PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF ENGINEERING

UNIVERSITY OF FLORIDA


2003

































Copyright 2003

by

Devin N. Lebrun

































Dedicated to my father, Denis Lebrun, for his love and support.















ACKNOWLEDGMENTS

I would like to thank my advisor, Dr. Johannes van Oostrom, for all the advice he

has given me on realizing my career goals. I knew what I wanted, but he helped me get

from there to here and now onto somewhere great. He endured barrages of

scatterbrained, unintelligent questions and always surfaced with the answers I needed.

I am grateful to Dr. Charlene Krueger for allowing me to share in her research and

express my gratitude to Dr. Krueger, Dr. van Oostrom, and the College of Nursing for

providing me with funding throughout my research. I would also like to thank Kelly

Spaulding, John Hardcastle, Dr. JS Gravenstein, Dr. Richard Melker, and Dr. John Harris

for their support in various areas.

Words cannot describe my appreciation for my parents. They have taught me to be

determined, self-motivated, and not afraid of failure. I thank Linnea, for her

companionship and assistance throughout the past 6 years. Finally, thanks go to all my

friends for keeping me from doing something stupid (well, more stupid).
















TABLE OF CONTENTS
Page

A C K N O W L E D G M E N T S ................................................................................................ IV

LIST OF TABLES ...... ........................... ............. ........... ............... VIII

LIST OF FIGURES ......... ....................... ............ ........ ........... IX

ABSTRACT .............. .................... .......... .............. XI

CHAPTER

1 IN T R O D U C T IO N ............................................................................. .............. ...

2 AUTONOMIC NERVOUS SYSTEM AND HEART RATE VARIABILITY ..........4

A utonom ic N ervous System ...........................................................................4
D evelopm ental C changes ................................................................ ...................... 6
V agal Tone and H om eostasis ............................................... ............................ 7
C ardiac O renting R response ............................................................................. 8
Respiratory Sinus Arrhythm ia ............................................................. ........ .....8

3 PRELIMINARY RESULTS AND STUDY PROCEDURES............................... 10

P relim inary F finding s ......................................................... .. ...... .. ...... ............10
Proposed Study ............................................................... .............. ........ 11
Recitation Sessions .................. .......................... .. .. .. .............. ... 12
T e st S e ssio n s ................................................................................................. 1 3

4 T E S T IN G .......................................................................................16

5 QRS DETECTION .............................................................................. 21

N oise/B aselin e W an d er ..............................................................................................2 3
S lo p e D ete ctio n ..................................................................................................... 2 7
O u tlie rs ..............................................................................3 2
D direct RRi Analysis ............................................................................................ 34
Slope Analysis ........................................................................ ......... .......... ........ 36
D differential A naly sis ........................................................................... 37
Interp olation ........... .. ......... ............. .. .. ....................................... 39



v









6 STA TISTIC A L A N A L Y SIS ............................................................. ....................44

T im e B a se d ....................................................................................................4 4
First O rder Statistics ............................................................................ 44
P o in car P lots ...............................................................4 5
Q u adrant A naly sis ............................. .... ................ .. ......... .... ............4 8
Frequency Based............................................ ........ 50
Fast Fourier Spectral A nalysis......................................... .......................... 51
L om b Spectral A naly sis ........................................ .........................................52
Empirical Mode Decomposition...................... ..... .......................... 54
Statistical A analysis of Param eters ........................................ .......................... 57

7 DISCU SSION AND RESULTS .......................................... ............................. 59

E engineering R results .............................................. ... ............ .............. 59
QRS D election Algorithm ...........................................................................59
H RV Plots and Statistics .................................. ......................................60
HR Plot .................. ..... ................... 62
Q u a d ran t P lo t........................................................................................... 6 3
P oin care P lot ............................................................64
F ourier and L om b P lots......... .............................................. ..... .......... 65
EMD Plots ................ ......... .......... ........67
Clinical Results ............. ...... ..... .. ......................................... 68

8 C O N C L U SIO N ......... ......................................................................... ......... ........72

9 FU TU R E W O R K ........................ ............................ .. .... ................ 74

Record Baseline A dult HRV D ata.................................... ........................ 75
Creating the Index of ANS Maturation .............. ....................................76

APPENDIX

A TIME-DOMAIN HRV ANALYSIS .............................. 77

W e e k O n e .......................................................................... 7 7
W e e k T w o ........................................................................ 7 8
W eek T h re e ...................................................................... 7 9
W eek F ou r ............................................................................ 80
W eek Five .................................... .................................. .......... 81
W e e k S ix ......................................................................... 8 2
W eek S ev en ........................................................................... 83

B FREQUENCY-DOMAIN HRV ANALYSIS ................................... ............... 84

W e e k O n e .......................................................................... 8 4
W e e k T w o ........................................................................ 8 5










W e e k T h re e ............................................................................................................ 8 5
W eek Four ................................................................86
W e e k F iv e .............................................................................................................. 8 6
W e e k S ix ................................................................................................................ 8 7
W eek Seven ................................................................87

LIST OF REFEREN CE S ................................................................. .. ......... 88

BIO GR A PH ICA L SK ETCH ........................................................................................ 92
















LIST OF TABLES


Table p

3-1 Organization of R ecitation and Test Sessions ...........................................................11

7-1 Shows the performance of the QRS detection algorithm ...........................................60

7-2 Results of the repeated measures ANOVA with independent variable (SV power)..69

7-3 Results of the repeated measures ANOVA with independent variable (HF power)..70
















LIST OF FIGURES


Figure p

2-1 The divisions of the central nervous system ........................................ ..................4

3-1 Time-line for the windowing of test data. .......... ..... ...........................14

4-1 Example of text file generated from Data Acquisition Program............................18

4-2 The setup for the test sessions. ........................................ ............................ 19

5-1 Q R S detection algorithm ............................ ............. .................................... 23

5-2 Magnitude and phase response of the 60Hz notch filter used to filter the raw ECG
signal. ..................................................................25

5-3 Magnitude and phase response of the bandpass filter used to filter the raw ECG
signal. ..................................................................26

5-4 Magnitude and phase response of the differentiator function.............................27

5-6 Example of detected QRS complexes and beat intervals...................................31

5-7 Histogram of RRi distribution from one subject................................. ...............33

5-8 Standard deviations of the Normal Distribution .......................................... 34

5-9 Illustrates the bouncing effect and upper and lower bound symmetry problems with
using the direct RRi analysis method for tagging outliers. .....................................35

5-10 Illustrates an RRi trend that would not be preserved using the slope analysis
m ethod for tagging outliers. ...... ........................... ......................................37

5-11 Outliers that have been tagged, using the differential method, are marked by a
green dot. Shows the efficiency of the method to preserve trend and remove
bouncing effects. ......................................................................39

5-12 LaGrange Polynomial Interpolation...................... .... .......................... 40

5-13 L near Interpolation ...... ...................................................................... .. .... ... ....4 1

5-14 Cubic Spline Interpolation ............................................... ............................ 41









5-15 C ub ic Interp olation .......................................................................... ................... 42

6-1 Poincare plot representation of RRi. Since data points overlap, colors, coordinated
with the visible spectrum, are used to indicate magnitude at each point ................47

6-2 Defines the quadrants used in Quadrant Analysis.................................. ............... 48

6-3 Example of Quadrant Analysis applied to an RRi dataset .................................49

7-1 HR plots of example RRi data. From top to bottom they represent: the (a)
developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus periods............62

7-2 Quadrant Analysis plots of example RRi data. From top to bottom they represent:
the (a) developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus periods.63

7-3 Poincare plots of example RRi data. From top to bottom they represent: the (a)
developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus periods............64

7-4 Fourier PSD plots of example RRi data. From top to bottom they represent: the (a)
developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus periods............65

7-5 Lomb PSD plots of example RRi data. From top to bottom they represent: the (a)
developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus periods............66

7-6 Hilbert-Haung spectrum of the COR windows of an example RRi dataset.............68

7-7 Plot of mean SV power with standard deviation bars across sessions ...................69

7-8 Plot of mean HF power with standard deviation bars across sessions....... ........ 70















Abstract of Thesis Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Master of Engineering

ANALYSIS OF NEONATAL HEART RATE VARIABILITY AND CARDIAC
ORIENTING RESPONSES


By

Devin N. Lebrun

August 2003

Chair: Johannes H. van Oostrom
Major Department: Electrical and Computer Engineering

We have applied signal-processing methods to evaluate the clinical significance of

heart rate variability (HRV) features, assess developmental changes in HRV, evaluate

vagal tone, and demonstrate a cardiac orienting response (COR).

The study is based on a sample of 28 low-risk premature newborns randomly

assigned to one of two groups. Group 1 is exposed to an auditory stimulation beginning

at 28 weeks GA and Group 2 begins exposure at 32 weeks GA. An eight-minute ECG

recording is taken once a week from 28-34 weeks GA. The recordings were analyzed in

the time and frequency domain using MATLAB. Extracted features are used to classify

developmental changes in HRV across individuals and subject groups. In addition, we

have used a staggered grouping approach to separate cardiac orienting responses from

reflexive responses.

At the conclusion of this thesis, the sample size (n=5) was not large enough to elicit

any developmental trends in GA. However, we were able to extract useful features to









describe developmental changes in HRV and we have acknowledged Empirical Mode

Decomposition as an effective measure of a cardiac orienting response. Also, we have

identified numerous exogenous interactions and perturbations which can interfere with

the baroreceptor and chemoreceptor responses controlling HRV. These interactions were

not considered in the original design of the study but will now be considered during final

analysis.














CHAPTER 1
INTRODUCTION

The premature infant faces many obstacles beyond the control of present day

medical science. Since the loss of the "in utero environment" strips the premature infant

of the most advantageous setting for growth and development, it is increasingly important

for medical science to reach a greater understanding of the developmental changes during

this final trimester. Survival rates of neonates will continue to increase as science begins

to elucidate more regarding this crucial period of development. In general, scientific

advancements are predicated upon a divide and conquer approach. It is important to

embrace the past discoveries as background for future research. Therefore, this study

expands upon a preliminary study conducted by the Principal Investigator (PI), Charlene

Krueger, Ph.D, ARNP. The findings of the preliminary study indicate that the analysis of

heart rate variability (HRV) can provide significant insight into the development of the

premature and fetus' autonomic nervous system (ANS). In collaboration with Dr.

Krueger, the current study moves from the comforts of the uterus to the premature

infant's continuing development in a Neonatal Intensive Care Unit (NICU) setting.

About a third of premature births occur for no apparent reason and with little or no

warning. The list below describes some known causes of premature birth [1,2]:

* Pre-eclampsia, occurring in about 1 in 14 pregnancies, causes around a third of all
premature births.

* Pregnancies involving multiple fetuses are likely to end early.

* Stressful events can elicit the mother into a premature labor.









* An infection involving the uterus may trigger an early delivery.

* The mother's water may break early, starting the delivery process.

* A shortage of oxygen exchange to the placenta may cause a doctor to advise a
caesarian section.

* Malnutrition or chronic diseases, including high blood pressure, diabetes, kidney
disease and hypothyroidism can cause premature births.

Regardless of the cause for an early birth, premature infants are commonly highly

dependant upon the continuation of breathing support, intravenous nutrition, intravenous

antibiotics, and phototherapy. These infants are frequently placed in NICUs, which

specialize in caring for the very young and/or very small. The neonate's high

dependency upon medical interventions demonstrates a need to establish criteria

regarding his or her efficiency at regulating homeostasis. Homeostasis is the

maintenance of equilibrium in a biological system by means of automatic mechanisms,

which counteract influences trending toward disequilibria. It involves constant internal

monitoring and regulating of numerous factors, including oxygen, carbon dioxide,

nutrients, hormones, and organic and inorganic substances. Due to homeostatic controls,

the concentrations of these substances in body fluid remain unchanged, within limits,

despite changes in the external environment. Efficient regulation of these parameters is

necessary to support life, on all levels. Therefore, the ability to evaluate homeostatic

controls would facilitate the decision process leading to the neonate's discharge from the

unit. The events surrounding the control of homeostasis will be discussed further in

Chapter 2.

We hypothesize that an index of HRV will provide sufficient insight into the

maturity level of the neonate's ANS. In turn, this will provide us with information

regarding the ability of the neonate to control events of homeostasis. In a study by









Chatow et al, increasing conceptual age, CA, or gestational age, GA, is correlated with

neonatal autonomic control. Chatow used a power ratio, discussed in Chapter 2, which

was assumed to reflect the sympathovagal balance of the ANS [3]. In a study of the fetal

baboon, Stark et al. suggest that "Sequential decreases in fetal heart rate, increases in RR

interval (RRi) variability and increases in changes in RRi and ARRi with age imply an

overall maturation in autonomic cardio-regulatory control processes. Increases with

gestation of high frequency components of variability are compatible with enhanced

parasympathetic modulation of the fetal heart rate" [4]. This thesis will attempt to

expand upon these studies to examine developmental changes, cardiac orienting

responses, and to introduce a maturation index based on time and frequency-domain

analysis of neonatal HRV. Current methods of frequency domain analysis tend to center

around one method of analysis. Therefore, new methods of analysis will be constructed

as a result of combining the most functional signal processing tools described in recent

literature. This thesis will explore and compare a variety of spectral analysis methods and

focus on those methods deemed most clinically significant. In addition, we will parallel

techniques used to examine the frequency spectrum of HRV with time-domain

techniques to provide a complete clinical evaluation of each patient and grouping.














CHAPTER 2
AUTONOMIC NERVOUS SYSTEM AND HEART RATE VARIABILITY

Autonomic Nervous System

The ANS is part of the human central nervous system (CNS). The CNS can

essentially be divided into three closely interacting sections: the motor nervous system,

sensory nervous system, and autonomic nervous system. The ANS originates from the

brainstem at the base of the spinal cord. It includes the nerves, which innervate the

smooth muscles of the heart, internal organs, and glands. The ANS is subdivided into the

sympathetic (SNS) and the parasympathetic (PNS) systems, which originate from

different areas in the brainstem and spinal cord.


Sensory
Nervous Syste


Figure 2-1. The divisions of the central nervous system.









Both systems work simultaneously in an antagonistic nature, with an increased

activity in one, inevitably causing a decrease in the other. However, neither substructure

of the ANS is ever completely disabled by the response of the other. While both the

branches of the ANS serve to regulate the internal organs by stimulating or inhibiting

their activity, these two substructures respond to different biochemical transmitter

substances [5]. The actions of the PNS are started by the release of acetylcholine.

Parasympathetic nerve fibers slow the heart rate, narrow the bronchioles, and cause the

release of insulin and digestive fluids. On the other hand, the actions of the SNS are

started by the release of norepinepherine. SNS responses include relaxing the tubes of

the bronchioles, speeding up the heart, and slowing the digestive tract muscles, while

increasing the release of glucose from the liver [5]. Given the broad spectrum of

physiological events controlled by the ANS, it is important to establish criteria to

evaluate its tone. Past research has indicated, through the measurement of heart rate

variability, an inference can be made regarding the functionality of the ANS [3,6,7].

HRV is defined as beat-to-beat fluctuations in heart rate attributed to electrical

impulses generated by the sino-atrial (SA) node. The SA node is commonly referred to

as the heart's internal pacemaker, and the vagus nerve, or the Xth cranial nerve, inhibits

this pacemaker. The PNS influences the SA and AV nodes via the vagus nerve. This

influence is produced by release of the neurotransmitter acetylcholine at the vagus nerve

endings, which result in the slowing of activity at the SA node and a slowing of the

cardiac impulse passing into the ventricles. The SNS has the opposite effect through the

release of norepinepherine at the sympathetic nerve endings, thereby increasing the heart

rate [5].









HRV reflects the response of the ANS to exogenous or endogenous perturbations.

Numerous factors contribute to HRV including blood pressure, temperature, respiration,

state of oxygenation, ventilation, biochemical influences of the acid-base balance, and

psychological parameters [3]. Continuous changes in sympathetic and parasympathetic

neural impulses induce changes in heart rate and cause oscillations around the mean,

defined as HRV. By studying heart rate variability, an opportunity exists to study cardiac

dynamic behavior and functionality of the ANS. Since the Autonomic Nervous System's

control of heart rate is dynamic and nonlinear, it is necessary to explore viable

parameters, which can be used to access the causes of these beat-to-beat fluctuations.

This thesis will address time-domain, linear frequency-domain, and nonlinear frequency-

domain characteristics of the HRV. Thus the framework for the study has been set, with

an attempt to answer each of the following questions:

1. What clinically significant features can be extracted from HRV?
2. What are the developmental changes in HRV?
3. Can we establish criteria to assess vagal tone?
4. Can the premature infant learn or demonstrate a cardiac orienting response, and
does HRV change as learning emerges?

Developmental Changes

Immediately following birth, the heart rate is extremely sensitive to physiological

state of the infant. The developmental changes in HRV are hypothesized to indicate a

change in the parasympathetic tone. While early fetal life indicates a strong dependence

on the SNS, these two systems trend toward an equal balance in the term infant [3]. This

thesis assumes a parallel between the fetus and the premature infant. The developmental

changes in heart rate variability can be measured through time and frequency domain

statistics. Developmental trends of HRV will be measured across subject groups and









gestational age. Time and frequency domain indexes will be used to access maturation of

the neonate's ANS, or ability to regulate homeostasis.

Vagal Tone and Homeostasis

Cardiac vagal tone is a construct that describes the functional relationship between

the brainstem and the heart. Cardiac responses to brief low intensity stimuli should be

deceleratory, indicating receptiveness to incoming stimuli. Cardiac responses to

sustained stimuli should be acceleratory, indicating attempts to shut it out [8]. Vagal tone

reflects individual differences in degrees of influence of the vagus on the heart rate. High

vagal tone is associated with the ability to attend selectively to stimuli and maintain

attention during a task, the ability to maintain homeostatic integrity, and infers healthier

ANS functioning over the long-term.

Homeostasis is primarily maintained by the PNS and provides a physiological

framework for the development of complex behaviors. Therefore infants with a high

parasympathetic tone are more efficient regulators of homeostasis than those dominated

by sympathetic responses [6]. The PNS is primarily responsible with coordinating the

life-support process, which intrinsically provides maximum growth and restoration of the

body's tissues and organs [9]. More effective regulation of homeostatis may be an

important factor underlying the relationship between high resting vagal tone and

favorable neurobehavioral development [6]. Moreover, vagal stimulation controls life-

support processes including sucking, swallowing, thermoregulation, and vocalization [9].

The search for an accurate measurement of parasympathetic and sympathetic tone has

only begun to express its value. This thesis is anticipated to contribute to the

comprehension of these physiological parameters.









Cardiac Orienting Response

The existence of a cardiac orienting response (COR) suggests that HRV includes a

neurongenic response. The deceleratory parasympathetic component is provided by the

bi-directional vagus nerve, which includes both efferent and afferent fibers. Efferent

fibers originating in the brain stem terminate on the SA node. Vagal afferent fibers

originate throughout the heart and project to the nucleus tractus solitarius [10]. These

fibers serve as a continuous response and feedback mechanism between the heart and the

brain, facilitating regulation of cardiac function. Porges et al. suggest that the time

course of the response, the effects of neural blockades, and studies with clinical

populations support this theory. As noted by Porges the heart rate deceleration associated

with the cardiac orienting response is rapid and usually returns rapidly to baseline. Also,

the latency characteristics of the cardiac orienting response are similar to the opto-vagal,

vaso-vagal, baroreceptor-vagal and chemoreceptor-vagal reflexes [8]. Since short latency

heart rate reactivity is mediated by the vagus nerve, the magnitude of the cardiac

orienting response may be an index of vagal regulation.

Respiratory Sinus Arrhythmia

Repiratory sinus arrhythmia (RSA) has long been used to estimate vagal tone

because of its association with PNS. RSA results from increases in vagal efference

during exhalation, which decelerates heart rate, and decreases in vagal efference during

inhalation, which accelerate heart rate. Recently, the accuracy of the correlation between

RSA and vagal tone has been questioned because it only partially accounts for the beat-

to-beat fluctuations of the heart [10]. This thesis will address RSA as a possible criterion

to extract vagal tone, but due to the nature of the study will not pharmacologically prove






9


the efficiency of this technique. However, since RSA accounts for a significant amount

of variability in heart rate, it is impossible to completely ignore this aspect.















CHAPTER 3
PRELIMINARY RESULTS AND STUDY PROCEDURES

Preliminary Findings

In a preliminary study by the principal investigator [7], cardiac activities during

two consecutive quiet periods of 16 healthy fetuses were measured. The exploratory

phase used a repeated measures design to follow the normative changes in fetal HRV

between 28-34 weeks. Changes in HRV were explained by normalized power shifts in the

Fast Fourier Transform (FFT) frequency spectrum. A multi-group pretest-posttest

experimental design determined whether 28-34 week-old fetuses would become familiar

with a rhyme, recited by the maternal parent. Developmental changes in the HRV

frequency spectrum occurred in a nonlinear fashion and individual differences were

apparent. The total power diminished significantly with increasing gestational age.

Furthermore, the study provided conclusive evidence that the 28-34 week old fetus was

capable of responding to a nursery rhyme their mother recited aloud. The response to the

nursery rhyme was described by a shift in the power of HRV frequency spectrum to the

high frequency bands, which is synonymous with an increase in parasympathetic activity.

This conclusion parallels all hypothesized learning capabilities of the premature infant.

The testing sessions for the preliminary study utilized a recording of a rhyme

spoken by an unfamiliar female, which insured the fetuses' reactions would be elicited by

familiarity with the acoustic properties of the rhyme. Once learning was detected,

ongoing interactions between learning, movement and HRV were detected. A









nonparametric signed test was used to compare differences across time. The preliminary

study revealed a significant downward trend in fetal movement and an upward trend in a

high (.17-.40 Hz) frequency once a cardiac orienting response or learning was identified.

Since the subject size was small, there exists a need for a more expanded and complete

study before any conclusions can be drawn. However, the findings of the fetal study

provide the framework for this expanded research, conducted in collaboration with the PI

of the preliminary study.

Proposed Study

A convenience-based sample of 28 low-risk premature newborns admitted to the

NICU at Shand's Teaching Hospital are recruited, in blocks of four, and randomly

assigned to one of two groups. The target age for admission into the study is 27-28

weeks post-conceptional age. Group 1 will begin exposure to the nursery rhyme at 28

weeks and those in Group 2 will begin exposure at 32 weeks. Staggering the groups shall

assist in determining whether an orienting response was a function of learning or simply a

reflexive response to presentation of the rhyme. All subjects will be excluded from

analysis if they are subject to one or more of the following:

* An abnormal head ultrasound
* A sensorineural hearing loss
* A confirmed prenatally transmitted virus/bacterial infection
* Cardiac abnormalities

Table 3-1. Organization of Recitation and Test Sessions

Group 28 weeks 29 weeks 30 weeks 31 weeks 32 weeks 33 weeks 34 weeks
1 Recit Recit Recit Recit Recit Recit Reie
____Test Test Test Test Test Test Test
2 Test Test Test Test Recit Recite Recit-
Only Onlnl y nly Only Test -Test Test









Recitation Sessions

The mother of the infant will complete a recording of a nursery rhyme, identical to

that presented in the preliminary study. The nursery rhyme, taking approximately 15

seconds to recite, will be played 3 consecutive times to provide the infant with 45s of

auditory stimulation. The maternal recitation sessions will occur once in the morning and

again in the evening, each day; Group 1 will begin at 28 weeks and Group 2 at 32 weeks

of age. The purpose of the recitation session is to introduce the neonate to the rhyme

using a familiar voice.

The nursery rhyme will be presented once the premature infant is judged to be in a

state of active sleep and at least fifteen minutes following a feeding [11]. Between 28-34

weeks post-conceptional age active sleep accounts for approximately 70% of the entire

sleep-wake cycle [12]. Learning in the normal newborn, can occur during periods of

active sleep [13,14] and have occurred in the fetus during quiet periods of their quiet-

activity cycle [7]. Criteria established by Thoman and Whitney [15] and Holditch-Davis

[12] will be used to detect active sleep. The subject is in an active sleep state when:

* Eyes closed
* Respiratory rate under 25 breaths per minute
* Limited movement of extremities

Subjects will be monitored at the same time each week because of potential

circadian influences on heart rate patterns and movement [16]. A hospital log of

recitation times and the time since the neonate's last meal will be kept because variations

related to unexpected changes in hospital routines are expected. These recitation sessions

will be completed entirely separate from the testing session described below.









Test Sessions

Each group will undergo a 480 second ECG recording during each test session.

Once more, the criteria established by Thoman and Whitney [15] and Holditch-Davis [12]

will be used to detect active sleep before the test session begins. The ECG recording will

be divided into two stages. The first stage of this study incorporates an exploratory, two-

group pretest-posttest subject set which will be used to methodologically replicate the

findings from the preliminary study in the premature newborn during a similar

developmental time period (28-34 weeks post-conceptional age). The exploratory design

will follow specific changes in HRV, measured by time series analysis, power spectral

analyses, and time dependant spectral analysis over 28-34 weeks post-conceptional age.

A 240-s window separated from the delivery of the auditory stimulus and commencement

of the second phase will be used to represent developmental changes in HRV. It is

hypothesized that developmental changes in HRV will occur in a decreasing, nonlinear

fashion over 28 to 34 weeks gestational age. Two interactions (gestational age vs. gender

and gestational age vs. group assignment) will be included in the mixed general linear

model. The mixed model allows for the individual random effects when estimating fixed

population effects [17].

During a second period of confirmed active sleep, a two-group pretest-posttest

experimental design will be used to document when a cardiac orienting response to a

presented nursery rhyme first occurs, and how movement and HRV change once it is

detected. Again, a nursery rhyme identical to that presented in the preliminary study will

be used. An unfamiliar, English-speaking female performs the recitation of the nursery

rhyme used in these test sessions. To avoid additional exposure to this unfamiliar voice,

this person will not be in attendance for presentations of the rhyme to the infant. The









nursery rhyme takes approximately 15 seconds to recite, and will be repeated 3

consecutive times to provide the 45s window of auditory stimulation. The testing

sessions using the unfamiliar recording of the nursery rhyme will occur once a week for

both groups from 28-34 weeks post-conceptional age.








Variable 240 s 45 s 45 s 45 s


Washout Developmental Prestimulus Stimulus Poststimulus
Analysis




Figure 3-1. Time-line for the windowing of test data.

This second stage will be separated into three 45s windows, in which statistical

analysis of heart rate variability will be individually analyzed. The first window will be

used to detect the baseline state of the infant before the presence of any auditory

stimulation. The second window, immediately following the first, chronologically

corresponds to the recitation of the nursery rhyme. This window will be used to observe

time and frequency domain changes in HRV, which may reflect the hypothesized COR.

A third consecutive window will be used to examine how the neonate responds after the

auditory stimulation has been removed from the environment. The third window will be

used to verify that an observed COR response is not merely a result of coincidental

changes in HRV. In addition, a count of movements will be recorded onto the hospital

log during each test session. Our hypotheses pertaining to a COR suggest, (a) over 28-34






15


weeks of gestational age, it will take 2-5 weeks of recitation to detect and maintain a

cardiac orienting response in each subject and (b) once a cardiac orienting response has

been observed, subsequent testing session will demonstrate a trend of movement counts

and specific HRV characteristics.














CHAPTER 4
TESTING

The rhymes will be played over a 2.5-inch Radio Shack speaker positioned 20 cm

from the infant's ear. Both versions of the rhyme (maternal and unfamiliar) will be

delivered at 55 dB (+ 5 dB) as measured by a Bruel & Kjaer sound level meter (2239)

using A-weighted sound pressure levels [18]. This decibel level is just below the normal

levels of conversational speech (58-60 dB) and is significantly lower than the background

sound levels reported in the literature (70-80 dB + 20 dB) for NICUs [19,20]. The

decision to present both versions of the rhyme (maternal and unfamiliar) at 55 dB (+ 5

dB) was based on preliminary findings in the premature fetus where a sound level of 75

dB was used [7]. This decibel level corrects for a 20-decibel attenuation due to passage

of sound through the uterine wall [21], bringing the sound level down to approximately

55 dB or to just below the level of normal conversational speech (58-60 dB). The study

will record decibel levels at the infant's ear prior to presentation of the maternal and

unfamiliar rhyme. In the preliminary study, this sound level (75 dB) did not elicit a

cardiac deceleration until a history of recitation of the nursery rhyme had occurred [7].

An RS232 interface from the Agilent Neonatal CMS 2001 neonatal monitor to an

IBM compatible computer (Dell Inspiron 8100 Laptop) was constructed to obtain the

ECG signal and respiratory rates. The data acquisition program is designed to run in

DOS and may be activated from the Windows Desktop of Dell Inspiron 8100 Laptop by

double clicking on the "Neonatal Data Acquisition" icon. The C program, originally

written by Hewlett Packard (HP), was modified from an earlier version to in order to









record the ECG and respiratory rate. The HP program was initially written to

demonstrate that information from their monitors could be transferred and displayed on

peripheral devices through RS-232 interfacing. The sample program was only available

in a DOS format. Although this format is not optimally user friendly, an extensive user

interface was not necessary to transfer the required signals into readable text files.

Furthermore, HP uses a complex communication protocol, which makes it desirable to

avoid writing an original program. The original program was capable of display the ECG

signal on a peripheral monitor, which implies that the ECG information was already

being transferred. The program modification involved the following:

* Requesting the respiratory rate signal and displaying the information on the screen.
* Creating a user interface to request the patient number.
* Transferring the requested information to a readable text file.


Before data acquisition can begin, the Agilent monitor must be set to a baud rate of

19200 and the transition must be set to a high to low trigger. If the baud rate or transition

is not correctly set, the program will time-out after eight seconds. When the data

acquisition program is started the user enters the patient number and is asked to confirm

this number before the recording process begins. The program passes a text string, which

signifies the rhyme state and allows the user to mark the data when the rhyme is

presented and upon completion. The three measurements: ECG, respiratory rate, and the

rhyme state are recorded in an ASCII text file labeled with the month, date, hour (in 24h

format), and minute (MMDDHHmm.txt) the program is initialized. ASCII was the

chosen format for data storage to permit the user to utilize other analysis programs (i.e.

Microsoft EXCEL). The parameters are stored in the format directly readable as a

MATLAB matrix and can later be separated into three separate vectors. Preceding the











matrix is an ASCII timestamp, which is recorded once the program has finished the


initialization process (DDHHmmss), and appended by another timestamp when the

program is terminated. These time stamps allow the user to verify the length of the

recording.

The program is designed to work with an ECG/Resp module, which will transfer


both the ECG and the respiratory rate signal to the Dell Inspiron 8100 Laptop. However,

in the absence of the ECG/Resp module it is possible to use only the singular ECG


module. This will not record the respiratory signal, but will still record the ECG signal.

The first column recorded in the ASCII text file represents the unsigned 12-bit digital


values of the ECG signal. The second column contains the information regarding the

rhyme state. The final column contains the respiratory rate. The resulting format of the


text file is illustrated in the figure below.




06101635
1985,0,51,
1996,0, 51,
1996,0, 51,
1997,0,51,
1997,0, 51,


2027, 1,58,
2028, 1 58,
2028, 1 58,
2028, 1, 58,
2049,1, 58,


2079,0,61,
2078,0, 61,
2089,0, 61,
2078,0, 61,
2078,0, 61,
06102452



Figure 4-1. Example of text file generated from Data Acquisition Program









All columns are recorded at a sampling rate of 500Hz, which was preserved from the HP

demo program as the rate information was sent to the screen. This sampling rate, much

higher than the 200Hz often described in literature, was preserved to prevent jagged QRS

complexes. The column corresponding to the state of the rhyme is initialized to '0', until

the 'r' key is pressed and the value changes to '1'. The value remains at '1' until the 'd'

key is pressed, which signifies the end of the presentation and returns the value to '0'.

When all desired data has been recorded, the user presses the 'q' key twice to exit the

program and return to the Windows desktop.


Figure 4-2. The setup for the test sessions.









Figure 4-2 illustrates the equipment setup during the test sessions. A male 25-pin

RS232 port in the back of the Agilent monitor is interfaced with the female 9-pin RS232

port Dell Inspiron 8100 Laptop. The RS232 cable provides means of serial

communication between the Agilent monitor and the Dell Laptop. The Agilent Neonatal

CMS 2001 monitor continually monitors the infant's vital signs, with a standard 3-lead

ECG setup. The speakers, connected to the CD CopyWriter Live, are placed 20cm from

the infant's ear. This setup is replicated for every test session, regardless of the subject's

group number.














CHAPTER 5
QRS DETECTION

Once all data from the testing session is acquired the focus shifts to the data

analysis. The first goal is to perform time-series operations, yielding discernable

characteristics of the HRV. The HRV is basically controlled by the effect of the vagal

and sympathetic influences of the sino-atrial node. Therefore, it seems the best

representation of the HRV variability would be expressed as a function of the time

intervals between consecutive P waves. However, due to the inconsistent positioning and

suboptimal attachment of the infants ECG leads, P waves characteristically have a low

signal-to-noise (SNR) ratio [22]. Therefore, in the absence of cardiac arrithymias,

addressed in the exclusion criteria, it seems more reliable to replace the detection of PP

intervals (PPi), with RRi. While this may include an extra variable in the analysis of the

HRV, it seems trivial compared to the accuracy of the detection of PPi. It is our

assumption that variations in the PR intervals will be equal or less than error of detection

of the PPi. When considering an error-percentage of P-wave detection greater than the

standard deviation of the PR interval and greater than the percent error in R-wave

detection, the replacement of PPi with RRi becomes trivial. The occurrence of the

following cardiac arrhythmias causes this assumption to fail [23]:

* ATRIAL FLUTTER the atrial flutter waves, known as F waves, are larger than
normal P waves and they have a saw-toothed waveform. Not every atrial flutter
wave results in a QRS complex because the AV node acts as a filter. Some flutter
waves reach the AV node when it is refractory and thus are not propagated to the
ventricles. The ventricular rate is usually regular but slower than the atrial rate.









* ATRIAL FIBRILLATION Atrial fibrillation occurs when the atria depolarize
repeatedly and in an irregular uncontrolled manner. As a result, there is no
concerted contraction of the atria. No P-waves are observed in the EKG due to the
chaotic atrial depolarization. The chaotic atrial depolarization waves penetrate the
AV node in an irregular manner, resulting in irregular ventricular contractions. The
QRS complexes have normal shape, due to normal ventricular conduction.
However, the RR intervals vary from beat to beat.

* VENTRICULAR TACHYCARDIA Ventricular tachycardia occurs when
electrical impulses originating either from the ventricles cause rapid ventricular
depolarization. Since the impulse originates from the ventricles, the QRS
complexes are wide and bizarre. Ventricular impulses can be sometimes conducted
backwards to the atria, in which case, P-waves may be inverted. Otherwise, regular
normal P waves may be present but not associated with QRS complexes. The RR
intervals are usually regular.

* THIRD DEGREE AV BLOCK Atrial rate is usually normal and the ventricular
rate is usually low, with the atrial rate always faster than the ventricular rate. P
waves are normal with constant P-P intervals, but not associated with the QRS
complexes. QRS may be normal or widened depending on where the escape
pacemaker is located in the conduction system. The result is unrelated atrial and
ventricular activities due to the complete blocking of the atrial impulses to the
ventricles.

Subjects whose electrocardiogram display any of the arrhythmias listed above are

removed from analysis on the basis of violating the exclusion criteria. In conclusion, all

ECG time-domain analysis of the HRV will consider RRi rather than the ideal analysis of

the PPi for the remainder of the thesis.

The first step in the analysis of the HRV is to extract the RRi from the recorded

ECG waveforms. The ASCII text formatted matrix addressed in the data acquisition

chapter is read into MATLAB. The result is an N by 3 matrix, with N being the number

of samples. This matrix is then separated into 3 vectors of length N. Vector one contains

the ECG information; vector two contains the respiratory rate; and vector three contains

information regarded the state of the rhyme. Only vector one expresses useful

information for analysis of RRi. Figure 5-1 describes the transform of the raw ECG


signal to a vector containing RRi.






























Figure 5-1. QRS detection algorithm

Noise/Baseline Wander

The difficulty in accurate R-wave or QRS complex detection is directly related to

the common noise characteristics existing in a raw ECG signal. Noise sources include,

but are not limited to: electrode motion, muscle movement, power-line interference,

baseline drift, and T-wave interference. Using a method similar to that described by

Jiapu Pan [24] we will attempt to improve the SNR of the raw ECG signal and detect the

QRS complexes. This process is split into two linear steps and one nonlinear

transformation. The first linear process requires filtering the signal with a bandpass

filter, which is split into a low-pass and high-pass filter. A low-pass filter with cutoff

frequency of approximately 11Hz is first applied the raw ECG signal to remove the low-

frequency baseline drift. In addition, the filtered ECG signal is sent through a high-pass

filter to remove muscle movement, powerline interference, and T-wave interference. The

cutoff frequency of the high-pass filter is approximately 5Hz. This yields a filtered ECG









signal (fECG) containing frequencies of 5-1 1Hz, which has been described to preserve

maximum information regarding the QRS complexes [24].

Low-pass filter:

Hz) Z-6 2
H ( 1)2

High-pass filter:

H(z) (- 132z-16+ -32)
1( z1)

Filtering with the Pan-Tompkins method described above, yielded only minor

improvements to the raw ECG signal. Problems with the above technique arise as a

result of the electrode placement differing from subject to subject. In some case, the

Pam-Tompkins algorithm removed substantial frequency components of the QRS

complex, thereby deteriorating the quality of the signal. Therefore, we need to replace

the low and high pass filters described above with a notch and bandpass filter. The

nature of the ECG signal calls for the use of an IIR filter. We will explore the

Butterworth and Chebyshev filters. Both are high order filter designs, which can be

realized by using simple first order stages cascaded together to achieve the desired order,

passband response, and cut-off frequency. The Butterworth filter yields no passband

ripple and a roll-off of -20dB per pole. However, the flatness of the passband response

comes at the expense of roll-off. The Chebyshev filter displays a much steeper roll-off,

but is characterized by a significant passband ripple. The Chebyshev filter thereby

optimizes roll-off at the expense of passband ripple. Due to the nature of the ECG signal

the roll-off frequency of the passband is not as important as the gain in the passband.

Therefore, we will use Butterworth filters for all filtering of the raw ECG signal.











Since 60Hz powerline interference is common among electrical signals, we need to

construct a notch filter to remove this noise. The notch filter used is a 5th order


Butterworth filter with center frequency at 60Hz. The magnitude and phase of this filter


is shown in Figure 5-2.


Notch Filter
50






.100

-150 -
0 10 20 30 40 50 60 70 80 90 100
Frequency (Hz)


.200 -
S -400


-ao
-1000 ---_
-1200 I i i
0 10 20 30 40 50 60 70 80 90 100
Frequency (Hz)



Figure 5-2. Magnitude and phase response of the 60Hz notch filter used to filter the raw
ECG signal.

The phase of this filter is essentially linear over the region of attenuation. After the raw


ECG signal is passed through the notch filter, it is passed through an 8th order bandpass


filter serving to remove baseline wander and high frequency noise; thereby improving the


SNR of the raw ECG signal. Both the notch and bandpass filters are implemented using


Butterworth filters. As illustrated in the Figure 5-3 below, the phase over the passband


region is approximately linear. Keeping the phase of the filter linear over the passband


prevents phase distortion in the resulting filtered ECG waveform.












Bandpass Filter

i--------- ---I----------- ------------- *
20




-40
II I I I I i
0 10 20 30 40 50 60 70 80 90 100
Frequency (Hz)

1500

S1000

500 .



-500
.500 --- I -- I -- I----------------1 --- I -
0 10 20 30 40 50 60 70 80 90 100
Frequency (Hz)



Figure 5-3. Magnitude and phase response of the bandpass filter used to filter the raw
ECG signal.

During the recording phase, we obtained ECG recordings that included significant


periods where the signal underwent high and low voltage clipping. The clipping resulted


in epochs of flat line, thereby making it impossible to detect subsequent QRS complexes.


For HRV analysis, a continuous ECG signal is needed. The undesired clipping resulted


in extensive signal losses that required considerable epochs of QRS interpolation. It was


later discovered that this problem could be overcome by setting the ECG module to filter


the signal before it was relayed to the laptop for text conversion. In addition, incorrect


electrode placement proved the monitors default selection of Lead II to be an errant


assumption. Lead selection was then delegated to the user before the acquisition program


was initialized. After realizing the need to select the "filter" option and correct lead, we


were able to attain signals that needed little or no secondary filtering and was not


necessary to bandpass filter the recordings. Powerline interference during the passage











from the monitor to the laptop was not completely alleviated and the 60Hz notch filter,


described above, was still applied. When the detection phase demonstrated a need for


secondary filtering, the bandpass Butterworth Filter, described above, was applied to the


signal to conclude the filtering stage of the detection algorithm.


Slope Detection

In parallel with the Pan-Tompkins algorithm, after filtering, the fECG signal is sent


through a differentiator yielding QRS complex slope information; this method results in a


differentiated version (dfECG) of the fECG. In accordance with Pan [24], we will use a


5-point transfer function to approximate the derivative.


H(z) =(1/8T)(-z2 -2z1 +2' +z2)


where T is the sampling period.


This transfer function is roughly linear in the region between dc and 100Hz. Therefore it


provides a quick approximation of the ideal derivative over the region of interest.


Differentiator
20 ---------F-- --I--------
10

0 /
-10


-30

0 10 20 30 40 50 60 70 80 90 100
Frequency (Hz)

-50

-100





-2500
0 10 20 30 40 50 60 70 60 90 100
Frequency (Hz)



Figure 5-4. Magnitude and phase response of the differentiator function









Following this step, the dfECG signal is passed through a nonlinear point-by-point

squaring operator yielding the sdfECG signal:

y(nT)= [x(nT)]2

This operation makes all data points positive and supplies a nonlinear amplification of the

signal at the higher frequencies. Finally, the sdfECG is smoothed using a ten point

averaging function, which produces a signal containing the power information of the

ECG waveform's slope, pECG. This reduces the presence of local minima in the sdfECG

signal. From this point, the method used for detection of QRS complexes deviates from

that used by Pan [24]. Since it is not necessary to perform QRS detection in an online

fashion, we have the luxury of determining thresholds using the entire information

contained in pECG signal. The threshold limit of the QRS slope detection is set in the

following fashion:

S= [max(dECG 2)+mean(dCG 2 )]/ 2

The selection of the current threshold (y) allows for detection of poorly defined QRS

complexes, while distinguishing the R-waves from occurrences of P and T-waves. The

threshold is then adjusted throughout the detection process. The following equations

illustrates how the threshold is updated:

77 = [max(dfCG (i: i + win))+mean(dfECG (i : i +win))]/2

where threshold parameter (rl) is updated in accordance with the current sample (i) over a

user-defined window (win)

S6=-y


;i y(g*3)









where 6 represents is the threshold difference and includes a memory constant (E) to

reduce high frequency noise interference during the update

Once this threshold has been set, it is now necessary to find maximum values

during each interval that the signal spends above threshold. The detection algorithm

searches the pECG until it reaches a y-value above the threshold. An array is

constructed, storing the waveform from this point until the y-value falls below the

threshold. The x-value corresponding to the maximum y-value inside this array is stored

as a possible QRS complex.

QRS = max(x[n,n + 1,n + 2, .,m]

where y[n] and y[m] are above the threshold, but y[n-1] and y[m+l] fall below the

threshold.


Figure 5-5. Visualization of QRS peak detection









These maximums, or peaks, define the occurrences of possible QRS complexes. If a

QRS complex is detected outside the expected range given below,

mean(RRi) < QRS(k) < (mean(RRi) A)
A

5. RR() = 0.375s (160bpm)
6. 1 7. QRS(k) is the detected QRS complex


the algorithm reverts back to the previous QRS complex, adjusts the threshold as follows,


/ threshold

Threshold is a user definable parameter (default = 1.5)

and searches for a new QRS complex. This searchback technique plays a crucial role in

preventing high frequency noise contamination, P-waves, and T-waves from being

marked as a QRS complex.

If a complex is not found within the expected range after 5 searchbacks, a

complex is interpolated as an average of the previous two intervals.

QRS(k) = [QRS(k 1) + QRS(k 2)/2

Subsequent detections are stored in a vector of increasing length until the algorithm has

operated on the entire pECG signal. In addition, the interpolated detections are stored in

a separate vector for further analysis. Units of these vectors are stored according to

sample numbers. The vectors are then divided by the sampling rate for unit conversion to

the time-domain. Finally, another vector is constructed, which represents a time-based

representation of the RRi in the following fashion:

RRi = {QRS(2)- QRSl), QRS3) QRS2),..., QRS(N) QRS(N 1)










where N = total number of QRS detections

Figure 5-6 illustrates the final results of the QRS detection algorithm where the

QRS complexes are marked in red. This figure also displays a linear interpolation of the

RRi waveform corresponding to the detected beats.




QRS Complexes Marked in Red
3500

3000

S2500-

S2000--

1500

1000 --
100 100.5 101 101.5 102 102.5 103 103.5 104
Time(s)
Beat Intervals


0.365


D 0 -355/
i =\ / /

0.35
0.35 5------------
100 100.5 101 101.5 102 102.5 103 103.5 104
Time (s)


Figure 5-6. Example of detected QRS complexes and beat intervals

The detection is offset from the peaks of the QRS complex due to the intrinsic phase

delays of the filtering process. As a result, the delays are consistent throughout each

signal and fall out after computing the RRi vector, as follows:

RRi= QRS(N)- ] [QRS(N -1)- rj}
= {QRS(N)-- QRS(N-1)+ ri
= {QRS(N) QRS(N 1)}









Outliers

Outliers are ill-defined byproducts of data acquisition systems. There is no

rigorous definition of an "outlier"; it is only generally conceived that selectively

eliminating inconvenient data points from analysis creates a more robust dataset.

Outliers are acknowledged to skew sample means, but tend to have less affect on the

median than the mean. For the purpose of this thesis, it is necessary to form a working

definition of an outlier.

An outlier is a data point that is an "unusual" observation or an extreme value in the
data set.

The lack of a universally accepted definition of an outlier inherently causes a deficiency

of objective mathematical processes for outlier recognition. Adding to the complication

of outlier detection, RRi datasets are time dependant vectors, rather than independent

samples. Therefore, it is essential to add to our working definition of an outlier.

A time dependant outlier is a data point that is an "unusual" observation or an
extreme value in the data set, which is not part of a trend in time.

Again, the definition of a "trend" is subjective and needs to be defined for the purposes of

this paper.

A trend is two or more consecutive data points, which move in the same general
direction within a given statistical range.

These definitions will be used to construct an algorithm to remove time dependant

outliers from the RRi datasets. False detections and missed detections result in possible

outliers in these datasets. The presence of outliers in the RRi dataset adversely effect

both time and frequency domain analysis of HRV. Possible skewed time domain

representations of HRV include the mean and standard deviation of RRi windows. In

addition, these outliers produce "ghost powers" in the high frequency range of the










spectrum. Not only does this increase the power in the high frequency bands, it skews

the power in the low frequency bands. Low frequency bands are of the foremost

importance during analysis of HRV in the frequency domain. This will become more

evident through discussion in Chapter 6.

Now that the importance of removing outliers has been established, we must

construct a mathematically sound algorithm for their removal. Three techniques are

explored as possible strategies for outlier detection: direct RRi analysis, slope analysis,

and differential analysis. For each of the aforementioned techniques, we will use standard

deviation as the exclusion criterion. Since literature is not conclusive for the distribution

of RRi datasets, the exclusion criterion is based on normal distribution of the data. The

figures below illustrate the RRi distribution and the normal distribution. Visually, the

normal distribution seems an accurate estimation of the RRi dataset. We will discuss the

validity of this assumption further in Chapter 6.


Distibuton of RRi
160 I



120

100

80

60






0.28 0.3 032 0.34 0,36 0.38


Figure 5-7. Histogram of RRi distribution from one subject.


0.4 0.42




























Figure 5-8. Standard deviations of the Normal Distribution

In a normal distribution, one standard deviation away from the mean in either

direction on the horizontal axis accounts for around 68 percent of the values in the

dataset. Two standard deviations away from the mean includes for roughly 95 percent of

the data. Finally, three standard deviations account for about 99 percent of the data [25].

Therefore, we will assume data outside of three standard deviations from the mean to be

statistically irrelevant and tag them as outliers unless they are part of a trend.

Direct RRi Analysis

The simplest method of removing outliers would be to examine the RRi dataset

directly and remove data that fall three standard deviations above or below the mean.

The following equation is used for outlier detection:

outliers = RRi(n)> mean(RRi)+ 3 std(RRij


for n = 1: length(RRi)











This method inadequately detects outliers and is not optimal to detect a trend, which

moves beyond the bounds of the acceptable range. In addition, the direct RRi analysis

method does not account for consecutive false or missed detection that exhibit a bouncing

effect around the actual QRS complex.






0.48

0.46

0.44

0.42

0.4

t 0.38

S0.36-

0.34

0.32

0.3-

0.28





Figure 5-9. Illustrates the bouncing effect and upper and lower bound symmetry
problems with using the direct RRi analysis method for tagging outliers.

The figure above demonstrates problems, which arise from direct RRi analysis

method for tagging outliers. The green lines represent the bounds for acceptable

datapoints. The bouncing effect is clearly present in this example. According to the

exclusion criteria, one of the datapoints would be tagged as an outlier, while the other

point would pass through the algorithm without being tagged. This figure also unveils

weaknesses with the computation of the mean. When the heart rate deviates to far from









the mean (marked by the black line), situations exist where the upper and lower bound

are not symmetrical about the data over a given window. Given all the abovementioned

weaknesses, it is necessary to explore a more robust technique to mark the presence of

an outlier.

Slope Analysis

In order to correct for the bouncing effect described above, slope information could

be used to isolate these occurrences and tag outliers. First, compute the derivative of the

RRi dataset, as follows:

aRRi = R t)

Then, tagging outliers using the same detection method as in the direct analysis:

outliers = 8RRi(n) > mean(8RRi)+ 3 std(RRi)

for n = 1: length(ORRi)

This algorithm adequately accounts for the removal of the bouncing effect, but we will

demonstrate that it is impossible to exclude trends from being tagged. Figure 5-10 is an

example of trend which seems to begin abruptly, but is sustained over subsequent

intervals. The slope algorithm solves the problems in accordance with the mean and

bounce effect, but introduces troubles in excluding trends from being tagged as outliers.

Over the course of this trend, only those datapoints whose slope exceeded the bound of

acceptance (marked by the green line) would be tagged. No mathematically simple

method exists to overcome the flaws in this algorithm. Merely modifying this algorithm

to remove the outliers and loop through the process again would serve to completely

eliminate the trend. Therefore, for a second time it is necessary to search for a more

robust algorithm to tag outliers.














RRi Plot
0.6


0.5


0.4


0.3


RRi Derivative Plot
0.15

0.1 -

0.05

0

-0.05

-0.1




Figure 5-10. Illustrates an RRi trend that would not be preserved using the slope analysis
method for tagging outliers.

Differential Analysis

To this point, three problems in the outlier detection process have been identified:

mean calculations that deteriorate the symmetry of the acceptance bound, the bouncing

effect, and the exclusion of trends. We introduce the differential analysis method as a

means, which circumvents all three of these problems. First, a difference vector is

constructed as follows:


RRidff(n) x(n) x(n 1)


for n =2 : m, where m = length(RRi)









RRid (1) = mean{RRi) RRi(l)

The end product is a vector, which is the same length as RRi. The resulting mean of the

difference vector is close to zero, which further implies the data is predominately

normally distributed. The difference vector also alleviates the first problem mentioned

above. Then, a vector is created using the same exclusion criterion discussed in the

previous two subsections. However, using this difference vector alone to tag outliers

outside the acceptable standard deviation range is not sufficient. A single outlier in a

dataset would produce two datapoints in this difference vector that would be tagged as

outliers (n and n+1). Being that it is not necessary to perform the outlier removal real-

time, we are able to create an additional vector, as follows:

RRi^ff (n)= x(n) x(n +1

for n = 1: m -1, where m = length(RRi)


RRidff (m) = mean(RRi)- RRi(m)

The same exclusion criterion is used to create a vector containing the outliers. As a

result of this computation, a single outlier in the dataset would produce two datapoints

tagged as outliers. Since taking the difference of RR intervals backwards through time

generates this vector, the two resulting datapoints tagged as outliers would be n and n-1.

Now, we compare the two difference vectors and note only those values for n, which

occur in both. These are the actual RRi datapoints that are separated from the previous

and subsequent datapoint by three standard deviations above or below the differential

mean. Figure 5-11 is a section of RRi data with trends that would normal be removed

using the Direct RRi Analysis method, described above, and a bounce effect.










Differential Analysis Outlier Removal


Figure 5-11. Outliers that have been tagged, using the differential method, are marked by
a green dot. Shows the efficiency of the method to preserve trend and remove
bouncing effects.

In figure above, the points that were tagged as outliers are denoted by green dots. It is

now evident that we have constructed a robust method to find and remove outliers, as

defined above.

Interpolation

The RRi vector described above is sufficient for time domain analysis of the HRV.

The RRi vector can also be converted to an instantaneous representation of the heart rate

(iBPM). However, there is a problem with using this vector to represent the HRV in the

frequency domain. When plotting the RRi on a time scale, the x-axis must include a

cumulative sum of all the preceding RRi plotted against the RRi or iBPM. This results in










an unequally sampled representation of the HRV. In order to apply a meaningful power

spectral representation of the HRV, the data must be evenly sampled at a known

frequency. Many papers using spectral analysis fail to address this topic and there are no

universal methods to solve this problem.

Resolving this issue of uneven sampling is essential for spectral analysis. Without

addressing this topic, may forms of spectral analysis applied to the RRi would produce

meaningless results. The choice of an interpolation method is very important, as it can

significantly compromise the integrity of the data. Therefore, we will explore four

methods of interpolating between successive data points of RRi:

* Linear Interpolation (LI)

* LaGrange Polynomial Interpolation (LPI)

* Cubic Spline (CS)

* Cubic Interpolation (CI)

LaGrange Polynormi Irtapola
0 38
0,375.





//


0.345[ -
J4- r /


TIme(s)


Figure 5-12. LaGrange Polynomial Interpolation









41




Linear Interpolation


--------------L- ------- -- L
1 2 3 4
Times)


5 6 7


Figure 5-13. Linear Interpolation


Cubic Spine
I I~ -


S 1 2 3 4
Times)


5 6 7


Figure 5-14. Cubic Spline Interpolation


*) 37


0365





S0 25S





n?F


n.,O


'A'







42



Cubic Interpolation
0 38

3.375 -

(137

0A





0,35 X
0 35 -4


0 34




0 1 2 3 4 5 6 7
Time(s)



Figure 5-15. Cubic Interpolation

The performance of the interpolation techniques is shown applied to a sample RRi

dataset marked by X's in the figures above. LI is continuous, but the slope changes at

the vertex points. LPI has an adjustable order of polynomial and reduces to LI when the

order equals one. As seen above, the LPI algorithm does not perform well at the

endpoints. CSI produces the smoothest results of all the interpolation methods, but

performs poorly if the data is highly non-uniform. CI results in both the interpolated data

and its derivative being continuous. The CI technique seems the most realistic

interpolation method because it does not force a sinusoidal interpolation. Therefore, we

will use CI to interpolate and resample RRi in preparation for spectral analysis. It is

necessary to resample the RRi at a frequency where the Nyquist rate is above the

frequency range of interest for spectral analysis, discussed further in Chapter 6. Also in






43


Chapter 6, spectral analysis techniques will be introduced that do not require resampling

of the RRi waveform prior examining HRV in the frequency domain.















CHAPTER 6
STATISTICAL ANALYSIS

Time Based

The simplest methods used to evaluate HRV are based in the time domain. Time

domain statistics can be divided into those derived from RRi or iHR and those derived

from differences in RRi, namely ARRi [26]. The following time-domain techniques will

be used:

* Mean RRi (mRRi)
* Median RRi (medRRi)
* Standard deviation of RR interval (SDRR)
* Poincare plots
* Quadrant analysis


First Order Statistics

Analysis of mRRi and medRRi yields some useful information about HRV. We

will use this information to analyze trends of the mRRi and medRRi with gestational age

and across genders. Also, the mRRi for the 45-s prestimulus period will be subtracted

from each heartbeat of the 45-s prestimulus and stimulus period. The result will be a

string of difference scores. The difference scores will be used to confirm that no

significant HRV trending has occurred during the prestimulus period, which may mask or

falsely demonstrate a COR. During the stimulus period, the difference scores will be

used to detect an orienting response. The mRRi is calculated in the following manner:

N
ZRRi(n)
mRRi = "
N









for N= number of RRi

To find the median, we must first sort the RR intervals by ascending values. Then, apply

the following formula to find medRRi:

medRRi = RRisted (N / 2)

for even values ofN


medRRi = RRsorted ([N -1]/ 2) + RRsorted ([N + 1]/2)
2

for odd values ofN

SDRR is highly dependant upon the length of the recording [26]. When using this

value as a comparison across datasets, we must be careful to uniformly define the length

of the recording. Therefore, we have kept all window lengths constant from subject-to-

subject and session-to-session. The variance (SDRR squared) is mathematically equal to

the total power of the spectrum; therefore SDRR will also be used to analyze trends

corresponding to increasing gestational age. SDRR is the primary first order statistical

analysis method we will use to discuss HRV in the time domain. SDRR is calculated in

the following manner:


Z(RRi(n) -mRRi)2
SDRR = -"
(N 1)

where N= number of RRi

Poincare Plots

In addition to numerical analysis in the time domain, the importance of geometric

representations, such as Poincare plots are useful tools for HRV analysis [26,27]. The

nonlinear Poincare plot produces a figure with each RRi plotted against the previous RRi.









According to Brennan et al., the plot provides summary information as well as detailed

beat-to-beat information on the behavior of the heart. The dispersion of points

perpendicular to the line of identity reflects the level of short-term variability [27].

In subjects with normal cardiac function, the geometry of Poincare plots tends to be

elliptical in nature. No literature pertaining to analysis of Poincare plots using an ellipse

were found, but we will explore an ellipse fitting method as a possible representation of

short-term and long-variability. MATLAB code written by Marcos Duarte [28] will be

used to fit an ellipse to the data. Duarte calculates the major and minor axis of the ellipse

using principal component analysis (PCA). The method results in 85.35% of the

datapoints lying inside of the ellipse. The datapoints that fall outside the ellipse can be

described by increasing the order of the ellipsoid, but it is only necessary to include the

first two eigenvalues when analyzing the short and long-term variability of the dataset.

PCA is a multivariate procedure, which rotates the data such that maximum

variabilities are projected onto the axes. Essentially through computing the eigenvalues

of the correlation matrix, a set of correlated variables are transformed into a set of

uncorrelated variables, which are ordered by reducing variability. The uncorrelated

variables are linear combinations of the original variables, and the last of these variables

can be removed with minimum loss of real data. The largest eigenvalue corresponds to

the principal component in the direction of greatest variance; the next largest eigenvalue

corresponds to the principal component in the perpendicular direction of next greatest

variance.

PCA can be viewed as a rotation of the existing axes to new positions in the space

defined by the original variables. In this new rotation, there will be no correlation







47


between the new variables defined by the rotation. The first new variable contains the

maximum amount of variation. The second new variable contains the maximum amount

of variation unexplained by, and orthogonal to, the first. The result is an ellipsoid in

multidimensional space and reduces to an ellipse for the two dimensional Poincare Plot.

When a Gaussian random vector has covariance matrix that is not diagonal, the

axes of the resulting ellipse are perpendicular to each other, but are not parallel to the

coordinate axes. As shown in the Figure 6-1, the principal components of the RR interval

Poincare representation result in an ellipse whose axes are rotated versions of the

coordinate axes. Therefore, as discussed in Chapter 5, the assumption that a Normal

distribution can describe RRi datasets is valid. The minor axis of the ellipse serves as an

interpretation of short-term variability; and the major axis serves as an index on long-

term variability.


Poincare Plot
0.45

0.44

0.43 *
0-2**
0.42


039,




0-38-

0.37

0.36
0,36 0.37 0.38 0.39 04 0.41 0.42 0.43 0,44 0,45
RR(n)(s)

Figure 6-1. Poincare plot representation of RRi. Since data points overlap, colors,
coordinated with the visible spectrum, are used to indicate magnitude at each
point.









Quadrant Analysis

Quadrant analysis is a graphical tool used to represent the number of sustained

increases, sustained decreases, and alternating changes in RRi. Two vectors are created

as follows:

ARRi(n) = RRi(n) RRi(n 1)

ARR(n +1)= RRi(n +1) -RRi(n)

Each difference value was then plotted against the previous difference, resulting in

the ARRi (n) versus ARRi (n+1) plot shown below [4]. The plot displays values

distributed over an area defined by four quadrants. Each quadrant indicates the direction

of two consecutive changes in interval length.








(-,+) (+,+)



A B ARR+1


C D



(-A -) (+, -)



ARR,


Figure 6-2. Defines the quadrants used in Quadrant Analysis.







49


We will use this analysis tool to examine sustained increases or decreases of RRi


and alternating RRi. The number of points in quadrants A and D measure high frequency


variation, or quick changes in RRi. Low frequency heart rate variation is characterized


by many sequential increases and decreases in interval length, which is indicated by


points in quadrants B and C. If changes in RRi occurred randomly and independently, all


quadrants would have an approximately equal number of data points. In addition, if the


heart rate tended to increase quickly and decrease slowly, more points would fall in


quadrants A and B than in C and D and visa versa. Figure 6-3 is an example of a


quadrant analysis plot extract from real RRi data.


Quadrant Analysis
0.3



0.2



0.1



0
-n


-0.1



-0.2



-0.3
-0.3 -0.2 -0.1 0 0.1 0.2 0.3
deltRR(n)



Figure 6-3. Example of Quadrant Analysis applied to an RRi dataset.









Frequency Based

While analytical methods were first restricted to the time domain, they have

quickly shifted to the frequency domain. Spectral analysis is now a common

representation of HRV. Transforming HRV to the frequency domain provides a means of

separating the parasympathetic and sympathetic responses of the ANS. Heart rate

fluctuations in the premature will be examined in the .00033-1Hz, while the frequencies

below this will be discarded to reduce the influence of slow trends or artifacts. Research

has suggested that the power spectral decomposition of the ECG yields frequency bands

essential to analysis of the ANS. These frequency bands are very inconsistent throughout

the literature referenced for this subject. This thesis will follow the frequency divisions

described by Stein et al. [29]:

* The ultra low frequency band (ULFB) <.0033Hz

* The very low frequency band (VLFB) 0.0033 0.04Hz

* The low frequency band (LFB) 0.04 0.15Hz

* The high frequency band (HFB) 0.15 0.4Hz

* The total power (TP) .0033 1Hz

The low frequency band of the power spectrum is composed of two peaks. The

mid-range peak (.09 .15Hz) corresponds to the Mayer waves, which have been related

to the baroreceptor reflex [3]. In addition, there exist a second peak between .02 .09Hz,

which has been linked to thermoregulatory fluctuations and variations in vasomotor tone.

The high-frequency band reflects the parasympathetic response, while the low frequency

band is influenced by the PNS and SNS [3]. For each frequency band listed above, we

will extract the normalized power over the described range. Also, the power ratio of the









high and low frequency bands has been assumed to reflect the sympathovagal balance

[3]. This power ratio is observed as follows:

LFP
SV
rto (LFP + HFP)


The SV ratio should provide us with insight into the ANS balance and should trend

downward in accordance with increasing gestational age. Lastly, a range just above the

high-frequency band contains information regarding the RSA, which is not represented

by a well-defined peak, but is usually centered about the .5-.6Hz range. We will explore

three methods of power spectral analysis in the following subsections:

* Fast Fourier Spectral Analysis
* Lomb-Scargle Spectral Analysis
* Empirical Mode Decomposition


Fast Fourier Spectral Analysis

The simplest and most commonly exploited technique for spectral analysis is the

Fast Fourier Transform (FFT). The FFT is commonly used because it is computationally

simple algorithm for spectral analysis. Fourier analysis decomposes a time-series into

global sinusoidal components with fixed amplitudes. As a result, Fourier spectral

analysis is limited to systems that are linear and strictly stationary. However, a high

number of variables affect the heart rate, which suggests that an RRi signal is highly

nonlinear and non-stationary. Most research pertaining to HRV utilizes the FFT to

describe to the spectral components of the RRi signal. We will utilize the FFT spectral

analysis method to obtain results, which can be compared with past and future research.

The FFT is computed using the internal FFT function provided by MATLAB.

Since the FFT method is commonly used in signal processing, we will not discuss this









method in detail, rather we will explain how we computed the transform. All RRi

windows for the RRi dataset (developmental, prestimulus, stimulus, and poststimulus) are

analyzed separately based on a N-point FFT, where N is equal to the total number of

point in the intervals. Often FFT windows are padded with zeros or with repeated copies

of the window. Neither of these techniques adds any information to the resulting FFT.

Padding only serves to increase the resolution of the transformation. Setting N to the

length to the window abolishes the need for padding. Each of the aforementioned

windows is extracted from the original recording based on cumulative time, and the

resulting windows are unevenly sampled signals. The cubic spline interpolation

technique, described in Chapter 5, is applied before computing the FFT and resampling

the interpolated signal at 4Hz. We chose 4Hz to prevent aliasing above of the Nyquist

frequency (2Hz) from falling into the frequency bands of interest. Finally, the resampled

signal is passed through the FFT algorithm and the DC component of the FFT

transformation vector is removed by dropping the first term. Results from this algorithm

can be found in the following chapter.

Lomb Spectral Analysis

In order to avoid the consequence of interpolation, we explore the Lomb method to

estimate the power spectral density (PSD) of unevenly sampled signals. Laguna et al.

concluded that for PSD of these unevenly sampled signals the Lomb method is more

suitable than FFT estimates requiring linear or cubic interpolation [30]. The Lomb

estimate avoids the low pass effect of resampling and avoids aliasing up to the mean

Nyquist frequency:

1
Nf =(2 mRR)









The Lomb method is based on the minimization of the squared differences between

the projection of the signal onto the basis function and the signal under study [30].

X(t )+ 6, = c(i)b,(t )

Minimizing the variance of sE (mean squared error), is synonomous to the minimization

of

N
L x (t,) c (i)b,(t, 2
n=1

which results in [30]

1N
cO(i-= )x(t,)b:(t,) k= b,(tn2
kn n=l

In addition, the signal power at index i of the transformation as described by Lomb

[31]:

N
PJ (i) =c(i)f x (t )b (t= k c(i) 2 22 where 4(i)= c(i)k
n=l

Setting the basis function to that of the Fourier transform

b, (t )= e 2'tn

Then,


k = e2ftn 2 = N
n=1

And

P(i) = PA)= 2f)

Lomb presented a simplification that permits one to solve for each coefficient

independently, known as the Lomb normalized periodogram. First, we solve for the

mean and variance [31]:









x= mRR
o 2 = SDRR)2

Then, we generate the Lomb normalized periodogram using the ith angular frequency

IN -2 )coN 2)sin(ot
(x(t) x)cos(o (t r,)) Y (x(t- x)sin (a), -r)
P(, = 2cos2 ( )) sin2, (n -))

n=l n=l

where r, is defined as,



Y1 sin(2 at)
r =- arctan 1-I
20, zcos(2atn)
Sn=l

MATLAB code written by David Glover is used to compute this Lomb normalized

periodogram [32]. We apply this method of spectral analysis to the identical windows

used to compute the Fourier PSD; and we will compare the results in the following

chapter.

Empirical Mode Decomposition

Approximations are the key to understanding complex natural phenomena and

often there is a significant advantage to these methods. However, science is continuing to

develop ways to describing complex system and behaviors. Why approximate the RRi

signal as a linear, stationary signal when tools exist to analyze nonstationary, nonlinear

systems? The answer is rhetorical because the only reason describes the simplicity of the

computation of the FFT. Given all the assumptions needed to compute the FFT, it is

necessary to search for a more reliable spectral transformation. We have already

observed the Lomb method, which still assumes the signal is linear and stationary. Now,









using a technique developed by Haung et al. [33], this thesis will compare the

effectiveness of the Empirical Mode Decomposition (EMD) as an analysis tool of HRV.

The goal of EMD is to decompose a time series into superposition of Intrinsic

Mode Functions (IMF). EMD uses a shifting process to produce IMFs, which enables

complicated data to be reduced into a form with defined instantaneous frequencies. The

IMFs has the appearance of a generalized Fourier analysis, but with variable amplitudes

and frequencies. Haung et al. describes EMD as the first local and adaptive method in

frequency time analysis. Furthermore, Haung suggests that the advantage of EMD and

Hilbert spectral analysis is that it gives the best-fit local sine or cosine form to the local

data. The frequency resolution for any point is uniformly defined by the local phase

method and is especially effective in extracting low frequencies oscillations [33]. Since

the HRV is characterized by low frequency oscillations, this method should produce

informative frequency domain representation of HRV.

Computation of the EMD algorithm is done in MATLAB using modified code

originally written by G. Rilling [34]. The EMD algorithm [33] coded by Rilling begins

with the sifting process. First, all local extrema of the original signal are identified. The

local maxima are then connected via a cubic spline and the procedure is repeated for the

local minima, respectively forming the upper and lower envelopes. Then, we find the

first component of the sifting process [33], h1.

h, = X(t) = RRi(t)
h, = h0 m,

where m, is the mean of the envelope

The shift process is repeated k times, as follows, until h,, is an IMF (c,).









hik = hl(k-) k
C, = hlk

Then, a residual is found

r = X(t) -c,

Finally the IMF procedure is repeated n times until the residual, r, is a trend or constant.

rn = h, ) c(n+ )

It is then necessary to determine if all the calculated IMF are actually orthogonal to one

another, as the theory suggests. Huang suggests using the follow equation as a check for

orthogonality [33]:


C +C2
10fg -c c2+c2 0 for all IMFs
f g

In order to visualize the result of EMD, we must find the Hilbert transforms of the

IMF components and represent the RRi signal in terms of time, frequency, and amplitude.

We compute the Hilbert transforms using the Cauchy principal value integral to find the

magnitude, a (t), and phase, (t), and instantaneous frequencies, CO [33].



j t td

Z (t) = X (t)+ iY, (t)= a (t)e' (t)

a, (t) = 2X (t) + Y, (t)

Yc (t)
,j (t)= arctan YJ


ast

Finally, we can express the RRi dataset as follows










RRi(t) = a,(t)exp if o, (t\t
j=1

Cleary, the equation above represents the generalized Fourier expansion, but the

amplitude and instantaneous frequency are variable with respect to time [33]. At this

point, we will use the EMD method to provide us with a visual representation of the

power spectrum changes over a prestimulus -* stimulus -* poststimulus window. In

accordance with the description of a COR, detailed in Chapter 2, we should see

significant shift in the magnitude of frequency spectrum over this window.

Statistical Analysis of Parameters

Statistical analysis of the parameters extracted above is performed using a

technique for longitudinal data analysis, called mixed general liner models (MixMod) or

random regression models [35]. The MixMod model allows the concurrent identification

of both group and individual patterns in longitudinal sets with covariates that vary over

time, inconsistently timed data, irregularly timed data, and randomly missed values.

Given the vulnerable state of the infants in the NICU, we must perform test sessions in

accordance with the health status of the infant. These stability considerations can result

in the testing schedules varying from subject to subject and session to session.

Fortunately, problems arising from this effect are alleviated using the MixMod analysis.

In addition, subjects are often transferred or released before reaching 34 weeks GA. The

MixMod model allows us to include the data from these subjects, thereby reducing the

amount of data, which would otherwise be lost. Strengths of the MixMod model include

the acceptance of multiple characteristics and the exploration of patterns across subjects

and groups. Using multiple characteristics, MixMod leaves the statistical correlation and

inclusive of parameters up to the investigator [35]. All MixMod analysis is performed by






58


Dr. Kruger using readily available SAS Proc Mixed software. The results of the MixMod

analysis can be found in the Clinical Results section of Chapter 7.















CHAPTER 7
DISCUSSION AND RESULTS

Engineering Results

We will explore the results of the signal-processing phase of this thesis in

accordance with ECG recordings from one subject. The example subject is not intended

to prove the clinical significance of the hypothesis described in Chapter 3. For clinical

evaluation of aforementioned hypotheses, please refer to the "Clinical Results"

subsection of this chapter.

QRS Detection Algorithm

The following table is a representation of the performance of the QRS detection

algorithm, detailed in Chapter 5. The percent of QRS complexes is computed, as follows,

D det sections
QR detected (# det ec on
Sdetectnons missed

where the interpolated detections are tagged as missed detections and are stored in a

separate array to determine the QRS detection ratio defined above. The interpolations are

a result of no detected QRS complexes falling within the user defined expectancy range,

described in Chapter 5.

Table 7-1 clearly illustrates that on average the algorithm only misses 0.85% of the

QRS complexes in the signal. This detection algorithm is a robust system, which can

effectively detected QRS complexes in the presence of noise and muscle artifact and

produce reliable RRi datasets. Beyond the detection process, we improve quality of the










RRi signal through the implementation of the outliers removal algorithm designed in

Chapter 5.

Table 7-1. Shows the performance of the QRS detection algorithm


Patient


Recording Number QRS detections (%)


3 12130842 99.25
3 12200930 98.8
3 12130842 98.65
3 12060959 93.73
3 12020916 100
3 11221027 99.94
3 11151124 95.11
3 11120913 100
5 2191432 99.11
5 1311022 97.94
5 2071118 99.86
5 2160627 99.36
5 1241004 98.09
5 1211028 99.83
6 3110940 99.88
6 2281027 99.77
6 2211035 99.78
6 2160711 99.94
6 2050923 99.94
6 2071032 99.24
7 3281006 99.67
7 3210950 99.91
7 3141005 100
7 2281002 99.86
7 3071003 99.61
7 2190951 99.92
7 2211017 99.63
Average: 99.14148148


HRV Plots and Statistics

Scientific findings are based on equally important qualitative and quantitative

measures. Quantitative methods observe distinguishing characteristics and elemental

properties resulting in numerical comparative analyses and statistical analyses. For each

recording, an excel file is created, which contains all the pertinent quantitative measures


for each window:









* mRR

* medRR

* SDRR

* Number of sustained increase and decreases in HR derived through Quadrant
Analysis

* The length, width, and (W/L) ratio of the PCA fit elliptical representation of the
Poincare plot

* ULF, VLF, LF, and HF power of the Lomb PSD

* TP and SVrao of Lomb PSD

* ULF, VLF, LF, and HF power of the Fourier PSD

* TP and SVroao of Fourier PSD

These quantitative measures are used for MixMod analysis in the "Clinical Results"

section.

Qualitative measures are those associated with interpretative approaches, from the

investigators point of view, rather than discrete observations. While qualitative

methodologies can lead to biased results, there is no doubt that they should be included in

medical research. The power of visualization and quick interpretation of data possessed

by humans cannot be replaced by quantitative measures. Therefore, for each recording,

six plots are created for qualitative analysis:

* HR with outlier removal

* Quadrant analysis

* Poincare

* Fourier PSD

* Lomb PSD

* Hilbert-Haung spectrum of EMD







62


HR Plot



160

A 150


0 50 100 150 200
160

B 2 150-

140

0 5 10 15 20 25 30 35 40 45
160 ..-

) 150-
C)a IL
140

0 5 10 15 20 25 30 35 40


D) i 150-
140

0 5 10 15 20 25 30 35 40 45
Time (s)


Figure 7-1. HR plots of example RRi data. From top to bottom they represent: the (a)
developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus periods.




The HR plot shown above illustrates the original RRi (blue) waveform displayed

as instantaneous heart rates and the resulting RRi (red) after outliers have been removed.

The figure above is divided into four plots representing the four windows of analysis.

Plot (b) illustrates the removal of an outlier, which causes a high frequency spike in the

signal. Plot (a) demonstrated the algorithm efficiency of preserving trends in RR

intervals. Figure 7-1 is used to visualize cardiac accelerations and decelerations,

quantitatively described through Quadrant Analysis.








63



Quadrant Plot


The Quadrant plot shown in Figure 7-2, below, is organized in the same manner as


Figure 7-1 with each plot representing a particular analysis window. Quadrant Analysis


is a more useful tool for quantitative analysis, but the visual representation can provide


insight into the amount of high and low frequency variability. The excel file described


above contains some useful information regarding this plot that will later be used for


MixMod analysis.


0.





-0.

0.
s,

A) s

a,



-0.

0.
B) -


-0.

0.


C) S
-

-D -0.


Quadrant Analysis
05 I, I, ,I,-,-,
05



0
05 L L L L L---
-0.05 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 0.C
05


0


05 L L -
-0.05 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 0.C
05





05
0.05 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 0.C
-0.05 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 0.C


0.05
-)

D) S 0
0
"o -0.05------


0
deltRR(n)


Figure 7-2. Quadrant Analysis plots of example RRi data. From top to bottom they
represent: the (a) developmental, (b) prestimulus, (c) stimulus, and (d)
poststimulus periods.











Poincare Plot

The Poincare plot is a valuable tool for both quantitative and qualitative analysis.


The Poincare method of RRi interpretation is observed to prove normal cardiac function


and analytically evaluated by the fitted ellipse to described short and long-term HRV.


0.46
u 0.44
0.42
A) +
S0.4
of 0.38
0.36
0.
0.46
0.44
B) 0.42
S 0.4
rO
of 0.38
0.36
0.
0.46
0.44
C) 0.42
S 0.4
0.38
0.36
0.
0.46
0.44
D) 0.42
S 0.4
of 0.38
0.36
0.


Poincare Plot

+-



36 0.37 0.38 0.39 0.4 0.41 0.42 0.43 0.44 0.45 0.'






36 0.37 0.38 0.39 0.4 0.41 0.42 0.43 0.44 0.45 0.,


36



36


0.37 0.38 0.39 0.4 0.41
RR(n)(s)


0.42 0.43 0.44 0.45 0.46


Figure 7-3. Poincare plots of example RRi data. From top to bottom they represent: the
(a) developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus periods.



Again the magnitude corresponding of a given point is represented by the basic


colors of the visible light spectrum with red being the next to least frequent (black is the


least) and violet being the most frequent. Plot (a) represents the dataset with the most


long-term variability or low-frequency changes inferred by the length of the major axis of


the ellipse. This is expected because plot (a) is a 240s window, while plots (a,b,c) are 45s










windows. We can also see that during the stimulus period the major and minor axis of

the ellipse become significantly longer, perhaps signifying an increase the LF power, HF

power, and total power of the HRV spectrum. We will next look at the Fourier and Lomb

PSD to confirm this hypothesis.

Fourier and Lomb Plots

The Fourier and Lomb PSD plots essentially describe the same information;

therefore we will observe them together.


0.03

A) 0.02
0.01
0

0.06
B) 0.04
0.02
0

0.06
C) 0.04
0.02


FFT Power Spectrum

--



0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1





0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
I-


0.(
D) 0.(
0.(


0.5
Frequency (Hz)


Figure 7-4. Fourier PSD plots of example RRi data. From top to bottom they represent:
the (a) developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus
periods.












Lomb Power Spectrum
I I I


100-

A) 50


0 1 I -/L I I I I I I I I
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0

0

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
5


10-

5

0
0 0.1 0.2


0.4 0.5 0.6 0.7 0.8 0.9 1


t 0.5
Frequency (Hz)


Figure 7-5. Lomb PSD plots of example RRi data. From top to bottom they represent:
the (a) developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus
periods.



First, we will address the observational hypothesis from the Poincare Plot section.


Inspecting the Fourier and Lomb stimulus plots, we note that the area under the curve


increases from the prestimulus period to the stimulus period, which means there is an


increase in total power. Also, power increases in the LF and HF bands correspond to


increases in the length of the major and minor axes of the ellipse, respectively. In


conclusion, through qualitative analysis, we have shown that the Poincare plot is a time


domain representation of the PSD.


Now we will compare the PSD representations of the Fourier and Lomb


Spectrum. The Fourier and Lomb plots are comparable over the developmental window,


1









but the Fourier plot is choppier than the Lomb plot over the windows corresponding to

the COR. The choppy Fourier plot is a result of the low frequency resolution of the N-

point Fourier transformation. Therefore, the frequency resolution is lower for the short

windows analyzed in plots b, c, and d. However, the Lomb computation allows the user

to set the frequency resolution. For the purpose of obtaining a smooth transformation, we

have set the Lomb PSD resolution to 0.001Hz. To increase the frequency resolution of

the FFT, we must pad the data before performing the Fourier transform. For reasons

discussed in Chapter 6, we have refrained from doing so. The apparent equivalent Lomb

and Fourier PSD plots for the developmental window suggest that the Lomb method is at

least as successful as the Fourier for PSD estimation. For clinical interpretation, we will

use the Lomb Spectrum when computing the cumulative power in the frequency bands of

interest.

EMD Plots

At this point, we use the Hilbert-Haung Spectrum (HHS) only as a qualitative

analysis tool. If the clinical findings indicate the usefulness of this time-frequency

domain representation of HRV, we will expand to use EMD as a quantitative analysis

tool. Figure 7-6 is the HHS of the example RRi dataset used in this section. Black lines

separate the prestimulus, stimulus, and poststimulus periods. We can relate the HHS

power shift of the stimulus period to those described by the Poincare and PSD

assessments of HRV. All four representations are in agreement that during the stimulus

period the LF power, HF power, and TP increase. However, the HHS shows us that this

increase occurs just after the rhyme begins only to tail-off significantly shortly thereafter.

Initially the EMD seems to be a robust method to describe the HRV characteristics of the

COR with respect to time. However, since this method has not often been used for









clinical interpretation, the inclusion of the method in the analysis of the hypothesized

COR will be left up to the PI.


6b
Time


Figure 7-6. Hilbert-Haung spectrum of the COR windows of an example RRi dataset.

Clinical Results

At the time of submission, we only have four complete datasets (recordings for the

28-34 weeks GA) and one incomplete dataset (recordings for 28-33 weeks GA) that can

be used for clinical interpretation. The current sample size (n=5) is well below the

minimum sample size (n=28) described in the design of this study [7]. Therefore, we will

not be able to perform MixMod analysis of the data at this time. Instead, we will use a

repeated measures ANOVA analysis of the PSD. Repeated measure ANOVA allows us

to compare the variance caused by the independent variables to a more accurate error

term, which includes the variance effects of individual subjects. Thus, this method










requires fewer participants to yield an adequate interpretation of the data than the

MixMod method.

Early analysis based on the Poincare plots, Quadrant Analysis, and EMD will not

be performed due to the lack of comparable literature. We will exploit the

aforementioned methods upon completion of the study, so the PSD results can be used to

confirm the clinical significance of these models. The repeated measure ANOVA model

returned only two independent variables (HF power and SV power), which passed the

equal variance test. We will first examine results of the repeated measures ANOVA

analysis of the SV power defined in Chapter 6.

Table 7-2. Results of the repeated measures ANOVA with independent variable (SV
power).


Session n p a
1 5 0.807 0.0857
2 5 0.729 0.089
3 5 0.846 0.121
4 5 0.878 0.0501
5 5 0.742 0.239
6 5 0.866 0.0836
7 4 0.613 0.186


SV Power

1.2


S0.8 --\-

0.6
| 0.4
0.2
0 ------------------------
1 2 3 4 5 6 7
Session

Figure 7-7. Plot of mean SV power with standard deviation bars across sessions










Now we examine repeated measures ANOVA analysis of the HF power, also defined in

Chapter 6.

Table 7-3. Results of the repeated measures ANOVA with independent variable (HF
power).


Session n p a
1 5 68.2 25.925
2 5 133.6 144.921
3 5 48.46 26.837
4 5 41.32 23.877
5 5 178.1 172.465
6 5 47.22 58.952
7 4 150.75 126.818


HF Power


400
350
300
250
200
150
100
50
0
-50


S7
1) 13 A r_ n


Session


Figure 7-8. Plot of mean HF power with standard deviation bars across sessions

While Figure 7-7 shows a possible non-linear trend with respect to GA, Figure 7-

8 does not elicit any significant trend in time. First, the small sample size doesn't allow

for the development of a trend. Secondly, there exist many exogenous variables, which

conceal a developing trend corresponding to GA. These variables include, but are not

limited to:


\J r-









* Drug interventions and interactions

* Respiratory support interventions

* Physiological discrepancies in ANS maturation

* Changes in NICU environment

* Incorrect evaluations in GA

Patient logs will be analyzed following the completion of the study for drug

interventions. The PI will attempt to parallel the interventions with the recordings that

contribute most to the standard deviation of each variable. It is known that certain drugs

inhibit or even block parasympathetic and sympathetic responses. Furthermore, many of

our subjects are on multiple drugs whose interactions with respect to the ANS are not

known.

Many subjects intermittently require respiratory interventions such as continuous

positive airway pressure (CPAP). A closed-loop study showed that adequate application

of long-term CPAP therapy produced substantial alterations in the major physiological

mechanisms that influence HRV and blood pressure variability [36]. Based upon

preliminary results, see Figures 7-7 and 7-8, final results need to include an analysis of

how drug and respiratory interventions hinder the ability to identify trends in the ANS

landscape. The significance of all exogenous variables mentioned above will emerge as

the number of subject in the dataset increases. Unfortunately, these results are not

included in this thesis, but we refer you to further publications by the PI and the author of

this paper. For a complete time and frequency domain analysis of a randomly selected

subject, please refer to Appendix A and B, respectively.














CHAPTER 8
CONCLUSION

We have successfully modified code written by HP to communicate with the

Agilent/HP NICU monitor. This program records the ECG, respiratory rate, and rhyme

state in a text file at a 500 Hz sampling rate. The text file can be read into MATLAB for

QRS detection and HRV analysis. In MATLAB, we have created an algorithm to

perform accurate detection of QRS complexes. The average percentage of missed

detections is approximately 0.8%, which we feel represents a near optimal detection

algorithm. For those QRS complexes that are undetectable, we have constructed an

interpolation method to replace the missing data. Using our novel method of tagging

outliers, we are able to remove outliers, but preserve trends across time.

We have found that the most useful measurement of HRV in the time-domain is

through the analysis of Poincare plots. Fitting an ellipse to the Poincare representation of

the RRi is usefully in relating time-domain analysis to the frequency domain. Currently,

most spectral representations described in literature use the Fourier transform. We have

found that the Lomb normalized periodogram is a better representation of the power

spectrum of RRi waveforms. The Lomb algorithm avoids the negative effects of

interpolation and resampling required for FFT computation. Also, Empirical Mode

Decomposition is a favorable method for measuring changes in the frequency spectrum

across time. The EMD spectrum preserves more information regarding shifts in the

frequency spectrum during the stimulus period than either the Lomb or Fourier spectrum.

EMD has been proven to be a useful analysis tool to depict a COR.






73


We have discussed the importance of both qualitative and quantitative measures

of developmental changes in HRV and in identifying a COR; we have extracted

numerous features using time and frequency domain HRV analyses, which will be used

to analyze and test the hypotheses of the clinical study still in progress. We feel we have

completed the design phase of the study and need only to increase our sample size. As

the sample size increases, the clinical significance of each feature will be evaluated and

we will move to interpret our findings and test our hypotheses.














CHAPTER 9
FUTURE WORK

A need exists in the NICU and beyond to determine the state of development of the

ANS. Before an infant is discharged from the NICU, it is important to determine their

ability to perform certain life support tasks. In Chapter 2, these life-supporting tasks

were linked to the maturity of the ANS, specifically the PNS. The list below is compiled

from the 1998 suggested guidelines for neonatal discharge established by the American

Academy of Pediatrics (AAP) [37]. In the judgment of the responsible physician there

has been:

* A sustained pattern of weight gain of sufficient duration.

* Adequate maintenance of normal body temperature with the infant fully clothed in
an open bed with normal ambient temperature (240C to 250C).

* Competent suckle feeding, breast or bottle, without cardiorespiratory compromise.

* Physiologically mature and stable cardiorespiratory function of sufficient duration.

* Appropriate immunizations have been administered.

* Appropriate metabolic screening has been performed.

* Hematologic status has been assessed and appropriate therapy instituted as
indicated.

* Nutritional risks have been assessed and therapy and dietary modification instituted
as indicated.

* Sensorineural assessments, hearing and funduscopy, have been completed as
indicated.

* Review of hospital course has been completed, unresolved medical problems
identified, and plans for treatment instituted as indicated.









The AAP has established these guidelines, which include subjective observations of

cardiorespiratory function and sensorineural assessments. It is our hypothesis that

cardiac function, as well as, neural assessments can be objectively realized using an index

of ANS maturation. We will create this index of ANS functionality based on HRV

comparisons with mature young adults. Below, we suggest a pilot study, which shall

provide evidence that an index of maturation can be realized.

Record Baseline Adult HRV Data

In order to establish a baseline for a mature ANS with a synonymously high vagal

tone, it is important to analyze the HRV of the young, healthy adult. ECG recordings

will be obtained for a sample of 28 male and female adults meeting the following criteria:

* Between the ages of 18-30

* No history of congestive heart failure, myocardial infarction, coronary artery
disease, head trauma, or diabetes

* Normal sinus ECG rhythm


The length of the recordings should be equal to those used in the study of the fetus. This

is important we comparing both time and frequency domain features of HRV and the

COR.

Each subject should be oriented in the supine position to ensure consistency with

ECG recordings of the neonates. Before recording, all subjects should remain in this

position without auditory and visual simulation, until they are verified to be in REM or

paradoxical sleep. Consistency with the neonatal recording described in the paper is a

significant issue for accurate interpretation of HRV data and the adult period of REM

sleep corresponds to the neonatal state of "active sleep" [38]. Siegel suggests that the

phasic motor activation that characterizes REM sleep in the adult is also present in the









neonate. The amplitude of the twitching, described during REM sleep, is more intense in

neonates than in the adult, but a developmental continuity can be observed between

"active sleep" in the neonate and REM sleep in the adult.

Both active sleep in the neonate and REM sleep in the adult can be defined by

purely behavioral criteria [38]. Therefore, the recordings should be taken in the presence

of a medical professional qualified to identify phases of paradoxical sleep. An auditory

stimulation that is widely familiar and contains no political or nationalistic propaganda

should be used to prevent to prevent undesirable neural responses due to stress. The ECG

recordings during the auditory stimulation can later be used to compare the COR of the

neonate and adult.

Creating the Index of ANS Maturation

We suggest using the EMD method, detailed in Chapter 6, to compare the COR of

the HRV groupings described in the previous subsection with the mature adult. The

EMD data provides a time-based representation of the COR. The examiner should

compare differences in the COR of the adult and neonate to determine an efficient

method for relating these changes. We suggest exploring the use of power-ratios, neural

networks, windowing functions, etc. as possible criteria to represent the EMD data in the

form of an index. In addition, the development data can be compared using PCA and a

neural network to classify the Lomb spectrum into different level of ANS maturation.






















APPENDIX A

TIME-DOMAIN HRV ANALYSIS


Week One




155 --


3 25 30 35 40
Time (s)


( 0 4

0.35
0,37 0-38 0.39


0.37 0.38 0.39
0,44
S0.42 -
S 0.4
n 0.38
0 37 038 0.39


Poincare Plot



* *


0.4 0.41 042 0.43 0.44


0.4 0.41 042 0.43 0.44






0.4 0.41 042 0.43 0.44
RR(nXs)


5 10








78



Week Two


0 50 100


150


200


0 5 10 15 20 25 30 35 40 45
150

5 145 -

140',
A A~ /


Time (s)


Poincare Plot
0.46 -
0.44 *
0.42- -
S0.4 -c- -
0.38
0.38 0.39 0.4 0.41 0.42 0.43 0.44 0.45 0.46 0.47

0.46 -
0.44 -
0 44
S0.42 --
o 0.4 -.----
0.38
0.38 039 0.4 0.41 0.42 0.43 0.44 0.45 046 047

0.46 -
0.44- -
0.42 -* t
o4 04
0.38
0.38 039 0.4 0.41 0.42 0.43 0.44 0.45 046 047

S0.46 -
0.44-2
S0.42 t
W 0.4
0.38
0.38 0.39 0.4 0.41 0.42 0.43 0.44 0.45 0.46 0.47
RR(n)(s)








79



Week Three


0 5 10 15 20 25 30 35 40


0 5 10 15 20 25 30 35 40


0 5 10 15 20 25 30 35 40 45
Time (s)


Poincare Plot







0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38







0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38







0.31 0.32 0.33 0.34 0.35 0.36 0.37 038


0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38
RR(n)(s)


0.38
. 0.36
+ 0.34
j 0.32
0.3
0.3
0.38
.. 0.36
+ 0.34
f 0.32
0.3
0.3
0.38
. 0.36
S0.34
r 0.32
0.3
0.3
0.38
. 0.36
0.34
of 0.32
0.3
0.3


-





-








80



Week Four


50 100 150 200


182
180
I 178
- 176
m 174-
172
0
182
180
r 178
0 176
m 174
172
0
182
180
r 178
a 176
m 174-
172
0
182
180
1 178
a0 176
m 174-
172
0





0.38
0.36
S0.34
"f 0.32
0.3
0.3
0.38
S0.36
0.34
0.32
0.3
0.3
0.38
0.36
S0.34
If 0.32
0.3
0.3
0.38
." 0.36
g 0.36-
0.34
"i 0.32
0.3
0.3


5 10 15 20 25 30 35 40


5 10 15 20 25 30 35 40
Time (s)


Poincare Plot







0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38


0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38

0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38


0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38
RR(n)(s)


5 10 15 20 25 30 35 40








81



Week Five


180

S170
m 160
0 160


170
0..
m 160

150

180

170
0-..
m 160

150

180

170

m 160


Poincare Plot
0.4
0.38- .
+ 0.36 t-. *
S0.34 -
0.32
0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.4
0.4 p l?.ii
S0.38-
0.36 "
' 0.34-
0.32
0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.4
0.4 pp;
S0.38-
S0.36 -
O 0.34
0.32 I I
0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.4
0.4 p

X 0.38

0.34
0.32
0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.4
RR(n)(s)








82



Week Six




180

S170 r s

160

0 50 100 150 200
180 Tim

170

160

0 5 10 15 20 25 30 35 40
180

S170

160

0 5 10 15 20 25 30 35 40
180

170-

160

0 5 10 15 20 25 30 35 40




Poincare Plot
0.38 -

0.36 *-

" 0.34 -

0.32
0.32 0.33 0.34 0.35 0.36 0.37 0.38
0.38

0.36 -

S0.34

0.32I
0.32 0.33 0.34 0.35 0.36 0.37 0.38
0.38 PP'

0.36-

0.34-

0 .3 2 L L
0.32 033 0.34 0.35 0.36 0.37 038
0.38

0.36

' 0.34

0.32 I
0.3 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.4
RR(n)(s)








83



Week Seven


1I-(U- Ir -' I

2 160o

150

0 50 100 150 200
170- i ,(,

160 55 A

150

0 5 10 15 20 25 30 35 40 45


0 5 10 15 20 25 30 35 40
170 Tima


0 5 10 15 20 25
Time (s)


30 35 40 45


Poincare Plot


--0.42
0.4
S0.38
a 0.36
aO 0.34
0.32
0.:

S0.42
0.4
0.38
" 0.36
n, 0.34
0.32
0.:

-0.42
0.4
0.38
a 0.36
Of 0.34


w 0.42
S0.4
+0.38
- 0.36
rO 0.34
0.32
0.32


0.34 0.36


0.38
RR(n)(s)


0.42 0.44



















APPENDIX B
FREQUENCY-DOMAIN HRV ANALYSIS


Week One


Lomb Power Spectrum
40 II------ I I -------


40
20


0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
20


10


0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
60

40

20 -

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9


Frequency (Hz)


0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9













Week Two



Lomb Power Spectrum
100


50



0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1


30

20

10

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9


20 -

10

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Frequency (Hz)



Week Three




150



50









S 0.2 0.3 0.4 0.5 0. 0.7 08 0.9
60 r r i I i



20

Cl 02 03 0C4 Cb Cr3 06 0 C.S '


0.1 0.2 03 0.4 0.5 0.6 0.7 08 0,9 1
Frequency (Hz)








86



Week Four



1501





0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 08 09 1
30

20


I0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
40



20

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
30


20 I I
20


10
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Frequency (Hz)


Week Five



150

100

50

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
40 1 1



20__ I

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
30 ,

20 -

10

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
40




0 Ii
20

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Frequency (Hz)








87



Week Six



60 Lqmb Powr Spoectrurm-n


-' I ,


0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
20"





0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
20





0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
20 I I




0. --
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Frequency (Hz)



Week Seven


Lomb Power Spectrum
100




0 ^ ^ I I-- ,4,^ I A- I--.--------------
0 0.1 02 03 0.4 0.5 0.6 0.7 08 09 1
15

10

5

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1




10




20

10
30 i 0I 0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
30 -----------------------------------


Frequency (Hz)















LIST OF REFERENCES


1. BLISS, Causes of premature birth (2003),
http://www.bliss.org.uk/support/causes.asp, Accessed 4/4/2003.

2. Thomas, K., Causes of premature birth (2001),
http://www.kerri.thomas.btintemet.co.uk/causes.htm, Accessed 4/4/2003.

3. Chatow, U., Davidson, S., Reichman, B. & Akselrod, S. (1995). Development and
maturation of theautonomic nervous system in premature and term infants using
spectral analysis of heart rate fluctuations. Pediatric Research, 37, 294-302.

4. Stark, R., Myers, M., Daniel, S.,Garland, M., & Kim, Y. (1999). Gestational age
related changes in cardiac dynamics of the fetal baboon. Early Human
Development, 53, 219-237.

5. Boon, R., Learning Discoveries Psychological Services (2002),
http://home.iprimus.com.au/rboon/HeartRhythmsandHRV.htm, Accessed
1/12/2003.

6. Groome, L., Mooney, D., Holland, S., Smith, L., Atterbury, J., & Loizou, P.
(1999). High Vagal Tone is Associated with More Effiecient Regulation of
Homeostasis in Low-Risk Human Fetuses. Dev. Psychophysiology, 36, 25-34.

7. Krueger, C. (2002). Fetal heart rate variability and learning in the 28-34 week fetus.
(Doctoral dissertation, University of North Carolina at Chapel Hill, 2001).
Dissertation Abstracts International, 62, 11B.

8. Porges, S. (1995). Orienting in a defensive world: Mammalian modifications of
our evolutionary heritage. A Polyvagal Theory. Psychophysiology., 32, 301-318.

9. Porges, S., Doussard-Roosevelt, J., Stifter, C., McClenny, B., & Riniolo, T.
(1999). Sleep state vagal regulation of heart period patterns in the human newborn:
An extension of the polyvagal theory. Psychophysiology, 36, 14-21.

10. Beauchaine, T. (2001). Vagal tone, development, and Gray's motivational theory:
Toward an integrated model of autonomic nervous system functioning in
psychopathology. Development and Psychopathology, 13, 183-214.

11. Gagnon, R. (1992). Fetal behavioral states: Biological alteration. Seminars in
Perinatology, 16, 234-238.




Full Text

PAGE 1

ANALYSIS OF NEONATAL HEART RATE VARIABILITY AND CARDIAC ORIENTING RESPONSES By DEVIN N. LEBRUN A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF ENGINEERING UNIVERSITY OF FLORIDA 2003

PAGE 2

Copyright 2003 by Devin N. Lebrun

PAGE 3

Dedicated to my father, Denis Lebrun, for his love and support.

PAGE 4

ACKNOWLEDGMENTS I would like to thank my advisor, Dr. Johannes van Oostrom, for all the advice he has given me on realizing my career goals. I knew what I wanted, but he helped me get from there to here and now onto somewhere great. He endured barrages of scatterbrained, unintelligent questions and always surfaced with the answers I needed. I am grateful to Dr. Charlene Krueger for allowing me to share in her research and express my gratitude to Dr. Krueger, Dr. van Oostrom, and the College of Nursing for providing me with funding throughout my research. I would also like to thank Kelly Spaulding, John Hardcastle, Dr. JS Gravenstein, Dr. Richard Melker, and Dr. John Harris for their support in various areas. Words cannot describe my appreciation for my parents. They have taught me to be determined, self-motivated, and not afraid of failure. I thank Linnea, for her companionship and assistance throughout the past 6 years. Finally, thanks go to all my friends for keeping me from doing something stupid (well, more stupid). iv

PAGE 5

TABLE OF CONTENTS Page ACKNOWLEDGMENTS................................................................................................IV LIST OF TABLES..........................................................................................................VIII LIST OF FIGURES..........................................................................................................IX ABSTRACT......................................................................................................................XI CHAPTER 1 INTRODUCTION........................................................................................................1 2 AUTONOMIC NERVOUS SYSTEM AND HEART RATE VARIABILITY...........4 Autonomic Nervous System.........................................................................................4 Developmental Changes...............................................................................................6 Vagal Tone and Homeostasis.......................................................................................7 Cardiac Orienting Response.........................................................................................8 Respiratory Sinus Arrhythmia......................................................................................8 3 PRELIMINARY RESULTS AND STUDY PROCEDURES....................................10 Preliminary Findings..................................................................................................10 Proposed Study...........................................................................................................11 Recitation Sessions..............................................................................................12 Test Sessions.......................................................................................................13 4 TESTING....................................................................................................................16 5 QRS DETECTION.....................................................................................................21 Noise/Baseline Wander..............................................................................................23 Slope Detection...........................................................................................................27 Outliers.......................................................................................................................32 Direct RRi Analysis.............................................................................................34 Slope Analysis.....................................................................................................36 Differential Analysis...........................................................................................37 Interpolation................................................................................................................39 v

PAGE 6

6 STATISTICAL ANALYSIS......................................................................................44 Time Based.................................................................................................................44 First Order Statistics............................................................................................44 Poincar Plots......................................................................................................45 Quadrant Analysis...............................................................................................48 Frequency Based.........................................................................................................50 Fast Fourier Spectral Analysis.............................................................................51 Lomb Spectral Analysis......................................................................................52 Empirical Mode Decomposition..........................................................................54 Statistical Analysis of Parameters..............................................................................57 7 DISCUSSION AND RESULTS.................................................................................59 Engineering Results....................................................................................................59 QRS Detection Algorithm...................................................................................59 HRV Plots and Statistics.....................................................................................60 HR Plot.........................................................................................................62 Quadrant Plot................................................................................................63 Poincare Plot................................................................................................64 Fourier and Lomb Plots................................................................................65 EMD Plots....................................................................................................67 Clinical Results...........................................................................................................68 8 CONCLUSION...........................................................................................................72 9 FUTURE WORK........................................................................................................74 Record Baseline Adult HRV Data.......................................................................75 Creating the Index of ANS Maturation...............................................................76 APPENDIX A TIME-DOMAIN HRV ANALYSIS...........................................................................77 Week One...................................................................................................................77 Week Two...................................................................................................................78 Week Three.................................................................................................................79 Week Four..................................................................................................................80 Week Five...................................................................................................................81 Week Six.....................................................................................................................82 Week Seven................................................................................................................83 B FREQUENCY-DOMAIN HRV ANALYSIS............................................................84 Week One...................................................................................................................84 Week Two...................................................................................................................85 vi

PAGE 7

Week Three.................................................................................................................85 Week Four..................................................................................................................86 Week Five...................................................................................................................86 Week Six.....................................................................................................................87 Week Seven................................................................................................................87 LIST OF REFERENCES...................................................................................................88 BIOGRAPHICAL SKETCH.............................................................................................92 vii

PAGE 8

LIST OF TABLES Table page 3-1 Organization of Recitation and Test Sessions............................................................11 7-1 Shows the performance of the QRS detection algorithm...........................................60 7-2 Results of the repeated measures ANOVA with independent variable (SV power)..69 7-3 Results of the repeated measures ANOVA with independent variable (HF power)..70 viii

PAGE 9

LIST OF FIGURES Figure page 2-1 The divisions of the central nervous system..............................................................4 3-1 Time-line for the windowing of test data.................................................................14 4-1 Example of text file generated from Data Acquisition Program..............................18 4-2 The setup for the test sessions..................................................................................19 5-1 QRS detection algorithm..........................................................................................23 5-2 Magnitude and phase response of the 60Hz notch filter used to filter the raw ECG signal........................................................................................................................25 5-3 Magnitude and phase response of the bandpass filter used to filter the raw ECG signal........................................................................................................................26 5-4 Magnitude and phase response of the differentiator function..................................27 5-6 Example of detected QRS complexes and beat intervals.........................................31 5-7 Histogram of RRi distribution from one subject......................................................33 5-8 Standard deviations of the Normal Distribution......................................................34 5-9 Illustrates the bouncing effect and upper and lower bound symmetry problems with using the direct RRi analysis method for tagging outliers.......................................35 5-10 Illustrates an RRi trend that would not be preserved using the slope analysis method for tagging outliers......................................................................................37 5-11 Outliers that have been tagged, using the differential method, are marked by a green dot. Shows the efficiency of the method to preserve trend and remove bouncing effects.......................................................................................................39 5-12 LaGrange Polynomial Interpolation.........................................................................40 5-13 Linear Interpolation..................................................................................................41 5-14 Cubic Spline Interpolation.......................................................................................41 ix

PAGE 10

5-15 Cubic Interpolation...................................................................................................42 6-1 Poincar plot representation of RRi. Since data points overlap, colors, coordinated with the visible spectrum, are used to indicate magnitude at each point.................47 6-2 Defines the quadrants used in Quadrant Analysis....................................................48 6-3 Example of Quadrant Analysis applied to an RRi dataset.......................................49 7-1 HR plots of example RRi data. From top to bottom they represent: the (a) developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus periods............62 7-2 Quadrant Analysis plots of example RRi data. From top to bottom they represent: the (a) developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus periods.63 7-3 Poincar plots of example RRi data. From top to bottom they represent: the (a) developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus periods............64 7-4 Fourier PSD plots of example RRi data. From top to bottom they represent: the (a) developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus periods............65 7-5 Lomb PSD plots of example RRi data. From top to bottom they represent: the (a) developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus periods............66 7-6 Hilbert-Haung spectrum of the COR windows of an example RRi dataset.............68 7-7 Plot of mean SV power with standard deviation bars across sessions.....................69 7-8 Plot of mean HF power with standard deviation bars across sessions.....................70 x

PAGE 11

Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Engineering ANALYSIS OF NEONATAL HEART RATE VARIABILITY AND CARDIAC ORIENTING RESPONSES By Devin N. Lebrun August 2003 Chair: Johannes H. van Oostrom Major Department: Electrical and Computer Engineering We have applied signal-processing methods to evaluate the clinical significance of heart rate variability (HRV) features, assess developmental changes in HRV, evaluate vagal tone, and demonstrate a cardiac orienting response (COR). The study is based on a sample of 28 low-risk premature newborns randomly assigned to one of two groups. Group 1 is exposed to an auditory stimulation beginning at 28 weeks GA and Group 2 begins exposure at 32 weeks GA. An eight-minute ECG recording is taken once a week from 28-34 weeks GA. The recordings were analyzed in the time and frequency domain using MATLAB. Extracted features are used to classify developmental changes in HRV across individuals and subject groups. In addition, we have used a staggered grouping approach to separate cardiac orienting responses from reflexive responses. At the conclusion of this thesis, the sample size (n=5) was not large enough to elicit any developmental trends in GA. However, we were able to extract useful features to xi

PAGE 12

describe developmental changes in HRV and we have acknowledged Empirical Mode Decomposition as an effective measure of a cardiac orienting response. Also, we have identified numerous exogenous interactions and perturbations which can interfere with the baroreceptor and chemoreceptor responses controlling HRV. These interactions were not considered in the original design of the study but will now be considered during final analysis. xii

PAGE 13

CHAPTER 1 INTRODUCTION The premature infant faces many obstacles beyond the control of present day medical science. Since the loss of the in utero environment strips the premature infant of the most advantageous setting for growth and development, it is increasingly important for medical science to reach a greater understanding of the developmental changes during this final trimester. Survival rates of neonates will continue to increase as science begins to elucidate more regarding this crucial period of development. In general, scientific advancements are predicated upon a divide and conquer approach. It is important to embrace the past discoveries as background for future research. Therefore, this study expands upon a preliminary study conducted by the Principal Investigator (PI), Charlene Krueger, Ph.D, ARNP. The findings of the preliminary study indicate that the analysis of heart rate variability (HRV) can provide significant insight into the development of the premature and fetus autonomic nervous system (ANS). In collaboration with Dr. Krueger, the current study moves from the comforts of the uterus to the premature infants continuing development in a Neonatal Intensive Care Unit (NICU) setting. About a third of premature births occur for no apparent reason and with little or no warning. The list below describes some known causes of premature birth [1,2]: Pre-eclampsia, occurring in about 1 in 14 pregnancies, causes around a third of all premature births. Pregnancies involving multiple fetuses are likely to end early. Stressful events can elicit the mother into a premature labor. 1

PAGE 14

2 An infection involving the uterus may trigger an early delivery. The mother's water may break early, starting the delivery process. A shortage of oxygen exchange to the placenta may cause a doctor to advise a caesarian section. Malnutrition or chronic diseases, including high blood pressure, diabetes, kidney disease and hypothyroidism can cause premature births. Regardless of the cause for an early birth, premature infants are commonly highly dependant upon the continuation of breathing support, intravenous nutrition, intravenous antibiotics, and phototherapy. These infants are frequently placed in NICUs, which specialize in caring for the very young and/or very small. The neonates high dependency upon medical interventions demonstrates a need to establish criteria regarding his or her efficiency at regulating homeostasis. Homeostasis is the maintenance of equilibrium in a biological system by means of automatic mechanisms, which counteract influences trending toward disequilibria. It involves constant internal monitoring and regulating of numerous factors, including oxygen, carbon dioxide, nutrients, hormones, and organic and inorganic substances. Due to homeostatic controls, the concentrations of these substances in body fluid remain unchanged, within limits, despite changes in the external environment. Efficient regulation of these parameters is necessary to support life, on all levels. Therefore, the ability to evaluate homeostatic controls would facilitate the decision process leading to the neonates discharge from the unit. The events surrounding the control of homeostasis will be discussed further in Chapter 2. We hypothesize that an index of HRV will provide sufficient insight into the maturity level of the neonates ANS. In turn, this will provide us with information regarding the ability of the neonate to control events of homeostasis. In a study by

PAGE 15

3 Chatow et al, increasing conceptual age, CA, or gestational age, GA, is correlated with neonatal autonomic control. Chatow used a power ratio, discussed in Chapter 2, which was assumed to reflect the sympathovagal balance of the ANS [3]. In a study of the fetal baboon, Stark et al. suggest that Sequential decreases in fetal heart rate, increases in RR interval (RRi) variability and increases in changes in RRi and RRi with age imply an overall maturation in autonomic cardio-regulatory control processes. Increases with gestation of high frequency components of variability are compatible with enhanced parasympathetic modulation of the fetal heart rate [4]. This thesis will attempt to expand upon these studies to examine developmental changes, cardiac orienting responses, and to introduce a maturation index based on time and frequency-domain analysis of neonatal HRV. Current methods of frequency domain analysis tend to center around one method of analysis. Therefore, new methods of analysis will be constructed as a result of combining the most functional signal processing tools described in recent literature. This thesis will explore and compare a variety of spectral analysis methods and focus on those methods deemed most clinically significant. In addition, we will parallel techniques used to examine the frequency spectrum of HRV with time-domain techniques to provide a complete clinical evaluation of each patient and grouping.

PAGE 16

CHAPTER 2 AUTONOMIC NERVOUS SYSTEM AND HEART RATE VARIABILITY Autonomic Nervous System The ANS is part of the human central nervous system (CNS). The CNS can essentially be divided into three closely interacting sections: the motor nervous system, sensory nervous system, and autonomic nervous system. The ANS originates from the brainstem at the base of the spinal cord. It includes the nerves, which innervate the smooth muscles of the heart, internal organs, and glands. The ANS is subdivided into the sympathetic (SNS) and the parasympathetic (PNS) systems, which originate from different areas in the brainstem and spinal cord. CNS Sensory Nervous System SNS PNS Autononmic Nervous System Motor Nervous System Figure 2-1. The divisions of the central nervous system. 4

PAGE 17

5 Both systems work simultaneously in an antagonistic nature, with an increased activity in one, inevitably causing a decrease in the other. However, neither substructure of the ANS is ever completely disabled by the response of the other. While both the branches of the ANS serve to regulate the internal organs by stimulating or inhibiting their activity, these two substructures respond to different biochemical transmitter substances [5]. The actions of the PNS are started by the release of acetylcholine. Parasympathetic nerve fibers slow the heart rate, narrow the bronchioles, and cause the release of insulin and digestive fluids. On the other hand, the actions of the SNS are started by the release of norepinepherine. SNS responses include relaxing the tubes of the bronchioles, speeding up the heart, and slowing the digestive tract muscles, while increasing the release of glucose from the liver [5]. Given the broad spectrum of physiological events controlled by the ANS, it is important to establish criteria to evaluate its tone. Past research has indicated, through the measurement of heart rate variability, an inference can be made regarding the functionality of the ANS [3,6,7]. HRV is defined as beat-to-beat fluctuations in heart rate attributed to electrical impulses generated by the sino-atrial (SA) node. The SA node is commonly referred to as the heart's internal pacemaker, and the vagus nerve, or the Xth cranial nerve, inhibits this pacemaker. The PNS influences the SA and AV nodes via the vagus nerve. This influence is produced by release of the neurotransmitter acetylcholine at the vagus nerve endings, which result in the slowing of activity at the SA node and a slowing of the cardiac impulse passing into the ventricles. The SNS has the opposite effect through the release of norepinepherine at the sympathetic nerve endings, thereby increasing the heart rate [5].

PAGE 18

6 HRV reflects the response of the ANS to exogenous or endogenous perturbations. Numerous factors contribute to HRV including blood pressure, temperature, respiration, state of oxygenation, ventilation, biochemical influences of the acid-base balance, and psychological parameters [3]. Continuous changes in sympathetic and parasympathetic neural impulses induce changes in heart rate and cause oscillations around the mean, defined as HRV. By studying heart rate variability, an opportunity exists to study cardiac dynamic behavior and functionality of the ANS. Since the Autonomic Nervous Systems control of heart rate is dynamic and nonlinear, it is necessary to explore viable parameters, which can be used to access the causes of these beat-to-beat fluctuations. This thesis will address time-domain, linear frequency-domain, and nonlinear frequency-domain characteristics of the HRV. Thus the framework for the study has been set, with an attempt to answer each of the following questions: 1. What clinically significant features can be extracted from HRV? 2. What are the developmental changes in HRV? 3. Can we establish criteria to assess vagal tone? 4. Can the premature infant learn or demonstrate a cardiac orienting response, and does HRV change as learning emerges? Developmental Changes Immediately following birth, the heart rate is extremely sensitive to physiological state of the infant. The developmental changes in HRV are hypothesized to indicate a change in the parasympathetic tone. While early fetal life indicates a strong dependence on the SNS, these two systems trend toward an equal balance in the term infant [3]. This thesis assumes a parallel between the fetus and the premature infant. The developmental changes in heart rate variability can be measured through time and frequency domain statistics. Developmental trends of HRV will be measured across subject groups and

PAGE 19

7 gestational age. Time and frequency domain indexes will be used to access maturation of the neonates ANS, or ability to regulate homeostasis. Vagal Tone and Homeostasis Cardiac vagal tone is a construct that describes the functional relationship between the brainstem and the heart. Cardiac responses to brief low intensity stimuli should be deceleratory, indicating receptiveness to incoming stimuli. Cardiac responses to sustained stimuli should be acceleratory, indicating attempts to shut it out [8]. Vagal tone reflects individual differences in degrees of influence of the vagus on the heart rate. High vagal tone is associated with the ability to attend selectively to stimuli and maintain attention during a task, the ability to maintain homeostatic integrity, and infers healthier ANS functioning over the long-term. Homeostasis is primarily maintained by the PNS and provides a physiological framework for the development of complex behaviors. Therefore infants with a high parasympathetic tone are more efficient regulators of homeostasis than those dominated by sympathetic responses [6]. The PNS is primarily responsible with coordinating the life-support process, which intrinsically provides maximum growth and restoration of the bodys tissues and organs [9]. More effective regulation of homeostatis may be an important factor underlying the relationship between high resting vagal tone and favorable neurobehavioral development [6]. Moreover, vagal stimulation controls life-support processes including sucking, swallowing, thermoregulation, and vocalization [9]. The search for an accurate measurement of parasympathetic and sympathetic tone has only begun to express its value. This thesis is anticipated to contribute to the comprehension of these physiological parameters.

PAGE 20

8 Cardiac Orienting Response The existence of a cardiac orienting response (COR) suggests that HRV includes a neurongenic response. The deceleratory parasympathetic component is provided by the bi-directional vagus nerve, which includes both efferent and afferent fibers. Efferent fibers originating in the brain stem terminate on the SA node. Vagal afferent fibers originate throughout the heart and project to the nucleus tractus solitarius [10]. These fibers serve as a continuous response and feedback mechanism between the heart and the brain, facilitating regulation of cardiac function. Porges et al. suggest that the time course of the response, the effects of neural blockades, and studies with clinical populations support this theory. As noted by Porges the heart rate deceleration associated with the cardiac orienting response is rapid and usually returns rapidly to baseline. Also, the latency characteristics of the cardiac orienting response are similar to the opto-vagal, vaso-vagal, baroreceptor-vagal and chemoreceptor-vagal reflexes [8]. Since short latency heart rate reactivity is mediated by the vagus nerve, the magnitude of the cardiac orienting response may be an index of vagal regulation. Respiratory Sinus Arrhythmia Repiratory sinus arrhythmia (RSA) has long been used to estimate vagal tone because of its association with PNS. RSA results from increases in vagal efference during exhalation, which decelerates heart rate, and decreases in vagal efference during inhalation, which accelerate heart rate. Recently, the accuracy of the correlation between RSA and vagal tone has been questioned because it only partially accounts for the beat-to-beat fluctuations of the heart [10]. This thesis will address RSA as a possible criterion to extract vagal tone, but due to the nature of the study will not pharmacologically prove

PAGE 21

9 the efficiency of this technique. However, since RSA accounts for a significant amount of variability in heart rate, it is impossible to completely ignore this aspect.

PAGE 22

CHAPTER 3 PRELIMINARY RESULTS AND STUDY PROCEDURES Preliminary Findings In a preliminary study by the principal investigator [7], cardiac activities during two consecutive quiet periods of 16 healthy fetuses were measured. The exploratory phase used a repeated measures design to follow the normative changes in fetal HRV between 28-34 weeks. Changes in HRV were explained by normalized power shifts in the Fast Fourier Transform (FFT) frequency spectrum. A multi-group pretest-posttest experimental design determined whether 28-34 week-old fetuses would become familiar with a rhyme, recited by the maternal parent. Developmental changes in the HRV frequency spectrum occurred in a nonlinear fashion and individual differences were apparent. The total power diminished significantly with increasing gestational age. Furthermore, the study provided conclusive evidence that the 28-34 week old fetus was capable of responding to a nursery rhyme their mother recited aloud. The response to the nursery rhyme was described by a shift in the power of HRV frequency spectrum to the high frequency bands, which is synonymous with an increase in parasympathetic activity. This conclusion parallels all hypothesized learning capabilities of the premature infant. The testing sessions for the preliminary study utilized a recording of a rhyme spoken by an unfamiliar female, which insured the fetuses reactions would be elicited by familiarity with the acoustic properties of the rhyme. Once learning was detected, ongoing interactions between learning, movement and HRV were detected. A 10

PAGE 23

11 nonparametric signed test was used to compare differences across time. The preliminary study revealed a significant downward trend in fetal movement and an upward trend in a high (.17-.40 Hz) frequency once a cardiac orienting response or learning was identified. Since the subject size was small, there exists a need for a more expanded and complete study before any conclusions can be drawn. However, the findings of the fetal study provide the framework for this expanded research, conducted in collaboration with the PI of the preliminary study. Proposed Study A convenience-based sample of 28 low-risk premature newborns admitted to the NICU at Shands Teaching Hospital are recruited, in blocks of four, and randomly assigned to one of two groups. The target age for admission into the study is 27-28 weeks post-conceptional age. Group 1 will begin exposure to the nursery rhyme at 28 weeks and those in Group 2 will begin exposure at 32 weeks. Staggering the groups shall assist in determining whether an orienting response was a function of learning or simply a reflexive response to presentation of the rhyme. All subjects will be excluded from analysis if they are subject to one or more of the following: An abnormal head ultrasound A sensorineural hearing loss A confirmed prenatally transmitted virus/bacterial infection Cardiac abnormalities Table 3-1. Organization of Recitation and Test Sessions Group 28 weeks 29 weeks 30 weeks 31 weeks 32 weeks 33 weeks 34 weeks 1 Recite Test Recite Test Recite Test Recite Test Recite Test Recite Test Recite Test 2 Test Only Test Only Test Only Test Only Recite Test Recite Test Recite Test

PAGE 24

12 Recitation Sessions The mother of the infant will complete a recording of a nursery rhyme, identical to that presented in the preliminary study. The nursery rhyme, taking approximately 15 seconds to recite, will be played 3 consecutive times to provide the infant with 45s of auditory stimulation. The maternal recitation sessions will occur once in the morning and again in the evening, each day; Group 1 will begin at 28 weeks and Group 2 at 32 weeks of age. The purpose of the recitation session is to introduce the neonate to the rhyme using a familiar voice. The nursery rhyme will be presented once the premature infant is judged to be in a state of active sleep and at least fifteen minutes following a feeding [11]. Between 28-34 weeks post-conceptional age active sleep accounts for appproximately 70% of the entire sleep-wake cycle [12]. Learning in the normal newborn, can occur during periods of active sleep [13,14] and have occurred in the fetus during quiet periods of their quiet-activity cycle [7]. Criteria established by Thoman and Whitney [15] and Holditch-Davis [12] will be used to detect active sleep. The subject is in an active sleep state when: Eyes closed Respiratory rate under 25 breaths per minute Limited movement of extremities Subjects will be monitored at the same time each week because of potential circadian influences on heart rate patterns and movement [16]. A hospital log of recitation times and the time since the neonates last meal will be kept because variations related to unexpected changes in hospital routines are expected. These recitation sessions will be completed entirely separate from the testing session described below.

PAGE 25

13 Test Sessions Each group will undergo a 480 second ECG recording during each test session. Once more, the criteria established by Thoman and Whitney [15] and Holditch-Davis [12] will be used to detect active sleep before the test session begins. The ECG recording will be divided into two stages. The first stage of this study incorporates an exploratory, two-group pretest-posttest subject set which will be used to methodologically replicate the findings from the preliminary study in the premature newborn during a similar developmental time period (28-34 weeks post-conceptional age). The exploratory design will follow specific changes in HRV, measured by time series analysis, power spectral analyses, and time dependant spectral analysis over 28-34 weeks post-conceptional age. A 240-s window separated from the delivery of the auditory stimulus and commencement of the second phase will be used to represent developmental changes in HRV. It is hypothesized that developmental changes in HRV will occur in a decreasing, nonlinear fashion over 28 to 34 weeks gestational age. Two interactions (gestational age vs. gender and gestational age vs. group assignment) will be included in the mixed general linear model. The mixed model allows for the individual random effects when estimating fixed population effects [17]. During a second period of confirmed active sleep, a two-group pretest-posttest experimental design will be used to document when a cardiac orienting response to a presented nursery rhyme first occurs, and how movement and HRV change once it is detected. Again, a nursery rhyme identical to that presented in the preliminary study will be used. An unfamiliar, English-speaking female performs the recitation of the nursery rhyme used in these test sessions. To avoid additional exposure to this unfamiliar voice, this person will not be in attendance for presentations of the rhyme to the infant. The

PAGE 26

14 nursery rhyme takes approximately 15 seconds to recite, and will be repeated 3 consecutive times to provide the 45s window of auditory stimulation. The testing sessions using the unfamiliar recording of the nursery rhyme will occur once a week for both groups from 28-34 weeks post-conceptional age. Figure 3-1. Time-line for the windowing of test data. This second stage will be separated into three 45s windows, in which statistical analysis of heart rate variability will be individually analyzed. The first window will be used to detect the baseline state of the infant before the presence of any auditory stimulation. The second window, immediately following the first, chronologically corresponds to the recitation of the nursery rhyme. This window will be used to observe time and frequency domain changes in HRV, which may reflect the hypothesized COR. A third consecutive window will be used to examine how the neonate responds after the auditory stimulation has been removed from the environment. The third window will be used to verify that an observed COR response is not merely a result of coincidental changes in HRV. In addition, a count of movements will be recorded onto the hospital log during each test session. Our hypotheses pertaining to a COR suggest, (a) over 28-34

PAGE 27

15 weeks of gestational age, it will take 2-5 weeks of recitation to detect and maintain a cardiac orienting response in each subject and (b) once a cardiac orienting response has been observed, subsequent testing session will demonstrate a trend of movement counts and specific HRV characteristics.

PAGE 28

CHAPTER 4 TESTING The rhymes will be played over a 2.5-inch Radio Shack speaker positioned 20 cm from the infants ear. Both versions of the rhyme (maternal and unfamiliar) will be delivered at 55 dB ( + 5 dB) as measured by a Bruel & Kjaer sound level meter (2239) using A-weighted sound pressure levels [18]. This decibel level is just below the normal levels of conversational speech (58-60 dB) and is significantly lower than the background sound levels reported in the literature (70-80 dB + 20 dB) for NICUs [19,20]. The decision to present both versions of the rhyme (maternal and unfamiliar) at 55 dB ( + 5 dB) was based on preliminary findings in the premature fetus where a sound level of 75 dB was used [7]. This decibel level corrects for a 20-decibel attenuation due to passage of sound through the uterine wall [21], bringing the sound level down to approximately 55 dB or to just below the level of normal conversational speech (58-60 dB). The study will record decibel levels at the infants ear prior to presentation of the maternal and unfamiliar rhyme. In the preliminary study, this sound level (75 dB) did not elicit a cardiac deceleration until a history of recitation of the nursery rhyme had occurred [7]. An RS232 interface from the Agilent Neonatal CMS 2001 neonatal monitor to an IBM compatible computer (Dell Inspiron 8100 Laptop) was constructed to obtain the ECG signal and respiratory rates. The data acquisition program is designed to run in DOS and may be activated from the Windows Desktop of Dell Inspiron 8100 Laptop by double clicking on the Neonatal Data Acquisition icon. The C program, originally written by Hewlett Packard (HP), was modified from an earlier version to in order to 16

PAGE 29

17 record the ECG and respiratory rate. The HP program was initially written to demonstrate that information from their monitors could be transferred and displayed on peripheral devices through RS-232 interfacing. The sample program was only available in a DOS format. Although this format is not optimally user friendly, an extensive user interface was not necessary to transfer the required signals into readable text files. Furthermore, HP uses a complex communication protocol, which makes it desirable to avoid writing an original program. The original program was capable of display the ECG signal on a peripheral monitor, which implies that the ECG information was already being transferred. The program modification involved the following: Requesting the respiratory rate signal and displaying the information on the screen. Creating a user interface to request the patient number. Transferring the requested information to a readable text file. Before data acquisition can begin, the Agilent monitor must be set to a baud rate of 19200 and the transition must be set to a high to low trigger. If the baud rate or transition is not correctly set, the program will time-out after eight seconds. When the data acquisition program is started the user enters the patient number and is asked to confirm this number before the recording process begins. The program passes a text string, which signifies the rhyme state and allows the user to mark the data when the rhyme is presented and upon completion. The three measurements: ECG, respiratory rate, and the rhyme state are recorded in an ASCII text file labeled with the month, date, hour (in 24h format), and minute (MMDDHHmm.txt) the program is initialized. ASCII was the chosen format for data storage to permit the user to utilize other analysis programs (i.e. Microsoft EXCEL). The parameters are stored in the format directly readable as a MATLAB matrix and can later be separated into three separate vectors. Preceding the

PAGE 30

18 matrix is an ASCII timestamp, which is recorded once the program has finished the initialization process (DDHHmmss), and appended by another timestamp when the program is terminated. These time stamps allow the user to verify the length of the recording. The program is designed to work with an ECG/Resp module, which will transfer both the ECG and the respiratory rate signal to the Dell Inspiron 8100 Laptop. However, in the absence of the ECG/Resp module it is possible to use only the singular ECG module. This will not record the respiratory signal, but will still record the ECG signal. The first column recorded in the ASCII text file represents the unsigned 12-bit digital values of the ECG signal. The second column contains the information regarding the rhyme state. The final column contains the respiratory rate. The resulting format of the text file is illustrated in the figure below. Figure 4-1. Example of text file generated from Data Acquisition Program

PAGE 31

19 All columns are recorded at a sampling rate of 500Hz, which was preserved from the HP demo program as the rate information was sent to the screen. This sampling rate, much higher than the 200Hz often described in literature, was preserved to prevent jagged QRS complexes. The column corresponding to the state of the rhyme is initialized to , until the r key is pressed and the value changes to . The value remains at until the d key is pressed, which signifies the end of the presentation and returns the value to . When all desired data has been recorded, the user presses the q key twice to exit the program and return to the Windows desktop. Figure 4-2. The setup for the test sessions.

PAGE 32

20 Figure 4-2 illustrates the equipment setup during the test sessions. A male 25-pin RS232 port in the back of the Agilent monitor is interfaced with the female 9-pin RS232 port Dell Inspiron 8100 Laptop. The RS232 cable provides means of serial communication between the Agilent monitor and the Dell Laptop. The Agilent Neonatal CMS 2001 monitor continually monitors the infants vital signs, with a standard 3-lead ECG setup. The speakers, connected to the CD CopyWriter Live, are placed 20cm from the infants ear. This setup is replicated for every test session, regardless of the subjects group number.

PAGE 33

CHAPTER 5 QRS DETECTION Once all data from the testing session is acquired the focus shifts to the data analysis. The first goal is to perform time-series operations, yielding discernable characteristics of the HRV. The HRV is basically controlled by the effect of the vagal and sympathetic influences of the sino-atrial node. Therefore, it seems the best representation of the HRV variability would be expressed as a function of the time intervals between consecutive P waves. However, due to the inconsistent positioning and suboptimal attachment of the infants ECG leads, P waves characteristically have a low signal-to-noise (SNR) ratio [22]. Therefore, in the absence of cardiac arrithymias, addressed in the exclusion criteria, it seems more reliable to replace the detection of PP intervals (PPi), with RRi. While this may include an extra variable in the analysis of the HRV, it seems trivial compared to the accuracy of the detection of PPi. It is our assumption that variations in the PR intervals will be equal or less than error of detection of the PPi. When considering an error-percentage of P-wave detection greater than the standard deviation of the PR interval and greater than the percent error in R-wave detection, the replacement of PPi with RRi becomes trivial. The occurrence of the following cardiac arrhythmias causes this assumption to fail [23]: ATRIAL FLUTTER the atrial flutter waves, known as F waves, are larger than normal P waves and they have a saw-toothed waveform. Not every atrial flutter wave results in a QRS complex because the AV node acts as a filter. Some flutter waves reach the AV node when it is refractory and thus are not propagated to the ventricles. The ventricular rate is usually regular but slower than the atrial rate. 21

PAGE 34

22 ATRIAL FIBRILLATION Atrial fibrillation occurs when the atria depolarize repeatedly and in an irregular uncontrolled manner. As a result, there is no concerted contraction of the atria. No P-waves are observed in the EKG due to the chaotic atrial depolarization. The chaotic atrial depolarization waves penetrate the AV node in an irregular manner, resulting in irregular ventricular contractions. The QRS complexes have normal shape, due to normal ventricular conduction. However, the RR intervals vary from beat to beat. VENTRICULAR TACHYCARDIA Ventricular tachycardia occurs when electrical impulses originating either from the ventricles cause rapid ventricular depolarization. Since the impulse originates from the ventricles, the QRS complexes are wide and bizarre. Ventricular impulses can be sometimes conducted backwards to the atria, in which case, P-waves may be inverted. Otherwise, regular normal P waves may be present but not associated with QRS complexes. The RR intervals are usually regular. THIRD DEGREE AV BLOCK Atrial rate is usually normal and the ventricular rate is usually low, with the atrial rate always faster than the ventricular rate. P waves are normal with constant P-P intervals, but not associated with the QRS complexes. QRS may be normal or widened depending on where the escape pacemaker is located in the conduction system. The result is unrelated atrial and ventricular activities due to the complete blocking of the atrial impulses to the ventricles. Subjects whose electrocardiogram display any of the arrhythmias listed above are removed from analysis on the basis of violating the exclusion criteria. In conclusion, all ECG time-domain analysis of the HRV will consider RRi rather than the ideal analysis of the PPi for the remainder of the thesis. The first step in the analysis of the HRV is to extract the RRi from the recorded ECG waveforms. The ASCII text formatted matrix addressed in the data acquisition chapter is read into MATLAB. The result is an N by 3 matrix, with N being the number of samples. This matrix is then separated into 3 vectors of length N. Vector one contains the ECG information; vector two contains the respiratory rate; and vector three contains information regarded the state of the rhyme. Only vector one expresses useful information for analysis of RRi. Figure 5-1 describes the transform of the raw ECG signal to a vector containing RRi.

PAGE 35

23 Figure 5-1. QRS detection algorithm Noise/Baseline Wander The difficulty in accurate R-wave or QRS complex detection is directly related to the common noise characteristics existing in a raw ECG signal. Noise sources include, but are not limited to: electrode motion, muscle movement, power-line interference, baseline drift, and T-wave interference. Using a method similar to that described by Jiapu Pan [24] we will attempt to improve the SNR of the raw ECG signal and detect the QRS complexes. This process is split into two linear steps and one nonlinear transformation. The first linear process requires filtering the signal with a bandpass filter, which is split into a low-pass and high-pass filter. A low-pass filter with cutoff frequency of approximately 11Hz is first applied the raw ECG signal to remove the low-frequency baseline drift. In addition, the filtered ECG signal is sent through a high-pass filter to remove muscle movement, powerline interference, and T-wave interference. The cutoff frequency of the high-pass filter is approximately 5Hz. This yields a filtered ECG

PAGE 36

24 signal (fECG) containing frequencies of 5-11Hz, which has been described to preserve maximum information regarding the QRS complexes [24]. Low-pass filter: 212611zzzH High-pass filter: 132161321zzzzH Filtering with the Pan-Tompkins method described above, yielded only minor improvements to the raw ECG signal. Problems with the above technique arise as a result of the electrode placement differing from subject to subject. In some case, the Pam-Tompkins algorithm removed substantial frequency components of the QRS complex, thereby deteriorating the quality of the signal. Therefore, we need to replace the low and high pass filters described above with a notch and bandpass filter. The nature of the ECG signal calls for the use of an IIR filter. We will explore the Butterworth and Chebyshev filters. Both are high order filter designs, which can be realized by using simple first order stages cascaded together to achieve the desired order, passband response, and cut-off frequency. The Butterworth filter yields no passband ripple and a roll-off of -20dB per pole. However, the flatness of the passband response comes at the expense of roll-off. The Chebyshev filter displays a much steeper roll-off, but is characterized by a significant passband ripple. The Chebyshev filter thereby optimizes roll-off at the expense of passband ripple. Due to the nature of the ECG signal the roll-off frequency of the passband is not as important as the gain in the passband. Therefore, we will use Butterworth filters for all filtering of the raw ECG signal.

PAGE 37

25 Since 60Hz powerline interference is common among electrical signals, we need to construct a notch filter to remove this noise. The notch filter used is a 5 th order Butterworth filter with center frequency at 60Hz. The magnitude and phase of this filter is shown in Figure 5-2. Figure 5-2. Magnitude and phase response of the 60Hz notch filter used to filter the raw ECG signal. The phase of this filter is essentially linear over the region of attenuation. After the raw ECG signal is passed through the notch filter, it is passed through an 8 th order bandpass filter serving to remove baseline wander and high frequency noise; thereby improving the SNR of the raw ECG signal. Both the notch and bandpass filters are implemented using Butterworth filters. As illustrated in the Figure 5-3 below, the phase over the passband region is approximately linear. Keeping the phase of the filter linear over the passband prevents phase distortion in the resulting filtered ECG waveform.

PAGE 38

26 Figure 5-3. Magnitude and phase response of the bandpass filter used to filter the raw ECG signal. During the recording phase, we obtained ECG recordings that included significant periods where the signal underwent high and low voltage clipping. The clipping resulted in epochs of flat line, thereby making it impossible to detect subsequent QRS complexes. For HRV analysis, a continuous ECG signal is needed. The undesired clipping resulted in extensive signal losses that required considerable epochs of QRS interpolation. It was later discovered that this problem could be overcome by setting the ECG module to filter the signal before it was relayed to the laptop for text conversion. In addition, incorrect electrode placement proved the monitors default selection of Lead II to be an errant assumption. Lead selection was then delegated to the user before the acquisition program was initialized. After realizing the need to select the filter option and correct lead, we were able to attain signals that needed little or no secondary filtering and was not necessary to bandpass filter the recordings. Powerline interference during the passage

PAGE 39

27 from the monitor to the laptop was not completely alleviated and the 60Hz notch filter, described above, was still applied. When the detection phase demonstrated a need for secondary filtering, the bandpass Butterworth Filter, described above, was applied to the signal to conclude the filtering stage of the detection algorithm. Slope Detection In parallel with the Pan-Tompkins algorithm, after filtering, the fECG signal is sent through a differentiator yielding QRS complex slope information; this method results in a differentiated version (dfECG) of the fECG. In accordance with Pan [24], we will use a 5-point transfer function to approximate the derivative. 21122281zzzzTzH where T is the sampling period. This transfer function is roughly linear in the region between dc and 100Hz. Therefore it provides a quick approximation of the ideal derivative over the region of interest. Figure 5-4. Magnitude and phase response of the differentiator function

PAGE 40

28 Following this step, the dfECG signal is passed through a nonlinear point-by-point squaring operator yielding the sdfECG signal: 2nTxnTy This operation makes all data points positive and supplies a nonlinear amplification of the signal at the higher frequencies. Finally, the sdfECG is smoothed using a ten point averaging function, which produces a signal containing the power information of the ECG waveforms slope, pECG. This reduces the presence of local minima in the sdfECG signal. From this point, the method used for detection of QRS complexes deviates from that used by Pan [24]. Since it is not necessary to perform QRS detection in an online fashion, we have the luxury of determining thresholds using the entire information contained in pECG signal. The threshold limit of the QRS slope detection is set in the following fashion: 2/max22dfECGmeandfECG The selection of the current threshold () allows for detection of poorly defined QRS complexes, while distinguishing the R-waves from occurrences of P and T-waves. The threshold is then adjusted throughout the detection process. The following equations illustrates how the threshold is updated: 2/):():(max22winiidfECGmeanwiniidfECG where threshold parameter () is updated in accordance with the current sample (i) over a user-defined window (win)

PAGE 41

29 where represents is the threshold difference and includes a memory constant () to reduce high frequency noise interference during the update Once this threshold has been set, it is now necessary to find maximum values during each interval that the signal spends above threshold. The detection algorithm searches the pECG until it reaches a y-value above the threshold. An array is constructed, storing the waveform from this point until the y-value falls below the threshold. The x-value corresponding to the maximum y-value inside this array is stored as a possible QRS complex. mnnnxQRS,,2,1,max where y[n] and y[m] are above the threshold, but y[n-1] and y[m+1] fall below the threshold. Figure 5-5. Visualization of QRS peak detection

PAGE 42

30 These maximums, or peaks, define the occurrences of possible QRS complexes. If a QRS complex is detected outside the expected range given below, *)()()(RRimeankQRSRRimean 5. RR(0) = 0.375s (160bpm) 6. 1<<2 7. is the detected QRS complex )(kQRS the algorithm reverts back to the previous QRS complex, adjusts the threshold as follows, threshold threshold is a user definable parameter (default = 1.5) and searches for a new QRS complex. This searchback technique plays a crucial role in preventing high frequency noise contamination, P-waves, and T-waves from being marked as a QRS complex. If a complex is not found within the expected range after 5 searchbacks, a complex is interpolated as an average of the previous two intervals. 2)2()1()( kQRSkQRSkQRS Subsequent detections are stored in a vector of increasing length until the algorithm has operated on the entire pECG signal. In addition, the interpolated detections are stored in a separate vector for further analysis. Units of these vectors are stored according to sample numbers. The vectors are then divided by the sampling rate for unit conversion to the time-domain. Finally, another vector is constructed, which represents a time-based representation of the RRi in the following fashion: 1,,23,12 NQRSNQRSQRSQRSQRSQRSRRi

PAGE 43

31 where N = total number of QRS detections Figure 5-6 illustrates the final results of the QRS detection algorithm where the QRS complexes are marked in red. This figure also displays a linear interpolation of the RRi waveform corresponding to the detected beats. Figure 5-6. Example of detected QRS complexes and beat intervals The detection is offset from the peaks of the QRS complex due to the intrinsic phase delays of the filtering process. As a result, the delays are consistent throughout each signal and fall out after computing the RRi vector, as follows: 111 NQRSNQRSNQRSNQRSNQRSNQRSRRi

PAGE 44

32 Outliers Outliers are ill-defined byproducts of data acquisition systems. There is no rigorous definition of an outlier; it is only generally conceived that selectively eliminating inconvenient data points from analysis creates a more robust dataset. Outliers are acknowledged to skew sample means, but tend to have less affect on the median than the mean. For the purpose of this thesis, it is necessary to form a working definition of an outlier. An outlier is a data point that is an "unusual" observation or an extreme value in the data set. The lack of a universally accepted definition of an outlier inherently causes a deficiency of objective mathematical processes for outlier recognition. Adding to the complication of outlier detection, RRi datasets are time dependant vectors, rather than independent samples. Therefore, it is essential to add to our working definition of an outlier. A time dependant outlier is a data point that is an "unusual" observation or an extreme value in the data set, which is not part of a trend in time. Again, the definition of a trend is subjective and needs to be defined for the purposes of this paper. A trend is two or more consecutive data points, which move in the same general direction within a given statistical range. These definitions will be used to construct an algorithm to remove time dependant outliers from the RRi datasets. False detections and missed detections result in possible outliers in these datasets. The presence of outliers in the RRi dataset adversely effect both time and frequency domain analysis of HRV. Possible skewed time domain representations of HRV include the mean and standard deviation of RRi windows. In addition, these outliers produce ghost powers in the high frequency range of the

PAGE 45

33 spectrum. Not only does this increase the power in the high frequency bands, it skews the power in the low frequency bands. Low frequency bands are of the foremost importance during analysis of HRV in the frequency domain. This will become more evident through discussion in Chapter 6. Now that the importance of removing outliers has been established, we must construct a mathematically sound algorithm for their removal. Three techniques are explored as possible strategies for outlier detection: direct RRi analysis, slope analysis, and differential analysis. For each of the aforementioned techniques, we will use standard deviation as the exclusion criterion. Since literature is not conclusive for the distribution of RRi datasets, the exclusion criterion is based on normal distribution of the data. The figures below illustrate the RRi distribution and the normal distribution. Visually, the normal distribution seems an accurate estimation of the RRi dataset. We will discuss the validity of this assumption further in Chapter 6. Figure 5-7. Histogram of RRi distribution from one subject.

PAGE 46

34 Figure 5-8. Standard deviations of the Normal Distribution In a normal distribution, one standard deviation away from the mean in either direction on the horizontal axis accounts for around 68 percent of the values in the dataset. Two standard deviations away from the mean includes for roughly 95 percent of the data. Finally, three standard deviations account for about 99 percent of the data [25]. Therefore, we will assume data outside of three standard deviations from the mean to be statistically irrelevant and tag them as outliers unless they are part of a trend. Direct RRi Analysis The simplest method of removing outliers would be to examine the RRi dataset directly and remove data that fall three standard deviations above or below the mean. The following equation is used for outlier detection: RRistdRRimeannRRioutliers*3 for n RRilength:1

PAGE 47

35 This method inadequately detects outliers and is not optimal to detect a trend, which moves beyond the bounds of the acceptable range. In addition, the direct RRi analysis method does not account for consecutive false or missed detection that exhibit a bouncing effect around the actual QRS complex. Figure 5-9. Illustrates the bouncing effect and upper and lower bound symmetry problems with using the direct RRi analysis method for tagging outliers. The figure above demonstrates problems, which arise from direct RRi analysis method for tagging outliers. The green lines represent the bounds for acceptable datapoints. The bouncing effect is clearly present in this example. According to the exclusion criteria, one of the datapoints would be tagged as an outlier, while the other point would pass through the algorithm without being tagged. This figure also unveils weaknesses with the computation of the mean. When the heart rate deviates to far from

PAGE 48

36 the mean (marked by the black line), situations exist where the upper and lower bound are not symmetrical about the data over a given window. Given all the abovementioned weaknesses, it is necessary to explore a more robust technique to mark the presences of an outlier. Slope Analysis In order to correct for the bouncing effect described above, slope information could be used to isolate these occurrences and tag outliers. First, compute the derivative of the RRi dataset, as follows: tRRiRRi Then, tagging outliers using the same detection method as in the direct analysis: RRistdRRimeannRRioutliers*3 for n RRilength:1 This algorithm adequately accounts for the removal of the bouncing effect, but we will demonstrate that it is impossible to exclude trends from being tagged. Figure 5-10 is an example of trend which seems to begin abruptly, but is sustained over subsequent intervals. The slope algorithm solves the problems in accordance with the mean and bounce effect, but introduces troubles in excluding trends from being tagged as outliers. Over the course of this trend, only those datapoints whose slope exceeded the bound of acceptance (marked by the green line) would be tagged. No mathematically simple method exists to overcome the flaws in this algorithm. Merely modifying this algorithm to remove the outliers and loop through the process again would serve to completely eliminate the trend. Therefore, for a second time it is necessary to search for a more robust algorithm to tag outliers.

PAGE 49

37 Figure 5-10. Illustrates an RRi trend that would not be preserved using the slope analysis method for tagging outliers. Differential Analysis To this point, three problems in the outlier detection process have been identified: mean calculations that deteriorate the symmetry of the acceptance bound, the bouncing effect, and the exclusion of trends. We introduce the differential analysis method as a means, which circumvents all three of these problems. First, a difference vector is constructed as follows: 1)( nxnxnRRidiff for n, where m:2 RRilengthm

PAGE 50

38 1)1(RRiRRimeanRRidiff The end product is a vector, which is the same length as RRi. The resulting mean of the difference vector is close to zero, which further implies the data is predominately normally distributed. The difference vector also alleviates the first problem mentioned above. Then, a vector is created using the same exclusion criterion discussed in the previous two subsections. However, using this difference vector alone to tag outliers outside the acceptable standard deviation range is not sufficient. A single outlier in a dataset would produce two datapoints in this difference vector that would be tagged as outliers (n and n+1). Being that it is not necessary to perform the outlier removal real-time, we are able to create an additional vector, as follows: 1)(nxnxnRRidiff for n, where 1:1m RRilengthm mRRiRRimeanmRRidiff)( The same exclusion criterion is used to create a vector containing the outliers. As a result of this computation, a single outlier in the dataset would produce two datapoints tagged as outliers. Since taking the difference of RR intervals backwards through time generates this vector, the two resulting datapoints tagged as outliers would be n and n-1. Now, we compare the two difference vectors and note only those values for n, which occur in both. These are the actual RRi datapoints that are separated from the previous and subsequent datapoint by three standard deviations above or below the differential mean. Figure 5-11 is a section of RRi data with trends that would normal be removed using the Direct RRi Analysis method, described above, and a bounce effect.

PAGE 51

39 Figure 5-11. Outliers that have been tagged, using the differential method, are marked by a green dot. Shows the efficiency of the method to preserve trend and remove bouncing effects. In figure above, the points that were tagged as outliers are denoted by green dots. It is now evident that we have constructed a robust method to find and remove outliers, as defined above. Interpolation The RRi vector described above is sufficient for time domain analysis of the HRV. The RRi vector can also be converted to an instantaneous representation of the heart rate (iBPM). However, there is a problem with using this vector to represent the HRV in the frequency domain. When plotting the RRi on a time scale, the x-axis must include a cumulative sum of all the preceding RRi plotted against the RRi or iBPM. This results in

PAGE 52

40 an unequally sampled representation of the HRV. In order to apply a meaningful power spectral representation of the HRV, the data must be evenly sampled at a known frequency. Many papers using spectral analysis fail to address this topic and there are no universal methods to solve this problem. Resolving this issue of uneven sampling is essential for spectral analysis. Without addressing this topic, may forms of spectral analysis applied to the RRi would produce meaningless results. The choice of an interpolation method is very important, as it can significantly compromise the integrity of the data. Therefore, we will explore four methods of interpolating between successive data points of RRi: Linear Interpolation (LI) LaGrange Polynomial Interpolation (LPI) Cubic Spline (CS) Cubic Interpolation (CI) Figure 5-12. LaGrange Polynomial Interpolation

PAGE 53

41 Figure 5-13. Linear Interpolation Figure 5-14. Cubic Spline Interpolation

PAGE 54

42 Figure 5-15. Cubic Interpolation The performance of the interpolation techniques is shown applied to a sample RRi dataset marked by Xs in the figures above. LI is continuous, but the slope changes at the vertex points. LPI has an adjustable order of polynomial and reduces to LI when the order equals one. As seen above, the LPI algorithm does not perform well at the endpoints. CSI produces the smoothest results of all the interpolation methods, but performs poorly if the data is highly non-uniform. CI results in both the interpolated data and its derivative being continuous. The CI technique seems the most realistic interpolation method because it does not force a sinusoidal interpolation. Therefore, we will use CI to interpolate and resample RRi in preparation for spectral analysis. It is necessary to resample the RRi at a frequency where the Nyquist rate is above the frequency range of interest for spectral analysis, discussed further in Chapter 6. Also in

PAGE 55

43 Chapter 6, spectral analysis techniques will be introduced that do not require resampling of the RRi waveform prior examining HRV in the frequency domain.

PAGE 56

CHAPTER 6 STATISTICAL ANALYSIS Time Based The simplest methods used to evaluate HRV are based in the time domain. Time domain statistics can be divided into those derived from RRi or iHR and those derived from differences in RRi, namely RRi [26]. The following time-domain techniques will be used: Mean RRi (mRRi) Median RRi (medRRi) Standard deviation of RR interval (SDRR) Poincar plots Quadrant analysis First Order Statistics Analysis of mRRi and medRRi yields some useful information about HRV. We will use this information to analyze trends of the mRRi and medRRi with gestational age and across genders. Also, the mRRi for the 45-s prestimulus period will be subtracted from each heartbeat of the 45-s prestimulus and stimulus period. The result will be a string of difference scores. The difference scores will be used to confirm that no significant HRV trending has occurred during the prestimulus period, which may mask or falsely demonstrate a COR. During the stimulus period, the difference scores will be used to detect an orienting response. The mRRi is calculated in the following manner: NnRRimRRiNn1)( 44

PAGE 57

45 for N= number of RRi To find the median, we must first sort the RR intervals by ascending values. Then, apply the following formula to find medRRi: )2/(NRRimedRRisorted for even values of N 2)2/1()2/1( NRRiNRRimedRRisortedsorted for odd values of N SDRR is highly dependant upon the length of the recording [26]. When using this value as a comparison across datasets, we must be careful to uniformly define the length of the recording. Therefore, we have kept all window lengths constant from subject-to-subject and session-to-session. The variance (SDRR squared) is mathematically equal to the total power of the spectrum; therefore SDRR will also be used to analyze trends corresponding to increasing gestational age. SDRR is the primary first order statistical analysis method we will use to discuss HRV in the time domain. SDRR is calculated in the following manner: 1)(12NmRRinRRiSDRRNn where N= number of RRi Poincar Plots In addition to numerical analysis in the time domain, the importance of geometric representations, such as Poincar plots are useful tools for HRV analysis [26,27]. The nonlinear Poincar plot produces a figure with each RRi plotted against the previous RRi.

PAGE 58

46 According to Brennan et al., the plot provides summary information as well as detailed beat-to-beat information on the behavior of the heart. The dispersion of points perpendicular to the line of identity reflects the level of short-term variability [27]. In subjects with normal cardiac function, the geometry of Poincar plots tends to be elliptical in nature. No literature pertaining to analysis of Poincar plots using an ellipse were found, but we will explore an ellipse fitting method as a possible representation of short-term and long-variability. MATLAB code written by Marcos Duarte [28] will be used to fit an ellipse to the data. Duarte calculates the major and minor axis of the ellipse using principal component analysis (PCA). The method results in 85.35% of the datapoints lying inside of the ellipse. The datapoints that fall outside the ellipse can be described by increasing the order of the ellipsoid, but it is only necessary to include the first two eigenvalues when analyzing the short and long-term variability of the dataset. PCA is a multivariate procedure, which rotates the data such that maximum variabilities are projected onto the axes. Essentially through computing the eigenvalues of the correlation matrix, a set of correlated variables are transformed into a set of uncorrelated variables, which are ordered by reducing variability. The uncorrelated variables are linear combinations of the original variables, and the last of these variables can be removed with minimum loss of real data. The largest eigenvalue corresponds to the principal component in the direction of greatest variance; the next largest eigenvalue corresponds to the principal component in the perpendicular direction of next greatest variance. PCA can be viewed as a rotation of the existing axes to new positions in the space defined by the original variables. In this new rotation, there will be no correlation

PAGE 59

47 between the new variables defined by the rotation. The first new variable contains the maximum amount of variation. The second new variable contains the maximum amount of variation unexplained by, and orthogonal to, the first. The result is an ellipsoid in multidimensional space and reduces to an ellipse for the two dimensional Poincar Plot. When a Gaussian random vector has covariance matrix that is not diagonal, the axes of the resulting ellipse are perpendicular to each other, but are not parallel to the coordinate axes. As shown in the Figure 6-1, the principal components of the RR interval Poincar representation result in an ellipse whose axes are rotated versions of the coordinate axes. Therefore, as discussed in Chapter 5, the assumption that a Normal distribution can describe RRi datasets is valid. The minor axis of the ellipse serves as an interpretation of short-term variability; and the major axis serves as an index on long-term variability. Figure 6-1. Poincar plot representation of RRi. Since data points overlap, colors, coordinated with the visible spectrum, are used to indicate magnitude at each point.

PAGE 60

48 Quadrant Analysis Quadrant analysis is a graphical tool used to represent the number of sustained increases, sustained decreases, and alternating changes in RRi. Two vectors are created as follows: 1 nRRinRRinRRi nRRinRRinRR 11 Each difference value was then plotted against the previous difference, resulting in the RRi (n) versus RRi (n+1) plot shown below [4]. The plot displays values distributed over an area defined by four quadrants. Each quadrant indicates the direction of two consecutive changes in interval length. Figure 6-2. Defines the quadrants used in Quadrant Analysis.

PAGE 61

49 We will use this analysis tool to examine sustained increases or decreases of RRi and alternating RRi. The number of points in quadrants A and D measure high frequency variation, or quick changes in RRi. Low frequency heart rate variation is characterized by many sequential increases and decreases in interval length, which is indicated by points in quadrants B and C. If changes in RRi occurred randomly and independently, all quadrants would have an approximately equal number of data points. In addition, if the heart rate tended to increase quickly and decrease slowly, more points would fall in quadrants A and B than in C and D and visa versa. Figure 6-3 is an example of a quadrant analysis plot extract from real RRi data. Figure 6-3. Example of Quadrant Analysis applied to an RRi dataset.

PAGE 62

50 Frequency Based While analytical methods were first restricted to the time domain, they have quickly shifted to the frequency domain. Spectral analysis is now a common representation of HRV. Transforming HRV to the frequency domain provides a means of separating the parasympathetic and sympathetic responses of the ANS. Heart rate fluctuations in the premature will be examined in the .00033-1Hz, while the frequencies below this will be discarded to reduce the influence of slow trends or artifacts. Research has suggested that the power spectral decomposition of the ECG yields frequency bands essential to analysis of the ANS. These frequency bands are very inconsistent throughout the literature referenced for this subject. This thesis will follow the frequency divisions described by Stein et al. [29]: The ultra low frequency band (ULFB) <.0033Hz The very low frequency band (VLFB) 0.0033 0.04Hz The low frequency band (LFB) 0.04 0.15Hz The high frequency band (HFB) 0.15 0.4Hz The total power (TP) .0033 1Hz The low frequency band of the power spectrum is composed of two peaks. The mid-range peak (.09 .15Hz) corresponds to the Mayer waves, which have been related to the baroreceptor reflex [3]. In addition, there exist a second peak between .02 .09Hz, which has been linked to thermoregulatory fluctuations and variations in vasomotor tone. The high-frequency band reflects the parasympathetic response, while the low frequency band is influenced by the PNS and SNS [3]. For each frequency band listed above, we will extract the normalized power over the described range. Also, the power ratio of the

PAGE 63

51 high and low frequency bands has been assumed to reflect the sympathovagal balance [3]. This power ratio is observed as follows: HFPLFPLFPSVratio The SV ratio should provide us with insight into the ANS balance and should trend downward in accordance with increasing gestational age. Lastly, a range just above the high-frequency band contains information regarding the RSA, which is not represented by a well-defined peak, but is usually centered about the .5-.6Hz range. We will explore three methods of power spectral analysis in the following subsections: Fast Fourier Spectral Analysis Lomb-Scargle Spectral Analysis Empirical Mode Decomposition Fast Fourier Spectral Analysis The simplest and most commonly exploited technique for spectral analysis is the Fast Fourier Transform (FFT). The FFT is commonly used because it is computationally simple algorithm for spectral analysis. Fourier analysis decomposes a time-series into global sinusoidal components with fixed amplitudes. As a result, Fourier spectral analysis is limited to systems that are linear and strictly stationary. However, a high number of variables affect the heart rate, which suggests that an RRi signal is highly nonlinear and non-stationary. Most research pertaining to HRV utilizes the FFT to describe to the spectral components of the RRi signal. We will utilize the FFT spectral analysis method to obtain results, which can be compared with past and future research. The FFT is computed using the internal FFT function provided by MATLAB. Since the FFT method is commonly used in signal processing, we will not discuss this

PAGE 64

52 method in detail, rather we will explain how we computed the transform. All RRi windows for the RRi dataset (developmental, prestimulus, stimulus, and poststimulus) are analyzed separately based on a N-point FFT, where N is equal to the total number of point in the intervals. Often FFT windows are padded with zeros or with repeated copies of the window. Neither of these techniques adds any information to the resulting FFT. Padding only serves to increase the resolution of the transformation. Setting N to the length to the window abolishes the need for padding. Each of the aforementioned windows is extracted from the original recording based on cumulative time, and the resulting windows are unevenly sampled signals. The cubic spline interpolation technique, described in Chapter 5, is applied before computing the FFT and resampling the interpolated signal at 4Hz. We chose 4Hz to prevent aliasing above of the Nyquist frequency (2Hz) from falling into the frequency bands of interest. Finally, the resampled signal is passed through the FFT algorithm and the DC component of the FFT transformation vector is removed by dropping the first term. Results from this algorithm can be found in the following chapter. Lomb Spectral Analysis In order to avoid the consequence of interpolation, we explore the Lomb method to estimate the power spectral density (PSD) of unevenly sampled signals. Laguna et al. concluded that for PSD of these unevenly sampled signals the Lomb method is more suitable than FFT estimates requiring linear or cubic interpolation [30]. The Lomb estimate avoids the low pass effect of resampling and avoids aliasing up to the mean Nyquist frequency: mRRNf*21

PAGE 65

53 The Lomb method is based on the minimization of the squared differences between the projection of the signal onto the basis function and the signal under study [30]. ninntbictx)( Minimizing the variance of n (mean squared error), is synonomous to the minimization of Nnnintbictx12 which results in [30] Nnnintbtxkic1*1 Nnnitbk12 In addition, the signal power at index i of the transformation as described by Lomb [31]: 221*)(*icicktbtxiciPNnninx where kicic Setting the basis function to that of the Fourier transform nitfjnietb2 Then, NekNntfjni122 And fcfPiPXx2 Lomb presented a simplification that permits one to solve for each coefficient independently, known as the Lomb normalized periodogram. First, we solve for the mean and variance [31]:

PAGE 66

54 22SDRRmRR x Then, we generate the Lomb normalized periodogram using the angular frequency thi NniniNnininNniniNnininiXttxtxttxtxP122112212sinsincoscos21 where i is defined as, NnniNnniiitt112cos2sinarctan21 MATLAB code written by David Glover is used to compute this Lomb normalized periodogram [32]. We apply this method of spectral analysis to the identical windows used to compute the Fourier PSD; and we will compare the results in the following chapter. Empirical Mode Decomposition Approximations are the key to understanding complex natural phenomena and often there is a significant advantage to these methods. However, science is continuing to develop ways to describing complex system and behaviors. Why approximate the RRi signal as a linear, stationary signal when tools exist to analyze nonstationary, nonlinear systems? The answer is rhetorical because the only reason describes the simplicity of the computation of the FFT. Given all the assumptions needed to compute the FFT, it is necessary to search for a more reliable spectral transformation. We have already observed the Lomb method, which still assumes the signal is linear and stationary. Now,

PAGE 67

55 using a technique developed by Haung et al. [33], this thesis will compare the effectiveness of the Empirical Mode Decomposition (EMD) as an analysis tool of HRV. The goal of EMD is to decompose a time series into superposition of Intrinsic Mode Functions (IMF). EMD uses a shifting process to produce IMFs, which enables complicated data to be reduced into a form with defined instantaneous frequencies. The IMFs has the appearance of a generalized Fourier analysis, but with variable amplitudes and frequencies. Haung et al. describes EMD as the first local and adaptive method in frequency time analysis. Furthermore, Haung suggests that the advantage of EMD and Hilbert spectral analysis is that it gives the best-fit local sine or cosine form to the local data. The frequency resolution for any point is uniformly defined by the local phase method and is especially effective in extracting low frequencies oscillations [33]. Since the HRV is characterized by low frequency oscillations, this method should produce informative frequency domain representation of HRV. Computation of the EMD algorithm is done in MATLAB using modified code originally written by G. Rilling [34]. The EMD algorithm [33] coded by Rilling begins with the sifting process. First, all local extrema of the original signal are identified. The local maxima are then connected via a cubic spline and the procedure is repeated for the local minima, respectively forming the upper and lower envelopes. Then, we find the first component of the sifting process [33], h. 1 1010)(mhhtRRitXh where is the mean of the envelope 1m The shift process is repeated k times, as follows, until is an IMF (). kh1 1c

PAGE 68

56 kkkkhcmhh111111 Then, a residual is found 11ctXr Finally the IMF procedure is repeated n times until the residual, is a trend or constant. nr 11 nnnchr It is then necessary to determine if all the calculated IMF are actually orthogonal to one another, as the theory suggests. Huang suggests using the follow equation as a check for orthogonality [33]: 022gfgffgccccIO for all IMFs In order to visualize the result of EMD, we must find the Hilbert transforms of the IMF components and represent the RRi signal in terms of time, frequency, and amplitude. We compute the Hilbert transforms using the Cauchy principal value integral to find the magnitude, and phase, tan tn and instantaneous frequencies, n [33]. tttXtYttYtXtaetatiYtXtZtdtttctYjjjjjjjtijjjjjjj arctan122 Finally, we can express the RRi dataset as follows

PAGE 69

57 ttitatRRijnjjexp1 Cleary, the equation above represents the generalized Fourier expansion, but the amplitude and instantaneous frequency are variable with respect to time [33]. At this point, we will use the EMD method to provide us with a visual representation of the power spectrum changes over a prestimulus stimulus poststimulus window. In accordance with the description of a COR, detailed in Chapter 2, we should see significant shift in the magnitude of frequency spectrum over this window. Statistical Analysis of Parameters Statistical analysis of the parameters extracted above is performed using a technique for longitudinal data analysis, called mixed general liner models (MixMod) or random regression models [35]. The MixMod model allows the concurrent identification of both group and individual patterns in longitudinal sets with covariates that vary over time, inconsistently timed data, irregularly timed data, and randomly missed values. Given the vulnerable state of the infants in the NICU, we must perform test sessions in accordance with the health status of the infant. These stability considerations can result in the testing schedules varying from subject to subject and session to session. Fortunately, problems arising from this effect are alleviated using the MixMod analysis. In addition, subjects are often transferred or released before reaching 34 weeks GA. The MixMod model allows us to include the data from these subjects, thereby reducing the amount of data, which would otherwise be lost. Strengths of the MixMod model include the acceptance of multiple characteristics and the exploration of patterns across subjects and groups. Using multiple characteristics, MixMod leaves the statistical correlation and inclusive of parameters up to the investigator [35]. All MixMod analysis is performed by

PAGE 70

58 Dr. Kruger using readily available SAS Proc Mixed software. The results of the MixMod analysis can be found in the Clinical Results section of Chapter 7.

PAGE 71

CHAPTER 7 DISCUSSION AND RESULTS Engineering Results We will explore the results of the signal-processing phase of this thesis in accordance with ECG recordings from one subject. The example subject is not intended to prove the clinical significance of the hypothesis described in Chapter 3. For clinical evaluation of aforementioned hypotheses, please refer to the Clinical Results subsection of this chapter. QRS Detection Algorithm The following table is a representation of the performance of the QRS detection algorithm, detailed in Chapter 5. The percent of QRS complexes is computed, as follows, missedectionsectionsectedQRS###detdetdet where the interpolated detections are tagged as missed detections and are stored in a separate array to determine the QRS detection ratio defined above. The interpolations are a result of no detected QRS complexes falling within the user defined expectancy range, described in Chapter 5. Table 7-1 clearly illustrates that on average the algorithm only misses 0.85% of the QRS complexes in the signal. This detection algorithm is a robust system, which can effectively detected QRS complexes in the presence of noise and muscle artifact and produce reliable RRi datasets. Beyond the detection process, we improve quality of the 59

PAGE 72

60 RRi signal through the implementation of the outliers removal algorithm designed in Chapter 5. Table 7-1. Shows the performance of the QRS detection algorithm HRV Plots and Statistics Scientific findings are based on equally important qualitative and quantitative measures. Quantitative methods observe distinguishing characteristics and elemental properties resulting in numerical comparative analyses and statistical analyses. For each recording, an excel file is created, which contains all the pertinent quantitative measures for each window:

PAGE 73

61 mRR medRR SDRR Number of sustained increase and decreases in HR derived through Quadrant Analysis The length, width, and (W/L) ratio of the PCA fit elliptical representation of the Poincar plot ULF, VLF, LF, and HF power of the Lomb PSD TP and of Lomb PSD ratioSV ULF, VLF, LF, and HF power of the Fourier PSD TP and of Fourier PSD ratioSV These quantitative measures are used for MixMod analysis in the Clinical Results section. Qualitative measures are those associated with interpretative approaches, from the investigators point of view, rather than discrete observations. While qualitative methodologies can lead to biased results, there is no doubt that they should be included in medical research. The power of visualization and quick interpretation of data possessed by humans cannot be replaced by quantitative measures. Therefore, for each recording, six plots are created for qualitative analysis: HR with outlier removal Quadrant analysis Poincar Fourier PSD Lomb PSD Hilbert-Haung spectrum of EMD

PAGE 74

62 HR Plot Figure 7-1. HR plots of example RRi data. From top to bottom they represent: the (a) developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus periods. The HR plot shown above illustrates the original RRi (blue) waveform displayed as instantaneous heart rates and the resulting RRi (red) after outliers have been removed. The figure above is divided into four plots representing the four windows of analysis. Plot (b) illustrates the removal of an outlier, which causes a high frequency spike in the signal. Plot (a) demonstrated the algorithm efficiency of preserving trends in RR intervals. Figure 7-1 is used to visualize cardiac accelerations and decelerations, quantitatively described through Quadrant Analysis.

PAGE 75

63 Quadrant Plot The Quadrant plot shown in Figure 7-2, below, is organized in the same manner as Figure 7-1 with each plot representing a particular analysis window. Quadrant Analysis is a more useful tool for quantitative analysis, but the visual representation can provide insight into the amount of high and low frequency variability. The excel file described above contains some useful information regarding this plot that will later be used for MixMod analysis. Figure 7-2. Quadrant Analysis plots of example RRi data. From top to bottom they represent: the (a) developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus periods.

PAGE 76

64 Poincare Plot The Poincare plot is a valuable tool for both quantitative and qualitative analysis. The Poincar method of RRi interpretation is observed to prove normal cardiac function and analytically evaluated by the fitted ellipse to described short and long-term HRV. Figure 7-3. Poincar plots of example RRi data. From top to bottom they represent: the (a) developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus periods. Again the magnitude corresponding of a given point is represented by the basic colors of the visible light spectrum with red being the next to least frequent (black is the least) and violet being the most frequent. Plot (a) represents the dataset with the most long-term variability or low-frequency changes inferred by the length of the major axis of the ellipse. This is expected because plot (a) is a 240s window, while plots (a,b,c) are 45s

PAGE 77

65 windows. We can also see that during the stimulus period the major and minor axis of the ellipse become significantly longer, perhaps signifying an increase the LF power, HF power, and total power of the HRV spectrum. We will next look at the Fourier and Lomb PSD to confirm this hypothesis. Fourier and Lomb Plots The Fourier and Lomb PSD plots essentially describe the same information; therefore we will observe them together. Figure 7-4. Fourier PSD plots of example RRi data. From top to bottom they represent: the (a) developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus periods.

PAGE 78

66 Figure 7-5. Lomb PSD plots of example RRi data. From top to bottom they represent: the (a) developmental, (b) prestimulus, (c) stimulus, and (d) poststimulus periods. First, we will address the observational hypothesis from the Poincar Plot section. Inspecting the Fourier and Lomb stimulus plots, we note that the area under the curve increases from the prestimulus period to the stimulus period, which means there is an increase in total power. Also, power increases in the LF and HF bands correspond to increases in the length of the major and minor axes of the ellipse, respectively. In conclusion, through qualitative analysis, we have shown that the Poincar plot is a time domain representation of the PSD. Now we will compare the PSD representations of the Fourier and Lomb Spectrum. The Fourier and Lomb plots are comparable over the developmental window,

PAGE 79

67 but the Fourier plot is choppier than the Lomb plot over the windows corresponding to the COR. The choppy Fourier plot is a result of the low frequency resolution of the N-point Fourier transformation. Therefore, the frequency resolution is lower for the short windows analyzed in plots b, c, and d. However, the Lomb computation allows the user to set the frequency resolution. For the purpose of obtaining a smooth transformation, we have set the Lomb PSD resolution to 0.001Hz. To increase the frequency resolution of the FFT, we must pad the data before performing the Fourier transform. For reasons discussed in Chapter 6, we have refrained from doing so. The apparent equivalent Lomb and Fourier PSD plots for the developmental window suggest that the Lomb method is at least as successful as the Fourier for PSD estimation. For clinical interpretation, we will use the Lomb Spectrum when computing the cumulative power in the frequency bands of interest. EMD Plots At this point, we use the Hilbert-Haung Spectrum (HHS) only as a qualitative analysis tool. If the clinical findings indicate the usefulness of this time-frequency domain representation of HRV, we will expand to use EMD as a quantitative analysis tool. Figure 7-6 is the HHS of the example RRi dataset used in this section. Black lines separate the prestimulus, stimulus, and poststimulus periods. We can relate the HHS power shift of the stimulus period to those described by the Poincar and PSD assessments of HRV. All four representations are in agreement that during the stimulus period the LF power, HF power, and TP increase. However, the HHS shows us that this increase occurs just after the rhyme begins only to tail-off significantly shortly thereafter. Initially the EMD seems to be a robust method to describe the HRV characteristics of the COR with respect to time. However, since this method has not often been used for

PAGE 80

68 clinical interpretation, the inclusion of the method in the analysis of the hypothesized COR will be left up to the PI. Figure 7-6. Hilbert-Haung spectrum of the COR windows of an example RRi dataset. Clinical Results At the time of submission, we only have four complete datasets (recordings for the 28-34 weeks GA) and one incomplete dataset (recordings for 28-33 weeks GA) that can be used for clinical interpretation. The current sample size (n=5) is well below the minimum sample size (n=28) described in the design of this study [7]. Therefore, we will not be able to perform MixMod analysis of the data at this time. Instead, we will use a repeated measures ANOVA analysis of the PSD. Repeated measure ANOVA allows us to compare the variance caused by the independent variables to a more accurate error term, which includes the variance effects of individual subjects. Thus, this method

PAGE 81

69 requires fewer participants to yield an adequate interpretation of the data than the MixMod method. Early analysis based on the Poincar plots, Quadrant Analysis, and EMD will not be performed due to the lack of comparable literature. We will exploit the aforementioned methods upon completion of the study, so the PSD results can be used to confirm the clinical significance of these models. The repeated measure ANOVA model returned only two independent variables (HF power and SV power), which passed the equal variance test. We will first examine results of the repeated measures ANOVA analysis of the SV power defined in Chapter 6. Table 7-2. Results of the repeated measures ANOVA with independent variable (SV power). Session n 1 5 0.807 0.0857 2 5 0.729 0.089 3 5 0.846 0.121 4 5 0.878 0.0501 5 5 0.742 0.239 6 5 0.866 0.0836 7 4 0.613 0.186 SV Power00.20.40.60.811.21234567SessionMagnitude Figure 7-7. Plot of mean SV power with standard deviation bars across sessions

PAGE 82

70 Now we examine repeated measures ANOVA analysis of the HF power, also defined in Chapter 6. Table 7-3. Results of the repeated measures ANOVA with independent variable (HF power). Session n 1 5 68.2 25.925 2 5 133.6 144.921 3 5 48.46 26.837 4 5 41.32 23.877 5 5 178.1 172.465 6 5 47.22 58.952 7 4 150.75 126.818 HF Power-500501001502002503003504001234567SessionMagintude Figure 7-8. Plot of mean HF power with standard deviation bars across sessions While Figure 7-7 shows a possible non-linear trend with respect to GA, Figure 7-8 does not elicit any significant trend in time. First, the small sample size doesnt allow for the development of a trend. Secondly, there exist many exogenous variables, which conceal a developing trend corresponding to GA. These variables include, but are not limited to:

PAGE 83

71 Drug interventions and interactions Respiratory support interventions Physiological discrepancies in ANS maturation Changes in NICU environment Incorrect evaluations in GA Patient logs will be analyzed following the completion of the study for drug interventions. The PI will attempt to parallel the interventions with the recordings that contribute most to the standard deviation of each variable. It is known that certain drugs inhibit or even block parasympathetic and sympathetic responses. Furthermore, many of our subjects are on multiple drugs whose interactions with respect to the ANS are not known. Many subjects intermittently require respiratory interventions such as continuous positive airway pressure (CPAP). A closed-loop study showed that adequate application of long-term CPAP therapy produced substantial alterations in the major physiological mechanisms that influence HRV and blood pressure variability [36]. Based upon preliminary results, see Figures 7-7 and 7-8, final results need to include an analysis of how drug and respiratory interventions hinder the ability to identify trends in the ANS landscape. The significance of all exogenous variables mentioned above will emerge as the number of subject in the dataset increases. Unfortunately, these results are not included in this thesis, but we refer you to further publications by the PI and the author of this paper. For a complete time and frequency domain analysis of a randomly selected subject, please refer to Appendix A and B, respectively.

PAGE 84

CHAPTER 8 CONCLUSION We have successfully modified code written by HP to communicate with the Agilent/HP NICU monitor. This program records the ECG, respiratory rate, and rhyme state in a text file at a 500 Hz sampling rate. The text file can be read into MATLAB for QRS detection and HRV analysis. In MATLAB, we have created an algorithm to perform accurate detection of QRS complexes. The average percentage of missed detections is approximately 0.8%, which we feel represents a near optimal detection algorithm. For those QRS complexes that are undetectable, we have constructed an interpolation method to replace the missing data. Using our novel method of tagging outliers, we are able to remove outliers, but preserve trends across time. We have found that the most useful measurement of HRV in the time-domain is through the analysis of Poincar plots. Fitting an ellipse to the Poincar representation of the RRi is usefully in relating time-domain analysis to the frequency domain. Currently, most spectral representations described in literature use the Fourier transform. We have found that the Lomb normalized periodogram is a better representation of the power spectrum of RRi waveforms. The Lomb algorithm avoids the negative effects of interpolation and resampling required for FFT computation. Also, Empirical Mode Decomposition is a favorable method for measuring changes in the frequency spectrum across time. The EMD spectrum preserves more information regarding shifts in the frequency spectrum during the stimulus period than either the Lomb or Fourier spectrum. EMD has been proven to be a useful analysis tool to depict a COR. 72

PAGE 85

73 We have discussed the importance of both qualitative and quantitative measures of developmental changes in HRV and in identifying a COR; we have extracted numerous features using time and frequency domain HRV analyses, which will be used to analyze and test the hypotheses of the clinical study still in progress. We feel we have completed the design phase of the study and need only to increase our sample size. As the sample size increases, the clinical significance of each feature will be evaluated and we will move to interpret our findings and test our hypotheses.

PAGE 86

CHAPTER 9 FUTURE WORK A need exists in the NICU and beyond to determine the state of development of the ANS. Before an infant is discharged from the NICU, it is important to determine their ability to perform certain life support tasks. In Chapter 2, these life-supporting tasks were linked to the maturity of the ANS, specifically the PNS. The list below is compiled from the 1998 suggested guidelines for neonatal discharge established by the American Academy of Pediatrics (AAP) [37]. In the judgment of the responsible physician there has been: A sustained pattern of weight gain of sufficient duration. Adequate maintenance of normal body temperature with the infant fully clothed in an open bed with normal ambient temperature (24C to 25C). Competent suckle feeding, breast or bottle, without cardiorespiratory compromise. Physiologically mature and stable cardiorespiratory function of sufficient duration. Appropriate immunizations have been administered. Appropriate metabolic screening has been performed. Hematologic status has been assessed and appropriate therapy instituted as indicated. Nutritional risks have been assessed and therapy and dietary modification instituted as indicated. Sensorineural assessments, hearing and funduscopy, have been completed as indicated. Review of hospital course has been completed, unresolved medical problems identified, and plans for treatment instituted as indicated. 74

PAGE 87

75 The AAP has established these guidelines, which include subjective observations of cardiorespiratory function and sensorineural assessments. It is our hypothesis that cardiac function, as well as, neural assessments can be objectively realized using an index of ANS maturation. We will create this index of ANS functionality based on HRV comparisons with mature young adults. Below, we suggest a pilot study, which shall provide evidence that an index of maturation can be realized. Record Baseline Adult HRV Data In order to establish a baseline for a mature ANS with a synonymously high vagal tone, it is important to analyze the HRV of the young, healthy adult. ECG recordings will be obtained for a sample of 28 male and female adults meeting the following criteria: Between the ages of 18-30 No history of congestive heart failure, myocardial infarction, coronary artery disease, head trauma, or diabetes Normal sinus ECG rhythm The length of the recordings should be equal to those used in the study of the fetus. This is important we comparing both time and frequency domain features of HRV and the COR. Each subject should be oriented in the supine position to ensure consistency with ECG recordings of the neonates. Before recording, all subjects should remain in this position without auditory and visual simulation, until they are verified to be in REM or paradoxical sleep. Consistency with the neonatal recording described in the paper is a significant issue for accurate interpretation of HRV data and the adult period of REM sleep corresponds to the neonatal state of active sleep [38]. Siegel suggests that the phasic motor activation that characterizes REM sleep in the adult is also present in the

PAGE 88

76 neonate. The amplitude of the twitching, described during REM sleep, is more intense in neonates than in the adult, but a developmental continuity can be observed between "active sleep" in the neonate and REM sleep in the adult. Both active sleep in the neonate and REM sleep in the adult can be defined by purely behavioral criteria [38]. Therefore, the recordings should be taken in the presence of a medical professional qualified to identify phases of paradoxical sleep. An auditory stimulation that is widely familiar and contains no political or nationalistic propaganda should be used to prevent to prevent undesirable neural responses due to stress. The ECG recordings during the auditory stimulation can later be used to compare the COR of the neonate and adult. Creating the Index of ANS Maturation We suggest using the EMD method, detailed in Chapter 6, to compare the COR of the HRV groupings described in the previous subsection with the mature adult. The EMD data provides a time-based representation of the COR. The examiner should compare differences in the COR of the adult and neonate to determine an efficient method for relating these changes. We suggest exploring the use of power-ratios, neural networks, windowing functions, etc. as possible criteria to represent the EMD data in the form of an index. In addition, the development data can be compared using PCA and a neural network to classify the Lomb spectrum into different level of ANS maturation.

PAGE 89

APPENDIX A TIME-DOMAIN HRV ANALYSIS Week One 77

PAGE 90

78 Week Two

PAGE 91

79 Week Three

PAGE 92

80 Week Four

PAGE 93

81 Week Five

PAGE 94

82 Week Six

PAGE 95

83 Week Seven

PAGE 96

APPENDIX B FREQUENCY-DOMAIN HRV ANALYSIS Week One 84

PAGE 97

85 Week Two Week Three

PAGE 98

86 Week Four Week Five

PAGE 99

87 Week Six Week Seven

PAGE 100

LIST OF REFERENCES 1. BLISS, Causes of premature birth (2003), http://www.bliss.org.uk/support/causes.asp Accessed 4/4/2003. 2. Thomas, K., Causes of premature birth (2001), http://www.kerri.thomas.btinternet.co.uk/causes.htm Accessed 4/4/2003. 3. Chatow, U., Davidson, S., Reichman, B. & Akselrod, S. (1995). Development and maturation of theautonomic nervous system in premature and term infants using spectral analysis of heart rate fluctuations. Pediatric Research, 37, 294-302. 4. Stark, R., Myers, M., Daniel, S.,Garland, M., & Kim, Y. (1999). Gestational age related changes in cardiac dynamics of the fetal baboon. Early Human Development, 53, 219-237. 5. Boon, R., Learning Discoveries Psychological Services (2002), http://home.iprimus.com.au/rboon/HeartRhythmsandHRV.htm Accessed 1/12/2003. 6. Groome, L., Mooney, D., Holland, S., Smith, L., Atterbury, J., & Loizou, P. (1999). High Vagal Tone is Associated with More Effiecient Regulation of Homeostasis in Low-Risk Human Fetuses. Dev. Psychophysiology, 36, 25-34. 7. Krueger, C. (2002). Fetal heart rate variability and learning in the 28-34 week fetus. (Doctoral dissertation, University of North Carolina at Chapel Hill, 2001). Dissertation Abstracts International, 62, 11B. 8. Porges, S. (1995). Orienting in a defensive world: Mammalian modifications of our evolutionary heritage. A Polyvagal Theory. Psychophysiology., 32, 301-318. 9. Porges, S., Doussard-Roosevelt, J., Stifter, C., McClenny, B., & Riniolo, T. (1999). Sleep state vagal regulation of heart period patterns in the human newborn: An extension of the polyvagal theory. Psychophysiology, 36, 14-21. 10. Beauchaine, T. (2001). Vagal tone, development, and Grays motivational theory: Toward an integrated model of autonomic nervous system functioning in psychopathology. Development and Psychopathology, 13, 183-214. 11. Gagnon, R. (1992). Fetal behavioral states: Biological alteration. Seminars in Perinatology, 16, 234-238. 88

PAGE 101

89 12. Holditch-Davis, D. (1990). The development of sleeping and waking states in high-risk preterm infants. Infant Behavior and Development 13, 513-53 1. 13. Brackbill, Y. (1977). Behavioral state and temporal conditioning of heart rate in infants. In G. Etzel, J. LeBlanc, D. Baer (Eds.), New developments in behavioral research: Theory, method and application. Hillsdale, NJ: Erlbaum. 14. Kisilevsky, B. & Muir, D. (1984). Neonatal habituation and dishabituation to tactile stimulation during sleep. Developmental Psychology, 20, 367-373. 15. Thoman, E.B., & Whitney, M.P. (1990). Behavioral states in infants. Individual differences and individual analyses. In LJ. Colombo & J. Fagen (Eds.) Individual differences in infancy: Reliability, stability, and prediction, (pp. 113-136). Hillsdale, HJ: Erlbaum. 16. Arduini, D., Rizzo, G. & Romanini, C. (1995). Fetal behavioral states and behavioral transitions in normal and compromised fetuses. In J.-P. Lecanuet, W. Fifer, N. Krasnegor and W. Smotherman (Eds.), Fetal development: A psychobiological perspective, (pp. 83-100), Hillsdale, NJ: Erlbaum. 17. Faircloth, D., & Helms, R. (1986). A mixed model with linear covariance structure: A sensitivity analysis of the maximum likelihood estimates. Journal of Statistical Computer Simulations, 25, 205-206. 18. Gray, L. & Philbin, M. (2000). Measuring sound in hospital nurseries. Journal of Perinatology, 20(8), S100-S104. 19. Gottfried, A. & Hodgman, J. (1984). How intensive is newborn intensive care? Pediatrics, 74(8), 292-294. 20. Morris, B. Philbin, M. & Bose, C. (2000). Physiological effects of sound on the newborn. Journal of Perinatology, 20(8), S55-S60. 21. Gerhardt, K. (1989). Characteristics of the fetal sheep sound environment. Seminars in Perinatology, 13, 362-370. 22. Laguna, P., Caminal, P., Raimon, J., & Rix, H. (1991). Evaluation of HRV by PP and RR Interval Analysis Using a New Time Delay Estimate. IEEE Trans. Biomedical Eng., 63, 63-66. 23. American Heart Association (2002). Handbook of emergency cardiovascular care. AHA, 2002: Dallas, TX. 24. Pan, J., & Tompkins, W. (1985). Evaluation of HRV by PP and RR Interval Analysis Using a New Time Delay Estimate. IEEE Trans. Biomedical Eng., 32 n.3, 230-236.

PAGE 102

90 25. Niles, R., Stats lessons (2003), http://www.robertniles.com/stats/stdev.shtml Accessed 5/4/2003. 26. Malik, M. (1996). Heart rate variability: Standards of measurement, physiological interpretation, and clinical use. European Heart Journal, 17, 354-381. 27. Brennan, M., Palaniswami, M., & Kamen, P. (2002). Poincar plot interpretation using physiological model of HRV based on a network of oscillators. Am. J. Physiol. Heart. Circ. Physiol., 283, H1873-H1886. 28. Duarte, M. Software (1999), http://www.usp.br/eef/lob/md/software Accessed 6/4/2003. 29. Stein, P.K., & Kleiger, R.E. (1999). Insights from the study of heart rate variability. Annu. Rev. Med., 50, 249-261. 30. Laguna, P., Moody, G.B., & Mark, R.G. (1998). Power spectral density of unevenly sampled data by least-squares analysis: Performance and application to heart rate signals. IEEE Trans. Biomedical Eng., 45 n.6, 698-715. 31. Lomb N.R. (1976). Least-squares frequency analysis of unequally spaced data. Astrophysical and Space Science, 39, 447. 32. Glover, D.M. Sequence analysis II: Optimal filtering and spectral analysis (2000), http://w3eos.whoi.edu/12.747/notes/lect07/l07s05.html Accessed 5/20/2003. 33. Huang, N., Shen, Z., Long, S., Wu, M., Shih, H., Zheng, Q., Yen, N., Tung, C., & Liu, H. (1998). The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. Lond. A., 454, 903-995. 34. Rilling, G., Flandrin, P., & Gonalvs, P. Empirical mode decomposition (2002), http://www.ens-lyon.fr/~flandrin/emd.html Accessed March 25 th 2003. 35. Holditch-Davis, D., Edwards, L.J., &Helms, R.W. (1998). Modeling development of sleep-wake behaviors: I. Using the Mixed General Linear Model. Physiology & Behavior. 63 n.3, 6311-318. 36. Belozeroff, V., Berry R.B., Sassoon C.S.H., & Khoo M.C.K. (2002). Effects of CPAP therapy on cardiovascular variability in obstructive sleep apnea: A closed-loop analysis. Am. J. Physiol. Heart. Circ. Physiol., 282, H110-H121. 37. Lemons, J.A., & The American Academy of Pediatrics Committee on Fetus and Newborn (1998). Hospital discharge of the high-risk neonate: Proposed guidelines (RE9812). Pediatrics, 102 n.2, 411-417.

PAGE 103

91 38. Siegel, J (1999). The evolution of REM sleep. Handbook of Behavioral State Control. Lydic, R and Baghdoyan (Eds.), (pp. 87-100), Boca Raton, FL: CRC Press.

PAGE 104

BIOGRAPHICAL SKETCH Devin Lebrun received his bachelors degree in electrical engineering with a focus in communications from the University of Florida in 2001. He received his Master of Engineering in electrical engineering from the University of Florida in August 2003. His research interests include biomedical imaging, digital signal processing, and neural networks. 92