Improving the accuracy of depth of anaesthesia using modified detrended fluctuation analysis method

Biomedical Signal Processing and Control 5 (2010) 59–65

Contents lists available at ScienceDirect

Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc

Improving the accuracy of depth of anaesthesia using modified detrended fluctuation analysis method T. Nguyen-Ky *, Peng Wen, Yan Li University of Southern Queensland, Toowoomba, QLD 4350, Australia

A R T I C L E I N F O

A B S T R A C T

Article history: Received 3 October 2008 Received in revised form 23 January 2009 Accepted 4 March 2009 Available online 8 April 2009

This paper presents a modified detrended fluctuation analysis (MDFA) to improve the monitoring accuracy of the depth of anaesthesia (DoA). We first use MDFA to classify anaesthesia state levels into awake, light, moderate, deep and very deep states. Then we build up five zones using linear regression method from very deep anaesthesia state to awake state, corresponding with different box sizes. Finally, the Lagrange method is applied to compute the DoA. Comparing with the most popular Bispectral Index (BIS) method, our modified DFA method extends the ranges of the moderate anaesthesia, deep anaesthesia and very deep anaesthesia to provide more information about the DoA. This extension is very significant in the clinical perspective as these states are within the ranges for operations and need more attention. Simulation results demonstrate that the new technique monitors the DoA in all anaesthesia states accurately. ß 2009 Elsevier Ltd. All rights reserved.

Keywords: Depth of anaesthesia Modified detrended fluctuation analysis Linear regression Lagrange method Electroencephalogram

1. Introduction Anaesthesia can be defined as a lack of response and recall to noxious stimuli. It includes the triad of paralysis, unconsciousness and analgesia. Depth of anaesthesia (DoA) depends on: equilibration of drug concentrations in plasma, the relationship between drug concentration and drug effects, and the influence of noxious stimuli [1]. Clinical assessment of depth of anaesthesia which uses autonomic signs, such as pulse, blood pressure, sweating and lacrimation is popular but is far from reliable. Several devices have been developed in recent years to monitor the DoA. Currently BIS based on electroencephalogram (EEG) is the most common one in hospitals. The BIS is based on EEG power spectrum and phase spectrum. It quantifies the coupling of phase angles of different frequencies. The BIS integrates several descriptors of the EEG into a single variable [2–4]. A BIS value of 0 indicates an isoelectric EEG signal. A BIS value nearing 100 indicates an ‘‘awake’’ clinical state. The anaesthesia state levels can be divided as the following: very deep anaesthesia state has BIS value from 0 to 20, deep anaesthesia state from 20 to 40, moderate anaesthesia state from 40 to 60, light anaesthesia state from 60 to 80 and awake state has BIS value from 80 to 100. Maintaining BIS values in the range of 40–60 provides

* Corresponding author. E-mail addresses: [email protected] (T. Nguyen-Ky), [email protected] (P. Wen), [email protected] (Y. Li). 1746-8094/$ – see front matter ß 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.bspc.2009.03.001

adequate hypnotic effect during balanced general anaesthesia and helps improve emergence and recovery. Therefore, the determination of anaesthesia state levels is very important. Clinicians can use this information to help guide anaesthetic dosing. There are also other methods developed recently. In [5], the authors presented a model that generalizes the autoregressive class of poly-spectral models by having a semi-parametric description of the residual probability density. Their results indicate that in two out of three anaesthetic agents, better classification can be achieved with higher order spectral. Their analysis also suggested that propofol had the strongest effect on the EEG signal recordings. Another new approach for quantifying the relationship between brain activity patterns and depth of anaesthesia is presented in [6] by analyzing the spatio-temporal patterns in the EEG using Lempel-Ziv complexity analysis. Compared with other indexes, such as approximate entropy, spectral entropy and median frequency, their results not only demonstrate better performance across all of the patients, but also an algorithm which is easy to implement for real-time use. In [7], a fully automated system is developed for the DoA estimation. The system determines the anaesthesia depth by assessing the characteristics of the mid-latency auditory evoked potentials (MLAEP). The discrete time wavelet transformation is used for compacting the MLAEP which localizes the time and the frequency of the waveform. The authors in [8] use an artificial neural network (ANN) to integrate different EEG variables in order to estimate the depth of anaesthesia. Comparison with correspondent BIS confirms

