Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Noninvasive monitoring of human circadian phase using model-based particle filter estimation Mott, Christopher Grey 2006

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

Item Metadata

Download

Media
831-ubc_2006-0740.pdf [ 8.48MB ]
Metadata
JSON: 831-1.0064896.json
JSON-LD: 831-1.0064896-ld.json
RDF/XML (Pretty): 831-1.0064896-rdf.xml
RDF/JSON: 831-1.0064896-rdf.json
Turtle: 831-1.0064896-turtle.txt
N-Triples: 831-1.0064896-rdf-ntriples.txt
Original Record: 831-1.0064896-source.json
Full Text
831-1.0064896-fulltext.txt
Citation
831-1.0064896.ris

Full Text

Noninvasive Monitoring of H u m a n Circadian Phase using Model-Based Particle Filter Estimation by Christopher Grey Mott B.A.Sc, University of British Columbia, 2000 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Applied Science in THE FACULTY OF GRADUATE STUDIES (Electrical and Computer Engineering) The University of Brit ish Columbia December 2006 © Christopher Grey Mott, 2006 Abst rac t This thesis presents a novel method for monitoring the state of an individual's circadian physiology. While scientific understanding of circadian physiology has advanced rapidly in recent years, the difficulty of measuring an individual's circadian phase in a non-invasive and portable manner remains a barrier to applications in workplace safety and medical fields. Existing approaches to circadian monitoring use measurements from a single physiological marker such as core body temperature. A new modelling framework is developed here that allows the integration of multiple physiological measurements and a predictive model of cir-cadian physiology dynamics. On top of this framework a particle filter estimation algorithm is implemented that accurately captures the nonlinear dynamics of the circadian system and allows real-time incorporation of sensory input using Bayesian statistics. A human study in which subjects spent five days in an laboratory while continuously monitored with an array of sensors generated data from which the algorithm is demonstrated. Tuning parameters and the theoretical performance limits of the algorithm are also explored through a series of simulations. The novel conceptual aspects of the modelling framework wil l require further validation however results indicate that this system presents a set of capabilities leading towards application for circadian monitoring beyond the laboratory. ii Table of Contents Abstract ii Table of Contents iii List of Tables v List of Figures vi Acknowledgements x Dedication x i i 1 Introduction 1 1.1 Motivation and Goals 1 1.2 Thesis Organization 2 2 Background on Circadian Rhythms 3 2.1 Circadian Rhythm Physiology 3 2.1.1 Circadian Pacemaker Mechanism 3 2.1.2 Effects of Circadian Rhythms 4 2.2 Circadian Monitoring Applications 8 2.2.1 Alertness and Safety 8 2.2.2 Medical Treatment 9 3 Circadian Monitoring: State of the Art and Research Objectives 10 3.1 Measuring Circadian Phase 10 3.1.1 Constant Routine 11 3.1.2 Purification 12 3.1.3. Actigraphy 12 3.2 Modeling and Predicting Circadian Dynamics 13 3.2.1 Mathematical Models 13 3.2.2 Application of Model-Based Prediction 16 3.3 Objectives 18 4 Circadian Model Framework 20 4.1 Parameter Transformation of Circadian Pacemaker Model 21 4.1.1 Characteristics of existing states x and xc 22 4.1.2 Parameter transformation to new states <p and A 23 4.1.3 State-Space Model 30 4.2 Physiological Markers with Probability Distributions 32 4.3 Integrated model with multiple outputs 34 iii C O N T E N T S 5 Particle Filter Phase Estimator 37 5.1 Estimation Algorithm Selection 37 5.1.1 Recursive Bayesian Estimation Overview 38 5.1.2 Requirements for Circadian State-Space Model Estimation 40 5.1.3 Particle Filter Algorithm Selection 44 5.2 Particle Filter Design 45 5.2.1 Particle Filtering Algorithm 46 5.3 Probability Density Reconstruction 50 5.3.1 Probability Distribution Comparison Metrics 51 5.3.2 Kernel Density Estimation 52 6 Simulation Results 55 6.1 Scenario 1 55 6.2 Scenario 2 59 6.3 Scenario 3 65 6.4 Scenario 4 68 6.5 Summary 70 7 H u m a n Study Results 71 7.1 Equipment 71 7.1.1 Sensors 7 2 7.1.2 Laboratory 75 7.2 Protocol 75 7.3 Data Processing 77 7.3.1 Artifact Removal and Denoising 77 7.3.2 Particle Filter Configuration 79 7.4 Results 82 7.4.1 Cerebral Temperature Sensor 82 7.4.2 Particle Filter Estimation Test 84 7.5 Summary 87 8 Conclusion 92 8.1 Significance of Results ; 92 8.2 Summary of Contributions 93 Bibliography 95 A Research Ethics Board Certificates of Approval 103 L i s t of Tables 3.1 Example Illumination Levels 15 5.1 Filter Selection Comparison 45 6.1 Simulation 1 Parameter Value - Case with No Process Noise 56 6.2 Simulation 1 Parameter Value - Case With Process Noise 57 6.3 Simulation 2 Parameter Values 60 6.4 Simulation 3 Parameter Values 66 6.5 Simulation 4 Parameter Values 69 v L i s t of F igu res 2.1 Circadian pacemaker a n a t o m y '. 5 3.1 Trajectories of oscillator states x and xc i n the circadian pacemaker model . . . 15 3.2 Mathematical model o f the circadian pacemaker as proposed b y Kronauer et al 16 4.1 Original circadian physiology model with core body temperature o u t p u t 21 4.2 Typical values o f circadian states x and xc 23 4.3 Trajectory o f states x and xc i n an x-y plot 23 4.4 Values o f transformed state variables 8 and A over a t w e n t y four hour period . 24 4.5 Comparison of phase offset in simulation of constant and shifted l i g h t / d a r k schedules 26 4.6 Comparison of phase estimates calculated from x minima times and from phase offset 28 4.7 Comparison of original and transformed circadian pacemaker m o d e l 29 4.8 State-space propagation function and the output function 32 4.9 Example of point phase estimates from melatonin and C o r t i s o l 33 4.10 Example o f PDF estimates from melatonin and C o r t i s o l 34 4.11 Block diagram showing the generalized physiological marker m o d e l 35 4.12 Block diagram showing the general model f o r multiple physiological markers connecting t o a central circadian phase 35 vi LIST OF FIGURES 5.1 Recursive Bayesian estimation steps 40 5.2 Block diagram of the circadian phase estimator. 41 5.3 Parameter spaces of the state variables n, A, and <t> 42 5.4 Illustration of two light randomization methods 43 5.5 Gaussian vs. point mass probability representation 44 5.6 Particle filter diagram 47 5.7 Kernel probability density reconstruction 51 5.8 Bimodal probability distribution 52 5.9 Kullback-Liebler divergence with variable bandwidth and number of particles. 53 5.10 Illustration of bandwidth selection using MISE Kernel Density Estimation. . . 54 6.1 Light schedule for simulation 1, following a regular eight hour sleep schedule. 56 6.2 Particle trajectories for the three state variables during Simulation 1A 57 6.3 Comparison of multiple simulations with different numbers of particles 58 6.4 Light schedule for simulation 2 59 6.5 Particle probability distributions for the three state variables during Simula-tion 2A 61 6.6 Graphical probability density representation of the phase state variable during Simulation 2A 62 6.7 Particle probability distributions for the three state variables during Simula-tion 2B 62 6.8 Graphical probability density representation of the phase state variable during Simulation 2B 63 6.9 Particle probability distributions for the three state variables during Simula-tion 2C • 63 vii L I S T O F F I G U R E S 6.10 Graphical probability density representation of the phase state variable dur ing Simulat ion 2C 64 6.11 Comparison of final phase probability distributions after 21 days of known l ight level wi th different noise settings. . 64 6.12 Ligh t schedule for simulation 2 65 6.13 Particle trajectories for the three state variables dur ing Simulat ion 3 66 6.14 Graphical probability density representation of the phase state variable dur ing Simulat ion 3 67 6.15 L igh t schedule for simulation 4 68 6.16 Particle trajectories for the three state variables dur ing Simulat ion 4 69 6.17 Graphical probability density representation of the phase state variable dur ing Simulat ion 4 70 7.1 Summary of physiological sensors dur ing human study. 72 7.2 Cerebral temperature sensor diagram 74 7.3 Cerebral temperature sensor equipment picture 74 7.4 Laboratory controlled-environment room 75 7.5 Act iv i ty and l ight scheduling protocol for laboratory study. 76 7.6 Temperature data signal before (top) and after (bottom) artifact removal and denoising 78 7.7 Temperature data from subject 1 showing a l l sk in and core temperatures. The visible floor effect on the hand temperature channel is a result of a sensor range l imitat ion. This channel is not used i n an analysis 79 7.8 Configuration of particle filter estimator for processing data from subject 1. . . 81 7.9 Cerebral temperature sensor measurements and rectal core body temperature measurements plotted for each of the five subjects w i th major artifacts removed. 83 viii LIST OF FIGURES 7.10 Study protocol schedule for subject 1 84 7.11 Phase probability distribution generated only from known light levels 85 7.12 Study protocol schedule for subject 1 with raw and fitted H R V data superimposed. 86 7.13 Phase probability distribution generated from light levels and HRV measures . 87 7.14 Core body temperature and salivary melatonin 88 7.15 Physiological markers collected from Subject 1 during the constant routine . . 89 7.16 Physiological markers assessed during the constant routine and shifted to gen-erate estimates of central circadian phase 90 7.17 Comparison between phase estimates generated by constant routine (CR) core body temperature and melatonin measurements, and particle filter estimates . 91 Acknowledgements I extend my gratitude to Dr. Guy Dumont arid Dr. Mihai Huzmezan from The University of British Columbia who served as excellent supervisors and mentors through different stages of this project. Acknowledgement is due to Dr. Huzmezan for his role in proposing the col-laboration between U B C and the industrial partner Pulsar Informatics, and my thanks to him for providing a very warm welcome to my entry into graduate school. I felt privileged to be working with Dr. Dumont as he shared his experience with me, and whose support in navigating the challenges of starting a new research program and openness to opportunities that arose made this process very rewarding. Working in collaboration with Pulsar Informatics I was inspired by the innovation and ambitious spirit brought by Daniel Mollicone and Matthew van Wollen, and appreciative of their personal support which I frequently drew upon for motivation. The studies which we performed would not have been possible without the invaluable contributions of Dr. Diane Boivin at the Centre for the Study and Treatment of Circadian Rhythms and McGill University. My thanks are also extended to Dr. Boivin and her staff at Dougals Hospital for their hospitality during my stay in Montreal. An additional collaborator who I would like to thank is Dr. Andre Dittmar from the Laboratory of Materials Physics, INSA-Lyon, France. Financial support from the British Columbia Innovation Council, the National Science and Engineering Research Council of Canada (NSERC), and Pulsar Informatics Inc. is grate-fully acknowledged. x Acknowledgements A personal thank you is extended to Bronwyn Maxwell who I value dearly for her tremen-dous support and wisdom, friends at U B C who enriched the breadth of my experience, and of course to my wonderful family whose love continues to provide me a great strength. xi In memory ofJAunt (Dede xii Chapter 1 I n t roduc t ion Circadian rhythms are twenty-four hour oscillations in physiological processes that occur in most living organisms. In humans, circadian rhythms are observed in processes ranging from rates of cellular activity to our ability to perform cognitive performance tasks. This thesis presents a novel set of engineering methods for monitoring human circadian rhythm physiology. 1.1 M o t i v a t i o n a n d Goals When in a state of good health and maintaining a regular diurnal sleep schedule, our internal clock operates seamlessly in the background as it synchronizes our body's physiology with daily routines. However, when daily activities become desynchronized, such as occurs during shift work, or when acute medical care is being administered to a body on the edge of health, such as during chemotherapy treatment, the impact of our internally driven physiological oscillations becomes a factor necessary to consider in decision-making. The capability to accurately monitor and predict the state of an individual's circadian state is therefore of significant interest to a number of medical and workplace safety applications. Understanding of the specific impacts of circadian rhythms and the underlying physiolog-ical mechanisms have become increasing well known. The molecular basis of the clock has been specifically identified, and mathematical models characterizing the behaviour of the 1 Chapter 1. Introduction circadian pacemaker have been developed. Despite increased scientific understanding there are a number of significant challenges limiting the application of this knowledge beyond the laboratory. One is the lack of an effective method of measuring the internal state of an in-dividual's circadian physiology in real-time with non-invasive sensors that allow ambulatory movement. A second is that while there exist detailed mathematical models of circadian physiology with significant predictive capability, the complexity of these models limits the ease with which they can be practically applied. Motivated by the goal of enabling circadian rhythm science to be applied to health and human performance applications, this project was initiated with the aim of applying engineering modelling and signal processing techniques to the development of new circadian monitoring tools address these challenges. 1.2 T h e s i s O r g a n i z a t i o n As this thesis is an engineering project embedded in a field of human physiology, Chapter 2 first provides a background of the fundamental science of circadian rhythm physiology and discusses some of the motivating applications in more detail. Subsequently, Chapter 3 re-views the current state of the art in the area of circadian rhythm monitoring, and the specific technical objectives for new work are laid out in Chapter 3.3. Detailed development of a novel method for circadian monitoring is then presented in two parts: the development of a modelling framework in Chapter 4, and a particle filter estimation algorithm in Chapter 5. Analysis of the performance of the method is based on both simulations and on data from a human study which was conducted as part of the project. The simulated results are presented in Chapter 6, and the design of the human study and the results from its data are presented in Chapter 7. Finally Chapter 8 concludes with an analysis of the significance of the results and suggestions for future work. 2 Chapter 2 B a c k g r o u n d on C i r c a d i a n Rhy thms The word "circadian" is derived from the Latin words circa, meaning about, and dies, or day. Circadian physiological rhythms are present in organisms across the animal and plant kingdoms. This section provides a background on the physiology of circadian rhythms in hu-mans, and the applications in which monitoring of human circadian rhythms is of significant interest. 2.1 Circadian Rhythm Physiology Circadian rhythms are driven by an internal pacemaker which maintains a self-regulating oscillation with a twenty-four hour period. Recent research has revealed the molecular struc-tures which form the core of the human circadian pacemaker. The pacemaker serves as a central timing mechanism which synchronizes the rhythms of a wide array of physiological systems. 2.1.1 Circadian Pacemaker Mechanism Daily fluctuations in human physiology, such as sleeping and body temperature changes, have long been observed, however it was not clear in modern Western science whether these were internally driven or simply responses to environmental stimuli. It was not until the 1970s 3 Chapter 2. Background on Circadian Rhythms that strong experimental evidence of the existence in humans of an endogenous circadian pacemaker emerged. Subsequently, in the 1990s the specific molecular basis of a central human circadian pacemaker was identified. Located in the suprachiasmatic nucleus (SCN) in the hypothalamus region of the brain [25] a molecular clock maintains an approximate twenty four hour rhythm. While evidence of additional peripheral oscillators exists, such as in the liver, the SCN pacemaker is believed to play the central role in regulating circadian timing signals for other physiological systems. Although the pacemaker has an intrinsic period close to twenty four hours [12], precise synchronization to the external environment is maintained by external stimulii called "zeit-gebers"1. The strongest zeitgeber for most organisms, including humans, is light[13]. The daily transitions between light and dark caused by the earth's rotation in front of the sun create a strong environmental stimulus which organisms synchronize to. Synchronization of the circadian pacemaker to light occurs through photoreceptors in the retina which have a neural pathway to the SCN distinct from the visual system [26]. Signals arriving to the SCN modify both the phase and amplitude of the pacemaker's oscillations. The duration, intensity, timing relative to circadian phase, and pattern of light exposure are all factors in determining the specific impact. The location of the S C N and the light pathways are shown in Fig. 2.1. 2.1.2 Effects of Circadian Rhythms Studying the effects of circadian rhythms is not as straight forward as just looking for twenty-four hour physiological oscillations. Daily patterns of physical activity and sleep-wake also generally occur on a twenty-four hour schedule, so a distinction must be made between be-haviour induced rhythms (e.g. body temperature rising during the day because of walking) and endogenously driven rhythms (e.g. body temperature rising based on internal circadian 'From German Zeit, time and Geber, giver which is "time-giver". 4 Chapter 2. Background on Circadian Rhythms Suprachiasmatic Alertness Hemoglobin F i g u r e 2.1: C i r c a d i a n p a c e m a k e r a n a t o m y s h o w i n g t he l oca t i on o f s u p r a c h i a s m a t i c n u c l e u s ( S C N ) , the l i g h t p a t h w a y to t he S C N , a n d the c i r c a d i a n s i g n a l s w h i c h o r i g i n a t e at t h e S C N a n d a re d i s t r i b u t e d t h o u g h o u t t h e body. t h e r m o r e g u l a t o r y s i gna l s ) . T w o p r e d o m i n a n t e x p e r i m e n t a l t e c h n i q u e s for i s o l a t i n g t he c i r c a d i a n ef fects a r e t h e " fo rced d e s y n c h r o n y " a n d " c o n s t a n t r o u t i n e " protoco ls . B o t h occu r i n t i m e i s o l a t i o n l abo ra to r i es . T h e first forces a n i n d i v i d u a l ' s s leep a n d w a k e schedu les to d e s y n c h r o n i z e f r o m t h e i r i n t e r n a l c i r -c a d i a n p a c e m a k e r , t he second e l i m i n a t e s s l eep /wake effects by k e e p i n g i n d i v i d u a l s a w a k e i n a c o n s t a n t e n v i r o n m e n t fo r m o r e t h a n t w e n t y fou r h o u r s . B a s e d o n s t u d i e s c o n d u c t e d w i t h t hese pro toco ls , k n o w l e d g e of t he speci f ic r e l a t i o n s h i p s b e t w e e n the c i r c a d i a n p a c e m a k e r a n d v a r i o u s p h y s i o l o g i c a l s y s t e m s h a s been i den t i f i ed . A s u m m a r y of c i r c a d i a n ef fects o n some k e y p h y s i o l o g i c a l s y s t e m s i s p r o v i d e d below. 2.1.2.1 Core Body Temperature C o r e body t e m p e r a t u r e e x h i b i t s one o f the s t ronges t o b s e r v a b l e c i r c a d i a n r h y t h m . W h i l e h u m a n t e m p e r a t u r e r e g u l a t i o n i s c o m m o n l y t h o u g h t o f as t he h o m e o s t a t i c m a i n t e n a n c e o f a c o n s t a n t t e m p e r a t u r e , i n fac t o u r bod ies m a i n t a i n a cons tan t o s c i l l a t i o n w i t h a n a m p l i t u d e o f app rox 1 °C . T h e m i n i m u m o f the a p p r o x i m a t e l y s i n u s o i d a l t ra jec to ry occu rs a t a r o u n d 5 Chapter 2. Background on Circadian Rhythms 4 : 0 0 a m f o r m o s t peop le [33] , a n d t h e p rec i se t i m i n g of t h e m i n i m u m i s u s e d as a n i n d i c a t o r o f c i r c a d i a n p h a s e . T h e m o s t d e t a i l e d a n a l y s i s o f t he c i r c a d i a n ef fects o n core t e m p e r a t u r e i n c l u d e m a t h e m a t i c a l m o d e l s o f t h e c i r c a d i a n t e m p e r a t u r e r e g u l a t i o n a n d the i n t e r a c t i o n of s leep a n d p h y s i c a l a c t i v i t y [61. 2.1.2.2 Melatonin T h e c o n c e n t r a t i o n o f t h e h o r m o n e m e l a t o n i n p rov i des a v e r y r e l i a b l e m a r k e r o f c i r c a d i a n p h a s e e x p r e s s i o n u n d e r l o w l i g h t e x p e r i m e n t a l cond i t i ons . M e l a t o n i n i s a h o r m o n e assoc i -a t e d w i t h s l e e p i n g , a n d for m o s t o f t h e d a y m e l a t o n i n h a s neg l i b l e l eve l s , b u t p r e c e e d i n g the a n t i c i p a t e d b e d t i m e , t he re i s a n onset o f a m e l a t o n i n s p i k e . T h e t i m i n g o f t h i s s p i k e i s r e f e r r e d to as t h e D i m L i g h t M e l a t o n i n O n s e t ( D L M O ) t i m e . [53]. 2.1.2.3 Cortisol C o r t i s o l i s a h o r m o n e p r i m a r i l y r e s p o n s i b l e for r e g u l a t i n g m e t a b o l i s m a n d the body 's re -sponse to i n f l a m m a t i o n a n d s t ress . H u m a n s tud ies h a v e s h o w n t h a t sec re t i on o f C o r t i s o l h a s a s i g n i f i c a n t fluctuation w h i c h i s g o v e r n e d b y t h e c i r c a d i a n p a c e m a k e r [11] w i t h p e a k l eve l s o c c u r i n g b e t w e e n 8 :00am a n d 1 0 : 0 0 a m a n d d e c l i n i n g t h r o u g h o u t t he day. C o r t i s o l a c c o m p a -n i e s m e l a t o n i n a n d core b o d y t e m p e r a t u r e as t h e m o s t w i d e l y u s e d p h y s i o l o g i c a l m a r k e r s o f c i r c a d i a n p h a s e . 2.1.2.4 Cellular Activity A t a c e l l u l a r l e v e l , t he r a t e of ce l l p r o l i f e r a t i o n i n bone m a r r o w a n d h a s been i d e n t i f i e d to h a v e c i r c a d i a n v a r i a t i o n s [52]. A d d i t i o n a l l y t u m o r s a n d t u m o r - b e a r i n g hos t s e x h i b i t s t r o n g c i r c a d i a n r h y t h m s [41]. 2.1.2.5 Heart Rate Variability T h e h u m a n c a r d i a c r e g u l a t o r y s y s t e m r e s p o n d s to changes f r o m the r e s p i r a t o r y cyc le , t he r -m o r e g u l a t i o n , a n d m o s t s i g n i f i c a n t l y f r o m t h e a u t o n o m i c n e r v o u s s y s t e m . T h e r e i s deba te 6 Chapter 2. Background on Circadian Rhythms about the extent to which endogenous circadian rhythms impact cardiac measures. A cir-cadian variation, distinct from sleep-wake patterns, in parasympathetic activation has been observed under controlled conditions [36]. In studies of day and night shift operators however, twenty four hour cardiac variations are observed but primarily attributed to sleep/wake behaviour rather than an endogenous circadian component [18; 28]. These studies had a limitation in that they did not record circadian phase so it is possible that night shift workers just had their circadian phases aligned. 2.1.2.6 Respiration Circadian variation in the chemoreceptive respiratory feedback system has been shown to have a minima that occurs six to eight hours earlier than core body temperature minima [55]. Not all respiratory measures exhibit an endogenous circadian component, however detectable rhythms in tidal volume and respiratory frequency have been observed [55]. Under ambula-tory conditions respiratory rate does show a nocturnal decrease in rate [49] however this is likely due to a corresponding decrease metabolic rate related to lower activity levels. 2.1.2.7 Performance and Alertness The circadian pacemaker coordinates the endocrine and metabolic systems in such a way that the body is alert during daytime hours and promoting sleep during evening hours. Studies of cognitive performance in individuals who are kept awake for extended durations or on shifted schedules have revealed a strong circadian modulation which is independent from sleep schedules. During periods in which the internal circadian clock is triggering a sleep response individuals experience a negative dip in their alertness and performance [5; 31]. 7 Chapter 2. Background on Circadian Rhythms 2.2 Circadian Monitoring Applications Circadian rhythms significantly affect human health and safety, and there are a number of applications where knowledge of the exact phase of an individual's circadian phase would prove useful in mitigating risks or delivering therapeutic treatment. Although the phase of an individual's circadian rhythm can typically be correlated with the average wake-up time, inter-individual differences can be substantial. Furthermore, in situations where an individual has had irregular sleep schedule (such as shift work) or there has been some desynchrony between and individual's sleep/wake pattern and environmental light patterns (such as time zone travel), accurate estimates of circadian phase are difficult to obtain. A way of estimating circadian phase in real-time using ambu-latory, non-invasive methods would enable a new range of applications in two major fields: alertness and safety, and medical treatment. 2.2.1 Alertness and Safety Individuals required to perform critical tasks when their circadian pacemaker is coordinat-ing the body's metabolic and endocrine systems for sleep experience significant impairments in performance which increases the risk of accidents. M a n y industrial workplace scenarios involve tasks demanding a high degree of alertness and reliability in addition to requiring workers to operate on irregular or shifting schedules. As a result of schedule variations, in-dividuals experience reduced levels of alertness and cognitive performance due to both sleep loss and circadian rhythm desynchrony. The associated risk of accidents is a concern in op-erational settings, such as transportation [10], health care [60], and emergency response [51]. Accidents at the Chernobyl and Three Mile Island nuclear power plants, the N A S A Challenger disaster, and the Exxon oil spill in Valdez, Alaska, were all partially caused by decrements in performance due to the effects of circadian phase and length of time awake. Strategies to mitigate the adverse effects of circadian phase desynchrony and sleep loss 8 Chapter 2. Background on Circadian Rhythms in shift work environments include education programs, the design of shift work rotation schedules, and the use of light exposure to shift the endogenous circadian pacemaker [46]. Strategies to reduce the effects of jet lag travel also include the specific timing of light expo-sure to accelerate adaptation to a new time zone [4]. Implementing reliable means to shift an individual's circadian pacemaker however requires accurate knowledge of the individual's circadian phase, which is not currently possible outside of the laboratory. 2.2.2 Medical Treatment The traditional medical concept of "homeostasis", that the human body maintains a constant internal state, is giving way to the recognition of continuous time-varying fluctuations [9]. Applying this information, medical treatments are being developed to deliver treatment at not only the right location but at the right time. For instance, the cycles present in human circadian physiological systems bring about pre-dictable changes in the body's tolerance to anticancer agents and a tumour's responsiveness to them. Indeed, the tolerability and the efficacy of chemotherapeutic drugs can vary by 50% or more as a function of dosing time in mice or rats [42]. Circadian-modulation of continu-ous drug delivery systems for chemotherapy patients has been demonstrated in randomized, multi-centre trials [40; 44] to have enhanced tumor response and increased survival rates. Administering chronomodulated therapy requires synchronized timing to the endogenous circadian phase of the individual being treated. In recent trials, physical activity patterns measured with wrist worn accelerometers, have been used to infer sleep/wake schedules and thus circadian phase. While the activity measures can lead to a general estimate of cir-cadian phase i f an individual has maintained a consistent schedule, there are known vari-abilities between individuals [32] and uncertainty in the method. It can be expected that an enhanced method of accurately monitoring circadian phase would translate to increased efficacy of chronotherapies. 9 Chapter 3 C i r c a d i a n M o n i t o r i n g : State of the A r t and Resea rch Object ives This chapter describes the current methods for monitoring and predicting human circadian rhythm states. The existing methods are largely confined to use in laboratory settings, which leads to the primary goal of this work, which is to make steps towards enabling monitoring of circadian rhythms in non-laboratory conditions and thus enable use in the applications discussed in Sec. 2.2. The primary desired characteristics of a circadian monitoring method are that it allows ambulatory movement, has noninvasive sensors, and provides real-time feedback. 3 .1 Measuring Circadian Phase Since the central circadian pacemaker mechanism is inaccessibly located in the brain, its state cannot be measured directly and must be inferred from observations of downstream physiological effects. The arising complication is that systems with an observable circadian modulation, such as core body temperature and melatonin secretion, are also responsive to other physiological systems or environmental stimulus. From the perspective of circadian rhythm analysis, the additional signals are considered to mask the circadian component. A primary challenge in direct circadian measurement is thus to "demask" the circadian signal components from the other non-circadian components. 10 Chapter 3. Circadian Monitoring: State of the Art and Research Objectives Two demasking approaches are the physical el imination of a l l t ime-varying exogenous st imulus, and the removal of the exogenous factors by signal processing techniques. Hold ing a l l exogenous s t imul i constant is accomplished wi th a laboratory protocol called the "constant routine", and the signal processing methods are commonly referred to as "purification". Both are described i n further detail below 3.1.1 Constant Routine Core body temperature [7] and melatonin [53] are the two physiological systems that exhibit the most consistent and observable circadian rhythms, however core temperature also re-sponds to physical activity, posture, ambient temperature, and sleep. Melatonin secretion is also highly responsive to ambient l ight exposure. The constant routine (CR) procedure de-veloped by Czeisler [11] is designed to eliminate the majority of confounding effects on core temperature and melatonin, by placing subjects i n a strictly controlled laboratory environ-ment. To eliminate the effects of sleep-wake transitions and posture changes, the subjects are kept awake for a period of up to 40 hours i n a semi-recumbent position. L igh t exposure is set very low to 10 lux, meals are introduced at one hour intervals, and physical activity is l imi ted. D u r i n g the constant routine core body temperature is measured continuously, and the roughly sinusoidal circadian oscillation component wi th an amplitude of approximately 2°C is monitored. The t iming of the lowest point typical ly occurs between 4:00 a.m. and 5:00 a.m. and is the most frequently used indicator of phase. The natural circadian melatonin cycle consists of an onset i n secretion approximately at one's typical bed time. The t iming of the onset is driven by the circadian pacemaker, however melatonin secretion is also affected by exposure to ambient l ight. The d im light conditions of the constant routine allow accurate measurements of the D i m Ligh t Melatonin Onset ( D L M O ) t ime [8]. The constant routine is the gold standard method of assessing circadian phase and is the 11 Chapter 3. Circadian Monitoring: State of the Art and Research Objectives primary method by which data about the phase shifting effects of light were collected. Despite the success of the constant routine, its application is limited to laboratory environments as a research tool. 3.1.2 Purification The "purification" approach is another method of circadian phase estimation which attempts to remove the masking signals using signal analysis rather than the restrictive physical con-straints of the constant routine. Physical activity and sleep are the two primary masking factors of core body temperature so current purification methods have focused on the sep-aration of their effects from the circadian component. In contrast to the constant routine, participants in purification studies have been allowed to follow regular sleep/wake schedules with free ambulatory movement during waking periods. Waterhouse has developed statistical methods of purification U t i l i z i n g data from activity sensors. One method involves c a t e g o r i z i n g activity d u r i n g waking and sleep periods and then calculating an association temperature effect from each category[58]. A second method uses a linear regression based on direct mean scores from activity sensors[58]. Subsequent develop-ment of these purification methods has introduced some basic thermoregulatory models[59]. There is contention about the accuracy of purification approaches, however as arguments citing the unique accuracy of the constant routine[17] have been made on one hand, purifica-tion results comparable to constant routine results have also been shown[57]. A significant limitation of the statistical purification approach identified by Klerman et. al. [35] is that during periods of significant desynchrony between sleep-wake times and circadian phase, linear methods to separate the two effects from temperature data are inherently unreliable. 3.1.3 Actigraphy One method of circadian phase measurement which is currently used operationally outside laboratory conditions is the recording of rest-activity patterns using actigraphs. Actigraphs 12 Chapter 3. Circadian Monitoring: State of the Art and Research Objectives are sensors which record gross physical movement, most commonly with a wrist-worn ac-celerometer. This method assumes a direct correlation between an individual's restactivity rhythm and their sleep-wake rhythm and thus their circadian phase. Actigraphy has been used to the determine the circadian phase of cancer patients for timing the delivery of chronomodulated chemotherapy drugs [43]. The type of circadian vari-ation present in actigraph measurements has also been shown to provide an indicator of 'health status' of cancer patients. However, there has not been strong evidence demonstrating a precise link in humans be-tween actigraphy and more direct measures of circadian phase such as core body temperature or melatonin. Furthermore this approach has been applied only to individual's following a regular diurnal schedule. The enabling feature which has allowed the use of actigraphy in field applications is that it meets the criteria of non-invasiveness, portability. However confounding factors such as inter-individual variations in circadian phase, differences in behavioural patterns, and irreg-ular schedules such as arise with shiftwork limit its accuracy and precision. 3.2 Modeling and Predicting Circadian Dynamics In addition to direct measurement of an individual's circadian phase, predictive models of circadian pacemaker physiology are another set of tools which are currently applied to human circadian physiology analysis. Mathematical models describing the dynamic response of the circadian pacemaker have been used to predict the behaviour of the circadian pacemaker under specific light exposure scenarios. 3.2.1 Mathematical Models Studies of the human circadian pacemaker have led to the development of mathematical models at both a molecular and macro level that describe both the self-oscillating and light 13 Chapter 3. Circadian Monitoring: State of the Art and Research Objectives responsive dynamics. At the molecular level, light responsive gene transcription reactions have been characterized [24], and at the macro level, changes in the timing of human core body temperature rhythms in response to light exposure have been used to characterize the circadian system [13]. The macro level models were created with a black-box approach and no mathematical connections yet exist between the molecular level and macro level models. The macro-level models are the only type that can link to the other physiological systems of interest in this project, so molecular models are not considered further here. The most widely accepted model of the circadian pacemaker was developed by Kronauer et al in 1987 [37] based on macro observations of dose-response relationships between light exposure and circadian phase shifts. It was inferred from the data that the model should have both a self-regulating oscillator representing the internal circadian pacemaker, and a light input response system representing the pathway from-light in the eye to a synchronizing input on the oscillator. Subsequent discovery of the molecular functioning of the circadian pacemaker have supported the general physiological basis of the model. A refined version of the model was published in 1999 [30] and is subsequently used throughout this work referred to as the Kronauer model. In the Kronauer model, the self-sustaining rhythm of the circadian pacemaker is mod-elled by a modified Van der Pol oscillator which maintains a steady state oscillation with a stable period and amplitude. A light input term describes how light intensity observed in the retina causes shifts in the phase and amplitude. The Van der Pol oscillator consists of two interacting states described by the equations where fx - 0.13, q = 1/3, TX = 24.2, k - 0.55, and B is a driving input due to light [30]. The state variables x and xc follow trajectories which are close to two sinusoids out of phase by (3.2.2) (3.2.1) 14 Chapter 3. Circadian Monitoring: State of the Art and Research Objectives Table 3.1: Example Illumination Levels Illuminance (lux) a Example 0.00005 n/a Starlight 1 n/a Moonlight 10 0.003 Candle one foot away 400 0.02 Brightly lit office 400 0.02 Sunrise or sunset on a clear day 1,000 0.04 Typical TV studio lighting 3,000 0.08 Light therapy lamp (min.) 10,000 0.2 Light therapy lamp (max.) 32,000 0.3 Sunlight on an average day (min.) 100,000 0.7 Sunlight on an average day (max.) 9 0 ° as shown in Fig. 3.1. 1.5 0.5 0 -0.5 -1.5 0 8 16 24 32 40 48 56 64 72 Time [hours] Figure 3.1: Trajectories of oscillator states x and xc in the circadian pacemaker model during steady state conditions. The response (a) of the human eye to light is modelled first by a logarithmic function (3.2.3) where I is the ambient light intensity in units of lux, a0 = 0.05, and p=0.5 [30]. The log-arithmic function accounts for the very wide responsiveness range of the eye as shown by Table 3.1. The second component to the circadian light response is a dynamic filter which 15 Chapter 3. Circadian Monitoring: State of the Art and Research Objectives r e l a tes a to t he d r i v i n g i n p u t (B) on the o s c i l l a t o r s ta te v a r i a b l e s n = 60 [a ( l -n) -fin] B = Ga(\-n){\- mi)(1 - mxc) (3.2.4) (3.2.5) w h e r e /3=0.0075 a n d G=19 .875 [30]. E q u a t i o n (3.2.4) m o d e l s a f i l t e r (n) a c t i n g u p o n (a) a n d e q u a t i o n (3.2.5) m o d e l s t he m o d u l a t i o n of t h e d r i v i n g i n p u t (B) b y t h e c u r r e n t s ta te o f t h e c i r c a d i a n p a c e m a k e r a n d t h e f i l ter . F i n a l l y , t h e o u t p u t f r o m the m o d e l r e l a tes t h e s ta te v a r i a b l e x to a p h y s i o l o g i c a l c i r c a d i a n p h a s e m a r k e r . T h e t i m e a t w h i c h h u m a n core body t e m p e r a t u r e ( C B T ) r e a c h e s a m i n i m u m i s one of t he m o s t r e l i a b l e m a r k e r s o f c i r c a d i a n p h a s e a n d i s de f i ned as a re fe rence po in t . T h e t i m e o f t he m i n i m u m of v a r i a b l e x i s de f i ned to occu r 0.8 h o u r s before t h e t i m e of core body t e m p e r a t u r e m i n i m u m . CBTmin = x m i n + 0.8 h o u r s A d i a g r a m s h o w i n g the f u l l c i r c a d i a n p a c e m a k e r m o d e l i s s h o w n i n F i g . 3.2. C i r c a d i a n P a c e m a k e r Light (3.2.6) 0minjCBT} Tlinit Xmil Xc mil F i g u r e 3.2: M a t h e m a t i c a l m o d e l o f t he c i r c a d i a n p a c e m a k e r as p roposed b y K r o n a u e r et a l [30]. 3.2.2 Application of Model-Based Prediction T h e u s e f u l n e s s o f a d e s c r i p t i v e m o d e l i s d e p e n d e n t u p o n t h e degree to w h i c h i t c a n be u s e d to gene ra te m e a n i n g f u l p red i c t i ons . T h e K r o n a u e r c i r c a d i a n p a c e m a k e r m o d e l h a s been u s e d 16 Chapter 3. Circadian Monitoring: State of the Art and Research Objectives with differential equation solving algorithms to generate simulations predicting the phase shift of the circadian pacemaker, starting from a known initial condition, in response to a given light exposure pattern. This predictive capability has been successfully used to confirm the design of experimental protocols [4], and interpret experimental observations of circadian phase shifts [50]. Despite the apparent usefulness of the Kronauer model, it has not actually been widely applied in broader contexts. Current models of human alertness incorporate both a sleep-related component and a circadian component, however five out of the seven most widely used models assume a fixed circadian phase [45]. In those models circadian variation in alertness is simply modelled by a series of sinusoidal harmonics with a predetermined and constant phase [5]. Any scenarios in which the circadian phase may be non-stationary such as shift work or transmeridian travel cannot be accurately modelled. Of the two models that incorporate changing circadian phase, one uses the Kronauer model [31] and the other just uses a rule of thumb for shifting in response to a time-zone change [1]. In a comparative review of human performance models, the lack of circadian dynamic modelling in performance models was identified as a critical area for future development [14]. One of the challenges in applying a model with detailed circadian pacemaker dynamics to real scenario predictions is that current simulation methods require precise specification of initial conditions, and complete knowledge of light levels during the course of the simulation. In uncontrolled environments, such as in an actual workplace, it is difficult to assess both the circadian phase and ambient light levels for a specific individual. It may be this reason that the Kronauer model has found application in simulating laboratory environment scenarios, where circadian phase and light levels can be controlled, but has not been widely used in operational scenarios. This inability to apply circadian predictions to real environments is partially responsible for the fact that despite a well-established model of the circadian pace-maker, it can still be difficult for scientists to provide definitive advice concerning specific 17 Chapter 3. Circadian Monitoring: State of the Art and Research Objectives circadian adjustment countermeasures [46]. 3.3 Objectives The current state of the art in circadian monitoring provides reliable estimates of an indi-vidual's phase through the use of invasive sensors in controlled laboratory environments. Models of the human circadian pacemaker also provide the capability for predictions of circa-dian phase response to light exposure, but have not been widely applied and do not provide for measures of uncertainty. While existing models are useful, no systematic application of model-based signal process-ing or data fusion techniques to extract information for more than one measurement source have been conducted to date. Many of the methods also focus on statistical analysis of cir-cadian rhythms over the course of a full study, rather than on real-time estimation of an evolving state. Additionally, a model has been developed which describes circadian dynam-ics, but has primarily been used as an illustrative tool or to design experimental protocols. It has not been used in an integrated way to develop real-time measurements or for practical predictions. With the goal of making circadian rhythm monitoring feasible outside laboratory envi-ronments, a number of opportunities were identified for the application of new modelling and signal processing techniques, then, guided by the primary objective of this project - to develop tools to enable application of circadian rhythms science to practical uses in both workplace safety and medical treatment fields - the following technical objectives were set for this project: 1. Estimate an individual's circadian phase state in real-time by integrating existing mod-els of circadian pacemaker dynamics with sensor measurements; 2. Combine information from multiple sensors to detect circadian phase; 18 Chapter 3. Circadian Monitoring: State of the Art and Research Objectives 3 . Evaluate the potential of new, non-invasive, ambulatory sensors for measuring circa-dian phase; and 4. Present estimates and predictions of circadian phase in a statistical framework, allow-ing the use of probability distributions in scenario modeling and analysis. 19 Chapter 4 C i r c a d i a n M o d e l F r a m e w o r k In order to develop a circadian monitoring system capable of real-time estimation using inte-grated sensor data and model predictions, a coherent modelling framework is required as a foundation. With the aim of applying recursive estimation algorithms, the modelling frame-work should allow continuous tracking of an underlying circadian phase state, and define relationships between the system states and physiological measurements. To meet the objec-tive of allowing multiple sensor measurements, the framework should also define statistical relationships between multiple sensor measurements. The well-established circadian pacemaker model by Kronauer et al [30] provides a basis for developing a modelling approach, however it does not of itself meet all the criteria to serve as a general circadian modelling framework. As described in Sec. 3.2.1, the Kronauer model consists of a dynamic model of a circadian pacemaker with a correlation to a measurable property of core body temperature. Its elements can be grouped into a circadian pacemaker component and a physiological marker component as shown in Fig. 4.1. Three limitations preventing this model's use as a general modelling framework are that the circadian phase is not presented as a continuous variable which can be monitored or updated, the correlation between the circadian phase state and the physiological marker is not defined in a statistical manner, and it only specifies a specific correlation to core body temperature. Using the Kronauer model as a starting point, a general modelling framework 20 Chapter 4. Circadian Model Framework C i r c a d i a n P a c e m a k e r O s c i l l a t o r Light (pmmjCBT} riniil Xinit Xc 'mil C i r c a d i a n P a c e m a k e r P h y s i o l o g i c a l M a r k e r Figure 4.1: Original Kronauer model [30] of the circadian pacemaker, with input light model, state transition model, and definition of output phase estimate from core body temperature (CBT) marker with a phase shift of 0.8h. is developed by resolving each of the identified limitations. First, the Kronauer circadian pacemaker model is augmented with a parameter transformation to create continuous phase and amplitude states. Second, the physiological marker system is generalized to provide circadian phase probability distributions rather than point estimates. Finally a circadian physiology architecture is designed to correlate multiple physiological markers to a central circadian phase state. 4.1 Parameter Transformation of Circadian Pacemaker Model The Kronauer circadian pacemaker model, described in detail in Ch. 2, is the most accepted physiological-level characterization of circadian pacemaker phase and amplitude dynamics. This model is therefore the preferred choice for use in a state estimation, however phase and amplitude are not accessible as state variables. The amplitude is defined as a nonlinear function of the state variables x and xc and the phase is defined as the time of the nadir (local minimum) of xc. Thus the location of the nadir provides the reference point for connecting the model's phase to the phase of observable physiological systems. Most commonly, the xc 21 Chapter 4. Circadian Model Framework nadir is connected to the nadir of core body temperature as described in Eq. 3.2.6. To overcome this limitation a state transformation is proposed to map x and xc into di-rectly useful states <f> and A. To develop this transformation the properties of the existing state variables are examined, a state transformation is derived, and an analysis of the char-acteristics of the new states is completed. 4.1.1 Characteristics of existing states x and xc The differential equations used in the Kronauer model contain three states x, xc, and n. The variables x and xc interact to create a modified van der Pol oscillator which self-oscillates at a period of approximately twenty-four hours. Both x and xc follow nearly sinusoidal trajecto-ries, in which xc is phase lagging by 90°. The variable n is a state which is part of a light filter system as described in Sec. 3.2.1. The reference point that Kronauer et al. used to define the phase of the circadian pacemaker model is the timing of the nadir (minimum point) of x. For example, for an individual with sleep occurring regularly between 12:00am and 8:00am, sim-ulations of the Kronauer model will show that the nadir of x occurs around 4.3h (4:22am), and in real experiments subjects will have a body temperature minimum occurring at 5.1h (5:06am) as described in Eq. 3.2.6. Based on a variety of experiments with different light exposure settings, the model has been refined so the nadir timing of xc maintains a constant timing shift relative to observed core body temperature nadirs. The nadir phase reference method, however, presents two limitations in the context of a circadian estimation system. The first is that the measurement of the model's phase can only occur at a single point every twenty four hours. The second limitation is the lack of an inverse relationship from which the system's state variables can be updated to match a given phase and amplitude. Currently, the initial conditions for model simulations are just set based on a look-up table of values for a typical situations (e.g. for a habitual schedule of 8h sleep and 16h awake in 160 lux the initial conditions at bedtime are x = -0.17, x c = —1.22, 22 Chapter 4. Circadian Model Framework a n d n = 0.50 [30]). 4.1.2 Parameter transformation to new states 0 and A B a s e d on a n u n d e r s t a n d i n g of t h e e x i s t i n g m o d e l s ta tes x a n d xc a t r a n s f o r m a t i o n i s de-v e l o p e d to c rea te n e w p h a s e a n d a m p l i t u d e s ta tes t h a t c a n be c o n t i n u o u s l y m o n i t o r e d a n d u p d a t e d con t i nuous l y . T h e s ta tes x a n d xc f o l l ow n e a r s i n u s o i d a l t ra jec to r ies w i t h a c o n s t a n t p h a s e d i f fe rence , w i t h a t y p i c a l e x a m p l e s h o w n i n F i g . 4 .2 . I f x a n d xc a re c o n s i d e r e d as 1.5 1 0.5 -o -' S -0.5 --1 --1.5 ^ 4 12 Time [h] 16 20 24 F i g u r e 4 .2 : T y p i c a l v a l u e s o f c i r c a d i a n s ta tes x (b lack) a n d xc (grey) ove r a t w e n t y fou r h o u r p e r i o d , w i t h one -hou r i n t e r v a l s i n d i c a t e d b y s o l i d dots . A e x a m p l e set o f po i n t s i s h i g h l i g h t e d a t t = 15h. x -y coo rd ina tes on a C a r t e s i a n p l a n e t h e y desc r i be a n e a r - c i r c u l a r l ocus . F i g . 4.3 s h o w s the s a m e d a t a po i n t s f r o m F i g . 4.2 t r a n s p o s e d i n t o x - y coo rd ina tes . G i v e n the c i r c u l a r x - y p lo t , x 0 F i g u r e 4 .3 : T r a j e c t o r y o f s ta tes x a n d xc i n a n x - y p lo t s h o w i n g the t r a n s f o r m a t i o n to po la r coo rd ina tes w i t h v a r i a b l e s 6 a n d A. T h e e x a m p l e p o i n t a t t = 15ft, i s a lso h i g h l i g h t e d for c o m p a r i s o n to F i g . 4 .2 , i s s h o w n he re . 23 Chapter 4. Circadian Model Framework the trajectory can then be described naturally in polar coordinate system using amplitude, A, and phase angle, 9. Considering the amplitude, A, as the distance of the vector from the origin to the point (x, xc), transformation from x and xc to amplitude, A, is given by y/xTTxl (4.1.1) as described in [30]. To calculate an instantaneous phase angle, the vector from the origin to (x, xc) is considered to form an angle, 9, from the horizontal axis. The equation for 9 is given by 0 i f xc > 0 and x > 0 (quadrant I), 0 = tan_1(^) + < 7 r i f x c < 0 (quadrant II or III), (4.1.2) 2VT i f xc > 0 and x < 0 (quadrant IV), such that as time increases, the point (x, xc) follows a circular trajectory with a period of twenty-four hours and 9 increases from 0 to 2TT radians. The inverse transformation, T - 1 , is accomplished with x = Asin(9) (4.1.3) x c = ;4.cos(0). (4.1.4) In practice, it is most meaningful to describe the circadian system in terms of a phase offset 1.5 1 0.5 0 360 . ; . . . , . . (D — ^ . < i i i 12 Time [h] Figure 4.4: Values of transformed state variables 9 and A over a twenty four hour period, using the same data with the point at time t = I5h highlighted. 24 Chapter 4. Circadian Model Framework which reflects the phase shift relative to defined reference phase, rather than in terms of a phase angle, 9, which reflects the continually moving angle from 0 to 2ir. Denoting the phase offset as 0 in units of hours, the relationship could be described in general as 2TT e = uj(t + to)-<j>— 24 <b = — {UJ (t + 1 0 ) - 9) (4.1.5) where ui is the baseline frequency, and t0 is a baseline offset parameter. Thus the phase offset, </> is the difference between the phase angle and a reference line. Since the baseline frequency corresponds to a constant twenty-four hour day UJ = 2TT/24 then Eq. 4.1.5 can be rewritten as 24 <t> = t +10 - —9 (4.1.6) To illustrate this equation, i f an individual follows a consistent twenty-four light exposure schedule then their circadian phase angle, #24/(2?T), will increase at the same rate as t, so the difference will be constant. Alternatively, i f an individual experiences a shift in light exposure, then their circadian phase angle will increase or decrease relative to t, and a change in the phase offset will occur. Examples of these two situations are shown in Fig. 4.5, with a constant phase difference on the left, and a shifting phase offset on the right. The baseline offset term, to, in Eq. 4.1.6 can be treated as a calibration constant, as it shifts the mean value of the phase offset <f>. To determine the appropriate calibration value, a comparison is made to the original phase offset assessment method defined as the times of the minima of x. For consistency, the new phase assessment of <p should be calibrated to match the values which would be obtained by the original method. Since the original x m i n method provides only a point measurement every twenty-four hours while the new method with d> provides continuous values, there could be multiple calibration methods. The two calibration possibilities considered were to adjust the baseline offset of <f> so that: • at the points where x m i n is measured, the values of <j> match the phase estimates from x m i n exactly; or 25 Chapter 4. Circadian Model Framework 0 24 48 72 96 120 24 •» 20 | 16 £ 12 48 72 Time [h] 24 48 72 96 120 24 48 72 96 120 24 20 | 16 0) 2 12 O 48 72 Time [h] F i g u r e 4 .5 : C o m p a r i s o n of c i r c a d i a n s ta te v a r i a b l e s x a n d xc ( top row) , p h a s e 6 ( m i d d l e row) , a n d p h a s e offset <j> (bo t tom row) i n t w o cases. F o r a case w i t h a c o n s t a n t s l e e p / w a k e s c h e d u l e ( left s ide) t he p h a s e offset r e m a i n s cons tan t , a n d for a case w i t h s l e e p / w a k e s c h e d u l e s h i f t e d f o r w a r d b y s i x h o u r s ( r igh t s ide) t he p h a s e offset i n c r e a s e s as t h e c i r c a d i a n p a c e m a k e r e n -t r a i n s to t h e n e w l i g h t schedu le . • t h e t w e n t y - f o u r h o u r m e a n o f <j> m a t c h e s the p h a s e e s t i m a t e s f r o m x m i n . W h i l e a r g u m e n t s c o u l d be m a d e for b o t h , t h e c a l i b r a t i o n to t h e m e a n w a s c h o s e n as t h e p r e f e r r e d m e t h o d i n o r d e r to e l i m i n a t e a b i a s i n t h e a v e r a g e e s t i m a t e s . U s i n g t h e t w e n t y - f o u r h o u r m e a n o f <j> a spec i f i c c a l i b r a t i o n v a l u e for t0 w a s d e t e r m i n e d f r o m a s i m u l a t i o n o f m o d e l u n d e r c o n d i t i o n s w h i c h a re c o n s i d e r e d as s t a n d a r d s c h e d u l e a n d l i g h t exposu re scenar io . R e p e a t e d d a y s w i t h a n e igh t h o u r s leep ep isode (0 L u x ) f o l l owed b y a s i x t e e n h o u r a w a k e ep isode w i t h c o n s t a n t l i g h t l eve l o f 150 L u x , w e r e s i m u l a t e d u n t i l t h e m o d e l a c h i e v e d a s teady s ta te cond i t i ons . U n d e r t h i s c o n d i t i o n , t h e v a l u e o f t0 w h i c h of fsets t h e m e a n of <j> to m a t c h t h e p h a s e e s t i m a t e s f r o m x m i n w a s 17.1. T h e r e s u l t i n g c a l i b r a t e d 26 Chapter 4. Circadian Model Framework v e r s i o n o f E q . 4.1.6 i s 24 <j) = t + 17.1 - — 6. (4.1.7) Z7T w h e r e 0 i s g i v e n by E q . 4 .1 .2 . T h e c a l i b r a t e d p h a s e e s t i m a t e s f r o m t h e s teady s ta te c o n d i t i o n d e s c r i b e d above a r e s h o w n i n t he l e f t - h a n d p lo ts o f F i g . 4 .6 . T h e c o n t i n u o u s p h a s e e s t i m a t e s f r o m </> e x h i b i t a v a r i a t i o n o f abou t ± 0 . 5 h d u r i n g e a c h t w e n t y - f o u r h o u r p e r i o d . T h i s f l u c t u a t i o n i s c l e a r l y v i s i b l e i n t h e b o t t o m p lo ts of F i g . 4 .6 w h i c h s h o w the d i f fe rence b e t w e e n <fi a n d t h e s i ng l e p o i n t es t i -m a t e s from x m i n . T h e cause of the f l u c t u a t i o n i s t h a t t h e V a n de r P o l e q u a t i o n s ( E q . 3 .2 .1 , a n d 3.2.2) do no t fo l l ow a per fec t s i n u s o i d a l t ra jec tory . A s d i s c u s s e d s u b s e q u e n t l y i n C h . 5 , t h e m a g n i t u d e o f t h e f l u c tua t i ons i s s m a l l e n o u g h to a l l o w a c c u r a t e u s e i n a p h a s e e s t i m a t i o n a l g o r i t h m . To i l l u s t r a t e t h e p e r f o r m a n c e o f t he c o n t i n u o u s p h a s e e s t i m a t i o n i n a case w i t h a c i r c a d i a n p h a s e sh i f t , a s i m u l a t i o n w a s r u n i n w h i c h a n i n d i v i d u a l ' s s leep i s d e l a y e d b y s i x h o u r s f r o m a n e n t r a i n e d r h y t h m . T h i s case i s s h o w n i n the r i g h t - h a n d p lo ts o f F i g . 4 .6 . A s i x h o u r sh i f t i n t he c i r c a d i a n p h a s e as i t e n t r a i n s to t he n e w l i g h t / d a r k s c h e d u l e i s e v i d e n t i n t h e m i d d l e r i g h t p lo t . D u r i n g t h e cou rse o f t h e sh i f t b e t w e e n e n t r a i n e d c o n d i t i o n s , i t i s n o t e d t h a t cf> h a s a pos i t i ve b i a s w i t h respec t to t h e p h a s e a s s e s s m e n t f r o m x m i n . T h i s b i a s i s h i g h l i g h t e d i n t h e b o t t o m r i g h t p lo t o f F i g . 4 .6 . S i n c e 4> r e t u r n s to a n u n b i a s e d s ta te d u r i n g t h e e n t r a i n e d pe r i ods , t h e b i a s d u r i n g t he n o n - e n t r a i n e d pe r i ods i s a s m a l l t r a n s i e n t effect t h a t does no t s i g n i f i c a n t l y af fect s u b s e q u e n t u s e i n a p h a s e e s t i m a t i o n a l g o r i t h m . I n c o n c l u s i o n , one- to-one t r a n s f o r m a t i o n s b e t w e e n t he o r i g i n a l s ta tes i n t he K r o n a u e r m o d e l [30], x a n d xc, a n d t h e n e w s ta tes , a m p l i t u d e A a n d p h a s e of fset <j>, h a v e been de-v e l o p e d . W h i l e p h a s e m e a s u r e m e n t w a s poss ib le a t p o i n t v a l u e s e v e r y t w e n t y - f o u r h o u r s w i t h x m i n , t h e n e w s ta te v a r i a b l e cf> a l l o w s c o n t i n u o u s p h a s e m e a s u r e m e n t . W i t h t h e n o n -l i n e a r t r a n s f e r f u n c t i o n a p p e n d e d to t h e e x i s t i n g m o d e l t h e o r i g i n a l d y n a m i c s a re r e t a i n e d , 27 Chapter 4. Circadian Model Framework 0 24 48 72 96 120 144 168 0 24 48 72 96 120 144 168 0 24 48 72 96 120 144 • 168 0 24 48 72 96 120 144 168 Time[h] Time[h] F i g u r e 4 .6 : C o m p a r i s o n o f p h a s e e s t i m a t e s c a l c u l a t e d f r o m x m i n i m a t i m e s a n d f r o m p h a s e offset. b u t p h a s e a n d a m p l i t u d e c a n be t r e a t e d d i r e c t l y as c o n t i n u o u s s ta te v a r i a b l e s . T h e o r i g i n a l c i r c a d i a n m o d e l a n d the t r a n s f o r m e d c i r c a d i a n m o d e l a re s h o w n i n F i g . 4 .7 , i n t h e top a n d b o t t o m respec t i ve ly . T h e f o r w a r d t r a n s f o r m a t i o n b lock T i n F i g . 4.7 cons i s t s o f E q . 4.1.1 a n d 4 .1 .7 , a n d t h e i n v e r s e t r a n s f o r m a t i o n b lock T _ 1 cons is t s o f E q . 4 .1.3 a n d 4.1 .4 . T h e t r a n s -f o r m e d m o d e l c a n be d i r e c t l y i n c o r p o r a t e d i n t o a r e c u r s i v e e s t i m a t i o n a l g o r i t h m i n w h i c h t h e p h a s e a n d a m p l i t u d e a r e t r e a t e d as the s y s t e m s ta tes . 28 Chapter 4. Circadian Model Framework a) Light C i r c a d i a n P a c e m a k e r O s c i l l a t o r limit Xinit Xc mil b) Light C i r c a d i a n P a c e m a k e r L i g h t F i l t e r O s c i l l a t o r Tiinit l\ h B ri K 1 ' \\ 1 n X Xinit Xc Xc inii (pinit Ah Figure 4.7: Comparison of the the original Kronauer model of the circadian pacemaker (a) and the modified circadian pacemaker model with transformations to continuous amplitude, A, and phase offset, <j> states (b). 29 Chapter 4. Circadian Model Framework 4.1.3 S ta te -Space M o d e l To enable subsequent use of the transformed Kronauer model in estimation algorithms, it is helpful to cast the model into a state-space form. State-space models are defined by a state vector that contains the time varying properties of the system, a state transition function that describes how the state vector evolves in time, and an output function that describes how the states can be observed. For the case of a discrete time system the states evolve at incremental time intervals with sampling period T. Where x denotes the state vector we can use the notation x f e = x(t) where t = t0 + kT, k e N (4.1.8) . To create a state-space model of the circadian pacemaker a state vector is defined with the three variables used in the transformed Kronauer model A I (4.1.9) where A is circadian amplitude, ri> is circadian phase, and n is the light filter state. State-space models contain two functions, a state transition function and an output function. The state transition function can be generally written for discrete-time systems as x f e + 1 = f(xjt,u f c, v f c) (4.1.10) where x f e is the current state vector at time step k, Ufc is the measured input, v f e is the unmeasured input, and x f e + i is the predicted state vector at the future time step k + 1. The unmeasured input v f c is often called "process noise" and modelled as a random variable. To cast the Kronauer model in this form we can first see that the transformed model as shown in F i g 4.7 describes the transition function x^+i = ffc(xfc,Ifc) where I is the light input and x f c is given by E q . 4.1.9. To complete the generalized to a state space form the process noise v remains to be added. There is no precedent for specifying process noise however however 30 Chapter 4. Circadian Model Framework as a state-space formulation of this circadian model has not previously been developed. A choice made here is to assume that v is an additive Gaussian noise which is the predominant process noise model used in many other systems. Therefore, the state-space state transition equation for the transformed Kronauer model is Xfc+1 f(Xfe.Ifc) + Vfc (4.1.11) A fi(xfe,Ife) vl <f> = f 2(x f e,I f e) + v2 n Mxfc.Ifc) v3 The general form of a state-space output function is z f c = h f e(x f c ; w f c) (4.1.12) where w i s a random variable considered the "measurement noise". The output function for the circadian model is chosen based on the states to which measurement information could potentially be correlated. Of the three states, A, cj>, and n, the phase d> is the key state which will be observed. Therefore <j> is chosen as the single system output using the linear transfer function given by Zfc = h f c(xfc,Wfc) Hx f e + Wfc 0 1 0 (4.1.13) + Wfc + Wfc where w is an independent random Gaussian variable with variance R. A diagram illus-trating the functional relationship between the state propagation function and the output function is shown in Fig. 4.8. In situations in which circadian amplitude is being measured, this output matrix could be modified so that both phase and amplitude are outputs. 31 Chapter 4. Circadian Model Framework X 0 . f ( x , i ) h ( x ) - D e l a y +-O u t p u t S t a t e P r o p a g a t i o n F i g u r e 4 .8 : B l o c k d i a g r a m s h o w i n g t h e i n t e r c o n n e c t i o n o f t he s ta te -space s ta te p r o p a g a t i o n f u n c t i o n a n d t h e o u t p u t f u n c t i o n . 4.2 Physiological Markers with Probability Distributions T h e m e a s u r a b l e p h y s i o l o g i c a l m a r k e r o f c i r c a d i a n p h a s e i n t h e K r o n a u e r m o d e l i s t he t i m e a t w h i c h core body t e m p e r a t u r e r e a c h e s a m i n i m u m v a l u e (nad i r ) . T h e c e n t r a l c i r c a d i a n p h a s e i s t h e n d e f i n e d r e l a t i v e to t h a t po in t [30]. T h e r e a r e a v a r i e t y o f d i f fe ren t p h y s i o l o g i c a l s y s t e m s w h i c h i n d i c a t e c i r c a d i a n p h a s e , a n d t h e y e a c h h a v e c o r r e s p o n d i n g f ea tu re de tec t i on a l g o r i t h m s . F o r core b o d y t e m p e r a t u r e the m a r k e r i s t he m i n i m a . F o r o the r s y s t e m s t h i s c o u l d be t h e m a x i m u m , a t h r e s h o l d c r o s s i n g , o r some o t h e r a l g o r i t h m . C o m m o n a m o n g a l l t h e p h y s i o l o g i c a l m a r k e r s i s t h a t t h e y a re t y p i c a l l y on l y u s e d to p r o v i d e p o i n t e s t i m a t e s . A s a n e x a m p l e w e c a n l ook a t t he m e l a t o n i n a n d C o r t i s o l s a m p l e s t a k e n f r o m one i n d i -v i d u a l d u r i n g a l a b o r a t o r y s tudy. M e l a t o n i n a n d C o r t i s o l c a n b o t h be a n a l y z e d for a p h a s e m a r k e r b y f i t t i n g t h e d a t a to a F o u r i e r ser ies a n d i d e n t i f y i n g t h e t i m e a t w h i c h the m a x i m a occurs . S h o w n i n F i g . 4.9 a r e t h e C o r t i s o l a n d m e l a t o n i n d a t a co l l ec ted . T y p i c a l l y one o f t h e p h a s e re fe rence p o i n t s g e n e r a t e d f r o m the m e l a t o n i n a n d C o r t i s o l d a t a w o u l d be t a k e n as t h e s i ng l e s t r o n g p h a s e i n d i c a t o r for a n i n d i v i d u a l . S t a t i s t i c a l a n a l y s i s c o u l d t h e n be done c o m p a r i n g t h e i n te r - sub jec t v a r i a b i l i t y o r o the r p o p u l a t i o n d i s t r i b u t i o n m e t r i c s [34]. H o w e v e r , for i n t e g r a t i o n i n t o a m u l t i - p a r a m e t e r e s t i m a t i o n a l g o r i t h m w h e r e t h e i n t e r e s t i s o n t h e s ta te o f a g i v e n i n d i v i d u a l ' s c i r c a d i a n p h a s e a t a spec i f ic t i m e a p r o b a b i l i t y d i s t r i -b u t i o n i s s i g n i f i c a n t l y m o r e u s e f u l t h a t a po in t e s t i m a t e . A B a y e s i a n p r o b a b i l i t y d i s t r i b u t i o n 32 Chapter 4. Circadian Model Framework - © — Melatonin measurements 1 Fourier2 fit - G — Cortisol measurements 1 Fourier2 fit 16 20 0 4 Clock Time (hours) 16 20 0 4 Clock Time (hours) Figure 4.9: Example melatonin and C o r t i s o l measurements taken from one subject and used as circadian physiological phase markers with references points defined as point estimates (dashed arrows) at the peaks of the fitted curves. associated with a measurement indicates a degree of belief in the accuracy of the measure-ment. In a recursive estimation algorithm, the degree of belief from a new measurement can indicate to what degree the overall state estimation should be updated. As part of the modelling framework for the circadian monitoring system, physiological marker feature detection algorithms are augmented to generate probability distributions rather than point estimates. Specific implementation depends on the nature of the algorithm so there is not a single method, but the common case of locating the maxima or minima of a fitted Fourier is discussed here. With a first order Fourier fit (a single sinusoid) using least-mean-squares, the phase parameter directly indicates peak location, and a confidence bound on the phase estimate can be directly calculated. For higher order Fourier fits, the likeli-hood of the peak location is not linearly related to the Fourier parameters, so non-parametric approaches such as a grid search algorithm [56] can be applied to generate probability distri-butions for the peak locations. Reexamining the melatonin and C o r t i s o l example, the point estimates are replaced with probability distribution functions of the peak location generated using a grid search as shown in Fig. 4.10. The width and shape of the probability distributions is based upon both the peak amplitude and the degree of noise in the measured data. Data sets with smaller amplitude 33 Chapter 4. Circadian Model Framework Clock Time (hours) Clock Time (hours) Figure 4.10: Example melatonin and C o r t i s o l measurements taken from one subject with probability distributions for the peak locations of Fourier fits used as circadian physiological phase markers. peaks and higher noise will have wider probability distributions convey an increased uncer-tainty in the location of the phase marker. 4.3 Integrated model with multiple outputs With the objective to integrate multiple physiological sensors, the final piece of the modelling framework is a method to connect the data from multiple physiological markers to a common circadian phase estimate. The Kronauer model specifies an explicit connection between the circadian pacemaker phase and a single physiological marker of core body temperature, so it serves as a starting point from which a generalization is needed. The physiological marker component in the Kronauer model (refer back to Fig. 4.1) con-sists a time shift of 0.8h and a minimum point feature detection algorithm. This can be 34 Chapter 4. Circadian Model Framework abstracted to the general components shown in Fig. 4.11. CBT M a r k e r C B T S h i f t TCIJT ' . '- - - - • . . . P h y s i o l o g i c a l P h a s e E s t i m a t e Figure 4.11: Block diagram showing the generalized physiological marker model. Allowing arbitrary feature detection and time shifts for any physiological function, we can make a conceptual model which connects multiple physiological markers to a central circa-dian phase as shown in Fig. 4.12 The additional marker could include core body temperature, H u m a n C i r c a d i a n P a c e m a k e r 0 C B T CBT M a r k e r , . S h i f t T(;BT P n m„ M a r k e r , S h i f t 0 i J P h y s i o l o g i c a l P h a s e E s t i m a t e Figure 4.12: Block diagram showing the general model for multiple physiological markers connecting to a central circadian phase. melatonin, and any other system which has a measurable phase, and a known phase shift to the central pacemaker. As there are no known mechanisms for measuring a 'true' circadian phase from the central pacemaker in the SCN, there is no absolute reference value for the circadian phase state variable <f>. Therefore <j> can be arbitrarily calibrated. For the purposes of the subsequent work, the central phase <j> will be defined with a phase equal to that given by the minima of the x variable of the circadian model. Thus the calibration described in Sec. 4.1.2 will be 35 Chapter 4. Circadian Model Framework maintained consistently. Note that this modelling structure does not necessarily suggest internal physical mecha-nisms in the human body, but rather it just provides a means to describe the observable phys-iological behaviour. The integrated connection between multi-parameter output and a single pacemaker sets the necessary foundation for implementing a recursive state estimation al-gorithm, which simultaneously utilizes multiple measurements and central state transition updates. 36 Chapter 5 Par t i c le F i l t e r Phase Es t ima to r In C h . 4 a modelling framework was developed to describe the dynamics of circadian physiol-ogy in a state-space model, and to integrate multiple measurements and model predictions in a coherent manner. A n estimation algorithm is now required to generate circadian state es-timates within this framework. There are a number of options for estimation algorithms, so we begin by stepping through a selection process that leads to the choice of a particle filter as the appropriate choice. The detailed design of the particle is then presented along with a dis-cussion of the tuning parameters. Final ly the implementation of a kernel probability density construction algorithm, which is a useful complement to the particle filter, is detailed. 5.1 Estimation Algorithm Selection Methods for calculating an underlying property of a system based on observed measurements and a system model include multiple variants of recursive filtering and parametric curve fitting. The general circadian estimation problem involves a dynamic system model with process noise (rather than an algebraic system model), incremental processing of measure-ments (rather than batch data analysis), and statistical distributions of measurements and parameter estimates. In response the focus of estimation algorithm selection is narrowed to recursive filtering methods using Bayesian statistics, which will be referred to as 'Bayesian 37 Chapter 5. Particle Filter Phase Estimator e s t i m a t i o n ' f r o m he re . T h e d e f i n i n g f ea tu re o f B a y e s i a n e s t i m a t i o n i s t h e i n t e g r a t i o n o f p re -d i c t i o n u p d a t e s a n d m e a s u r e m e n t s u p d a t e s u s i n g i n f e r e n t i a l s ta t i s t i c s . A n o v e r v i e w o f r e c u r s i v e B a y e s i a n e s t i m a t i o n i s p r e s e n t e d b e l o w w i t h t h e g e n e r a l m a t h e -m a t i c a l e q u a t i o n s as t h e y a re a p p l i e d to s ta te -space m o d e l s . T h e r e q u i r e m e n t s fo r d e v e l o p i n g a s o l u t i o n to t he B a y e s i a n e q u a t i o n s spec i f i ca l l y for t he c i r c a d i a n e s t i m a t i o n fo l low, a n d t h e n a c o m p a r i s o n o f d i f fe ren t s o l v i n g a l g o r i t h m s i s c o n d u c t e d w h i c h l e a d s to t he se lec t i on o f t h e p a r t i c l e filter. 5.1.1 Recursive Bayesian Estimation Overview B a y e s i a n s ta t i s t i c s p rov i des a n a p p r o a c h to on - l i ne e s t i m a t i o n i n w h i c h t h e p robab i l i t y , o r be l ie f , o f a sys tem 's p r o p e r t y i s u p d a t e d b a s e d on a n i n i t i a l bel ie f , n e w m e a s u r e m e n t s , a n d p r e d i c t i o n s o f t h e sys tem 's i n t e r n a l d y n a m i c s . A s t a t i s t i c a l p r o b a b i l i t y i s a s s o c i a t e d w i t h e a c h o p e r a t i o n w h i c h a l l o w s for n a t u r a l e x p r e s s i o n o f t h e u n c e r t a i n t y t h a t i s i n h e r e n t i n m o s t r e a l s y s t e m s , a n d a l l o w s for o p t i m a l a n a l y s i s o f n o i s y m e a s u r e m e n t s . A p p l i e d to s ta te -space m o d e l s , t he v a l u e s o f a s ta te vec to r a r e e s t i m a t e d s t a r t i n g f r o m a p r i o r p r o b a b i l i t y d i s t r i b u t i o n a n d s e q u e n t i a l l y u p d a t e d w i t h p r e d i c t i o n s f r o m a s ta te t r a n s i t i o n e q u a t i o n a n d a d j u s t m e n t s f r o m a m e a s u r e m e n t e q u a t i o n . L e t ' s c o n s i d e r t h e d i s c re te - t ime s ta te -space m o d e l w i t h a s ta te vec to r x t h a t i s e v a l u a t e d a t t i m e s ti = iT, w h e r e T i s t h e s a m p l i n g p e r i o d . D e n o t i n g m e a s u r e m e n t s a t t i m e tt as zt, t he set o f a l l obse rva t i ons u p to t i m e k i s d e f i n e d as Zk = {zi} i = 1 , . . . , k}. T h e s ta te -space m o d e l cons i s t s o f a s ta te t r a n s i t i o n m o d e l a n d a m e a s u r e m e n t m o d e l w h i c h a re e a c h u s e d i n B a y e s i a n u p d a t e e q u a t i o n s to u p d a t e the s ta te e s t i m a t e s . T h e s ta te t r a n s i t i o n m o d e l x f c = f f c ( x f c _ i , U f c _ i , V f e _ i ) (5.1.1) desc r i bes t h e e v o l u t i o n of s ta tes f r o m a p r i o r s ta te x f c _ i to a f u t u r e s ta te x f c w i t h m e a s u r e d 38 Chapter 5. Particle Fi lter Phase Estimator i n p u t s Uk a n d no ise Vk- T h e m e a s u r e m e n t m o d e l z f c = h f c ( x f e , w f c ) (5.1.2) desc r i bes a r e l a t i o n s h i p b e t w e e n a n obse rved m e a s u r e m e n t z f c a n d t h e sys tem 's i n t e r n a l s ta te Xfc a n d m e a s u r e m e n t no i se wk-T o a p p l y B a y e s i a n s t a t i s t i c a l log ic , b o t h t h e s ta te t r a n s i t i o n a n d m e a s u r e m e n t e q u a t i o n s a r e e x t e n d e d to i n c l u d e p r o b a b i l i t y d i s t r i b u t i o n s . A s s u m i n g a k n o w n p r o b a b i l i t y d i s t r i b u t i o n o f t h e s ta tes a t t i m e t k - \ i s p(xfc_i |Zfc_i) t h e n t h e p r o b a b i l i t y d i s t r i b u t i o n o f s ta tes a t a f u t u r e t i m e s tep p(xfc|Zfc_i) i s d e f i n e d b y t h e C h a p m a n - K o l m o g o r o v e q u a t i o n [48] /oo p(xfc |x fc_ 1 )p(x f e _i |Zfc_ 1 ) ( ix fc_ 1 . (5.1.3) -oo T h i s prediction update e q u a t i o n a l l ows t h e b e l i e f i n t he s y s t e m s ta te vec to r to evo lve i n t i m e b a s e d o n a p r e d i c t i o n s f r o m t h e s ta te t r a n s i t i o n m o d e l . A p r o b a b i l i s t i c u p d a t e o f s ta tes b a s e d i n m e a s u r e m e n t s i s p r o v i d e d by B a y e s t h e o r e m w h i c h s ta tes t h a t p(zfc|x f c ) -p(xfc|Zfc_i) p(xfc|z f c ) = — — (5.1.4) p ( z f c l z f c - i ) w h e r e p(xfc|Zfc_i) r e p r e s e n t s t h e p r i o r k n o w l e d g e of t he s ta tes , p(zfc|xfc) i s t h e l i k e l i h o o d f r o m t h e o b s e r v e d d a t a , a n d p(xfc |zfc) i s the pos te r i o r p r o b a b i l i t y o f t h e s ta te v a r i a b l e s . T h e d e n o m -i n a t o r p(zfc|Zfc_i) i s a n o r m a l i z a t i o n c o n s t a n t , the re fo re E q . 5.1.4 c a n be e x p r e s s e d as p(xfc|zfc) = Cp(zfc|xfc)p(xfc|Zfc_i). (5.1.5) w h e r e C i s a n o r m a l i z a t i o n cons tan t . T h i s measurement update e q u a t i o n a l l o w s a n u p d a t e of t h e sys tem 's s ta te vec to r b a s e d on n e w m e a s u r e m e n t s . G i v e n a n i n i t i a l i z a t i o n of t h e p r i o r p r o b a b i l i t y d i s t r i b u t i o n s o f t h e s ta tes , t h e m e a s u r e -m e n t u p d a t e a n d p r e d i c t i o n u p d a t e s teps a re a p p l i e d i t e r a t i v e l y as s h o w n i n F i g . 5 .1 . T h i s g e n e r a l sequence o f ope ra t i ons cons t i t u tes B a y e s i a n filtering. T h e spec i f i c i m p l e m e n t a t i o n o f t h e m e a s u r e m e n t u p d a t e a n d p r e d i c t i o n u p d a t e s teps v a r y w i d e l y h o w e v e r b a s e d o n t h e c h a r a c t e r i s t i c s o f t he s y s t e m o f i n t e r e s t a n d o n d e s i g n c r i t e r i a . 39 Chapter 5. Particle Fi lter Phase Estimator INITIALIZE PREDICTIONTJPDATE M E A S U R E M E N T U P D A T E Figure 5.1: Recursive Bayesian estimation steps. 5.1.2 Requirements for Circadian State-Space Model Estimation To implement Bayesian filtering for circadian phase estimation, a specific solving algorithm must be developed. Integrating the state space circadian pacemaker model and the physio-logical marker components into the general recursive Bayesian filtering scheme from leads to the specific system architecture illustrated in Fig. 5.2. An appropriate selection of algorithms for the prediction and measurement update blocks is then based on the properties of the cir-cadian model. Important characteristics include the degree and types of nonlinearity, and the properties of the noise. Analysis of the circadian state-space model developed in Sec. 4.1.3 highlighted three important features: 1. the phase state, <f>, has a nonlinear parameter space; 2. the state transition equation is nonlinear in a way that may lead to bimodal probability densities; 3. the unknown variability of light input should be treated parametrically rather than a simple additive, independent Gaussian random variable. Each of these points is discussed in detail below. 40 Chapter 5. Particle Filter Phase Estimator Phys io log ica l Phase Es t ima te Measurement U p d a t e State Estimate Pred i c t i on U p d a t e Figure 5.2: Block diagram of the circadian phase estimator. A typical assumption for many systems is that the state vector exists in independent linear parameter spaces of continuous real number such that x e R™. In examining the state vector of the transformed circadian model Eq. 4.1.9, A n we can see that the amplitude A and light response n have typical parameter spaces however the phase offset <f> is an exception. By its definition, <f> indicates the phase point within a twenty-four hour day, which like the hour hand of a clock is constrained to the range 0 < 4> < 24 and is in a circular parameter space since the time of OhOO is equivalent to 24h00. The parameter space of the phase state variable can be thus defined using the modulo operator as 4> e (R mod 24). A graphical representation of the parameter spaces for each variable in the state vector is illustrated in Fig. 5.3. 41 Chapter 5. Particle Fi l ter Phase Estimator A <p n F i g u r e 5.3: R e p r e s e n t a t i o n o f t h e p a r a m e t e r spaces o f t he s ta te v a r i a b l e s f r o m E q . 4 .1 .9 , s h o w i n g the p h a s e s ta te b e l o n g i n g to a c i r c u l a r r a n g e <j> e (R m o d 24), a n d the a m p l i t u d e a n d l i g h t filter s ta tes w i t h t y p i c a l d o m a i n s o f A e R a n d n € R. T h e second i m p o r t a n t f ea tu re o f t he c i r c a d i a n p a c e m a k e r m o d e l i s t h a t t h e s ta te t r a n s i -t i o n e q u a t i o n s a re n o n l i n e a r . T h e V a n de P o l osc i l l a t o r e q u a t i o n s , E q . 3.2.1 a n d E q . 3 .2 .2 , c o n t a i n h i g h o r d e r t e r m s o f x a n d xc, a n d the l i g h t i n p u t e q u a t i o n s , E q . 3.2.3 a n d E q . 3 .2.5, c o n t a i n e x p o n e n t i a l a n d m u l t i p l i c a t i v e t e r m s . P r e v i o u s w o r k s h o w e d t h a t l i n e a r a p p r o x i m a -t i ons c o u l d be m a d e for t h e m o d e l [47]. T h e l i n e a r i z e d m o d e l w a s u s e d to e n a b l e a m o d e l -p r e d i c t i v e con t ro l a l g o r i t h m to d e t e r m i n e l i g h t i n p u t s for c h a n g i n g t h e c i r c a d i a n p h a s e . I n the c u r r e n t con tex t however , o u r spec i f i c ob jec t ive i s h i g h a c c u r a c y t r a c k i n g o f t h e c i r c a d i a n p h a s e so l i n e a r i z a t i o n w a s no t c o n s i d e r e d a good s o l u t i o n u n l e s s i t w a s s h o w n t h a t a l l t h e n o n l i n e a r effects w e r e m a r g i n a l . O n e n o n l i n e a r p r o p e r t y of t he c i r c a d i a n s y s t e m t h a t w a s w e l l k n o w n q u a l i t a t i v e l y i n t h e l i t e r a t u r e w a s t h a t t he t i m i n g o f l i g h t a r o u n d t h e p h a s e m i n -i m a c o u l d l e a d to a d i v e r g e n c e o f p h a s e sh i f t d i r ec t i ons . I f l i g h t i s i n t r o d u c e d s l i g h t l y before t h e p h a s e m i n i m a i t w i l l c a u s e a sh i f t b a c k w a r d , a n d i f s l i g h t l y a f te r t h e n i t w i l l c a u s e a sh i f t f o r w a r d . I t w a s i n f e r r e d t h a t t h i s d i ve rgence b e h a v i o u r c o u l d l e a d to b i - m o d a l p r o b a b i l i t y d i s t r i b u t i o n s . T h e final n o n l i n e a r i t y w h i c h w a s obse rved w a s t h a t v a r i a b i l i t y i n l i g h t i n p u t w o u l d no t be m o d e l l e d w e l l b y a s i m p l e a d d i t i v e G a u s s i a n r a n d o m v a r i a b l e i n t he f o r m o f I = I + v. T w o p r i m a r y d e t e r m i n i n g fac to rs o f i n d i v i d u a l l i g h t e x p o s u r e a r e t h e t i m i n g o f l i g h t expo-su re c h a n g e s f r o m s leep a n d w a k e t r a n s i t i o n s a n d the a m b i e n t l i g h t l e ve l s r e l a t e d to l o c a t i o n (e.g. s u n l i g h t , d i m r o o m , b r i g h t room) . F o r t he p u r p o s e s o f t h i s pro ject a n a s s u m p t i o n w a s 42 Chapter 5. Particle Filter Phase Estimator m a d e t h a t l i g h t i n p u t c a n be d e t e r m i n e d as a f u n c t i o n o f a t i m i n g p a r a m e t e r , a n d l e v e l p a -rame te r . I t t h e n fo l l ows t h a t u n c e r t a i n t y s h o u l d be i n t r o d u c e d as r a n d o m v a r i a b i l i t y to t h e p a r a m e t e r v a l u e s . A M o n t e C a r l o s i m u l a t i o n i l l u s t r a t i n g t h e d i f fe rence b e t w e e n l i g h t v a r i -a b i l i t y i n t r o d u c e d as a a d d i t i v e r a n d o m v a r i a b l e a n d v a r i a b i l i t y i n t r o d u c e d to t i m i n g a n d l e v e l p a r a m e t e r s i s s h o w n i n F i g . 5.4. W h i l e no t b a s e d o n q u a n t i t a t i v e s t ud ies , i t w a s as -400 300 ~ 200 _1 100 0 0 8 16 24 32 40 48 56 64 72 80 Time [hours] 400 300 £ 200 CT) _J 100 0 0 8 16 24 32 40 48 ' 56 64 72 80 Time [hours] F i g u r e 5.4: I l l u s t r a t i o n o f t w o m e t h o d s o f r a n d o m i z i n g l i g h t i n p u t d u r i n g w a k i n g pe r i ods u s i n g M o n t e C a r l o s i m u l a t i o n w i t h t e n p a r t i c l e s : a n a d d i t i v e G a u s s i a n r a n d o m v a r i a b l e (top) a n d r a n d o m v a r i a b l e s i n t r o d u c e d for t r a n s i t i o n t i m i n g a n d l e v e l p a r a m e t e r s (bot tom). s u m e d t h a t a p a r a m e t e r i z e d l i g h t v a r i a b i l i t y m e t h o d s w o u l d be m o r e r e p r e s e n t a t i v e o f a t yp -i c a l s cena r i os t h a n t h e s i m p l e a d d i t i v e no i se as p a r a m e t e r s w o u l d m o r e c lose ly m o d e l h u m a n b e h a v i o u r a l c h a r a c t e r i s t i c s . B a s e d on t h i s a n a l y s i s t h e spec i f i c p a r a m e t e r i z a t i o n w a s no t i m -p o r t a n t , h o w e v e r i t w a s d e t e r m i n e d t h a t t h e B a y e s i a n filtering a l g o r i t h m w o u l d r e q u i r e t h e g e n e r a l c a p a c i t y to m o d e l n o n l i n e a r l i g h t i n p u t v a r i a b i l i t y . 43 Chapter 5. Particle Filter Phase Estimator 5.1.3 Particle Filter Algorithm Selection Implementing a recursive Bayesian filter for a given system requires developing solutions to equations 5.1.5 and 5.1.3. For linear systems with Gaussian noise, analytical solutions exist resulting in the classical Kalman filter, however in cases with nonlinearities, approximate solutions must usually be used. There are variety of well-developed general algorithms, so implementation for the circadian model first involved selecting an appropriate class of solu-tions to meet the requirements identified in Sec. 5.1.2. Three common classes of Bayesian filter algorithms are the Kalman Filter, the Extended [21 and Unscented Kalman Filters [54], and Particle Filters [3]. The common feature of all Kalman filters is that they represent probabilities with Gaussian distributions. Each proba-bility is completely specified by two parameters: the mean and the standard deviation. Parti-cle filters however represent probability distributions with a set of point masses (see Fig. 5.5) and use sequential Monte Carlo estimation. A particle representation is an approximation 0 5 10 0 5 10 x x Figure 5.5: Comparison of a Gaussian probability density function specified by its mean and standard deviation (left) and a representation of the same probability distribution using point masses (right). so the accuracy will increase with the number of particles. Particle filters can be compu-tationally intensive and are complex to tune, however their strength lies in their ability to represent probability densities of arbitrary shapes. Selection between the methods is dependent on a number of factors including whether 44 Chapter 5. Particle Filter Phase Estimator the model is linear or nonlinear, whether the noise is Gaussian or not, whether the system needs to have multi-modal PDFs to accurately represent its behaviour, and any computa-tional constraints. Table 5.1 summarizes the capabilities of the three filter classes against these criteria. The Kalman filter can only handle linear, Gaussian noise systems, while the Table 5.1: Filter Selection Comparison Kalman Filter Extended / Un- Particle Filter scented Kalman Filter Linear Model V V V Gaussian Noise V V V Nonlinear Model x y/ y/ NonGaussian Noise x V V Multimodal PDFs x x yj Computation Time 0(n) 0(n) 0{m),m^>n n =number of states, in =number of particles Extended and Unscented Kalman Filters extend their capabilities to some nonlinear systems, and the particle filter is the most flexible method. There is a cost in terms of complexity and computation time for the more flexible algorithms therefore the simplest method that can meet the estimation requirements is preferable. The requirements for the circadian system, identified in the previous section, were the capability to accurately track nonlinearities in the state parameter space, a state transition function that could lead to bimodal probability densities, and a nonlinear input noise model. Given the high number of nonlinear elements the particle filter was selected as the Bayesian filter implementation method. 5.2 Particle Filter Design There are many variations to particle filtering, and the method selected here is Sequential Importance Resampling with a Markov chain Monte Carlo move step. The algorithm is dis-cussed in detail below, followed by a discussion of the noise models and tuning parameters. 45 Chapter 5. Particle Fi lter Phase Estimator 5.2.1 Particle Filtering Algorithm Particle filters consist of a sequence of operations that propagate a set of particles through recursive Bayesian filters. As with all Bayesian filtering algorithms a prediction update and measurement update are necessary operations (see Fig. 5.1). In particle filters additional steps such as resampling are required to ensure stability and optimize performance [15; 16]. There are a large number of proposed particle filter variants which can be chosen to meet the needs of a particular problem. For the circadian tracking application, a four step particle filter was implemented consisting of Sequential Importance Resampling (SIR) [20] and a Markov Chain Monte Carlo (MCMC) Move step [19]. The full sequence of steps therefore consists of a prediction update to project each par-ticle forward in time based on predictions from the dynamic circadian pacemaker model, a measurement update to weight the particles based on a likelihood determined from phase measurements, then resampling to remove particles in areas with low probabilities and in-crease the particles in areas with high probability, and finally the move step redistributes the particles to avoid convergence problems. This sequence of operations is shown schematically in Fig. 5.6 and listed in Algorithm 1. Algorithm 1 Particle Filter l : procedure [xfc] = P A R T I C L E F l L T E R ( x f c _ i , u^_i.z*.) 2: for i = 1 : N do 3: x\ = l N C R E M E N T ( x ^ _ j . U f c _ j , Zfc, Qfc, Rfc) > x\ ~ p ( x k . **) 4: wl = l M P O R T A N C E W E I G H T ( x j i , z f e , Q f c , R f c ) t> wi oc pCzj^xj.) 5: xlk = R E S A M P L E ( x ^ ) t> returns evenly weighted particles 6: xi = M C M C - M O V E ( x ^ u f c _ 1 , Q f c , R f c ) 7: end for 8: end procedure 5.2.1.1 Increment Operation Beginning with the particle distribution represented in the state's prior PDF at time k - 1 which is p(xj._j |Zfc_i), the I N C R E M E N T step propagates each particle based on the system's 46 Chapter 5. Particle Fi lter Phase Estimator 9 I N C R E M E N T 0 0 I M P O R T A N C E W E I G H T i 4 4 i S /'(x:!x: ;.z. ;.) 0 M C M C - M O V E 6 p ( x i i z j Figure 5.6: Diagram showing progression of particles through one cycle of the particle filter operations. state transition model. The algorithm is listed in Algorithm 2. Algorithm 2 Prediction Update l: procedure [x'fc] = i N C R E M E N T f x j ^ . U f c - i . Z f c . Q j ^ R f c ) 2: H = V h ( x f e l ) 3: S = H Q f c H ' + R f c 4: S = Q f c - Q f c H ' S ^ H Q f c 5: for i = 1 : N do 6: xj, — f (x]i._ j , U fc_ i ) > prediction update 7: 4-1 =h(x l f e ,u f e _ 1 ) 8: x\ = x\ + £ H ' R ^ 1 M O D M l N U S ( z f c . z ^ , 24) i> measurement update 9: Draw xj: ~ W(x^.. d i a g E ) > draw x!fc from p(xk\xlk_1,zk) 10: x' f c(2) = mod (x [ ( 2 ) , 24) > keep phase within [0,24) 11: end for 12: end procedure 5.2.1.2 Importance Weight Operation Next the I M P O R T A N C E W E I G H T operation weights each particle based on the likelihood func-tion p(zk\xlk). The weighted particles now represent the posterior probability density p(xk\Zk). The algorithm is listed in Algorithm 3 47 Chapter 5. Particle Filter Phase Estimator Algorithm 3 Measurement Update 1: procedure [wjp = IMPORTANCE W E I G H T ^ . , zk. Qk. Rk) 2: S = H Q f e H ' + R f c 3: for i = 1 : N do 4: ±i = h(xik:uk) 5: = MODMlNUS(zfc. i\, 24) o phase difference using modulo 6: wj. = Af(dj,-,0. S(2.2)) > set particle weight to p(zfc|x'A.) 7: end for 8: for i = \ : N Ao 9: wjj. = w ^ / E ^ W j > normalize the weights 10: end for 11: end procedure 5.2.1.3 Resampling The R E S A M P L E operation then discards particles in areas of low probability and multiples the particles in area of high probability. While the particles now represent the desired prior density function a practical problem is that the resampling step reduces the diversity of particle locations and the particles would eventually collapse to a single point. The algorithm is listed in Algorithm 4. Algorithm 4 Systematic Resampling l 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 procedure [x7 f c,wjp = RESAMPLE(xj.,wj.) 3 = 1 C = N~1 CSW = 0 while C •< 1 do if CSW>C then C = C + N-1 else select i randomly from [j, N] CSW = CSW + w[ x l = xi 3=3 + 1 end if end while for i = 1 : N do w i = A T - 1 end for end procedure > increment cumulative sum of weights i> swap new particle > set weights to uniform values 48 Chapter 5. Particle Fi l ter Phase Estimator 5.2.1.4 Markov chain Monte Carlo Move To redistribute the particles while still maintaining the statistical distribution characteristics the particles are moved using a Markov chain Monte Carlo in the M C M C - M O V E step. The M C M C is implemented using the Metropolis-Hastings algorithm. The algorithm for move operation is listed in Algorithm 5. Algorithm 5 Markov chain Monte Carlo Move l : procedure [xjp = M C M C - M O V E ( x J . . uk, Q f c , R f e ) 2: S = H Q f c H ' + Rfe 3: for i = 1 : N do 4: Z = h(x f c ,U f c ) 5: dj, = M O D M l N U S ( z f c , z{. 24) 6: 1/ = Af(d^: 0, S) > likelihood of current particle p(zfc|xfc) 7: draw xjj. ~ U({(xk_1,uk), Qfc) > generate a new random proposal 8: z = h(xji,Ufc) 9: = M O D M l N U S ( z f c , z{, 24) 10: L z = N(dj,\ 0, S) > likelihood of proposed particle p(zfc|xfc) 11: PA = min(l, \J1 t> probability of acceptance of new particle 12: draw r ~ [0,1] o random number from 0 to 1 13: i f Pa > r then 14: x f c = xlk > replace x with new proposal 15: end if 16: end for 17: end procedure 49 Chapter 5. Particle Filter Phase Estimator 5.3 Probability Density Reconstruction While particle locations sufficiently represent probability distributions for the use in the esti-mation algorithm, they are difficult to graphically interpret. Continuous probability density functions are much better suited as communication tools, therefore a method for converting from particle distributions to continuous probability density functions is desired. Approaches for generating probability density functions from point samples fall into the parametric and nonparametric classes. Parametric methods include maximum likelihood fitting of a known distribution functions, while non-parametric methods do not assume specific distribution shape and include histograms and kernel density functions [29] approaches. Since we want to allow arbitrarily shaped distributions for circadian phase, parametric methods are not an appropriate choice. Among the non-parametric methods histograms provide a simple method however they are sensitive to placement of bin locations. The kernel density approach while more computationally demanding offers superior consistency and smoothness therefore is selected for use with the particle filter. The concept of the kernel density method is to replace each particle with a continuous 'Kernel' function. Given a set of particles then given a kernel function w(x), an approximation of the continuous probability density function is given by N p(x) = ^2w{xi) (5.3.2) A Gaussain PDF is a often selected as the kernel function due to its well understood statis-tical properties and smooth window. An example of a Gaussian kernel replacement is shown in Fig. 5.7. To achieve accurate probability density reconstructions using kernels two inter-related parameters must be selected. One is the shape of the kernel function, and another N (5.3.1) 50 Chapter 5. Particle Filter Phase Estimator F i g u r e 5.7: A p r o b a b i l i t y d e n s i t y f u n c t i o n r e p r e s e n t e d b y d i sc re te p o i n t m a s s e s (a) i s r econ -s t r u c t e d b y r e p l a c i n g e a c h po in t w i t h a G a u s s i a n k e r n e l (b) a n d t h e n s u m m i n g a l l t h e k e r n e l s to c rea te a s i ng l e p r o b a b i l i t y d i s t r i b u t i o n (c). i s t h e n u m b e r o f pa r t i c l es . T h e p r i m a r y s h a p e p a r a m e t e r o f t h e k e r n e l i s t he b a n d w i d t h , o r t h e w i d t h . F o r G a u s s i a n k e r n e l s t h e s t a n d a r d d e v i a t i o n i s c o n s i d e r e d t h e c o r r e s p o n d i n g b a n d w i d t h pa rame te r . 5.3.1 Probability Distribution Comparison Metrics T o a s s i s t w i t h the se lec t i on o f a p p r o p r i a t e b a n d w i d t h p a r a m e t e r s , tes t cases c a n be e v a l u -a t e d i n w h i c h t h e u n d e r l y i n g t r u e p r o b a b i l i t y d i s t r i b u t i o n i s k n o w n . T h e a c c u r a c y o f a r e c o n -s t r u c t e d p r o b a b i l i t y d e n s i t y f u n c t i o n , p(x), c a n be d e t e r m i n e d by c a l c u l a t i n g t h e s i m i l a r i t y to t h e k n o w n d i s t r i b u t i o n p(x). T w o m e t h o d s a re the K u l l b a c k - L i e b l e r d i ve rgence [39] g i v e n b y f°° v(x) DKL(P,P)= / p(x)log^-Ldx (5.3.3) . / - o o P\X) a n d t h e m e a n i n t e g r a t e d s q u a r e d e r r o r g i v e n by /oo E[{p{x) - p{x)}2)dx. (5.3.4) -oo A s a n e x a m p l e le t u s t a k e a b i m o d a l d i s t r i b u t i o n s h o w n i n F i g . 5.8 as a ' t rue ' d i s t r i b u t i o n a n d t h e n tes t r e p r e s e n t a t i o n s o f i t u s i n g v a r i o u s c o m b i n a t i o n s o f p a r t i c l e s n u m b e r s a n d ke r -n e l b a n d w i d t h s . F o r e a c h r e c o n s t r u c t i o n t h e K u l l b a c k - L i e b l e r d i v e r g e n c e i s c a l c u l a t e d a n d a su r f ace m a p i n F i g . 5.9 i s g e n e r a t e d s h o w i n g the o p t i m a l r e c o n s t r u c t i o n zone i n w h i t e . 51 Chapter 5. Particle Filter Phase Estimator 0.14 x Figure 5.8: Bimodal probability distribution used as an example to test particle representa-tion. 5.3.2 Kernel Density Estimation In situations where the true probability density is unknown, which occurs in most cases, the challenge of selecting the optimal number particles and bandwidths is not straightforward. Based on the observations from the previous test case we can start with some general prin-ciples. As can be seen from Fig. 5.9, the accuracy always improves with a greater number of particles. So more particles are better. A practical constraint in that dimension however is computational power. The selection of bandwidth does not follow a simple rule however since it case be seen in Fig. 5.9 that high divergence can occur with both too large a bandwidth and too small a bandwidth. To address this challenge a class of optimization solutions called Kernel Density Esti-mators have been developed for automatically selected bandwidths based on minimization of the mean integrated squared error (MISE) for a given set of particles. A Kernel Density Estimation algorithm based on a dual-tree algorithm [21] and selected from a Matlab K D E toolbox [27] is used here to implement K D E . To illustrate its performance we can repeat the bimodal fitting test. For each set of N particles drawn from bimodal distribution, rather than testing an array of different bandwidths, a single 'optimal' bandwidth can be selected using the K D E algorithm. The bandwidths selected by the K D E are shown in Fig. 5.10 overlayed 52 Chapter 5. Particle Fi lter Phase Estimator 600 500 400 300 200 100 KL Divergence A) B) C) M PA AA E) F) M M 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 X X X Kernel Bandwidth Figure 5.9: Surface map of Kullback-Liebler divergence for various kernel bandwidths and numbers of particles (left) where lowest values, shown as white, correspond to the best matches. Examples of reconstructed distributions (right), corresponding to the points from the Kullback-Liebler surface map, illustrate the different types of probability distributions reconstructed, with the best fit to the original bimodal distribution occurring at point B. on the Kullback-Liebler divergence surface. Despite the fact that the K D E does not use any knowledge of the true P D F , it can be seen that it selects bandwidths which minimize the divergence. This K D E method is used for all subsequent translations from particle represen-tations to continuous probability density function. 53 Chapter 5. Particle Filter Phase Estimator KL Divergence 0.5 1 1.5 2 Bandwidth Figure 5.10: The bandwidth selected by the Kernel Density Estimation algorithm for a given number of particles using the bimodal probability distribution, is shown by the black line superimposed on the K L Divergence surface map. Note that the K D E algorithm has no knowledge of the true probability distribution but selects near-optimal bandwidth choices. 54 Chapter 6 S imu la t i on Resu l ts The particle filter approach to circadian phase estimation, in conjunction with the new mod-elling framework, introduce the capability for real-time tracking of a circadian phase with Bayesian probability distributions. To explore the theoretical capabilities and limitations of this approach, and to determine appropriate tuning parameters for the particle filter, simu-lations are performed for series of scenarios with different particle filter parameter choices, light input assumptions, and noise models. 6.1 Scenario 1 A preliminary simulation scenario is designed to explore the tuning parameter selection of the number of particles and the level of process noise. We will take a typical scenario and test the performance of the particle filter with various combinations of the tuning parameters. As a typical scenario a light input setting is set to a pattern representative of an individual consistently following a standard sleep schedule for seven days. Lights are turned off at midnight for an eight hour sleep episode and then turned on at eight a.m. For this simulation we will also assume a consistent light intensity of600 Lux throughout the day. The light/dark exposure pattern for is shown in Fig. 6.1. A good test the behaviour of the particle filter is to use a uniformly distributed initial particle distribution and observe the convergence after multiple simulation days. For this 55 Chapter 6. Simulation Results Table 6.1: Simulation 1 Parameter Value - Case with No Process Noise Parameter Value No. Particles 12 24 48 72 240 480 960 Initial Distribution Uniform over [0,24] Process Noise, Qg (negligible) Light Levels Dark 0 Lux Bright 600 Lux Measurements (none) simulation it means we are considering the case in which no a priori information is known about the individual's circadian phase, and would like to know the final probability distribu-tion to indicate a degree of belief in the location of the circadian phase after one week of a consistent schedule. What is expected is that an entrainment effect will occur and the phase distribution will converge with a mean of 0 = 4. To map the effects of the number of particles the simulation is run with seven different choices of particle numbers, and in a first case, is run with no process noise. The settings listed in Table 6.1 are used. Another set of simulations is run, repeating the choice of particle numbers, but this time with the phase state process noise set to a non-negligible value of Q = 3 x 10~ 4 as listed in Table 6.2. Since we currently have no quantitative data from which a process noise level can be selected, this value was chosen based on qualitative physical assumptions about variability inherent in the circadian pacemaker system. After running Clock Time • Light Figure 6.1: Light schedule for simulation 1, following a regular eight hour sleep schedule. Dark periods are 0 lux and light periods are 600 lux. 56 Chapter 6. Simulation Results Table 6.2: Simulation 1 Parameter Value - Case With Process Noise Parameter Value No. Particles 12 24 48 72 240 480 960 Initial Distribution Uniform over [0,24] Process Noise, Qq 3 X 10~ 4 Light Levels Dark 0 L u x Bright 600 L u x Measurements (none) the particle filter with this array of parameters we have results from 14 trials to examine. Illustrating the results from a single trial, the particle trajectories for each of the three states in the simulation with 24 particles and no process noise is shown in Fig. 6.2. With only twenty-four particles and negligible process noise the precise path of each particle. The phase of the particles become further entrained to the timing of the sleep/wake schedule after each light exposure period. The most useful comparative analysis is an examination of the final Amplitude State, A 1.5i 1 1 1 r— 0 1 2 3 4 5 6 7 Phase State, <(> Time {hours) Figure 6.2: Particle trajectories for the three state variables during Simulation 1A. probability distributions predicted by each algorithm after the end of day seven. The selection 57 Chapter 6. Simulation Results of an appropriate number of particles can be chosen based on a point at which increasing numbers of particles still converge at the same distribution. Referring back to the Kernel Density Estimate Figure 5.10, it is the converging behaviour that indicates a good match to the 'true' probability. The left panel of Fig 6.1 shows the predictions from each number of particles in the case with no process noise and the right panel shows the predictions in case with process noise on the phase state. In all simulations a peak in the distribution consistently converges at approximately 4h, which is the expected phase for entrainment to a sleep schedule between 0 and 8 hours. When there is no process noise the sharpness of the probability distribution increases as the number of particles increases from 24 to 960 without an apparent conver-gence, indicating significantly larger numbers of particles would be required to accurately capture the true PDF. On the right panel of Fig 6.1, a convergence can be seen with the predictions from 240, 480, and 960 particles all producing similar results. An explanation a 0.15 0.05 0 2 4 6 8 10 12 14 16 18 20 22 24 Circadian Fhase, *> (hours) 0.25 0.15 0.05 -A— 12,BW=2.8663 -e—24,BW=1.6307 -6)— 48,BW=0S7551 -0— 72,BW=1.0787 -#— 240,BW=0.52047 —I— 480,BW=0.42443 960, BW=038727 2 4 6 8 10 12 14 16 18 20 22 24 Circadiai Fhase, $ (hours) Figure 6.3: Comparison of multiple simulations with different numbers of particles in case of no process noise (left) and the case with process noise (right). of this behaviour is that the process noise is limiting the minimum width of the probability distribution and introducing a tendency towards a Gaussian shape. 58 Chapter 6. Simulation Results W e a s s u m e t h a t t he i n t r o d u c i n g a n o n - n e g l i g i b l e p rocess no ise i s a r e a l i s t i c a s s u m p t i o n a l t h o u g h a spec i f i c v a l u e c a n n o t be j u s t i f i e d a t t h i s po in t . T o p roceed w i t h f u r t h e r a n a l y s i s t h e cho ice o f Q = 3 x 1 0 ~ 4 w i l l be se lec ted as t h e de fau l t l e v e l , a n d the re fo re a cho ice o f N > = 240 s h o u l d be se lec ted to ach ieve convergence to a n a c c u r a t e P D F . I n f u r t h e r s i m u l a t i o n s a d e f a u l t v a l u e of N = 240 w i l l be u s e d to m i n i m i z e c o m p u t a t i o n a l comp lex i t y . T h e o p t i m a l n u m b e r o f p a r t i c l e s m a y v a r y for d i f fe ren t s i m u l a t i o n scena r i os , a n d c h o o s i n g a h i g h n u m b e r o f p a r t i c l e s c a n p r o v i d e a bu f fe r for those cases w i t h no p e n a l t y excep t for c o m p u t a t i o n a l comp lex i t y . A s p r o c e s s i n g t i m e i s s i g n i f i c a n t (20 m i n u t e s to r u n t h i s s i m u l a t i o n w i t h N = 9 6 0 ) t h e cho ice o f N = 240 w i l l be u s e d as a d e f a u l t for s u b s e q u e n t s i m u l a t i o n s . 6.2 Scenario 2 T h i s s c e n a r i o f u r t h e r e x a m i n e s t he p r e d i c t i o n s of e n t r a i n m e n t to a fixed s c h e d u l e a n d tes t s d i f f e ren t no i se mode l s . W e choose to s i m u l a t e t he case o f a n i n d i v i d u a l f o l l o w i n g a cons i s -t e n t l y t i m e d e igh t h o u r s leep r e g i m e n for a d u r a t i o n o f t h r e e w e e k s as s h o w n i n F i g . 6.4. T h e w a k i n g l i g h t l eve l s a re chosen to be t y p i c a l of i n d o o r b r i g h t l i g h t of 3 8 0 L u x . C o n s i d e r i n g Clock Time • Light F i g u r e 6.4: L i g h t schedu le for s i m u l a t i o n 2. a case i n w h i c h no a p r i o r i i n f o r m a t i o n i s k n o w n abou t t h e i n d i v i d u a l ' s c i r c a d i a n p h a s e , t h e final p r o b a b i l i t y d i s t r i b u t i o n w i l l i n d i c a t e a degree o f b e l i e f i n t h e l o c a t i o n o f t he c i r c a d i a n p h a s e a f te r t h e t h r e e w e e k s o n t h i s schedu le . T h i s a n a l y s i s h a s a s i g n i f i c a n t c o n n e c t i o n to a c l i n i c a l s cena r i o s ince a t h r e e - w e e k m o n i t o r i n g p e r i o d i s t y p i c a l l y u s e d to e n s u r e t h a t sub jec ts 59 Chapter 6. Simulation Results Table 6.3: Simulation 2 Parameter Values Parameter Value Case 2A Case 2B Case 2C No. Particles 240 240 240 Initial Distribution Uniform over [0,24] Process Noise, Qg 1 x 10- 9 3 x 10~4 3 x 10~4 Light Levels Dark O L u x Bright 380 L u x Light Noise Qt 0 0 2 QL 0 0 50 Measurements (none) have a well entrained circadian phase prior to entering a laboratory study. The distribution of the final phase probability distribution will give an indication as to the strength of that assumption. To illustrate the effects of noise parameters, three different simulations are run using the same baseline schedule. In the first case no noise sources are included, in the second case process noise is added to the phase state, and the third case a further random variability is introduced to the light input. The light variability is introduced by a Gaussian random vari-ation in the timing of each light to dark transition with variance Qt, and a second Gaussian random variation in the light level with variance Qt- The parameters for each of the three simulations are given in Table 6.3. To plot the results from these simulations that have 240 particles, the probability dis-tributions will be shown rather than each particle's trajectory as done with 24 particles in Fig. 6.2. Over the course of the 21 day period, the particle probability distributions for simu-lation with no noise (2A) are shown in Fig. 6.5 as contour plots. The same phase probability distribution is shown in Fig. 6.6 in an alternate representation. As in the previous scenario the most significant interpretation here can be gained from the final probability distribution. A t day 21 the maximum likelihood of the phase occurs around 4:00, which corresponds to the 60 Chapter 6. S imula t ion Results Amplitude State, A 8 10 12 14 Phase State, if 18 20 -~i n -i r _i i_ 8 10 12 14 Light Response State, n 20 Figure 6.5: Particle probability distributions for the three state variables during Simulation 2A. expected result for an individual synchronized to a regular sleep-wake cycle. When process noise is introduced, in Case 2B, there is a broadening of the probability distribution which can be seen in Fig. 6.7. The particle probability distributions from the final simulation case, with additional variability in the light input, is shown in Fig. 6.9. This case shows a further widening of the distribution which represents additional uncertainty in circadian phase location. A comparison of the final phase probability distributions from each of the three cases is shown in Fig. 6.2. The successive decrease of the certainty with the introduction of each noise term is immediately apparent. One interesting observation is that the introduction of process noise in case IB reduced the peak probability but maintained the same most likely value as the case with no noise. The introduction of light variability however had the interesting effect of skewing the most likely value forward. 61 to w ct> Cfi> CO T3 o cr w cr a . cn' c-l-3. o" c r+ M • o 3 cn 3> "! rt-3 d CD 02 *j E. 3 05 re (5 5. w re CO a 3 3. 3 (K5 a : —. 0 3 Phase (Hours) Amplitude 0 3 NO > 05 05 o CO •O 3" — ft £. T3 -! o c r c r 3 fB T) "! CD cn rt> 3 O 3 3 d CD •3 D" pa CO re re < 5. po re 2 -3 Probability, p(<t>) 3" -ft -! 05 GO 3 c re —• o 3 s o to Peg CD o> co TJ 9-T3 "( o cr » cr If M CD 3. cr c r>-I-I • O 3 01 o> 1 in-=r n ri-cr ~. CD CD •< 09 3. 03 cr CO a. c 3. D W — • 3 Phase (Hours) -* -* ro O Oi O CJl o Amplitude 3 -O 3 to t o bo a *i 09 T3 jT o" SL -i o er 09 cr to 3 CD TJ *1 CD CO n> a o 3 cr CD =r | CO It < 09 5. 09 p£ CtT E-Probability, p(<t>) o 3-BP 13 rc -. 05 — • 3 c_ s S' s H n ai s o c 35 to o 3 3 Chapter 6. Simulation Results a. 0.3 J . . : ' : : : : Phase, <|> (hours Figure 6.10: Graphical probability density representation of the phase state variable during Simulation 2C. Simulation 1A Simulation IB Simulation 1C • -f I : j ' Iv-. h \ 01 \ \ 0 2 4 6 8 10 12 14 16 18 20 22 24 Circadian Phase, <(> (hours) Figure 6.11: Comparison of final phase probability distributions after 21 days of known light level, in the case with no noise (2A), process noise on the phase state (2B), and process noise plus light input noise (2C). 64 Chapter 6. Simulation Results 6.3 Scenario 3 To explore phase shifting properties of the circadian pacemaker model, we can create a sce-nario of an individual making an abrupt eight hour advance in their sleep/wake schedule as can happen with transmeridian airplane travel. As shown experimentally [4], and as predicted by the circadian model [38], the introduction of strong light pulses at appropriate circadian phases will increase the rate of adjustment to a shifted schedule. The greatest ef-fect is achieved by introducing light near the nadir of the circadian when the light response is most sensitive, however the nadir is also the critical point at which shifts occur in opposite directions on either side. Therefore, by introducing light near the nadir there is a risk of causing a shift in the opposite direction to the one desired. With the particle filter we can test the probabilistic outcomes of phase shift direction. A bright light pulse is chosen with a duration of five hours, starting at time 3:00. Clock Time Figure 6.12: Light schedule for simulation 2. The simulation parameters in Table 6.4 are identical to Simulation 2B except that the ini-tial condition is set to a known prior distribution that represents an individual well entrained to the initial schedule. As seen in the trajectories of the phase state in Figure 6.13 there is a split in the probabil-ity distribution after the first introduction of the light pulse. The optimal adjustment to the 65 Chapter 6. Simulation Results T a b l e 6.4: S i m u l a t i o n 3 P a r a m e t e r V a l u e s P a r a m e t e r V a l u e N o . P a r t i c l e s 2 4 0 I n i t i a l D i s t r i b u t i o n G a u s s i a n W(4 .38,1.3) P r o c e s s N o i s e , Qe 3 x L i g h t L e v e l s D a r k T y p i c a l B r i g h t M e a s u r e m e n t s O L u x 380 L u x 10000 L u x (none) new schedu le occurs on t h e p h a s e t ra jec to ry w h i c h dec reases to ze ro a n d t h e n w r a p s a r o u n d to s t a b i l i z e a t 2 2 h by day t w e l v e ( ten d a y s a f te r t he sh i f t ) . T h e o the r sh i f t , i n w h i c h phase i nc reases doesn ' t e n t r a i n to t h e n e w schedu le u n t i l a b o u t d a y 23 ( twen ty -one d a y s a f te r t he sh i f t ) , w h i c h i s t w i c e as l o n g . W h i l e t h e d i v e r g e n t sh i f t b e h a v i o u r h a s b e e n c a p t u r e d i n t h e c i r c a d i a n p a c e m a k e r m o d e l , the c a p a b i l i t y to a n a l y z e p r o b a b i l i s t i c s c e n a r i o s i n t h i s m a n n e r h a s no t been poss ib le p rev ious l y . Amplitude State. A 10 12 14 Phase State,* 18 20 8 10 12 14 Light Response State, n 16 18 20 22 F i g u r e 6 .13 : P a r t i c l e t ra jec to r ies for t h e t h ree s ta te v a r i a b l e s d u r i n g S i m u l a t i o n 3 . 66 Chapter 6. Simulation Results Phase, (j> (hour: Figure 6.14: Graphical probability density, representation of the phase state variable during Simulation 3. 67 Chapter 6. Simulation Results 6.4 Scenario 4 The previous two scenarios involved predicting the effect of known light input on the circa-dian pacemaker, however the additional capability to incorporate feedback from measure-ments has not been illustrated. We can revisit the previous scenario in which the model predicted a bimodal probability distribution resulting from the introduction of bright light stimulus around the circadian nadir. That situation would lead to difficulty making decisions based on the individual's state since the two trajectories lead to two possible scenarios which are diametrically opposed in their physiological circadian effects. The opportunity for physi-ological measurements to enhance state estimates is highlighted in this scenario as a small set of measurements could resolve the ambiguity between two possible trajectories. To test this hypothesis, we can run a simulation with the same conditions as Scenario 3 but with the addition of measurement points. Clock Time Figure 6.15: Light schedule for simulation 4. Let's choose the situation in which the individual in question is actually phase advancing (on the lower trajectory) then add two measurements at day 6 and day 8 and assume accurate mean measurements, but allowing relatively low measurement precision a2 = 3. Comparing the trajectories of the phase state in Figure 6.16 to that of the scenario 2, 6.16, it is apparent that the measurements have provided sufficient information to determine that 68 Chapter 6. Simulation Results Table 6.5: Simulation 4 Parameter Values Parameter Value No. Particles 240 Initial Distribution Gaussian W(4.38.1.3) Process Noise, Qg 3 x 10~'1 Light Levels Dark Typical Bright Measurements O L u x 380 Lux 10000 Lux 9{M) = 1 ± 3 6»(8<i) = 22 ± 3 the individual has a high probability of being on the lower trajectory. Looking to ambulatory, non-laboratory settings where sensors will not have the same degree of precision as in con-trolled environments, the capability to integrate multiple somewhat noisy measurements in this manner will be advantageous. Amplitude State. A 10 12 14 Phase State,» 8 10 12 14 Light Response State, n 16 18 20 22 10 12 14 16 18 20 Time (days) Figure 6.16: Particle trajectories for the three state variables during Simulation 4. 69 Chapter 6. Simulation Results Phase, $ (hour: Figure 6.17: Graphical probability density representation of the phase state variable during Simulation 4. 6.5 Summary-Through a series of four simulated scenarios, characteristics of the particle filter and the properties of the circadian pacemaker model have been demonstrated. A n appropriate num-ber of particles was selected by looking for convergence in the probability distributions. Based on an assumption of a process noise level a default value of 240 particles was selected. Differ-ent noise models were also shown to affect the width of the phase probability distributions. A confirmation that the particle filter can represent nonlinear aspects of the circadian pace-maker model was shown in a scenario designed to elicit a bimodal probability distribution. Finally, the capability to introduce measurements to enhance the predictions was tested. 70 Chapter 7 H u m a n S tudy Resu l ts The majority of previous human studies of circadian rhythm physiology have been designed with the objective of obtaining accurate phase estimates through core body temperature or melatonin measurements. To investigate the feasibility of using the novel circadian moni-toring algorithms in conjunction with non-invasive sensors to monitor circadian phase a new human study was conducted to build a data set with simultaneous measurements from a wide array of non-invasive sensors in addition to the standard invasive sensors for compari-son purposes. The study was conducted at the Center for Study and Treatment of Circadian Rhythms and Douglas Hospital in Montreal, Canada. Subjects were monitoring for a multi-day period with continuous monitoring from sensors on their bodies. To process the collected data, arti-facts are removed and then a particle filter is configured for the specific sensor channels that will be analyzed. The results demonstrate the application of the circadian monitor approach to an real individual and comparisons are made to the traditional circadian estimation ap-proaches. 7.1 Equipment A large ensemble of sensor equipment was utilized to physiologically monitor subjects during the period of time they spent in the laboratory. To meet the data collection objectives, sensors 71 Chapter 7. Human Study Results from a variety of different sources were used simultaneously. Each of the sensor components is described below followed by an explanation of the laboratory facilities. 7.1.1 Sensors The locations of the temperature sensors and a list of the additional non-temperature sensors used for each subject are shown in Fig. 7.1. To achieve simultaneous measurements on all Temperature Sensors Other Sensors Respiration E C G Posture Activity SP02 Figure 7.1: Summary of all physiological sensors installed during the human study. of the physiological processes four different sensory system were utilized. The MiniLogger and DAS systems were standard equipment for the laboratory and the Lifeshirt and Cerebral temperature sensor had not previously been used at this lab. 72 Chapter 7. Human Study Results MiniLogger The MiniLogger (Mini-Mitter Inc., Sunriver, Oregon USA) is a commonly used device for monitoring skin and core temperatures. One MiniLogger was worn by each subject to record four skin temperature locations. D A S A core body temperature measurement using a YSI-400 rectal probe was recorded by the DAS computer installed in the laboratory. The DAS also measured ambient environmental conditions of light intensity and air temperature in each subjects' room. LifeShirt The Lifeshirt-400 (Vivometrics Inc., Ventura, California, USA) is an ambulatory monitoring system consisting of a light weight vest with an integrated respiratory inductive plethysmo-graph, three-probe cardiac E C G sensor, a two-axis motion sensor, a skin temperature sensor and blood oxygen sensor. Data is stored on a digital computer stored in a hip worn pouch with data cards that were swapped out at twelve hour intervals. Cerebral Temperature Sensor The cerebral temperature sensor is an active thermal sensor which was a single prototype developed by Dr. Dittmar from Lyons, France. It is a noninvasive temperature sensor which consists of a rubber pad with a surface areas of approximately 2" x 3" which is placed on the side of the head. An active heating element creates an isothermal boundary condition and thus "draws-out" temperatures 2 to 3cm deep up to the surface. A skin surface probe can then be used to measure temperature at the surface as shown in Fig. 7.2 Located on the side of the temple this results in the measurement of the brain temperature. A photograph of the cerebral sensor is shown in Fig. 7.3. Since this sensor was a prototype it was not suited for ambulatory use therefore measurements were only possible during the constant routine 73 Chapter 7. Human Study Results Figure 7.2: Diagram i l lus t ra t ing the placement of the cerebral temperature sensor on the side of the head and the isothermal region which i t creates wi th the brain. Figure 7.3: Picture o f the cerebral temperature sensor equipment, wi th data acquisition unit (above) and skin mounted sensor (below). period when the subjects were i n a stationary location. The cerebral temperature sensor had previously been used f o r short studies o f durat ion o f less than a f e w hours, so this protocol produced the first data set f o r this sensor which exceeded twenty f o u r hours. Hormonal Assays Salivary melatonin and C o r t i s o l assays were collected periodically dur ing the course o f the study f o r subsequent chemical assays. 74 Chapter 7. Human Study Results 7.1.2 Laboratory The facilities of the Centre for the Study and Treatment of Ci rcadian Rhythms include time isolation rooms which allow subjects to be maintained for multi-day periods inside a com-pletely controlled environment. L igh t levels and air temperature are accurately controlled, and the subjects are continuously monitored by staff us ing closed circuit television. A photo-graph of one of the rooms is shown i n Fig . 7.4. Figure 7.4: Controlled-environment room at the Center for the Study and Treatment of Ci r -cadian Rhythms, Douglas Hospi tal , i n which subjects l ived dur ing the five day laboratory protocol. 7.2 Protocol Six healthy young subjects were selected as participants i n the study. Toxicological analysis of urine was performed on each subject to confirm that subjects were drug-free at the t ime of study. None had a history of t ransmeridian travel across more than one time zone or of shift work i n the 3 months preceding the start of the study. Subjects gave informed consent i n accordance wi th the guidelines of the Douglas Hospi ta l Research Ethics Board and the U B C 75 Chapter 7. Human Study Results R e s e a r c h E t h i c s B o a r d . T o co l lec t d a t a se ts o f su f f i c ien t l e n g t h to e n a b l e i d e n t i f i c a t i o n o f t h e c i r c a d i a n p rocesses a five d a y l a b o r a t o r y pro toco l w a s se lec ted as s h o w n i n F i g . 7.5. T h e a m b u l a t o r y d a y s a l l o w for the a n a l y s i s o f d a t a i n a m b u l a t o r y c o n d i t i o n s so s u c h t h a t some degree o f ' m a s k i n g ' fac to rs w o u l d be p r e s e n t w h i l e t h e final c o n s t a n t r o u t i n e a l l o w s for a s u b s e q u e n t p r e c i s e e s t i m a t e o f c i r c a d i a n p h a s e . Clock Time 4 8 12 16 20 24 _ l l l I l l l l I i _ a 2 ] Ambulatory, 150 Lux ^| Constant Routine, 10 Lux I Sleep, 0 Lux F i g u r e 7.5: A c t i v i t y a n d l i g h t s c h e d u l i n g p ro toco l for l a b o r a t o r y s tudy. Preparatory Monitoring (Days -21 to 0) I t i s i m p o r t a n t to e n s u r e t h a t sub jec ts h a v e been e n t r a i n e d to a r e g u l a r s l e e p - w a k e cyc le p r i o r to t h e i r e n t r a n c e i n t o the labora to ry . F o r t h r e e w e e k s p r i o r to e n r o l m e n t i n t he s tudy, e a c h sub jec t w a s a s k e d to m a i n t a i n a r e g u l a r schedu le o f 8 -hou rs of s l e e p / d a r k n e s s a n d con-s i s t e n t l y m a i n t a i n t h e i r h a b i t u a l s leep onse t t i m e . T h e y m a i n t a i n e d a s l e e p / w a k e l o g a n d c o m p l i a n c e to t h e i r s chedu le w a s be ve r i f i ed the w e e k p r i o r to a d m i s s i o n by a w r i s t a c t i v i t y m o n i t o r i n g w i t h a so l i d s ta te , po r t ab l e d a t a co l l ec t ion dev i ce a n d a l i g h t de tec to r ( A c t i w a t c h w - L , f r o m M i n i - M i t t e r Inc. , S u n r i v e r , O r e g o n U S A ) . Ambulatory Laboratory (Days 0 to 3) Sub jec t s e n t e r e d t h e l a b o r a t o r y o n the e v e n i n g o f D a y 1 a n d w e r e e q u i p p e d w i t h t he senso rs t h a t reco rded c o n t i n u o u s l y t h r o u g h the e x p e r i m e n t a l p r o c e d u r e . D u r i n g w a k i n g ep isodes the r o o m l i g h t w a s m a i n t a i n e d a t a t y p i c a l i n d o o r l eve l o f 150 L u x , a n d the sub jec ts w e r e a l l o w e d to f ree l y m o v e i n t h e r o o m a l t h o u g h w e r e not p e r m i t t e d to engage i n s t r e n u o u s ac t i v i t y . T h r e e 76 Chapter 7. Human Study Results consecutive nights of sleep occurred. Constant Routine Laboratory (Days 4 to 5) To acquire a reliable measurement of the subjects circadian phase against which the previ-ous measurements can be correlated a 32-hour constant routine (CR) was performed in the final two days of the study, starting on the morning of day 3 and end in the evening of day 4. The CR procedure is designed to limit the masking effects that light exposure, postu-ral changes, activity, and meal intake may have on the endogenous expression of measured circadian markers. In this procedure, each subject was asked to remain awake, in a semi-recumbent posture under dim light conditions ( 10 lux). Throughout the CR, they remained in the required posture and were relatively inactive. The cerebral temperature sensor was added to the existing sensors at the start of the constant routine. The CR was followed by an ad libitum recuperative sleep episode. 7 . 3 Data Processing To prepare for the analysis of the sensor measurements collected from the study, the data is first cleaned to remove artifacts and major problems and then the particle filtering system is configured. 7.3.1 Artifact Removal and Denoising On all of the data signals artifacts from known causes and outliers were first manually iden-tified. Known sources of outlier patterns primarily consisted of sensor dislocation or removal or external stimulus such as a shower. Other variability in the data that was unlikely to reflect true physiological measurements but could not be attributed to a likely cause was treated more conservatively. Sections of outliers that were removed were either left blank or linearly interpolated depending on the situation. Following manual cleaning the data signals were denoised using a wavelet filter from the Matlab wavelet toolbox. An example of a signal 77 Chapter 7. Human Study Results before and after denoising is shown in Fig. 7.6. 38 I 1 1 36 I—! 1 ! 1 U 1 ! 1 L_J 00:00 00:00 00:00 00:00 00:00 00:00 00:00 00:00 00:00 00:00 Time (Clock Hours) Figure 7.6: Temperature data signal before (top) and after (bottom) artifact removal and denoising. To provide a general illustration of the data collected during the study, the temperature signals from one of the subjects over the five day period are shown in Fig. 7.7. A challenge which presented itself in the data analysis was that in five out of six subjects there were significant artifacts which rendered large segments unusable. An important signal which was impacted was core body temperature signal which was intended to provide a key base-line. The cause of the noisy signals was not known, but one likely possibility is that the high encumbrance on each subject due to the unusually high number of sensors may have contributed to sensor displacements. 78 Chapter 7. Human Study Results Head Temp. Foot Temp. Leg Temp. Rectal Temp. — Cerebral Temp. Chest Temp. 00:00 00:00 Figure 7.7: Temperature data from subject 1 showing all skin and core temperatures. The visible floor e f f e c t on the hand temperature channel is a result of a sensor range limitation. This channel is not used in an analysis. 7.3.2 Particle Filter Configuration The circadian pacemaker model at the core of particle filtering algorithm does not have any customized tuning parameters so it remains a fixed component. For a given set measure-ments however the physiological marker components must be selected, and there are a num-ber of different choices to consider. To set up a particle filter for the analysis of data from subject 1 physiological signals are chosen from three categories for comparison. A set of invasive measurements includes rectal core body temperature, salivary melatonin, and salivary C o r t i s o l . A non-invasive but non am-bulatory signal is from the cerebral temperature sensor, and a non-invasive and ambulatory signal is heart rate variability as collected by the ECG on the Lifeshirt. Phase markers for the temperatures and hormone levels can all be determined with direct 79 Chapter 7. Human Study Results application of a Fourier fit combined with either a maxima or minima location and time shift. Assuming that the cerebral temperature follows the same core body temperature as the rectal sensor, both can be analyzed with a second-order Fourier fit and a minimum point. The necessary time shift to correlate to the central circadian phase is defined by the Kronauer model to be 0.8 hours. Therefore the shift from core body temperature minima to circadian phase is TCBT = -0.8h. With melatonin and Cortisol the maxima of a second order Fourier fit is selected however there are a variety of different phase analysis options [34] which could be equally well se-lected. The relative time shifts between melatonin and Cortisol and core body temperature have not been extensively studied, however based on published results in [34] values can be determined. The melatonin maxima occurs approximately 2.3 hours before core body tem-perature minima, therefore a relationship the central pacemaker is TMET = 2.3 — 0.8 = 1.5. The Cortisol maxima occurs approximately 2.2 hours after core body temperature minima therefore its central shift is rcort = —2.2 — 0.8 = —3.0. The selection of heart rate variability (HRV) parameters is less clearly defined than the other measures. Heart rate variability itself has multiple calculation methods and interpre-tation, and the degree to which H R V reflects an endogenous circadian component is also not definitive (see Background). Our experimental protocol does not enable separation of sleep-wake effects during the first days, however in the interest of including a fully non-invasive, non-ambulatory sensor into the analysis, H R V is analyzed for phase markers based on the waking E C G segments for each of the two preliminary study days. Respiratory and cardiac measures were analyzed by the Vivologic analysis software (Vivometrics) associated with the Lifeshirt. This software performs HR calculation and peak detection algorithms. Respiratory sinus arrythmia (RSA) is used as a HRV metric. The RSA is quantified in the time domain using the peak-valley algorithm to estimate amplitude of RSA [22]. This method has been 80 Chapter 7. Human Study Results previously validated as highly correlated with spectral analysis estimations [23]. The min-ima of the Fourier fit to the calculated heart rate is then determined and introduced with a phase shift O(THRV = -12 - 0.8 = -12.8. The configuration of the particle filter estimator with all of the physiological marker com-ponents is shown in Fig. 7.8. CBT Cerebral Temp Melatonin CoriAsol ECG Human fourier2 min Shift -0.8 h <P<:UT fourier2 min Shift -0.8 h <P,, fourier2 max Shift 4-1.5h <pu,, fourier3 max Shift -3.0 h <Pc„ H R V fourier2 Shift •f. min —*> -12.8 h Physiological Phase Estimate Measurement Update Prediction Update State Estimate Figure 7.8: Configuration of the particle filter estimator for processing the data collected from subject 1. Time shifts and physiological marker detection functions are selected for each signal. 81 Chapter 7. Human Study Results 7.4 Results A highlight from the data was the performance of the cerebral temperature sensor which had not previously been tested in a long duration study. Further discussion of the cerebral sensor results is followed by results from the application of the particle filtering algorithm to one of the subject data sets. 7.4.1 Cerebral Temperature Sensor Results from the first subject showed a strong linear correlation between the measurement produced by the cerebral temperature sensor and the core body temperature measured by the rectal probe. In this trial the relationship was continuous over the 32 hour constant routine. Subsequent trials had significantly higher artifacts in both the rectal and cerebral temperature. To facilitate observation of the cerebral temperature sensor response across all subjects, large artifacts were removed from both the cerebral and rectal temperature signals and plot-ted without any further denoising in Fig. 7.9. Due to the large number of data segments removed in an ad-hoc manner a statistical analysis was not conducted to compare the two signals. One interesting observation is that the cerebral temperature signal for subjects 3 and 4 actually appears more consistent with the expected core body temperature trajectory than the rectal temperature probe. Due to the invasiveness of the rectal temperature sensor and the current lack of accurate, non-invasive alternatives suitable for extended use, the initial results shown here strongly support further development and testing of the cerebral temperature sensor. 82 Chapter 7. Human Study Results 1) 2) 37 00:00 06:00 12:00 18:00 00:00 06:00 12:00 18:00 00:00 06:00 12:00 18:00 00:00 3) 4) 38.5 38 37.5 37 36.5 36 35.5 Cerebral Temp. Rectal Temp. 5) 38.5 38 37.5 37 36.5 36 35.5 Cerebral Temp. Rectal Temp. Figure 7.9: Cerebral temperature sensor measurements and rectal core body temperature measurements plotted for each of the five subjects with major artifacts removed. 83 Chapter 7. Human Study Results 7.4.2 Particle Filter Estimation Test T h e d a t a set f r o m sub jec t 1 i s v e r y c l e a n the re fo re i s se r ves as a good case for e x a m i n i n g t h e p e r f o r m a n c e of t he pa r t i c l e f i l te r a l g o r i t h m on a r e a l d a t a set . A m o t i v a t i n g ob jec t ive i s t he d e t e r m i n a t i o n o f p h a s e e s t i m a t e s b a s e d s t r i c t l y on n o n - i n v a s i v e a n d a m b u l a t o r y m e a s u r e -m e n t s so a c o m p a r i s o n w i l l be m a d e b e t w e e n the i n v a s i v e a n d n o n i n v a s i v e m e a s u r e s . T h e c o m p a r i s o n w i l l be m a d e a t t h e e n d o f D a y 5. 7.4.2.1 Phase estimates from noninvasive, ambulatory data T h e a v a i l a b l e i n f o r m a t i o n for Sub jec t 1 w h i c h c a n be a p p l i e d to gene ra te p h a s e e s t i m a t e s i n c l u d e s the s leep h i s t o r y p r i o r to the sub jec t e n t e r i n g the l a b , t he k n o w n l i g h t e x p o s u r e d u r i n g the l a b p ro toco l , a n d the n o n i n v a s i v e d a t a co l lec ted d u r i n g a m b u l a t o r y d a y s 2 a n d 3. T h e pro toco l i s s h o w n i n F i g . 7.10. 0 4 8 12 16 20 24 F i g u r e 7.10: S t u d y p ro toco l s c h e d u l e for sub jec t 1. F i r s t w e c a n gene ra te a n e s t i m a t e b a s e d s t r i c t l y o n k n o w i n g the l i g h t e x p o s u r e d u r i n g the l a b study. In t h i s case t h e subject 's c i r c a d i a n p h a s e p r i o r to e n t e r i n g t h e l ab i s c o n s i d e r e d comp le te l y u n k n o w n . T h e r e s u l t i s s h o w n o n t h e le f t o f F i g . 7.11 A s e c o n d e s t i m a t e c a n be 84 Chapter 7. Human Study Results generated by adding knowledge of the subject's sleep history which was kept in a log prior to arriving at the lab. A number of assumptions are made in this analysis since the exact light levels are not known, and uncertainty surrounding the accuracy of the sleep log is unknown. Light was assumed to be constant at 380 Lux for the duration of waking. The new estimate, shown on the right in Fig. 7.13, has a narrower distribution, reflecting the additional information. Time (days) Phase, <|> (hours) Time (days) Phase, <]) (hours) Figure 7.11: Phase probability distribution generated only from known light levels during the study (left), and with the addition of known sleep history prior to entering the lab (right). To build further information into the estimate a noninvasive measurement of heart rate variability phase is extracted from the two ambulatory days as shown in Fig. 7.12. The phase estimates generated just from using just the two H R V measurements are shown in the left of Fig. 7.13, and the combined estimate using all of the sleep history and the HRV measures is shown on the right in Fig. 7.13. 7.4.2.2 Phase estimates from constant routine data During the constant routine monitoring from all the sensors used during the ambulatory periods was maintained with the addition of the cerebral temperature sensor and melatonin, 85 Chapter 7. Human Study Results 0 4 8 12 16 20 24 0 1 2 3 4 i>! * if F i g u r e 7.12: S t u d y p ro toco l s chedu le for sub ject 1 w i t h r a w a n d f i t ted H R V d a t a s u p e r i m -posed . a n d C o r t i s o l assays . A n i l l u s t r a t i o n o f t he r e c t a l t e m p e r a t u r e a n d m e l a t o n i n s i g n a l s on top of the m u l t i - d a y schedu le i s s h o w n i n F i g . 7.14. P r e c i s e p h a s e m a r k e r s w e r e c a l c u l a t e d f r o m the two core t e m p e r a t u r e s a n d t w o h o r m o n e a s s a y s as s h o w n i n F i g . 7.15. A f t e r t h e d e r i v i n g t he p h a s e m a r k e r s , e a c h p r o b a b i l i t y d i s t r i b u t i o n i s s h i f t e d a c c o r d i n g to t he c a l i b r a t i o n se t t i ngs p rec ious d e t e r m i n e d . T h e p h a s e m a r k e r s before a n d a f te r s h i f t i n g a r e s h o w n toge the r i n F i g . 7.16. T h e t e m p e r a t u r e m e a s u r e m e n t s a n d m e l a t o n i n m a r k e r s a r e h i g h l y cons i s ten t , w i t h d i f fe rence obse rved i n t he Cor t i so l . Comparison A c o m p a r i s o n b e t w e e n the f i n a l p h a s e e s t i m a t e s g e n e r a t e d a t t h e end of t he of s t u d y b y the s t a n d a r d core body t e m p e r a t u r e a n d m e l a t o n i n m e t h o d s a n d t h e e s t i m a t e s g e n e r a t e d b y t h e pa r t i c l e f i l t e r a n d t h e n o n i n v a s i v e H R V m e a s u r e m e n t i s s h o w n i n F i g . 7.17 T h e p r e d i c t i o n s g e n e r a t e d b y the p a r t i c l e f i l t e r e s t i m a t i o n t h a t u s e d pas t l i g h t e x p o s u r e a n d h e a r t ra te v a r i -a b i l i t y l i es b e t w e e n the core body t e m p e r a t u r e , m e l a t o n i n a n d Cortisol m a r k e r s . 86 Chapter 7. Human Study Results 0.6 v . -Time (days) o Phase, <|> (hours) i e (days) o Phase, <|> (hours) Figure 7.13: Phase probability distribution generated from light levels during the study and H R V measures (left), and with the further addition of known sleep history prior to entering the lab (right). 7.5 Summary A human study was designed and conducted to collect data sets suitable for evaluating the feasibility of noninvasive circadian phase assessment in conjunction with new estimation methods. While artifacts in the data limited the potential for analyzing the data from all subjects, the particle filtering approach was testing on the data from one of the subjects. Two notable observations resulted. First, the phase estimates from the prototype of a new non-invasive cerebral temperature sensor provided a nearly identical estimate to the that generated from the standard core body temperature assessment. This indicates a promising avenue for further sensor development. Second, a comparison showing the phase assessment that can be obtained from strictly non-invasive means is compared to that generated from the invasive measurements under controlled constant routine conditions. The results from the controlled conditions were more precise, however the particle filtering option still shows significant estimates which could be a desired trade-off for the advantage of non-invasiveness. 87 Chapter 7. H u m a n Study Results F i g u r e 7.14: C o r e body t e m p e r a t u r e (left) a n d s a l i v a r y m e l a t o n i n ( r igh t ) r e c o r d e d d u r i n g t h e cons tan t r ou t i ne . Chapter 7. Human Study Results 8 12 16 20 0 4 8 12 16 12 16 20 0 4 8 12 Clock Time (hours) Clock Time (hours) 8 12 16 20 0 4 8 12 8 12 16 20 0 4 8 12 Clock Time (hours) Clock Time (hours) 8 12 16 20 0 4 8 12 8 12 16 20 0 4 8 12 Clock Time (hours) Clock Time (hours) Figure 7.15: Physiological markers collected from Subject 1 during the constant routine shown with the probability of peak locations shown below plots of the raw and Fourier fitted for each. Clockwise from the top left are rectal temperature, cerebral temperature, Cortisol and melatonin. 89 Chapter 7. Human Study Results Peak Marker 0.6 m 0.5 m Q-o 0.4 Cerebral temp. Rectal temp. Cortisol Melatonin I A I \ i hi l i J ^1 12 16 20 0 4 8 12 Clock Time (hours) Figure 7.16: Physiological markers assessed during the constant routine as determined from the preliminary peak distributions (top) which are then adjusted with corresponding time shifts to generate estimates of central circadian phase (bottom). 90 Chapter 7. Human Study Results Shifted Circadian Phase Marker Clock Time (hours) F i g u r e 7 .17: C o m p a r i s o n b e t w e e n p h a s e e s t i m a t e s g e n e r a t e d by c o n s t a n t r o u t i n e ( C R ) core body t e m p e r a t u r e a n d m e l a t o n i n m e a s u r e m e n t s , a n d p a r t i c l e f i l t e r e s t i m a t e s b a s e d on A ) n o n e B ) H R V O P r i o r H i s t o r y D ) H R V + P r i o r h is to ry . E ) C e r e b r a l t e m p e r a t u r e 91 Chapter 8 Conclusion The goal of this project was to apply modelling and signal processing tools to enhance hu-man circadian physiology monitoring. A model framework was developed that allows the integration of physiological measurements with prediction models, which were two compo-nents previously treated independently. A complementary particle filter algorithm was then developed that generates real-time circadian phase estimates that presents rich probability distributions. While substantial validation of the algorithm against human data has not yet be performed, initial results from a subject in the study that was conducted showed promising results for the algorithm and for a noninvasive cerebral temperature sensor. 8.1 Significance of Results The model-based Bayesian approach to circadian physiology monitoring presents a new set of capabilities to those who are interested in addressing applied problems related to human circadian physiology. Even without the use of physiological sensors, the simulation capabil-ity alone provides the capability for analyzing the evolution of circadian phase probability distributions through any configuration of arbitrary light exposure scenarios. The particle filter implementation allows for full modelling of nonlinearities present in the human circa-dian system, and as such allows interesting emergent properties such as bimodal probability distributions to be accurately captured. The probabilistic framework furthermore provides 92 Chapter 8. Conclusion the capability to evaluate scenarios in which there is not perfect information about the initial conditions of an individual's state, as is often the case is practice. This work also introduces a method for combining physiological circadian marker mea-surements in an integrated way with the predictive model. Through the recursive Bayesian filtering framework, predictions in the evolution of an individual's circadian phase can be made using the existing circadian pacemaker model, and updated with information resulting from any direct measurements. The significance of this capability is that it directly addresses the challenge of operating in environments with sparse and noisy measurements. Finally, a sensor that shows promise as a new method for noninvasively assessing core body temperature has been tested for the first time in an extended duration trial. While noisy data and a relatively small sample size do not enable a conclusive statement, the cere-bral brain temperature sensor showed a very strong correlation to the standard rectal tem-perature sensor. 8.2 Summary of Contributions The key novel contributions that have been made by this project are: • development of a nonlinear model transformation for a commonly used circadian pace-maker model which allows the circadian phase and amplitude to be treated as continu-ous states in a state space model • introduction of a recursive Bayesian filtering technique to the field of circadian phys-iology which is the first method to allow direct measurements of circadian state to be combined with model predictions • introduction of the capability for a rigorous treatment of uncertainty and probability distributions to circadian modeling and prediction. 93 Chapter 8. Conclusion • testing of a wide array of noninvasive sensors for circadian monitoring and shown re-sults from a cerebral temperature sensor which indicate the potential for an accurate noninvasive alternative for core body temperature measurement 94 Bibliography [1] T. Akerstedt, S. Folkard, and C. Portin, Predictions from the three process model of alert-ness, AVIATION SPACE A N D E N V I R O N M E N T A L MEDICINE 75 (2004), no. 3, Suppl., A75-A83. [2] B.D. Anderson and J.B. Moore, Optimal filtering, Prentice-Hall, New Jersey, 1979. [3] S. Arulampalam, S. Maskell, N . Gordon, and T. Clapp, A tutorial on particle filters for non-linear / non-Gaussian Bayesian tracking, I E E E Trans. Signal Processing 50 (2002), 174-188. [4] D.B. Boivin and F.O. James, Phase-dependent effect of room light exposure in a 5-h ad-vance of the sleep-wake cycle: Implications for jet lag, 17 (2002), no. 3, 266-267. [5] A. Borbly and P. Achermann, Sleep homeostasis and models of sleep regulation, Journal of Biological Rhythms 14 (1999), no. 6, 557-568. [6] E .N. Brown, Y. Choe, H . Luithardt, and C A . Czeisler, A statistical model of the human core-temperature circadian rhythm, American Journal of Physiology - Endocrinology and Metabolism 279 (2000), E669-E683. [7] E.N. Brown and C A . Czeisler, The statistical analysis of circadian phase and amplitude in constant-routine core-temperature data, Journal of Biological Rhythms 7 (1992), no. 3, 177-202. 95 B I B L I O G R A P H Y [8] H.J. Burgess and C.I. Eastman, The dim light melatonin onset following fixed and free sleep schedules, Journal of Sleep Research 14 (2005), no. 3, 229-237. [9] R.E. Burns, Biological time and in vivo research: A field guide to pitfalls, The Anatomical Record 261 (2000), 141-152. [10] J.A. Caldwell, Fatigue in aviation, Travel Medicine and Infectious Disease 3 (2005), 85— 96. [11] C A . Czeisler, J.S. Allan, S.H. Strogatz, E.J. Ronda, R. Sanchez, C D . Rios, G.S. Freitag, G.S. Richardson, and R.E. Kronauer, Bright light resets the human circadian pacemaker independent of the timing of the sleep-wake cycle, Science 233 (1986), no. 4764, 667-671. [12] C A . Czeisler, J.F. Duffy, T.L. Shanahan, E .N. Brown, D.W. Rimmer, E.J. Ronda, E.J. Silva, J.S. Allan, J.S. Emens, D.J. Dijk, and R.E. Kronauer, Stability, precision, and near-24-hour period of the human circadian pacemaker, Science 284 (1999), no. 5423, 2177-2181. [13] C A . Czeisler, G.S. Richardson, J. Zimmerman, M . C Moore-Ede, and E.D. Weitzman, Entrainment of human circadian rhythms by light-dark cycles: a reassessment, Photo-chemistry and Photobiology 34 (1981), no. 2, 239-247. [14] D.F. Dinges, Critical research issues in development ofbiomathematical models of fatigue and performance, AVIATION SPACE A N D E N V I R O N M E N T A L MEDICINE 75 (2004), no. 3, A181-A191 Suppl. [15] A . Doucet, N . de Freitas, and N.J. Gordon (eds.), Sequential Monte Carlo methods in practice, Springer, New York, 2001. [16] A . Doucet, S. Godsill, and Andrieu C , On sequential Monte Carlo sampling methods for Bayesian filtering, Statistics and Computing 10 (2000), no. 3, 197-208. 96 BIBLIOGRAPHY [17] J.F. Duffy and D.J. Dijk, Getting through to circadian oscillators:why use constant rou-tines?, Journal of Biological Rhythms 17 (2002), no. 4-13. [18] J. Freitas, P. Lago, J. Puig, M . Carvalho, O. Costa, and A. Falcao de Freitas, Circadian heart rate variability rhythm in shift workers, Journal of Electrocardiology 30 (1997), no. 1, 39-44. [19] W.R. Gilks and C. Berzuini, Following a moving target - Monte Carlo inference for dy-namic Bayesian models, Journal of the Royal Statistical Society, B 63 (2001), 127-146. [20] N.J. Gordon, D.J. Salmond, and A.F .M. Smith, Novel approach to nonlinear/non-Gaussian Bayesian state estimation, I E E E Proceedings-F 140 (1993), no. 2, 107—113. [21] A. Gray and Moore, Very fast multivariate kernel density estimation using via computa-tional geometry, Proceedings of Joint Stat. Meeting, 2003. [22] P. Grossman, J. Karemaker, and W. Wieling, Prediction of tonic parasympathetic cardiac control using respiratory sinus arrhythmia: the need for respiratory control, Psychophys-iology 28 (1991), 201-216. [23] P. Grossman, J. van Beek, and C. Wientjes, A comparison of three quantification methods for estimation of respiratory sinus arrhythmia, Psychophysiology 27 (1990), 702-714. [24] R. Gunawan and F.J. Doyle, Phase sensitivity analysis of a circadian rhythm gene net-work, Proceedings of the 44th IEEE Conference on Decision and Control, 2005, pp. 3687-3692. [25] S.L. Harmer, S. Panda, and S.A. Kay, Molecular bases of circadian rhythms, Annual Review of Cell and Developmental Biology 17 (2001), 215-253. 97 BIBLIOGRAPHY [26] S. Hattar, H . Liao, M . Takao, D.M. Berson, and K. Yau, Melanopsin-containing reti-nal ganglion cells: Architecture, projections, and intrinsic phtotosensiticity., Science 295 (2002), 1065-1070. [27] A. Ihler, Kernel density estimation toolbox for MATLAB, http://www.ics.uci.edu/ ih-ler/code/kde.shtml, 2003, (Aug 2006). [28] H . Ito, M . Nazaki, T. Maruyama, Y. Kaji, and T. Tsuda, Shift work modified the circa-dian patterns of heart rate variability in murses, International Journal of Cardiology 79 (2001), 231-236. [29] A. J. Izenman, Recent developments in nonparametric density estimation, Journal of the American Statistical Association 86 (1991), no. 413, 205-224. [30] M . E . Jewett, D.B. Forger, and R.E. Kronauer, Revised limit cycle oscillator model of human circadian pacemaker, Journal of Biological Rhythms 14 (1999), no. 6, 493-499. [31] M . E . Jewett and R.E. Kronauer, Interactive mathematical models of subjective alertness and cognitive throughput in humans, Journal of Biological Rhythms 14 (1999), no. 6, 588-597. [32] G.A. Kerkhof and H.P. van Dongen, Morning-type and evening-type individuals differ in the phase position of their endogenous circadian oscillator, Neuroscience Letters 218 (1996), no. 3, 153-156. [33] S.B. Khalsa, M . E . Jewett, J .F Duffy, and C A . Czeisler, The timing of the human circa-dian clock is accurately represented by the core body temperature rhythm following phase shifts to a three-cycle light stimulus near the critical zone, Journal of Biological Rhythms 16(1999), no. 6. [34] E.B. Herman, H.B. Gershengorn, J.F. Duffy, and R.E. Kronauer, Comparisons of the 98 BIBLIOGRAPHY variability of three markers of the human circadian pacemaker, Journal of Biological Rhythms 17 (2002), no. 2, 181-193. [35] E.B. Herman, Y. Lee, C.A. Czeisler, and R.E. Kronauer, Linear demasking techniques are unreliable for estimating the circadian phase of ambulatory temperature dat, Journal of Biological Rhythms 14 (1999), no. 4, 260-274. [36] K. Krauchi and A. Wirz-Justice, Circadian rhythm of heat production, heart rate, and skin and core temperature under unmasking conditions in men, American Journal of Physiology - Regulatory, Integrative, and Comparative Physiology 267 (1994), no. 3, R819-R829. [37] R.E. Kronauer, A model for the effect of light on the human 'deep' circadian pacemaker, Sleep Research 16 (1987), no. 621. [38] R.E. Kronauer, D.B. Forger, and M . E . Jewett, Quantifying human circadian pacemaker response to brief, extended, and repeated light stimuli over the photopic range, Journal of Biological Rhythms 16 (1999), no. 6. [39] S. Kullback and R. A. Leibler, On information and sufficiency, Annals of Mathematical Statistics 22 (1951), 79-86. [40] F. Levi, From circadian rhythms to cancer chronotherapeutics, Chronobiology Interna-tional 19 (2002), no. 1, 1-19. [41] F. Levi and M.C. Mormont, Circadian-system alterations during cancer processes: A re-view, International Journal of Cancer 70 (1997), 241-247. [42] F. Lvi , Handbook of experimental pharmacology: Physiology and pharmacology of bio-logical rhythms, ch. Chronopharmacology of anticancer agents, pp. 299-331, Springer-Verlag, Berlin, 1997. 99 BIBLIOGRAPHY [43] , Chronotherapeutics: the relevance of timing in cancer therapy, Cancer Causes Control 17 (2006). [44] F. Lvi , R. Zidani, and J.L. Misset, International organization for cancer chronotherapy: randomized multicentre trial of chronotherapy with oxaliplatin, fluorouracil, andfolinic acid in metastatic colorectal cancer., Lancet 350 (1997), 681-686. [45] M . Mallis, S. Mejdal, T.T. Nguyen, and D.F. Dinges, Summary of the key features of seven biomathematical models of human fatigue and performance, AVIATION SPACE A N D E N V I R O N M E N T A L MEDICINE 75 (2004), no. 3, A4-A14 Suppl. [46] T.H. Monk, What can the chronobiologist do to help the shift worker?, Journal of Biologi-cal Rhythms 15 (2000), no. 2, 86-94. [47] C. Mott, M . Huzmezan, D. Mollicone, and M . van Wollen, Modifying the human circa-dian pacemaker using model-based predictive control, Proceedings of the 2003 American Control Conference, June 2003, pp. 453-458. [48] A . Papoulis, Probability, random variables, and stochastic processes, 2 ed., p. 531, McGraw-Hill, New York, 1984. [49] M . Raciti, P. Pisani, C. Emdin, C. Carpeggiani, S. Ruschi, G. Kraft, R. Francesconi, G. Membretti, and C. Marchesi, Circadian dynamics of respiratory parameters from am-bulatory monitoring, Computers in Cardiology 1994 (1994), 581-584. [50] D.W. Rimmer, D.B. Boivin, T.L. Shanahan, R.E. Kronauer, J .F Duffy, and C A . Czeisler, Dynamic resetting of the human circadian pacemaker by intermittent bright light, A M E R I C A N J O U R N A L OF PHYSIOLOGY-REGULATORY INTEGRATIVE A N D COMPARATIVE PHYSIOLOGY 279 (2000), no. 5, R1574-R1579. [51] A. Samel, M . Vejvoda, and H . Maass, Sleep deficit and stress hormones in helicopter 100 BIBLIOGRAPHY pilots on 7-day duty for emergency medical services, Aviation, Space, and Environmental Medicine 75 (2004), no. 11, 935-940. [52] L . E . Scheving, E.R. Burns, J.E. Pauly, and T.H. Tsai, Circadian variation in cell division of the mouse alimentary tract, bone marrow and corneal epithelium, Anat Rec 191 (1978), 479-486. [53] T.L. Shanahan, R.E. Kronauer, J.F. Duffy, G.H. Williams, and C A . Czeisler, Melatonin rhythm observed throughout a three-cycle bright-light stimulus designed to reset the hu-man circadian pacemaker, Journal of Biological Rhythms 14 (1999), no. 3, 237-253. [54] A . Sitz, U. Schwarz, and J. Kurths, The unscented kalman filter, a powerful tool for data analysis, International Journal of Bifurcation and Chaos 14 (2004), no. 6, 2093-2105. [55] C. Spengler, C A . Czeisler, and S.A. Shea, A n endogenous circadian rhythm of respiratory control in humans, The Journal of Physiology 526 (2000), no. 3, 683-694. [56] Y. Wang and M . Brown, A flexible model for human circadian rhythms, Biometrics 52 (1996), 588-596. [57] J. Waterhouse, S. Kao, D. Weinert, B. Edwards, G. Atkinson, and T. Reilly, Measuring phase shifts in humans following a simulated time-zone transition: Agreement beween constant routine and purification methods, Chronobiology International 22 (2005), no. 5, 829-858. [58] J. Waterhouse, D. Weinert, D. Minors, S. Folkard, D. Owens, G. Atkinson, D. Macdonald, N . Sytnik, P. Tucker, and T. Reilly, A comparison of some different methods for purifying core temperature data from humans, Chronobiology International 17 (2000), no. 4, 539— 566. 101 BIBLIOGRAPHY [59] D. W e i n e r t , A . N e v i l l , R . W e i n a n d y , a n d J . W a t e r h o u s e , The development of new purifi-cation methods to assess the circadian rhythm of body temperature in mongolian gerbils, C h r o n o b i o l o g y I n t e r n a t i o n a l 20 (2003), no. 2 , 2 4 9 - 2 7 0 . [60] M . B . W e i n g e r a n d S . A n c o l i - I s r a e l , Sleep deprivation and clinical performance, J o u r n a l o f t h e A m e r i c a n M e d i c a l A s s o c i a t i o n 287 (2002), no. 8, 9 5 5 - 9 5 7 . 102 Appendix A Research Ethics Board Certificates of Approval 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0064896/manifest

Comment

Related Items