- Email: [email protected]
Assessment of anaesthetic depth by clustering analysis and autoregressive modelling of electroencephalograms
Computer Methods and Programs in Biomedicine, 34 (1991) 125-138 © 1991 Elsevier Science Publishers B.V. 0169-2607/91/$03.50
125
COMMET 01137
Assessment of anaesthetic depth by clustering analysis and autoregressive modelling of electroencephalograms C . E . T h o m s e n a, A. R o s e n f a l c k 1 a n d K . N o r r e g a a r d C h r i s t e n s e n 2 1 Department of Medical Informatics and Image Analysis, Aalborg University, DK 9220 Aalborg, Denmark, and 2 Department of Anaesthesia, Aalborg Hospital, DK 9000 Aalborg, Denmark
The brain activity electroencephalogram(EEG) was recorded from 30 healthy women scheduled for hysterectomy. The patients were anaesthetized with isoflurane, halothane or etomidate/fentanyl. A multiparametric method was used for extraction of amplitude and frequency information from the EEG. The method applied autoregressive modelling of the signal, segmented in 2 s fixed intervals. The features from the EEG segments were used for learning and for classification. The learning process was unsupervised and hierarchical clustering analysis was used to construct a learning set of EEG amplitude-frequency patterns for each of the three anaesthetic drugs. These EEG patterns were assigned to a colour code corresponding to similar clinical states. A common learning set could be used for all patients anaesthetized with the same drug. The classification process could be performed on-line and the results were displayed in a class probability histogram. This histogram reflected in all patients the depth of anaesthesia, when the concentration of the anaesthetic agent was adjusted either based on clinical signs or according to the protocol. This uniform display, where colours in a class probability histogram indicate the depth of anaesthesia, may in the future serve as on-line advice for the administration of anaesthetics. A comparison of multiparametric with single parametric methods, based on calculation of median, spectral edge and peak frequencies, questions the reliability of the single parametric methods in monitoring anaesthetic depth. Anesthesia, depth; Unsupervised learning; Hierarchical clustering; Isoflurane; Etomidate; Halothane; Human interface; On-line computer methods; Electroencephalogram amplitude-frequency patterns
1. Introduction C o m p u t e r i z e d analysis of the electroencephalogram ( E E G ) proved to b e valuable for long time m o n i t o r i n g of patients with head i n j u r y a n d during surgery [1]. T h e r e are two m a j o r reasons for using o n - l i n e computerized methods: (i) the imm e n s e a m o u n t of d a t a can be reduced in a n objective way, a n d (ii) signal analysis m e t h o d s c a n extract i n f o r m a t i o n on the a m p l i t u d e - f r e q u e n c y
Correspondence: C.E. Thomsen, Department of Medical Informatics and Image Analysis,Aalborg University, Frederik Bajers Vej 7 D, DK 9220 Aalborg, Denmark.
p a t t e r n of the E E G , which c a n n o t be detected b y visual i n s p e c t i o n of the E E G trace [2]. The use of E E G m o n i t o r i n g as a guide for the d e p t h of anaesthesia has b e e n controversial. D u r ing surgery m a n y factors m a y i n f l u e n c e the E E G a n d it becomes difficult to detect changes that can solely be related to the d e p t h of anaesthesia, because fluctuations in the arterial c a r b o n dioxide tension, b o d y temperature, p a i n stimuli a n d cardiovascular changes are present. I n emergency cases hypoxia will interfere with changes caused b y anaesthetics [1]. C u r a r i z a t i o n a n d the a d v e n t of new anaesthetics leave the anaesthetist with limited i n d i c a t i o n s of the anaesthetic depth. T o o light anaesthesia
126
might lead to increased need for analgesic and sedative drugs in the post-surgical period and to increased post-surgical stress in general. Awareness of pain or recall of other events during general anaesthesia has become a current medico-legal problem in the U.S.A. and in the U.K. [3]. The cerebral function analysing monitor (CFAM) has proven to be an appropriate EEG monitor in clinical practice [1]. The EEG activity is displayed in the traditional frequency bands (delta, theta, alpha and beta) together with the mean peak to peak value and the 10th and 90th percentiles of amplitude. The use of EEG monitoring systems still requires an extensive experience in interpretation of changes in amplitude and frequency content of the EEG. Efforts have been made to extract single parameters of the EEG activity to optimize the administration of anaesthetics. The median frequency (MF) of the power spectrum has been used as a single parameter characterizing the degree of central nervous system depression [4,5]. Another parameter is the spectral edge frequency (SEF) introduced by Rampil [6]. The development of the present method aimed towards the application of multiparametric techniques, to extract specific classes of amplitudefrequency patterns related to different levels of anaesthesia [7]. By means of this analysis, uniform information concerning the depth of anaesthesia could be obtained for three different anaesthetic techniques (halothane, isoflurane and etomidate/ fentanyl) despite the agent specificity of electroencephalographic changes.
2. Patients We investigated 30 healthy women (29-55 years old, body weight not exceeding 10% of ideal) undergoing simple hysterectomy. There was no history of neurological disorders or administration of analgesics and tranquillizers in the weeks before operation. The day before operation, informed consent was obtained. The investigation was approved by the local Ethics Committee. Three groups with ten patients in each group were investigated:
Group A was anaesthetized with halothane. Induction was made with thiopentone. For tracheal intubation a priming dose of gallamine and suxamethonium was administered. Anaesthesia was maintained with halothane 0.75-1.5 vol% in oxygen-nitrous oxide 1 : 2. In group B induction of anaesthesia was performed as described for group A, and anaesthesia was maintained with isoflurane 1.0-2.5 vol% in oxygen-nitrous oxide 1 : 2. In group C anaesthesia was induced and maintained with etomidate and fentanyl. Maintenance infusion was 0.5-1.0 m g / m i n etomidate and 1.63.5 m g / m i n fentanyl. All patients were mechanically ventilated by MC 801 Dameca anaesthetic ventilator. The duration of anaesthesia was 70-175 min. Monitoring of EEG from both hemispheres (P3-Fpa and P4-Fp2) started before induction of anaesthesia. In all patients temperature, end-tidal carbon dioxide concentration and mean arterial blood pressure were kept within physiological limits. All information was stored on a Racal Store 8 frequency modulated tape recorder. Further analysis of signals was performed off-line. From induction of anaesthesia the mechanical ventilator was adjusted to an end-tidal carbon dioxide concentration of 5.0-5.5 vol% corresponding to normoventilation. 20-25 rain after start of surgery, when the clinical condition was considered stable, a slight hyperventilation was induced corresponding to an end-tidal carbon dioxide concentration of 4.0-4.5 vol%. After 20 min hyperventilation was discontinued. Approximately 30 min after normoventilation was reestablished, the dosage of anaesthetics was increased by approximately 50% to study the E E G at a deeper anaesthetic level. After 20 min the maintenance concentration was reestablished.
3. Signal processing The signals were analysed off-line using an Intel 310 computer system equipped with a high resolution colour graphic display (1024 x 512 pixels in 256 colours) and analog to digital converter (12 bit resolution).
127 The EEG was filtered with a 25 Hz 4th order antialiasing filter, sampled at 100 Hz and preemphasized digitally with a first order high-pass filter at 4.2 Hz (Fig. 1). This filter gives approximately 10 d B / d e c a d e preemphasizing up to 25 Hz. The digitized EEG signal was divided into 2 s fixed segments using rectangular window function. The length of the window must be greater than 1 / 4 of the wavelength of the lowest frequency, and short enough to ensure stationarity within the window. A window length of 1-5 s is suitable, but 2 s preferable [8].
3.1. Feature extraction From each segment 11 features were extracted (ten normalized autocorrelation coefficients (ACCs) and RMS amplitude). These features were transformed to amplitude or power spectra by autoregressive modelling using the Durbin algorithm [9]. ACCs: N-i--1
r(i)=
~
x(n)x(n+i)
f o r i = 0 . . . . . 10
n=0
10o
o
y
g~
/ ,/
I02L 10 I
/
Normalized ACCs:
R(i)=r(i)/r(O)
for i = 1', . . . . 10
Amplitude: RMS = c r ( ( ~
c = scaling constant
During experiments ACCs were compared to reflection coefficients [9,10] and prediction coefficients [9] combined with either the mean squared power or the RMS amplitude. Prediction coefficients, reflection coefficients and ACCs contain exactly the same information, but encoded quite differently. The normalized ACCs combined with the RMS amplitude were chosen as features because they proved to be the least sensitive set of features with respect to interpatient variability.
4. Learning/clustering analysis Based on feature vectors extracted from a representative set of segments the learning was performed unsupervised by hierarchical clustering. The aim of the unsupervised approach was to find a natural grouping in the EEG pattern (not hmited by the definition of delta, theta, alpha and beta activities). But how do we define 'natural' and how can we express that data units within a cluster are more similar than data units taken from two different clusters? These questions can be expressed as two separate problems: - how is similarity between data units measured? and - how can the partition of data units in clusters be evaluated?
4.1. Similarity measures 100
101
10z
Frequency lHz] Fig. 1. Frequency response of the high-pass filter used to preemphasize the EEG signal. In the range from 2-20 Hz the filter characteristic is similar to that used by Prior and Maynard [1] in the CFAM. In the range below 2 Hz the filter has a slope of 20 dB/decade, whereas Prior and Maynard use filtering with 60 riB/decade.
An obvious choice of similarity (or dissimilarity) measure between data units is the Euclidean or geometric distance (1). A generalization of the Euclidean distance is the Mahalanobis distance (2). The Mahalanobis distance includes the covariance matrix and therefore compensates for the correlation between components of the feature vector.
128 The optimal clustering is the one which minimizes L [ll].
Euchdean distance:
de = ~ ( x - _ ~ ) r ( x - ~ ) =/l_x-~ll
(1) 4.3. Hierarchical clustering
Mahalanobis distance:
dm = ~/(x - ~) T~- l(X --_/£)
(2)
4.2. Criteria for clustering An essential problem is to evaluate the quality of the clustering analysis, and to define the optimal division of the n design units into an unknown number c of groups. An often used criterion function is the error sum of squares (or the minimum variance method). Let nj be the number of data units and x # the k-th sample both for the j-th cluster, and let ~j be the mean value for that cluster. Then the error sum of squares is defined using the Euclidean distance:
nj Je = E E (X--jk--I£j)T(x--jk--~j) j~l k~l c
(3/
A matrix of similarities between all feature vectors is the usual base for hierarchical clustering analysis. However this matrix contains n * ( n + 1 ) / 2 elements, where n denotes the total number of data units. A considerable reduction in m e m o r y could be obtained when, for each input vector, only the identifier and the distance to the nearest vector is stored. This matrix only contains 2n elements. The result from hierarchical clustering can be represented in a dendrogram (Fig. 2). The branches in the bottom represent one entity or one data unit, while the root at the top represents the entire collection of entities. The attention is here directed to agglomerative hierarchical clustering building the tree from branches to the root. Notice that in hierarchical methods clusters are nested. When two clusters are merged, they are joined permanently and become a building block for later merges. Herein lie both the strength and the weakness of hierarchical clustering analysis.
0 m
I 40 .~
60~
i00 Fig. 2. A dendrogram is the graphic description of the results from hierarchical clustering analysis. The Figure illustrates how clusters are merged as a function of similarity (ordinate). The spectra shown are mean values of the spectra for merging branches in the tree
from the subtrees below. This dendrogram is only part of the dendrogram for the group of patients anaesthetized with isoflurane.
129
The general procedure for agglomerative hierarchical methods can be described in Pseudo Pascal:
a)
PROCEDURE HIERARCHICALCLUSTERING BEGIN FOR i:=l TO n_unitDO c := n_umi; WHILE c > 1 DO BEGIN Find the two clusters,which have a minimum inter-distance.Let's call them •i and t)j; Merge f/i and Dj in a new cluster Di; Delete Gj; Update similaritymatrix for D~ Update similaritymatrix for clusterswith Di or ~j as nearestneighbours c := c-1; END; END;
b)
The number of data units is denoted n_uni t and ~2, is the i-th cluster. In the procedure the term 'distance' is used. In the previous section we discussed the choice of similarity measure. Between which points within the two clusters are we measuring, and shall the distance be weighted according to the size of the clusters? Three often used cluster distances are:
dmin(~'~i,~Qj)=
min
ll_x-_x'll
dma×(#,, ~j)=
max
I I x - x'll
dmean(~i,
~))=
Fig. 3. T o illustrate t h r e e d i f f e r e n t s i m i l a r i t y m e a s u r e s f o r c l u s t e r i n g of E E G p a r a m e t e r s f o r the c o n s t r u c t i o n o f a c o m m o n l e a r n i n g set f o r a specific a n a e s t h e t i c d r u g . E l e m e n t s f r o m c l u s t e r Y~l a n d [22 a r e m a r k e d w i t h (o) a n d (x). T h e m e a n values f o r the t w o clusters a r e m a r k e d w i t h ( O ) a n d (X). I n ( a ) , ( b ) a n d ( c ) the m i n i m u m , t h e m a x i m u m a n d the m e a n o c c u r r i n g d i s t a n c e b e t w e e n e l e m e n t s f r o m the c l u s t e r s ~ a n d 122 are u s e d as i n t e r - c l u s t e r d i s t a n c e .
II_bti--~j]l
dmin(#i,
~j) is the minimum occurring distance from data units in cluster 12~ to data units in cluster 12j (Fig. 3a). Analogous dmax (12i, 9j) is the maximum occurring distance between samples within the two clusters (Fig. 3b). In the third method the distance is measured between mean values or centroids of clusters (Fig. 3c). Combining the minimum variance methods and the centroid methods, it was possible to obtain stepwise optimization of the error sum of squares Je (3) by using a weighted dmean(Oi, ~2j) [12,13]: d~ean(~i,
c)
~j)= ~/(rl i * rlj)/(rl i + tl2)
*
[[_~i --_PL, H
(4) n i and nj are the number of data units within cluster 12i and f2j. This measure could be
interpreted as the square root of the increase in the error sum of squares if cluster/2, and Dj were merged. Even though the clustering was stepwise optimal, there was no quarantee for global minimization of this criteria function.
5. Construction of learning sets
During learning we selected a set of 504 E E G segments from each patient representing typical EEG patterns from different levels of anaesthesia. During hierarchical clustering segments were merged together - - step by step starting with the 504 segments ending with one. The modified Euclidean distance (4) between feature vectors was used as similarity measure. The clustering sequence was represented graphically in a dendro-
130
of EEG patterns in all patients anaesthetized with that particular agent.
6. Classification
Fig. 4. To illustrate 'linked' or stepwise clustering analysis, when data from five or ten patients were included in the learning process. Below: Clustering sequences from each patient, starting with 504 feature vectors. The arrow indicates the level where 50 classes from each dendrogram are selected and merged in a dataset representing c o m m o n E E G patterns. Above: The c o m m o n set of patterns was used in an additional clustering analysis to define a c o m m o n set of classes (to be used as a learning set during classification), specific for one anaesthetic agent.
gram. The dendrogram from hierarchical clustering could be inspected interactively. This inspection was guided by displaying mean spectra from subtrees below selected nodes (Fig. 2). During inspection special attention was paid to nodes where the distance increased (i.e. similarity decreased) rapidly. Using this method a set of classes typical for this patient could be defined. Classes containing equal or less segments than the number of features were removed from the class definition. To generate a common learning set for a specific anaesthetic agent, data from up to ten patients were included in the clustering sequence using 'linked' clustering analysis (Fig. 4). From the dendrograms, for each of the patients anaesthetized with the same drug, we selected a level where 50 clusters were represented. The 50 clusters from each patient were merged and clustering analysis repeated. After this procedure 12-16 classes were defined to be used as learning set for classification
The probability p(~2~ I_x) for the object belonging to the class ~2~ given the input feature vector _x was calculated from Bayes theorem. Calculating p (12~ I_x) for all classes and choosing the class with the largest probability gave a minimum probability for misclassifications [11]. We assumed that the probability density function for each class had a multidimensional normal distribution and that all classes had the same a priori probability. The individual class centroids were determined by the mean value vectors it, and the size and shape by the pooled covariance matrix Z. The decision problem was expressed in a set of linear discriminant functions [11] g~(x), i = 1 . . . . . c. The classifier assigned the feature vector to the class 12~, if g , ( x ) > g j ( x ) for all j v~ i.
The set of linear discriminant functions used in the system were g~ (x_) = Wlr~x + WO~ --1
where: W 1 i = Y, #, and: W0, = - ( 1 / 2 ) ~ Z --
-
T
-I
~,
The results from the classification were presented as 'raw' class probability histograms displayed together with density spectral array (DSA) plots as a function of time (Fig. 5). DSA plots are the classical presentation of the amplitude frequency pattern of the EEG [14]. Each line represents one amplitude spectrum, where the
Fig. 5. Example of probability scores for EEG amplitude frequency classes derived from a patient anaesthetized with isoflurane. Vertical axis is time and the total duration is 2 h and 50 min. Events are marked along the time axis at the left: (T) thiopentone; (x), (X) start and stop of anaesthesia; (o) onset of surgery; ( 1"), ( - ) start and stop of hyperventilation; ( > ), ( < ), increased and decreased induction of anaesthesia. Note that for this patient hyperventilation was induced before the change in induction of anaesthesia.
0~ 4-J
DSA spectrum
Class 13 histogram
Class 11+12+13 histogram
Class 4 histogram
Class 2+3+4+8 histogram
Class 6 histogram
Class 6+7+10 histogram
W-
1 -I ..m
I
--! P
i
I :
q
4 4
"t
P
iii~Li,ill i---!!
b X
ip
132 anaesthesia (e.g. s t a r t / s t o p of isoflurane induction and i n / d e c r e a s e d concentration of inspired isoflurane. For the example in Fig. 5 it was clear that class 6, 7 and 10 corresponded to increased anaeshetic depth, and that class 11, 12 and 13 occurred during induction of anaesthesia and in the awakening period, corresponding to light anaesthesia. Using this procedure classes or sums of classes were assigned colours according to the level of anaesthesia (Fig. 6 [7]) for each anaesthetic agent
frequency axis is the abscissa and the amplitude is expressed in a grey-scale, black denotes maximal amplitude. This presentation is used in Figs. 5 and 7 to support the evaluation of the class probability histogram. The DSA traces are globally scaled to give a true impression of amplitude changes during surgery. The class probability histograms for some of the classes and for summated classes are shown in black and white print. The occurrence of each class was then related to clinical statements, especially statements related to the depth of
-- "-----i~.:I--F:ISS 1 6
DROWSY (dark blue) ,
.
'
,
,
n
if: R S S
i2
VERY LIGHT (light b l u e )
N O R M A L
0
20
30
50
0
20
30
50
°I~"" C-L A S S ~
2
0
20
30
50
30
50
~ o b" -C L'A ~ S S ,~, 3, ,
"
--
-- ~
-
CLASS
8
[yellow}
6
' 26
3b
'
5b
o
2o
3o
50
30
5o
o
20
~o
26 3~
I 50
[purple}
0 ~0 30
----Pre.
.
' ~
3b
'
5b
o
.
.
2o
.
.
Fig. 6. Set of averaged amplitude spectra for classes from the group of patients anaesthetized with isoflurane. One to four classes were assigned to the same colour code because they all appear at the same depth of anaesthesia.
133 (blue = 'drowsy', light blue = 'very light anaesthesia', green ='light anaesthesia', yellow = ' n o r m a l surgical anaesthetic depth', purple = ' d e e p anaesthesia' and red = ' very deep anaesthesia'). The sets of classes and colour codings were used during on-line pattern recognition. Classes representing awake states, artifacts and movements were assigned to background colour (black on the screen and white on the paper print out). Even though the whole configuration of the amplitude spectrum was essential in classification of anaesthetic depth, the peak frequency was found to be one of the dominating single parameters in the learning set. For isoflurane a peak frequency of 10-12 Hz was found during light anaesthesia. Surgical depth of anaesthesia was characterized by a peak frequency of 4 - 6 Hz and deep levels of anaesthesia by a peak frequency of 1-3 Hz.
7. Results from monitoring of anaesthetic depth
The three learning sets of classes and colour codings, obtained for each of the three different anaesthetic drugs, were used during on-line classification. The information from the EEG analysis was displayed as a function of time [7]: (a) as density spectral array plots (DSA) (Fig. 7, left column). Each p l o t / l i n e represents the averaged amplitude spectrum during the last 10 s; (b) as class probability histograms (Fig. 7, right column). Each line represents the relative occurrence of any class in percent during the last 10 s. The class probability histograms were smoothed by a 2nd order low pass filter with two real poles. In all patients this histogram reflected the depth of anaesthesia when the concentration of the anaesthetic agent was adjusted either based on clinical signs or according to the protocol. During steady-state the class probability histogram often indicated an increased anaesthetic level, when pain stimuli were less, towards the end of surgery.
8. Results when outliers were excluded from the learning set
The disturbing effect of outliers in the learning sets was clearly demonstrated by clustering analysis from the three groups of patients. During isoflurane anaesthesia, all recordings were very 'clean' and contained practically no artifacts, resulting in a reliable set of classes, and a well defined anaesthetic scale. For the two other anaesthetic agents the sets of classes based on data from all ten patients gave some inconsistency in the class probability histogram, and a number of classes could not be assigned according to the anaesthetic depth. This problem was minimized by reducing the number of patients contributing to the learning set to five. The five patients were selected as those who responded most uniformly to surgical depth of anaesthesia. When these learning sets were used for all ten patients in each group, there was a better consistency and a more reliable configuration of the class probability histogram during online analysis. The improvement was not only related to patients inside the learning group, but was also observed for the five patients not included in the learning material. 8.1. Single parametric methods
Assuming that the anaesthetic scale used in the class probability histogram (Fig. 7, right column) was reliable we evaluated a number of single parametric methods against the present multiparametric method. The anaesthetic scale was converted to an analog scale (0-5) giving each anaesthetic depth from Fig. 6 a value: 'drowsy' = 0, ' v e r y light anaesthesia' = 1, 'light anaesthesia' = 2, 'normal surgical anaesthetic d e p t h ' = 3, 'deep anaesthesia' = 4, and ' very deep anaesthesia' = 5. From each 10 s epoch of the EEG, the median frequency (MF), the peak frequency (PF) and the spectral edge frequency (SEF) were calculated. The values of MF, SEF and PF were calculated for one patient solely and for the ten patients in group B together. The values for approximately 700 epochs in each patient were plotted against
30~U-% i v
134
" ~ , ' ~ - - .~¢v-'~ -,.',,* w' r l p ~ m " l F w
0
DSA spectrum 15 30[Hz]
0
Class prob. histogram 50 1001%1
T
t
<:
X
CLAS:~
CI3D
I I'JL~ :
135
the analogue scale for depth of anaesthesia, and the linear regression lines were calculated (Fig. 8). MF, SEF and PF frequencies decreased as expected with the depth of anaesthesia (Fig. 9). In evaluating the performance of these single EEG parameters one should keep in mind that the misclassification rate corresponded to the squared correlation coefficients. The scatter around the regression line for the median frequency was illustrated by selecting the values encountered in the range of the anaesthetic scale between 2.5 to 3.5. This corresponds to periods in the class probability histogram, where normal surgical anaesthetic depth (yellow) was scored by the multiparametric method. The median frequencies for one patient ranged from 5 to 10 Hz and for ten patients from 4 to 11 Hz (Fig. 9).
9. Discussion
If EEG monitoring shall serve as decision support during administration of anaesthetic drugs, a 'normal' population of patients has to be investigated during standard conditions. The aim of this study was to use a multiparametric method to analyse the EEG from 'normal' patients anaesthetized with different anaesthetic techniques and to facilitate the implementation of this method in an 'anaesthetic level monitoring system'. The multiparametric method used features from autoregressive modelling of the EEG signal, segmented in fixed time intervals. Before the feature extraction the signal was preemphasized or prewhitened by a filter (Fig. 1) similar to that used by Prior and Maynard [1]. Preemphasizing was done to reduce the redundant information in the signal caused by the inverse proportionality between the amplitude and frequency of EEG changes during anaesthesia [15,16].
The learning process was unsupervised and applied clustering analysis to define a learning-set of classes for three anaesthetic drugs. Each set of basic classes were assigned colours according to the level of anaesthesia. The l e a r n i n g s e t w a s c o m m o n for each anaesthetic agent but different for the three agents. It was used for classification of the EEG activity. The result of the classification was presented in class probability histograms and demonstrated a high degree of simplicity. To suppress short periods of disturbances from artifacts in the EEG trace the class probability histograms were smoothed by a 2nd order low pass filter with two real poles. The effect was that the histograms became more clear. At different levels of anaesthesia more than two classes were seldom represented. The simple display, where the same colour is related to the same level of anaesthesia for three different anaesthetic agents may support the anaesthesiologist. It has been suggested that single parametric methods: the median and the spectral edge frequencies can be used to control of the depth of anaesthesia. Stoeckel and Schwilden [4,17] have shown that the time course of the median frequency parallels the visual staging of anaesthesia and is closely related to the time course of the declining blood concentration of etomidate. They also found changes in the median frequency when the patient was anaesthetized with halothane and isoflurane but the variation with the depth of anaesthesia was different for the three drugs. Rampil and Matteo [6] found consistent changes in SEF related to haemodynamic responses during anaesthesia. White and Boyle [18] were not able to confirm these findings but found good correlation between propofol blood levels and EEG changes recorded by the CFAM [19]. A comparison of multiparametric and single parametric methods was attempted. From the re-
Fig. 7. Example of a colour copy of results from EEG analysis from a patient anaesthetized with isoflurane. Vertical axis is time and the total duration is 2 h and 20 rain. Events are marked along the time axis at the left; (T) thiopentone; (x), (X) start and stop of anaesthesia; (o) onset of surgery; ( > ) , ( < ) , increased and decreased induction of anaesthesia; (1"), ( - ) start and stop of hyperventilation. Left column is the DSA plot. Right column is the class probability histogram, colour coded as illustrated in Fig. 6. The trace above is a 2 s segment of the raw EEG signal, which is continuously displayed for control.
136
30
30
Ten paUents
Pa0ent 2 y : -I.5; ÷ 12.74, r2 =
y = -2,11x + 13.79, r2 = .88
.52
25 ~'o
2~-
~_.
20
~
15
o
o
o
o o
°
oo
o o ooo
o
o
o co ° o O
o o ¢
o ~~
o
~
o
o o
o
~
o qo~ooo
o
a
o
o
~
tO
oo o
o o
o
5
•5
l
1.5
2 2.5 3 Anaesthetic depth
y = -1.28x
Patient 2
3.5
+ 24.42,
r2
4
4.5
0
45
= .16
.5
oo
o
I
1.5
2 2.5 3 Anaesthetiedepth
"Ten pa~ents
4
+ 26.33,
r2
4.5
5
o y = -1.89x
o
= .38
o o
O
40
3.5
40
o
o
o
o
~o
o
O
o
o
35
~ o
°o o
o°
30
o
•~
oo o
o e
o
°
o
o
o
o #°°0
ff
o I
30
~
25
o
0
0
Oo
20
o 15} 0
I .5
oo
o
oo
2
o
~
°
:
o
o
25~,~._ o 8 ~ O o
o
o o%
°
~'
o
8
35
1.5
2 2.5 Anaesthetic
Patient 2
3 depth
3.5
4
4.5
5
.5
I
1.5
. . . 2 2.5 Anaesthetic
. 3 depth
3.5
q
.
.
4
4.5
5
4.5
5
30 Ten patients
y = -2.86x + 13.39, r2 = .61
y = -2.57x + 12.24, r2 = .66
25
2o
O
2O
O
0
o
oo
~ o o
~ . ~ ° o
o
5
®
o
° °/
°'u o o
o
o
o - - - ~ t q j ~ S o , , d t s o ¢oo o - 0 o ~ 8 o~8o o
o
o~'
_
°
# oo
oo
,~
o oo
o
o
oo
-~
o
oo
a
o
o
o
o
o
o
°
o
o
o
o
ooa
10 o
o o
oo
8
OoOO~
• o
°o
o
o
°
-
oo
.
o 0
.5
I
1.5
2 2.5 Anaesthetic
3 depth
3.5
4
4.5
0
" ~ .5
I
"0"~~
1.5 2 2.5 3 Anaesthetic depth
~ 3.5
° 4
137
4c~
Pa~ent2
2
N
N
2,
2
3
4
5
6
7 8
9 101112
Median f r e q u e n c y [Hz]
2
3 4
5
6
7 8
9 1 0 1 1 12
Median frequency [Hz]
Fig. 9. Histogram illustrating the scatter in median frequency in the range of the anaesthetic scale between 2.5 and 3.5 for one patient (left) and for ten patients (right) anaesthetized with isoflurane. This is the region in the class probability histogram dominated by yellow indicating that the patient is under normal surgical anaesthesia.
suits obtained by the multiparametric method, we constructed a scale from 0 - 5 for the level of anaesthesia and related it to calculations of single parameters (Figs. 8 and 9). The parameters were significantly correlated to anaesthetic depth (Fig. 8). The scatter around the regression line was however great indicating that the performance of single parameters is poor compared with the multiparametric descriptions of spectral information. The absolute values for median frequencies obtained with our technique are higher than those determined by Schwilden and Stoechel [4,17]. They used an F F T algorithm on a non-filtered signal. In
our method the signal was preemphasized to counteract the filtering effect when the E E G activity is transmitted through the skull. With the multiparametric method described in this paper it was shown that the three drugs modified the EEG in such a way that specific amplitude/frequency classes were found at the same level of anaesthesia. The learning sets for the three drugs were different, but they are similar for all patients anaesthetized with the same drug. In all patients a change to lighter or deeper levels of anaesthesia caused a stereotyped shift in class histogram to profiles representing the same
Fig. 8. Scatter plot of the median frequency (MF), the spectral edge frequency (SEF) and the peak frequency (PF) as a function of depth of anaesthesia (D) for one patient (left) and for ten patients (right) anaesthetized with isoflurane. The regression lines are MF (Hz) = - 1.5 (D) + 12.7, r = 0.5 (one patient), M F (Hz) = - 2.1 (D) + 13.8, r = 0.7 (ten patients), SEF (Hz) = - 1.3 (D) + 24.4, r = 0.2 (one patient), SEF (Hz) = - 1.9 (D) + 26.3, r = 0.4 (ten patients), PF (Hz) = - 2.9 (D) + 13.4, r = 0.6 (one patient), PF (Hz) = - 2.6 ( D ) + 12.2, r = 0.7 (ten patients).
138
trend. It was an important finding that the c o m mon learning sets based on five patients could be used for patients outside the learning material. This demonstrates that general anaesthesia tends to decrease the interpatient variability of the EEG trace. The improvement in the performance of the method, when only the five patients were included in the clustering analysis, indicated that the method as many other parametric statistical methods, e.g. linear regression, was sensitive to outhers. Aberrant classes interfered with the clustering process because features from these classes deviated heavily from the mean values and decreased the performance of the clustering process. This result emphasizes that it is essential to remove outliers, when a reference material is collected [20]. The improvement in performance was most likely due to the exclusion of periods with artifacts in the signal. The study showed that clustering analysis is useful for investigation of changes in the electroencephalogram caused by specific alterations of physiological conditions. The administration of general anaesthetics was used as an example for development of a universal anaesthetic depth monitor. The decrease in cost and increase in power of personal computers now allow for a clinical test of the method.
Acknowledgement This study was supported by the Danish Council for Technical and Scientific Research, by the Danish Medical Research Council, and by KUSIN.
References [1] P.F. Prior and D. Maynard, Monitoring Cerebral Function. Long Term Monitoring of EEG and Evoked Potentials, p. 441 (Elsevier, Amsterdam, 1987). [2] F.H. Lopes da Silva, Computerized EEG analysis: a tutorial overview, in: A Textbook of Clinical Neurophysiology, eds. A.M. Halliday, S.R. Butler and R. Paul, pp. 61-102 (John Wiley and Sons, New York, 1987). [3] M. Rosen and J.N. Lunn, Consciousness Awareness and
Pain in General Anaesthesia, p. 195 (Butterworth, Essex, 1987). [4] H. Schwilden and H. Stoeckel, Quantitative EEG analysis during anaesthesia with isoflurane in nitrous oxide at 1.3 and 1.5 MAC, Br. J. Anaesth. 59 (1987) 738-745. [5] H. Schwilden, J. Schtittler and H. Stoeckel, Closed loop feed back control of methohexital anaesthesia by quantitative EEG analysis in humans, Anesthesiology 67 (1987) 341-347. [6] I.J. Rampil and R.S. Matteo, Changes in EEG spectral edge frequency correlate with haemodynamic response in laryngoscopy and intubation, Anesthesiology 67 (1987) 139-42. [7] C.E. Thomsen, K.N. Christensen and A. Rosenfalck, Computerized monitoring of depth of anaesthesia with isoflurane. Br. J. Anaesth. 63 (1989) 36-43. [8] B.H. Jansen, EEG Segmentation and Classification, Thesis (Vrije Universiteit, Amsterdam, 1979). [9] J. Makhoul, Linear prediction: a tutorial review, Proc. IEEE 63 (1975) 561-580. [10] J.D. Markel and A.H. Grey, Linear Prediction of Speech (Springer-Verlag, New York, 1976). [11] R.O. Duda and P.E. Hart, Pattern Classification and Scene Analysis (John Wiley and Sons, New York, 1973). [12] M.R. Anderberg, Cluster Analysis for Apphcations (Academic Press, New York, 1973). [13] H. Sp~th, Cluster Analysis Algorithms (John Wiley and Sons, New York, 1980). [14] R.G. Bickford, J. Brimm, L. Berger and M. Aung, Application of compressed spectral array in clinical EEG, in: Automation of Chnical Electroencephalography, eds. P. KeUaway and I. Petersen, pp. 55-64 (Raven Press, New York, 1973). [15] J.R. Cox, F.M. Nolle and R.M. Arthur, Digital analysis of the electroencephalogram, the blood pressure wave, and the electrocardiogram, Proc. IEEE 60 (1972) 1137-1164. [16] A. Kumar, A real-time system for pattern recognition of human sleep stages by fuzzy system analysis, Pattern Recognit. 9 (1977) 43-46. [17] H. Stoeckel and H. Schwilden, Median EEG frequency, in: Consciousness Awareness and Pain in General Anaesthesia, eds. M. Rosen and J.N. Lunn, pp. 53-60 (Butterworth, Essex, 1987). [18] F. White and W.A. Boyle, Relationship between hemodynamic and electroencephalographic changes during general anesthesia, Anesth. Analg. 68 (1989) 177-181. [19] P.M. Yate, D.E. Maynard, E. Major, M. Frank, A.J.W. Verniquet, H.K. Adams and E. Douglas, Anaesthesia with ICI 35 868 monitored by the cerebral function analysing monitor (CFAM), Eur. J. Anaesthesiol. 3 (1987) 159-166. [20] B. Falck, S. Andreassen, T. Groth, H. Lang, M. Melander, A. Nurmi, A. Puusa, A. Rosenfalck, E. Sthlberg and M. Suojanen, The development of a multicenter database for reference values in chnical neurophysiology - - principles and examples, Comput. Methods Programs Biomed. 34 (1991) 145-162.
125
COMMET 01137
Assessment of anaesthetic depth by clustering analysis and autoregressive modelling of electroencephalograms C . E . T h o m s e n a, A. R o s e n f a l c k 1 a n d K . N o r r e g a a r d C h r i s t e n s e n 2 1 Department of Medical Informatics and Image Analysis, Aalborg University, DK 9220 Aalborg, Denmark, and 2 Department of Anaesthesia, Aalborg Hospital, DK 9000 Aalborg, Denmark
The brain activity electroencephalogram(EEG) was recorded from 30 healthy women scheduled for hysterectomy. The patients were anaesthetized with isoflurane, halothane or etomidate/fentanyl. A multiparametric method was used for extraction of amplitude and frequency information from the EEG. The method applied autoregressive modelling of the signal, segmented in 2 s fixed intervals. The features from the EEG segments were used for learning and for classification. The learning process was unsupervised and hierarchical clustering analysis was used to construct a learning set of EEG amplitude-frequency patterns for each of the three anaesthetic drugs. These EEG patterns were assigned to a colour code corresponding to similar clinical states. A common learning set could be used for all patients anaesthetized with the same drug. The classification process could be performed on-line and the results were displayed in a class probability histogram. This histogram reflected in all patients the depth of anaesthesia, when the concentration of the anaesthetic agent was adjusted either based on clinical signs or according to the protocol. This uniform display, where colours in a class probability histogram indicate the depth of anaesthesia, may in the future serve as on-line advice for the administration of anaesthetics. A comparison of multiparametric with single parametric methods, based on calculation of median, spectral edge and peak frequencies, questions the reliability of the single parametric methods in monitoring anaesthetic depth. Anesthesia, depth; Unsupervised learning; Hierarchical clustering; Isoflurane; Etomidate; Halothane; Human interface; On-line computer methods; Electroencephalogram amplitude-frequency patterns
1. Introduction C o m p u t e r i z e d analysis of the electroencephalogram ( E E G ) proved to b e valuable for long time m o n i t o r i n g of patients with head i n j u r y a n d during surgery [1]. T h e r e are two m a j o r reasons for using o n - l i n e computerized methods: (i) the imm e n s e a m o u n t of d a t a can be reduced in a n objective way, a n d (ii) signal analysis m e t h o d s c a n extract i n f o r m a t i o n on the a m p l i t u d e - f r e q u e n c y
Correspondence: C.E. Thomsen, Department of Medical Informatics and Image Analysis,Aalborg University, Frederik Bajers Vej 7 D, DK 9220 Aalborg, Denmark.
p a t t e r n of the E E G , which c a n n o t be detected b y visual i n s p e c t i o n of the E E G trace [2]. The use of E E G m o n i t o r i n g as a guide for the d e p t h of anaesthesia has b e e n controversial. D u r ing surgery m a n y factors m a y i n f l u e n c e the E E G a n d it becomes difficult to detect changes that can solely be related to the d e p t h of anaesthesia, because fluctuations in the arterial c a r b o n dioxide tension, b o d y temperature, p a i n stimuli a n d cardiovascular changes are present. I n emergency cases hypoxia will interfere with changes caused b y anaesthetics [1]. C u r a r i z a t i o n a n d the a d v e n t of new anaesthetics leave the anaesthetist with limited i n d i c a t i o n s of the anaesthetic depth. T o o light anaesthesia
126
might lead to increased need for analgesic and sedative drugs in the post-surgical period and to increased post-surgical stress in general. Awareness of pain or recall of other events during general anaesthesia has become a current medico-legal problem in the U.S.A. and in the U.K. [3]. The cerebral function analysing monitor (CFAM) has proven to be an appropriate EEG monitor in clinical practice [1]. The EEG activity is displayed in the traditional frequency bands (delta, theta, alpha and beta) together with the mean peak to peak value and the 10th and 90th percentiles of amplitude. The use of EEG monitoring systems still requires an extensive experience in interpretation of changes in amplitude and frequency content of the EEG. Efforts have been made to extract single parameters of the EEG activity to optimize the administration of anaesthetics. The median frequency (MF) of the power spectrum has been used as a single parameter characterizing the degree of central nervous system depression [4,5]. Another parameter is the spectral edge frequency (SEF) introduced by Rampil [6]. The development of the present method aimed towards the application of multiparametric techniques, to extract specific classes of amplitudefrequency patterns related to different levels of anaesthesia [7]. By means of this analysis, uniform information concerning the depth of anaesthesia could be obtained for three different anaesthetic techniques (halothane, isoflurane and etomidate/ fentanyl) despite the agent specificity of electroencephalographic changes.
2. Patients We investigated 30 healthy women (29-55 years old, body weight not exceeding 10% of ideal) undergoing simple hysterectomy. There was no history of neurological disorders or administration of analgesics and tranquillizers in the weeks before operation. The day before operation, informed consent was obtained. The investigation was approved by the local Ethics Committee. Three groups with ten patients in each group were investigated:
Group A was anaesthetized with halothane. Induction was made with thiopentone. For tracheal intubation a priming dose of gallamine and suxamethonium was administered. Anaesthesia was maintained with halothane 0.75-1.5 vol% in oxygen-nitrous oxide 1 : 2. In group B induction of anaesthesia was performed as described for group A, and anaesthesia was maintained with isoflurane 1.0-2.5 vol% in oxygen-nitrous oxide 1 : 2. In group C anaesthesia was induced and maintained with etomidate and fentanyl. Maintenance infusion was 0.5-1.0 m g / m i n etomidate and 1.63.5 m g / m i n fentanyl. All patients were mechanically ventilated by MC 801 Dameca anaesthetic ventilator. The duration of anaesthesia was 70-175 min. Monitoring of EEG from both hemispheres (P3-Fpa and P4-Fp2) started before induction of anaesthesia. In all patients temperature, end-tidal carbon dioxide concentration and mean arterial blood pressure were kept within physiological limits. All information was stored on a Racal Store 8 frequency modulated tape recorder. Further analysis of signals was performed off-line. From induction of anaesthesia the mechanical ventilator was adjusted to an end-tidal carbon dioxide concentration of 5.0-5.5 vol% corresponding to normoventilation. 20-25 rain after start of surgery, when the clinical condition was considered stable, a slight hyperventilation was induced corresponding to an end-tidal carbon dioxide concentration of 4.0-4.5 vol%. After 20 min hyperventilation was discontinued. Approximately 30 min after normoventilation was reestablished, the dosage of anaesthetics was increased by approximately 50% to study the E E G at a deeper anaesthetic level. After 20 min the maintenance concentration was reestablished.
3. Signal processing The signals were analysed off-line using an Intel 310 computer system equipped with a high resolution colour graphic display (1024 x 512 pixels in 256 colours) and analog to digital converter (12 bit resolution).
127 The EEG was filtered with a 25 Hz 4th order antialiasing filter, sampled at 100 Hz and preemphasized digitally with a first order high-pass filter at 4.2 Hz (Fig. 1). This filter gives approximately 10 d B / d e c a d e preemphasizing up to 25 Hz. The digitized EEG signal was divided into 2 s fixed segments using rectangular window function. The length of the window must be greater than 1 / 4 of the wavelength of the lowest frequency, and short enough to ensure stationarity within the window. A window length of 1-5 s is suitable, but 2 s preferable [8].
3.1. Feature extraction From each segment 11 features were extracted (ten normalized autocorrelation coefficients (ACCs) and RMS amplitude). These features were transformed to amplitude or power spectra by autoregressive modelling using the Durbin algorithm [9]. ACCs: N-i--1
r(i)=
~
x(n)x(n+i)
f o r i = 0 . . . . . 10
n=0
10o
o
y
g~
/ ,/
I02L 10 I
/
Normalized ACCs:
R(i)=r(i)/r(O)
for i = 1', . . . . 10
Amplitude: RMS = c r ( ( ~
c = scaling constant
During experiments ACCs were compared to reflection coefficients [9,10] and prediction coefficients [9] combined with either the mean squared power or the RMS amplitude. Prediction coefficients, reflection coefficients and ACCs contain exactly the same information, but encoded quite differently. The normalized ACCs combined with the RMS amplitude were chosen as features because they proved to be the least sensitive set of features with respect to interpatient variability.
4. Learning/clustering analysis Based on feature vectors extracted from a representative set of segments the learning was performed unsupervised by hierarchical clustering. The aim of the unsupervised approach was to find a natural grouping in the EEG pattern (not hmited by the definition of delta, theta, alpha and beta activities). But how do we define 'natural' and how can we express that data units within a cluster are more similar than data units taken from two different clusters? These questions can be expressed as two separate problems: - how is similarity between data units measured? and - how can the partition of data units in clusters be evaluated?
4.1. Similarity measures 100
101
10z
Frequency lHz] Fig. 1. Frequency response of the high-pass filter used to preemphasize the EEG signal. In the range from 2-20 Hz the filter characteristic is similar to that used by Prior and Maynard [1] in the CFAM. In the range below 2 Hz the filter has a slope of 20 dB/decade, whereas Prior and Maynard use filtering with 60 riB/decade.
An obvious choice of similarity (or dissimilarity) measure between data units is the Euclidean or geometric distance (1). A generalization of the Euclidean distance is the Mahalanobis distance (2). The Mahalanobis distance includes the covariance matrix and therefore compensates for the correlation between components of the feature vector.
128 The optimal clustering is the one which minimizes L [ll].
Euchdean distance:
de = ~ ( x - _ ~ ) r ( x - ~ ) =/l_x-~ll
(1) 4.3. Hierarchical clustering
Mahalanobis distance:
dm = ~/(x - ~) T~- l(X --_/£)
(2)
4.2. Criteria for clustering An essential problem is to evaluate the quality of the clustering analysis, and to define the optimal division of the n design units into an unknown number c of groups. An often used criterion function is the error sum of squares (or the minimum variance method). Let nj be the number of data units and x # the k-th sample both for the j-th cluster, and let ~j be the mean value for that cluster. Then the error sum of squares is defined using the Euclidean distance:
nj Je = E E (X--jk--I£j)T(x--jk--~j) j~l k~l c
(3/
A matrix of similarities between all feature vectors is the usual base for hierarchical clustering analysis. However this matrix contains n * ( n + 1 ) / 2 elements, where n denotes the total number of data units. A considerable reduction in m e m o r y could be obtained when, for each input vector, only the identifier and the distance to the nearest vector is stored. This matrix only contains 2n elements. The result from hierarchical clustering can be represented in a dendrogram (Fig. 2). The branches in the bottom represent one entity or one data unit, while the root at the top represents the entire collection of entities. The attention is here directed to agglomerative hierarchical clustering building the tree from branches to the root. Notice that in hierarchical methods clusters are nested. When two clusters are merged, they are joined permanently and become a building block for later merges. Herein lie both the strength and the weakness of hierarchical clustering analysis.
0 m
I 40 .~
60~
i00 Fig. 2. A dendrogram is the graphic description of the results from hierarchical clustering analysis. The Figure illustrates how clusters are merged as a function of similarity (ordinate). The spectra shown are mean values of the spectra for merging branches in the tree
from the subtrees below. This dendrogram is only part of the dendrogram for the group of patients anaesthetized with isoflurane.
129
The general procedure for agglomerative hierarchical methods can be described in Pseudo Pascal:
a)
PROCEDURE HIERARCHICALCLUSTERING BEGIN FOR i:=l TO n_unitDO c := n_umi; WHILE c > 1 DO BEGIN Find the two clusters,which have a minimum inter-distance.Let's call them •i and t)j; Merge f/i and Dj in a new cluster Di; Delete Gj; Update similaritymatrix for D~ Update similaritymatrix for clusterswith Di or ~j as nearestneighbours c := c-1; END; END;
b)
The number of data units is denoted n_uni t and ~2, is the i-th cluster. In the procedure the term 'distance' is used. In the previous section we discussed the choice of similarity measure. Between which points within the two clusters are we measuring, and shall the distance be weighted according to the size of the clusters? Three often used cluster distances are:
dmin(~'~i,~Qj)=
min
ll_x-_x'll
dma×(#,, ~j)=
max
I I x - x'll
dmean(~i,
~))=
Fig. 3. T o illustrate t h r e e d i f f e r e n t s i m i l a r i t y m e a s u r e s f o r c l u s t e r i n g of E E G p a r a m e t e r s f o r the c o n s t r u c t i o n o f a c o m m o n l e a r n i n g set f o r a specific a n a e s t h e t i c d r u g . E l e m e n t s f r o m c l u s t e r Y~l a n d [22 a r e m a r k e d w i t h (o) a n d (x). T h e m e a n values f o r the t w o clusters a r e m a r k e d w i t h ( O ) a n d (X). I n ( a ) , ( b ) a n d ( c ) the m i n i m u m , t h e m a x i m u m a n d the m e a n o c c u r r i n g d i s t a n c e b e t w e e n e l e m e n t s f r o m the c l u s t e r s ~ a n d 122 are u s e d as i n t e r - c l u s t e r d i s t a n c e .
II_bti--~j]l
dmin(#i,
~j) is the minimum occurring distance from data units in cluster 12~ to data units in cluster 12j (Fig. 3a). Analogous dmax (12i, 9j) is the maximum occurring distance between samples within the two clusters (Fig. 3b). In the third method the distance is measured between mean values or centroids of clusters (Fig. 3c). Combining the minimum variance methods and the centroid methods, it was possible to obtain stepwise optimization of the error sum of squares Je (3) by using a weighted dmean(Oi, ~2j) [12,13]: d~ean(~i,
c)
~j)= ~/(rl i * rlj)/(rl i + tl2)
*
[[_~i --_PL, H
(4) n i and nj are the number of data units within cluster 12i and f2j. This measure could be
interpreted as the square root of the increase in the error sum of squares if cluster/2, and Dj were merged. Even though the clustering was stepwise optimal, there was no quarantee for global minimization of this criteria function.
5. Construction of learning sets
During learning we selected a set of 504 E E G segments from each patient representing typical EEG patterns from different levels of anaesthesia. During hierarchical clustering segments were merged together - - step by step starting with the 504 segments ending with one. The modified Euclidean distance (4) between feature vectors was used as similarity measure. The clustering sequence was represented graphically in a dendro-
130
of EEG patterns in all patients anaesthetized with that particular agent.
6. Classification
Fig. 4. To illustrate 'linked' or stepwise clustering analysis, when data from five or ten patients were included in the learning process. Below: Clustering sequences from each patient, starting with 504 feature vectors. The arrow indicates the level where 50 classes from each dendrogram are selected and merged in a dataset representing c o m m o n E E G patterns. Above: The c o m m o n set of patterns was used in an additional clustering analysis to define a c o m m o n set of classes (to be used as a learning set during classification), specific for one anaesthetic agent.
gram. The dendrogram from hierarchical clustering could be inspected interactively. This inspection was guided by displaying mean spectra from subtrees below selected nodes (Fig. 2). During inspection special attention was paid to nodes where the distance increased (i.e. similarity decreased) rapidly. Using this method a set of classes typical for this patient could be defined. Classes containing equal or less segments than the number of features were removed from the class definition. To generate a common learning set for a specific anaesthetic agent, data from up to ten patients were included in the clustering sequence using 'linked' clustering analysis (Fig. 4). From the dendrograms, for each of the patients anaesthetized with the same drug, we selected a level where 50 clusters were represented. The 50 clusters from each patient were merged and clustering analysis repeated. After this procedure 12-16 classes were defined to be used as learning set for classification
The probability p(~2~ I_x) for the object belonging to the class ~2~ given the input feature vector _x was calculated from Bayes theorem. Calculating p (12~ I_x) for all classes and choosing the class with the largest probability gave a minimum probability for misclassifications [11]. We assumed that the probability density function for each class had a multidimensional normal distribution and that all classes had the same a priori probability. The individual class centroids were determined by the mean value vectors it, and the size and shape by the pooled covariance matrix Z. The decision problem was expressed in a set of linear discriminant functions [11] g~(x), i = 1 . . . . . c. The classifier assigned the feature vector to the class 12~, if g , ( x ) > g j ( x ) for all j v~ i.
The set of linear discriminant functions used in the system were g~ (x_) = Wlr~x + WO~ --1
where: W 1 i = Y, #, and: W0, = - ( 1 / 2 ) ~ Z --
-
T
-I
~,
The results from the classification were presented as 'raw' class probability histograms displayed together with density spectral array (DSA) plots as a function of time (Fig. 5). DSA plots are the classical presentation of the amplitude frequency pattern of the EEG [14]. Each line represents one amplitude spectrum, where the
Fig. 5. Example of probability scores for EEG amplitude frequency classes derived from a patient anaesthetized with isoflurane. Vertical axis is time and the total duration is 2 h and 50 min. Events are marked along the time axis at the left: (T) thiopentone; (x), (X) start and stop of anaesthesia; (o) onset of surgery; ( 1"), ( - ) start and stop of hyperventilation; ( > ), ( < ), increased and decreased induction of anaesthesia. Note that for this patient hyperventilation was induced before the change in induction of anaesthesia.
0~ 4-J
DSA spectrum
Class 13 histogram
Class 11+12+13 histogram
Class 4 histogram
Class 2+3+4+8 histogram
Class 6 histogram
Class 6+7+10 histogram
W-
1 -I ..m
I
--! P
i
I :
q
4 4
"t
P
iii~Li,ill i---!!
b X
ip
132 anaesthesia (e.g. s t a r t / s t o p of isoflurane induction and i n / d e c r e a s e d concentration of inspired isoflurane. For the example in Fig. 5 it was clear that class 6, 7 and 10 corresponded to increased anaeshetic depth, and that class 11, 12 and 13 occurred during induction of anaesthesia and in the awakening period, corresponding to light anaesthesia. Using this procedure classes or sums of classes were assigned colours according to the level of anaesthesia (Fig. 6 [7]) for each anaesthetic agent
frequency axis is the abscissa and the amplitude is expressed in a grey-scale, black denotes maximal amplitude. This presentation is used in Figs. 5 and 7 to support the evaluation of the class probability histogram. The DSA traces are globally scaled to give a true impression of amplitude changes during surgery. The class probability histograms for some of the classes and for summated classes are shown in black and white print. The occurrence of each class was then related to clinical statements, especially statements related to the depth of
-- "-----i~.:I--F:ISS 1 6
DROWSY (dark blue) ,
.
'
,
,
n
if: R S S
i2
VERY LIGHT (light b l u e )
N O R M A L
0
20
30
50
0
20
30
50
°I~"" C-L A S S ~
2
0
20
30
50
30
50
~ o b" -C L'A ~ S S ,~, 3, ,
"
--
-- ~
-
CLASS
8
[yellow}
6
' 26
3b
'
5b
o
2o
3o
50
30
5o
o
20
~o
26 3~
I 50
[purple}
0 ~0 30
----Pre.
.
' ~
3b
'
5b
o
.
.
2o
.
.
Fig. 6. Set of averaged amplitude spectra for classes from the group of patients anaesthetized with isoflurane. One to four classes were assigned to the same colour code because they all appear at the same depth of anaesthesia.
133 (blue = 'drowsy', light blue = 'very light anaesthesia', green ='light anaesthesia', yellow = ' n o r m a l surgical anaesthetic depth', purple = ' d e e p anaesthesia' and red = ' very deep anaesthesia'). The sets of classes and colour codings were used during on-line pattern recognition. Classes representing awake states, artifacts and movements were assigned to background colour (black on the screen and white on the paper print out). Even though the whole configuration of the amplitude spectrum was essential in classification of anaesthetic depth, the peak frequency was found to be one of the dominating single parameters in the learning set. For isoflurane a peak frequency of 10-12 Hz was found during light anaesthesia. Surgical depth of anaesthesia was characterized by a peak frequency of 4 - 6 Hz and deep levels of anaesthesia by a peak frequency of 1-3 Hz.
7. Results from monitoring of anaesthetic depth
The three learning sets of classes and colour codings, obtained for each of the three different anaesthetic drugs, were used during on-line classification. The information from the EEG analysis was displayed as a function of time [7]: (a) as density spectral array plots (DSA) (Fig. 7, left column). Each p l o t / l i n e represents the averaged amplitude spectrum during the last 10 s; (b) as class probability histograms (Fig. 7, right column). Each line represents the relative occurrence of any class in percent during the last 10 s. The class probability histograms were smoothed by a 2nd order low pass filter with two real poles. In all patients this histogram reflected the depth of anaesthesia when the concentration of the anaesthetic agent was adjusted either based on clinical signs or according to the protocol. During steady-state the class probability histogram often indicated an increased anaesthetic level, when pain stimuli were less, towards the end of surgery.
8. Results when outliers were excluded from the learning set
The disturbing effect of outliers in the learning sets was clearly demonstrated by clustering analysis from the three groups of patients. During isoflurane anaesthesia, all recordings were very 'clean' and contained practically no artifacts, resulting in a reliable set of classes, and a well defined anaesthetic scale. For the two other anaesthetic agents the sets of classes based on data from all ten patients gave some inconsistency in the class probability histogram, and a number of classes could not be assigned according to the anaesthetic depth. This problem was minimized by reducing the number of patients contributing to the learning set to five. The five patients were selected as those who responded most uniformly to surgical depth of anaesthesia. When these learning sets were used for all ten patients in each group, there was a better consistency and a more reliable configuration of the class probability histogram during online analysis. The improvement was not only related to patients inside the learning group, but was also observed for the five patients not included in the learning material. 8.1. Single parametric methods
Assuming that the anaesthetic scale used in the class probability histogram (Fig. 7, right column) was reliable we evaluated a number of single parametric methods against the present multiparametric method. The anaesthetic scale was converted to an analog scale (0-5) giving each anaesthetic depth from Fig. 6 a value: 'drowsy' = 0, ' v e r y light anaesthesia' = 1, 'light anaesthesia' = 2, 'normal surgical anaesthetic d e p t h ' = 3, 'deep anaesthesia' = 4, and ' very deep anaesthesia' = 5. From each 10 s epoch of the EEG, the median frequency (MF), the peak frequency (PF) and the spectral edge frequency (SEF) were calculated. The values of MF, SEF and PF were calculated for one patient solely and for the ten patients in group B together. The values for approximately 700 epochs in each patient were plotted against
30~U-% i v
134
" ~ , ' ~ - - .~¢v-'~ -,.',,* w' r l p ~ m " l F w
0
DSA spectrum 15 30[Hz]
0
Class prob. histogram 50 1001%1
T
t
<:
X
CLAS:~
CI3D
I I'JL~ :
135
the analogue scale for depth of anaesthesia, and the linear regression lines were calculated (Fig. 8). MF, SEF and PF frequencies decreased as expected with the depth of anaesthesia (Fig. 9). In evaluating the performance of these single EEG parameters one should keep in mind that the misclassification rate corresponded to the squared correlation coefficients. The scatter around the regression line for the median frequency was illustrated by selecting the values encountered in the range of the anaesthetic scale between 2.5 to 3.5. This corresponds to periods in the class probability histogram, where normal surgical anaesthetic depth (yellow) was scored by the multiparametric method. The median frequencies for one patient ranged from 5 to 10 Hz and for ten patients from 4 to 11 Hz (Fig. 9).
9. Discussion
If EEG monitoring shall serve as decision support during administration of anaesthetic drugs, a 'normal' population of patients has to be investigated during standard conditions. The aim of this study was to use a multiparametric method to analyse the EEG from 'normal' patients anaesthetized with different anaesthetic techniques and to facilitate the implementation of this method in an 'anaesthetic level monitoring system'. The multiparametric method used features from autoregressive modelling of the EEG signal, segmented in fixed time intervals. Before the feature extraction the signal was preemphasized or prewhitened by a filter (Fig. 1) similar to that used by Prior and Maynard [1]. Preemphasizing was done to reduce the redundant information in the signal caused by the inverse proportionality between the amplitude and frequency of EEG changes during anaesthesia [15,16].
The learning process was unsupervised and applied clustering analysis to define a learning-set of classes for three anaesthetic drugs. Each set of basic classes were assigned colours according to the level of anaesthesia. The l e a r n i n g s e t w a s c o m m o n for each anaesthetic agent but different for the three agents. It was used for classification of the EEG activity. The result of the classification was presented in class probability histograms and demonstrated a high degree of simplicity. To suppress short periods of disturbances from artifacts in the EEG trace the class probability histograms were smoothed by a 2nd order low pass filter with two real poles. The effect was that the histograms became more clear. At different levels of anaesthesia more than two classes were seldom represented. The simple display, where the same colour is related to the same level of anaesthesia for three different anaesthetic agents may support the anaesthesiologist. It has been suggested that single parametric methods: the median and the spectral edge frequencies can be used to control of the depth of anaesthesia. Stoeckel and Schwilden [4,17] have shown that the time course of the median frequency parallels the visual staging of anaesthesia and is closely related to the time course of the declining blood concentration of etomidate. They also found changes in the median frequency when the patient was anaesthetized with halothane and isoflurane but the variation with the depth of anaesthesia was different for the three drugs. Rampil and Matteo [6] found consistent changes in SEF related to haemodynamic responses during anaesthesia. White and Boyle [18] were not able to confirm these findings but found good correlation between propofol blood levels and EEG changes recorded by the CFAM [19]. A comparison of multiparametric and single parametric methods was attempted. From the re-
Fig. 7. Example of a colour copy of results from EEG analysis from a patient anaesthetized with isoflurane. Vertical axis is time and the total duration is 2 h and 20 rain. Events are marked along the time axis at the left; (T) thiopentone; (x), (X) start and stop of anaesthesia; (o) onset of surgery; ( > ) , ( < ) , increased and decreased induction of anaesthesia; (1"), ( - ) start and stop of hyperventilation. Left column is the DSA plot. Right column is the class probability histogram, colour coded as illustrated in Fig. 6. The trace above is a 2 s segment of the raw EEG signal, which is continuously displayed for control.
136
30
30
Ten paUents
Pa0ent 2 y : -I.5; ÷ 12.74, r2 =
y = -2,11x + 13.79, r2 = .88
.52
25 ~'o
2~-
~_.
20
~
15
o
o
o
o o
°
oo
o o ooo
o
o
o co ° o O
o o ¢
o ~~
o
~
o
o o
o
~
o qo~ooo
o
a
o
o
~
tO
oo o
o o
o
5
•5
l
1.5
2 2.5 3 Anaesthetic depth
y = -1.28x
Patient 2
3.5
+ 24.42,
r2
4
4.5
0
45
= .16
.5
oo
o
I
1.5
2 2.5 3 Anaesthetiedepth
"Ten pa~ents
4
+ 26.33,
r2
4.5
5
o y = -1.89x
o
= .38
o o
O
40
3.5
40
o
o
o
o
~o
o
O
o
o
35
~ o
°o o
o°
30
o
•~
oo o
o e
o
°
o
o
o
o #°°0
ff
o I
30
~
25
o
0
0
Oo
20
o 15} 0
I .5
oo
o
oo
2
o
~
°
:
o
o
25~,~._ o 8 ~ O o
o
o o%
°
~'
o
8
35
1.5
2 2.5 Anaesthetic
Patient 2
3 depth
3.5
4
4.5
5
.5
I
1.5
. . . 2 2.5 Anaesthetic
. 3 depth
3.5
q
.
.
4
4.5
5
4.5
5
30 Ten patients
y = -2.86x + 13.39, r2 = .61
y = -2.57x + 12.24, r2 = .66
25
2o
O
2O
O
0
o
oo
~ o o
~ . ~ ° o
o
5
®
o
° °/
°'u o o
o
o
o - - - ~ t q j ~ S o , , d t s o ¢oo o - 0 o ~ 8 o~8o o
o
o~'
_
°
# oo
oo
,~
o oo
o
o
oo
-~
o
oo
a
o
o
o
o
o
o
°
o
o
o
o
ooa
10 o
o o
oo
8
OoOO~
• o
°o
o
o
°
-
oo
.
o 0
.5
I
1.5
2 2.5 Anaesthetic
3 depth
3.5
4
4.5
0
" ~ .5
I
"0"~~
1.5 2 2.5 3 Anaesthetic depth
~ 3.5
° 4
137
4c~
Pa~ent2
2
N
N
2,
2
3
4
5
6
7 8
9 101112
Median f r e q u e n c y [Hz]
2
3 4
5
6
7 8
9 1 0 1 1 12
Median frequency [Hz]
Fig. 9. Histogram illustrating the scatter in median frequency in the range of the anaesthetic scale between 2.5 and 3.5 for one patient (left) and for ten patients (right) anaesthetized with isoflurane. This is the region in the class probability histogram dominated by yellow indicating that the patient is under normal surgical anaesthesia.
suits obtained by the multiparametric method, we constructed a scale from 0 - 5 for the level of anaesthesia and related it to calculations of single parameters (Figs. 8 and 9). The parameters were significantly correlated to anaesthetic depth (Fig. 8). The scatter around the regression line was however great indicating that the performance of single parameters is poor compared with the multiparametric descriptions of spectral information. The absolute values for median frequencies obtained with our technique are higher than those determined by Schwilden and Stoechel [4,17]. They used an F F T algorithm on a non-filtered signal. In
our method the signal was preemphasized to counteract the filtering effect when the E E G activity is transmitted through the skull. With the multiparametric method described in this paper it was shown that the three drugs modified the EEG in such a way that specific amplitude/frequency classes were found at the same level of anaesthesia. The learning sets for the three drugs were different, but they are similar for all patients anaesthetized with the same drug. In all patients a change to lighter or deeper levels of anaesthesia caused a stereotyped shift in class histogram to profiles representing the same
Fig. 8. Scatter plot of the median frequency (MF), the spectral edge frequency (SEF) and the peak frequency (PF) as a function of depth of anaesthesia (D) for one patient (left) and for ten patients (right) anaesthetized with isoflurane. The regression lines are MF (Hz) = - 1.5 (D) + 12.7, r = 0.5 (one patient), M F (Hz) = - 2.1 (D) + 13.8, r = 0.7 (ten patients), SEF (Hz) = - 1.3 (D) + 24.4, r = 0.2 (one patient), SEF (Hz) = - 1.9 (D) + 26.3, r = 0.4 (ten patients), PF (Hz) = - 2.9 (D) + 13.4, r = 0.6 (one patient), PF (Hz) = - 2.6 ( D ) + 12.2, r = 0.7 (ten patients).
138
trend. It was an important finding that the c o m mon learning sets based on five patients could be used for patients outside the learning material. This demonstrates that general anaesthesia tends to decrease the interpatient variability of the EEG trace. The improvement in the performance of the method, when only the five patients were included in the clustering analysis, indicated that the method as many other parametric statistical methods, e.g. linear regression, was sensitive to outhers. Aberrant classes interfered with the clustering process because features from these classes deviated heavily from the mean values and decreased the performance of the clustering process. This result emphasizes that it is essential to remove outliers, when a reference material is collected [20]. The improvement in performance was most likely due to the exclusion of periods with artifacts in the signal. The study showed that clustering analysis is useful for investigation of changes in the electroencephalogram caused by specific alterations of physiological conditions. The administration of general anaesthetics was used as an example for development of a universal anaesthetic depth monitor. The decrease in cost and increase in power of personal computers now allow for a clinical test of the method.
Acknowledgement This study was supported by the Danish Council for Technical and Scientific Research, by the Danish Medical Research Council, and by KUSIN.
References [1] P.F. Prior and D. Maynard, Monitoring Cerebral Function. Long Term Monitoring of EEG and Evoked Potentials, p. 441 (Elsevier, Amsterdam, 1987). [2] F.H. Lopes da Silva, Computerized EEG analysis: a tutorial overview, in: A Textbook of Clinical Neurophysiology, eds. A.M. Halliday, S.R. Butler and R. Paul, pp. 61-102 (John Wiley and Sons, New York, 1987). [3] M. Rosen and J.N. Lunn, Consciousness Awareness and
Pain in General Anaesthesia, p. 195 (Butterworth, Essex, 1987). [4] H. Schwilden and H. Stoeckel, Quantitative EEG analysis during anaesthesia with isoflurane in nitrous oxide at 1.3 and 1.5 MAC, Br. J. Anaesth. 59 (1987) 738-745. [5] H. Schwilden, J. Schtittler and H. Stoeckel, Closed loop feed back control of methohexital anaesthesia by quantitative EEG analysis in humans, Anesthesiology 67 (1987) 341-347. [6] I.J. Rampil and R.S. Matteo, Changes in EEG spectral edge frequency correlate with haemodynamic response in laryngoscopy and intubation, Anesthesiology 67 (1987) 139-42. [7] C.E. Thomsen, K.N. Christensen and A. Rosenfalck, Computerized monitoring of depth of anaesthesia with isoflurane. Br. J. Anaesth. 63 (1989) 36-43. [8] B.H. Jansen, EEG Segmentation and Classification, Thesis (Vrije Universiteit, Amsterdam, 1979). [9] J. Makhoul, Linear prediction: a tutorial review, Proc. IEEE 63 (1975) 561-580. [10] J.D. Markel and A.H. Grey, Linear Prediction of Speech (Springer-Verlag, New York, 1976). [11] R.O. Duda and P.E. Hart, Pattern Classification and Scene Analysis (John Wiley and Sons, New York, 1973). [12] M.R. Anderberg, Cluster Analysis for Apphcations (Academic Press, New York, 1973). [13] H. Sp~th, Cluster Analysis Algorithms (John Wiley and Sons, New York, 1980). [14] R.G. Bickford, J. Brimm, L. Berger and M. Aung, Application of compressed spectral array in clinical EEG, in: Automation of Chnical Electroencephalography, eds. P. KeUaway and I. Petersen, pp. 55-64 (Raven Press, New York, 1973). [15] J.R. Cox, F.M. Nolle and R.M. Arthur, Digital analysis of the electroencephalogram, the blood pressure wave, and the electrocardiogram, Proc. IEEE 60 (1972) 1137-1164. [16] A. Kumar, A real-time system for pattern recognition of human sleep stages by fuzzy system analysis, Pattern Recognit. 9 (1977) 43-46. [17] H. Stoeckel and H. Schwilden, Median EEG frequency, in: Consciousness Awareness and Pain in General Anaesthesia, eds. M. Rosen and J.N. Lunn, pp. 53-60 (Butterworth, Essex, 1987). [18] F. White and W.A. Boyle, Relationship between hemodynamic and electroencephalographic changes during general anesthesia, Anesth. Analg. 68 (1989) 177-181. [19] P.M. Yate, D.E. Maynard, E. Major, M. Frank, A.J.W. Verniquet, H.K. Adams and E. Douglas, Anaesthesia with ICI 35 868 monitored by the cerebral function analysing monitor (CFAM), Eur. J. Anaesthesiol. 3 (1987) 159-166. [20] B. Falck, S. Andreassen, T. Groth, H. Lang, M. Melander, A. Nurmi, A. Puusa, A. Rosenfalck, E. Sthlberg and M. Suojanen, The development of a multicenter database for reference values in chnical neurophysiology - - principles and examples, Comput. Methods Programs Biomed. 34 (1991) 145-162.