非線性時間序列分析沮喪引起的心臟疾病

非線性時間序列分析沮喪引起的心臟疾病(第 3 年)

研究成果報告(完整版)

計 畫 類 別 ： 個別型

計 畫 編 號 ： NSC 95-2112-M-004-001-MY3

執 行 期 間 ： 97 年 08 月 01 日至 98 年 12 月 31 日

執 行 單 位 ： 國立政治大學應用物理研究所

計 畫 主 持 人 ： 蕭又新

計畫參與人員： 碩士班研究生-兼任助理人員：黃俊傑

碩士班研究生-兼任助理人員：金宛貞

碩士班研究生-兼任助理人員：謝佳宏

碩士班研究生-兼任助理人員：李恩慈

碩士班研究生-兼任助理人員：蔡羽青

碩士班研究生-兼任助理人員：方玠人

碩士班研究生-兼任助理人員：張雁茹

博士班研究生-兼任助理人員：吳宜家

處 理 方 式 ： 本計畫可公開查詢

中 華 民 國 99 年 03 月 29 日

5

死亡人數為 11441 人，高居第三位，不僅如此，二十年來心臟疾病一直都在十大

死因「榜上有名」，總是佔據著第三、四位，並有向上攀升的趨勢，這樣的訊息

顯示出若心臟疾病的問題不盡快獲得改善，將繼續主宰著我們生命的健康與否。

除此之外，近年來由研究學者提出憂鬱沮喪的心情會導致心臟疾病已引起廣泛的

注意與討論，例如在「第三屆國際心臟血管病大會」發表的一項美日合作研究結

果顯示，憂鬱沮喪的心情可能也是促使心血管病患發作的隱形殺手（如圖一）。

反之，當病人接受冠狀動脈繞道手術時，醫師在手術前多半會告訴病人與家屬，

有些病人在手術之後會產生抑鬱症狀，但若想知道哪些人會有症狀，嚴重情形如

何，就很少有人能夠回答了，因為相關的研究並不多，不是研究中病人取樣不夠

多，就是難以判斷心臟疾病與抑鬱症乃至於其他致病因子之間的關聯。此外，研

究人員亦發現在手術前後具有中度到重度憂鬱的病患，或者輕到重度的憂鬱情形

持續達六個月以上，病患的死亡機率比沒有憂鬱症狀的病患高了兩倍，而像這樣

從生理影響心理及從心理影響生理反應已成為新興的研究主題。

圖一

在醫學上已將心臟疾病的問題針對不同的情況有了特定的定義，也有一些方法作

為心臟疾病的診斷及治療，而最常見的就是心電圖的診斷了【2】

；其他領域方面，

生物研究者也做了許多動物的電生理量測【3】，以期能有效對應於人體心臟系

統。也正因深入研究後，接踵而來的卻是隱藏在心臟疾病問題之後的複雜心臟系

統問題。此乃因心跳為人類最基本且最靈敏的生命現象，其複雜的行為，至今仍

無合理解釋。其原因在於心臟的跳動是由節律點、baro-及 chemo 的感覺受體、

體液，以及中樞神經系統所組成迴饋系統加以控制。此系統形成的複雜心跳動力

學，可藉由雜亂的心率變化得知。更進一步來講，心跳、血壓及呼吸皆屬於心血

管系統的一部份，目前人們已從臨床數據得知，此三者的動力行為有同步現象產

生，而此同步行為可藉由測量身體各個不同部位的生理訊號來加以驗證。目前，

6

analysis）

、赫氏指數（Hurst exponent）

、多重碎形（multifractal analysis）

、隨機共

振（stochastic resonance）

、波元分析（wavelet analysis）及混沌理論（chaos theory）

，

這些統計及非線性的方法（統稱非線性時間序列分析法）目前已被廣泛地運用在

各種物理系統，我們整合這些想法並發展出新的分析工具於複雜的心跳訊號，目

前已有數篇文章發表（請參閱附件）

。

本計畫已訓練數名研究生及專題生，了解相關知識，及習得數值計算上的一些技

巧，培養具對心理與生理學有一定認識的跨領域研究人才，可供以後國內各領域

研究的新血注入。除此之外，在整個計畫的執行中，參與研究人員對平行環境的

建立、平行程式的撰寫及系統的管理等在過程中也習得經驗與技術，在資訊爆炸

的時代，相信也可高度應用於各學產業。

參考資料

1. http://www.doh.gov.tw/statistic/index.htm

2. Malcolm S. Thaler 原著，呂嘉陞編譯，“心電圖學必備（第三版）＂，

【合記

圖書出版社，2002】

3. Marcus L. Koller, Mark L. Riccio and Robert F. Gilmour, Jr., Am. J. Physiol. 275

(Heart Circ. Physiol. 44), H1635-H1642 (1998).

4. Jeffrey J. Fox, Mark L. Riccio, Fei Hua, Eberhard Bodenschatz, Robert F. Gilmour,

Jr., Circ Res. 90, 289-296 (2002).

5. Niels F. Otani, Chaos,12, 3, 829-842 (2002).

6. A. V. Panfilov, Phys Rev Lett. 88, 118101 (2002).

Other uses, including reproduction and distribution, or selling or

licensing copies, or posting to personal, institutional or third party

websites are prohibited.

In most cases authors are permitted to post their version of the

article (e.g. in Word or Tex form) to their personal website or

institutional repository. Authors requiring further information

regarding Elsevier’s archiving and manuscript policies are

encouraged to visit:

Does well-harmonized homeostasis exist in heart rate

fluctuations? Time series

analysis and model simulations

Yuo-Hsien Shiau

Graduate Institute of Applied Physics, National Chengchi University, Taipei 11605, Taiwan, Republic of China

a b s t r a c t

a r t i c l e i n f o

Article history:

Received 11 March 2008

Received in revised form 3 October 2008 Accepted 24 November 2008 Keywords: Heart Time series Homeostasis Mathematical models Diagnosis

Analyzing heart rate variability from electrocardiographic recordings has been an important method for assessing cardiovascular autonomic regulation. Researchers have conducted extensive analyses on normal as well as pathological hearts, however, it is still unclear whether increasing or decreasing the complexity of heart rate variability is a characteristic of healthy systems. In this study, wefind the existence of well-harmonized homeostasis in heart ratefluctuations, in particular, the evidence is verified among different individuals including healthy subjects, ICU patients, and one child with brainstem dysfunction. The methodology we used is composed of two parts, in which one is the consideration of reduction of cardiorespiratoryfluctuations inherited in the original R–R intervals and the other is based upon the concept of nonlinear dynamics to construct the low-dimensional trajectory in the angle plot. The cross-correlation measure between the theoretical angle map and the numerically derived angle trajectory is used to separate recovery (0.73 ± 0.13) from deterioration (0.25 ± 0.08) of ICU patients. In addition, a simple physiologic (deterministic) model of the interaction between the cardiovascular system and baroreceptor control of arterial pressure is used to explain why homeostasis can exist in heart ratefluctuations. Our study provides a potential link between the clinical data and circulatory system.

© 2008 Elsevier B.V. All rights reserved.

1. Introduction

It is known that both an epilepsy seizure and a heart attack may be considered as dynamic diseases used to describe diseases highlighted by a change in normal body rhythms. What are normal body rhythms? Conceptually, they are a healthy body's simple rhythms. From this point of view, the different body parts will tend towards homeostasis, and the interrelated systems reach a balance or have simple occasional behavior. Nevertheless, researchers suggested that chaos is the natural way to join different situations in the body.Goldberger et al. (1994)

have conducted an extensive analysis of normal and pathological hearts, and claim that chaos provides the body with theflexibility to respond to different stimuli. Thus, healthy systems want to exhibit chaotic/complicatedfluctuations. More recently,Andrés et al. (2006)