60

T. Nguyen-Ky et al. / Biomedical Signal Processing and Control 5 (2010) 59–65

that neural network could be trained to predict anaesthesia depth directly from the raw EEG. Another approach based on the analysis of a single-channel EEG signal using stationary wavelet transform is applied to study the cortical activity during general anaesthesia [9]. The wavelet coefficients calculated from the EEG are pooled into a statistical representation. The results show that the waveletbased anaesthetic value and BIS are well correlated during periods of steady-state. Recently, entropy and complexity of the EEG have been proposed as measures of depth of anaesthesia and sedation [10]. In [11], the detrended fluctuation analysis (DFA) is used to study the scaling behavior of the EEG as a measure of the level of consciousness. This method allows real-time implementation and enables its application in monitoring devices. The DFA algorithm finds a trend on the integrated EEG signals. The integrations on EEG signals are divided in segments of a length. The length of the segments is called box size. This trend is iteratively searched on different sizes of boxes over the signals. When the trend is found in each box, it is subtracted from the integrated signals. To investigate the effect of box sizes, [12] gives comprehensive reviews about the fractal and self-similarity properties of the EEG signals. It managed to find the optimum fractal-scaling exponent by selecting the best domain of box sizes, which have meaningful changes with different depth of anaesthesia. However, their results could not discriminate between light anaesthesia and moderate anaesthesia states. In this paper, we propose to monitor the DoA using a modified detrended fluctuation (MDFA), linear regression and Lagrange methods. We use the MDFA method to identify the EEG signals in different states of DoA: awake, light anaesthesia, moderate anaesthesia, deep anaesthesia and very deep anaesthesia. Using linear regression, we build five zones from very deep anaesthesia to awake state, corresponding with different box sizes. The Lagrange method is applied to compute the DoA. This paper is organized as follows. In Section 2, we introduce the DFA and MDFA algorithms to classify EEG signals. Section 3 presents the linear regression that is applied to choose box sizes. Implementation and simulation results are presented in Section 4. Finally, conclusions are drawn in Section 5. 2. Methodology 2.1. Detrended fluctuation analysis method Peng et al. [13] developed the DFA algorithm to estimate scaling exponents with local detrending to remove the nonstationary components. The time series is firstly integrated and divided in segments of length s. Each segment is then detrended by subtracting the best linear fit pn(i). Finally, the fluctuation function F(s) is calculated as the root mean square of the detrended time series as a function of the segment size s. If the time series is self-similar, the fluctuation function F(s) increases a by a power-law: F(s) a s . The scaling exponent a can be estimated by a linear fit on the log–log plot of F(s) versus s. The scaling exponent a represents the correlation properties of the signal: if a = 0.5, the signal is uncorrelated (white noise); if a < 0.5, the signal is anticorrelated (negative correlations); if a > 0.5, there are positive correlations in the signals [14–16]. For a given discrete-time sequence xn of length L, the DFA is calculated as the following:  Determine the profile YðiÞ ¼

i X ½xn  hxi; n¼1

i ¼ 1; . . . ; L;

(1)

where hxi is the average of the time series hxi ¼

L 1X xn : L n¼1

(2)

 Cut the profile Y(i) into Ls = [L/s], calculate the local trend for each segment n by a least-square fit pn(i) Y s ðiÞ ¼ YðiÞ  pn ðiÞ;

(3)

where pn(i) is the fitting polynomial in the nth segment, n = 1,. . ., 2Ls. Linear, cubic, or higher order polynomials can also be used in the fitting procedure (DFA1, DFA2, and higher order DFA). Since the detrending of the time series is done by the subtraction of the fits from the profile, these methods differ from their capabilities of eliminating trends in the data. In mth order DFA, trend’s order of m in the profile and order of m  1 in the original record are eliminated.  Calculate the variance for each of the 2Ls segments Fs2 ðnÞ ¼ hYs2 ðiÞi ¼

