- Email: [email protected]
Application of simulated annealing for estimating BAEPs in endocochlear pathologies
Medical Engineering & Physics 24 (2002) 385–392 www.elsevier.com/locate/medengphy
Application of simulated annealing for estimating BAEPs in endocochlear pathologies A. Naı¨t-Ali ∗, P. Siarry Universite´ Paris 12, LERISS 61, avenue du Ge´ne´ral de Gaulle, 94010 Cre´teil, France Received 19 July 2001; received in revised form 8 March 2002; accepted 19 March 2002
Abstract In this paper, we propose a new approach aimed at handling the temporal Brainstem Auditory Evoked Potentials (BAEPs) nonstationarity. It is pointed out that for some endocochlear pathologies, BAEPs could be randomly delayed from one response to another. This non-stationarity leads to smoothed BAEPs when applying ensemble averaging or any other technique based on BAEPs stationarity. In that case, waves identification is very difficult, sometimes impossible. The problem consists in estimating time delays. Knowing the distribution of delays allows subsequent study of the dynamic of the cochlea and, perhaps, identification of the nature of its pathology. The approach suggested in this paper is based on Simulated Annealing, used to minimize a non-linear criterion involving delays. This technique is advantageously compared to the non-corrected ensemble averaging method, using a set of simulated data based on a realistic model. As an illustration, results based on real signals recorded from two patients are presented and discussed at the end of the paper. 2002 IPEM. Published by Elsevier Science Ltd. All rights reserved. Keywords: Evoked potentials; Auditory system; Simulated annealing; Endocochlear
1. Introduction Brainstem Auditory Evoked Potentials (BAEPs) are electrical signals generated in response to acoustical stimulations of the auditory system [1,2]. BAEPs are recorded on the patient’s vertex and mastoid. They are used particularly for early diagnosis of acoustic neuroma [3–5]. The electroencephalogram (EEG) corresponding to a natural physiological activity corrupts the response to a single stimulation. Three main reasons cause BAEPs extraction difficulties. First, the very low signal-to-noise ratio (SNR) (its value is sometimes around ⫺20 dB); second, problems in EEG non-stationarity; and third, BAEP non-stationarity. Direct observation or extraction of BAEPs is not usually possible. One of the most classic techniques used to estimate BAEPs is the ‘ensemble averaging’, that performs well, especially when EEG and BAEPs are stationary. Various other techniques have been developed by taking into
Corresponding author. Tel: +33-1-4517-14-90; fax: +33-1-45171492. E-mail address: [email protected] (A. Naı¨t-Ali). ∗
account the non-stationarity properties. A well-known method, called the LCA technique (‘Latency Corrected Average’), has been proposed by McGillem et al. [6]. Latencies are detected at each sweep; signals having the same latencies are gathered. Yu et al. [7] have improved this technique by working out the PC-LCA (‘Peak Component Latency Corrected Average’): sweeps are filtered using an adaptive filter, before detecting latencies. Results are, of course, strongly correlated with the SNR. Other techniques are proposed in the literature to process late evoked potentials, as an adaptive estimation of latency change using a DLMS (‘Direct Least Mean Square’) adaptive TDE (‘Time Delay Estimation’) algorithm. Here the EEG is considered as Gaussian [8,9]. Another variant of this method was recently proposed by Kong et al. [10]. It is called DLMP (‘Direct LeastMean p-norm’) adaptive TDE, that works well under αstable noise conditions. In this paper, we are interested only in BAEPs temporal non-stationarity. BAEPs may be randomly delayed from one stimulation to another. The averaged signal will then be smoothed, as shown in Fig. 1. The idea here consists in keeping the classical technique of averaging, which is no longer applied to signals as they are
1350-4533/02/$22.00 2002 IPEM. Published by Elsevier Science Ltd. All rights reserved. PII: S 1 3 5 0 - 4 5 3 3 ( 0 2 ) 0 0 0 3 4 - 6
386
A. Naı¨t-Ali, P. Siarry / Medical Engineering & Physics 24 (2002) 385–392
Fig. 1.
Effects of BAEPs time delay on the averaged signal.
recorded, but after an optimal correction. Our aim is to estimate the delays which enhance the averaged signal. The estimated delays could later be analysed from the angle of distribution to understand the dynamic of the cochlea. However, filtering EEG will not be considered; the single assumption made on the EEG is that its energy tends to zero for M given number of sweeps. The paper comprises four sections. In section 2, we present the mathematical basis of our technique, that uses Simulated Annealing to estimate the delays and to improve the averaged signal. In section 3, we apply the method to simulated BAEPs, perturbed in two different ways; the results are then discussed. Results based on real signals recorded from two patients are finally presented and discussed. In section 4 a conclusion is given relating to the proposed method.
2. Mathematical basis of BAEPs estimating 2.1. Time delayed BAEP model Let xi(n) be the sampled signal recorded in response to the ith stimulus. One can write: xi(n) ⫽ s(n ⫹ Di) ⫹ vi(n) n ⫽ 0,1,…,N⫺1
the BAEP delay, randomly distributed according a given distribution; vi(n)is the EEG signal; N is the number of samples in a given response. The averaged signal x¯(n)over M sweeps is easily calculated as follows:
冘
M⫺1
x¯(n) ⫽
where : n is a discrete time variable; s(n) is the BAEP having a stationary waveform; Di is an integer measuring
冘
M⫺1
(2)
If we consider that M sweeps are sufficient to make the residual noise energy negligible in comparison with BAEP averaged energy, one can write:
冘
M⫺1
x¯(n)⬇
1 s(n ⫹ Di) Mi ⫽ 0
(3)
The technique that we propose consists in estimating delays maximizing the energy of the averaged signal. An endocochlear pathology is supposed to be the single cause of the observed delays. To make a diagnosis, it is necessary to study and classify those delays. More explanations are reported in section 3.3. The criterion to be optimized is formulated as follows :
冘冉 冘
N⫺1
(1)
冘
M⫺1
1 1 1 xi(n) ⫽ s(n ⫹ Di) ⫹ v (n) Mi ⫽ 0 Mi ⫽ 0 Mi ⫽ 0 i
fD ⫽
n⫽0
M⫺1
1 x (n⫺Di) Mi ⫽ 0 i
冊
2
with D ⫽ [D0,D1,D2,...DM⫺1].
(4)
A. Naı¨t-Ali, P. Siarry / Medical Engineering & Physics 24 (2002) 385–392
The first response is considered as a reference, then D0 ⫽ 0. Actually, this hypothesis ensures the unicity of the solution. Maximizing fD leads to minimizing over D the following criterion:
冘冉 冘
N⫺1
JD ⫽ ⫺
n⫽0
M⫺1
1 x (n⫺Di) Mi ⫽ 0 i
冊
2
(5)
The criterion (5) is not quadratic or convex, therefore optimization is driven using Simulated Annealing algorithm (SA) [11], which was successfully applied to numerous multiminima optimization problems [12] [13]. 2.2. Simulated Annealing algorithm The problem of finding the minimum of a given function J of M variables has received much attention for many years. Several algorithms that give good results for functions having only one minimum have been proposed. However, those search methods do not work satisfactorily when the cost function comprises several (‘global’ or ‘local’) minima within the domain of interest, because they tend to stop at the first minimum encountered. In many research and development fields, engineers have to deal with minimization problems of ever increasing dimensionality, leading to a combinatorial-like explosion of the number of local minima. For handling such global optimization problems, new approaches, called ‘metaheuristics’, have emerged, most of them being inspired by biological or physical phenomena : in particular, Simulated Annealing, Genetic Algorithms, Tabu Search and Ant Colonies. We propose here the use of the Simulated Annealing method, for which efficiency has been demonstrated in the solution of difficult optimization problems, when only the objective function is available. Simulated Annealing aims at solving such problems ‘at the best’ : in theory, it avoids the classical trapping into objective function local minima. It was inspired by the ‘annealing’ experimental technique, currently used by metallurgists to obtain a well-ordered solid, corresponding to the minimal energy state. The most significant features of the method are conditionally guaranteed convergence towards a global optimum, but in practice rather sensitive tuning of parameters. A pseudo-code of Simulated Annealing is given in Table 1. Starting from an initial point Dini, SA generates a succession of points D1, D2, etc.... A new point D’, randomly generated from the current point D (step 3 of Table 1), is submitted for acceptance to the Metropolis criterion, defined at the bottom of Table 1. The Metropolis algorithm allows SA to estimate first (at high temperature) the ‘coarse grain’ behaviour of the function under test and then to examine (at low temperature) its ‘finer grain’ behaviour.
387
Table 1 Pseudo-code of SA algorithm and Metropolis rule 1
Randomly choose randomly an initial solution Diniand evaluate the criterion value JD; 2 Choose an initial temperature T; 3 Perturb the solution to have a new one D’ ⫽ D ⫹ ⌬D; 4 Calculate ⌬J ⫽ JD⬘⫺JD; 5 The solution is accepted or rejected according to the Metropolis rule (see below) ; 6 If the thermodynamic equilibrium is reached, then decrease the temperature, else go to step 3; 7 If the system is ‘frozen’, go to step else go to step 3; 8 Solution = the best point found. Metropolis rule If ⌬J ⫽ JD⬘⫺JDⱕ0 then keep this solution and D’; D←D’; JD⬘⫺JD Else calculate P ⫽ exp{⫺ }; T Generate a random number R uniformly distributed within [0,1]; If RⱕP, then accept the new solution D’; D←D’, else reject the solution D’.
The sequence of points generated by SA at a fixed temperature T is stopped as soon as the thermodynamic equilibrium, detected by some adequate condition, is reached (step 6 of Table 1). Then the temperature is suitably adjusted, before beginning a new sequence of moves, and so on. This process is finally stopped (step 7 of Table 1) when the system is ‘frozen’, i.e. at least one of several independent stopping tests on downhill moves number, temperature and number of evaluations of J is satisfied. These tests are basically equivalent in the sense that they all aim at terminating the annealing process through estimation of J local behaviour: the SA search is stopped if J is sufficiently flat in the neighbourhood under exploration, as detected by at least one of the tests. To sum up, SA allows ‘uphill’ moves under control of the T temperature parameter. It is clear that, at higher temperature, only the coarse grain behaviour of the objective function J is relevant to the search. As temperature decreases, the finer grain behaviour of the objective function J is examined, finally leading to a solution theoretically near to the global minimum of J. A pseudo-code of SA and Metropolis criterion is given in Table 1. Several criteria exist to stop SA (step 7 of Table 1); for example, one may decide to stop SA when no more tested solutions are accepted, or when the temperature T is lower than a threshold temperature close to 0. The parameters of SA are fixed as follows: 1. Initial temperature T0 : it can be evaluated using the following empirical technique: 앫 Do 100 random perturbations: evaluate the average aute energy variation ⬍ 兩⌬E兩 ⬎ 앫 Choose and initial rate of acceptance t0, according to the confidence in the ‘quality’ of the initial configuration, for example:
388
A. Naı¨t-Ali, P. Siarry / Medical Engineering & Physics 24 (2002) 385–392
– Bad quality: t0 ⫽ 50% (beginning at high temperature), – Good quality: t0 ⫽ 20% (beginning at low temperature), 앫 DeduceT0 from therelation exp(⫺ ⬍ 兩⌬E兩 ⬎ / T0) ⫽ t0 2. Temperature decreasing law : Tn ⫹ 1 ⫽ 0.98Tn, 3. Thermodynamic equilibrium criterion : the temperature is changed when one of the two following conditions are satisfied: – 12M accepted perturbations, – 100M tried perturbations (M being the number of parameters). 4. Stopping criterion: three successive temperature stages without any accepted perturbation.
3. Results and discussion The method has been tested through simulated signals. A realistic BAEP model [14] has been taken, shown in Fig. 2. Five waves were identified, the first three are due to the cochlear action potential, the cochlear nuclei and the protuberantial olive, respectively. The complex IV/V is generated by the lateral lemniscus and the inferior colliculus. The simulation has been applied on 100 identical BAEP waveforms, each BAEP was randomly delayed in two different ways : Gaussian delay distribution. Pre-determined delay curve.
3.1. Gaussian delay distribution BAEPs are perturbed according to a Gaussian distribution having a zero mean and a standard deviation of 0.3 ms. This standard deviation has been chosen by analysing its influence on BAEP distortion. It could be demonstrated that it is high enough to make the complex IV/V and waves II, III smoothed; latencies and conduction time are also modified. In this case it is clear that a direct interpretation by clinicians fails and an optimal correction is necessary. Fig. 3 shows BAEPs surface before the correction. BAEPs alignment achieved by SA is visually clear on Fig. 4, the waves are easily distinguished. The BAEP averaged before any correction is represented in Fig. 5, showing the distortion of the signal. In this case, we cannot state whether the second wave corresponds to wave II or to wave III. Latencies and conduction times couldn’t be properly measured. The BAEP extraction has been advantageously improved using the SA algorithm. As shown in Fig. 6, clinical parameters could be easily measured. In Table 2, we compare the averaged BAEP before correction and after correction, the BAEP model being taken as a reference. The difference could reach sometimes as much as 0.3 ms. 3.2. Pre-determined delay curve In the same simulation conditions, SA is applied on a delay curve having particular variations, shown in Fig. 7. It is not a typical case, since one does not know exactly how BAEPs moved from one stimulation to
The application of this method to real signals recorded from patients suffering from endocochlear pathologies is presented at the end of the section.
Fig. 2. Realistic BEAP model.
Fig. 3.
One hundred BAEPs randomly delayed before correction.
A. Naı¨t-Ali, P. Siarry / Medical Engineering & Physics 24 (2002) 385–392
Fig. 6.
Fig. 4.
BAEP alignment using SA technique.
389
Improvement of the BAEP using the SA algorithm.
Table 2 Comparison of latency and conduction time before and after correction BAEPModel
Latency I (ms) Latency II (ms) Latency III (ms) Latency V (ms) Conduction time I-III (ms) Conduction time I-V (ms) a
Before correction After correction
1.8 2.8 3.6 5.5 1.9
1.5
3.8
3.6
a a
5.1 a
1.7 2.8 3.7 5.4 1.7 3.7
Wave not identified
Fig. 5. Distortion of the BAEP signal averaged before any correction.
another, since a direct observation is impossible, due to the low SNR. Our particular aim in this simulation is to show how the optimal solution can be found with the SA technique independently from the delay curve. The same remark as in the previous case could be made concerning the problem of waves identification, latencies and conduction times related to non-corrected averaged BAEP. Signals before and after correction are shown respectively in Fig. 8 and Fig. 9. The correction performed using SA clearly allows better results to be obtained.
Fig. 7. Predetermined delay curve relative to the first signal, positive for right delay and negative for left delay.
390
A. Naı¨t-Ali, P. Siarry / Medical Engineering & Physics 24 (2002) 385–392
Fig. 8.
Distortion of the averaged BAEP.
Fig. 10. BAEP obtained by the classical averaging method for patient no. 1, suffering from an endocochlear pathology.
affirm if non-stationarity is of the type delay or not even if delays could be observed in some examples. A special study of signals coming from patients suffering from Me´ nie`re’s disease should be used. The same procedure should be employed for other pathologies. The BAEP shown in Fig. 10 (signal recorded from a patient 70 yr old with a presbyacousis and vertigo) has been obtained through the classical averaging method. We notice in this case that wave III is completely lacking and wave II has a low amplitude. The BAEP obtained using the SA method is shown in Fig. 11. The BAEP has been enhanced with a view of clinical interpretation, especially with the apparition of wave III. It is clear that the aim of this approach is not only to obtain a better BAEP shape, but to explore and understand the non-stationnarity that leads to a random delay of the response to stimulation. This allows us to establish a functional Fig. 9. Enhancement of the BAEP using the SA algorithm.
3.3. Application on real signals The SA method has been applied on real signals recorded from patients suffering from endocochlear pathologies1. The used database has been recorded by taking into account only patient symptoms as clinical parameters, no pathology was identified a priori, our aim for the moment being to verify the accuracy of some non-stationarity models. That is why we can’t actually associate our model to a specific pathology. To do this, we need to perform a statistical study for a long period and take into account the final result diagnosis. For example, in the case of Me´ nie`re’s disease, we can’t 1 Data recorded at ‘Centre d’Explorations Fonctionnelles Oto-neurologiques (CREFON)’, 10 rue Falguie`re 75015 Paris, France.
Fig. 11.
BAEP enhanced using SA method for patient no. 1.
A. Naı¨t-Ali, P. Siarry / Medical Engineering & Physics 24 (2002) 385–392
Fig. 12.
Histogram of delays for patient no. 1.
Fig. 14.
391
BAEP enhanced using SA method for patient no. 2.
exploration at the level of the cochlea. From the ‘optimal’ delay vector a histogram has been constructed in Fig. 12. This result is very important in the sense that it allows classification of the pathology. The classification problem is not studied in this paper. A second example is presented in Fig. 13 (signal recorded from a patient 55 yr old with a progressive deafness, tinnitus and recruitment) where the complex IV/V is not well estimated when the classical averaging method is used. The optimal correction shows a better BAEP shape in Fig. 14. By analysing the corresponding histogram in Fig. 15, we notice that the histogram is not Gaussian, delays are more concentrated on the right. We can conclude that in this case the cochlea has a dynamic different from that of the first example. Fig. 15. Histogram of delays for patient no. 2.
4. Conclusion
Fig. 13. BAEP obtained by the classical averaging method for patient no. 2, suffering from an endocochlear pathology.
The SA algorithm was applied on a set of randomly delayed BAEPs. Simulations have shown that using this technique allows improvement of the averaged BAEP extraction. Waves can be easily identified whatever sweeps are randomized. The enhancement of the signals means that we can know how BAEPs have been delayed; in consequence we are able to identify the dynamic of the cochlea. A classification step is very important in order to study the cochlea pathologies. We have shown that our method is thus suitable to estimate BAEPs recorded from ‘endocochlear patients’. In this paper, we did not study the performances of the technique versus the SNR. Indeed we considered, for a given M, that noise energy is negligible in comparison with BAEP energy. The simulation has been performed on 100 signals using
392
A. Naı¨t-Ali, P. Siarry / Medical Engineering & Physics 24 (2002) 385–392
a Pentium II computer at 350 MHz: few seconds are necessary for the optimization. Of course, the computing time non-linearly increases when increasing the number of signals to be processed. Ten years ago, SA algorithms were sometimes inefficient because of computing limitation, now computers are becoming faster, thus making SA significantly of interest for estimating BAEPs in the case of some pathologies.
[5] [6]
[7]
[8] [9]
References [1] Naı¨t-Ali A, Adam O, Motsch JF. The brainstem auditory evoked potentials estimation using a bayesian deconvolution method. In: Proceedings of the 18th Annual International Conference of IEEE Engineering in Medicine and Biology Society (CD-ROM), Amsterdam, November 1996. [2] Naı¨t-Ali A, Adam O, Motsch JF. On optimal aperiodic stimulation for brainstem auditory evoked potentials estimation. In: Proceedings of the 19th Annual International Conference of IEEE Engineering in Medicine and Biology Society (CD-ROM), Chicago, November 1997. [3] Abramson M, Stein B, Pedley T, Emerson R, Wazen J. Intraoperative BAER monitoring and hearing preservation in the treatment of acoustic neuromas. Laryngoscope 1985;95:1318–22. [4] Harper M, Harner S, Slavit D, Litchy W, Daube R, Beatty C,
[10]
[11] [12]
[13]
[14]
Ebersold M. Effect of BAEP monitoring on hearing preservation during acoustic neuroma resection. Neurology 1992;42:1551–3. Ohresser M, Toupet M, Vialla P, Strekers P. Les neurinomes a` BERA normaux. Ann Oto-laryng 1986;103:215–21. McGillem CD, Aunon JI, Pomalaza CA. Improved waveform estimation procedures for event related potentials. IEEE Trans. Biomed Eng 1985;32:371–9. Yu XH, Zhang YS, He ZY. Peak component latency-corrected average: method for evoked potential waveform. IEEE Trans. Biomed Eng 1994;41:1072–82. Kong X, Thakor NV. Adaptive estimation of latency changes in evoked potentials. IEEE Trans. Biomed Eng 1996;43:189–97. Kong X, Zhu L, Kaufman B, Rog A. Adaptive estimation of evoked potential latency under low SNR conditions. Proceedings of the 16th International Conference IEEE Engineering in Medicine and Biology Society, Baltimore, USA, November 1994, pp. 213–4. Kong X, Qiu T. Adaptive estimation of latency change in evoked potentials by direct least mean p-norm time delay estimation. IEEE Trans. Biomed. Eng 1999;99:4–1003. Kirkpatrick S, Gelatt CD, Vecchi MP. Optimization by simulated annealing. Science 1983;220(4598):671–80. Courat JP, Raynaud G, Mrad I, Siarry P. Electronic component model minimization based on Log Simulated Annealing. IEEE Trans Circuits Syst Part 1 1994;41(12):790–5. Siarry P, Berthiau G, Durbin F, Haussy J. Enhanced Simulated Annealing for globally minimizing functions of many continuous variables. ACM Trans Math Software 1997;23:209–28. Motsch JF. La dynamique du tronc ce´ re´ bral. The`se d’Etat, Universite´ Paris 12, 1987.
Application of simulated annealing for estimating BAEPs in endocochlear pathologies A. Naı¨t-Ali ∗, P. Siarry Universite´ Paris 12, LERISS 61, avenue du Ge´ne´ral de Gaulle, 94010 Cre´teil, France Received 19 July 2001; received in revised form 8 March 2002; accepted 19 March 2002
Abstract In this paper, we propose a new approach aimed at handling the temporal Brainstem Auditory Evoked Potentials (BAEPs) nonstationarity. It is pointed out that for some endocochlear pathologies, BAEPs could be randomly delayed from one response to another. This non-stationarity leads to smoothed BAEPs when applying ensemble averaging or any other technique based on BAEPs stationarity. In that case, waves identification is very difficult, sometimes impossible. The problem consists in estimating time delays. Knowing the distribution of delays allows subsequent study of the dynamic of the cochlea and, perhaps, identification of the nature of its pathology. The approach suggested in this paper is based on Simulated Annealing, used to minimize a non-linear criterion involving delays. This technique is advantageously compared to the non-corrected ensemble averaging method, using a set of simulated data based on a realistic model. As an illustration, results based on real signals recorded from two patients are presented and discussed at the end of the paper. 2002 IPEM. Published by Elsevier Science Ltd. All rights reserved. Keywords: Evoked potentials; Auditory system; Simulated annealing; Endocochlear
1. Introduction Brainstem Auditory Evoked Potentials (BAEPs) are electrical signals generated in response to acoustical stimulations of the auditory system [1,2]. BAEPs are recorded on the patient’s vertex and mastoid. They are used particularly for early diagnosis of acoustic neuroma [3–5]. The electroencephalogram (EEG) corresponding to a natural physiological activity corrupts the response to a single stimulation. Three main reasons cause BAEPs extraction difficulties. First, the very low signal-to-noise ratio (SNR) (its value is sometimes around ⫺20 dB); second, problems in EEG non-stationarity; and third, BAEP non-stationarity. Direct observation or extraction of BAEPs is not usually possible. One of the most classic techniques used to estimate BAEPs is the ‘ensemble averaging’, that performs well, especially when EEG and BAEPs are stationary. Various other techniques have been developed by taking into
Corresponding author. Tel: +33-1-4517-14-90; fax: +33-1-45171492. E-mail address: [email protected] (A. Naı¨t-Ali). ∗
account the non-stationarity properties. A well-known method, called the LCA technique (‘Latency Corrected Average’), has been proposed by McGillem et al. [6]. Latencies are detected at each sweep; signals having the same latencies are gathered. Yu et al. [7] have improved this technique by working out the PC-LCA (‘Peak Component Latency Corrected Average’): sweeps are filtered using an adaptive filter, before detecting latencies. Results are, of course, strongly correlated with the SNR. Other techniques are proposed in the literature to process late evoked potentials, as an adaptive estimation of latency change using a DLMS (‘Direct Least Mean Square’) adaptive TDE (‘Time Delay Estimation’) algorithm. Here the EEG is considered as Gaussian [8,9]. Another variant of this method was recently proposed by Kong et al. [10]. It is called DLMP (‘Direct LeastMean p-norm’) adaptive TDE, that works well under αstable noise conditions. In this paper, we are interested only in BAEPs temporal non-stationarity. BAEPs may be randomly delayed from one stimulation to another. The averaged signal will then be smoothed, as shown in Fig. 1. The idea here consists in keeping the classical technique of averaging, which is no longer applied to signals as they are
1350-4533/02/$22.00 2002 IPEM. Published by Elsevier Science Ltd. All rights reserved. PII: S 1 3 5 0 - 4 5 3 3 ( 0 2 ) 0 0 0 3 4 - 6
386
A. Naı¨t-Ali, P. Siarry / Medical Engineering & Physics 24 (2002) 385–392
Fig. 1.
Effects of BAEPs time delay on the averaged signal.
recorded, but after an optimal correction. Our aim is to estimate the delays which enhance the averaged signal. The estimated delays could later be analysed from the angle of distribution to understand the dynamic of the cochlea. However, filtering EEG will not be considered; the single assumption made on the EEG is that its energy tends to zero for M given number of sweeps. The paper comprises four sections. In section 2, we present the mathematical basis of our technique, that uses Simulated Annealing to estimate the delays and to improve the averaged signal. In section 3, we apply the method to simulated BAEPs, perturbed in two different ways; the results are then discussed. Results based on real signals recorded from two patients are finally presented and discussed. In section 4 a conclusion is given relating to the proposed method.
2. Mathematical basis of BAEPs estimating 2.1. Time delayed BAEP model Let xi(n) be the sampled signal recorded in response to the ith stimulus. One can write: xi(n) ⫽ s(n ⫹ Di) ⫹ vi(n) n ⫽ 0,1,…,N⫺1
the BAEP delay, randomly distributed according a given distribution; vi(n)is the EEG signal; N is the number of samples in a given response. The averaged signal x¯(n)over M sweeps is easily calculated as follows:
冘
M⫺1
x¯(n) ⫽
where : n is a discrete time variable; s(n) is the BAEP having a stationary waveform; Di is an integer measuring
冘
M⫺1
(2)
If we consider that M sweeps are sufficient to make the residual noise energy negligible in comparison with BAEP averaged energy, one can write:
冘
M⫺1
x¯(n)⬇
1 s(n ⫹ Di) Mi ⫽ 0
(3)
The technique that we propose consists in estimating delays maximizing the energy of the averaged signal. An endocochlear pathology is supposed to be the single cause of the observed delays. To make a diagnosis, it is necessary to study and classify those delays. More explanations are reported in section 3.3. The criterion to be optimized is formulated as follows :
冘冉 冘
N⫺1
(1)
冘
M⫺1
1 1 1 xi(n) ⫽ s(n ⫹ Di) ⫹ v (n) Mi ⫽ 0 Mi ⫽ 0 Mi ⫽ 0 i
fD ⫽
n⫽0
M⫺1
1 x (n⫺Di) Mi ⫽ 0 i
冊
2
with D ⫽ [D0,D1,D2,...DM⫺1].
(4)
A. Naı¨t-Ali, P. Siarry / Medical Engineering & Physics 24 (2002) 385–392
The first response is considered as a reference, then D0 ⫽ 0. Actually, this hypothesis ensures the unicity of the solution. Maximizing fD leads to minimizing over D the following criterion:
冘冉 冘
N⫺1
JD ⫽ ⫺
n⫽0
M⫺1
1 x (n⫺Di) Mi ⫽ 0 i
冊
2
(5)
The criterion (5) is not quadratic or convex, therefore optimization is driven using Simulated Annealing algorithm (SA) [11], which was successfully applied to numerous multiminima optimization problems [12] [13]. 2.2. Simulated Annealing algorithm The problem of finding the minimum of a given function J of M variables has received much attention for many years. Several algorithms that give good results for functions having only one minimum have been proposed. However, those search methods do not work satisfactorily when the cost function comprises several (‘global’ or ‘local’) minima within the domain of interest, because they tend to stop at the first minimum encountered. In many research and development fields, engineers have to deal with minimization problems of ever increasing dimensionality, leading to a combinatorial-like explosion of the number of local minima. For handling such global optimization problems, new approaches, called ‘metaheuristics’, have emerged, most of them being inspired by biological or physical phenomena : in particular, Simulated Annealing, Genetic Algorithms, Tabu Search and Ant Colonies. We propose here the use of the Simulated Annealing method, for which efficiency has been demonstrated in the solution of difficult optimization problems, when only the objective function is available. Simulated Annealing aims at solving such problems ‘at the best’ : in theory, it avoids the classical trapping into objective function local minima. It was inspired by the ‘annealing’ experimental technique, currently used by metallurgists to obtain a well-ordered solid, corresponding to the minimal energy state. The most significant features of the method are conditionally guaranteed convergence towards a global optimum, but in practice rather sensitive tuning of parameters. A pseudo-code of Simulated Annealing is given in Table 1. Starting from an initial point Dini, SA generates a succession of points D1, D2, etc.... A new point D’, randomly generated from the current point D (step 3 of Table 1), is submitted for acceptance to the Metropolis criterion, defined at the bottom of Table 1. The Metropolis algorithm allows SA to estimate first (at high temperature) the ‘coarse grain’ behaviour of the function under test and then to examine (at low temperature) its ‘finer grain’ behaviour.
387
Table 1 Pseudo-code of SA algorithm and Metropolis rule 1
Randomly choose randomly an initial solution Diniand evaluate the criterion value JD; 2 Choose an initial temperature T; 3 Perturb the solution to have a new one D’ ⫽ D ⫹ ⌬D; 4 Calculate ⌬J ⫽ JD⬘⫺JD; 5 The solution is accepted or rejected according to the Metropolis rule (see below) ; 6 If the thermodynamic equilibrium is reached, then decrease the temperature, else go to step 3; 7 If the system is ‘frozen’, go to step else go to step 3; 8 Solution = the best point found. Metropolis rule If ⌬J ⫽ JD⬘⫺JDⱕ0 then keep this solution and D’; D←D’; JD⬘⫺JD Else calculate P ⫽ exp{⫺ }; T Generate a random number R uniformly distributed within [0,1]; If RⱕP, then accept the new solution D’; D←D’, else reject the solution D’.
The sequence of points generated by SA at a fixed temperature T is stopped as soon as the thermodynamic equilibrium, detected by some adequate condition, is reached (step 6 of Table 1). Then the temperature is suitably adjusted, before beginning a new sequence of moves, and so on. This process is finally stopped (step 7 of Table 1) when the system is ‘frozen’, i.e. at least one of several independent stopping tests on downhill moves number, temperature and number of evaluations of J is satisfied. These tests are basically equivalent in the sense that they all aim at terminating the annealing process through estimation of J local behaviour: the SA search is stopped if J is sufficiently flat in the neighbourhood under exploration, as detected by at least one of the tests. To sum up, SA allows ‘uphill’ moves under control of the T temperature parameter. It is clear that, at higher temperature, only the coarse grain behaviour of the objective function J is relevant to the search. As temperature decreases, the finer grain behaviour of the objective function J is examined, finally leading to a solution theoretically near to the global minimum of J. A pseudo-code of SA and Metropolis criterion is given in Table 1. Several criteria exist to stop SA (step 7 of Table 1); for example, one may decide to stop SA when no more tested solutions are accepted, or when the temperature T is lower than a threshold temperature close to 0. The parameters of SA are fixed as follows: 1. Initial temperature T0 : it can be evaluated using the following empirical technique: 앫 Do 100 random perturbations: evaluate the average aute energy variation ⬍ 兩⌬E兩 ⬎ 앫 Choose and initial rate of acceptance t0, according to the confidence in the ‘quality’ of the initial configuration, for example:
388
A. Naı¨t-Ali, P. Siarry / Medical Engineering & Physics 24 (2002) 385–392
– Bad quality: t0 ⫽ 50% (beginning at high temperature), – Good quality: t0 ⫽ 20% (beginning at low temperature), 앫 DeduceT0 from therelation exp(⫺ ⬍ 兩⌬E兩 ⬎ / T0) ⫽ t0 2. Temperature decreasing law : Tn ⫹ 1 ⫽ 0.98Tn, 3. Thermodynamic equilibrium criterion : the temperature is changed when one of the two following conditions are satisfied: – 12M accepted perturbations, – 100M tried perturbations (M being the number of parameters). 4. Stopping criterion: three successive temperature stages without any accepted perturbation.
3. Results and discussion The method has been tested through simulated signals. A realistic BAEP model [14] has been taken, shown in Fig. 2. Five waves were identified, the first three are due to the cochlear action potential, the cochlear nuclei and the protuberantial olive, respectively. The complex IV/V is generated by the lateral lemniscus and the inferior colliculus. The simulation has been applied on 100 identical BAEP waveforms, each BAEP was randomly delayed in two different ways : Gaussian delay distribution. Pre-determined delay curve.
3.1. Gaussian delay distribution BAEPs are perturbed according to a Gaussian distribution having a zero mean and a standard deviation of 0.3 ms. This standard deviation has been chosen by analysing its influence on BAEP distortion. It could be demonstrated that it is high enough to make the complex IV/V and waves II, III smoothed; latencies and conduction time are also modified. In this case it is clear that a direct interpretation by clinicians fails and an optimal correction is necessary. Fig. 3 shows BAEPs surface before the correction. BAEPs alignment achieved by SA is visually clear on Fig. 4, the waves are easily distinguished. The BAEP averaged before any correction is represented in Fig. 5, showing the distortion of the signal. In this case, we cannot state whether the second wave corresponds to wave II or to wave III. Latencies and conduction times couldn’t be properly measured. The BAEP extraction has been advantageously improved using the SA algorithm. As shown in Fig. 6, clinical parameters could be easily measured. In Table 2, we compare the averaged BAEP before correction and after correction, the BAEP model being taken as a reference. The difference could reach sometimes as much as 0.3 ms. 3.2. Pre-determined delay curve In the same simulation conditions, SA is applied on a delay curve having particular variations, shown in Fig. 7. It is not a typical case, since one does not know exactly how BAEPs moved from one stimulation to
The application of this method to real signals recorded from patients suffering from endocochlear pathologies is presented at the end of the section.
Fig. 2. Realistic BEAP model.
Fig. 3.
One hundred BAEPs randomly delayed before correction.
A. Naı¨t-Ali, P. Siarry / Medical Engineering & Physics 24 (2002) 385–392
Fig. 6.
Fig. 4.
BAEP alignment using SA technique.
389
Improvement of the BAEP using the SA algorithm.
Table 2 Comparison of latency and conduction time before and after correction BAEPModel
Latency I (ms) Latency II (ms) Latency III (ms) Latency V (ms) Conduction time I-III (ms) Conduction time I-V (ms) a
Before correction After correction
1.8 2.8 3.6 5.5 1.9
1.5
3.8
3.6
a a
5.1 a
1.7 2.8 3.7 5.4 1.7 3.7
Wave not identified
Fig. 5. Distortion of the BAEP signal averaged before any correction.
another, since a direct observation is impossible, due to the low SNR. Our particular aim in this simulation is to show how the optimal solution can be found with the SA technique independently from the delay curve. The same remark as in the previous case could be made concerning the problem of waves identification, latencies and conduction times related to non-corrected averaged BAEP. Signals before and after correction are shown respectively in Fig. 8 and Fig. 9. The correction performed using SA clearly allows better results to be obtained.
Fig. 7. Predetermined delay curve relative to the first signal, positive for right delay and negative for left delay.
390
A. Naı¨t-Ali, P. Siarry / Medical Engineering & Physics 24 (2002) 385–392
Fig. 8.
Distortion of the averaged BAEP.
Fig. 10. BAEP obtained by the classical averaging method for patient no. 1, suffering from an endocochlear pathology.
affirm if non-stationarity is of the type delay or not even if delays could be observed in some examples. A special study of signals coming from patients suffering from Me´ nie`re’s disease should be used. The same procedure should be employed for other pathologies. The BAEP shown in Fig. 10 (signal recorded from a patient 70 yr old with a presbyacousis and vertigo) has been obtained through the classical averaging method. We notice in this case that wave III is completely lacking and wave II has a low amplitude. The BAEP obtained using the SA method is shown in Fig. 11. The BAEP has been enhanced with a view of clinical interpretation, especially with the apparition of wave III. It is clear that the aim of this approach is not only to obtain a better BAEP shape, but to explore and understand the non-stationnarity that leads to a random delay of the response to stimulation. This allows us to establish a functional Fig. 9. Enhancement of the BAEP using the SA algorithm.
3.3. Application on real signals The SA method has been applied on real signals recorded from patients suffering from endocochlear pathologies1. The used database has been recorded by taking into account only patient symptoms as clinical parameters, no pathology was identified a priori, our aim for the moment being to verify the accuracy of some non-stationarity models. That is why we can’t actually associate our model to a specific pathology. To do this, we need to perform a statistical study for a long period and take into account the final result diagnosis. For example, in the case of Me´ nie`re’s disease, we can’t 1 Data recorded at ‘Centre d’Explorations Fonctionnelles Oto-neurologiques (CREFON)’, 10 rue Falguie`re 75015 Paris, France.
Fig. 11.
BAEP enhanced using SA method for patient no. 1.
A. Naı¨t-Ali, P. Siarry / Medical Engineering & Physics 24 (2002) 385–392
Fig. 12.
Histogram of delays for patient no. 1.
Fig. 14.
391
BAEP enhanced using SA method for patient no. 2.
exploration at the level of the cochlea. From the ‘optimal’ delay vector a histogram has been constructed in Fig. 12. This result is very important in the sense that it allows classification of the pathology. The classification problem is not studied in this paper. A second example is presented in Fig. 13 (signal recorded from a patient 55 yr old with a progressive deafness, tinnitus and recruitment) where the complex IV/V is not well estimated when the classical averaging method is used. The optimal correction shows a better BAEP shape in Fig. 14. By analysing the corresponding histogram in Fig. 15, we notice that the histogram is not Gaussian, delays are more concentrated on the right. We can conclude that in this case the cochlea has a dynamic different from that of the first example. Fig. 15. Histogram of delays for patient no. 2.
4. Conclusion
Fig. 13. BAEP obtained by the classical averaging method for patient no. 2, suffering from an endocochlear pathology.
The SA algorithm was applied on a set of randomly delayed BAEPs. Simulations have shown that using this technique allows improvement of the averaged BAEP extraction. Waves can be easily identified whatever sweeps are randomized. The enhancement of the signals means that we can know how BAEPs have been delayed; in consequence we are able to identify the dynamic of the cochlea. A classification step is very important in order to study the cochlea pathologies. We have shown that our method is thus suitable to estimate BAEPs recorded from ‘endocochlear patients’. In this paper, we did not study the performances of the technique versus the SNR. Indeed we considered, for a given M, that noise energy is negligible in comparison with BAEP energy. The simulation has been performed on 100 signals using
392
A. Naı¨t-Ali, P. Siarry / Medical Engineering & Physics 24 (2002) 385–392
a Pentium II computer at 350 MHz: few seconds are necessary for the optimization. Of course, the computing time non-linearly increases when increasing the number of signals to be processed. Ten years ago, SA algorithms were sometimes inefficient because of computing limitation, now computers are becoming faster, thus making SA significantly of interest for estimating BAEPs in the case of some pathologies.
[5] [6]
[7]
[8] [9]
References [1] Naı¨t-Ali A, Adam O, Motsch JF. The brainstem auditory evoked potentials estimation using a bayesian deconvolution method. In: Proceedings of the 18th Annual International Conference of IEEE Engineering in Medicine and Biology Society (CD-ROM), Amsterdam, November 1996. [2] Naı¨t-Ali A, Adam O, Motsch JF. On optimal aperiodic stimulation for brainstem auditory evoked potentials estimation. In: Proceedings of the 19th Annual International Conference of IEEE Engineering in Medicine and Biology Society (CD-ROM), Chicago, November 1997. [3] Abramson M, Stein B, Pedley T, Emerson R, Wazen J. Intraoperative BAER monitoring and hearing preservation in the treatment of acoustic neuromas. Laryngoscope 1985;95:1318–22. [4] Harper M, Harner S, Slavit D, Litchy W, Daube R, Beatty C,
[10]
[11] [12]
[13]
[14]
Ebersold M. Effect of BAEP monitoring on hearing preservation during acoustic neuroma resection. Neurology 1992;42:1551–3. Ohresser M, Toupet M, Vialla P, Strekers P. Les neurinomes a` BERA normaux. Ann Oto-laryng 1986;103:215–21. McGillem CD, Aunon JI, Pomalaza CA. Improved waveform estimation procedures for event related potentials. IEEE Trans. Biomed Eng 1985;32:371–9. Yu XH, Zhang YS, He ZY. Peak component latency-corrected average: method for evoked potential waveform. IEEE Trans. Biomed Eng 1994;41:1072–82. Kong X, Thakor NV. Adaptive estimation of latency changes in evoked potentials. IEEE Trans. Biomed Eng 1996;43:189–97. Kong X, Zhu L, Kaufman B, Rog A. Adaptive estimation of evoked potential latency under low SNR conditions. Proceedings of the 16th International Conference IEEE Engineering in Medicine and Biology Society, Baltimore, USA, November 1994, pp. 213–4. Kong X, Qiu T. Adaptive estimation of latency change in evoked potentials by direct least mean p-norm time delay estimation. IEEE Trans. Biomed. Eng 1999;99:4–1003. Kirkpatrick S, Gelatt CD, Vecchi MP. Optimization by simulated annealing. Science 1983;220(4598):671–80. Courat JP, Raynaud G, Mrad I, Siarry P. Electronic component model minimization based on Log Simulated Annealing. IEEE Trans Circuits Syst Part 1 1994;41(12):790–5. Siarry P, Berthiau G, Durbin F, Haussy J. Enhanced Simulated Annealing for globally minimizing functions of many continuous variables. ACM Trans Math Software 1997;23:209–28. Motsch JF. La dynamique du tronc ce´ re´ bral. The`se d’Etat, Universite´ Paris 12, 1987.