provided a study on cardiac dynamics, and they found that the existence of premature ventricular contraction increases the embed-ding dimension of heart rate variability (HRV). In addition, patients with congestive heart failure also show an increase in the embedding dimension of HRV. Therefore, it is still unclear whether increasing or decreasing the complexity of HRV is a characteristic of healthy systems.

Analyzing HRV from electrocardiographic (ECG) recordings has been an important method for assessing cardiovascular autonomic

regulation (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996). The most commonly used prognostic HRV index has been the standard deviation in N_{–N intervals (SDNN; one of linear methods), which is} analyzed over a 24-h period for risk stratification. Spectral analysis of HRV allows assessment of frequency-specific fluctuations in heart rate and provides prognostic information beyond the SDNN measure. Although all the measures of HRV differ in their manner of computation and analysis, these methods are fundamentally based on moment statistics and describe the magnitude of HRV. It is therefore not surprising that SDNN and spectral analysis all have a relatively close mutual correlation, and that there are only minor differences in prognostic power between them. More recently, nonlinear dynamics has opened new approaches for studying and understanding the characteristics of heart rate (Goldberger, 1996). These methods differ from the above-mentioned measures of HRV in that they are not designed to assess the magnitude of variability. Rather, they estimate the correlation properties and complexity of HRV and other features in heart rate dynamics that are not uncovered by methods based on variance and mean. One of typical nonlinear measures is Poincaré plot analysis that allows visual and quantitative analysis of instantaneous and continuous R–R interval variability and also provides more powerful prognostic information on patients with heart failure and on arrhythmic risk (Huikuri et al., 1996).

In the clinical data, it is known that sinus rhythm in very sick patients often varies in a regular way that appears much simpler

E-mail address:[email protected].

1566-0702/$– see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.autneu.2008.11.017

Contents lists available atScienceDirect

Autonomic Neuroscience: Basic and Clinical

j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / a u t n e u

in form than the variability observed in normal, healthy indi-viduals (Goldberger et al., 1990). Besides, deterministic homoclinic orbits illustrated in the Poincaré plot can be directly observed from symptomatic sinus node dysfunction (Bergfeldt and Haga, 2003). Nevertheless, Kantz and Schreiber (1998) concluded that there is no clear evidence for determinism in R–R intervals. Therefore, the conservative working hypothesis would be that the process which governs the initiation of new cardiac cycles is effectively stochastic, superimposed by the regulations of auto-nomic nervous system.

Recently, studies were made to create either deterministic or stochastic models describing the behavior of R–R intervals based on the understanding of the physiological mechanisms underlying their variations (Rosenblum and Kurths, 1995; Seidel and Herzel, 1998; Ivanov et al., 1998).Suder et al. (1998) paid their attention on the evolution of angular component of the R–R interval map in humans undergoing paced respiration at a frequency close to 0.1 Hz, and the deterministic structure was found in the angle map. Their point of view is that the complexity of spontaneous respiratory movements obviously precludes the identification of finite dimensional attractors in heart rate fluctuations within the frequency range of breathing, i.e., so-called respiratory sinus arrhythmia (RSA). Therefore, one way to overcome this problem is to introduce paced breathing. The main result of their paper is that HRV during voluntarily induced slow-paced breathing obeys a one-dimensional, nonlinear law of motion. In addition, the one-dimensional law of motion breaks down for cycle lengths close to that of spontaneous breathing. More recently,

Balocchi et al. (2004)derived/studied RSA from the heartbeat time series using empirical mode decomposition. In their study, the spontaneous respiratory signal (i.e., respirogram) was recorded simultaneously with ECG, in which they found the existence of phase and frequency synchronization between the R–R intervals associated breathing and the respiratory signal itself. Moreover, it can be verified that the phase distribution of the respiratory signal displays a noise-like homogeneous profile in a range from −π to π (Wu et al., 2006). According to these findings, it could be hypothesized that RSA fluctuations include a noisy part which is from the condition of spontaneous breathing.

In the present study, we report the existence of well-harmonized homeostasis in heart ratefluctuations, in particular, the evidence is verified among different individuals including healthy subjects, ICU patients, and one child with brainstem dysfunction. Our analyzing method is a nonlinear measure which allows visual and quantitative analyses of instantaneous and continuous R–R interval variability. Basically, our calculations are based upon the traditional concept of fluctuations in nature, in which a mixture of deterministic and stochastic factors should be considered. In order to extract the deterministic factor embedded in R–R intervals, a noise-reduction method is used tofind the final as well as deterministic R–R intervals displaying a one-dimensional and well-harmonized nonlinear angle map derived from the Poincaré plot, in which the so-called home-ostasis can be realized. In this study, we also provide a clinical research on risk stratification for ICU patients in terms of the angle plot. Surprisingly, the outcome of analyzed 25 patients is well described by the present method. Moreover, a simple physiologic (deterministic) model of the interaction between the cardiovascular system and baroreceptor control of arterial pressure is used to explain the status of ICU patients, which could be the result from the influence of time delay in the human baroreceptor-mediated reflex. Thus, our study provides a potential link between the clinical data and the circulatory system.

The remainder of this paper is organized as follows. Subjects and methods are given in Section 2. Section 3 contains the central part of our paper including analyzing the clinical data, surrogate data analysis, and model simulations. Discussions and concluding remarks are given in Sections 4 and 5, respectively.

2. Subjects and methods 2.1. Subjects

A group of 25 healthy normal subjects were recruited/recorded in the Taipei Veterans General Hospital, where sex (10 females and 15 males), age (30 ± 1 years), and body mass index (21.5 ± 0.4 kg/m2_).

And 25 patients with diseases including those requiring admission to intensive care unit such as myocardial infarction, multiple organ dysfunction syndrome, sepsis, and heart failure were recruited in this study. The study protocol was approved by the local ethics committee and all participants gave their informed consent. The study was conducted according to the principles of the Helsinki declaration.

2.2. Study protocol

For all healthy normal subjects, no alcoholic or caffeine-containing drinks were taken for at least 24 h before the study. The examination was performed in a quiet room during the daytime. Subjects received ECG measurement in the supine position afterfive min rest. During the ECG measurement, subjects were instructed to fully relax, stay awake, and not to speak.

2.3. Data acquirement

The raw R–R intervals were deduced from the adjacent normal sinus beats, which were then transferred to a personal computer and post-processed by a program. The missing intervals (due to extra-systoles) were linearly interpolated and the resulting R–R intervals were resampled at 4 Hz by the Berger method (Berger et al., 1986). All analyzed data were checked by a qualified medical doctor. Therefore, possible artifacts appeared in R–R intervals were deleted.

2.4. Nonlinear noise reduction

Suppose we have a scalar time series {xi}, i = 1,…,T, where the xi

are composed of a clean signal yiwith some noise wiadded, xi= yi+

wi. Thebw2N is called the absolute noise level. The reduction scheme

is to replace each noisy xiby the average value of this coordinate

over points in a suitably chosen neighborhood. The neighborhoods are defined in a phase space reconstructed by delay coordinates. In order to define the neighborhoods, one has to fix positive integers k and l and construct embedding vectorsxi_{= (x}

i− k,…, xi + l) as usual.

Note that past and future coordinates are involved. Further, choose a radius r for the neighborhoods. For each value xifind the set Uirof all

neighbors xj for which ||xj−xi||supbr, i.e., all segments of the

trajectory which are close during a time lasting from k iterations in the past to l iterations in the future. Then, replace the“present” coordinate xiby its mean value inUir,

xcorri = 1 jUr ij ∑ Ur i xj: ð1Þ

The implementation of the algorithm is straightforward. To obtain optimal results it is essential to choose r, the size of the neighbor-hoods, appropriately. In addition, the procedure can be iterated. If one takes the rms of the correction made as a new value for r, r will decrease exponentially with the number of iterations until eventually no neighbors are found for any point and no further correction is made (Schreiber, 1993).

To estimate the embedding dimension d for the reconstruction of the phase space, the false nearest neighbor method proposed by