s 1X Y 2 ½ðn  1Þs þ i; s i¼1 s

 Average over all segments and take the square root " #1=2 2Ls 1 X Fs2 ðnÞ : FðsÞ ¼ 2Ls n¼1

(4)

(5)

For the EEG time series, the collected data have five states of anaesthetic depth. In Fig. 1, states A, B, C, D and E have BIS values of 15, 30, 50, 80 and 97, which are corresponding to very deep anaesthesia, deep anaesthesia, moderate anaesthesia, light anaesthesia and awake states. The source of the signals used in this paper is from and discussed in [12]. The raw EEG signals were recorded by BIS XP monitor (Aspect Medical System Inc.) through contact electrodes placed on a patient’s forehead. The length of the data is 10 s with 128 Hz sampling rate [12] (see [12] for more details). The relationship between F(s), the average fluctuation as a function of segment size (box size), and the segment size s is shown in Fig. 2. Typically, F(s) will increase with s. We define sA, sB, sC, sD and sE as box sizes for states A, B, C, D and E, respectively. In Fig. 2, we can distinguish the different states of anaesthesia. With the same value s (sA = sB = sC = sD = sE), when BIS values increase from state A (BIS = 15) to state E (BIS = 97), DFA values decrease from FA to FE. These results show the relation between BIS and DFA values and are consistent with the results obtained in [11,12]. However, we can see that DFA values of states C and D are very close in Fig. 2. It is impossible to clearly distinguish the two states C and D from the DFA values. In the next section, we use modified DFA to classify five states and separate states C, D. 2.2. Modified detrended fluctuation algorithm (MDFA) to classify anaesthesia states In DFA method, we do not plot FDFA in log–log representation and use the power-law scaling because they cannot indicate the DoA values. Moreover, many F(s) values in DFA results are overlapped in states C and D (see Fig. 2) and are also overlapped in log–log representation. We, therefore, apply modified DFA to detect and separate states C and D. 2.2.1. Classify five anaesthesia states We propose modified DFA to classify five states A, B, C, D and E. When calculating the variance for each of the 2Ls segments of the detrended time series Ys(i), we modify Eq. (4) as: Fs ðnÞ ¼ hYs2 ðiÞi ¼

s 1 X Y 2 ½ðn  1Þs þ i; Ls  s i¼1 s

(6)

T. Nguyen-Ky et al. / Biomedical Signal Processing and Control 5 (2010) 59–65

61

Fig. 1. (a) The collected EEG time series have five states of anaesthetic deep. (b) The integration on EEG samples of the time series.

and compute the value F*(s) over all segments by defining "

# 2Ls 1 X  log Fs ðnÞ : F ðsÞ ¼ exp Ls  s n¼1 

(7)

We define: P j ¼ log

F  ðsÞ ; s

j ¼ A; B; C; D; E:

(8)

Fig. 3 presents the values of Pj for five states A, B, C, D and E. We have PA = 114.47, PB = 92.3175, PC = 67.9175, PD = 64.3995 and PE = 59.9935. We can classify the five states of anaesthesia using the above values from Pj. 2.2.2. Separate fluctuation function F(s) in states C and D As shown in Fig. 2, the F(s) values from the DFA fluctuation function in states C and D are close to each other. In order to separate them, we redefine Eq. (5) as: "

Fig. 2. The DFA fluctuation function F(s) versus s.

L 1 X F MDFA ðsÞ ¼ F 2 ðsÞ q j Ls s¼3

#1=2 ;

q j > 0; j ¼ A; B; C; D; E:

(9)

T. Nguyen-Ky et al. / Biomedical Signal Processing and Control 5 (2010) 59–65

62

Fig. 4. Using modified DFA to separate states C and D.

Fig. 3. Using modified DFA to classify five states.

where parameter qj is introduced to change the amplitude of FMDFA(s). We denote qj as: qj ¼



14 2

if if

P j  PO ; P j > PO

PO ¼ 66:

(10)

