Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Dynamical system approach to heart rate variability: QT versus RR interval
(USC Thesis Other)
Dynamical system approach to heart rate variability: QT versus RR interval
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
DYNAMICALSYSTEMAPPROACHTOHEARTRATEVARIABILITY: QT VERSUSRRINTERVAL by FaribaAriaei ADissertationPresentedtothe FACULTYOFTHEGRADUATESCHOOL UNIVERSITYOFSOUTHERNCALIFORNIA InPartialFulfillmentofthe RequirementsfortheDegree DOCTOROFPHILOSOPHY (ELECTRICALENGINEERING) December2008 Copyright 2008 FaribaAriaei ii Dedication To my husband, Hossein andtomyparents, Ahmad and Mehrangiz iii Acknowledgments Iwouldliketosincerelythankmythesisadvisor,ProfessorEdmondJonckheerefor his valuable and expert academic guidance. This research could not be completed without his constant support and encouragement. This thesis would never have materialized without the constant financial support of the Global Gateway Foundation. In particular, I thank Mr. James Turner and Mrs. Betsy Lehrfeld for their constant support and encouragement. Iwouldalsoliketoexpressmyappreciationtomyfamilyandmyhusband,Hossein, for their patience and encouragement. iv Table of Contents Dedication ii Acknowledgements iii List of Tables vii List of Figures viii Abstract xii Chapter 1: Introduction 1 1.1 The RRInterval ................................. 1 1.2 The QT versus The RRInterval ........................ 8 1.3 Summary ..................................... 11 Chapter 2: Heart Physiology 12 2.1 Introduction.................................... 12 2.2 DynamicsofExcitableCells........................... 12 2.3 ActionPotential ................................. 18 2.3.1 Threshold for Initiation of the Action Potential . . . . . . . . . . . . 21 2.3.2 PlateauinSomeActionPotentials ................... 22 2.3.3 TheHodgkin-HuxleyEquations..................... 23 2.4 ElectricityandtheHeart ............................ 24 2.4.1 Leads,rate,rhythm,andcardiacaxis ................. 26 2.5 Basicterminologies................................ 33 2.6 HeartDisease ................................... 40 2.6.1 Tachyarrhytmias ............................. 40 2.6.2 Disorders of Impulse Conduction . . . . . . . . . . . . . . . . . . . . 41 2.7 Cardiac Arrest and Sudden Cardiac Death . . . . . . . . . . . . . . . . . . . 42 2.7.1 Long QTSyndrome ........................... 42 2.7.2 L QT - RRRelationship ......................... 45 Chapter 3: Abstract Dynamical System 48 3.1 Introduction.................................... 48 3.2 SymbolicDynamicalSystems .......................... 51 3.3 HyperbolicDynamic ............................... 53 v 3.4 MeasureTheory ................................. 54 3.4.1 Invariant Measures and Measure Preserving Maps . . . . . . . . . . 55 3.4.2 ExistenceofInvariantMeasures..................... 55 3.5 Ergodicity..................................... 56 3.5.1 UniquelyErgodicMaps ......................... 59 3.5.2 MixingDynamicalSystems ....................... 60 3.6 Recurrence .................................... 64 3.6.1 ReturnTimesofDynamicalSystem .................. 65 3.6.2 ReturnTimesStatistics ......................... 66 Chapter 4: Chaos 68 4.1 WhatisChaos?.................................. 68 4.1.1 Lyapunov Exponents . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.1.2 Kolmogorov-Sinai Entropy. . . . . . . . . . . . . . . . . . . . . . . . 74 4.1.3 Kaplan-Yorke(Lyapunov)Dimension ................. 75 4.1.4 Scatter Plot Analysis. . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.2 Stationarity .................................... 77 4.2.1 TestingMethods ............................. 79 Chapter 5: Prelude to Novel Heart Dynamics: Computational Informa- tion Theory of Abstract Dynamical Systems 86 5.1 Introduction.................................... 86 5.2 Kolmogorov-Sinai Axiomatization. . . . . . . . . . . . . . . . . . . . . . . . 88 5.3 Directed Mutual Information . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.4 CanonicalCorrelationAnalysis ......................... 93 5.4.1 LinearCanonicalCorrelationAnalysis ................. 95 5.4.2 NonlinearCanonicalCorrelationAnalysis ............... 97 5.4.3 AggregatedProcesses .......................... 99 5.5 ModelSelection(LagWindow) ......................... 100 5.5.1 AkaikeInformationCriteria....................... 100 Chapter 6: Dynamical System Theory Analysis of Clinical Data 102 6.1 Introduction.................................... 102 6.2 Clinical Data Collection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.3 DynamicalSystemApproachtoHeartPhysiology............... 103 6.3.1 From Time Series to Kolmogorov Abstract Dynamical Systems . . . 104 6.3.2 Kolmogorov-Sinai Axiomatization . . . . . . . . . . . . . . . . . . . 105 6.3.3 PoincaréRecurrence ........................... 106 6.3.4 Chaos ................................... 109 6.3.5 StatisticsofHeartDynamics ...................... 113 6.4 InformationTheoreticAnalyses......................... 122 6.4.1 CanonicalCorrelationAnalysis ..................... 124 Chapter 7: Conclusion 132 vi Bibliography 136 Appendix A: Lyapunov Exponent 143 A.1 Algorithm[89].................................. 143 Appendix B: Hénon Map 147 B.1 StationarityTest ................................. 149 B.2 Chaos-Lyapunov Exponent . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 B.3 Kolmogorov-Sinai Mutual Information . . . . . . . . . . . . . . . . . . . . . 155 vii List of Tables Table2.1:Calculatingthecardiacaxis......................... 32 Table 6.1:Lyapunov exponents, Kolmogorov-Sinai entropy and Kaplan-Yorke di- mension.................................... 113 Table 6.2:Average of mutual information of RR-QT and QT-RR (in bold). . . . 126 viii List of Figures Figure 1.1: A typical computer generated trajectory of heart dynamics. . . . . . 2 Figure 1.2: Hear rate for di fferentsubjects. ..................... 6 Figure 2.1: Transport pathways through the cell membrane, and the basic mech- anismsoftransport. ........................... 13 Figure 2.2: Di ffusion of a fluid molecule during a thousandth of a second. . . . . 14 Figure 2.3: Transport of sodium and potassium ions through protein channels. Also shown are conformational changes in the protein molecules to openorclose"gates"guardingthechannels............... 16 Figure 2.4: E ffect of concentration of a substance on rate of di ffusion through a membrane by simple di ffusion and facilitated di ffusion. This shows that facilitated di ffusion approaches a maximum rate called the V max .16 Figure 2.5: Postulated mechanism for facilitated di ffusion.............. 17 Figure 2.6: Action potential (in millivolts) from a Purkinje fiber of the heart, showinga"plateau." ........................... 22 Figure 2.7: The His-Purkinje conduction system. . . . . . . . . . . . . . . . . . . 28 Figure2.8: TheECGstrip............................... 28 Figure 2.9: Position of the six chest electrodes for standard 12 lead electrocar- diography. V1: right sternal edge, 4th intercostal space; V2: left sternal edge, 4thintercostal space; V3: betweenV2and V4; V4: mid- clavicular line, 5th space; V5: anterior axillary line, horizontally in line with V4; V6: mid-axillary line, horizontally in line with V4. . . 29 Figure2.10:Hexaxial diagram (projection of six leads in vertical plane) showing eachlead’sviewoftheheart. ...................... 32 Figure2.11:ComplexshowingPwavehighlighted. ................. 33 Figure2.12:Normal duration of PR interval is 0.12-0.20 s (three to five small squares)................................... 35 ix Figure2.13:CompositionofQRScomplex....................... 36 Figure2.14:The ST segment should be in the same horizontal plane as the TP segment; the J point is the point of inflection between the S wave and STsegment................................. 38 Figure2.15:ComplexshowingTwavehighlighted. ................. 39 Figure2.16:The QT interval is measured in lead aVL as this lead does not have prominentUwaves(diagramisscaledup)................ 40 Figure3.1: Poincarérecurrencetime. ........................ 65 Figure 4.1: Geometrical definition for Lyapunov exponents. . . . . . . . . . . . . 72 Figure 4.2: A typical scatter plot. . . . . . . . . . . . . . . . . . . . . . . . . . . 77 Figure 5.1: Relationship between Information Theory and the other fields. . . . 87 Figure5.2: Closedloopsystemdiagram. ...................... 89 Figure 5.3: A feedbacked system: Antegrade system A and retrograde system R.91 Figure 6.1: Feedback as a schematic way to exhibit complex relation. . . . . . . 104 Figure 6.2: Probability distribution and best Erlang fit for a patient under drug consumption; solid line: Probability distribution of RR intervals; dot- dash lines: Best Erlang distribution fit.................. 109 Figure 6.3: Left/right fit to the experimental probability distribution of RR in- terval; left tail: Erlang distribution; right tail: Normal distribution. . 110 Figure 6.4: Left/right fit to the experimental probability distribution of RR in- terval; left tail: Erlang distribution; right tail: Weibull distribution. . 110 Figure 6.5: Lyapunov exponents: (a) Baseline, (b) Placebo, (c) Low Dose, (d) MediumDose,(e)HighDose. ...................... 112 Figure 6.6: Poincaré meta recurrence stationarity test for the control subject: (a) Auto prediction error; (b) Cross prediction error. . . . . . . . . . . . 114 Figure 6.7: Poincarémetarecurrencestationaritytestfortheplacebosubject: (a) Auto prediction error; (b) Cross prediction error. . . . . . . . . . . . 114 Figure 6.8: Poincaré meta recurrence stationarity test for subject under low drug doage: (a) Auto prediction error; (b) Cross prediction error. . . . . . 115 Figure 6.9: Poincarémetarecurrencestationaritytestforthesubjectundermedium drug dosage: (a) Auto prediction error; (b) Cross prediction error. . 115 x Figure6.10:Poincaré meta recurrence stationarity test for the subject under high drug dosage: (a) Auto prediction error; (b) Cross prediction error. . 116 Figure6.11:Space-time separation plot stationarity test for the control subject: (a) Embedding dimension m=3; (b) Embedding dimension m=8. . . 117 Figure6.12:Space-time separation plot stationarity test for the placebo subject: (a) Embedding dimension m=3; (b) Embedding dimension m=8. . . 117 Figure6.13:Space-time separation plot stationarity test for the subject under low drug dosage: (a) Embedding dimension m=3; (b) Embedding dimen- sionm=8. ................................. 118 Figure6.14:Space-time separation plot stationarity test for the subject under medium drug dosage: (a) Embedding dimension m=3; (b) Embed- dingdimensionm=8............................ 118 Figure6.15:Space-timeseparationplotstationaritytestforthesubjectunderhigh drug dosage: (a) Embedding dimension m=3; (b) Embedding dimen- sionm=8. ................................. 119 Figure6.16:Scatter plot for the control subject: (a) RR interval; (b) QT − RR intervals................................... 120 Figure6.17:Scatter plot for the placebo subject: (a) RR interval; (b) QT − RR intervals................................... 120 Figure6.18:Scatter plot for subject under low drug dosage: (a) RR interval; (b) QT − RRintervals............................. 121 Figure6.19:Scatter plot for subject under medium drug dosage: (a) RR interval; (b) QT − RRintervals........................... 121 Figure6.20:Scatter plot for subject under high drug dosage: (a) RR interval; (b) QT − RRintervals............................. 122 Figure6.21:Power spectral plots: (a) Control subject; (b) Placebo subject; (c) Subjectunderlowdrugdosage;(d)Subjectundermediumdrugdosage; (e) Subject under high drug dosage. . . . . . . . . . . . . . . . . . . 123 Figure6.22:Mutual Kolmogorov-Sinai information for the control subject: (a) 24- Hourvariation;(b)Hourlyvariation. .................. 127 Figure6.23:Mutual Kolmogorov-Sinai information for the placebo subject: (a) 24-Hourvariation;(b)Hourlyvariation. ................ 128 Figure6.24:Mutual Kolmogorov-Sinai information for the subject under low drug dosage: (a) 24-Hour variation; (b) Hourly variation. . . . . . . . . . 129 xi Figure6.25:Mutual Kolmogorov-Sinai information for the subject under medium drug dosage: (a) 24-Hour variation; (b) Hourly variation. . . . . . . 130 Figure6.26:MutualKolmogorov-Sinaiinformationforthesubjectunderhighdrug dosage: (a) 24-Hour variation; (b) Hourly variation. . . . . . . . . . 131 Figure B.1: Hénon map for a=1 .4,and b=0 .3. .................. 148 Figure B.2: Closed-loop discrete feedback system of Hénon map. . . . . . . . . . 149 Figure B.3: Stationaritytestforthesignal Y inHénonmap: Auto-predictionerror test. .................................... 151 Figure B.4: Stationarity test for the signal Y in Hénon map: Cross-prediction errortest. ................................. 151 Figure B.5: Stationaritytestforthesignal U inHénonmap: Auto-predictionerror test. .................................... 152 Figure B.6: Stationarity test for the signal U in Hénon map: Cross-prediction errortest. ................................. 152 Figure B.7: Space-time separation plot for Hénon map, signal Y .......... 153 Figure B.8: Space-time separation plot for Hénon map, signal U.......... 154 Figure B.9: Lyapunov exponents of Signal Y - Hénon map; embedding dimension m=3.................................... 156 FigureB.10:Lyapunov exponents of signal U - Hénon map; embedding dimension m=3.................................... 156 FigureB.11:Hénonmap: Kolmogorov-Sinaimutualinformationbetweensignals U and Y;segmentedtimeseriesanalysis. ................. 157 FigureB.12:Hénonmap: Kolmogorov-Sinaimutualinformationbetweensignals U and Y;fulltimeseriesanalysis. ..................... 157 xii Abstract Heart dynamics and heart rate variability, as well as their characterization, have been of interest to researchers for years. Sudden death due to heart failure, which could be the result of heart diseases, drug consumptions, or even emotional shocks, is one of the most important issues in public health. Henceforth, prediction and/or awareness of such conditions has always been a top priority. In the present work, a time domain approach, as opposedtothetraditionalfrequency-domainapproach,toheartratedynamicsisintroduced. The Electro-Cardio-Gram (ECG) signals, more specificallytheRRandQTintervalsas characteristics of heart dynamics, are investigated. It is hypothesized that the Erlang distribution of the RR interval is a manifestation of the Poincaré return time phenomenon, and based on clinical results available to us, it is conjectured that for the healthy heart, it is necessary but not su fficient that the RR interval would be the k-fold Poincaré return time. The higher frequency, k times that of normal heartbeat, is hypothesized to be related to the synchronization of the array of pacemaker cells in the sinoatrial (SA) node. In this work, chaoticbehaviorofheart dynamicsisobservedbymeansofsuchnonlineartechniques such as calculation of Lyapunov exponents, Kolmogorov-Sinai (KS) entropy, and Kaplan- Yorke (K-Y) dimensions. These techniques measure chaos, the amount of uncertainty and xiii complexity of the system, respectively. The existence of positive Lyapunov exponent and non-integerK-Ydimensionintheheartrateconfirmsitschaoticdynamics. Heartdynamics is also analyzed via statistical techniques such as linear and nonlinear canonical correlation analyses and mutual KS information. These techniques are implemented to unravel the at point subtle relationship between RR and QT intervals. It is conjectured that past RR intervals are nonlinear regression variables on the future QT intervals for the healthy heart, and in case of any abnormalities, such as prolongation of QT interval, the cause to e ffect relationship from past RR to future QT might be compromised. 1 Chapter 1 Introduction 1.1 The RR Interval In heart physiology, the sinus node (also called sinoatrial node)(SA)isknown as the dominant cardiac pacemaker, which is under the influence of the vagal nerve. The SA node cells have the capability of self-excitation and they behave like a large population of electrically coupled oscillators (pacemakers) with di ffering intrinsic frequencies [64], [85]. Hence, nosinglecell inthe SAnodeserves as pacemaker, andby Michaels et al. hypothesis [64], the pacing rate is not that of the intrinsically fastest pacemaker, but it may actually be from a "democratic consensus" dynamic process of contribution of all cells to establish the rhythm. The SA node depolarizations are transmitted to both the AV node and the Purkinje fibers. The new impulse from the SA node discharges both the AV node and Purkinje fibers before their own threshold for self-excitation can occur in either of these [28]. Thus, the SA node actually controls the heart beats and it is virtually always the pacemaker of the normal heart. 2 Figure 1.1: A typical computer generated trajectory of heart dynamics. Heartratevariability(HRV)referstothe beat-to-beat alternation in heart rate (HR) [57]. Heart rate and heart rate variability are important dynamic measures of the state of the cardiovascular system and the autonomic nervous system (ANS). Heart rate is perhapsthemostaccessiblecardiovascularsignalforanalysisofvariability. Itistraditionally estimatedastheaverageofthereciprocalofthe RR intervalwithinaspecifiedtimewindow [6]. The RR interval is defined as the time di fference between the two peaks (beats) in Electro-Cardo-Gram (ECG) signals. Figure 1.1 shows typical trajectory generated by the dynamical model in [63] in the 3-D space given by (x, y, z). The dashed line reflects the limit cycle of unit radius while the small circles show the positions of the P, Q, R, S,and T events. A healthy heart of a young subject has a large HRV, while decreased or absence of variability may indicate cardiac disease or advanced age. Several investigators have demonstrated that HRV is essentially stable over time despite significant inter-individual HRV variation. This is also true in patients who are on drug therapy, which suggests that 3 analysisofHRVcanbeusedfortheassessmentofthee ffectofpharmacologicalinterventions on cardiac autonomic activity. Another major reason for the interest in measuring HRV stemsfromitsabilitytopredictsurvivalafterheartattack. Overhalfofadozenprospective studies have shown that reduced HRV predicts sudden death in some patients [57]. Sudden cardiac death (also called sudden arrest) is death resulting from an abrupt loss of heart function (cardiac arrest). The victim may or may not have a diagnosed heart disease. Hence, the time and the mode of the death are unexpected; death occurs within minutes after symptoms appear, leaving almost no time for emergency intervention. Abnormalities of conduction, drug e ffects and miscellaneous disorders including noncardiac diseases can also cause dramatic changes in heart dynamics. For example, medications such as digitalis, antidepressants, and antiarrhythmic agents can profoundly alter the heart rhythm observed from ECG signals. Until recently, the time series methodologies applied in the analysis of HRV have mainly been standard frequency domain and time domain techniques. During the past decades methods mainly in statistical physics, including chaotic dynamics, have been used to study the human heart rate. Historically, the idea of the existence of “chaotic” behavior in biological systems in general, and in HRV in particular, is not a new one [21], [49]; whether chaos might exist in the pattern of HRV has been the subject of several investigations. However, the study of chaos in dynamical physiological systems is di fficult because of limitations in data length andintegrityandthecomplexityofmathematicalmodelstothesedata[22]. Therehasbeen increasing evidence to support the case that chaos plays a positive role in the physiology of 4 the organism. It is widely expected that the presence of chaos should be healthy because it providestheorganismwithan“information-rich(broadband)state”and“spectralreserve.” Most of the research seems to indicate that chaotic behavior is a necessary ingredient of normal functioning. Healthy heart responds easily to external signals, and is able to adapt to a wide range of beating rates. This responsiveness gives rise to the high variability of the heart rate. Conversely, high heart rate variability promotes responsiveness to external stimuli. On the other hand, severe heart diseases decrease the responsiveness of the heart withrespecttothewholespectrumofsignalsarrivingfromtheautonomousnervoussystem; this leads to the loss of apparent complexity of the HRV signal. According to classical concepts of physiologic control, healthy systems are self-regulated to reduce variability and maintain physiologic constancy. However, the output of a wide variety of systems, such as thenormalhumanheartbeat,fluctuatesinacomplexmanner,evenunderrestingconditions [21], [88], [101]. This contributes tophysiology beingan inexhaustible source of new control paradigms. Dr. Chi-Sang Poon, principal research scientist in the Harvard-MIT Division of Health Sciences and Technology (HST), has been investigating the seemingly indistinguish- ablerhythmsof healthyanddiseasedheartstoidentifypeopleatriskforcardiacdisease. In his most recent findings, he demonstrates that subtle heart-rate variabilities are not simply random fluctuations, but are complex yet deterministic patterns that are probably chaotic. In earlier studies published in 1997, Poon and C.K. Merrill discovered that heart disease actually makes the heart beat less chaotically. Paradoxically, that means an orderly and regular heartbeat might indicate heart disease. Poon confirmed his previous finding that a 5 weakly chaotic signal might be more susceptible to noise than a strongly chaotic signal and that it would be more likely to indicate congestive heart failure. Poon and Merrill proposed theuseoftheseanalyticalmethodsinnonlineardynamicstodetectanddiagnosecongestive heart failure. Poon also pointed out that when one fails to detect chaos, or even nonlinear dynamics, inexperimentaldata, itdoesnotnecessarilymeanthatthesystemisnotchaotic. The background noise might be so high that chaos in the signal is totally covered by the noise. The chaotic behavior of the heart rate variability can be explained in terms of dynamics of the complex of sinoatrial (SA) and atrioventricular (AV)nodes,whichthe former has been modeled as a system of non-linear coupled oscillators, responsible for the heart rhythm [88]. In the case of human heart, the “noise” comes from the autonomous nervoussystemintheformofvagus nerveinputs regulatingtheheartrate. Byobservations itisunderstoodthattheheartrateiscontrolledbynon-deterministicande ffectivelyrandom signals arriving from the autonomous nervous system [88]. Further support of chaos being inherent to the physiological control of heart rate is provided by Narayanan et al. [67]. They discovered unstable periodic orbits (UPO) in RR timeseries. Periodicorbits reflect thedeterministic dynamics ofthe underlyingsystem, but the UPOs are unstable, which means that they have a positive Lyapunov exponent and are therefore typical for chaotic systems [101]. Figure (1.2) shows the heart rate for di fferent subjects. Parts A and C in Figure (1.2) are from patients with severe congestive heart failure; part B is an illustration of a healthyheartrate; andpartDisfromapatientwithaconditionknownasatrialfibrillation, 6 Figure 1.2: Hear rate for di fferent subjects. in which the heart beats erratically [101]. Loss of complexity in the heartbeat is a result of illness or aging. Hashida and Tasaki in 1984 [32] have investigated the nature of the RR interval during the atrial fibrillation in man and proposed a function for the atrioventricular ( AV) node. Theyhaveclaimedthatthepatternofthe RRintervalinatrialfibrillationisunimodal andskewedtotherightandsocanbe fitted to an Erlang distribution. This conclusion was based on the assumption that the atrial inputs can be considered to be summated to a certain threshold for ventricular activation in the AV node; it was assumed that inputs to 7 the AV node from the atria occur at times determined by a Poisson process and the AV node firing is due to a summed e ffect of k inputs, and as soon as it reaches a threshold. Hence the function of the AV node during the atrial fibrillation would be to transform an exponentially distributed input into an Erlang distributed output. This is the issue to be discussedinthisworkinmoredetailslater,asitis flawed because of the failure to take into consideration the e ffect of digitalis on which the patients had been put. It has been also claimed that the sequence of the RR intervals are independent and hence belongs to a renewal process (last paragraph on page 672 in ref. [32]). "... Since the stationarity and independency of the sequence of RR intervals in atrial fibrillation can be proved as shown below, this sequence belongs to a renewalprocess...." Also on page 679 of the same paper it is written that, "... Atrial inputs arrive at the AV node in accord with a Poisson processes. Thispossibilityisnotquitesoartificialasitmayappearbecauseoftheexistence of a general theorem in the theory of queueing which states that, under fairly general conditions, the superposition of several independent series of simpler forms is indistinguishable from a Poisson process ...." It is well-known that the superposition of a finite number of renewal process is not renewal anymore. Under fairly general conditions the distribution of the waiting time is approximately exponential and so the combined process is close to Poisson process. Based on the Limit theorem this is true asymptotically when infinitely large number of renewal processes are superposed. Also from the physiology point of view, the mechanism of the AV node is as a bridge or passage for the SA node conduction and there is no threshold potential in AV node. In this work, we propose a new hypothesis regarding the dynamics of the Sinoa- trial (SA) node, which itself generates the RR intervals. We proceed from the reasonable 8 hypothesis that the SA node is a dynamical attractor. An earlier hypothesis was that a heartbeat is emitted when some abstract state of the SA attractor enters some finite region A, so that the RRintervalwouldhavebeenthePoincaréreturntimeto A. Unfortunately, the theoretically exponential distribution of the Poincaré return time is not a fitofthe experimental RR histograms. However, the theoretical Erlang distribution of the k-fold return time is a surprisingly tight fitoftheexperimental RR histograms, of some but not all research subjects. 1.2 The QT versus The RR Interval As it will appear clear soon, it is necessary to develop a complex system para- digm specific to processes that can be observed through multiple sensing ports, but the data collected on such ports leaves the experimenter with a bewildering array of variables (e.g., RR, QT, QRS, ST level,...), with the di fficulty to determine which variables could be considered as causes and which one should be considered as e ffects. A good example is provided by weather phenomena: the geostrophic wind is usually considered to be a conse- quence of the atmospheric pressure, although on frontal systems the causality is reversed, as it is the wind that is the cause of the pressure drop ahead of the cold front. Electrical circuits have long been known to embody the same puzzle, as in the impedance formulation the current is the input and the voltage is the output, while in the admittance formulation, the voltage is the input and the current is the output. This puzzle persists under electro- mechanical analogy, as in robotics impedance control deals with the issue of whether the force is the input and the velocity the output, or the other way round, or some kind of 9 halfway compromise. It is anticipated that this puzzle will soon pervade neurophysiology. Indeed, it is usually admitted that nerve transmission is electrochemical and regulated by the Hodgkin-Huxley equations. However, recently, Heimburg and Jackson [42] formulated evidence that nerve transmission is through acoustic wave in the fiber membrane [90, 15]. Their (still disputed) argument is based on the fact that chemically di fferent anesthetics share the common physical property of lowering the temperature of the chain melting tran- sition of the lipid membrane so that the shock wave that normally brings the membrane 85% through chain melting no longer happens under anesthetics [44]. Whether nerve trans- mission is electrochemical or acoustic might not be the relevant issue, as the two aspects might coexist, with the real issue being the causality relation between the acoustic and the electrochemical phenomena. In fact, more recently, Heimburg and Jackson [43, Section 4.4] have argued that the shock wave is the cause while the electric activity is the e ffect. Petrov [73] holds the same stand. Tasaki [99] on the other hand has the more traditional standthatthemechanicaldeformationistheeffectoftheactionpotential. Controlprovides one motivation for the need to determine which variables are causes and which variables are e ffects, as the variables that are causes have enough control authority over the e ffect variablesandshouldhencebechosenascontrolinputs. Butotherapplicationsincludeover- all system monitoring, especially for systems that can switch among di fferent attractors, as some bifurcations are accompanied by causality flips. Probably the most popular approach thus far to determine which variables are “inputs” and which variables are “outputs” is Willems’ behavioral theory [102]. However, the latter purely algebraic approach applies in the realm of linear processes subject to 10 perfect measurements. Here we are rather dealing with nonlinear, event driven phenomena, observable through corrupted measurements, so that, instead of an algebraic approach, a “softer,” information theoretic one is called for. More specifically, this work develops a heart dynamics case study of the above concepts. Whilethereareagreatmanyvariablesthatcanbeassociatedwiththecardiovas- cular system, here we retain two parameters: the RR and QT sequences. The RR interval is the amount of time between consecutive heart beats. The QT intervalistheamountof time between the beginning of the Qwaveandtheendofthe T wave. In a few words, the QT intervalistheamountoftimetheventriclestaketodepolarizeandrepolarizeduringthe heart cycle. Under normal conditions, the RR sequence is the pacemaker and should hence control the QT interval. However, under some pathological condition typically created by Sodium channels blockers, the QT interval is prolonged, which could lead to ventricular fibrillation, a life threatening situation. Under prolonged QT interval, the RR interval no longer completely controls the QT interval; in fact, a causality is compromised: the QT interval could sporadically become the cause and the RR the e ffect. Arelatedapproachwas devisedin[74] forbrainwaves, thesizabledi fferencebeing that the latter does not relate causality to feedback. The terminology of antegrade/retrograde system is borrowed from cardiology [53, 1], where antegrade refers to the normal conduction from the atria to the ventricles and retrograde refers to conduction from the ventricles to the atria. The latter is sometimes referred to as reentry [13]. 11 1.3 Summary The main focus of the present work is on the measure theoretic and statistical properties of dynamical systems. As a summary, we begin with the Poincaré recurrence and Ergodic theorems, and hypothesize that the RR interval is a candidate for the return time of the heart dynamical system. We also investigate chaos in heart dynamics. Then we proceed to study the statistical properties of a heart dynamics by looking at the correlation functions and their long term behavior. An outline of the rest of this work is as follows: Section 2 is an introduction to Heart physiology and contains some general facts about HRV. Sections 3 to 5 contain the mathematical background and the methods used to analyze the HRV time series, such as Lyapunov Exponents, Kolmogorov-Sinai entropy, and linear/nonlinear Canonical Correla- tion analyses. Section 6 includes discussions on the analyses of clinical data and simulation results; conclusion and future studies are followed thereafter in Section 7. The algorithms and examples are presented in Appendices. 12 Chapter 2 Heart Physiology 2.1 Introduction In this chapter, we briefly review the heart physiology and some fatal heart con- ditions [28]. 2.2 Dynamics of Excitable Cells Cell physiology and structure has been extensively studied during past decades. Transmembranepotential and electrical conduction of cells have been investigated by many scientists such as Young in 1936, Cole, Curtis, Marmont, Hodgkin and Huxley in 1940s and 1950s (see ref. in Cole 1968, Hodgkin 1964). Hodgkin and Huxley were awarded the Nobel Prize in Physiology or Medicine in 1963 for their work on modeling the ionic basis of action potential. 13 Figure 2.1: Transport pathways through the cell membrane, and the basic mechanisms of transport. The life of the cells is extremely dependent on the substances in the intercellular rand extracellular fluid and their di fferences. Of the most important substances are sodium and potassium; for example, the extracellular fluid contains a large amount of sodium but only asmallamountof potassium, whereas, exactlytheoppositeistruefor the intracellular fluid. Also, the extracellular fluid contains a large amount of chloride ions, whereas the intracellular fluid contains very little. The membrane consists almost entirely of a lipid bilayer, with large numbers of protein molecules in the lipid, many of which penetrate all the way through the membrane, asshowninFigure2.1. Thelipidbilayerisabarrieragainsttheextracellularandintracellu- lar fluid compartments; however, a few substances, mainly lipid-soluble ones, can penetrate it by di ffusing directly through the lipid substance itself. Theproteinmoleculesinthemembranehaveentirelydi fferentpropertiesfortrans- porting substances. Di fferent proteins function di fferently. Two main of these proteins are channel proteins and carrier proteins. Both the channel proteins and the carrier proteins are usually highly selective in the types of molecules or ions that are allowed to cross the membrane. 14 Figure 2.2: Di ffusion of a fluid molecule during a thousandth of a second. Transport through the cell membrane, either directly through the lipid bilayer or through the proteins, occurs by one of two basic processes: di ffusion or active transport. Briefly speaking, di ffusion means continual and random molecular movement of substances molecule by molecule, either through intermolecular spaces in the membrane or in combi- nation with a carrier protein. By contrast, active transport is the movement of ions or other substances across themembraneincombinationwithacarrierproteinthatcausesthesubstancetomovefrom a low-concentration state to a high-concentration state. Di ffusionthroughthecellmembraneisdividedintotwosubtypescalledsimple dif- fusion and facilitated di ffusion.Simpledi ffusion means that kinetic movement of molecules or ions occurs through a membrane opening or through intermolecular spaces without any interaction with carrier proteins in the membrane. The rate of di ffusion is determined by the amount of substance available, the velocity of kinetic motion, and the number and sizes of openings in the membrane through which the molecules or ions can move. Facilitated di ffusion requires interaction of a carrier protein. The carrier protein aids passage of the 15 molecules or ions through the membrane by binding chemically with them and shuttling them through the membrane in this form. Many of the protein channels are highly selective for transport of one or more specific ions or molecules. This results from the characteristics of the channel itself, such as its diameter, its shape, and the nature of the electrical charges and chemical bonds along its inside surfaces. As an example, one of the most important of the protein channels, the so-called sodium channel, specifically selective for passage of sodium ions, is only 0.3 by 0.5 nanometer in diameter, but more important, the inner surfaces of this channel are strongly negatively charged, as shown by the negative signs inside the channel proteins in the top panel of Figure 2.3. These strong negative charges can pull small dehydrated sodiumions intothesechannels, actuallypullingthesodiumionsawayfromtheirhydrating water molecules. Once in the channel, the sodium ions di ffuse in either direction according to the usual laws of di ffusion. The same mechanism is for the so called potassium channel, specifically for passage of potassium. Althoughtherateofsimpledi ffusion through an open channel increases propor- tionately with the concentration of the di ffusing substance, in facilitated di ffusion the rate of di ffusion approaches a maximum, called Vmax, as the concentration of the di ffusing substance increases. This di fference between simple di ffusion and facilitated di ffusion is demonstrated in Figure 2.4. The figure shows that as the concentration of the di ffusing substance increases, the rate of simple di ffusion continues to increase proportionately, but in the case of facilitated di ffusion, the rate of di ffusion cannot rise greater than the Vmax level. 16 Figure 2.3: Transport of sodium and potassium ions through protein channels. Also shown are conformational changes in the protein molecules to open or close "gates" guarding the channels. Figure 2.4: E ffect of concentration of a substance on rate of di ffusion through a membrane bysimpledi ffusionandfacilitateddi ffusion. Thisshowsthatfacilitateddi ffusionapproaches a maximum rate called the V max . 17 Figure 2.5: Postulated mechanism for facilitated di ffusion. The reason is the mechanism illustrated in Figure 2.5. Specific molecules are transported through the carrier protein by means of a binding "receptor" on the inside of the protein. The binding receptor interacts with the molecules by a weak binding force that the thermal motion of the attached molecule causes it to break away and the molecule is released on the opposite side of the membrane. The rate at which molecules can be transported by this mechanism can never be greater than the rate at which the carrier protein molecule can undergo change back and forth between its two states. The molecules can diffuse in both ways in this mechanism. The opening and closing of gates are controlled in two principal ways: Voltage gating and Chemical (ligand) gating [28]. Perhaps the most direct evidence for the exis- tence of single channels in the membranes of cells comes from the patch-clamp technique, introducedbyNeherandSakmannin1991[28]. Therewillbeonlyonechannelinthepatch of membrane subtended by the rim of the electrode, and one will pick up a signal, showing 18 thatthechannelopensandclosesinanapparentlyrandomfashion, allowingafixedamount of current (on the order of picoamperes) to flow through it when it is in the open state [27]. 2.3 Action Potential Electrical potentials exist across the membranes of virtually all cells of the body. In addition, some cells, such as nerve and muscle cells, are capable of generating rapidly changing electrochemical impulses at their membranes, and these impulses are used to transmitsignalsalongthenerveormusclemembranes. Wediscussthemembranepotentials generated both at rest and during action by nerve and muscle cells. Di ffusion Potential caused by an ion concentration di fference on the two sides of the membrane. The rate at which the substance di ffuses inward is proportional to the concentration of molecules on the outside, because this concentration determines how many molecules strike the outside of the membrane each second. Conversely, the rate at whichmoleculesdi ffuseoutwardisproportionaltotheirconcentrationinsidethemembrane. Therefore, the rate of net di ffusion into the cell is proportional to the concentration on the outside minus the concentration on the inside, i.e., C o − C i ,where C o is concentration outside and C i is concentration inside. The Nernst potential for an ion is defined as the exact opposite of the di ffusion potentiallevelacrossamembrane. Itdeterminesatwhatpointthetwoforcesfromchemical and electrical gradients balance each other and an ionic equilibrium happens. The magni- tude of this Nernst potential (also called equilibrium potential) is determined by the ratio of the concentrations of that specific ion on the two sides of the membrane. The greater 19 this ratio, the greater the tendency for the ion to di ffuse in one direction, and therefore the greater the Nernst potential required to prevent additional net di ffusion. The following equation, called the Nernst equation, can be used to calculate the Nernst potential for any univalent ion is: E = V in − V out = RT zF log( Concentration Outside = C o Concentration inside = C i ) (2.1) where E is electromotive force (in millivolts), R is Rydberg gas constant, T is the tempera- ture in degrees Kelvin, z is the charge on the ion, and F is Faraday’s constant. When using this formula, it is usually assumed that the potential in the extracellular fluid outside the membrane remains at zero potential, and the Nernst potential is the potential inside the membrane. Also, the sign of the potential is positive (+) if the ion di ffusing from inside to outside is a negative ion, and it is negative (-) if the ion is positive [28]. For a body temperature (2.1) is simplified to E =±61log( C 0 C i ). Example 1 Nernst Potential for Potassium. Assume the extracellular concentration of potassium [ K + ] o =4(mM) and its intracellular concentration [ K + ] i = 120(mM).If the cell were to become freely permeable to potassium (only), what would be the equilibrium potential. Solution 2 By (2.1), in a body temperature, the equilibrium potential for potassium di ffu- sion across the cell membrane is E=61log( [ K + ] o [ K + ] i )=61log( [ K + ] o [ K + ] i )= −90 .1(mV ol t) As we expected the equilibrium potential is negative since potassium is a positive ion and di ffusingfrominsidetooutside. 20 The principal ions involved in an action potential are sodium and potassium ions; sodium ions enter the cell, and potassium ions leave, restoring equilibrium. Relatively few ions need to cross the membrane for the membrane voltage to change drastically. The ions exchanged during an action potential, therefore, make a negligible change in the interior and exterior ionic concentrations. The few ions that do cross are pumped out again by the continual action of the sodium—potassium pump, which, with other ion transporters, maintains the normal ratio of ion concentrations across the membrane. During the early portion of the action potential, the ratio of sodium to potassium conductance increases more than 1000-fold. Therefore, far more sodium ions flow to the interior of the fiber than do potassium ions to the exterior. This is what causes the mem- brane potential to become positive at the action potential onset. Then the sodiumchannels begin to close and the potassium channels to open, so that the ratio of conductance shifts far in favor of high potassium conductance but low sodium conductance. This allows very rapid loss of potassium ions to the exterior but virtually zero flow of sodium ions to the interior. Consequently, the action potential quickly returns to its baseline level. The membranes of almost all cells of the body have a calcium pump similar to the sodium pump, and calcium serves along with (or instead of) sodium in some cells to cause most of the action potential. The calcium pump acts similar to the sodium pumps. The calcium channels are slow to become activated, requiring 10 to 20 times as long for activation as the sodium channels. Therefore, they are called slow channels, in contrast to the sodium channels, which are called fast channels. 21 2.3.1 Threshold for Initiation of the Action Potential Enough initial rise in the membrane potential from -90 millivolts toward the zero level, the rising voltage itself causes many voltage-gated sodium channels to begin opening. This allows rapid inflow of sodium ions, which causes a further rise in the membrane po- tential, thus openingstill more voltage-gatedsodiumchannels andallowingmorestreaming of sodium ions to the interior of the fiber. This process is a positive-feedback vicious cy- cle that, once the feedback is strong enough, continues until all the voltage-gated sodium channels have become activated (opened). Then, within another fraction of a millisecond, the rising membrane potential causes closure of the sodium channels as well as opening of potassium channels, and the action potential soon terminates. An action potential will not occur until the initial rise in membrane potential is great enough to create the vicious cycle described in the preceding paragraph. This occurs when the number of Na + ions entering the fiber becomes greater than the number of K + ions leaving the fiber. A sudden rise in membrane potential of 15 to 30 millivolts usually is required. Therefore, a sudden increase in the membrane potential in a large nerve fiber from -90 millivolts up to about -65 millivolts usually causes the explosive development of an action potential. This level of -65 millivolts is said to be the threshold for stimulation. An action potential elicited at any one point on an excitable membrane usually excites adjacent portions of the membrane, resulting in propagation of the action potential along the membrane. The action potential travels in all directions away from the stimulus- even along all branches of a nerve fiber-until the entire membrane has become depolarized. 22 Figure 2.6: Action potential (in millivolts) from a Purkinje fiber of the heart, showing a "plateau." 2.3.2 Plateau in Some Action Potentials In some instances, the excited membrane does not repolarize immediately after depolarization; instead, the potential remains on a plateau near the peak of the spike potential for many milliseconds, and only then does repolarization begin. Such a plateau is shown in Figure 5-13; one can readily see that the plateau greatly prolongs the period of depolarization. This type of action potential occurs in heart muscle fibers, where the plateau lasts for as long as 0.2 to 0.3 second and causes contraction of heart muscle to last for this same long period. Thecauseoftheplateauisacombinationofseveralfactors. First,inheartmuscle, two types of channels enter into the depolarization process: (1) the usual voltage-activated sodium channels, called fast channels, and (2) voltage-activated calcium-sodium channels, which are slow to open and therefore are called slow channels. Opening of fast channels causes the spike portion of the action potential, whereas the slow, prolonged opening of the slowcalcium-sodiumchannelsmainly allowscalciumions to enter thefiber, whichis largely responsible for the plateau portion of the action potential as well. 23 A second factor that may be partly responsible for the plateau is that the voltage- gated potassium channels are slower than usual to open, often not opening very much until the end of the plateau. This delays the return of the membrane potential toward its normal negative value of -80 to -90 millivolts. 2.3.3 The Hodgkin-Huxley Equations The Hodgkin—Huxley model is a scientific model that describes how action poten- tials in neurons are initiated and propagated. It is a set of nonlinear ordinary di fferential equations that approximates the electrical characteristics of excitable cells such as neurons and cardiac myocytes. Alan Lloyd Hodgkin and Andrew Huxley described the model in 1952 to explain the ionic mechanisms underlying the initiation and propagation of action potentials in the squid giant axon [46]. They received the 1963 Nobel Prize in Physiology or Medicine for this work. The Hodgkin-Huxley equation can be expressed as a system of four coupled non- linear ordinary di fferential equations, dV dt = − 1 C £ ¯ g Na m 3 h( V − E Na )+¯ g K n 4 ( V − E K )+¯ g L ( V − E L )+ I stim ¤ , dm dt = α m (1 − m) − β m m, dh dt = α h (1 − h) − β h h, (2.2) dn dt = α n (1 − n) − β n n, where ¯ g ∗ is the maximal conductance, E ∗ is the equilibrium potential, I stim is the total stimulus current, n is a "gating" variable, m is the activation variable, h is inactivation variable, and α ∗ and β ∗ are voltage-dependent rate constants. I stim might be a periodic 24 pulse train or a constant ("bias") current. n is actually controlling the opening and closing of potassium channel. The rate constants, α ∗ and β ∗ control the transitions between the closed and open states of the gate. Since the full Hodgkin-Huxley equations are four-dimensional system of ordinary di fferential equations, it is di fficult to obtain a visual picture of trajectories in this system. In Hodgkin-Huxley equations, variables m and V tracked each other during the action potential, so that one could be expressed as an algebraic function of the other (this is also hold true for h and n). FitzHugh and Nagumo et al. could simplify Hodgkin-Huxley equations to a simple two dimensional equation at the same time. dx dt = c( x − x 3 3 + y + S( t)) , dy dt = − ( x − a+ by) c , where x is a variable representing transmembrane potential and excitability, y is variable for refractoriness and accommodation. The function S( t) is the stimulus, and a, b,and c are parameters. Although the Hodgkin-Huxley model has been a great success, there are some dis- crepancies betweenexperimentandmodel. Thishasled toupdatedversionoftheHodgkin- Huxley model being produced to account for these discrepancies [27]. 2.4 Electricity and the Heart Electricity, an innate biological electricity, is what makes the heart go. The ECG is a recorder of the heart’s electrical activity; through the perturbations in the normal electrical patterns many di fferent cardiac disorders are diagnosed. Cardiac cells are electri- 25 cally polarized, that is, their insides are negatively charged with respect to their outsides by means of membrane pumps that ensure the appropriate distribution of ions (primar- ily potassium, sodium, chloride, and calcium) necessary to keep the insides of these cells relatively electronegative [100]. cardiac cells are losing their internal negativity in depolar- ization, which is a fundamental electrical event of the heart. Depolarization is propagated from cell to cell, producing a wave of depolarization that can be transmitted across the entire heart. This wave represents a flow of electricity, an electrical current that can be detected by electrodes placed on the surface of the body. After depolarization is complete, thecardiaccellsareabletorestoretheirrestingpolaritythroughaprocessof repolarization. This, too, can be sensed by recording electrodes. From the standpoint of the electrocardiographer, the heart consists of three types of cells [100]: • Pacemaker cells — the electrical power source of the heart • Electrical conducting cells — the hard wiring of the heart • Myocardial cells — the contractile machinery of the heart. sinoatrial(SA)node: Thegroupofcells,dominantpacemakercells,locatedhigh up in the right atrium. This is alsocalledsinus nodefor short. Thesecellstypicallyfire at a rate of 60 to 100 times per minute, varying tremendously depending on the activity of the autonomic nervous system (sympathetic stimulation from adrenaline accelerates the sinus node,whereasvagalstimulationslowsit)andthedemandsofthebodyforincreasedcardiac output. In fact, everycellsinthehearthastheabilitytobehavelikea pacemaker cell. This 26 so-called automatic ability is normally suppressed unless the dominant cells of the sinus node fail or if something in the internal or external environment of a cell (sympathetic stimulation, disease, etc.) stimulates its automatic behavior. Atrioventricular (AV)node: The AV nodeislocatedclosetothewallbetween right and left atriums. The electrical conduction from SA node must be funneled along the interventricular septum, the wall that separates the right and left ventricles. Here, the AV node slows conduction to a crawl. This pause lasts only a fraction of a second. This physiologic delay in conduction is essential to allow the atria to finish contracting before the ventricles begin to contract. Like sinus node, the AV node is also under the influence of the autonomic system [100]. The waves that appear on an ECG primarily reflect the electrical activity of the myocardial cells, which compose the vast bulk of the heart. 2.4.1 Leads, rate, rhythm, and cardiac axis Electrocardiography is a fundamental part of cardiovascular assessment. It is an essential tool for investigating cardiac arrhythmias and is also useful in diagnosing cardiac disorders such as myocardial infarction. Familiarity with the wide range of patterns seen in theelectrocardiogramsofnormalsubjectsandanunderstandingofthee ffectsofnon-cardiac disorders on the trace are prerequisites to accurate interpretation. The contraction and relaxation of cardiac muscle results from the depolarization and repolarization of myocardial cells. These electrical changes are recorded via electrodes placed on the limbs and chest wall and are transcribed on to graph paper to produce an electrocardiogram (commonly known as an ECG). The sinoatrial node acts as a natural 27 pacemaker and initiates atrial depolarization. The impulse is propagated to the ventricles by the atrioventricular node and spreads in a coordinated fashion throughout the ventricles via the specialized conducting tissue of the His-Purkinje system. Thus, after delay in the atrioventricular mode, atrial contraction is followed by rapid and coordinated contraction of the ventricles. The word "Lead" sometimes causes confusion. Sometimes it is used to mean the piecesofwirethatconnectthepatienttotheECGrecorder. Properly, aleadisanelectrical picture of the heart. The electrical signal from the heart is detected at the surface of the body through five electrodes, which are joined to the ECG recorder by wires [31]. The electrocardiogram is recorded on to standard paper travelling at a rate of 25 mm/s. The paper is divided into large squares, each measuring 5 mm wide and equivalent to 0.2s. Each large square is five small squares in width, and each small square is 1 mm wide and equivalent to 0.04 s. The electrical activity detected by the electrocardiogram machine is measured in millivolts. Machines are calibrated so that a signal with an amplitude of 1 mV moves the recording stylus vertically 1 cm. Throughout this text, the amplitude of waveforms will be expressed as: 0.1 mV = 1 mm = 1 small square. The amplitude of the waveform recorded inanyleadmaybeinfluenced by the myocardial mass, the net vector of depolarization, the thickness and properties of the intervening tissues, and the distance between the elec- trode and the myocardium. Patients with ventricular hypertrophy have a relatively large myocardial mass and are therefore likely to have high amplitude waveforms. In the pres- 28 Figure 2.7: The His-Purkinje conduction system. Figure 2.8: The ECG strip. 29 Figure 2.9: Position of the six chest electrodes for standard 12 lead electrocardiography. V1: right sternal edge, 4th intercostal space; V2: left sternal edge, 4th intercostal space; V3: between V2 and V4; V4: mid-clavicular line, 5th space; V5: anterior axillary line, horizontally in line with V4; V6: mid-axillary line, horizontally in line with V4. ence of pericardial fluid, pulmonary emphysema, or obesity, there is increased resistance to current flow, and thus waveform amplitude is reduced. The direction of the deflection on the electrocardiogram depends on whether the electrical impulse is travelling towards or away from a detecting electrode. By convention, an electrical impulse travelling directly towards the electrode produces an upright (“positive”) deflection relative to the isoelectric baseline,whereas an impulse moving directly away from an electrode produces a downward (“negative”) deflection relative to the baseline. When the wave of depolarization is at right angles to the lead, an equiphasic deflection is produced. The six chest leads (V1 to V6) “view” the heart in the horizontal plane. The information from the limb electrodes is com- bined to produce the six limb leads (I, II, III, aVR, aVL, and aVF), which view the heart in the vertical plane. The information from these 12 leads is combined to form a standard electrocardiogram. 30 The arrangement of the leads produces the following anatomical relationships: leads II, III, and aVF view the inferior surface of the heart; leads V1 to V4 view the anterior surface; leads I, aVL, V5, and V6 view the lateral surface; and leads V1 and aVR look through the right atrium directly into the cavity of the left ventricle. Rate The term tachycardia is used to describe a heart rate greater than 100 beats/min. A bradycardia is defined as a rate less than 60 beats/min (or < 50 beats/min during sleep). One large square of recording paper is equivalent to 0.2 seconds; there are five large squares per second and 300 per minute. Thus when the rhythm is regular and the paper speed is running at the standard rate of 25 mm/s, the heart rate can be calculated by counting the number of large squares between two consecutive R waves, and dividing this number into 300. Alternatively, the number of small squares between two consecutive R waves may be divided into 1500. Some countries use a paper speed of 50 mm/s as standard; the heart rate is calcu- lated by dividing the number of large squares between R waves into 600, or the number of small squares into 3000. “Rate rulers” are sometimes used to calculate heart rate; these are used to measure two or three consecutive R-R intervals, of which the average is expressed as the rate equivalent. When using a rate ruler, take care to use the correct scale according to paper speed (25 or 50 mm/s); count the correct numbers of beats (for example, two or three); and restrict the technique to regular rhythms. When an irregular rhythm is present, the heart rate may be calculated from the rhythmstrip(seenextsection). Ittakesonesecondtorecord2.5cmoftrace. Theheartrate 31 per minute can be calculated by counting the number of intervals between QRS complexes in 10 seconds (namely, 25 cm of recording paper) and multiplying by six. Rhythm To assess the cardiac rhythm accurately, a prolonged recording from one lead is used to provide a rhythm strip. Lead II, which usually gives a good view of the P wave, is most commonly used to record the rhythm strip. The term “sinus rhythm” is used when the rhythm originates in the sinus node and conducts to the ventricles. Young, athletic people may display various other rhythms, particularly during sleep. Sinus arrhythmia is the variation in the heart rate that occurs during inspiration and expiration. There is “beat to beat” variation in the RR interval, the rate increasing with inspiration. It is a vagally mediated response to the increased volume of blood returning to the heart during inspiration. cardiac axis The cardiac axis refers to the mean direction of the wave of ventricular depolar- ization in the vertical plane, measured from a zero reference point. The zero reference point looks at the heart from the same viewpoint as lead I. An axis lying above this line is given a negative number, and an axis lying below the line is given a positive number. Theoretically, the cardiac axis may lie anywhere between 180 and -180 ◦ . The normal range for the cardiac axis is between -30 ◦ and 90 ◦ . An axis lying beyond -30 ◦ is termed left axis deviation, whereas an axis > 90 ◦ is termed right axis deviation. Several methods can be used to calculatethe cardiac axis, though occasionally itis extremely di fficult todetermine. 32 Figure2.10: Hexaxialdiagram(projectionofsixleads inverticalplane)showingeachlead’s view of the heart. Table 2.1: Calculating the cardiac axis. Normal Axis Right Axis Deviation Left Axis Deviation Lead I Positive Negative Positive Lead II Positive Positive / Negative Negative Lead III Positive / Negative Positive Negative The simplest method is by inspection of leads I, II, and III, as shown in Table 2.1. A more accurate estimate of the axis can be achieved if all six limb leads are examined. The hexaxial diagram shows each lead’s view of the heart in the vertical plane. The direction of current flow is towards leads with a positive deflection, away from leads with a negative deflection, and at 90 ◦ toaleadwithanequiphasic QR S complex. The axis is determined as follows: Choose the limb lead closest to being equiphasic. The axis lies about 90 ◦ to the right or left of this lead x With reference to the hexaxial diagram, inspect the QR S complexes in the leads adjacent to the equiphasic lead. If the lead to the left side is positive, then the axis is 90 ◦ to the equiphasic lead towards the left. If the lead to the right side is positive, then the axis is 90 ◦ to the equiphasic lead towards the right. 33 Figure 2.11: Complex showing P wave highlighted. 2.5 Basic terminologies Thisarticleexplainsthegenesisofandnormalvaluesfortheindividualcomponents of the wave forms that are seen in an electrocardiogram. To recognize electrocardiographic abnormalities the range of normal wave patterns must be understood. P wave- The sinoatrial node lies high in the wall of the right atrium and initiates atrial depolarization, producing the P wave on the electrocardiogram. Although the atria are anatomically two distinct chambers, electrically they act almost as one. They have relatively little muscle and generate a single, small P wave. P wave amplitude rarely exceeds two and a half small squares (0.25 mV). The duration of the P wave should not exceed three small squares (0.12 s). The wave of depolarization is directed inferiorly and towards the left, and thus the P wave tends to be upright in leads I and II and inverted in lead aVR. Sinus P waves are usually most prominently seen in leads II and V1. A negative P wave in lead I may be due to incorrect recording of the electrocardiogram (that is, with transposition of the left 34 and right arm electrodes), dextrocardia, or abnormal atrial rhythms. The P wave in V1 is often biphasic. Early right atrial forces are directed anteriorly, giving rise to an initial positive deflection; these are followed by left atrial forces travelling posteriorly, producing a later negative deflection. A large negative deflection (area of more than one small square) suggests left atrial enlargement. Normal P waves may have a slight notch, particularly in the pericardial (chest) leads. Bifid P waves result from slight asynchrony between right and left atrial depolarization. A pronounced notch with a peak-to-peak interval of >1mm (0.04s)isusually pathological, and isseeninassociation withaleftatrial abnormality–for example, in mitral stenosis. PR interval- After the P wave there is a brief return to the isoelectric line, resulting inthe “PR interval.”Duringthis timethe electrical impulse is conductedthrough the atrioventricular node, the bundle of His and bundle branches, and the Purkinje fibers. The PR interval is the time between the onset of atrial depolarization and the onset of ventricular depolarization, and it is measured from the beginning of the P wave to the first deflection of the QRS complex (see next section), whether this be a Q wave or an R wave. It therefore measures the time from the start of atrial depolarization to the start of ventricular depolarization [100]. The normal duration of the PR interval is three to five small squares (0.12-0.20 s). Abnormalities of the conducting system may lead to transmission delays, prolonging the PR interval. QRS complex- The QR S complex represents the electrical forces generated by ventricular depolarization. With normal intraventricular conduction, depolarization occurs inane fficient,rapidfashion. Thedurationofthe QR S complexismeasuredintheleadwith 35 Figure 2.12: Normal duration of PR interval is 0.12-0.20 s (three to five small squares). the widest complex and should not exceed two and a half small squares (0.10 s). Delays in ventriculardepolarization–forexample,bundlebranchblock–giverisetoabnormallywide QRS complexes ( >0.12 s). The depolarization wave travels through the interventricular septum via the bundle of His and bundle branches and reaches the ventricular myocardium viathePurkinjefibernetwork. Theleftsideoftheseptumdepolarizesfirst,andtheimpulse then spreads towards the right. Lead V1 lies immediately to the right of the septum and thusregistersaninitialsmallpositivedeflection( Rwave)asthedepolarizationwavetravels towards this lead. When the wave of septal depolarization travels away from the recording electrode, the first deflection inscribed is negative. Thus small “septal” Q waves are often present in the lateral leads, usually leads I, aVL, V5, and V6. These non-pathological Q waves are less than two small squares deep and less than one small square wide, and should be < 25% of the amplitude of the corresponding R wave. The wave of depolarization reaches the endocardium at the apex of the ventricles, and then travels to the epicardium, spreading 36 Figure 2.13: Composition of QRS complex. outwards in all directions. Depolarization of the right and left ventricles produces opposing electrical vectors, but the left ventricle has the larger muscle mass and its depolarization dominates the electrocardiogram. In the precordial leads, QR S morphology changes de- pending on whether the depolarization forces are moving towards or away from a lead. The forces generated by the free wall of the left ventricle predominate, and therefore in lead V1 asmall R wave is followed by a large negative deflection ( S wave). The R wave in the precordial leads steadily increases in amplitude from lead V1 to V6, with a corresponding decrease in S wave depth, culminating in a predominantly positive complex in V6. Thus, the QR S complex gradually changes from being predomi- nantly negative in lead V1 to being predominantly positive in lead V6. The lead with an equiphasic QRS complex is located over the transition zone; this lies between leads V3 and V4, but shifts towards the left with age. The height of the R wave is variable and increases progressively across the precordial leads; it is usually < 27 mm in leads V5 and V6. The R wave in lead V6, however, is often smaller than the R wave in V5, since the V6 electrode is further from the left ventricle. 37 The S waveisdeepestintherightprecordialleads;itdecreasesinamplitudeacross the precordium, and is often absent in leads V5 and V6. The depth of the S wave should not exceed 30 mm in a normal individual, although S waves and R waves > 30 mm are occasionally recorded in normal young male adults. ST segment- The QR S complex terminates at the J point or ST junction. The ST segment lies between the J point and the beginning of the T wave, and represents the period between the end of ventricular depolarization and the beginning of repolarization. The ST segment should be level with the subsequent “TP segment” and is normally fairly flat, though it may slope upwards slightly before merging with the T wave. In leads V1 to V3 the rapidly ascending S wave merges directly with the T wave, making the J point indistinct and the ST segment di fficult to identify. This produces elevation of the ST segment,andthisisknownas“hightake-o ff.”Non-pathologicalelevationofthe ST segment isalsoassociatedwithbenignearlyrepolarization(seearticleonacutemyocardialinfarction later in the series), which is particularly common in young men, athletes, and black people. Interpretation of subtle abnormalities of the ST segment is one of the more di fficult areas of clinical electrocardiography;nevertheless, any elevation or depression of the ST segment must be explained rather than dismissed. T wave- Ventricular repolarization produces the T wave. The normal T wave is asymmetrical, the first half having a more gradual slope than the second half. T wave orientation usually corresponds with that of the QRS complex, and thus is inverted in lead aVR, and may be inverted in lead III. T wave inversion in lead V1 is also common. It is occasionally accompanied by T wave inversion in lead V2, though isolated T wave 38 Figure 2.14: The ST segment should be in the same horizontal plane as the TP segment; the J point is the point of inflection between the S wave and ST segment. inversioninleadV2isabnormal. T wave inversion in two or more of the right precordial leads is known as a persistent juvenile pattern; it is more common in black people. The presence of symmetrical, inverted T waves is highly suggestive of myocardial ischemia, though asymmetrical inverted T waves are frequently a non-specific finding. No widely accepted criteria exist regarding T wave amplitude. As a general rule, T wave amplitude corresponds with the amplitude of the preceding R wave, though the tallest T waves are seen in leads V3 and V4. Tall T waves may be seen in acute myocardial ischemia and are a feature of hyperkalaemia. QT interval- The QT interval is measured from the beginning of the QRS com- plex to the end of the T wave and represents the total time taken for depolarization and repolarizationoftheventricles. Thereforeitincludesalloftheelectricaleventstakeplacein the ventricles. From the stand point of time, more of the QT interval is devoted to ventric- ular repolarization than depolarization (i.e. the T wave is wider than the QRS complex). The duration of the QT interval is proportional to the heart rate [100]. The QT interval lengthens as the heart rate slows, and thus when measuring the QT interval the rate must 39 Figure 2.15: Complex showing T wave highlighted. betakenintoaccount. Asageneralguidethe QT interval shouldbe0.35-0.45 s,andshould not be more than half of the interval between adjacent R waves ( R-R interval). In general, the QT interval composes about 40% of the normal cardiac cycle, the R- R interval [100]. The QT interval increases slightly with age and tends to be longer in women than in men. Bazett’s correction is used to calculate the QT interval corrected for heart rate (QT c): QT c = QT / √ RR (seconds) . Prominent U waves can easily be mistaken for T waves, leading to overestimation of the QT interval. This mistake can be avoided by identifying a lead where U waves are not prominent–for example, lead aVL. U wave- The U wave is a small deflection that follows the T wave. It is generally upright except in the aVR lead and is often most prominent in leads V2 to V4. U waves result from repolarization of the mid-myocardial cells–that is, those between the endo- cardium and the epicardium–and the His-Purkinje system. Many electrocardiograms have no discernible U waves. Prominent U waves may be found in athletes and are associated 40 Figure 2.16: The QT interval is measured in lead aVL as this lead does not have prominent U waves (diagram is scaled up). with hypokalaemia and hypercalcaemia. 2.6 Heart Disease Normal sinus rhythm is arbitrarily limited to impulse formation beginning in the sinusnodeatratesbetween60and100beats/min. Ratesbelow50beats/minareconsidered to be bradycardia and rates greater than 100 beats/min are considered to be tachycardia. The normal sequence of electrical activation of the heart is from the sinus node through the atria to the atrioventricular ( AV) node and His-Purkinje system and to ventricular myocardium [85]. Arrhythmias resulting in bradycardia or tachycardia can be thought of as specific disorders of each of these components [69]. 2.6.1 Tachyarrhytmias Tachyarrhytmias are broadly characterized as being superventricular tachycardia (AVT), defined as a tachycardia in which the driving circuit or focus originates in tis- sue above the ventricle (i.e., sinus node, atria, AV node, or His bundle) or a ventricular 41 tachycardia (VT), defined as a tachycardia in which the driving circuit solely originates in ventricular tissue or Purkinje fibers [69]. In general, VT implies the presence of significant heart disease, whereas SVT is usually not lethal. Usually a narrow QRS complex is associ- ated withSVT, but a wide complex tachycardiacan be SVT or VT. It is oftenvery di fficult to detect VT. 2.6.2 Disorders of Impulse Conduction Delay or block in conduction can result in bradyarrhythmias or tachyarrhythmia. Bradyarrhythmias occur when propagating impulse is blocked and is followed by asystole or slow scope. Tachycardia occurs due to reentrant excitation due to delay and block. More commonly, impulses are blocked at rapid rates as a result of incomplete recovery of refractoriness caused by time- or voltage dependent excitability. For example, such incomplete recovery is the usual mechanism responsible for a non conducted premature P wave or one that conducts with a functional bundle branch block [85]. Reentry Ecentrical activity during each normal cardiac cycle begins in the sinus node and continues until the entire heart has been activated. Then heart fibers are completely refrac- tory,thenextcellexcitationstartsatsinusnode. It,however,agroupoffibersnotactivates during the initial wave of depolarization, recover excitability in time to be discharged be- fore the impulse dies out. They may reexcite areas that have just been discharged. Such a processiscalledreentry. Reentryisprobablythecauseofmanytachyarrhythmias,including various types of SVT and VT and fibrillation [69]. 42 2.7 Cardiac Arrest and Sudden Cardiac Death Sudden CardiacDeath (SCD) isdefined as natural death from cardiac causes such ascardiacarrestandcardiovascularcollapse. Electrophysiologicalabnormalities,diseasesof the AV nodeandHis-Purkinjesystemandthepresenceofaccessorypathwaysofconduction are two groups of structural abnormalities of specialized conduction that may be associated with SCD [66]. Substantial evidence now shows that genetically determined abnormalities of proteins that control electrical activity of the heart can cause cardiac arrest in the struc- turallyintactheart. Atleast14geneshavebeenassociatedwithinheritedarrythmo-genetic disease and are linked to sudden death in the apparently normal heart [76]. 2.7.1 Long QT Syndrome Long QT Syndrome is an inherited arrythmogenetic disease characterized by sus- ceptibility to life-threatening arrhythmias. TwomajorformsofLQTShavebeenidentified, one transmitted as an autosomal dominant trait (Romano Ward Syndrome), the second transmitted as an autosomal recessive disease in which the cardiac phenotype is associated with neurosensory deafness (Jervell and Lange-Nelsen syndrome) [76]. The electrocardio- graphic marker of LQTS consists of prolonged repolarization (i.e. prolonged QT interval), abnormal morphology of the T wave and of a characteristic polymorphic ventricular tachy- cardia called "torsades de pointes" that is usually induced by activation of the sympathetic nervoussystembutalsobyasuddenpause[92]. Torsadedepointesiscommonlythespecific arrhythmia that triggers or degenerates into Ventricular Fibrillation (VF). The congenital LQTS is a functional abnormality caused by hereditary defects of molecular structure of 43 ion channel protein and is associated with environmental or neurogenic triggers that can initiate symptomatic or letal arrhythmias. More than 35 mutations in four cardiac ion channel genes—KVLQT1 (voltage-gated K channel gene causing one of the autosomal domi- nantformsofLQTS)(LQT1),HERG(humanether-a-go-gorelatedgene.) (LQT2),SCN5A (LQT3), and KCNE1 (minK, LQT5)—have been identified in LQTS. These genes encode ion channels responsible for three of the fundamental ionic currents in the cardiac action potential. Clinical manifestations of LQTS are usually syncope, fainting, cardiac arrest, and also sudden death [76]. Their occurrences often depend on the physical or emotional stress (e.g. fear, anger, loud noises, sudden awakening, central nervous system injury), but they may occur at rest. Higher levels of risk are associated with female gender, specifically, greater degree of QT prolongation or QT alternans, unexplained syncope etc. Patients with the syndrome require avoidance of drugs that are associated with QT lengthening and careful medical management, which may include implantable defibrillators [66]. From a physiologic standpoint, dispersion occurs with repolarization between 3 layersoftheheart,andtherepolarizationphasetendstobeprolongedinmyocardium. This iswhytheTwaveisnormallywideandtheintervalfromT peak to T end (T p-e )represents the transmural dispersion of repolarization(TDR).InLQTS,TDRincreasesandcreatesa functional substrate for transmural reentry. In LQTS, QT prolongation can lead to polymorphic ventricular tachycardia, or torsadedepointes,whichitselfmayleadtoventricularfibrillationandsuddencardiacdeath. 44 Torsade de pointes is widely thought to be triggered by reactivation of calcium channels, reactivation of a delayed sodium current, or a decreased outward potassium current that results in early after depolarization (EAD) in a condition with enhanced TDR usually associated with a prolonged QT interval. TDR serves as a functional reentry substrate to maintain torsade de pointes. BoththecongenitalandacquiredLQTSareduetoabnormalities(intrinsicand/or acquired)oftheioniccurrentsunderlyingrepolarization. Prolongationoftherepolarization phase acts as a primary step for the generation of EADs. EAD-induced triggered beats arise predominantly from the Purkinje network. In the LQTS, prolonged repolarization is associated with increased spatial dispersion of repolarization. The focal EAD-induced triggered beat(s) can infringe on the underlying substrate of inhomogeneous repolarization to initiate polymorphic reentrant ventricular tachyarrhythmia [20]. Currently, practically all major noninvasive methods of assessing an individual’s susceptibility to cardiac arrhythmias include some analysis of the QT and/or RR interval spatio-temporal distribution. Indeed, the QT interval dispersion is based on the assessment of myocardial repolarization inhomogeneity. Recent advances in computer technology have led to improvements in automatic analysisofheartrateand QT intervalvariations. Itiswellknownnowthatthe QT interval’s spatial variability (QT dispersion) observations performed separately or in combination with the heart rate (or RR interval) variability analysis provides a tool for the assessment of individual susceptibility to cardiac arrhythmias. 45 2.7.2 L QT - RR Relationship Currently no tests have high predictive value for sudden cardiac death. Long QT syndrome (LQTS) has been known as a cause of sudden cardiac death. LQTS can also lead to an abnormal heart rhythm (arrhythmia) [56]. A prolonged QT interval can increase the risk for atypeof arrhythmiacalledtorsade de pointes. When torsadedepointesoccurs, the heart cannot pump enough oxygen-rich blood to the rest of the body, especially the brain. Torsade de pointes can also lead to ventricular fibrillation, a dangerous form of arrhythmia that causes rapid, uncoordinated contractions in the muscle fibers of the ventricles. With ventricular fibrillation, the heart cannot pump oxygen-rich blood to the rest of the body, which can lead to death. Heart Rate Variability and prolonged QT interval have been the subject of active investigation in many research schools such as John Hopkins School of Medicine Heart Institute’s Arrhythmia Service, and Hospital’s Department of Medicine, and California Institute of Technology (Caltech). Dr. Ronald Berger, a Professor of Biomedical Engineering and Medicine, and the Director of Cardiac Electro-Physioloy Training Program at John Hopkins University, has focused his research on mechanisms of Sudden Cardiac Death (SCD). Dr. Berger’s patent, "Methodology for automated QT variability measurement," [7] in 1996, is related to analyzing electrocardiograph signals to determine risk of malignant arrhythmias. His work had been directed toward identifying electrical abnormalities associated with Dilated Cardiomyopathy (DCM). His studies on DCM have revealed significant prolongation of ac- tion potential duration (APD) in cells isolated from failing hearts, regardless of etiology, 46 compared to those taken from normal hearts. QT prolongation alone is believed to lack adequate specificity to serve as a useful predictor of malignant arrhythmias, and mechanis- tically may not be important without accompanying electrical abnormalities. Berger et al. [8] also determined the coherence between heart rate and QT interval fluctuations using spectral analysis. They have shown that the action potential duration is normally closely coupled with the RR interval through a variety of mechanisms, which may trigger an alert in the failing heart in the failing heart. They used an algorithm to quantify beat-to-beat fluctuations in the QT interval, and analyzed the relations of QT variability with heart rate variability in patients with DCM. In their method, the heart rate mean (HRm) and variance (HRv) and QT interval mean (QTm) and variance (QTv) were computed from the respective time series. A normalized QT variability index (QTVI) was then derived according to the equation QT V I=log 10 h QT v /Q T 2 m HR v /H R 2 m i . The QTVI represents the log ratio between the QT interval and heart rate vari- abilities. The HR exhibits substantial beat-to-beat variability, which is mirrored in the instantaneous QT interval; that is, when HR rises QT falls and vice versa. In this case, not surprisingly, heart rate variability is small. However, in contrast to the findings in the normal subject, in the DCM patient the QT interval varies widely and erratically. They found that a coupling between HR and QT interval variations exists in control subjects and is reduced in DCM patients. This suggests either that the mechanisms driving heart rate changes also modulate QT interval or that the heart rate fluctuations themselves sec- ondarily drive QT interval oscillations. The latter mechanism relates to the well-known phenomenon of electrical restitution, that is, the dependence of action potential duration 47 on the diastolic interval of the preceding beat. Hence, electrical restitution must be altered in these cells since there is no longer a unique action potential duration for a given diastolic interval. There are other researches on QT prolongation and QT- RR relationship. Fossa et al. [23] has introduced a technique to demonstrate the dynamic assessment of the beat-to- beat QT- RR interval relationship to di fferentiate QT interval prolongation e ffects. Their techniquemayprove usefultostudy QT interval heterogeneity, hysteresis, andbeat-to-beat HRV. QT- RR clouds with hysteresis during acceleration of HR not following the Bazett and Fridericia correction relationships were interpreted as QTc prolongation. What distinguishes our approach from the previous two is that, here, we claim that the causality relationship between RR and QT is the criterion to determine whether the heart is at risk for Torsade de Pointes, ventricular arrhythmia, or other life threatening diseases. Marsden and his research group at Caltech have been working on the measure of the separation rate of dynamical systems. They have been using the concept of the Poincaré homoclinic tangle (through the Lagrange Coherent Structure technique) to view thetransientchaosoverintermediatetimescaleandinvestigatethemixingpropertiesoftime dependent dynamical systems and their behavior in the boundary of stable and unstable manifolds. Theybelievethattheconceptsofdynamicalsystemsmayexplainthecomplexity of biological systems. "Maybe in the future, Poincaré homoclinic tangle will be substituted for a type of heart disease," said Marsden smiling in his visit to the University of Southern California. Clearly, Marsden’s point of view is the one that comes closest to ours. 48 Chapter 3 Abstract Dynamical System 3.1 Introduction The theory of dynamical systems is a major mathematical discipline closely in- tertwined with all main areas of mathematics. It has greatly stimulated research in many sciences and has given rise to the vast new areavaryinglycalledapplied dynamics, non- linear science, or chaos theory. Kolmogorov, Sinai, Arnold, Kouchnirenko, and Rohlin [5], [54], [81], [96], [95] developed the notion of "abstract dynamical system" as the unifying framework encompassing all systems, whether they are deterministic, chaotic or stochastic. A dynamical system generally includes the following ingredients: • A "phase space" Ω, which contains the possible states of the system, • A "time," which could be continuous or discrete. It could either extend only in the future (positive, irreversible, or noninvertible processes), or to the past as well (pos- itive/negative, reversible, or invertible processes). For discrete processes we assume 49 thatthis is anaturalcorrespondencetothesetofallintegers, andforcontinuous time to the set of real numbers. • A the time-evolution law. In general this is rule that allows us to predict the state of the system in any time t, knowing the previous states. Thus, the most general time-evolution law is time dependent and has infinite memory. In continuous systems there is a there is a one parameter family of time shifts. ϕ t : Ω7−→ Ω x7−→ T(x, t) These transformations are related to each other for di fferent time t. As for any time t+ s, ϕ t+ s evolves to the new state T(x, t + s) such that ϕ t+ s = ϕ t ◦ ϕ s . A reversible discrete- time dynamical systems is represented by a cyclic group © T n =( ϕ 1 ) n | n ∈Z ª of one-to-one transformation of Ω onto itself [52]. More precisely, an abstract dynamical system is a quadruple ( Ω , Σ,μ,T) where Σ is a σ − al g ebr a of subsets of Ω,thatis ∀ A i ∈ Σ =⇒∪ A i ∈ Σ ∩ A i ∈ Σ A c ∈ Σ where A c isthecomplementaryset. μisameasure(willbedefinedlater)andforprobability spaces μ( Ω)=1. T is the time shift operator or "flow," a measure-preserving map (will also be defined later in this chapter). 50 Atthispoint,weperceivethatabstractdynamicalsystemsandclassicaldynamical systems can be envisioned within the same conceptual framework. Recall that a classical dynamicalsystemisdefinedoveradi fferentiableRiemannmanifold. Thisisthemostsalient di fferencebetweenthetwo, intheabstractcase, Ωhasonlyameasurablestructure, whilein theclassicalcase Ωhas,inaddition,adi fferentiablestructure. Thequestionis,startingfrom an abstract dynamical system, with a measurable space structure, whether there are some nonlinear(measurable)transformationon Ωthatcouldendowitwithstrongerdi fferentiable properties, thereby leading to a classical dynamical system. This has been formulated by Arnold [5], and is called "classical realization of an abstract dynamical system."Clearly not all abstract systems have a classical realization. Whether this is possible depends on Kolmogorov-Sinai entropy of the system. For systems with infinite entropy the transition from "abstract" to "classical" is impossible. The concept of Kolmogorov-Sinai entropy and nonlinear transformation of dynamical systems will be discussed in subsequent chapters. In the theory of smooth dynamical systems or di fferentiable dynamics, the phase space possesses the structure of smooth manifold, and is endowed with di ffeomorphisms anditeratesofnoninvertibledi fferentiablemaps. Sincefinite-dimensionalsmoothmanifolds have a natural locally compact topology, the theory of smooth dynamical systems naturally draws upon the notions and results of topological dynamics, specifically in dealing with asymptotic behavior of these systems that are very complicated and nonsmooth. In fact, the concept of asymptotic behavior of these systems belongs to topological dynamics, that is, they are defined only in terms of topology not the di fferentiable structure [52]. In particular, some invariant sets for a smooth system, for example, attractors, may not have 51 any smooth structure. In such cases, a specific class of topological dynamical systems, symbolic dynamics, which occur as closed invariant subsets of the shift transformation in a sequence space, is studied. The part of the theory of smooth dynamical systems that concerns measure- theoretic properties of such systems are called smooth ergodic theory. The theory of "statistical properties of dynamical systems" becomes relevant as soonasthesystemisexcitedbeyondthesimplestbifurcation,sothattheprecisegeometrical information about the shape of the attractor or the motion on it is no longer available [19]. 3.2 Symbolic Dynamical Systems Symbolic dynamics is the practice of modelling a topological or smooth dynamical system by a discrete space consisting of infinite sequences of abstract symbols, each of which corresponds to a state of the system, with the dynamics (evolution) given by the shift operator. The symbolic dynamics typically gives information about the transients of the system which are often hard to study. The distinctive feature in symbolic dynamics is that time is measured in discrete intervals. So at each time interval the system is in a particular state. Each state is associated with a symbol and the evolution of the system isdescribedbyaninfinite sequence of symbols – represented e ffectively as strings. If the system states are not inherently discrete, then the state vector must be discretized, so as to get a coarse-grained description of the system [52]. Symbolic dynamical systems are a class of topological dynamical systems of par- ticularimportanceforthetheoryofsmoothdynamicalsystems. Onereasonisthatinmany 52 respects symbolic systems serve as model for smooth ones; it is often easier to see many properties in the symbolic case first and then to carry them over to the smooth cases. Sec- ond, restrictions of some smooth dynamical systems to various (often important) invariant sets look very much like symbolic systems. Furthermore, symbolic systems can be used to "code" some smooth systems. In symbolic dynamical systems, first we introduce the sequence spaces. Let us consider for each natural number N ≥ 2 the space Ω N ={ ω=(··· ,ω −1 ,ω 0 ,ω 1 ,···)| ω i ∈{0 ,1 ,··· ,N −1} for i ∈Z} of two-sided sequences of N symbols and a similar one-sided space Ω R N ={ ω=( ω 0 ,ω 1 ,ω 2 ,···)| ω i ∈{0 ,1 ,··· ,N −1} for i ∈N 0 } . We can define a topology by noting that each Ω N is the direct product ofZ copies of the finite set {0 ,1 ,··· ,N −1}, each with the discrete topology, and using the product topology. Then we define a cylinder subset as C n 1 ,··· ,n k α 1 ,··· ,α k ={ ω ∈ Ω N | ω n i = α i for i=1 ,··· ,k} where n 1 <n 2 <··· <n k and numbers α 1 ,··· ,α k ∈ {0 ,1 ,··· ,N −1}.Thenumber k of fix digits is called the rank of cylinder. We can also define topology by declaring that all cylinders are open sets and that they form a base for the topology. Since the complement of a cylinder is a finite union of cylinders, then every cylinder is also close [52]. The metric can be introduced as d λ (ω, ω 0 )= ∞ P n= −∞ | ω n − ω 0 n | λ |n| , for any fixed λ> 1. Cylinder sets are particularly useful in providing the base of the natural topology oftheproductofacountablenumberofcopiesofaset. If Ωisafiniteset,theneachelement 53 of Ω can be represented by a letter, and the countable product can be represented by the collection of strings of letters. 3.3 Hyperbolic Dynamic Hyperbolicdynamicsischaracterizedbythepresenceofexpandingandcontracting directions for the derivative. This is a situation where the di fferential alone provides strong local, semilocal or even global information about the dynamics. The hyperbolic structure is defined for a subset of a manifold, when its tangent bundle may be split into two invariant subbundles, one of which is contracting, and the other expanding. Let M be a compact smooth manifold, and let f : M → M be a di ffeomorphism. An f −invariant subset Ω of M is said to be hyperbolic (or to have a hyperbolic structure) if there is a splitting of the tangent bundle M restricted to Ω into a Whitney sum of two Df −invariant subbundle, E s and E u , the stable and unstable manifolds. Note that Df| E s and Df| E u representcontractionandexpansion, respectively. Therefore, foranyreal numbers λ and c,wehave ∃ 0 <λ< 1 and c>0: T Ω M = E s ⊕ E u and Df( x) E s x = E s f( x) , Df( x) E u x = E u f( x) , = ⇒ kDf n vk <cλ n k vk,v ∈ E s = ⇒ ° ° Df − n v ° ° <cλ n k vk,v ∈ E u ,n > 0 . 54 Indeed, under iteration the presence of these directions produces exponential rel- ative behavior of orbits on some sets, and this a ffords much insight into topological and measurable aspects of the orbit structure. This stretching and folding typically gives rise to complicated long-term behavior in these systems. The dynamics appears in many ways e ffectively random, even though these systems are completely deterministic. 3.4 Measure Theory Measureisadynamicsstatisticsinterplay. Themeasure μisdefendona σ −algebra of subsets of Ω, such that μ : A7−→ μ( A) ∈R + ∪{ ∞} : Σ −→R + ∪{ ∞} with properties: • μ( Ω)=1 • μ( ∅)=0 • μ( A)= μ ¡ T −1 A ¢ , T is measurable and μ is invariant under T Ameasurablemapisdefinedasanymap T : X −→ Y,where X and Y aretwosets with σ −algebrasA andB of subsets of X and Y, respectively, if T −1 ( B) ∈A, ∀ B ∈B.The triple (X,A,μ) with a σ −algebraA andameasure μ onAis ameasurespace, T : X −→ X is a measure preserving map and that μ is invariant under T,if μ( A)= μ( T −1 ( A)).Note that any continuous map is measurable [62]. The dynamic behavior of measure-preserving maps is the theme of ergodic theory. 55 3.4.1 Invariant Measures and Measure Preserving Maps An invariant measure reflects the asymptotic behavior of almost all points with respecttothatmeasure. Considermeasurablespaces(X,A,μ)and (Y,B,ν)andmeasurable transformation T : X −→ Y.Ifmoreover ν( B)= μ( T −1 ( B)) for every B ∈ B,then T is calledmeasurepreserving. Ifameasurepreservingmap T isinvertibleandtheinverse T −1 is measurable, then clearly T −1 is also measure preserving. Therefore T is an isomorphism in the category of measure spaces. We can also write ν = μ ◦ T −1 .If (X,A,μ)=(Y,B,ν),we call T ameasure preserving endomorphisms; wewillalsosaythatmeasure μis T −invariant, or that T preserves μ.Inthecaseof (X,A,μ)=(Y,B,ν) an isomorphism T is called automorphism. 3.4.2 Existence of Invariant Measures Considerthediscretetimedynamicalsystem Ωand T : Ω → Ω,when Ωiscompact and T is continuous. Ω ⊆R n implies that for any infinite sequence ω n , n ∈N 0 ,thereexists a subsequence ω n k , k ∈ {0 ,1 ,2, ...} that converges and lim k →∞ ω n k exists. Ω is endowed with a σ −algebra of subsets. By Krylov-Bogoliubov theorem [52] there exists an invariant measure μ for any subset A in the σ −algebra. Theorem 3 (Krylov-Bogoliubov Theorem) Any continuous map on a metrizable compact space has an invariant Borel probability measure. Proof. Detailed proof is found in [52]. Remark 4 The Krylov-Bogoliubov theorem does not claim that the measure is unique. 56 Remark 5 The invariant measure may be singular. Remark 6 If T is homeomorphism, μ a T-invariant measure, and A ⊂ Ω measurable then μ( T( A)) = μ( A). 3.5 Ergodicity Considerameasure-preservingmapofaspace (X,A,μ)andanintegrablefunction f : X →R.Wewantto find the condition under which the limit lim n →+ ∞ f( x)+ f( T( x))+···+ f( T n −1 ( x)) n (3.1) exists and is constant almost every where. In 1931 Birkho ff [9] proved that the above limit exists for T and f. He also showed its value to be constant a.e. if and only if there is no set A ∈A such that 0 <μ( A) < 1 and T −1 ( A)= A. Maps which satisfy this condition are called Ergodic. Theorem 7 (Birkho ff)Let (X,A,μ) be a probability space and T : X → X measure pre- serving map. a) If f ∈L 1 ( X),thelimit lim n →+ ∞ 1 n n −1 P j=0 f( T j ( x)) exists for almost every point x ∈ X. b) If f ∈L p ( X), 1 ≤ p ≤∞, the function ˜ f defined by ˜ f( x)= lim n→∞ 1 n n −1 P j=0 f( T j ( x)) 57 is in L p ( X) and satisfies lim ° ° ° ° ° ˜ f( x) − 1 n n −1 P j=0 f ◦ T j ° ° ° ° ° p =0 , ˜ f( T( x)) = ˜ f( x) a.e. c) For every f ∈L p ( X) we have R X ˜ fdμ = R X fdμ. Proof. Detailed proof is found in [62]. The ergodic theorem says that the long time behavior of the system is asymptoti- cally described by the behavior of the ergodic components of the space. It actually allows us to consider only the long-term behavior of a system and not worry about transients. This somehow simplifies the problems. For example, the physical long term behavior is on attractors, but the geometric study of attractors is mathematically di fficult. On the other hand, studying attractors using invariant measures makes computations easier. Thereby, the ergodic properties of systems are intertwined with the invariant measures [19]. In general, an invariant measure may be decomposable into several di fferent in- variant measures. However, the ergodic properties of dynamical systems is associated with theindecomposability. Hence, wecansayaninvariantmeasurecanbeuniquelyrepresented as a superposition of ergodic measures [19]. Now consider the orbital average of f by (3.1), called ˜ f.When f is the charac- teristic function f A of a set A ∈ A,thenumber ˜ f A ( x) is called the average time (sojourn time) spent by x in the set A, and is written τ A ( x) and is defined as τ A ( x)= lim n→∞ 1 n # © 0 ≤ j ≤ n −1| T j ( x) ∈ A ª . 58 Therefore, an immediate consequence of the ergodic theorem is that the measure of A is equal to the average sojourn time. μ( A)= R X τ A dμ=lim n →∞ 1 n n −1 P j=0 τ A ( T j x) . We deduce that for almost every x the frequency of hitting A by the forward trajectory equalstothemeasure(probability)of A,namely lim n →∞ #{0 ≤ j< n : Tj( x) ∈ A} /n,isequal to μ( A). This means for example that if we choose a point in X in a euclidean space at random its su fficiently long forward trajectory fills X with the density being approximately the density of μ with respect to the Lebesgue measure, provided μ is equivalent to the Lebesgue measure. Consider a measure-preserving map T of a probability space (X,A,μ) and a T- invariant set A ∈ A,sothat T −1 ( A)= A.Bydefinition, map T is ergodic if every T- invariant set has measure 0 or 1. Definition 8 A measurable transformation T : X → X of a measure space (X,A,μ) is said to be ergodic if for any measurable set A, μ( T −1 ( A)÷ A)=0 ⇒ μ( A)=0 or μ( X\ A)=0. (Recall the notation B÷ C=( B\ C) ∪( C\ B).) Note that we did not assume in the definition of ergodicity that μ is T-invariant (neither that μ is finite). Equivalently, for every A, B ∈A we have lim n →∞ 1 n n −1 P m=0 μ( T − m ( A) ∩ B)= μ( A) μ( B) (3.2) and that τ A ( x)= μ( A) almost everywhere. 59 Consider the "time average" of a well-behaved function f. The time average of function f is defined as the average (if it exists) over iterations of T starting from some initial point x. ˜ f( x)= lim n→∞ 1 n n −1 P j=0 f( T j x) Consider also the "space average" or "phase average" of f,defined as ¯ f = R fdμ In general, the time average and space average may be di fferent. But if the trans- formation is ergodic, and the measure is invariant then the time mean is equal to the space mean almost everywhere. This is the celebrated ergodic theorem. A general measure-preserving system is not necessarily ergodic, but we shall in- troduce the ergodic decomposition, which allows one to express any non-ergodic measure as an average of ergodic measures (generalizing the decomposition of a permutation into disjoint cycles). 3.5.1 Uniquely Ergodic Maps A continuous map T of a compact metric space X may happen to have exactly one invariant measure. Such maps are called uniquely ergodic [62]. Theorem 9 [62] Let X be a compact metric space and T : X −→ X a continuous map. The following properties are equivalent: a) T is uniquely ergodic. 60 b) For every f ∈ C 0 ( X) the limit lim n−→∞ 1 n+1 n P j=0 f( T j ( x)) exists for every x and does not depend on x. c) For every f ∈ C 0 ( X) the sequence of functions 1 n+1 n P j=0 f ◦ T j converges uniformly to a constant. 3.5.2 Mixing Dynamical Systems Inphysics,adynamicalsystemissaidtobemixingifthephasespaceofthesystem becomes strongly intertwined, according to at least one of several mathematical definitions. Definition 10 A topological dynamical system T : Ω7−→ Ω is called topologically mixing if for any two open sets A, B ⊂ Ω there exists a positive integer N = N(A, B) such that for any n>N the intersection of T n ( A) ∩ B is nonempty. Since every topological mixing map is also topologically transitive, no translation could be topologically mixing. This follows from the fact that translation preserves the natural metric on the torus induced by the standard Euclidean metric onR n and from the following general criteria. Lemma 11 If a topological dynamical system T : Ω 7−→ Ω preserves a metric on Ω which generates the topology of Ω,then T is not topologically mixing. 61 Themathematicaldefinitionsofmixingaremeanttocapturethenotionofphysical mixing. A canonical example is the Cuba libre: suppose that a glass initially contains 20% rum (the set A)and80%cola(theset B) in separate regions. After stirring the glass, any region of the glass contains approximately 20% rum. Furthermore, the stirred mixture is in a certain sense inseparable: no matter where one looks, or how small a region one looks at, one will find 80% cola and 20% rum. Mixing is a stronger notion that ergodicity. Every mixing transformation is er- godic. An equivalent definition can be given in the language of measure-preserving dynam- ical systems. Let ( Ω , Σ,μ, T) be a dynamical system, with T being the time-evolution or shift operator. Then, if for all A, B ∈A, one has lim k →∞ μ( A ∩ T − k B)= μ( A) μ( B) then the system is called strongly mixing. A dynamical system is said to be weakly mixing if lim n→∞ 1 n n P k=0 ¯ ¯ ¯ μ( A ∩ T − k B) − μ( A) μ( B) ¯ ¯ ¯=0 Strong mixing implies weak mixing, and every weakly-mixing system is ergodic. Now consider mixing dynamical systems in which the function φ determines the rate of mixing while the separation function f specifies a lower bound for the size of the “gap” m that is necessary to get the mixing property. A measure-preservative dynamical system (T, μ) is called (φ, f)-Mixing if there exits a sequence φ( n),n=0 ,1,. .. of positive numbers satisfying ∞ P n=1 φ 1 /2 ( n) < ∞ 62 such that for every cylinder sets A and B ¯ ¯ μ( A ∩ T − m − n B) − μ( A) μ( B) ¯ ¯ ≤ φ( m) μ( A) μ( B) Recall that in symbolic dynamics, a cylinder sets is defined as the subsets of the space X ={1 ,2,.. .,M} Z = Q Z ={ x =··· x −1 x 0 x 1 ··· ,x i ∈ Q} (ϕ, f)-mixing implies weak mixing [40]. Definition 12 [45] (Speed of mixing). Let (X,B,T,μ) be a dynamical system and ξ a finite or countable measurable partition of X.Weset ξ k = ∨ K −1 j=0 T − j ξ and σ( ξ k ) the σ −algebra generated by ξ k . • Uniform mixing. The partition is uniformly mixing with speed γ( n) going to zero for n going to infinity if for any n, γ( n)=sup k,l sup A ∈ σ( ξ k ) B ∈ T −( n+ k) σ( ξ l ) | μ( A ∩ B) − μ( A) μ( B)| . • α-mixing. The partition is α-mixing with speed α( n) going to zero for n going to infinity if for any n, α( n)=sup k,l sup A ∈ σ( ξ k ) B ∈ T −( n+ k) σ( ξ l ) ¯ ¯ ¯ ¯ μ( A ∩ B) μ( A) − μ( B) ¯ ¯ ¯ ¯ . 63 • ϕ-mixing. The partition ξ is ϕ-mixing with speed ϕ( n)goingtozerofor n going to infinity if for any n, ϕ( n)=sup k,l sup A ∈ σ( ξ k ) B ∈ T −( n+ k) ξ l ¯ ¯ ¯ ¯ μ( A ∩ B) μ( A) μ( B) −1 ¯ ¯ ¯ ¯ . • Weak-Bernoulli. The partition ξ is weak-Bernoulli with speed β( n) going to zero when n goes to infinity, if for any n, β( n)=sup k,l sup A ∈ σ( ξ k ) B ∈ T −( n+ k) ξ l | μ( A ∩ B) − μ( A) μ( B)| . Remark 13 We state some general implications and results verified by the preceding types of mixing. 1. ϕ-mixing implies α-mixing which implies uniform mixing. For any n, γ( n) ≤ α( n) ≤ ϕ( n). 2. ϕ-mixingimpliesweak-Bernoullimixingwhichimpliesuniformmixing. Forany n, γ( n) ≤ β( n) ≤ ϕ( n). 3. If ξ is a generating partition of an uniformly mixing dynamical system, then the system is mixing. 4. If ξ is a generating weak-Bernoulli partition then the system is metrically conjugated to a Bernoulli shift. 64 3.6 Recurrence The recurrence properties such as recurrence of an orbit, minimality, topological transitivity, and topological mixing can be expressed in a quantitative way by considering asymptotic frequencies. Recall the definition of topologically transitive and minimality as follows. Definition 14 [52] A topological dynamical system T : Ω 7−→ Ω is topologically transitive if there exist a point x ∈ Ω such that its orbit O T ( x)={ T n ( x)} n ∈Z is dense in Ω. Definition 15 [52] A topological dynamical system T : Ω 7−→ Ω is called minimal if the orbit of every point x ∈ Ω is dense in Ω, or, if T has no proper closed invariant set. Bothtopologicalmixingandminimalityarestrongernotionsthantopologicaltran- sitivity [33]. In order to show that for some orbits such asymptotic frequencies exist we have to appeal to measure theory. Another aspect of asymptotic behavior is related to the type of recurrence with respect to the time. Topological transitivity implies that iterates of any open set over the time intersects any other open set. A stronger version of recur- rence is reflected by the concept of mixing. Nontrivial recurrence is defined as irritates of a point coming back arbitrary close to the initial position without returning exactly to that position. 65 Set A Initial condition Iterates Return to set A Set A Initial condition Iterates Return to set A Figure 3.1: Poincaré recurrence time. 3.6.1 Return Times of Dynamical System The Poincaré Recurrence Theorem The Poincaré recurrence theorem states that a system having a finite amount of energy and confined to a finite spatial volume will, after a su fficiently long time, return to an arbitrarily small neighborhood of its initial state. The Poincaré recurrence time is the amount of time elapsed until the recurrence. Theorem 16 Poincaré Recurrence Theorem–For any S ∈ Σ, the set of those points of ω of S such that f n ( ω) / ∈S for all n> 0 has measure zero. That is, almost every point of S returns to S. In fact, almost every point returns infinitely often; i.e. μ({ ω ∈S : there exists N such that f n ( ω) / ∈S for all n>N})=0 Given a subset S ∈ Ω, the successive return times of the state ω in S are denoted as τ S . By a recently developed theoretical paradigm, the distribution of τ S reveals the dynamical properties of a system [38], [36], [40]. 66 3.6.2 Return Times Statistics Thestatisticsofreturntimes(alsoknownasPoincarérecurrences)hasbeenexten- sively studied over the last few years [40], [12], [45], [84], mainly to characterize the ergodic and statistical properties of dynamical systems. It has been proved that the asymptotic distribution of Poincaré recurrences is exponential for a wide class of mixing systems, even if they are not uniformly hyperbolic [84]. N. Haydn and S. Vaienti have proved that, for a large class of mixing dynamical systems, the return event is a Poisson process [40]; hence the return times are independent and exponentially distributed. We briefly recall some of the basic definitions about the statistics of return times from [84]. Consider a transformation T on a measurable space Ω with an invariant measure μ with respect to the dynamics T. For a given measurable subsetS ∈ Σ and a point ω ∈S, the first return time of ω into S is defined as: τ S ( ω)=min ³n k ∈N : T k ( ω) ∈S o ∪{+ ∞} ´ and the mean return time is given by h τ S i = Z S τ S ( ω)dμ S where μ S denotes the conditional measure with respect to S : μ S ³ ˜ S ´ = μ(S ∩ ˜ S) μ(S) for any measurable ˜ S ⊆ Ω. For ergodic dynamics, Kac’s theorem [75] says that the expectation of the return time to S starting from S,isjust μ(S) −1 . Theorem 17 (Kac’s Theorem) Let T : Ω → Ω be an ergodic transformation on the prob- ability space ( Ω , Σ,μ).Let S ∈ Σ have μ(S) > 0 then we definethereturntimefunction 67 τ S : S →Z + ∪{ ∞} (which is finite almost everywhere). The average return time (with respect to induced probability measure μ S )is Z S τ S ( ω) dμ S = 1 μ(S) Proof. The proof could be found in [75]. The statistics of the first return times then is defined as: F S ( t)= μ S ({ ω ∈S : τ S ( ω) /h τ S i >t}) The limit statistics F( t)=lim μ(S) →0 F S ( t) exists when the set S shrinks around apoint ω ∈ Ω. For a wide class of mixing systems (see [34], [35], [40]), it has been proved that the limit spectrum decays exponentially, that is F( t)= e − t ,if S is taken either as a ball or as a cylinder shrinking around μ-almost every point. ForaparticulardomainS thatintersectstheboundaryofthetwoinvariantregions and whose measure goes to zero, the limit statistics of first return times is simply ruled by the one of the mixing region, thus being exponential. In the case of systems that split into invariant regions, e.g. coupled dynamics, Haydn et al. [37] proved that the distributions for domains intersecting the boundaries are a linear superposition of the distributions char- acteristic of each region. These results may allow us to understand the dynamical behavior of systems of physical interest where di fferent dynamics exist. 68 Chapter 4 Chaos 4.1 What is Chaos? Chaos, in the technical sense, is used to denote a type of time evolution in which the di fference between two states that are initially closely similar grows exponentially over time. Systemsendowedwithsuchbehaviorrequireadissipativemechanismtopreventthem fromblowingapart[14]. Theyareaperiodicandboundedatleastfortimescalesoverwhich numerical computations are feasible [98]. Chaos is more easily understood through a comparison with randomness and peri- odicity. Random behavior never repeats itself and is inherently unpredictable and disorga- nized. Unlike random behavior, periodic behavior is highly predictable, because it always repeats itself over some finite time interval. In fact, chaos is distinct from periodicity and randomness,buthascharacteristicsofboth. Itlooksdisorganized,butisactuallyorganized. The most important criteria for chaotic behavior are summarized as follows [24]: 69 1. Chaos is deterministic and aperiodic and it never repeats itself exactly. There are no identifiable cycles that recur at regular intervals. 2. Chaotic systems have sensitive dependence on the initial conditions. In other words, very small di fferences in the initial conditions will later result in large di fferences in behavior. 3. Chaotic behavior is constrained. Although the system appears random, the behavior is bounded, and does not wander o ff to infinity. 4. Chaotic behavior has definite attributes. The behavior is constrained, and there is a particular pattern to the behavior. It was believed that nonlinearity is necessary and fundamental to chaos, which is not true in general [25], [77], [87]. The basic quantities to characterize chaotic behavior are the exponential diver- gence of nearby orbits (positive Lyapunov exponents), positive finite Kolmogorov entropy, and a non-integer dimension of the attractor. These quantities are invariant under smooth transformation of coordinates [89]. The following sections deal with classical computational techniques of chaos, level of uncertainty and complexity of the dynamical systems on a more mathematical level. 4.1.1 Lyapunov Exponents For a dynamical system, sensitivity to initial conditions is quantified by the Lya- punov exponents [82]. The Lyapunov exponent measures rate of orbital divergence or convergence of smooth dynamical systems [47]. For chaotic systems, it is a quantitative 70 measure of separation of the trajectories that diverge widely from their initial close posi- tions. There are as many Lyapunov exponents as there are dimensions in the state space of the system, but the largest is usually the most important. The magnitude of this exponent is related to how chaotic the system is; the larger the exponent, the more chaotic the sys- tem. Asdiscussed,thesignandthemagnitudeofLyapunovexponentsrevealtheunderlying dynamics of the system, although these exponents are not easily measured. For periodic signals, theLyapunovexponentiszero. Arandomsignalwillalsohaveanexponentof zero. A positive Lyapunov exponent indicates sensitive dependence on the initial conditions and is a diagnostic of chaos [82]. Any continuous time-dependent dynamical system without a fixed point will have at least one zero exponent [29], corresponding to the slowly changing magnitude of a principal axis tangent to the flow [25]. The sum of the Lyapunov exponents is the time-averaged divergence of the phase space velocity; hence any dissipative dynamical system will have at least one negative exponent, the sum of all of the exponents is negative, and the post transient motion of trajectories will occur on a zero volume limit set, an attractor. The exponential expansion indicated by a positive Lyapunov exponent is incompatible with motion on a bounded attractor unless some sort of folding process merges widely separated trajectories. Each positiveexponentreflectsadirectioninwhichthesystemexperiencestherepeatedstretching and folding that decorrelates nearby states on the attractor. Therefore, the long-term behavior of an initial condition that is specified with any uncertainty cannot be predicted; this is chaos. An attractor for a dissipative system with one or more positive Lyapunov exponentsissaidtobe strange or chaotic [103]. 71 The flow map Φ t :R n → R n x 7→ Φ t ( x) (4.1) describing the dynamical system acts on the n-dimensional state space M = R n and is generated by a vector field v ˙ x = v( x),x ∈R n ,t ∈R (4.2) To gather information about the time evolution of infinitesimally small perturbed initial states, the linearized flow map D x Φ t : T x M → T Φ t ( x) M u 7→ D x Φ t u (4.3) has to be considered. The linearized flow map D x Φ t is given by an invertible n× n matrix describing the time evolution of a vector u in the tangent space. For ergodic systems the Lyapunov exponents are defined as the logarithms of the eigenvalues μ i (1 ≤ i ≤ m) of the positive and symmetric limit matrix Λ x =lim t→∞ ¡ D x Φ t ∗ • D x Φ t ¢ 1 2 t asgivenbythetheoremofOseledec[70](*denotestransposition). TheLyapunovexponents are the logarithmic growth rates λ i =lim t→∞ 1 t ln ° ° D x Φ t e i ° ° , (1 ≤ i ≤ m) (4.4) where { e i :1 ≤ i ≤ m} are basis vectors that span the eigenspaces of Λ x [4]. 72 Figure 4.1: Geometrical definition for Lyapunov exponents. Recognizing that the length of the first principal axis is proportional to e λ 1 t ,the areadeterminedbythefirsttwoprincipalaxesisproportionalto e ( λ 1 + λ 2 ) t ,andconsequently, thevolumedeterminedbythefirst k principalaxesisproportionalto e ( λ 1 + λ 2 +···+ λ k ) t .Thus, the Lyapunov spectrum can be defined such that the exponential growth of a k-volume element is given by the sum of the k largest Lyapunov exponents. Note that information created by the system is represented as a change in the volume defined by the expanding principal axes [82]. From a geometrical point of view, to obtain the Lyapunov spectra, imagine an infinitesimal small ball with radius dr sitting on the initial state of a trajectory. The flow will deform this ball into an ellipsoid, such that after a finite time t all orbits that have started in that ball will be in the ellipsoid. Then, the i th Lyapunov exponent is defined by λ i =lim t→∞ 1 t µ dl i ( t) dr ¶ where dl i ( t) is the radius of the ellipsoid along its i th principal axis (Figure 4.1). TheseparationmustbemeasuredalongtheLyapunovdirectionswhichcorrespond to the principal axes of the ellipsoid previously considered. These Lyapunov directions are 73 dependent upon the system flow and are defined using the Jacobian matrix, i.e., the tan- gent map, at each point of interest along the flow which rules the growth of infinitesimal perturbations [82], [55]. Hence, it is important to preserve the proper phase space orien- tation by using a suitable approximation of the tangent map. This requirement, however, becomes unnecessary when calculating only the largest Lyapunov exponent. If we assume thatthereexistsanergodicmeasureforthesystem, thenthemultiplicativeergodictheorem of Oseledec justifies the use of arbitrary phase space directions when calculating the largest Lyapunov exponent with smooth dynamical systems. In ergodic systems most trajectories will yield the same Lyapunov exponent, as- ymptotically for long times [17]. The computation of the full Lyapunov spectrum requires considerably more e ffort than just the maximal exponent. As discussed before, an essential ingredient is some esti- mate of the local Jacobians, i.e., of the linearized dynamics. One either finds it from direct fitsoflocallinearmodelsofthetype s n+1 = a n s n + b n ,suchthatthefirstrowoftheJacobian is the vector a n ,and ( J) i, j = δ i −1 ,j for i=2,.. .m,where m is the embedding dimension. The a n isgivenbytheleastsquaresminimization σ 2 = P l ( s l+1 − a n s l − b n ) 2 where{ s l }isthe set of neighbors of s n [89], [18], [11]. Or one constructs a global nonlinear model and com- putes its local Jacobians by taking derivatives. In both cases, one multiplies the Jacobians one by one, following the trajectory, to as many di fferent vectors u k in tangent space as one wants to compute Lyapunov exponents. Every few steps, one applies a Gram-Schmidt orthonormalization procedure to the set of u k , and accumulates the logarithms of their rescaling factors. Their average, in the order of the Gram-Schmidt procedure, give the Lya- 74 punov exponents in descending order [97]. The routine used in this research, “lyap_spec”, uses this method, which goes back to [18] and [89], employing local linear fits. The code is borrowed from the website of Nonlinear Time Series Analysis [41]. The algorithm can be found in Appendix A. 4.1.2 Kolmogorov-Sinai Entropy The Kolmogorov-Sinai (KS) entropy measures how chaotic a dynamical system is. In the case of deterministic chaos, this quantity is positive and measures the average rate at which the information about the state of the system is lost over time. In other words, entropy is inversely proportional to the time interval over which the state of the system can be predicted. Assume that μ is an ergodic probability measure for a dynamical system. We definetheconceptofmeasure-theoretic entropyortheKolmogorov-Sinai invariantorsimply entropy as h( μ) [19]. Assume thatA=(A 1 ,.. .,A α ) be a finite ( μ-measurable) partition of thesupportof μ. Considerthe k th preimageofeachpartitionA j as T − k A j .Hence,by T − k A we denote the partition ( T − k A 1 ,. ..,T − k A α ). Finally, we refine the partitioning and define it as A ( n) = A ∨ T −1 A ∨··· ∨ T − n+1 A, which its pieces are A i 1 ∩ T −1 A i 2 ∩··· ∩ T − n+1 A i n with i j ∈ {1 ,2,. ..,α} [19]. Note that the partition T − k A is deduced from A by time evolution, whereas T k A is not necessarily a partition since in general the mapping T might be to many-to-one. We define the Kolmogorov-Sinai entropy as h( μ)= lim diam A →0 sup all measurable partitionings h( μ,A) 75 where the entropy rate is h( μ,A)=lim n→∞ £ H(A ( n+1) ) − H(A ( n) ) ¤ =lim n→∞ 1 n H(A ( n) ) and H(A)= − α P i=1 μ(A i )ln μ(A i ) (diam A=max i {diameter of A i }) [19], [48]. The relationship between the entropy and the Lyapunov exponents is very inter- esting. There is a general inequality introduced by the following theorem [19], [24], [86]. Theorem 18 (Ruelle 1978) Let T be a di fferentiable map of a finite-dimensional manifold and μ an ergodic measure with compact support. Then h( μ) ≤ X λ i >0 λ i . (4.5) Pesin[72]provedthattheequalityholdsif μisinvariantunderthedi ffeomorphism T,and μ has smooth density with respect to Lebesgue measure [82], [24]. The Kolmogorov-Sinai (KS) entropy can be evaluated quantitatively and is diag- nostic of chaos, whereassome other methods such as spectral analysis, time autocorrelation function and scatter plot construction are qualitative methods. 4.1.3 Kaplan-Yorke (Lyapunov) Dimension Lyapunov dimension is another Fractal dimension and a quantitative measure of chaos; a non-integer Lyapunov dimension is an indication of chaos. It was introduced by Kaplan and Yorke [51] based on the Lyapunov exponents. By definition [80], if Φ t is a map onR n and O + Φ ( x 0 ) is a bounded forward orbit having Lyapunov exponents λ j = λ j ( x 0 ; Φ), with, for the integer k such that λ 1 + λ 2 +···+ λ k ≥ 0 λ 1 + λ 2 +···+ λ k + λ k+1 < 0 76 then, the Lyapunov dimension of the orbit is dim L ( O + Φ ( x 0 )) = k + P k i=1 λ i | λ k+1 | (4.6) Notice that λ 1 + λ 2 +···+ λ k < | λ k+1 |,so dim L ( O + Φ ( x 0 )) <k+1. If the attractor has a positive Lyapunov exponent, then k ≥ 1. The fractal dimensions of chaotic flows are shown to be given D = m 0 + m + {1+ | λ + /λ − |},where m 0 and m + are the numbers of zero and positive Lyapunov characteristic exponents λ α and λ ± are the mean values of positive and negative λ α respectively [65]. 4.1.4 Scatter Plot Analysis The scatter plot, a technique taken from nonlinear dynamics, portrays the nature of fluctuations in the dynamical systems. It is a plot in which each state is plotted as a function of the previous one. The geometry of the scatter plot is essential and can be described by fitting an ellipse to the graph. The ellipse is fitted onto the so called line-of- identityat 45 ◦ tothenormalaxis. Theanalysescouldbedonequantitativelybycalculating the standard deviations of the distances of the state i to the lines y = x and y = − x+2 m, where m is the mean of all states i. The standard deviation of the points perpendicular and along to the line-of-identity (denoted by SD1 and SD2) describes short-term and long term variabilities, respectively [79]. Statistically, the plot displays the correlation between consecutive intervals in a graphical manner. Figure 4.2 is a typical illustration of the scatter plot. Here, nonlinear dynamicsconsidersthescatterplotasthetwodimensional(2-D)reconstructedphase-space, which is a projection of the reconstructed attractor describing the dynamics of the system. 77 Figure 4.2: A typical scatter plot. 4.2 Stationarity Ascientificmeasurementofanykindisinprinciplemoreusefulthemoreitis reproducible. In case of time series measurements, reproducibility is closely connected to di fferent notions of "Stationarity." Time series analysis methods give rise to algorithms which just compress data into a set of a few numbers. The results, however, can not be assumed to characterize the underlying system that produces non-stationary data. Hence, non-stationarity can not be simply ignored and since the stationarity can not be construc- tively proven for a given finite time series, the natural assumption would be that a given sequence of data is not stationary. The reason for testing stationarity is that we try to approach a dynamical phenomenon through a single finite time series, and hence is a re- quirement of almost all statistical tools for time series data, including the linear ones. There are many processes which are formally stationary when the limit of infi- nitely long observation times can be taken, but which behave e ffectively like non-stationary 78 processes when studied over finite time. A prominent phenomenon belonging to this class is called “intermittency.” It is thus necessary to devise not only tests for stationarity, but methods to determine over what time scale a time series can be considered stationary [58]. Itisneitheranecessarynorasu fficientconditionforstationaritythatthesystemproducing the time series must remainunchanged duringthe time of measurement. The reason is that there is no a posteriori distinction between a system parameter that remains constant and a variable (which may evolve in time). Thus a system with a rapidly fluctuating parameter may yield a stationary time series because these fluctuations can be averaged over, while a system with constant parameters can be considered nonstationary. An example for the lattercaseisgivenbyintermittencywherethetimeevolutionofsomevariablesmaybecome arbitrarily slow. For a stationary time series, the mean value is invariant under random shu ffling. If the process under observation is a probabilistic one, for a stationary process, these prob- abilities may not depend on time. Almost all the known methods and results on time series analysis assumes the validity of following conditions: • Parameters of the system remain constant, • The phenomenon is su fficiently sampled. Fromthedynamicalsystemspointofview,theimportanceoftestingstationarityis singled out by the Doob stochastic shift theorem [16, Chap. X]. It states that stationarity, in the strong sense, of signals u and/or y is necessary and su fficient for existence of an autonomous dynamics (W, S, m, τ) that generates u and/or y. 79 We should mention that we do not consider the notion of weak stationarity, which canbe foundinthe literature onlineartime seriesanalysis. Weakstationarityonlyrequires statistical quantities up to second order to be constant. In the nonlinear case, "weak stationarity" is certainly inadequate. 4.2.1 Testing Methods Asafirstrequirementfortestingstationarity,thetimeseriesshouldcoverastretch of time which is much longer than the longest characteristic time scale that is relevant for the identification of the system. The longest relevant time scale can be estimated as the inverse of the lowest frequency which still contains a significant fraction of the total power of the signal. The statistically most stable quantities are the mean and the variance. Other stable quantities such as spectral components, correlations or nonlinear statistics may be needed to detect less obvious non-stationarity such as a parameter drift which results in no visible drift in the mean or the distribution of values. Linear correlations and the spectrum may also be una ffected. Only the nonlinear dynamical relations and transition probabilities change appreciably. The non-stationarity also can be detected by comparing prediction errors with respect to a nonlinear model for di fferent sections of the data set. Power Spectral Density By the method of power spectral density, if there is too much power in the low frequencies, the time series must be considered "non-stationary," since the corresponding Fourier modes have only very few oscillations during the observation time. By Wiener- 80 Khinchin Theorem, the Fourier transform of the autocorrelation function equals the power spectrum. By Parseval’s theorem: P = N P n=1 | r n | 2 = N/2 P k= − N/2 ¯ ¯ ¯ ˜ R k ¯ ¯ ¯ 2 ˜ R k := 1 √ N N P n=1 r n e 2πik n/ N Recall that, dueto the linearity of theFouriertransform, thepowerspectrum P ( s) k of an additive superposition of a signal x and an independent noise η is the superposition of the signal P ( x) k and the noise P ( κ) k power, viz., P ( s) k = P ( x) k + P ( η) k As an example, the white noise, which is a prototype for a stationary process, has equalpowerinallfrequencybins. Hence,wehavetoexpectarealproblemonlyifthepower spectrumisdivergingaroundzero,suchas f − α withpositive α. TheBrownianmotionhasa f −2 spectrumandisknownnottoberecurrentintwoandmoredimensions. Thedecisionof whether data are stationary or not is much more di fficult if the power spectrum approaches anonzeroconstant for small frequencies whichis muchlargerthanthe background. We can only suggest that this be interpreted as a warning. Determinist chaotic signals may also have sharp spectral lines but even in the absence of noise, there will be a continuous part of the spectrum. This is an immediate consequence of the exponentially decaying autocorrelation function. Space-Time Separation Plot The embedding vectors at successive times are often also close in phase space due to the continuous time evolution. We call correlations of this kind Temporal. 81 The most important temporal correlations are caused by the fact that data close in time are also close in space, a fact which is not only true for purely deterministic systems but also for many stochastically driven processes. Temporal correlations are present in almost every data set. Virtually every process we find in nature is not stationary, since its parameters are depending on time, even if these dependences may be very weak. The question is, how to detect these temporal correlations and how to estimate a safe value of the correlation time t min . The problem was solved by Provenzale (1992) [78] by introducing the "space-time separation plot." The idea is that in the presence of temporal correlations the probability that a given pair of points has a distance smaller than does not only depend on but also on the time that has elapsed between the two measurements. Space-time separation plot [78] is a technique to detect temporal correlation. If { s n : n=0 ,1 ,2,...} is a time series and s i = { s i ,s i −1 ,...,s i − m+1 } is the delay embedding of dimension m, temporal correlation refers to pairs ( i, j) such that s i , s j are within an -neighborhood. Temporal correlations are present in almost every data set. Temporal correlation can be detected by plotting the number of pairs (i, j) within tolerance as a function of two variables, the time separation ∆ t and the spatial distance . N( ∆ t, )= N m −1 P i=1 j= i+1+ n ∆ t Θ( −k s i − s j k) where Θis the Heaviside step function, N m isthesizeoftheembeddedtimeseries,and n ∆ t is defined as the number of time steps separating s i and s j . Practically, we create for each time separation ∆ t an accumulated histogram of spatial distances.Thenwenormalize it properly and plot contour lines for a fraction of 1 /10, 2 /10, ···, of pairs are found. Temporal correlation are present as long as the contour curves do not saturate. For the 82 processes with significant power in low frequencies, such as 1 /f noise or Brownian motion, thecontourcurveswill neversaturate. Inthiscase, all pointsinthedatasetaretemporally correlated and there is no way of determining an “attractor” dimension from the sample. A similar situation arises if the data set is too short. This deficiency of the data can also be seen as a kind of non-stationarity of the measurement. For non-recurrent processes a more suitable characterization of such processes is the analysis in terms of Hurst exponents. The exponent H, Hurst parameter, reflects temporal correlations, as it can be found in the power spectrum. Indeed also the power spectrum of such signals shows a power law decay towards large frequencies given by β = 1+2 H. Quasi-Poincaré Meta-Recurrence Plot Non-stationarity through time dependent system parameters can cause a lack of recurrence. We define M ij = Θ( −| s i − s j |) where M ij is a symmetric matrix of 0’s and 1’s, and s i and s j are the embedded vectors. Θ is the heaviside function. By choosing a distance,we find the points in the embedded space which are within an -neighborhood. The recurrence plot is the plot of the matrix M ij . It is a manifestation of system’s dynamics. If the plot contains many isolated dots, it means that there is coincidental closeness and hence a strong noise component in the data, oraninsu fficientembeddingdimension. Ifitismostlyscattereddots,thenthedeterministic component is absent or at least weak. Irregularity of the arrangement of the line segments indicates chaos. 83 Stationarity of the full time series requires that the density of line segments be uniforminthe i − j plane. Thetemporalspacingbetweenparallellinesinrecurrenceplotcan also be given hints to the existence of unstable periodic orbits inside the chaotic attractor. Meta-Recurrence Plot The following approach is based on the similarity between parts of the time series themselves, rather than the similarity of parameters derived from the time series by local averages. In particular, the (nonlinear) cross-prediction error, that is, the predictability of one segment using another segment as a data base,isevaluated. Thisisveryappropriate in the case where the data possesses strong nonlinear correlations, and it yields meaningful results even for data from stochastic process. In meta-recurrence plot, the segmented time series is used and cross-prediction errors provides the measure of the distance between segments i and j. The idea of breaking a data set into parts and considering the fluctuations of a certain statistics computed on the single parts is widespread. Let { s n ,n=1,... ,N} be a time series which is split into adjacent segments of length l,the i-th segment being called S l i ={ s ( i −1) l+1 ,. ..,s il }. Traditionally, a statistics γ i is now computed for each such segment. It is then tested if the sequence { γ i } is constant up to statistical fluctuations. Here we rather utilize a cross-prediction error γ ij = γ( S l i ,S l j ). Statistical testing with nonlinear parameters γ is di fficult because we can assume very little about the statistical properties of γ. Also empirical errors are not expected to be independent which complicates the estimation of the variance of γ.Byusingstatistics γ ij on pairs of segments we increase the number of parameters computed from N/l to (N/l) 2 . 84 It can be argued that we gain largely redundant information for the purpose of statistical testing since the γ ij for di fferent i, j are not expected to be independent. However, we will be able to detect di fferent and more hidden kinds of nonstationarity. We can get a more detailed picture about the nature of the changes and, in particular, locate segments of a nonstationary sequence. We exploit the information contained in the relative statistics γ( S i ,S j ), in addition to that contained in the diagonal terms γ( S i ,S i ). Let X ≡{ x n ,n=1,.. .,N x }and Y ≡{ y n ,n=1,... ,N y } betwo time seriesand m bea small integer denoting an embedding dimension. Fromboth time series we can form embedding vectors { x n ,n = m,. ..,N x − 1} and { y n ,n = m, ...,N y − 1}, respectively, in the same m dimensional phase space, where x n =( x n − m+1 ,. ..,x n ).Furtherwe fixa length scale.Foreach y n we want to make a prediction one step into the future, that is, given y n =( y n − m+1 ,..., y n ) we want to estimate y n+1 ,using,however, X as a database. A locally constant approximation to the dynamics relating x n and x n+1 yields the estimate ˆ y X n+1 = 1 | U X ( y n )| P x ´ n ∈ U X ( y n) x ´ n+1 In the above formula, U X ( y n )= { x ´ n : k x ´ n − y n k < },isan -neighborhood of y n , formed, however, within the set X. ¯ ¯ U X ( y n ) ¯ ¯ denotes the number of elements in that neighborhood. For isolated points with empty neighborhoods we take the sample mean of the segment X as an estimate ˆ y X n+1 .Ifwetake X and Y to be the same but exclude from U Y ( y n ) all 2 m −1 vectors which share components with y n ,then ˆ y Y n+1 is an ordinary out-of-sample prediction of y n+1 . The root mean squared prediction error γ(X, Y) of the sequence Y,given X,isdefined by γ(X, Y)= s 1 N Y − m N Y −1 P n= m (ˆ y X n+1 − y n+1 ) 2 85 For a stationary time series, we expect that γ( S l i ,S l j ) depends on ( i − j) unless the coherence time of the process is longer than l. If there is variability in the sequence on time scales longer than l, be it due to a slow variable or due to a changing parameter, the diagonalterms γ( S l i ,S l i )willbetypicallysmallerthanthosewith i6= j. Notethatingeneral γ(X, Y) 6= γ(Y, X). In particular, if the attractor of Y is embedded within the attractor of X, for example, if Y forms a periodic orbit which is present as an unstable orbit in X, points in Y can be well predicted using X as a database. However, Y does not contain enough information to predict all points in X.Whiletheasymmetryof γ ij can provide valuable insights, it may also be confusing in some cases. One can then use a symmetrized statistic like γ ij + γ ji [91]. 86 Chapter 5 Prelude to Novel Heart Dynamics: Computational Information Theory of Abstract Dynamical Systems 5.1 Introduction Information theory is a branch of applied mathematics and electrical engineering involving the quantification of information. It has been believed to be a subset of commu- nication theory, but it also has fundamental contributions in other areas such as statistical physics (thermodynamics), computer science, statistical inference, and probability and sta- tistics (see Figure 5.1). Informationtheoryisgenerallyconsideredtohavebeenfoundedin1948byClaude Shannon in his seminal work, "A Mathematical Theory of Communication." [94] Generally, 87 Information Theory Limits of Communication Theory (Communication Theory) Thermodynamics Quantum Information Theory (Physics) Kolmogorov Complexity (Computer Science) Kelly Gambling (Economics) Inequalities (Mathematics) Fisher Information (Statistics) Limit Theorems (Probability Theory) Figure 5.1: Relationship between Information Theory and the other fields. the Shannon’s theory states that reliable communication is possible over noisy channels provided that the rate of communication is below a certain threshold called the channel capacity. Later on in the 1960’s, Kolmogorov, Chaitin, and Solomono ff presented the idea that complexity is the minimal description length of a string of data and introduced Kol- mogorov complexity theory, also known as "algorithmic information" theory. Gratifyingly, Kolmogorovcomplexity K isapproximatelyequaltotheShannonentropy H ifthesequence is drawn at random from a distribution that has entropy H. Boththeories aimatproviding ameans formeasuring"information." In the Shan- non approach we are interested in the minimum expected number of bits to transmit a message from a random source of known characteristics through an error-free channel. In Kolmogorov complexity we are interested in the minimum number of bits from which a particular message or file can e ffectively be reconstructed. They also provide a (distinct) 88 notion of mutual information that measures the information that one object gives about another object. In Shannon’s theory, this is the information that one random variable car- ries about another; in Kolmogorov’s theory ("algorithmic mutual information"), it is the information one sequence gives about another. In an appropriate setting the former notion can be shown to be the expectation of the latter notion [26]. Following in this chapter, we will discuss the Kolmogorov-Sinai entropy and the canonicalcorrelationanalysisasacomputationaltoolforcalculationofmutualinformation. 5.2 Kolmogorov-Sinai Axiomatization TheKolmogorov-Sinainotionofentropyisameasureof"howbizarre"adynamical systemlooks. Sinceitisdefinedforabstractdynamicalsystems,itappliestobothstochastic and deterministic systems. We utilize the Kolmogorov-Sinai entropy KS( u, y),relativeto (W, S, m, τ), as a measure of chaos. As a measure of chaos in, say y + ,weutilizethe Kolmogorov-Sinai entropy rate relative to (W, S, m, τ) KS( y + )= −sup Q lim m →∞ 1 m m −1 P i=0 ,j ∈ J m( T − i Q j )log m( T − i Q j ) where t j ∈ J Q j is a y + -measurable ( m mod 0) partitioning of W. The mutual Kolmogorov- Sinai information rate is defined as I( u − ,y + )= KS( y + ) − KS( y + | u − ) =sup Q,P lim m →∞ 1 2 m m −1 P i=0,j ∈ J m −1 P k=0 ,l ∈ L m( T − i Q j ∩ T k P l )log m( T − i Q j ∩ T k P l ) m( T − i Q j ) m( T k P l ) where KS( y + | u − ) is the conditional entropy and t j ∈ J P l is a u − -measurable ( m mod 0) partitioning of W. The preceding concepts cannot be confused with those built around the 89 G y u e Figure 5.2: Closed loop system diagram. Shannondi fferentialentropyrateSDE. ItshouldbeemphasizedthattheKolmogorov-Sinai entropy rate is not the same as the Shannon di fferential entropy SDE. The most striking manifestation of this discrepancy is as follows: Theorem 19 ([30]) For a general autonomous system y = A( u), KS( y) ≤ KS( u),whereas for a linear time-invariant stable system y( z)= A( z) u( z) , SDE( y) − SDE( u)= 1 2 π Z + π − π log| A(exp( jθ))| dθ . Using the Bode-like integral [30], it is possible to find a system y( z)= A( z) u( z) such that KS( y) <KS( u), as it should be, while SDE( y) >SDE( u). Example 20 ([30]) Consider the closed loop system in Figure (5.2). The system’s transfer function is given as G( z)= z+0 .1 ( Z −2)( Z −0 .2) By the Discrete Bode Theorem, the open-loop system has one unstable pole at 90 Z=2. which implies that SDE( y) − SDE( u)= 1 2 π Z + π − π log| A(exp( jθ))| dθ =2 hence, SDE( y) >SDE( u) This is a violation of Stratonovich’s conservation of information principle. 5.3 Directed Mutual Information Considerasystemonwhichweobservetwovariables: u, y. Supposethatsomefirst principles dictate that u should be the cause and y the e ffect. Thisleadsustoconjecture a relationship of the form y = A( u),where A,the antegrade system,isa causal operator. In the realm of complex systems where some known first principles might not capture the complete picture, it is imperative to allow for A to be as general as possible. In this spirit, causality cannot be defined by the absence of past input ( u − =0) implying the absence of future output ( y + =0), for this would rule out the case where u “modulates” y.One way out of this dilemma is to interpret causality using the concepts of information theory and more specifically, the mutual information between the past of y and the future of u, I( u + ,y − ). Note that the mutual information is defined relative to an invariant measure m of some autonomous dynamical system (W, S, m, τ). In the realm of complex systems with imperfect measurements, it is unlikely that the mutual information I( u + ,y − )=0, but causality could still be asserted if I( u + ,y − ) << 91 u(z) y(z) A R r(z) n(z) Figure 5.3: A feedbacked system: Antegrade system A and retrograde system R. I( u − ,y + ). However, consider a situation in which, under some conditions, under some mal- function, or during a certain time window, the system exhibits a causality flip, I( u + ,y − ) > I( u − ,y + ). In this case, the relationship between u and y can no longer be modeled as a causal system y = A( u). The problem with the (undirected, symmetric) mutual information is the second of the followings: I( u − ,y + )=0 ⇒ / ∃ feedback y = A( u) I( u − ,y + )6=0 6⇒∃ feedback y = A( u) Thus, unless the intuitive relationship I( u − ,y + ) >> I( u − ,y + ) is interpreted as the right hand side vanishing, not much can be asserted. Along the same line of idea, it is not true that I( u − ,y + ) accounts exclusively for an antegrade relationship u = R( y), unless y is a white noise. Indeed, if neither u nor y are white noises, I( u − ,u + ) > 0 and I( y − ,y + ) > 0; thus I( u − ,y + ) has some information in common with I( y − ,u + ), which is the retrograde relationship. 92 Whatweneedisaquantitydependingon u − and y + thatwouldexclusively indicate the strength of the antegrade relationship y = A( u). Therefore should this quantity be > 0, it would positively indicate a relationship y = A( u). Conversely, applied to y − and u + ,this quantity taking a nonvanishing value would positively indicate a retrograde relationship u = R( y). This quantity is the directed mutual information, introduced by Marko in 1973 [59]. Marko’s definition of directed information was later refined by Massey [60]. The two manifestations of directed mutual information are I( y − ⇒ u + ) and I( u − ⇒ y + ).Weusea notation slightly di fferent from the one of Massey [61], [60], I( u − → y + ), I( y − → u + ) as a mild warning that our definition is not quite the same as that of Massey. They are defined as follows: I( u − ⇒ y + )= I( u − ,y + | y − ) I( y − ⇒ u + )= I( y − ,u + | u − ) Consider the first equation. I( u − ,y + | y − ) is the mutual information between u − and y + , conditioned upon y − . This conditioning indeed removes the spurious correlation in the mutual information: The fact that y − and y + are correlated so that I( u − ,y + ) has some informationincommonwith I( u − ,y − ), which itself contains some information about the future of u and the past of y (after partitioning like (·) − =((·) −− ,(·) −+ )). This directed mutual information is defined as follows: I( u − ⇒ y + )= Z u − ,y + ,y − p( u − ,y + ,y − )log p( u − ,y + | y − ) p( u − | y − ) p( y + | y − ) du − dy + dy − 93 As already noted, the issue of mutual information versus directed mutual infor- mation is the removal of spurious correlations which insidiously bring correlation between the past of y and the future of u, even though we appear to correlate u − and y + .Our proposed I( u − ⇒ y + ) does not completely remove the spurious correlations. We note that even the Massey’s approach does not completely remove the spurious correlations. Within the proposed approach, there is an easy fix to remove the remaining spurious correlations by also conditioning on u + . We would therefore define the directed mutual information as I( u − ⇒ y + )= I ¡ u − ,y + |( u + ,y − ) ¢ In another way, y = A( u) only accounts for I( u − ⇒ y + ), whereas another causal relation u = R( y) accounts for I( u + ⇒ y − ).The retrograde mapping R can schematically be viewed as a feedback around A, as shown in Figure 5.3. This leads to the paradigm where feedback is no longer dictated by tracking or robustness specification, but as a way to exhibit a complex causality relation between two variables. Directed mutual information isacomputationallychallengingquantity, henceforth, weperformanalysesbasedonmutual information, which is simpler computationally. 5.4 Canonical Correlation Analysis Canonical Correlation Analysis (CCA) and mutual information techniques are meant to describe the statistical interface between two signals. The interface between the past and future, or on the other hand between the future and the past can be analyzed by these tools [48]. 94 CCAisasecondmomenttechnique. Thereforeitisnotsuitableinitslinearversion for systems with infinite variance, such as self-similar signals or highly noisy data; however, in its nonlinear version, this is no longer an issue. In the nonlinear CCA, the variance analysis is applied to a nonlinear distortion of the original process, which is restricted to result in a finite variance process [50]. Assumethatthetimeseries { y( k): k =··· , −1 ,0 ,1 ,···} is a centered process, bounded and viewed as weakly stationary process with finite covariance E( y( i) y( j)) = Λ i − j defined over the probably space. If the process is not stationary, we can simply compute z( k)= y( k) − y( k −1), which is usually stationary. The past and the future of the process are defined, respectively as y − ( k)= ( y( k),y( k −1),y( k −2),..., y( k − L+1)) T y + ( k)= ( y( k+1),y( k+2),y( k+3),..., y( k + L)) T where L is the lag. Interrelation between past and the future is a preliminary study of whether a recipe of the form y + = f( y − ) is likely to work. The ability to devise a good model can be gauged from the Kolmogorov-Sinai mutual information between the past and the future [50], [48]. The mutual information between the past y − and the future y + is the amount by which the entropy of the future decreases when the past is given [93]; that is, I( y − ,y + )= h( y + ) − h( y + | y − ) = RR log p( y − ,y + ) p( y − ) p( y + ) p( y − ,y + )dy − dy + In the above equation, h( y + )istheentropyofthefutureand h( y + | y − ) is the conditional entropy of the future given the past. 95 5.4.1 Linear Canonical Correlation Analysis In the linear version of CCA, the best linear regression between the past and the future data is sought. In this regard and to proceed froma numerical algebra point of view, the covariance of the past and the future are factored as E( y − ( k),y T − ( k)) = C −− = L T − L − E( y + ( k),y T + ( k)) = C ++ = L T + L + (5.1) where L − and L + are Cholesky factorization of the past and future, respectively. C −− and C ++ areToeplitzmatricesandmeasuresofstrengthofthepastandthefuture,respectively. These quantities are used for normalization to get the information interface, independently of the strength of the signals. Therefore the canonical correlation is defined as Γ( y − ,y + )= L − T − E( y −( k),y T − ( k)) L −1 + which asymptotically is a Hankel matrix in the case of large data records. The Singular Value Decomposition (SVD) of the canonical correlation matrix is given as Γ( y − ( k),y + ( k)) = U Σ V T (5.2) where U and V are orthogonal matrices and Σ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ σ 1 ··· 0 . . . . . . . . . 0 ··· σ L ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , 1 ≥ σ 1 ≥··· ≥ σ L ≥ 0 96 The σ i ’s are called canonical correlation coe fficients (CCC’s) and are invariant under the nonsingular linear transformations of the past and the future. For the Gaussian processes it is well known that I( y − ,y + )= − 1 2 log(det( I − Σ 2 )) (5.3) which is the maximum information can be achieved. It is claimed that there are only a restricted number D ≤ L of significant CCC’s, grouped in Σ 1 . Hence, matrices Σ, U,and V are partitioned as [50] U = ⎛ ⎜ ⎜ ⎝ U 1 U 2 ⎞ ⎟ ⎟ ⎠ , Σ = ⎛ ⎜ ⎜ ⎝ Σ 1 0 00 ⎞ ⎟ ⎟ ⎠ ,V = ⎛ ⎜ ⎜ ⎝ V 1 V 2 ⎞ ⎟ ⎟ ⎠ (5.4) The canonical past and the canonical future are defined as ¯ y − ( k)= U 1 L − T − y − ( k) ¯ y + ( k)= V 1 L −1 + y + ( k) State Space Model Thestateisdefinedastheminimumcollectionofpastmeasurablerandomvariables necessary to predict the future, that is, E(¯ y + ( k)|¯ y − ( k)). The state space model can be defined as x( k+1)= Fx( k)+ w( k) (5.5) where w( k) is a white noise and x( k)= Σ 1 ¯ y − ( k). F is the regression matrix of x( k+1) on x( k) and for the best prediction it is defined as F = E © x( k+1) x T ( k) ª¡ E © x( k) x T ( k) ª¢ −1 97 and the residual error is E n [ x( k+1) − Fx( k)][ x( k+1) − Fx( k)] T o ∼ ( I − Σ 2 1 ) 5.4.2 Nonlinear Canonical Correlation Analysis In nonlinear canonical correlation, a nonlinear processing is done on the past and thefuture. Thisisthecaseinwhichthemaximumpossiblemutualinformationisattempted to be reached. In general (non-Gaussian setup), sup f,g I( f( y − ),g( y + )) ≤ I( y − ,y + ) where f, g : < L → < L are measurable, bijective functions such that E( f)= E( g)=0, E(ff T ) < ∞× I,where I is the identity matrix. Equality is achieved if and only if f( y − ), g( y + )canbemadejointlyGaussian[104]. Inthiscasethelinearestimation ˆ g( y + )= Af( y − ) is optimum. In fact components of f( y − ), g( y + ) can be expressed as linear combination of functions p j ( y − ), q j ( y + ), j=1 ,2, ..., such that E − ( p j )= E + ( q j )=0, E − ( pp T ) < ∞× I and E + (qq T ) < ∞× I, and forming basis of the Lebesgue spaces of zero mean measurable functions such that E − ff T < ∞ and E + g T g< ∞, respectively [50]. For L< ∞,thesame procedure as the linear case is followed and canonical correlation matrix Γ( p( y − ),q( y + )) is computed. The motivation for nonlinear processing of the data is to gauge in comparison with the full dimensional linear case, how much increase in CCC’s is gained by going to the nonlinear analysis. Observation of increased information indicates nonlinearity of the process. For jointly Gaussian functions f( y − ), g( y + ), it can be found that [50] min A E ° ° g( y + ) − Af( y − ) ° ° 2 C −1 ++ = L − Trace( Γ T Γ( f( y − ),g( y + )) 98 for A = Σ f( y − ) ,g( y + ) . It is claimed that [48] f( y − )= U T 1 L − T − p( y − ) g( y + )= V 1 L −1 + q( y + ) where U 1 and V 1 are computed from factorization of the SVD of the nonlinear canonical correlation matrix. The state space model is constructed following the same procedure as in the linear case. In the nonlinear canonical correlation analysis, a nonlinear processing is done on the past of u and the future of y in an attempt to reach the mutual information in a non-Gaussian setup, sup f,g µ − 1 2 logdet ¡ I − Σ 2 ( f( u − n ),g( y + n )) ¢ ¶ ≤ I( u − n ,y + n ) where f, g : < L → < L are measurable, bijective functions such that E( f)= E( g)=0, E(ff T ) < ∞ × I,where I is the identity matrix. Equality is achieved if and only if f( u − n ), g( y + n ) canbemadejointlyGaussian[104]. Consideringthat f( u − n ), g( y + n )arejointly Gaussian, it can be found that [93] min F E ° ° g( y + n ) − Ff( u − n ) ° ° 2 C −1 ++ = L − Trace( Γ T Γ( f( u − n ),g( y + n )) for the regression matrix F = Σ( f( u − n ),g( y + n )). The components of f( u − n ), g( y + n ) can be expressed as linear combination of func- tions p j ( u − n ), q j ( y + n ), j=1 ,2,..., such that E − ( p j )= E + ( q j )=0 E − ( pp T ) < ∞× I E + (qq T ) < ∞× I 99 and forming bases of the Lebesgue spaces of zero-mean measurable functions such that E − ff T < ∞ and E + g T g< ∞, respectively [93]. For L< ∞, the same procedure as in the linear case is followed and the canonical correlation matrix Γ( p( u − n ),q( y + n )) is computed. It is claimed that [93] f( u − n )= U T 1 L − T − p( u − n ),g( y + n )= V 1 L −1 + q( y + n ) where U 1 and V 1 arecomputedfromtheSVDoftheoptimumnonlinearcanonicalcorrelation matrix. The state space model is constructed following the same procedure as in the linear case. The motivation for the nonlinear processing of the data ( u − n ,y + n ) is to gauge the amount of nonlinearity involved in the antegrade system A. 5.4.3 Aggregated Processes For a high resolution process ...,y( k −2),y( k −1),y( k),y( k+1),y( k+2), ... a lower resolution process can be generated as y ( m) k = 1 m m P i=1 y mk − i+1 If the high resolution process looks like the low resolution one, therefore the process is Self- Similar. The question is whether the canonical correlation matrix is changing under the aggregation, Γ( y − ,y + ) ? ≈ Γ( y ( m) − ,y ( m) + ) 100 or in the nonlinear case, Γ( f( y − ),g( y + )) ? ≈ Γ( f ( m) ( y − ),g ( m) ( y + ) . One can think of the invariance of the canonical correlation under linear (non-singular) transformations; however, for the nonlinear case, this might not be always true. 5.5 Model Selection (Lag Window) As is well known, a fundamental di fficulty in statistical analysis is the choice of an appropriate model, estimating and determining the order or dimension of a model. This is a common problem when a statistical model contains many parameters. The main purpose of model evaluation is to understand the observed data. Statistical data modeling is a field of statistical reasoning that seeks to fit models to data without knowing what the "true" model is or might be. Hence, the problem is posed how to choose the "best approximating" modelamongaclassofcompetingmodelswithdi fferentnumbersofparametersbyasuitable model selection criterion given a data set. The general principle is that for a given level of accuracy,asimpleroramoreparsimoniousmodelispreferredtoamorecomplexone,known as Occam’s Razor. Occam’s Razor emphasizes the desirability of selecting the accurate and parsimonious models of reality. Therefore, the best model is the one with the least complexity, or equivalently the highest information gain. 5.5.1 Akaike Information Criteria Theselection of an appropriate approximating model is critical to statistical infer- encefrommanytypesofempiricaldata. Informationtheoryand,inparticulartheKullback- 101 Leibler (Kullback-Leibler1951)"distance"or"information"formsthedeeptheoreticalbasis fordata-basedmodelselection. Akaikein1973foundasimplerelationshipbetweenexpected Kullback-LeiblerinformationandFisher’smaximizedlog-likelihoodfunction. Thisrelation- ship leads to a simple, e ffective, and very general methodology for selecting a parsimonious model for the analysis of empirical data. Akaike, in a very important series of papers, including Akaike (1973, 1974, 1977, and 1981) [71], was perhaps one of the first who laid the foundation of the modern field of statistical data modeling, statistical model identification or evaluation. He developed the information-theoretic, or entropic Akaike Information Criteria (AIC) for the identification of an optimal and a parsimonious model in data analysis from a class of competing models which takes model complexity into account [10]. The major development of AIC lies in the extension of an entropic or information-theoretic interpretation of the method of maximum likelihood. Its introduction is based on the minimization of Kullback-Leibler information quantity. Considerthepredictionofvariable Y fromregressionvariables X i ,i=1,... ,p,via the model θ( Y)= φ( X)+,where is the residual error, θ and φ are real values functions. By means of an appropriate regression method, variance for the best estimate of functions, ˆ σ 2 , and the order of the model (Lag), L, is computed. The Akaike Information Criteria can be found as AI C=log 2 ˆ σ 2 + 2 L N where N is the size of data. Minimization of the AIC leads to the choice of best estimation for functions θ and φ, and the best choice of order of the system [10]. 102 Chapter 6 Dynamical System Theory Analysis of Clinical Data 6.1 Introduction Inthischapter,wedevelopaheartdynamicscasestudyofthestatisticaldynamical systems concepts discussed in the previous chapters. While there are many variables that can be associated with the cardiovascular system, here we only retain two parameters: the RR and QT sequences. Recall that the RR interval is the amount of time between consecutive heart beats. The QT interval is the amount of time between the beginning of the Q wave and the end of the T wave. In a few words, the QT interval is the amount of time the ventricles are functioning during the heart cycle. Under normal conditions, the RR sequence is controlling the QT interval. How- ever, under some pathological condition typically created by Sodium channels blockers, the 103 QT interval is prolonged, which could lead to ventricular fibrillation, a life threatening situation. In the following, we discuss the results of nonlinear statistical analyses of RR and QT intervals as characteristics of heart dynamics. The results are categorized under two main approaches: statistical dynamical and information theoretic analyses. 6.2 Clinical Data Collection Thedataanalyzedinthisthesisisfromrandomized,doubleblind,5-waycrossover, Phase 1 clinical trial. There were 12 subjects included in the data set. Each subject had 5 recordings (one had four). The recordings were baseline, placebo, low-dose, medium dose and high-dose drug consumption in randomized order. All subjects gave informed consent to the pharmaceutical company that allowed the data to be used for research. The original use of the data was to evaluate a Sodium-channel blocker for arrhythmias and QT intervalprolongation. Thedrugisnolongerbeingdeveloped. Datasetsincludethe24-hour measurement of RR and QT intervals. There were also 60 data sets from normal subjects. The results of only one subject is discussed throughout this chapter. 6.3 Dynamical System Approach to Heart Physiology This section and the succeeding section are devoted to the statistical dynamical system analyses and their results. Theheartphysiologyisaverycomplexdynamicalsystem,especiallybecauseofthe influenceoftheAutonomousNervousSystem(ANS).Insteadofrelyingon“firstprinciples,” 104 (Atria SA node) Retrograde System R A Antegrade System (AV node Ventricles) u (RR Sequence) (QT sequence) y Figure 6.1: Feedback as a schematic way to exhibit complex relation. here, we adopt the point of view that the only objective reality about heart dynamics is the experimental data sequence (RR n ,QT n ). A model of the cardiovascular system would be a dynamical system, no matter how abstract, but provided it generates the experimental sequence. 6.3.1 From Time Series to Kolmogorov Abstract Dynamical Systems Goingfromtime-seriestodynamicalsystemsisformalizedunderthefamousDoob stochasticshifttheorem. ThedynamicalsystemconstructedintheDoobtheoremisformal- ized as a Kolmogorov abstract dynamical system, that is, a quadruple (W, S, m, τ),where W is the state-space, S is a Borel field of subsets of W, m is an invariant measure, and τ is a one-parameter family of dynamical shifts, that is, τ k w k 0 = w k 0 + k . The restriction broughtaboutbythisconceptisthatsuchasystemis autonomous, i.e., thedynamical shift τ k : w k 0 7→ w k 0 + k does not depend on k 0 . Having collected the RR n data, it is now treated as a time-driven system, and the resulting QT n ,asshowninFigure6.1,isalsotreatedastime-driven. If( RR, Q T)is(strictly) 105 stationary as a time-driven system, by the Doob stochastic shift theorem [16, Chap. X], there exists a dynamical system (W, S, m, τ) and a measurable mapping (u, y): W → R 2 such that ( u( w n ),y( w n )) = (RR n ,QT n ). In this formulation, experimental evidence indi- catesthattheinducedmeasure u ∗ ( m)fitstheErlangdistribution. In[3],(strict)stationarity of the collected {RR n } was considered and confirmed with varying degree of evidence. 6.3.2 Kolmogorov-Sinai Axiomatization Formally, the system that generates the RR interval is event-driven;morespecif- ically, it is driven by the R-wave event. The most likely scenario is that the R-wave is triggered by the state w( k) of some underlying dynamical system entering some subset S of the state space W. As such, the RR interval is the Poincaré return time (or compound return time) of the state to the subset S of the sate space W,with m, a measure invariant under the dynamical shift τ : W → W, guaranteed to exist under compactness of W by the Krylov-Bogoliubov theorem [68, Chap. VI, Th. 9.05]. Given a subset S ∈ S and w(0) ∈S, the first return time t 1 is defined as min © k>0: τ k ( w(0)) ∈S ª . The return time sequence t n is defined by induction. Under ergodicity, ( φ-)mixing [39] and some generic conditions onS,thesequence t n isstationaryinthesensethat t n isexponentiallydistributedindepen- dently of n. On the other hand, should the dynamics be integrable with vanishing entropy, then t n rather has a power law distribution. If the system is not ergodic with a partitioning W = W 1 ∪ W 2 ( mmod 0)and m( W 1 ∩ W 2 )=0intoinvariantsets, andforasubsetS across the boundary, m(S ∩ W 1 ) > 0, m(S ∩ W 2 ) > 0, a general result states that the distribution of t n is some weighted average of the distribution of t n in W 1 and the distribution of t n in W 2 . 106 The main point of [83] is that, conversely, the distribution of the return time providesinformationontheunderlyingdynamics. Clinicalevidencerevealsthat RR n (= t n ) is Erlang distributed, which in view of [39] implies that the dynamics is mixing. As it will be discussed, it occasionally occurs that the distribution is bimodal, indicating that the dynamics evolutes near the boundary between two invariant regions. 6.3.3 Poincaré Recurrence Poincaré recurrences are able to capture some of the fundamental properties of dynamical systems. It has been proved that the asymptotic distribution of k-fold Poincaré return times is Erlang for a wide class of mixing systems, even if they are not uniformly hyperbolic. We studied the RR intervals in the light of the k-fold Poincaré return time of the heart dynamical system. The main point is that the return time τ is related to the RR interval, so that the distribution of the latter would reveal heart dynamics. We propose a new hypothesis regarding the dynamics of the SA node, which itself generates the RR intervals. We proceed from the reasonable hypothesis that the SA node is a dynamical attractor. The next and bolder hypothesis is that a heartbeat is emitted when the state of the SA attractor has returned k consecutive times to some finite region S, so that the RR interval would be the k-fold Poincaré return time to S. It turns out that the theoretical Erlang distribution of the k-fold return time is a very good fitoftheexperimental RR histograms. More specifically, we developed a new heart dynamics paradigm to formulate a hypothesis that could justify the clinically observed fact that the RR interval is Erlang distributed. An explanation of this phenomenon has already been put forward for atrial 107 fibrillation patients [32]. The hypothesis was that the AV node would transform exponen- tiallydistributedarrivalsfromthemanyectopicsourcesontheatriatoanErlangdistributed process. However, our clinical data supports Erlang distribution on normal subjects. The discrepancy can be explained on the basis that [32] failed to take into consideration the e ffect of digitalis, on which the patients had been put. Our paradigm on the other hand is rooted in recent results on ergodic theory, saying that the k-fold Poincaré return time of a mixing system is Erlang distributed. The “carrier” k times faster than heart rate is hypothesized to be resulting from the synchronization of the pacemaker cells in the SA node. Under the influence of some drugs, the distribution becomes bimodal, which can be explained within the ergodic theory paradigm that the heart operates near the boundary between two attractors. Erlang fit versus multimodal distribution N. Haydn and S. Vaienti have proved that, for a large class of mixing dynamical systems, the return event is a Poisson process [40]; hence the return times are independent and exponentially distributed. Our observation from the experimental data shows that the RR interval for control subjects has an Erlang distribution. We reconcile the two by hypothesizing that the RR interval contains k> 1 return times. Recall that the Erlang distribution is a continuous distribution, and is parameter- ized by two parameters: the shape k, which is an integer, and the rate λ,whichisareal. The Erlang distribution is the distribution of the sum of k independent identically distrib- uted random variables each having an exponential distribution. The probability density 108 function of the Erlang distributed random variable τ is f( τ; k, λ)= λ k τ k −1 e − λτ ( k −1)! for τ> 0 where e is the base of the natural logarithm. The maps τ 1 , τ 2 ,defined over W 1 , W 2 , respectively, satisfy the following condi- tions: τ 1 = τ | W 1 \ W 1 ∩ W 2 ,τ 2 = τ | W 2 \ W 1 ∩ W 2 that is, they coincide with τ except possibly on the zero measure boundary W 1 ∩ W 2 [37]. Thus if W 1 is mixing with a faster return time than W 2 , it can be expected that the left “head” of the weighted distribution would be exponential for the first return time and Erlang for the k-fold return time. Our observation from the experimental data shows that the RR interval for some subjects (mainly normal) has an Erlang distribution. We hypothesize that the RR interval contain k> 1 return times. The RR sequence is weakly correlated, in full agreement with the theory developed by Haydn [40]. Simulation results for Erlang distribution of RR interval are shown in Figure 6.2. The data was collected during 24 hours on a subject under some drug with varying doses. For each subject, five data sets were collected for cases of baseline, placebo, low dose, medium dose and high dose drug. For the results depicted in the Figure 6.2, the subject’s baseline RR interval histogram has a good fit to an Erlang distribution. However, as the drug dosage increases, the distribution of RR interval deviates further from its best Erlang fit. Also,bimodalandmultimodaldistributionsareobservedforlowtohighdrugdoses. The left “head” of the bimodal distribution has a good fit to an Erlang distribution; however, 109 Figure 6.2: Probability distribution and best Erlang fit for a patient under drug consump- tion; solid line: Probability distribution of RR intervals; dot-dash lines: Best Erlang distri- bution fit. the right “tail” of the bimodal distribution does not appear to admit an easy fittoanyof the classical distributions, as shown in Figures 6.3, 6.4. 6.3.4 Chaos The chaotic behavior of the heart dynamics and its complexity are analyzed via Lyapunov exponents, Kolmogorov-Sinai entropy and Kaplan-Yorke (Lyapunov) dimension as follows. 110 400 600 800 1000 1200 1400 1600 0 1 2 3 4 5 6 7 8 x 10 -3 RR f RR (RR) Experimental Prob. Dist. Erlang Dist. Fit; f E (RR) Normal Dist. Fit; f N (RR) a * f E (RR) + b * f N (RR) Figure 6.3: Left/right fit to the experimental probability distribution of RR interval; left tail: Erlang distribution; right tail: Normal distribution. 400 600 800 1000 1200 1400 1600 0 0.5 1 1.5 2 2.5 3 3.5 x 10 -3 RR f RR (RR) Experimental Prob. Dist. Erlang Dist. Fit Weibull Dist. Fit a* f E + b* f W Figure 6.4: Left/right fit to the experimental probability distribution of RR interval; left tail: Erlang distribution; right tail: Weibull distribution. 111 Lyapunov Exponent The Lyapunov exponents as indication of chaos in dynamical system were inves- tigated for the sequence of RR interval among the group of 12 subjects under the drug consumption. The Lyapunov exponents were calculated for embedding dimension equal to 3 as depicted in Figure 6.5; this assumption is part of the phase space investigation of RR interval. Observe the need to compromise, as increasing the dimension of the system re- sults in excessively higher computation time. The outcome of positive Lyapunov exponents among all subjects ascertains that the dynamics of the heart rate is chaotic. Only one positive Lyapunov exponent was observed for each subject. Kolmogorov-Sinai Entropy The Kolmogorov-Sinai entropy of the subject is approximated by the summation over the positive Lyapunov exponents, using (4.5); there is only one positive Lyapunov exponent found for each subject. The results are tabulated in Table 6.1. Kaplan-Yorke (Lyapunov) Dimension The Lyapunov dimension is a measure of the complexity of the system. We cal- culate it using (4.6). It is about 2.6 (in average) for the selected subject. The noninteger Lyapunov dimension indicates the existence of chaos in RR interval and consequently in heart dynamics. The results are tabulated in Table 6.1. 112 (a) (b) (c) (d) (e) Figure 6.5: Lyapunov exponents: (a) Baseline, (b) Placebo, (c) Low Dose, (d) Medium Dose, (e) High Dose. 113 Table 6.1: Lyapunov exponents, Kolmogorov-Sinai entropy and Kaplan-Yorke dimension. Index (Patient No.) Lyapunov Exponents KS Entropy (estimated) KY Dimension (estimated) 52984 0.4248038 -0.01747596 -0.6580406 0.4248038 2.619001 38714 0.4291282 -0.02000201 -0.6633723 0.4291282 2.616737 50995 0.4239012 -0.0201440 -0.6744179 0.4239012 2.598675 72139 0.4100861 -0.0234355 -0.6571032 0.4100861 2.588417 94923 0.4102783 -0.02448697 -0.6585543 0.4102783 2.585815 6.3.5 Statistics of Heart Dynamics Stationarity In the following, we discuss the stationarity of RR interval, using the previously discussed techniques. It is desirable that the RR data be nonstationary, meaning that the heart dynamics can operate over a wider range of variation. Poincaré Meta-Recurrence Stationarity Figure 6.6 shows the meta-recurrence plot of the control subject, notunder any drug. Figures 6.7-6.10show thesame meta-recurrence plots, but for a subject under placebo, low, medium and high doses of a drug under clinical trial. Clearly, the drug a ffects that cardiovascular system by making RR more stationary. 114 (a) (b) Figure 6.6: Poincaré meta recurrence stationarity test for the control subject: (a) Auto prediction error; (b) Cross prediction error. (a) (b) Figure 6.7: Poincaré meta recurrence stationarity test for the placebo subject: (a) Auto prediction error; (b) Cross prediction error. 115 (a) (b) Figure 6.8: Poincaré meta recurrence stationarity test for subject under low drug doage: (a) Auto prediction error; (b) Cross prediction error. (a) (b) Figure 6.9: Poincaré meta recurrence stationarity test for the subject under medium drug dosage: (a) Auto prediction error; (b) Cross prediction error. 116 (a) (b) Figure 6.10: Poincaré meta recurrence stationarity test for the subject under high drug dosage: (a) Auto prediction error; (b) Cross prediction error. Space Time Separation Plot The space-time separation plots for the selected subject are shown in Figures 6.11-6.15. The stationarity test was performed with two embedding dimensions: m=3 and m=8. The results reveal the non-stationarity of RR interval. Scatter Plot Scatter plot analysis is traditional quantitative-visual technique whereby the shape of the plot is categorized into functional classes that indicate the degree of heart failureinasubject[55]. Theplotprovidessummaryinformationaswellasdetailedbeat-to- beat information on the behavior of the heart [17]. The RR interval scatter plot typically appears as an elongated cloud of points oriented along the line-of-identity. The dispersion ofpointsperpendiculartotheline-of-identityreflectsthelevelofshorttermvariability. The dispersion of points along the line-of-identity is thought to indicate the level of long-term variability [79]. The standard deviation of the points perpendicular to the line-of-identity denoted by SD1 describes short-term variability which is mainly caused by respiratory 117 (a) (b) Figure 6.11: Space-time separation plot stationarity test for the control subject: (a) Em- bedding dimension m=3; (b) Embedding dimension m=8. (a) (b) Figure 6.12: Space-time separation plot stationarity test for the placebo subject: (a) Em- bedding dimension m=3; (b) Embedding dimension m=8. 118 (a) (b) Figure 6.13: Space-time separation plot stationarity test for the subject under low drug dosage: (a) Embedding dimension m=3; (b) Embedding dimension m=8. (a) (b) Figure6.14: Space-timeseparationplotstationaritytestforthesubjectundermediumdrug dosage: (a) Embedding dimension m=3; (b) Embedding dimension m=8. 119 (a) (b) Figure 6.15: Space-time separation plot stationarity test for the subject under high drug dosage: (a) Embedding dimension m=3; (b) Embedding dimension m=8. sinus arrhythmia (RSA). The standard deviation along the line-of-identity denoted by SD2 describeslongtermvariability[79]. Thescatterplotof RRinterval(withlagone)ingeneral showstwotypesofdistribution: (I)almostoval(comet)shape;(II)non-oval(torpedo)shape (this is due to the Woo’s sisters; circa 1980). The first category (type I) was observed more amongthenormalsubjects,speciallythosewithanErlangdistributionfor RRinterval. The short time correlation exists for abnormal signals with varying degrees and this correlation is less for normal subjects and for normal subject the long term correlation is larger. The scatter plots of QT − RR intervals were also investigated. The more elliptical shapewereobservedamongthenormalsubjects. Fortheselectedsubject, Figures6.16-6.20 show the scatter plots for RR and QT − RR intervals. Strong relationship between good Erlang fit and type I shape of the scatter plots were observed. Power Spectral Analysis The results of power spectral analysis are depicted in Figure 6.21. They show an almost smooth distribution of power with frequency. 120 (a) (b) Figure 6.16: Scatter plot for the control subject: (a) RR interval; (b) QT − RR intervals. (a) (b) Figure 6.17: Scatter plot for the placebo subject: (a) RR interval; (b) QT − RR intervals. 121 (a) (b) Figure 6.18: Scatter plot for subject under low drug dosage: (a) RR interval; (b) QT − RR intervals. (a) (b) Figure 6.19: Scatter plot for subject under medium drug dosage: (a) RR interval; (b) QT − RR intervals. 122 (a) (b) Figure6.20: Scatterplotforsubjectunderhighdrugdosage: (a) RR interval; (b) QT − RR intervals. Too much power in the low frequencies indicates that the time series must be considered non-stationary. 6.4 Information Theoretic Analyses Simulationresultsfor24-hourandhourlyanalysesofthepast/futuremutualinfor- mation are shown in Figures 6.22-6.26. For each subject, the linear and nonlinear canonical correlation analyses were performed and the mutual informations for QT as past and RR as future and vice versa were computed and compared. Thereareseveraldi fferentmethodsforchoosingtheembeddingdimension, thelag L [97]. In theory, any value of the lag, L, is acceptable, but the shape of the embedding manifold will depend critically on the choice of L and it is wise to select a value of L that separates the data as much as possible. Typically, for the maximally spread data phase space, the vector field will be maximally smooth; that is, with minimum possible 123 (a) (b) (c) (d) (e) Figure 6.21: Power spectral plots: (a) Control subject; (b) Placebo subject; (c) Subject under low drug dosage; (d) Subject under medium drug dosage; (e) Subject under high drug dosage. 124 sharp changes in directions across the flow. Some of the methods for choosing the lag are the Akaike Information Criterion or the first zero of the autocorrelation function. For the latter, some authors recommend choosing the lag such that the autocorrelation has first dropped below a threshold (the so-called decorrelation time). For this analysis, the lag length was chosen to be convenient for running the simulations in a proper time. Hence, the values of the lag were selected from the set L ∈{25 ,50 ,100 ,500 ,1000 ,2000}. Di fferentlagshavethee ffectoflookingattheheartthroughdi fferentwindows. The RR and QT data are long empirical data sets, each containing about 75000 measurements inaverage. Thesamplingrateof RRand QT datasetsis0.01minuteinaverage. Therefore, a lag length of, for example, 500 is equivalent to 5 minutes of data recording. 6.4.1 Canonical Correlation Analysis In this section, we categorize the canonical correlation analysis mainly in four parts: 1) Linear Canonical Correlation Analysis (CCA), 2) Nonlinear Canonical Correlation Analysis, 3) Linear/Nonlinear Canonical Correlation Analysis of the Aggregated Data, 4) Kolmogorov-Sinai (KS) Mutual Information. Ineachcategorytheanalysesaredoneon RR, QT and RR − QT intervals. Weused the mutual information as a measure to investigate the nonlinearity of the RR interval; in caseof (RR − QT)analysis, itsvalueindicateswhichof RRor QT processesisthepredictor process for the other one. 125 For each lag choice, linear and nonlinear canonical correlation analyses were per- formed. The kth data (the first element of the past) was set to be in the middle of the processes in the 24-hour analysis, and the indices of the past and the future were in the interval [ k − L, k + L]. The choice of the lag length is e ffective in CCA, since a change of pattern in both linear and nonlinear cases was observed. Figures6.22-6.26depictthe24-hourandhourlyanalysesofKSmutualinformation fortheselectedsubjectundervaryingdoses. The24-houranalysisofKSmutualinformation, using both linear and nonlinear CCA, reveal a higher correlation for RR and RR/QT intervals under the nonlinear processing of the data. This reveals the nonlinearity in heart dynamics. Strong correlation between RR and QT intervals is observed as the lag length increases over 500 and KS mutual information goes to infinity. Figures 6.24(a)-6.26(a), show that I(QT − ,RR + ) >I( RR − ,QT + ) in 24-hour nonlinear CCA. In order to investigate the e ffect of aggregation on correlation analysis, the ag- gregation step, m, was chosen equal to 4. It is observed that the aggregated data has less canonical correlation in the nonlinear case, and keeps almost the same values in the linear case. Hence, nonlinear processing of aggregated data results in smaller mutual infor- mation than the linear one; a significant reduction of nonlinearity is observed specifically for L=300 ,500. Aggregation also results in significant reduction in canonical correlation coe fficients. Therefore, one may conclude that aggregation decreases the nonlinearity of the system. Hence, one can conclude that Γ( y − ,y + ) ' Γ( y ( m) − ,y ( m) + ), whereas in general Γ( f( y − ),g( y + ))6= Γ( f ( m) ( y − ),g ( m) ( y + ) [2]. 126 Table 6.2: Average of mutual information of RR-QT and QT-RR (in bold). No./Lag 25 50 100 500 52984 16.099 15.315 28.010 28.789 58.636 57.767 Inf Inf 38714 16.261 15.468 29.944 33.227 66.219 62.698 Inf Inf 50995 18.026 15.744 32.099 32.397 62.654 64.756 Inf Inf 72139 15.551 15.454 30.444 30.176 68.613 60.013 Inf Inf 94923 15.568 16.903 29.208 32.349 62.404 61.188 Inf Inf In the nonlinear canonical correlation analysis, the nonlinear operations on the past and the future are chosen to be exponential, e ( − x σ x ) . The rationale for this idea is that the exponential contains all powers. The general interpretation of the plots is that it is desirable to see the I( RR − ,QT + ) >I(QT − ,RR + ) . It appears that L = 500 is most relevant. Indeed, for this lag, the correct relationship is achieved in Figure 6.22, a baseline subject not under the drug. All the other cases are those of subjects under the drug (low, high, and medium dosage for Figures 6.24, 6.25, 6.26, respectively), and for those cases the desired I( RR − ,QT + ) >I( QT − ,RR + ) relation is clearly violated at L = 500. For hourly analysis of data, the average of mutual information for the selected subject has been tabulated in Table 6.2. 127 (a) (b) Figure 6.22: Mutual Kolmogorov-Sinai information for the control subject: (a) 24-Hour variation; (b) Hourly variation. 128 (a) (b) Figure 6.23: Mutual Kolmogorov-Sinai information for the placebo subject: (a) 24-Hour variation; (b) Hourly variation. 129 (a) (b) Figure 6.24: Mutual Kolmogorov-Sinai information for the subject under low drug dosage: (a) 24-Hour variation; (b) Hourly variation. 130 (a) (b) Figure 6.25: Mutual Kolmogorov-Sinai information for the subject under medium drug dosage: (a) 24-Hour variation; (b) Hourly variation. 131 (a) (b) Figure 6.26: Mutual Kolmogorov-Sinai information for the subject under high drug dosage: (a) 24-Hour variation; (b) Hourly variation. 132 Chapter 7 Conclusion Probably the most important point of this thesis is that the dynamical system paradigm that the Poincaré return reveals the fine structure of a system can be successfully applied to heart dynamics. The various cells making up the SA node would make a mixing attractor and their synchronization wouldexplain the Erlangdistribution. Theexplanation forthelattercanbefoundalongtherecentlydevelopedhypothesisthatthe SAnodeconsists of an array of self-exciting cells, of slightly di fferent intrinsic frequencies, f 1 ,f 2 ,... Through the network interactions among the various cells, their firings become synchronized. The nonlinearsynchronizationprocessgeneratesintegermultiplesof f 1 + f 2 +···, whichaccount for the k-fold faster frequency. Our observation from the experimental data shows that the RR interval for some subjects has an Erlang distribution. We reconcile the two by hypothesizing that the RR interval contain k> 1 return times. It appears that a heart beat is emitted after k returns of the state w of the underlying abstract dynamical system to some subset S.Thisis 133 consistent with the faster frequency hypothesis of the synchronization of the SA node cells. The Erlang distribution of RR interval was observed more among the normal subjects. An Erlang distributed RR interval corroborates normal RR versus QT relationship. The subjectsintheclassoflowentropyshowednogoodfitofErlangdistributionfor RRinterval. Bimodal distribution on the other hand is indicative of abnormal RR versus QT dynamics and the dynamical interpretation is that the heart operates near the boundary between a good and a bad attractor. The theoretical fact that if W 1 is mixing with a faster return time than W 2 , it can be expected that the left “head” of the weighted distribution would be exponential for the first return time and Erlang for the k-fold return time, which has been observed experimentally, as shown by the bimodal fit of Figure 6.3. The problem is that the theoretical distribution is relative to a great many initial conditions w ∈ S spreading across the boundary between the two regions. However, in the experimental setup, the average is relative to the time, so that we need some process to allow the system to jump from one invariant region to another. This process is probably the e ffect of the vagal nerve. Nevertheless, the left-sided Erlang fit would reveal an invariant set of normal QT − RR relationship, while the small variance of the right fitting distribution reveals less Heart Rate Variability, hence a pathology. Wealsoshowedthatthehuman’sheartisanonlinearchaoticdynamicalsystemby means of nonlinear transformation of heart rate and mutual information criteria. Nonlinear processing on data results in increment of mutual information, hence indicating that the system has nonlinearity. 134 We have developed an information theoretic technique to determine the causality relation between variables in a complex environment that has enough ergodic properties. The RR / Q T causality relation in heart dynamics has served as the experimental gateway to the concepts, the most important of which is that a mixed causality relation reveals a hidden feedback structure. Under normal conditions, the RR sequence is the pacemaker and should hence control the QT interval. However, under some pathological condition typically created by Sodium channels blockers, the QT interval is prolonged, which could lead to ventricular fibrillation, a life threatening situation. Under prolonged QT interval, the RR sequence does not completely control QT. Thee ffectof RRon QT,andsomefeedbackfrom QT to RR,havebeeninvestigated in [3]. Precisely, we utilized the Kolmogorov-Sinai mutual information to elucidate the antegraderelationshipfrom RRto QT andthefeedbackretrograderelationfrom QT to RR, asshowninFigure6.1. Unravelingsuchrelationshipisanarchetypicalcontrolproblem[102]. OurcorroboratingclinicalstudieshaverevealedacorrelationbetweendeparturefromErlang fit and abnormal RR versus QT relationship [3]. Probablythemostchallengingproblemleftforfurtherresearchistheutilizationof theproposedtechniquetodeterminethecausalityrelationbetweenmechanicaldeformation and action potential in nerve fibers [44]. Based on the observation of the simulation results and the theory of Poincaré return time of dynamical systems, wesummarize the conclusion in the following hypotheses and conjectures. 135 Hypotheses 1. Human’s heart is a nonlinear chaotic dynamical system. 2. SA node is a φ-mixing dynamical attractor. 3. RR interval contains k cycles of the SA oscillator, resulting from the nonlinear syn- chronization of the pacemaker cells in SA node. 4. RR interval would be the Poincaré return time. 5. Mixed causality relation between RR and QT intervals reveals a hidden feedback structure. 6. Fromtheclinicalstudies,clearlythedruga ffectsthatcardiovascularsystembychang- ing the underlying dynamics and the relation between RR and QT intervals. 7. Under the influence of some drugs/pathological conditions, the heart operates near the boundary between two attractors. Conjectures 1. For the normal (healthy) heart, the heart rate follows the Erlang probability distrib- ution function. 2. Heart dynamic operates in the boundary of (two) several attractors, if and only if the heart rate follows the (bi) multimodal probability distribution function. 3. If the KS entropy of the heart rate dynamics is low, then the subject is a strong candidate for having abnormal QT − RR relationship. 136 Bibliography [1] M.Akhtar,A.N.Damato,W.Batsford, J.Ruskin,andJ.B.Ogunkelu. Acomparative analysisofantegradeandretrogradeconductionpatternsinman. Circulation, Journal of the American Heart Association, 52:766—778, November 1975. [2] F. Ariaei. Canonical correlation analysis of the heart dynamics. Technical report, University of Southern California, 2005. [3] F. Ariaei, E. Jonckheere, W. Stuppy, and T. Callahan. Kolmogorov-sinai causality in chaotically intertwined dynamics: A heart rate variability case study. pages 980—985, Malta, Italy, March 2008. The 3rd International Symposium on Communications, Control and Signal Processing. [4] L.Arnold,H.Crauel,andJ.P.Ekmann. Lyapunov Exponents. Springer-Verlag,Berlin Heidelberg, 1991. [5] V.I. Arnold. Ergodic Problems of Classical Mechanics. W.A. Benjamin, New York, 1968. [6] R. Barbieri and E.N. Brown. Analysis of heart beat dynamics by point process adap- tive filtering. IEEE Transaction On Biomedical Engineering, 53(1), January 2006. [7] R.D. Berger. Methodology for automated QT variability measurement, 1996. United States Patent 5560368. [8] R.D. Berger, E.K. Kasper, K.L. Baughman, E. Marban, H. Calkins, and G.F. Tomaselli. Beat-to-beat QT interval variability: Novel evidence for repolarization lability in ischemic and nonischemic dilated cardiomyopathy. Circulation, 96:1557— 1565, 1997. [9] G.D. Birkho ff. Proof of the ergodic theorem. volume 17, pages 656—660. Proceedings of the National Academy of Sciences of the Unites States of America, 1931. [10] H. Bozdogan. Model selection and akaike’s information criterion (AIC): The general theory and its analytical extensions. Psychometrika, 52(3):345—370, September 1987. [11] R. Brown, P. Bryant, and H.D.I. Abarbanel. Computing the lyapunov spectrum of a dynamical system from an observed time series. Physical Review A, 43(6), 1991. 137 [12] H. Bruin, B. Saussol, S. Troubetzkay, and S. Vaienti. Return time statistics via inducing. Ergodic Theory and Dynamical Systems, 23:991—1013, 2003. [13] M. Courtemanche and A. Vinet. Reentry in excitable media. In A. Beuter, L. Glass, M.C. Mackey, and M. S. Titcombe, editors, Nonlinear Dynamics in Physiology and Medicine, Interdisciplinary Applied Mathematics; Mathematical Biology, chapter 7, pages 191—232. Springer, New York, 2003. [14] J.P. Crutchfield, J.D. Farmer, N.H. Packard, and R.S. Shaw. Chaos. Science Ameri- can, 255:38—49, 1986. [15] P.DasandW.H.Schwarz. Solitonsincellmembranes. Phys. Rev. E,51(4):3588—3612, Apr 1995. [16] J.L. Doob. Stochastic Processes. Wiley Classics Library, 1953. [17] B. Eckhardt and D. Yao. Local lyapunov exponents in chaotic systems. Physica D, 65:100—108, 1993. [18] J.P. Eckmann, S.O.Kamphorst, D.Ruelle, andS. Ciliberto. Lyapnovexponents from time series. Physical Review A, 34(6), 1986. [19] J.P.EckmannandD.Ruelle. Ergodic theoryofchaos andstrangeattractors. Reviews of Modern Physics, 57(3):617—656, July 1985. [20] N. El-Sherif. Mechanism of ventricular arrhythmias in the long QT syndrome: On hermeneutics. Journal of Cardiovascular Electrophysiology, 12(8):973—976, 2001. [21] T. Elbert. And physiology: Deterministic chaos in excitable cell assemblies. Physical Reviews, 74(1), 1994. [22] J.O. Fortrat, Y. Yamamoto, and R.L. Hughson. Respiratory influences on non-linear dynamics of heart rate variability in humans. Biol. Cybern., 77:1—10, 1997. [23] A.A. Fossa, T. Wisialowski, A. Magnano, E. Wolfgang, R. Winslow, W. Gorczyca, K. Crimin, and D.L. Raunig. Dynamic beat-to-beat modeling of the QT-RR interval relationship: AnalysisofQTprolongationduringalterationsofautonomicstateversus human ether a-go-go-related gene inhibition. The Journal OF Pharmacology and Experimental Therapeutics, 312(1), 2005. [24] P. Grassberger and I. Procaccia. The strangeness of strange attractors. Physica 9D, (9):189—208, 1983. [25] P. Grassberger and I. Procaccia. Dimensions and entropies of trange attractors from a fluctuating dynamics approach. Dimensions and Entropies of Strange Attractors from a Fluctuating Dynamics Approach, 13:34—54, 1984. [26] P. Grunwald and P. Vitanyi. Shannon information and kolmogorov complexity. IEEE Transactions on Information Theory, 2004. arXiv:cs/0410002v1 [cs.IT]. 138 [27] M.R. Guevara. Dynamicsof excitablecells. InA. Beuter, L. Glass, M.C. Mackey, and M.S. Titcombe, editors, Nonlinear Dynamics in Physiology and Medicine, volume 25, chapter 4. Springer, 2003. [28] A.C. Guyton and J.E. Hall. Textbook of Medical Physiology. Elsevier Saunders, eleventh edition, 2006. [29] H. Haken. At least one lyapunov exponent vanishes if the trajectory of an attractor does not contain a fixed point. Phys. Lett. A, 94:71, 1983. [30] A.Hammad. ControlofChaos.PhDthesis,Dept.ofElectricalEngineering,University of Southern California, 1994. [31] J.R. Hampton. The ECG Made Easy. Churchill Livingstone, sixth edition, 2003. [32] E.HashidaandT.Tasaki. Considerationsonthenatureofirregularityofthesequence of RR intervals and the function of the atrioventricular node in atrial fibrillation in man based on time series analysis. Japanese Heart Journal, 25(5):669—687—, Septem- ber 1984. [33] B. Hasselblatt and A.B. Katok. Handbook of Dynamical Systems. Elsevier, 2002. [34] N. Haydn. The distribution of the first return time for rational maps. Journal of Statistical Physics, 94(5-6):1027—1036, March 1999. [35] N. Haydn. Statistical properties of equilibrium states for rational maps. Ergodic Theory and Dynamical Systems, 20(5):1371—1390, October 2000. [36] N. Haydn, Y. Lacroix, and S. Vaienti. Hitting and return times in ergodic dynamical systems. The Annals of Probability, 33(5):2043—2050, 2005. [37] N. Haydn, E. Lunedei, L. Rossi, G. Turchetti, and S. Vaienti. Multiple returns for some regular and mixing maps. CHAOS, 15, August 2005. [38] N. Haydn, E. Lunedei, and S. Vaienti. Averaged number of visits. CHAOS, 17, September 2007. [39] N. Haydn and S. Vaienti. Fluctuations of the metric entropy for mixing measures. Stochastics and Dynamics, 4(4):595—627, 2004. [40] N. Haydn and S. Vaienti. The limiting distribution and error terms for return times of dynamical systems. Discrete and Continuous Dynamical Systems, 10(3):589—616, April 2004. [41] R. Hegger, H. Kantz, and T. Schreiber. Practical implementation of nonlinear time series methods: The TISEAN package. Chaos, 9:413, 1999. [42] T. Heimburg and A.D. Jackson. On soliton propagation in biomembranes and nerves. Proceedings of the National Academy of Sciences (PNAS), 102(28):9790—9795, July 2005. 139 [43] T. Heimburg and A.D. Jackson. On the action potential as a propagating density pulse and the role of anesthetics, 2006. http://www.citebase.org/abstract?id=oai:arXiv.org:physics/0610117. [44] T. Heimburg and A.D. Jackson. The thermodynamics of general anesthesia. Biophys- ical Journal, 92(9):3159—65, February 2007. [45] M. Hirata, B. Saussol, and S. Vaienti. Statisticsof return times: Ageneral framework and new applications. Communications in Mathematical Physics, 206:33—55, 1999. [46] A. Hodgkin and A. Huxley. A quantitative description of membrane current and its applicationtoconductionandexcitationinnerve. Journal of Physiology,117:500—544, 1952. [47] A.V. Holden. Chaos. Princeton University Press, Princeton, New Jersy, 1986. [48] E. Jonckheere. Chaotic systems. Technical report, University of Southern California, 2006. Lecture Notes. [49] E. Jonckheere, S. Musuvathy, and M. Stefanovic. A biologically inspired synchronization of lumped parameter oscillators through a distributed parame- ter channel. In IFAC Workshop on Control of Distributed Parameter Sys- tems (CDPS), University of Namur (FUNDP) Namur, Belgium, July 2007. http://www.fundp.ac.be/sciences/cdps07/. [50] E. Jonckheere, K. Shah, and S. Bohacek. Dynamic modeling of internet tra fficfor intrusion detection. [51] J.L. Kaplan and J.A. Yorke. Chaotic Behavior of Multidimensional Di fferential Equations, in Functional Di fferential Equations and Approximations of Fixed Points. Springer, Berlin, 1979b. [52] A.KatokandB.Hasselblatt.Introductiontothemoderntheoryofdynamicalsystems. In G.-C. Rota, editor, Encyclopedia of Mathematics and its Applications, volume 54. Cambridge University Press, 1995. [53] T. Kazakevicius, A. Puodziukynas, V. Sileikis, and V. Zabiela. Antegrade and retro- grade conduction properties after complete radiofrequency ablation of the slow atri- oventricular node pathway. Seminars in Cardiology, 12(3):111—113, 2006. [54] A.G. Kouchnirenko. An estimate from abouve for the entropy of a classical system. (in Russian) Dokl. Akad. Nauk. SSSR, 161(1):360—362, 1965. [55] C. Lerma, O. Infante, H. Perez-Grovas, and M.V. Jose. Poincare plot indexes of heart rate variability capture dynamic adaptations after haemodialysis in chronic renal failure patients.Clin.Physiol.andFunc.Im., 23:72—80, 2003. [56] B. London, A. Jeron, J. Zhou, P. Buckett, X. Han, C.F. Mitchell, and G. Koren. Long QT and ventricular arrhythmias in trangenic mice expressing the n terminus 140 and first transmembrane segment of a voltage-gated potassium channel. volume 95, pages 2926—2931. National Academy of Sciences, March 1998. [57] J.D.MacArthurandC.T.MacArthur. Researchnetworkonsocioeconomicstatusand health. 1986. [58] R. Manuca and R. Savit. Stationarity and nonstationarity in time series analysis. Physica D, 99:134—161, 1996. [59] H. Marko. The bidirectional communication theory - a general ization of information theory. IEEE Transactions on Communications, 21(12):1345—1351, December 1973. [60] J.L. Massey. Causality, feedback and directed information. In International Sympo- sium on Information Theory and its Applications, Waikiki, Hawaii, November 1990. [61] J.L.MasseyAnsP.C.Massey. Conservationofmutualanddirectedinformation. IEEE International Symposium on Information Theory, 1:157—158, 2005. [62] R. Ma˜ né. Ergodic Theory and Di fferential Dynamics. Springer-Verlag, 1987. [63] P.E. McSharry, G.D. Cli fford, L. Tarassenko, and L.A. Smith. A dynamical model for generating synthetic electrocardiogram signals. IEEETransactiononBiomedical Engineering, 50(3), March 2003. [64] D.C. Michaels, E.P. Matyas, and J. Jalife. Mechanisms of sinoatrial pacemaker syn- chronization: Anewhypothesis. Circulation Research,61(5):704—714,November1987. [65] H. Mori. Fractal dimensions of chaotic flows of autonomous dissipative systems. Progress of Theoretical Physics, 63(3):1044—1047. [66] R.J. Myerburg and A. Castellanos. Cardiac Arrest and Cardiac Death,volume1, chapter 36, pages 933—974. Saunders Elsevier, eighth edition, 2008. [67] K. Narayanan, R.B. Govindan, and M.S. Gopinathan. Unstable periodic orbits in human cardiac rhythms. Phys Rev E, 57:4594—4603, 1998. [68] V.V. Nemytskii and V.V. Stepanov. Qualitative Theory of Di fferential Equations. Dover, New York, 1989. [69] J.E. OlginandD.P. Zipes. Specific Arrhythmias: Diagnosis and Treatment,volume1, chapter 35, pages 863—931. Saunders Elsevier, eighth edition, 2008. [70] V.I. Oseledec. A multipicative ergodic theorem: Lyapunov characteristic numbers for dynamical systems. Transactions of the Moscow Mathematical Society, 19:197—231, 1968. [71] E.Parzen,K.Tanabe,andG.Kitagawa. Selected Papers of Hirotugu Akaike. Springer Verlag New York, Inc., 1998. [72] Ya.B. Pesin. Equations for the entropy of a geodesic flow on a compact riemannian manifold without conjugate points. Mathematical Notes, 24(4):553—570, 1978. 141 [73] A.G. Petrov. Liquid crystal physics and the physics of living matter. Mol. Cryst. Liq. Cryst., 332:577—584, 1999. [74] M.E. Pflieger and R.E. Greenblatt. Using conditional mutual information to approxi- matecausalityformultivariablephysiologicaltimeseries. IJBEM,2(5):285—288,2005. [75] M. Pollicott and M. Yuri. Dynamical Systems and Ergodic Theory. Cambridge Uni- versity Press, 1998. [76] S.G. Priori, C. Napolitano, and P.J. Schwartz. Genetics of Cardiac Arrhythmias, volume 1, chapter 9, pages 101—109. Saunders Elsevier, eighth edition, 2008. [77] I. Procaccia. Universal properties of dynamically complex systems: The organization of chaos. Nature, 333:618—623, 1988. [78] A. Provenzale, L.A. Smith, R. Vio, and G. Murante. Distinguishing between low- dimensional dynamics and randomness in measures time series. Physica D, 58:31—49, 1992. [79] A.U. Rajendra. Heart rate analysis in normal subjects of various age groups. Bio- Medical Engineering OnLine, 3:24 doi:10.1186/1475-925X-3-24, 2004. [80] R.C. Robinson. An Introduction to Dynamical Systems: Continuous and Discrete. Pearson Prentice Hall, 2004. [81] V.A.Rohlin. Entropyofmetricautomorphism. (inRussian)Dokl.Akad.Nauk.SSSR, 124:980—983, 1959. English summary in Math. Reviews, vol. 21, pp. 387, April 1960. [82] M.T.Rosenstein,J.J.Collins,andC.J. De Luca. A practical method for calculating largest lyapunov exponents from small data sets. Physica D, 65:117—134, 1993. [83] I.Rossi,G.Turchetti,andS.Vaiente. Poincarérecurrencesasatooltoinvestigatethe statistical properties of dynamical systems with integrable and mixing components. Journal of Physics, 7:94—100, 2005. [84] L.Rossi,G.Turchetti,andS.Vaienti. Poincarerecurrencesasatooltoinvestigatethe statistical properties of dynamical systems with integrable and mixing components. Chaotic Transport and Complexity in Fluids and Plasmas, 7:94—100, 2005. [85] M.RubartandD.P.Zipes. GenesisofCardiacArrhythmias: ElectrophysiologicalCon- siderations, volume 1, chapter 31, pages 727—762. Saunders Elsevier, eighth edition, 2008. [86] D. Ruelle. An inequality for the entropy of di fferentiable maps. Bulletin of the Brazilian Mathematical Society, 9(1):83—87, 1978. [87] D. Ruelle. Sensitive dependence on initial conditions and turbulent behavior of dy- namical systems. Annals of the New York Academy of Sciences, 316:408—416, 1979. 142 [88] M.Sakki, J.Kalda, M.Vainu, andM.Laan. Whatdoesmeasurethescalingexponent of the correlation sum in the case of human heart rate? arXiv:physics/0112031, 2(8), January 2003. [89] M. Sano and Y. Sawada. Measurement of the lyapunov spectrum from a chaotic time series.Phys.Rev.Lett., 55:1082, 1985. [90] M.F.Schneider, D.Marsh,B.Kloesgen, andT.Heimburg. Networkformationoflipid membranes: Triggeringstructural transitions by chainmelting. PNAS, 96(25):14312— 14317, December 1999. [91] T. Schreiber. Detecting and analysing nonstationarity in a time series with nonlinear cross-predictions. Physical Review Letters, 78(843), 1997. [92] P.J. Schwartz, S.G. Priori, and G. Napolitano. The Long QT Syndrome, pages 597— 615. WB Saunders, 3rd edition, 2000. [93] K. Shah, E. Jonckheere, and S. Bohacek. Dynamic modeling of internet tra fficfor intrusiondetection. EURASIPJournalon AppliedSignal Processing,2007:1—14,2007. [94] C.E. Shannon. A mathematical theory of communication. Bell System Technical Journal, 27:379—423 and 623—656, 1948. [95] Y. Sinai. Flows with finite entropy. (in Russian) Dokl. Akad. Nauk. SSSR, 125:1200— 1202, 1959. English summary in Math. Review, vol. 21, pp 386-387, April 1960. [96] Y.Sinai. Ontheconceptofentropyforadynamicalsystem. (in Russian) Dokl. Akad. Nauk. SSSR, 124:768—771, 1959. [97] M. Small. Applied Nonlinear Time Series Analysis,volume52of A.WorldScientific, 2005. [98] J.C. Sprott. Chaos and Time Series Analysis. Oxford University Press, 2003. [99] I. Tasaki and P.M. Byrne. Volume expansion of nonmyelinated nerve fibers during impulse conduction. Biophys. J., 57:633—635, March 1990. [100] M.S.Thaler. The Only EKG Book You’ll Ever Need. LippincottWilliamsandWilkins, 4th edition, 2003. [101] C.D. Wagner and P.B. Persson. Chaos in the cardiovascular system: An update. Cardivascular Research, 40:257—264, 1998. [102] J.C.Willems. Thebehavioralapproachtoopenandinterconnectedsystems: Modeling by tearing, zooming, and linking. IEEE Control Systems Magazine, 27(6):46—99, December 2007. [103] A.Wolf,J.B.Swift,H.L.Swinney,andJ.A.Vastano. Determininglyapunovexponents from a time series. Physica D, 16:285—317, 1985. [104] B.F. Wu. Identification and Control of Chaotic Processes—The Kolmogrov-Sini En- tropy Approach. PhD thesis, University of Southern California, 1992. 143 Appendix A Lyapunov Exponent A.1 Algorithm [89] The exponential divergence or convergence of nearby trajectories (Lyapunov ex- ponents) is conceptually the most basic indicator of deterministic chaos. Let us consider an observed trajectory x( t), which can be considered as a solution of a certain dynamical system: ˙ x = F( x) (A.1) defined in a d dimensional phase space. On the other hand, the evolution of a tangent vector in a tangent space at x( t) is represented by linearizing Eq. A.1 ˙ ξ = T( x( t)) ξ (A.2) where T = DF = ∂F /∂x is the Jacobian matrix of F. The solution of the linear non- autonomous Eq. A.2 can be obtained as ξ( t)= A T ξ(0) (A.3) 144 where A T is the linear operator which maps tangent vector ξ(0) to ξ( t).Themeanexpo- nential rate of divergence of the tangent vector ξ is defined as follows: λ( x(0),ξ(0)) = lim t→∞ 1 t ln k ξ( t)k k ξ(0)k (A.4) where k·k denotes a norm with respect to some Riemannian metric. Furthermore, there is a d dimensional basis { e i } of ξ(0),forwhich λ takes values λ i ( x(0)) = λ( x(0),e i ).These can be ordered by their magnitudes λ 1 ≥ λ 2 ≥··· ≥ λ d , and are the spectrum of Lyapunov characteristic exponents. These exponents are independent of x(0) if the system is ergodic. We often have no knowledge of the nonlinear equations of the system which pro- duces the observed time series. However, there is a possibility of estimating a linearized flow map A T of tangent space from the observed data. Let { x j }( j=1 ,2, ...) denote a time series of some physical quantity measured at the discrete time interval ∆ t, i.e., x j = x( t 0 +( j −1) ∆ t). Consider a small ball of radius E centered at the orbital point x j ,and find any set of points { x k i }( i=1 ,2,.. .,N) included in this ball, i.e., © y i ª ={ x k i − x j |k x k i − x j k ≤ ε} (A.5) where y i is the displacement vector between x k i and x j . We used a usual Euclidean norm defined as follows: k wk=( w 2 1 + w 2 2 +···+ w 2 d ) 1 /2 for some vector w =( w 1 ,w 2 , ...,w d ). After the evolution of a time interval τ = m ∆ t the orbital point x j will proceed to x j+ m and neighboring points x k i to { x k i + m }.The displacement vector y i = x k i − x j is thereby mapped to © z i ª ={ x k i + m − x j+ m |k x k i − x j k ≤ ε} (A.6) 145 If the radius ε is small enough for the displacement vectors © y i ª and © z i ª to be regarded as good approximation of tangent vectors in the tangent space, evolution of y i to z i can be represented by some matrix A j ,as z i = A j y i (A.7) The matrix A, is an approximation of the flow map A T at x j in Eq. A.3. Let us proceed to the optimal estimation of the linearized flow map A j from the data sets © y i ª and © z i ª . A plausible procedure for optimal estimation is the least-square-error algorithm, which minimizes the average of the squared error norm between z i and A j y i with respect to all components of the matrix A j as follows: min A j S=min A j 1 N N P i=1 ° ° z i − A j y i ° ° 2 (A.8) Denoting the (k, l) component of matrix A j by a kl ( j) and applying condition A.8, one obtains d× d equations to solve, ∂S/∂ a kl ( j)=0. One will easily obtain the following expression for A j : A j V = C, ( V) kl = 1 N N P i=1 y ik y il , ( C) kl = 1 N N P i=1 z ik y il (A.9) where V and C are d× d matrices, called covariance matrices, and y ik and z ik are the k components of vectors y i and z i , respectively. If N ≥ d and there is no degeneracy, Eq. A.9 has a solution for a kl ( j). Now that we have the variational equation in the tangent space along the experimentally obtained orbit, the Lyapunov exponents can be computed as λ i =lim n →∞ 1 nτ n P j=1 ln ° ° ° A j e j i ° ° ° (A.10) for i=1 ,2 .. .,d,where A j is the solution of Eq.A.9, and n e j i o ( i=1 ,2 .. .,d) is a set of basisvectorsofthetangentspaceat x j . Inthenumericalprocedure, chooseanarbitraryset 146 n e j i o . Operate with the matrix A j on n e j i o , and renormalize A j e j i to have length 1. Using the Gram-Schmidt procedure, maintain mutual orthogonality of the basis. Repeat this procedure for n iterations and compute Eq. A.10. The advantage of the present method is now clear, since we can deal with arbitrary vectors in a tangent space and trace the evolution of these vectors. In this method, these vectors are not restricted to observed data points, in contrast with the conventional method. The feature allows us to compute all exponents to good accuracy with great ease. Themethodwastestedonvariouschaoticdynamicalsystemstoseeifitcanbeused for characterizing experimental systems. For this purpose, a single variable of the model system (e.g., the x coordinate) was treated as experimental data, and then a d dimensional orbitwasreconstructedbyuseofdelaycoordinates,i.e., x i ={ x( iτ), ...,x( iτ+( d −1) t d )}, where t d is the delay time. We set an upper limit to the number N of points included. We proceed as follows: Firstwechooseanorbitalpoint x j and search the points included in the ball from the beginning of the data file { x j }( j=1 ..., M). When the number of points found in the ball exceeds the upper limit we stop the searching and proceed to solve Eq. A.9. In the case that the number does not exceed the upper limit, though the data file is exhausted, if the number satisfies the condition N ≥ d we proceed to solve Eq.A.9, but if the condition is not satisfied we abandon this point x j andgotothenextpoint x j+ m . For the value of the upper limit of N we chose 20 in this paper. It was confirmed that lower values of the limit, e.g., d ≤ N ≤ 5 for the system with d=3, gave similar results. 147 Appendix B Hénon Map The Hénon map is a benchmark discrete-time dynamical system. It is one of the moststudiedexamplesofdynamicalsystemsthatexhibitchaoticbehavior. TheHénonmap takes a point ( x, y) in the plane and maps it to a new point x k+1 = by k +1 − ax 2 k y k+1 = x k Themap depends ontwo parameters, a and b, which for the canonical Hénon map have values of a=1 .4 and b=0 .3. For the canonical values, the Hénon map is chaotic. For other values of a and b, the map may be chaotic, intermittent, or may converge to a periodic orbit. The map was introduced by Michel Hénon as a simplifiedmodelofthe Poincaré section of the Lorenz model. For the canonical map, an initial point of the plane will either approach a set of points known as the Hénon strange attractor, or diverge to infinity. The Hénon attractor is a fractal, smooth in one direction and a Cantor set in another. 148 -1.5 -1 -0.5 0 0.5 1 1.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 Henon Map Figure B.1: Hénon map for a=1 .4,and b=0 .3. ConsidertheclosedloopdiscretefeedbacksysteminterpretationoftheHénonmap showninFigureB.2. H( z), K( z) are transfer functions generated from the Hénon map. U, Y, r,and n are outputs of systems and input noises; respectively. The state-space representation of the Hénon map as follows: ⎛ ⎜ ⎜ ⎝ x k+1 y k+1 ⎞ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎝ 0 b 10 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ x k y k ⎞ ⎟ ⎟ ⎠ + ⎛ ⎜ ⎜ ⎝ 1 − ax 2 k 0 ⎞ ⎟ ⎟ ⎠ Assuming u k =1 − ax 2 k , we can rewrite the above formula as ⎛ ⎜ ⎜ ⎝ x k+1 y k+1 ⎞ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎝ 0 b 10 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ x k y k ⎞ ⎟ ⎟ ⎠ + ⎛ ⎜ ⎜ ⎝ 1 0 ⎞ ⎟ ⎟ ⎠ u k (B.1) wherewedefine the system matrices A = ⎛ ⎜ ⎜ ⎝ 0 b 10 ⎞ ⎟ ⎟ ⎠ , B = ⎛ ⎜ ⎜ ⎝ 1 0 ⎞ ⎟ ⎟ ⎠ , C = µ 10 ¶ ,and D=0. 149 U(z) Y(z) H( z) K r(z) n(z) Figure B.2: Closed-loop discrete feedback system of Hénon map. The transfer function of discrete systems can be defined as the z-transform of the impulse response of the system, H( z), ∞ P n=0 h( z) z − n = D+ ∞ P n=0 ¡ CA n −1 B ¢ z − n = D+ z −1 C ∙ ∞ P n=0 ¡ z −1 A ¢ n ¸ B Using the closed form sum of a matrix geometric series, we obtain: H( z)= D+ C(zI − A) −1 B = 1 z 2 − b The system defined by H( z) is always stable if | b| <1. We generate two time series as the output of the above systems for 50 ,000 points. B.1 Stationarity Test We first test the stationarity of the time series generated by systems H( z) and K, using the introduced methods of Meta-Recurrence and Space-Time Separation plots. Both methods are visual. Since both signals Y and U are generated by computer from a 150 model withconstantparameters, weexpectthemtobestationaryinthesenseofdynamical systems. Meta-Recurrence Plot Figures B.3-B.6 illustrate the variation of the auto-prediction and cross-prediction errors of signals Y and U. The two signals are divided to 50 equal segments. In auto- prediction test, the problem is simplified to an ordinary out-of-sample prediction. Figures B.3andB.5showthevariationoftheauto-predictionerrorsofsignals Y and U,respectively. The variation of the prediction errors is almost constant up to some statistical fluctuations. In the cross-prediction error test, each segment is used as data base to predict the others. Figures B.4 and B.6 illustrate the cross-prediction errors of signals Y and U, respectively. For both signals, the prediction error remains almost constant for all segments. In deter- ministic chaotic systems, this is an invariant property of the system. It is observed that the cross-prediction error for signal U is less than that of Y. Space-Time Separation Plot The space-time separation plots for signals Y and U are shown in Figures B.7 and B.8. Both signals are embedded in 3 −dimensional spaces. We can read the plot as follows: the contour lines shown indicate the distance we have to go to find a given fraction of pairs, depending on their temporal separation ∆ t. Only for the values of ∆ t where the contour lines are flat (at least on average, they may oscillate), temporal correlation has no more e ffect. We must be only concerned about the first 20 −50 time steps, where the level lines are changing significantly. 151 0 10 20 30 40 50 0.2 0.25 0.3 0.35 0.4 Stationarity Test; Henon Map − Y Segments Number Prediction Error Figure B.3: Stationarity test for the signal Y in Hénon map: Auto-prediction error test. 10 20 30 40 50 10 20 30 40 50 Data Base Stationarity Test; Henon Map − Y Predicted 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Figure B.4: Stationarity test for the signal Y in Hénon map: Cross-prediction error test. 152 0 10 20 30 40 50 0.25 0.3 0.35 0.4 Stationarity Test; Henon Map − U Segments Number Prediction Error Figure B.5: Stationarity test for the signal U in Hénon map: Auto-prediction error test. 10 20 30 40 50 5 10 15 20 25 30 35 40 45 50 Data Base Stationarity Test; Henon Map − U Predicted 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Figure B.6: Stationarity test for the signal U in Hénon map: Cross-prediction error test. 153 0 50 100 150 200 250 300 350 400 450 500 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 Separation in Time Separation in Space ( log 2 ε ) Space Time Separation Plot (m = 3) −− Henon Map − Y Figure B.7: Space-time separation plot for Hénon map, signal Y . The space-time separation plot tells us how many pairs have to be excluded to remove the temporal correlation. Temporal correlations are present as long as the contour curves do not saturate. The presence of temporal correlation means that the probability that a pair of points having a distance smaller than ε depends on both ε and the time elapsed between the two points. One can conclude that the flat curves in the space-time separation plots are indication of stationarity of the system. Hence, both signals Y and U are considered as stationary ones. B.2 Chaos-Lyapunov Exponent Figures B.9 and B.10 depict the Lyapunov exponents of signals Y and U,respec- tively. 154 0 100 200 300 400 500 −2.5 −1.5 −0.5 0.5 1.5 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 Separation in Time Separation in Space ( log 2 ε ) Space Time Separation Plot (m = 3) −− Henon Map − U Figure B.8: Space-time separation plot for Hénon map, signal U . The embedding dimension is m=3. Usually the embedding dimension is chosen greaterthantheactualdimensionofthesystemandhence,thespuriousLyapunovexponents may exist. To recognize these spurious exponents, the Lyapunov exponents are computed for both the forward and backward (reverse) time series. The spurious exponents of the forward and backward time series are not equal in magnitude. Figure B.9 illustrates the magnitude of the Lyapunov exponents of signal Y,com- puted in both forward and backward (reverse) iterations by solid and dash-dot lines, re- spectively. Only two of these exponents are equal in magnitude to their reverse time series’ exponents. The Lyapunov exponent λ ' 0 .79 is spurious. Meanwhile, the Lyapunov ex- ponent λ=0 .4 is exactly equivalent to the backward iteration’s. The same analogies are applied to the Lyapunov exponents of signal U, as shown by Figure B.10. The Lyapunov 155 exponent λ ' 0 .56 is in fact a spurious one. Signal U has a Lyapunov exponent λ = −1 .6 that is equivalent to the exponent of the reverse time series. The Hénon map is a two dimensional dissipative iterated map with chaotic solu- tions. The Kaplan-Yorke dimension is calculated as follows D KY = k + k P i=1 λ i /| λ k+1 | where k is the maximum integer such that the sum of k largest exponents is still non- negative. For both Y and U time series this dimension is equivalent to D KY =1 .25.The non-integer Lyapunov dimension is an indication of the chaotic dynamics of both signals Y and U. TheKolmogorov-Sinai entropy of both signals Y and U are approximatedby their positive Lyapunov exponents, which is bounded and positive, KS entropy| U,Y ≤ 0 .4. Since KS entropy of both signals Y and U are equivalent, both antegrade and retrograde systems H( z) and K in Figure B.2 are autonomous. B.3 Kolmogorov-Sinai Mutual Information Figures B.11 and B.12 show the variation of the Kolmogorov-Sinai (KS) mutual information of signals Y and U. Figure B.11 depicts the KS mutual information analysis over the segmented time series for the lag lengths of L=25 ,50 ,100 and 200,whereas, Figure B.12 illustrates the same analysis over the whole data sets. It is observed that the KS mutual information of the past and future of the signals Y and U vary in a tight manner with time, i.e. I( U − ,Y + ) ≈ I( Y − ,U + ). Hence, no really causally dominant system is observed for the feedback system in Figure B.2. 156 0 5 10 15 20 25 30 35 −2 −1.5 −1 −0.5 0 0.5 1 1.5 Iterations Lyapunov Exponent Lyapunov Exponents: Henon Map − Y + : Forward iterations; *: Backward iterations Figure B.9: Lyapunov exponents of Signal Y - Hénon map; embedding dimension m=3. 0 5 10 15 20 25 30 35 40 −2 −1.5 −1 −0.5 0 0.5 1 Iterations Lyapunov Exponent Lyapunov Exponents: Henon Map − U + : Forward iterations; *: Backward iterations Figure B.10: Lyapunov exponents of signal U - Hénon map; embedding dimension m=3. 157 1 11 21 31 41 51 61 71 81 91 100 0 20 40 Kolmogorov-Sinai Mutual Information of U/Y Solid line: Y Past/U Future; Dash-dot Line: U Past/Y Future. 1 11 21 31 41 49 0 20 40 1 6 11 16 21 24 40 50 60 1 3 5 7 9 11 12 80 100 120 Iterations Lag = 25 Lag = 50 Lag = 100 Lag = 200 Figure B.11: Hénon map: Kolmogorov-Sinai mutual information between signals U and Y; segmented time series analysis. 0 200 400 600 800 1,000 0 100 200 300 400 500 600 Lag Kolmogorov-Sinai Mutual Information; Solid Line: Y Past/U Future; Dash-dot Line: U Past/Y Future. ∞ ∞ Figure B.12: Hénon map: Kolmogorov-Sinai mutual information between signals U and Y; full time series analysis.
Abstract (if available)
Abstract
Heart dynamics and heart rate variability, as well as their characterization, have been of interest to researchers for years. Sudden death due to heart failure, which could be the result of heart diseases, drug consumptions, or even emotional shocks, is one of the most important issues in public health. Henceforth, prediction and/or awareness of such conditions has always been a top priority. In the present work, a time domain approach, as opposed to the traditional frequency-domain approach, to heart rate dynamics is introduced. The Electro-Cardio-Gram (ECG) signals, more specifically the RR and QT intervals as characteristics of heart dynamics, are investigated. It is hypothesized that the Erlang distribution of the RR interval is a manifestation of the Poincaré return time phenomenon, and based on clinical results available to us, it is conjectured that for the healthy heart, it is necessary but not sufficient that the RR interval would be the k-fold Poincaré return time. The higher frequency, k times that of normal heartbeat, is hypothesized to be related to the synchronization of the array of pacemaker cells in the sinoatrial (SA) node. In this work, chaotic behavior of heart dynamics is observed by means of such nonlinear techniques such as calculation of Lyapunov exponents, Kolmogorov-Sinai (KS) entropy, and Kaplan-Yorke (K-Y) dimensions. These techniques measure chaos, the amount of uncertainty and complexity of the system, respectively. The existence of positive Lyapunov exponent and non-integer K-Y dimension in the heart rate confirms its chaotic dynamics. Heart dynamics is also analyzed via statistical techniques such as linear and nonlinear canonical correlation analyses and mutual KS information. These techniques are implemented to unravel the at point subtle relationship between RR and QT intervals.
Linked assets
University of Southern California Dissertations and Theses
Asset Metadata
Creator
Ariaei, Fariba
(author)
Core Title
Dynamical system approach to heart rate variability: QT versus RR interval
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Electrical Engineering
Publication Date
12/16/2010
Defense Date
09/08/2008
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
Chaos,dynamical systems,ergodic theory,heart rate variability,information theory,Kolmogorov-Sinai entropy,OAI-PMH Harvest
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Jonckheere, Edmond A. (
committee chair
), Hsiai, Tzung K. (
committee member
), Nayak, Krishna S. (
committee member
)
Creator Email
ariaei@usc.edu,fariba.ariaei@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-m1929
Unique identifier
UC1125888
Identifier
etd-Ariaei-2542 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-148775 (legacy record id),usctheses-m1929 (legacy record id)
Legacy Identifier
etd-Ariaei-2542.pdf
Dmrecord
148775
Document Type
Dissertation
Rights
Ariaei, Fariba
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Repository Name
Libraries, University of Southern California
Repository Location
Los Angeles, California
Repository Email
cisadmin@lib.usc.edu
Tags
dynamical systems
ergodic theory
heart rate variability
information theory
Kolmogorov-Sinai entropy