Kennel et al. (1992)is used to determine the lower bound of the embedding dimension dl(Hegger and Kantz, 1999). The dlvalue for

within the range 7≤dl≤15. Thus, in the present study we choose d=21

(i.e., l = 10 and k = 10) as the embedding dimension, which is large enough to reconstruct the phase space for all analyzed subjects. In addition, the initial r values we used are 50 ms.

2.5. Angle map

Afirst attempt at identifying the underlying process generating the characteristic pattern of HRV is to plot the same data as a two-dimensional scatter (Poincaré) plot of successive R–R intervals. However, no sharply defined pattern is revealed in this plot; the points scatter around an elliptical structure. All R–R intervals of healthy subjects contain this feature, but they do not show a simple nonrandom structure. The elliptical structure reflects the long-term periodicity of R–R intervals under the control of the autonomic nervous system.

As a next step, we focus on the angular motion of the points in the scatter plot. This is done by introducing polar coordinates and by neglecting the_{fluctuations of the radii. The transformation needs the} definition of a center. The center is defined as the mean values of R–R intervals (RRa). For every pair of successive R–R intervals, an angle ϕ is

calculated. The angles are defined to vary between −π to π. n= arctan

RRn + 1−RRa

RRn−RRa

 

ð2Þ whereϕnand RRnrepresent the angle and the R–R interval at time

step n, respectively. It is deserved to note that testing for the angle map using surrogate data to preserve both the power spectrum and the histogram had been done in the work ofSuder et al. (1998). In the following we will give different examples to illustrate the determi-nistic one-dimensional angle map can be derived under different situations.

2.6. Surrogate data analysis

The purpose of surrogate data is to test for any nonlinearity in the original data (Theiler et al., 1992). Surrogate signal is produced by phase randomizing the original data. It has similar spectral properties as of the given data. The surrogate data sequence has the same mean, the same variance, the same autocorrelation function and therefore the same power spectrum as the original sequence, but phase relations are destroyed. In the case of data shuf_{fling, the histograms} of the surrogate sequence and the reference sequence are identical. The random phase spectrum is generated by using the method of phase shuf_{fle, where the phase values of the original spectrum are} used in random order.

2.7. Mathematical model

In this study, we use the classic three-element Windkessel model (Westerhof et al., 1971) to simulate HRV, which describes the interaction between the cardiovascular system and baroreceptor control of arterial pressure. On the basis of the Windkessel theory, the dynamic relationship between the arterial pressure P(t) and the cardiac output Q(t) is dP dt−r dQ dt = 1 RC½ðR + rÞQ−P; ð3Þ Q tð Þ =V_Tτð ÞP τð ÞP = V tð−τÞ T tð−τÞ; ð4Þ

where T and V mean, respectively, the period of the cardiac cycle and the stroke volume, andτ is the time delay in the human baroreceptor-mediated re_{flex. R, C, and r are, respectively, corresponding to a} peripheral resistance, a total arterial compliance, and an aortic

characteristic impedance. It is deserved to note that both T and V depend on the arterial pressure P. According to the sigmoidal law, T(P) can be described as

T Pð Þ = Ts+

Tm−Ts

1 +γe− αP=Pn; ð5Þ

where Tsand Tmestablish the shorter and longer cardiac period and

match, respectively, the maximal vasodepressor-induced sympathetic excitation and the maximal pressor-induced vagal activation (Franz, 1969). Pncorresponds to the steady level of arterial pressure.α and γ

are twofitting parameters which determine range and slope of the linear region of the T(P) curve. As for the stroke volume–pressure curve, the following expression is used.

V Pð Þ = Vmax

1 +_β P Pv−1

 −K ð6Þ

whereβ and K are fitting parameters, Vmaxis the maximum stroke

volume, and Pvis the threshold for cardiac output.Table 1illustrates

the detailed parameter values used for numerical simulations. It shall be noted that the value of these parameters was estimated by best-fitting data drawn from physiological literature (Milnor, 1989; Korner et al., 1974).

Accompanying with the steady-state T_{–P as well as V–P curves,} Eqs. (3) and (4) can be reduced to a delay-differential equation with one dynamical variable P(t). By putting the obtained P(t) into Eq. (5), then the period of the cardiac cycle is derived.

3. Results 3.1. Normal subjects

Fig. 1 illustrates time evolution of R–R intervals with/without nonlinear noise reduction. The original R_{–R intervals shown in the} panel I ofFig. 1are obtained from a healthy male voluntarily recorded under supine position. The recording process approximately lasts 2 h. It is obvious to see that this data exhibits low-frequency trend embedded in high-frequency fluctuations. Following the previous assumption of linear combination of a clean signal and a noise source, the deterministic R–R intervals can be derived after nonlinear noise reduction, which is shown in the panel III ofFig. 1. Evidently, high-frequency fluctuations are significantly reduced and the low-frequency variability is left. In particular, the histogram of high-frequency noisyfluctuations exhibits a symmetrical Gaussian-type distribution (Fig. 2). The underlying physiological process for noisy fluctuations can be revealed in terms of the well-known power spectrum analysis, from which it is quite straightforward to under-stand the contribution of noisyfluctuations being main from the respiratory band (see Section 4). Therefore, the variability of R–R intervals coming from RSA is treated as a stochastic source under our reduction scheme.

The associated angle maps corresponding to the R–R intervals shown in panels I and III ofFig. 1are illustrated in panels II and IV ofFig. 1, respectively. In the recurrence plot of rotation angles for the original R–R intervals, it exhibits quite scattered behavior and does not show the deterministic one-dimensional structure

Table 1

List of parameters for model simulations R 1.2 × 103 dyn s/cm5 r 52 dyn s/cm5 C 1.0 × 10− 3 cm5 /dyn Pn 89 mm Hg Pv 25 mm Hg Ts 0.66 s Tm 1.2 s Vmax 86 cm3 α 31 γ 6.7 × 1013 _{β 72 K 7}

(panel II). Nevertheless, the angle map for the deterministic R–R intervals displays a clear one-dimensional and well-harmonized nonlinear motion (panel IV). Of particular note is that the intersection points of the diagonal and the one-dimensional deterministic curve are located at (π/4,π/ 4) and (−3π/ 4,−3π/4). These points are corresponding to thefixed heart rate. The point (π/4,π/ 4) is composed of three successive fixed R–R intervals which are larger than RRa. In other words, three successivefixed

R–R intervals, smaller than RRa, are mapped to the point (−3π/ 4,

−3π /4). From the dynamical point of view, these two points can be realized as saddle points, from which _{“deterministic” HRV is} regarded as intrinsically unstable.

3.2. ICU patients

Extracting the deterministic structure from noisy R–R intervals is of importance in fundamental research. Nevertheless, it is an interesting issue to test the nonlinear angle plot applied to risk stratification. So far, there is no consensus about the best available index of HRV for clinical use. Before we demonstrate the diagnostic results for ICU patients, we shall admit that using the dynamic characteristic of healthy individuals to predict the outcomes of ICU patients seems to violate the basic logic. However, this contradictory thinking contrarily inspires the different point of view compared to intuition. There is no doubt that there are two possible outcomes for ICU patients, in which one is transferred to

Fig. 1. Panels I and II denote the original R–R intervals and its corresponding angle map. Panels II and IV denote the deterministic R–R intervals and its corresponding angle map. Of particular note is that the recording process approximately lasts 2 h. In order to clearly visualize the variation of R–R intervals, only truncated heart beats are shown in panels I and III. But the angle maps shown in panels II and IV are derived from the complete R–R intervals.

Fig. 2. The histogram of high-frequency noisyfluctuations is obtained from taking the difference of the original R–R intervals (panel I ofFig. 1) and the deterministic R–R intervals (panel III ofFig. 1).