The values of qj are determined by several facts below. Firstly, states C and D have PC = 67.9175, PD = 64.3995. These values are symmetrical through value 66. We choose PO = 66 as the break point to separate states C and D. Secondly, in Fig. 2, the fluctuation functions of F(s) in states A, B and C already separate the range of s  40, and the values of Pj in states A, B and C are higher than PO. We, therefore, choose qA = qB = qC = 2 as the original values in Eq. (5) of DFA method. Thirdly, in Eq. (9), the amplitude of FMDFA(s) decreases when the values of qj increase. We choose qD = qE = 14 as they provide the best value for selecting box size values and the range of FMDFA(s) as discussed in Section 3. They are, therefore, the optimal results based on our experiments. In this section, we define Eqs. (6)–(8) to classify the five states of the DoA in Fig. 3. By applying Eqs. (9) and (10), we can distinguish five states: very deep anaesthesia (state A), deep anaesthesia (state B), moderate anaesthesia (state C), light anaesthesia (state D) and awake (state E) as shown in Fig. 4.

is defined as ˆ e ¼ FðsÞ  FðsÞ:

(11)

The value of the best fit linear estimate is the sum of the squared distances. The linear equation can be written as ˆ ¼ a  s þ b: FðsÞ

(12)

where a is the slope of the line and b is called the F(s)-intercept. To find the derivatives of the Eq. (12) with respect to a and b, and let the derivatives equal to zero, we have: P a¼

¯ sFðsÞ  ns¯FðsÞ ; P 2 s  ns¯2

¯  as¯: b ¼ FðsÞ

¯ and s¯ are mean of F(s) and s. The values of a and b represent the FðsÞ straight line with the minimum sum of the squared distances. To set the range of FMDFA(s) from 0 to 200 and the range of the box size values from 0 to 50, we have the line Z: Z(s) = 4s + 200 in Fig. 5. This line is used to divide the zones of five states from very deep anaesthesia to awake state, corresponding to the DoA and box sizes. When the line Z is crossing every FˆMDFA ðsÞ line (BIS = 20, BIS = 40, BIS = 60, BIS = 80), we have five zones as follows:

3. State identification In Fig. 4, we can see that F(s) has an almost linear relationship with states A, B, C, D and E in the ranges of: 3  sA  20, 3  sB  25, 3  sC  35, 3  sD  45 and 3  sE  55, respectively; and their slope increases from the state E to the state A. These characteristics are used for the purpose of setting up the ranges of FMDFA and box size values. We assume that F(s) in the cases of BIS = 20, BIS = 40, BIS = 60 and BIS = 80 are lines in small ranges of s, 3  sA  60, and use the linear regression method to draw the lines ˆ from the data points F(s) in Fig. 5. These lines are used to divide FðsÞ five zones corresponding with five states from very deep anaesthesia to awake state. These results are used in the next step when we monitor the depth of anaesthesia by the Lagrange method. Linear regression is a technique that determines the best fit linear equation to a set of data points in terms of minimizing the ˆ sum of the squared distances between the line FðsÞ and the data points F(s) [17]. The distance from each point to its linear estimate

(13)

Fig. 5. The ranges of box sizes and FMDFA(s).

T. Nguyen-Ky et al. / Biomedical Signal Processing and Control 5 (2010) 59–65

    

Zone Zone Zone Zone Zone

1 2 3 4 5

has has has has has

the the the the the

space space space space space

formed formed formed formed formed

by by by by by

lines lines lines lines lines

63

s = 42, F(s) = 0 and Z(s), s = 35, F(s) = 30 and Z(s), s = 25, F(s) = 60 and Z(s), s = 12, F(s) = 100 and Z(s), s = 0, F(s) = 150 and Z(s).

In zone 1, line Z has FMDFA(s) values from 0 to 30 when s values decrease from 50 to 42. We define the range of FMDFA(s) (0–30) as range 1 to monitor the awake state. In zone 2, line Z has FMDFA(s) values from 30 to 60 when s values decrease from 42 to 35. We define the range of FMDFA(s) (30–60) as range 2 to monitor the light anaesthesia state. In zone 3, the line Z has the FMDFA(s) values from 60 to 100 when s values decrease from 35 to 25. We define the range of FMDFA(s) (60–100) as range 3 to monitor the moderate anaesthesia state. In zone 4, the line Z has FMDFA(s) values from 100 to 150 when s values decrease from 25 to 12. We define the range of FMDFA(s) (100–150) as range 4 to monitor the deep anaesthesia state. Finally in zone 5, the line Z has FMDFA(s) values from 150 to 200 when s values decrease from 12 to 0. We define the range of FMDFA(s) (150–200) as range 5 to monitor the very deep anaesthesia state in Fig. 5. We chose Z(s) depending on two parameters: box size values and FMDFA ranges. (a) Box size values: As shown in Fig. 1b, the values of Y(i) in five states have different fluctuations. Y(i) in state A has the fluctuation higher than those in states B, C, D and E. Similarly, Y(i) in state B has the highest fluctuation among those in states C, D and E; Y(i) in state C has the highest fluctuation among those in states D and E; and Y(i) in state D has the fluctuation higher than that in state E. With the Y(i) value having the higher fluctuation, we choose the small box size; and vice versa. Therefore, we decide to choose different box sizes for each state. We choose box sizes using the following rules: sA  sB  sC  sD  sE, with sA, sB, sC, sD and sE as box sizes for the states A, B, C, D and E. In Fig. 4, we can see that F(s) has an almost linear relationship with states A, B, C, D and E in the ranges of: 3  sA  20, 3  sB  25, 3  sC  35, 3  sD  45 and 3  sE  55, respectively; and their slope increases from the state E to the state A. This characteristic is used for the purpose of setting up the ranges of FMDFA and box size values. (b) FMDFA range: In order to expand the ranges of anaesthesia states in MDFA method, the ranges of FMDFA are chosen as below: awake state (0–30), light anaesthesia (30–60), moderate anaesthesia (60– 100), deep anaesthesia (100–150) and very deep anaesthesia (150–200). We finally obtain a relationship of range s and range of FMDFA as Z(s) = 4s + 200 as shown in Table 1 for the five zones.

4. Monitor the depth of anaesthesia using Lagrange method Let smaxT and sminT be the maximum and minimum box sizes in zone T (T = 1, 2, 3, 4, 5), FmaxT and FminT be the maximum and minimum of FMDFA, and sT be the value of the box size in the range

Fig. 6. The FMDFA(s) of EEG signal (BIS = 50) in zone 3.

of sminT < sT  smaxT. We have FMDFA value in the range of FminT  FMDFA < FmaxT. We propose to monitor DoA using Lagrange method [18]. We compute the smallest distance from data point F(s) to the Z(s) = 4s + 200 line using Lagrange method. For example, we consider FMDFA(s) of EEG signals in state C (BIS = 50) in Fig. 6. In this case, we have T = 3. In the range of (smin 3, smax 3), we have ten corresponding values for FMDFA(s) in the range (FminT, FmaxT). We need to find out the best value of FMDFA(s) for DoA. That is the value having the smallest distance from line Z. In Fig. 7, we choose one random point K(sK, FK) in the ellipse in Fig. 6, and compute the distance d from point K to line Z. The minimum d value is: minimize Subject to

dðsK Þ sminT < sK  smaxT ; K ¼ 1; 2; . . . ; 50 F minT  F K < F maxT :

(14)

Let e(sK) = DF = jFK  Z(sK)j be the deviation of FK and u be the angle between d and DF. We see that u equals the angle between line Z and line s = 0 (Fig. 7). The distance d is computed by d = DF cos u, with 0  u  90 and u = const. The optimization problem in Eq. (14) becomes: minimize Subject to

eðsK Þ sminT < sK  smaxT ; ZðsK Þ  F maxT ; K ¼ 1; 2; . . . ; 50:

(15)

The Lagrange of the problem in Eq. (15) is defined as LðsK ; lÞ ¼ eðsK Þ þ

50 X

lK ½ZðsK Þ  F maxT ;

(16)

K¼1

where lK is the Lagrange multiplier associated with the Kth inequality constraint Z(sK)  FmaxT. The optimization variable sK is called the primal variable. g(lK) is defined as the minimum value of the Lagrange over sK: gðlK Þ ¼ inf LðsK ; lÞ:

(17)

sK

Table 1 Define zone and range.

Zone Range s Range Z(s) BIS values

Awake

Light anaesthesia

Moderate anaesthesia

Deep anaesthesia

Very deep anaesthesia

1 50–42 0–30 80–100

2 42–35 30–60 80–60

3 35–25 60–100 60–40

4 25–12 100–150 40–20

5 12–0 150–200 20–10

T. Nguyen-Ky et al. / Biomedical Signal Processing and Control 5 (2010) 59–65

64

Table 2 The DoA values.

sK Z(sK) FK

Fig. 7. Choosing the optimum box size s.

The result is given by the optimality Karush–Kuhn–Tucker conditions that any primal-dual solution must satisfy [18]:

@L ¼0 @sK lK ½ZðsK Þ  F maxT  ¼ 0;

(18)

lK  0; F maxT  ZðsK Þ:

The above procedure is implemented in four steps to monitor the DoA in algorithm: Algorithm

Monitoring the DoA.

 Define the state using modified DFA. Select zone T and range s in Table 1.  Use modified DFA to compute the FK value with the box size sK.  Use Lagrange method to compute FK. Output FK as the DoA. The diagram for monitoring the DoA is shown in Fig. 8. The states of input EEG signals are defined by the modified DFA in Section 2.2. After identifying the states (from very deep anaesthesia to awake), we select zone T and range s in Table 1. We use the Lagrange method to compute the box size value s. With value s, we use the MDFA algorithms to compute the FK. The error feedback between the FK and Z(s) is used to adjust the box size s to the optimum value. With the optimum value s, we have the DoA value to equal the FK value. Table 2 shows the best values sK, FK of states A, B, C, D and E based on our simulation results. In zone T = 1, we compute the minimum of e(sK) in the range of s = (42–50) of

Fig. 8. Diagram for monitoring the DOA.

State E

State D

State C

State B

State A

44 24 21.514

42 32 33.356

32 72 76.747

22 112 111.598

11 156 147.848

state E. The DoA in state E is FK = 21.514, corresponding with sK = 44. Similarly, in zones T = 2, 3, 4, 5, we compute the minimum of e(sK) in the ranges of s = (35–41) of state D, s = (25–34) of state C, s = (13–24) of state B and s = (1–12) of state A. We have the DoA in state D (FK = 33.356), corresponding with sK = 42; state C (FK = 76.747), corresponding with sK = 32; state B (FK = 111.598), corresponding with sK = 22; and state A (FK = 147.848), corresponding with sK = 11. Fig. 9 shows the values of BIS and FK. In Fig. 9, we can see that the range of FK is extended in the states of moderate anaesthesia, deep anaesthesia and very deep anaesthesia. All consciousness monitoring technologies require processing time; and none is technically ‘‘real-time’’. There are two main tasks in obtaining a depth of anaesthesia value. The first one is to capture EEG data from a patient. The second one is to use signal processing techniques to process the data. In BIS Index, a patient’s BIS value is calculated using about 15– 30 s of captured EEG data. It is possible that a 10–15 s delay occurs between an EEG segment and their corresponding BIS value if the clinical situation changes rapidly. The BIS Index is updated every one second. The BIS trend is updated every 10 s. Time delay of index calculation for the Cerebral State Index, the Bispectral Index, and the Narcotrend Index was studied by Pilge et al. [19]. In this paper, we apply the Lagrange method for monitoring the DoA on real signals. The Lagrange algorithm is powerful in solving the original minimization problem in Eq. (14). The basic idea of the Lagrange method is to take the constraints of Eq. (14) into account by augmenting the objective function with a weighted sum of the constraint functions. With the Lagrange method, it is simpler to solve the dual problem rather than the primal one in some cases. Many traditional optimization methods are based on gradient descent algorithms, which suffer from slow convergence and sensitivity to the algorithm initialization and step-size selection. Comparing our method with BIS Index, we observed that: (a) Our method significantly reduces computational complexity The entire algorithm for the calculation of BIS Index has not been made public. The algorithm was only partially described in one review by Dr. Rampil [20]. According to this review, BIS Index is calculated from following four parameters:  Burst Suppression Ratio (BSR) by time domain analysis. BSR is a time domain EEG parameter to recognize those periods longer than 0.50 s, during which the EEG voltage does not exceed approximately 5.0 mV.