ordinary care units (OCUs) and the other is dying. The concept about OCU patients with stable signs of life seems reasonable, which definitely is the common point among healthy individuals. Therefore, using the dynamic characteristic of healthy individuals is to diagnose ICU patients with/without stable signs of life.

Fig. 3provides two typical Poincaré plots for ICU patients. Panel I displays normal sinus rhythm with a compact elliptical structure exception of some ectopic heartbeats. In panel III it is evident to see that ectopic heartbeats are predominant and, further, the compact elliptical structure disappears and is replaced by the fanned out pattern. To state more clearly, homoclinic orbits can be directly observed in the panel III, which is known to be resulted from various types of pathological states, e.g., symptomatic sinus node dysfunction (Bergfeldt and Haga, 2003). After the process of noise reduction, the associated angle maps corresponding to panels I and III ofFig. 3are illustrated in panels II and IV ofFig. 3, respectively. Panel II shows the deterministic one-dimensional structure analogous to that of healthy individuals, but in panel IV scattered points have a tendency tofill in the plane. Therefore, we may say that these two ICU patients have a distinct difference in mortality. Fig. 4 illustrates our dynamic risk stratification for ICU patients based upon the nonlinear map analysis. The last three days' ECG recordings of patient A are analyzed and shown inFig. 4(a), where stable signs of life are clearly demonstrated during these days. Therefore, we diagnose patient A was transferred to OCU. Our diagnosis is right according to the hospital records.Fig. 4(b) is the analyzed outcomes of patient B according to the last seven days' ECG recordings. The three angle plots from 21 Aug. to 23 Aug. exhibit deterministic one-dimensional curves. Nevertheless, the last four days' angle plots are quite scattered, from which we judge patient B was dying. Our prediction is proved again by the hospital records. It is deserved to explain why patient B had a critical transition from 23 August to 24 August. According to historical records, medical doctor suspected that the unexpected leakage after surgery could be the main reason. Therefore, using angle maps to exactly predict the unexpected leakage seems nontrivial, moreover, intensive studies in other populations performed by independent investigators are definitely necessary in that topic. In addition to patients A and B, the outcomes of other 23 patients are also well described by the present method. In order to quantify the relation between Eq. (3) and angle plots obtained from ICU patients, we also calculate the cross correlation as a measure which is 0.73 ± 0.13 for patients with stable signs of life (n = 10) and, however, the other is 0.25 ± 0.08 for critically ill patients (n = 15),

Fig. 3. Panels I and III denote two typical Poincaré plots for ICU patients. And the corresponding angle maps are shown in panels II and IV, respectively.

Fig. 4. The angle plot applied for the clinical diagnosis of two ICU patients, where different color symbols correspond to different measures. (a) It is evident to see the appearance of stable signs of life from patient A. Therefore, we judge this patient was transferred to OCU. (b) On the contrary, patient B lost stable signs of life, therefore, we judge this patient was dying. It shall be noted that the scattered points obtained from 24 August to 27 August have a tendency tofill in the plane. In order to clear visualization, we just show the data obtained from 27 August.

where the threshold separating the patient groups we used is 0.5. Therefore, there is a signi_{ficant difference between these two possible} outcomes of ICU patients. So far, it is known that a diagnosis of sepsis for ICU patients was associated with decreased total HRV, which strongly correlates with severity of illness (Garrard et al., 1993). However, in our study the degree of change in total HRV is not a good diagnostic index to describe illness severity of patient B. Thus, we might suggest that nonlinear map analysis on HRV applied to clinical utility could be promising. In order to give clear comparisons between different analyzed groups, the cross-correlation measure is listed in

Table 2. Of particular note is that analyzing healthy subjects from the physionet public website are also included inTable 2(www.physionet. org/physiobank/database). It is obvious to see that well-harmonized homeostasis still can be observed in this normal sinus rhythm R–R interval database.

3.3. Surrogate analysis

Surrogate analysis for the extracted deterministic data is per-formed and shown inTable 2, from which the extracted time series is statistically different from the surrogate. This rejects the null hypothesis and hence the fluctuations in the extracted time series have a deterministic/nonlinear structure. Moreover, surrogate analysis for the extracted noise is also performed in this study. The surrogate noise dataset was added to the extracted deterministic signal. Repeating the nonlinear noise-reduction procedure, wefind that the new-generated deterministic signal behaves the same angle plot as the old one. According to this finding, the hypothesis on noise embedded in R_{–R intervals is further verified.}

3.4. Model results

In these model simulations we only stress on the evolution of angle plots under the influence of time delay τ. When τ is smaller than 0.5 s, heart rate is not time dependent and the corresponding Poincaré plot shall exhibit afixed point. The spontaneous heart rate fluctuations appear whenτ is larger than 0.5 s, where τ=0.5 s is a threshold for occurring a supercritical Hopf bifurcation. In the range 0.5 s≤τ≤3 s the angle map displays the well-harmonized trajectory [Fig. 5(a)]. However, the well-harmonized profile will be gradually destructed due to the increase ofτ values [Fig. 5(b)], andfinally scattered points have a tendency tofill in the whole plane [Fig. 5(c)]. The basic reason for the angle plot shown inFig. 5(c) is the appearance of spiky-like heart ratefluctuations when τ is large enough, which is a well-known phenomenon for delay-differential equations. Spiky fluctuations indicate _{fixed heart rate cannot be temporally assembled, thus} scattered points in the angle plot cannot be accumulated at saddle points, i.e., (π/4,π/4) and (−3π/4,−3π/4). Instead, scattered points distribute in the whole plane. Compared to the map evolution of

patient B shown inFig. 4(b), we may say that the classic Windkessel model not only well describes why homeostasis can exist in heart rate fluctuations, but also provides a potential link between the clinical data and the circulatory system.

4. Discussions

Janson et al. (2001)introduced a model derived for the dynamics of angles of return times map of a periodic self-oscillatory system weakly forced by an arbitrary harmonic signal. The explicit map describing

Table 2

List of cross-correlation measures for different analyzed groups Data from Taipei Veterans General Hospital

Healthy subjects Male (n = 15) Female (n = 10) 0.75 ± 0.12 0.74 ± 0.15 0.01 ± 0.02a _{−0.01±0.03}a

ICU patients Recovery (n = 10) Deterioration (n = 15) 0.73 ± 0.13 0.25 ± 0.08 −0.03±0.04a

0.02 ± 0.03a

Data from the physionet public website

Healthy subjects Male (n = 15) Female (n = 13) 0.71 ± 0.11 0.72 ± 0.13 0.02 ± 0.03a

0.01 ± 0.03a a_{Means results from surrogate data.}

Fig. 5. Illustrations of the map evolution for differentτ values. τ=(a) 2.5 s, (b) 3.2 s, and (c) 4.3 s. In order to discriminate angle maps derived from the clinical data or model simulations,θnis used to represent model results.

the dynamics of angles for weak harmonic forcing is found to take the simple analytic form

n= arctan 2 cos 2ð πn− cot n−1Þ ð7Þ

whereξ is the rotation number, which is defined as the ratio of these two frequencies. Ifξ tends to a very small values, this map equation can well describe the angle plot shown in the panel IV of Fig. 1. Therefore, based upon this theoretical model, we can further realize that the deterministic one-dimensional structure we obtained might be mainly resulted from the interaction of two independent frequencies, which could be explained via the function of autonomic nervous system on HRV, i.e., a kind of well-harmonized homeostasis. In addition, Eq. (7) can be extended as the limiting case resulting from a larger number of oscillators, where the oscillating frequencies are located in two different narrow bands. Owing to this concept, the appearance of a little bit noisy structure in the panel IV ofFig. 1is expected.