Fig. 9. Comparison between BIS and Fk.

T. Nguyen-Ky et al. / Biomedical Signal Processing and Control 5 (2010) 59–65

 Quazi suppression index also by time domain analysis. It was designed to detect burst suppression in the presence of wandering baseline voltage.  Relative b ratio (log P30–47Hz/P11–20 Hz) by power spectral analysis, with Px–y is the sum of the spectral power in the band extending from frequency x to y.  SynchFastSlow (log B0.5–47 Hz/P40–47 Hz) by bispectral analysis, with Bx–y is the sum of the bispectrum activity in the area subtended from frequency x to y. Comparing with our method, we found that the computation of BIS value is very complex and requires more computation time than our method. EEG signals are essentially non-stationary time series that include artifacts, mainly due to muscle activities. The BIS Index is calculated focusing on the frequency domain. Although BIS Index requires the values of BSR and Quazi in time domain, as defined above, BSR value does not clearly reflect the variations of EEG signals in time domain. On the other hand, the DFA method has been established as an important fractal analysis technique for the detection of longrange temporal correlations in non-stationary time series. (b) Our method produces faster reaction to transients in patients’ consciousness levels The integration signals with higher amplitudes and lower fluctuations which correspond to less brain activities (state A in Fig. 1a) have resulted in a more wavy output Y(i) in Fig. 1b. The integration signals with high fluctuations and low amplitudes which correspond to high brain activities (state E in Fig. 1a) have smoothed output Y(i) in Fig. 1b. The rate of fluctuation Y(i) decreases when the states of a patient’s brain are changed from anaesthetized state to awaked situation. (c) Our method extends the ranges of DoA values which are very significant for monitoring the DoA accurately With a BIS Index value of less than 60, a patient has an extremely low probability of consciousness. BIS Index values lower than 40 signify a greater effect of the anaesthetic on the EEG. With low BIS values, the degree of EEG suppression is the primary determinant of the BIS value. Prospective clinical trials have demonstrated that maintaining BIS Index values in the range of 40–60 ensures adequate hypnotic effect during general anaesthesia while improving the recovery process. During this case, clinicians adjust the dose of medication to keep BIS values in the range of 40–60 [3]. Therefore, the DoA values in this range are very important. Comparing the range in our method and BIS Index, we can see that the ranges in our method (MDFA) are expanded in the states of moderate anaesthesia (60 < FMDFA < 100, 40 < BIS < 60), deep anaesthesia (100 < FMDFA < 150, 20 < BIS < 40) and very deep anaesthesia (150 < FMDFA < 200, 0 < BIS < 20) in Table 1. The expanded ranges will help clinicians having more DoA values during surgery. These results are very significant for monitoring the DoA accurately as the clinical perspective needs more attention and information about the DoA in these states. 5. Conclusion We have proposed a new method to monitor the DoA accurately using the modified detrended fluctuation algorithm method. In this approach, we first use MDFA to classify raw EEG signals into awake, light anaesthesia, moderate anaesthesia, deep anaesthesia and very deep anaesthesia states. Then we separate the diagrams of states C and D by using parameter qj in MDFA. By comparing

65

with the simple DFA method, the MDFA method can clearly identify one state from the other states. Using linear regression and Lagrange method, we build five zones from the very deep anaesthesia to the awake state, corresponding to different box sizes. Comparing the range in the BIS method to the range in the MDFA method, we find the MDFA range is expanded in the states of moderate anaesthesia, deep anaesthesia and very deep anaesthesia. These results are very significant for monitoring the DoA accurately as the clinical perspective needs more attention and information about the DoA in these states. Simulation results demonstrate that this new technique monitors the DoA values in all anaesthesia states accurately. The MDFA method proposed in this paper can provide a more accurate result for DoA monitoring.

Acknowledgments This research is supported by Australia Research Council (ARC) Discovery Program grant: DP0665216. The authors greatly appreciate the valuable comments and suggestions received from the reviewers of this paper.

References [1] M. Agarwal, R. Grifths, Monitoring the depth of anaesthesia, Applied Physics (2004) 343–344. [2] J. Bruhn, P.S. Myles, R. Sneyd, M.M.R.F. Struys, Depth of anaesthesia monitoring: what’s available, what’s validated and what’s next? British Journal of Anaesthesia 97 (June (1)) (2006) 85–94. [3] S.D. Kelly, Monitoring consciousness using the Bispectrum Index during anaesthesia, A pocket guide for clinicians, Aspect medical systems, 2007. [4] J.C. Drummond, Monitoring depth of anaesthesia, Anaesthesiology 93 (September (3)) (2000) 876–882. [5] I. Rezek, S.J. Roberts, R. Conrart, Increasing the depth of anaesthesia assessment, IEEE Engineering in Medicine and Biology Magazine (April 2007) 64–73. [6] X.-S. Zhang, R.J. Roy, E.W. Jensen, EEG complexity as a measure of depth of anaesthesia for patients, IEEE Transactions on Biomedical Engineering 48 (December (12)) (2001) 1424–1433. [7] J.W. Huang, Y.-Y. Lu, R.J. Roy, Depth of anaesthesia estimation and control, IEEE Transactions on Biomedical Engineering 46 (January (1)) (1999) 71–81. [8] O. Ortolani, A. Conti, A. Di Filippo, C. Adembri, E. Moraldi, A. Evangelisti, M. Maggini, S.J. Roberts, EEG signal processing in anaesthesia. Use of a neural network technique for monitoring depth of anaesthesia, British Journal of Anaesthesia 88 (January (5)) (2002) 644–648. [9] T. Zikov, S. Bibian, G.A. Dumont, M. Huzmezan, C.R. Ries, Quantifying cortical activity during general anaesthesia using wavelet analysis, IEEE Transactions on Biomedical Engineering 53 (April (4)) (2006) 71–81. [10] R. Ferents, T. Lipping, A. Anier, V. Ja¨ntti, S. Melto, S. Hovilehto, Comparison of entropy and complexity measures for the assessment of depth of a sedation, IEEE Transactions on Biomedical Engineering 53 (June (6)) (2006) 1067–1077. [11] M. Jospin, P. Caminal, E.W. Jensen, H. Litvan, M. Vallverdu´, M.M.R.F. Struys, H.E.M. Vereecke, D.T. Kaplan, Detrended fluctuation analysis of eeg as a measure of depth of anaesthesia, IEEE Transactions on Biomedical Engineering 54 (May (5)) (2007) 840–846. [12] P. Gifani, H.R. Rabiee, M.H. Hashemi, P. Taslimi, M. Ghabari, Optimal frac-scaling analysis of human eeg dynamic for depth of anaesthesia quantification, Journal of the Franklin Institute 344 (May–July (3–4)) (2007) 212–229. [13] C.K. Peng, J. Mietus, J.M. Hausdorff, S. Havlin, H.E. Stanley, A.L. Goldberger, Long range anticorrelations and non-gaussian behavior of the heartbeat, Physical Review Letters 70 (March) (1993) 1343–1346. [14] J.-M. Bardet, I. Kammoun, Asymptotic properties of the detrended fluctuation analysis of long range dependence processes, IEEE Transactions on Information Theory (July 2006) 1–9. [15] J.W. Kantelhardt, et al., Multifractal detrended fluctuation analysis of nonstationary time series, Physica A 316 (December (1–4)) (2002) 87–114. [16] A. Bunde, et al., Correlated and uncorrelated regions in heart-rate fluctuations during sleep, Physical Review Letters 85 (2000) 3736–3739. [17] R.I. Levin, D.S. Rubin, Statistics for management, Prentice-Hall, 1991. [18] S. Boyd, L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004. [19] S. Pilge, et al., Time delay of index calculation, Anesthesiology 104 (2006) 488– 494. [20] I.J. Rampil, A primer for EEG signal processing in anesthesia, Anesthesiology 89 (4) (1998) 981–1001.