The exhibition of autonomic nervous system can be traditionally realized via spectral analysis on HRV. In fact, interpretation of the spectrum itself is an active area of research. Usually the spectrum is broken into three regions for analysis. (a) The very low frequency (VLF) region covers from 0.000 to about 0.040 Hz. This region cannot usually be resolved but would be related with long-term factors such as thermoregulation of heart rate. (b) The low frequency range (LF, 0.040–0.150 Hz) often shows a peak at about 0.100 Hz, the origin of which is still unclear. Increased LF power may indicate sympathetic activation. (c) The high frequency region (HF, 0.150–0.400 Hz) covers rapid variations in heart rate due to vagal activity. In particular, human respiratory sinus arrhythmia is often seen between 0.180 and 0.400 Hz. And the central frequency for VLF, LF, and HF is corresponding to 0.003 Hz, 0.100 Hz, and 0.250 Hz, respectively.

Table 3illustrates linear measures, including SDNN as well as spectral analysis, for ICU patients. Obviously, most of these linear indices, exception of low-frequency power, are not good to identify the situations of ICU patients.

It is known that the spectral ratio LF/HF, i.e., the balance between sympathetic and vagal activities, was a classical homeostasis index for clinical diagnosis (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996). Nevertheless, Yien et al. (1997) clearly demonstrated that progressive increases in the power density values of both VLF and LF components appeared to be related to recovery for ICU patients. Conversely, progressive decreases in the power density values of these spectral components were indicative of deterioration and fatality. In the present study, the spectral power contributed by vagal tone is signi_{ficantly reduced by using nonlinear noise-reduction method. The} obtained deterministic one-dimensional angle map, in the point of view of autonomic tone, should be relevant to both VLF and LF. It is also deserved to note that the central-frequency ratio of VLF and LF is 0.03, which is close to the description by Eq. (7). To our best knowledge, it is little known that the correlation between VLF and LF

could be a prognostic homeostasis index for ICU patients. In other words, our _{findings are quite different from the classical risk} stratification (i.e., LF/HF) due to autonomic imbalance (Curtis and O'Keefe, 2002), and might be complementary to results reported by Yien et al.

Concerning the classic Windkessel model, it is known that the spectral power mainly falls into both VLF and LF bands when a longer time delay is used to simulate (greater than 2 s). However, HF component will be dominant in the power spectrum for a shorter time delay (smaller than 1 s). Berger et al. (1989) estimated distinctly different delays in response to vagal or sympathetic stimulations, where vagal mediated changes begin almost immediately (≈0.6 s) and sympathetic mediated changes may begin after 1.7–2 s. Therefore, comparisons between Figs. 4(b) and 5 we may suggest that the correlation between VLF and LF plays a critical role to determine whether well-harmonized homeostasis exists in heart rate fluctuations.

Till to now, several authors have studied nonlinear measures in order to test their feasibility to identify changes in autonomic nervous system. In particular, Hagerman et al. disrupted the autonomic nervous activity to the heart with propranolol and atropine and found a reduction in the largest Lyapunov exponent (Hagerman et al., 1996). This confirms the potential value of using measures of nonlinear dynamics as a tool for evaluating autonomic nervous system to the heart. Interestingly, Hagerman et al. could not totally eliminate the nonlinear structures in HRV by using combined blockade. Thus they concluded that other mechanisms like circulating hormones, preload, or afterload contribute to the nonlinearity in HRV. Here, we give an example to address the underlying mechanisms of HRV could be beyond the only consideration of autonomic nervous activities. It is no doubt that the function of brainstem is strongly related to autonomic nervous activities. However, we analyzed R–R intervals of a little boy (4 years old) with brainstem dysfunction, and found the cross-correlation measure is around 0.81, which has no significant difference with those of normal subjects (seeTable 1). This finding suggests that HRV could be from many factors rather than from autonomic nervous system, e.g., the classic Windkessel model. The detailed results of brainstem dysfunction will be published in elsewhere.

We would like to give more detailed remarks to address the possible difference between linear and nonlinear diagnoses. It is well known that linear SDNN method has a low accuracy for predicting the occurrence of life-threatening arrhythmias. What is the shortcoming of the SDNN method? We think it is a naturally born with problem. It is well accepted that cardiac dynamics is highly nonlinear. Bifurca-tions, chaos, and dynamical heterogenesis are all explored via experimental as well as theoretical studies in mammalian (Focus issue, 2002). A pronounced example of life-threatening arrhythmias is to consider ventricular tachycardia (VT), which is regarded as initiating abnormal/complicated spiral activities in ventricular tissue (Focus issue, 2002; Shiau et al., 2004). Intuitively, detecting spiral characteristics shall be analogous to HRV during VT period. These spirals are under the control of nonlinear evolution, therefore, the discovery of complexities buried in the nonlinear spiral could be beyond the capability of the linear method.

Finally, we would like to make our analyzed data available for analysis by other interested researchers. We wish that our data could be useful to the community of researchers.

5. Conclusions

Based upon the hypothesis of RSAfluctuations as a stochastic source, using the methods of noise reduction and angle map to extract the well-harmonized deterministic structure from R–R intervals is of great interest in both cardiac dynamics and nonlinear theory, where the concept of the interaction between different oscillators is suitable

Table 3

Results of SDNN and spectral analysis on ICU patients

ICU patients Recovery (n = 10) Deterioration (n = 15) P value

SDNN (ms) 20.9 ± 12.7 12.8 ± 8.9 ns LF (ms2 ) 107.5 ± 199.5 11.3 ± 16.3 b0.05 HF (ms2 ) 43.4 ± 72.7 14.8 ± 20.7 ns LF/HF 2.9 ± 1.9 2.3 ± 3.2 ns LFn(nu) 68.5 ± 13.2 50.7 ± 23.7 ns HFn(nu) 31.2 ± 14.2 48.3 ± 24.7 ns

SDNN, standard deviation of the R–R intervals; LF, Low-frequency power; HF, high-frequency power; LF/HF, the ratio of Low-high-frequency to high-high-frequency power; LFn, LF in

normalized units; HFn, HF in normalized units; and ns, no significance. All variables

were obtained from 15-minute recordings and the level of statistical significance was set at Pb0.05.

to interpret ourfindings. Significant reduction of the contribution of the HF component,_{fixed heart rate is verified as intrinsically unstable,} which could be influenced by both VLF and LF components. Particularly, the cross-correlation measure between the theoretical angle map and the numerically derived angle trajectory provides an additional non-invasive index in clinical research. Moreover, a simple physiologic model under the consideration of the influence of time delay in the human baroreceptor-mediated reflex is used to explain the clinical data. We wish that the discovery of well-harmonized homeostasis in heart rate fluctuations can raise more attractive studies in the future.

Acknowledgements

We really appreciate Dr. Sellier (Catania Univ., Italy) and Dr. Gong (Cornell Univ., USA) critically read the present manuscript and deliver encouraging comments. And we also thank for the research group of Dr. Yien (MD, Taipei Veterans General Hospital) providing HRV data This research was supported in part by the Institute of Physics, Academia Sinica, Taipei, Taiwan, the National Science Council, ROC, Grant No. NSC 93-2112-M-259-010, and Taipei Veterans General Hospital under Contract Nos. VGH 91-365-5.

References

Andrés, D.S., Irurzun, I.M., Mitelman, J., Mola, E.E., 2006. Increase in the embedding dimension in the heart rate variability associated with left ventricular abnormal-ities. Appl. Phys. Lett. 89, 144111.

Bergfeldt, L., Haga, Y., 2003. Power spectral and Poincaré plot characteristics in sinus node dysfunction. J. Appl. Physiol. 94, 2217–2224.

Balocchi, R., Menicucci, D., Santarcangelo, E., Sebastiani, L., Gemignani, A., Ghelarducci, B., Varanini, M., 2004. Deriving the respiratory sinus arrhythmia from the heartbeat time series using empirical mode decomposition. Chaos, Solitons Fractals 20, 171–177.

Berger, R.D., Akselrod, S., Gordon, D., Cohen, R.J., 1986. An efficient algorithm for spectral analysis of heart rate variability. IEEE Trans. Biomed. Eng. BME 33, 900–904. Berger, R.D., Saul, J.P., Cohen, R.J., 1989. Transfer function analysis of autonomic

regulation in canine atrial rate response. Amer. J. Physiol. 256, H142–H152. Curtis, B.M., O'Keefe, J.H., 2002. Autonomic tone as a cardiovascular risk factor: the

dangers of chronicfight or flight. Mayo. Clin. Proc. 77, 45–54 and references cited therein.

FOCUS ISSUE, 2002. Mapping and control of complex cardiac arrhythmias. Chaos 12, 732–981 and references cited therein.

Franz, G.N., 1969. Nonlinear rate sensitivity of the carotid sinus reflex as a consequence of static and dynamic nonlinearities in baroreceptor behavior. Ann. NY. Acad. Sci. 811–824.

Goldberger, A.L., 1996. Nonlinear dynamics for clinicians: chaos theory, fractals, and complexity at the bedside. Lancet 347, 1312–1314.

Goldberger, A.L., Rigney, D.R., West, B.J., 1990. Chaos, fractals and physiology. Sci. Am. 262, 42–49.

Goldberger, A.L., Mie, Tus Je., Rigney, D.R., Wood, M.I., Fortney, S.M., 1994. Effects of head-down bed rest on complex heart rate variability: response to LBNP testing. J. Appl. Physiol. 77, 2863–2869.

Garrard, C.S., Kontoyannis, D.A., Piepoli, M., 1993. Spectral analysis of heart rate variability in the sepsis syndrome. Clin. Auton. Res. 3, 5–13.

Huikuri, H.V., Seppanen, T., Koistinen, M.J., Airaksinen, K.E.J., Ikaheimo, M.J., Castellanos, A., Myerbutg, R.J., 1996. Abnormalities in beat-to-beat dynamics of heart rate before the spontaneous onset of life-threatening ventricular tachyarrhythmias in patients with prior myocardial infarction. Circulation 93, 1836–1844.

Hagerman, I., Berglund, M., Lorin, M., Nowak, J., Sylvén, C., 1996. Chaos-related deterministic regulation of heart rate variability in time and frequency domains: effects of autonomic blockade and exercise. Cardiovasc. Res. 31, 410–418. Hegger, R., Kantz, H., 1999. Improved false nearest neighbor method to detect

determinism in time series data. Phys. Rev. E 60, 4970–4973.

Ivanov, P.C., Amaral, L.A.N., Goldberger, A.L., Stanley, H.E., 1998. Stochastic feedback and the regulation of biological rhythms. Europhys. Lett. 43, 363–368.

Janson, N.B., Balanov, A.G., Anishchenko, V.S., McClintock, P.V.E., 2001. Modelling the dynamics of angles of human R–R intervals. Physiol. Meas. 22, 565–579. Kantz, H., Schreiber, T., 1998. Human EGG: nonlinear deterministic versus stochastic

aspects. IEE Proc.-Sci. Meas. Technol. 145, 279–284.

Kennel, M.B., Brown, R., Abarbanel, H.D.I., 1992. Determining embedding dimension for phase-space reconstruction using a geometrical construction. Phys. Rev. A 45, 3403–3411.

Korner, P.I., West, M.J., Shaw, J., Uther, J.B., 1974. Steady-state properties of baroreceptor heart rate reflex in essential hypertension in man. Clin. Exp. Pharmacol. Physiol. 1, 65–76.

Milnor, W.R., 1989. Hemodynamics. Williams & Wilkins, Baltimore.

Rosenblum, M., Kurths, J., 1995. A model of neural control of the heart rate. Physica A 215, 439–450.

Seidel, H., Herzel, H., 1998. Bifurcations in a nonlinear model of the baroreceptor– cardiac reflex. Physica D 115, 145–162.

Suder, K., Drepper, F.R., Schiek, M., Abel, H.H., 1998. One-dimensional, nonlinear determinism characterizes heart rate pattern during paced respiration. Am. J. Physiol. 275, H1092–H1102.

Schreiber, T., 1993. Extremely simple nonlinear noise-reduction method. Phys. Rev. E 47, 2401–2404.

Shiau, Y.H., Hsueh, M.P., Hseu, S.S., Yien, H.W., 2004. Complicated electrical activities in cardiac tissue. Int. J. Mod. Phys. B 18, 2645–2650.

Theiler, J., Eubank, S., Longtin, S., Galdrikian, B., Farmer, J., 1992. Testing for nonlinearity in time series: the method of surrogate data. Physca D 58, 77–94.

Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996. Heart rate variability: standards of measure-ment, physiological interpretation, and clinical use. Circulation 93, 1043–1065 and references cited therein.

See Fig. 11(c) inWu, M.C., Huang, M.C., Yu, H.C., Chiang, T.C., 2006. Phase distribution and phase correlation offinancial time series. Phys. Rev. E 73, 016118.

Westerhof, N., Elzinga, G., Sipkema, P., 1971. An artificial arterial system for pumping hearts. J. Appl. Physiol. 31, 776–781.

Yien, H.W., Hseu, S.S., Lee, L.C., Kuo, T.B.J., Lee, T.Y., Chan, S.H.H., 1997. Spectral analysis of systemic arterial pressure and heart rate signals as a prognostic tool for the prediction of patient outcome in the intensive care unit. Crit. Care Med. 25, 258–266.

Other uses, including reproduction and distribution, or selling or

licensing copies, or posting to personal, institutional or third party

websites are prohibited.

In most cases authors are permitted to post their version of the

article (e.g. in Word or Tex form) to their personal website or

institutional repository. Authors requiring further information

regarding Elsevier’s archiving and manuscript policies are

encouraged to visit:

Nonlinear measures on heart rate variability: A clinical tool or not?

To the Editor:

Spectral analysis on heart rate variability (HRV) has been a widely accepted linear method in the assessment of autonomic nervous system (ANS) (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996), in which HRV reveals a delicate balance between the two antagonistic parts of ANS: sympathetic and parasympathetic activities. Recently,Yeh et al. (2009)

used a nonlinear method, i.e., detrended fluctuation analysis (DFA) developed by Peng et al. (1995), to analyze HRV in late pregnant women. And their main finding, published in this journal, is that detrendedfluctuation scaling exponents might be new and indepen-dent measures of HRV in late pregnancy, in addition to those conventional linear measures. Hence, the underlying physiological meanings of these scaling exponents in pregnant women are still unclear. Owing to their results, we are motivated to reconsider this topic “Nonlinear measures on heart rate variability: A clinical tool or not?”.

It is well known that DFA is one of the powerful methods to detect self-similar scaling exponents (i.e., nonlinear indices) embedded in nonstationary biomedical signals. The use of nonlinear indices is usually justified by considering them as numerical indices, which are supported by clinical studies in large-scale database. However, in fact, the clinical interpretation of changes in these indices is quite difficult. One of the possible reasons is that they compress cardio-respiratory dynamics including heart rate, blood pressure, and respiration rate in a single numerical index, thus it is not easy to determine physiological relationships between the change in the index and the physiopath-ological causes. Nevertheless, to discover the clear physiphysiopath-ological meaning of these nonlinear indices is of fundamental importance.

Ursino and Magosso (2003)developed a detailed lumped parameter model (LPM) which contains the interconnection of ANS (sympatho-vagal activity) and non-autonomic components (peripheral resis-tance). According to numerical results of LPM (Rojo-Alvarez et al., 2007), it has been known that an increase of sympathetic tone or a decrease of vagal tone can significantly enhance the short-term α1 index (determined from 4 to 11 beats). In addition, the long-termα2 index (> 11 beats) can be significantly affected by the decrease of peripheral resistance.

There is a more natural way to realize the physiological meaning of nonlinear indices based upon traditional spectral analysis on HRV. Usually the spectrum is broken into three regions for analysis. (a) The very low-frequency (VLF) region covers from 0.000 to about 0.040 Hz. This region cannot usually be resolved but would be related with long-term factors such as thermoregulation of heart rate. (b) The low-frequency range (LF, 0.040–0.150 Hz) often shows a peak at about 0.100 Hz, the origin of which is still unclear. Increased LF power may indicate sympathetic activation. (c) The high-frequency region (HF, 0.150–0.400 Hz) covers rapid variations in heart rate due to vagal activity. It is known that in the rest status heart rate for normal subjects

is 60–100 bpm, then one heartbeat needs 0.6–1 s. Therefore, the time characteristic forα1 index covers from 2.4 s to 11 s, and the associated frequency range would be 0.417–0.090 Hz which is located in both of HF and LF bands. This is why the sympathovagal activity is strongly related to_{α1 index. Accordingly, α2 index covers LF and VLF bands.} This simple statement is consistent with the result byWillson et al. (2002). Compared to the previous physiological/theoreticalfindings of

Ursino and Magosso (2003), Rojo-Alvarez et al. (2007), and Willson et al. (2002), the conservative working hypothesis to reveal the physiological meaning of nonlinear indices would be thatα1 index should reflect the result of the sympathovagal balance (or imbalance), and_{α2 index is still unclear at present.}

It is well known that low-/high-frequency power ratio (LHR) is a well-accepted index to describe the sympathovagal activity (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996). Intuitively,α1 index should be strongly correlated with LHR. This concept we think is correct in principle. In postural-change experiments (Montano et al., 1994; Malliani et al., 1997) people have realized that sympathetic activity will be enhanced in the standing position because of gravitation. Hence, decreasing sympathetic activity would be expected in the supine position. Increasing (decreasing) sympathetic activity indicates low-frequency oscillations (LFOs) of HRV will be enhanced (reduced), thus LHR should be correspondingly increased (decreased). In the viewpoint of time series, LFOs describe the characteristic of positive correlations embedded in HRV data. In fact, DFA is an outstanding method to detect positive (or negative) correlations in biomedical signals, which can be expressed through the monofractal scaling exponent,α index. Therefore, it is not surprising that there is a mathematical relation between DFA and spectral analysis (Peng et al., 1995). It is known that the spectral density can be approximately expressed as a power-law distribution f− β, whereβ is an exponent andα and β indices obey the following relation:

β = 2α−1 ð1Þ

For example, when Gaussian white noise of the flat band is considered,β=0 and, therefore, α=0.5. LHR in mathematics should be equal to 0.44 (LF, 0.040–0.150 Hz; HF, 0.150–0.400 Hz). If β (α) is larger than zero (0.5), it indicates positive correlations are predominant in the original data, thus LHR will be increased over 0.44. However, it needs more careful treatments in clinical applications. As mentioned above,α1 index would be responsible for the sympathovagal balance (or imbalance) and the frequency characteristic for_{α1 index covers} from 0.090 Hz to 0.417 Hz, which does not occupy the entire LF band andfinally leads to the uncertainty for spectral-power calculations. Owing to that, LHR should be defined in a new frequency-domain rather than a traditional one. This procedure might give the correct relation betweenα1 index and LHR. Our interpretation is valid to address why the study of the correlation coefficient between α1 index and LHR in previous articles, see the text in the paper ofYeh et al. (2009), has not obtained consistent results. Of course, in practical aspects it still has a space to discuss the uncertainty embedded in LHR (see below).

1566-0702/$– see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.autneu.2009.09.016

Contents lists available atScienceDirect

Autonomic Neuroscience: Basic and Clinical

j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / a u t n e u

In mathematics spectral analysis is valid when time series exhibits the stationary characteristic. HRV, in fact, is a result from nonstation-ary cardio-respiratory processes, therefore, a crucial procedure to be simultaneously performed is to obtain respiratory rate in order to assess its synchronization with HF component. However, when the frequency of respiration is decreased and close to the LF rhythm in such a way that HF and LF merge into one single dominant oscillation. Thus, it is difficult to define LHR in this case. Although using controlled breathing is a way to maintain the frequency of respiration above the LF range, nevertheless, this is not physiological breathing and it can artificially shift the sympathovagal balance towards vagal predomi-nance. In short, it would not be surprising that previous studies, cited in the paper ofYeh et al. (2009), on correlation results between LHR and α1 index with/without controlled breathing have not given consistent conclusions. The possible way to overcome the present situation we suggest should measure respiration as well as HRV data simultaneously, which would be helpful to accurately define LHR rather than traditional usage in spectral analysis.

Although nonstationary cardio-respiratory processes will increase difficulties to realize the physiological basis of DFA and its relation to linear HRV measures, there is a mathematical relation between DFA and spectral analysis as simply described above. Thus, it is incorrect/ inappropriate to describe that DFA scaling exponents will not correlate with most HRV measures in the theoretical point of view, by which as well as previous inconsistent resultsYeh et al. (2009)claimed thatα1 index might be a new and independent measure of HRV, in addition to the conventional time and frequency-domain HRV measures.

In this letter, we try to give a better understanding about the relation between _{α1 index and LHR through various viewpoints,} including physiological/theoretical models and clinical symptoms during the postural change. We think nonlinear measures that resulted from sophisticated maths would become powerful tools for clinical applications. However, the understanding of the physiological basis of nonlinear measures would play a crucial role for future applications (Shiau, 2009), rather than considering them as numerical indices supported by clinical studies in large-scale database.

References

Malliani, A., Pagani, M., Furlan, R., et al., 1997. Individual recognition by heart rate variability of two different autonomic profiles related to posture. Circulation 96, 4143–4145.

Montano, N., Gnecchi Ruscone, T., Porta, A., Lombardi, F., Pagani, M., Malliani, A., 1994. Power spectrum analysis of heart rate variability to assess the changes in sympathovagal balance during graded orthostatic tilt. Circulation 90, 1826–1831. Peng, C.K., Havlin, S., Stanley, H.E., Goldberger, A.L., 1995. Quantification of scaling

exponents and crossover phenomena in nonstationary heartbeat time series. Chaos 5 (1), 82–87.

Rojo-Alvarez, J.L., Sanchez-Sanchez, A., Barquero-Perez, O., Goya-Esteban, R., Everss, E., Mora-Jimenez, I., Garcia-Alberola, A., 2007. Analysis of physiological meaning of detrendedfluctuation analysis in heart rate variability using a lumped parameter model. Comput. Cardiol. 34, 25–28.

Shiau, Y.H., 2009. Does well-harmonized homeostasis exist in heart ratefluctuations? Time series analysis and model simulations. Auton. Neurosci. Basic Clin. 146, 62–69. Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996. Heart rate variability: Standards of measure-ment, physiological interpretation, and clinical use. Circulation 93, 1043–1065, and references cited therein.

Ursino, M., Magosso, E., 2003. Role of short-term cardiovascular regulation in heart period variability: a modeling study. Am. J. Physiol. Heart Circ. Physiol. 284, H1479–H1493. Willson, K., Francis, D.P., Wensel, R., Coats, A.J., Parker, K.H., 2002. Relationship between

detrended fluctuation analysis and spectral analysis of heart-rate variability. Physiol. Meas. 23, 385–401.

Yeh, R.G., Shieh, J.S., Chen, G.Y., Kuo, C.D., 2009. Detrendedfluctuation analysis of short-term heart rate variability in late pregnant women. Auton. Neurosci. Basic Clin.

doi:10.1016/j.autneu.2009.05.241.

Yuo-Hsien Shiau Graduate Institute of Applied Physics, National Chengchi University, Taipei 11605, Taiwan, ROC E-mail address:[email protected].

journal homepage:www.elsevier.com/locate/physa

Concurrent sympathetic activation and vagal withdrawal in

hyperthyroidism: Evidence from detrended fluctuation analysis of heart

rate variability

Jin-Long Chen

a,b

, Yuo-Hsien Shiau

c

, Yin-Jiun Tseng

d

, Hung-Wen Chiu

e

, Tzu-Chien Hsiao

f

,

Niels Wessel

b,g

, Jürgen Kurths

g,h

, Woei-Chyn Chu

d,

∗

a_{Department of Medical Informatics, Tzu Chi University, Hualien, Taiwan} b_{Institute of Physics and Astronomy, University of Potsdam, Potsdam, Germany} c_{Graduate Institute of Applied Physics, National Chengchi University, Taipei, Taiwan} d_{Institute of Biomedical Engineering, National Yang-Ming University, Taipei, Taiwan} e_{Graduate Institute of Medical Informatics, Taipei Medical University, Taipei, Taiwan} f_{Department of Computer Science, National Chiao Tung University, Hsinchu, Taiwan} g_{Institute of Physics, Humboldt University, Berlin, Germany}

h_{Potsdam Institute for Climate Impact Research, Potsdam, Germany}

a r t i c l e

i n f o

Article history:

Received 29 August 2009 Available online 13 January 2010 Keywords:

Hyperthyroidism Autonomic nervous system Heart rate variability Detrended fluctuation analysis Nonlinear dynamics

a b s t r a c t

Despite many previous studies on the association between hyperthyroidism and the hy-peradrenergic state, controversies still exist. Detrended fluctuation analysis (DFA) is a well recognized method in the nonlinear analysis of heart rate variability (HRV), and it has phys-iological significance related to the autonomic nervous system. In particular, an increased short-term scaling exponentα1 calculated from DFA is associated with both increased sym-pathetic activity and decreased vagal activity. No study has investigated the DFA of HRV in hyperthyroidism. This study was designed to assess the sympathovagal balance in hyper-thyroidism. We performed the DFA along with the linear analysis of HRV in 36 hyperthy-roid Graves’ disease patients (32 females and 4 males; age 30±1 years, means±SE) and 36 normal controls matched by sex, age and body mass index. Compared with the normal controls, the hyperthyroid patients revealed a significant increase (P<0.001) inα1 (hy-perthyroid 1.28±0.04 versus control 0.91±0.02), long-term scaling exponentα2(1.05± 0.02 versus 0.90±0.01), overall scaling exponentα (1.11±0.02 versus 0.89±0.01), low frequency power in normalized units (LF%) and the ratio of low frequency power to high frequency power (LF/HF); and a significant decrease (P < 0.001) in the standard devia-tion of the R–R intervals (SDNN) and high frequency power (HF). In conclusion, hyperthy-roidism is characterized by concurrent sympathetic activation and vagal withdrawal. This sympathovagal imbalance state in hyperthyroidism helps to explain the higher prevalence of atrial fibrillation and exercise intolerance among hyperthyroid patients.

© 2010 Elsevier B.V. All rights reserved.

1. Introduction

Hyperthyroidism is characterized by clinical manifestations that resemble those of a hyperadrenergic state. Moreover, the fact that beta-adrenergic receptor blockers ameliorate these symptoms and signs has further suggested enhanced

∗_{Corresponding address: Institute of Biomedical Engineering, National Yang-Ming University, 155 Li-Nong St. Sec. 2, Beitou, Taipei 112, Taiwan. Tel.:}

+886 2 28267025; fax: +886 2 28210847. E-mail address:[email protected](W.-C. Chu).

0378-4371/$ – see front matter©2010 Elsevier B.V. All rights reserved.

receptors [7,8]; alterations in the quantity of guanine nucleotide-binding proteins [9]; an increase in catecholamine turnover at neural synapses [5,6]; and structural similarities between thyroid hormones and catecholamines [10]. Notwithstanding these possibilities, the notion that hyperthyroidism is associated with an increased sympathetic tone could not be fully ver-ified from these studies. On the other hand, a smaller increase in heart rate induced by atropine during the hyperthyroid state compared with the euthyroid state implied reduced vagal inhibition of heart rate in hyperthyroidism [11,12].

Analysis of the heart rate variability (HRV) provides a non-invasive and sensitive tool for the evaluation of autonomic regulation of the heart [13,14]. In clinical applications, reduced HRV is associated with increased cardiac mortality after acute myocardial infarction [15] and is an early warning sign of diabetic neuropathy among diabetic patients [16,17]. HRV analysis can be categorized into linear and nonlinear methods [14]. Conventionally, the beat-to-beat variation exhibited by the sinoatrial node is analyzed using linear methods. Previous studies have applied the linear analysis of HRV to investigate the autonomic nervous system of hyperthyroid patients; some have disclosed reduced [18–20] or normal [21] vagal activity, whereas others have shown both increased sympathetic and decreased vagal modulation of the heart rate in patients with hyperthyroidism [2,22].

However, multiple nonlinear mechanisms such as sympathetic nerves, vagal nerves, hormones and hemodynamics are involved in the regulation of the heart rate and these affecting factors interact mutually. Consequently, the heart rate regulating system appears to be a possible example where chaos theory can be applied [23,24]. These nonlinear phenomena could affect the genesis of heart rate fluctuation [25] and therefore a nonlinear analysis of HRV would be a more appropriate approach to interpreting the complex phenomena of heart rate dynamics. A recent study using a nonlinear analysis of HRV with the correlation dimension for hyperthyroidism has shown reduced complexity and impaired tolerance to cardiovascular stresses in hyperthyroid patients [26].

Nonlinear analysis of HRV can be quantified using parameters derived from chaos and fractal theory [27]. Detrended fluctuation analysis (DFA) is a nonlinear analysis of HRV [28]. Recent studies have indicated that DFA could be used not only to differentiate various patient groups from normal controls [29,30] but also to stratify high risk patient groups among post-myocardial infarction patients [31,32]. In addition, the short-term scaling exponent

α

1, which is calculated from the DFA, could be related to the state of the autonomic nervous system. An increased

α

1 is associated with concurrent increased sympathetic activity and decreased vagal activity. Conversely, a decreased

α

1 is related to co-activation of both the sympathetic and the vagal components of the autonomic nervous system [33].

At present, no study has investigated the DFA of HRV in hyperthyroidism. This study was designed to assess the autonomic nervous system in hyperthyroidism by the nonlinear analysis of HRV with DFA. We hypothesized that the autonomic dys-function in hyperthyroid patients is caused by the joint effect of increased sympathetic activity and decreased vagal activity.

2. Subjects and methods

2.1. Subjects

A group of 36, newly diagnosed, untreated hyperthyroid Graves’ disease patients from the outpatient clinic of a university hospital and a group of 36 healthy normal control subjects were recruited for this study. The hyperthyroid and control groups were matched for sex (32 females and 4 males versus 32 females and 4 males, hyperthyroid versus control), age (30

±

1 versus 29

±

1 years, means

±

SE) and body mass index (20.7

±

0.4 versus 21

.

8

±

0

.

5 kg

/

m2_{). The diagnosis of Graves’}

disease was established on the basis of clinical, biochemical, immunological, thyroid scintigraphic scanning and uptake data. Individuals with diabetes, cardiac arrhythmia, cardiovascular disease, pregnancy or those using medication were excluded. The study protocol was approved by the local ethics committee and all participants gave their informed consent. The study was conducted according to the principles of the Helsinki declaration.

2.2. Study protocol

The hyperthyroid patients were studied at the time of diagnosis before any medication was administered. For all participants, no alcoholic or caffeine-containing drinks were taken for at least 24 h before the study. The examination was performed in a quiet room during the daytime. Subjects received one-channel electrocardiogram (ECG) measurement for 30 min in the supine position after five minutes rest. During the ECG measurement, the subjects were instructed to fully relax, stay awake, breathe regularly, and not to speak.

2.3. Measurement of the ECG

The acquired analog ECG signals were transformed into digital signals by a 16-bit analog-to-digital converter with a sampling rate of 500 Hz. The digitized ECG signals were processed off-line. First, the R waves were detected and then