Single trial somatosensory evoked potential extraction with ARX filtering for a combined spinal cord intraoperative neuromonitoring technique
 Lorenzo Rossi^{1, 2}Email author,
 Anna Maria Bianchi^{2},
 Anna Merzagora^{2, 3},
 Alberto Gaggiani^{1},
 Sergio Cerutti^{2} and
 Francesco Bracchi^{1}
https://doi.org/10.1186/1475925X62
© Rossi et al; licensee BioMed Central Ltd. 2007
Received: 29 June 2006
Accepted: 04 January 2007
Published: 04 January 2007
Abstract
Background
When spinal cord functional integrity is at risk during surgery, intraoperative neuromonitoring is recommended. Tibial Single Trial Somatosensory Evoked Potentials (SEPs) and Hreflex are here used in a combined neuromonitoring method: both signals monitor the spinal cord status, though involving different nervous pathways. However, SEPs express a trialtotrial variability that is difficult to track because of the intrinsic low signaltonoise ratio. For this reason single trial techniques are needed to extract SEPs from the background EEG.
Methods
The analysis is performed off line on data recorded in eight scoliosis surgery sessions during which the spinal cord was simultaneously monitored through classical SEPs and Hreflex responses elicited by the same tibial nerve electrical stimulation. The single trial extraction of SEPs from the background EEG is here performed through AutoRegressive filter with eXogenous input (ARX). The electroencephalographic recording can be modeled as the sum of the background EEG, which can be described as an autoregressive process not related to the stimulus, and the evoked potential (EP), which can be viewed as a filtered version of a reference signal related to the stimulus. The choice of the filter optimal orders is based on the Akaike Information Criterion (AIC). The reference signal used as exogenous input in the ARX model is a weighted average of the previous SEPs trials with exponential forgetting behavior.
Results
The moving average exponentially weighted, used as reference signal for the ARX model, shows a better sensibility than the standard moving average in tracking SEPs fast intertrial changes. The ability to promptly detect changes allows highlighting relations between waveform changes and surgical maneuvers. It also allows a comparative study with Hreflex trends: in particular, the two signals show different fall and recovery dynamics following stressful conditions for the spinal cord.
Conclusion
The ARX filter showed good performances in single trial SEP extraction, enhancing the available information concerning the current spinal cord status. Moreover, the comparison between SEPs and Hreflex showed that the two signals are affected by the same surgical maneuvers, even if they monitor the spinal cord through anatomically different pathways.
Background
The monitoring of the functionality of vital parameters is known to be a basic aspect in uptodate surgical techniques. This is particularly true in the case of surgery performed on the vertebral column (i.e.: correction of serious scoliosis), when a continuous monitoring of the spinal cord functionality is usually required [1–4].
In fact, during surgery the integrity and functionality of the spinal cord may be compromised by surgical maneuvers to such an extent as to cause prolonged insufficient blood supply to the cord, or mechanical compression and, as a consequence, its temporary malfunctioning. Sometimes this malfunctioning, if undetected at an early stage, can lead to irreversible spinal cord damage [5–7].
However, whereas the Hreflex analysis does not require refined processing methods, the SEPs due to their poor signal to noise ratio (SNR), require proper tools for the analysis of the EEG in order to extract the evoked potential waveform.
The current clinical SEP extraction techniques usually average several trials in order to enhance the SNR. However, the SEP has a trialtotrial variability which is lost if it is extracted using the averaging technique. This would lead to the loss of information about the trialtotrial changes of the SEP signal.
When combining SEPs and Hreflex, a limit is introduced by the long recovery time of the Hreflex, so that the interstimulus interval (ISI) must be at least 10 s [15][16]. In case of a SEPs' sudden changes the averaging technique provides intolerably delayed information, for this reason, a single trial SEP analysis is needed in order to compare it with the Hreflex obtained from the same electrical stimulus.
Methods
Clinical protocol
Surgery sessions monitored
Age  Sex  Level of malformation  Surgery correction technique 

31  F  D5L4  CD 
13  F  D5L1  HR 
17  M  D3L3  HR 
20  F  D5L3  CD 
14  F  D4L4  CD 
17  F  D3L1  CD 
16  F  D3L4  CD 
53  F  D4L5  CD 
EEG trials for the evoked potentials and EMG for the Hreflex were acquired during surgery form each of eight patients. The signals were elicited stimulating the posterior tibial nerve with a 1 ms rectangular voltage pulse up to 150 V at 0.1 Hz repetition rate.
The EEG signal was recorded through needle electrodes placed on the scalp in Cz' with an auricular reference. The Hreflex was acquired using a surface electrode (RedDot, 3 M, USA) on the soleus muscle and a reference one over the Achilles tendon.
The EEG and EMG acquisition were performed by means of two analog amplifilters ICP 511AP by Grass Instruments, with different gains and filter settings (EEG: gain = 80dB, bandpass 1 Hz – 300 Hz; EMG: gain = 40 dB, bandpass 30 Hz – 3 kHz).
The stimulus triggered EMG recording lasts 62.5 ms, sampled at 10 kHz; the EEG recording lasts 125 ms, sampled at 2.5 kHz by PCI board (type PCIMIO16E4, National Instruments, USA). The acquisition board range is +/1.25 V with a 12bit resolution.
The analysis software used for offline processing the recorded data has been implemented by means of Virtual Instruments developed in LabView with integrated Matlab processing scripts. The software performs the single trial extraction of SEP waveform and the detection of amplitude and latency both of the SEP and of the HReflex simulating an online neuromonitoring system.
EEG signal processing
ARX model
Many algorithms have been proposed in order to obtain a single trial SEP from the EEG [17–22] or for improving the signal to noise ratio for neuromonitoring application.
An Autoregressive Model with an eXogenous Input (ARX) has been chosen, because it has been successfully applied and tested in the analysis of different kinds of evoked potentials (visual, acoustic, somatosensory [23–25]).
According to the ARX model, the signal y(k), recorded after a somatosensory stimulation, can be modeled as the sum of two different contributions: a former not related to the stimulus (the background EEG) that, according to the literature [26], can be described as an autoregressive process driven by a white noise e(k); a latter related to the stimulus (the EP), that can be viewed as a filtered version of a reference signal u(k). The general equation of an ARX process is:
$y(k)={\displaystyle \sum _{i=1}^{n}{a}_{i}}\cdot y(ki)+{\displaystyle \sum _{j=d}^{m+d1}bj\cdot u(kj)}+e(k)\left(1\right)$
The model orders (n, m, d) are defined as follows:

n and m, that are respectively the model orders of the autoregressive and of the moving average sections;

d, takes into account the possible temporal anticipation between the reference input and the output of the filter (the filter is not causal).
In the ztransform domain, where z^{1} is the unit delay operator, the ARX model is expressed as:
A (z)·Y = B (z)·U + E (2)
and
$Y=\frac{B(z)}{A(z)}\cdot U+\frac{E}{A(z)}\left(3\right)$
If u(k) and y(k) are known, then the a _{ i }and b _{ j }coefficients can be estimated together with the variance of the input white noise, for every prefixed model order set and delay (n, m, d), by means of a least squares approach, that minimizes the following figure of merit:
$J=\frac{1}{N}{\displaystyle \sum _{1}^{N}{\left[err(k)\right]}^{2}}\left(4\right)$
where err(k) is the prediction error of the model and N is the number of samples in the trial. The prediction error is the difference between the recorded signal and its prediction according the model expressed in eq. 1.
Choice of the model order set (n, m, d)

Choice of an initial set of model orders (n, m, d);

Minimization of a figure of merit J in order to find the best ARX coefficients for the chosen order.

Check the whiteness of prediction error err(k): the set of identified parameters is accepted for the following calculations when the err(k) is a white sequence with a confidence level of 95%; the whiteness of the residual is tested using the cumulative Anderson test.

If the whiteness test is not verified, the model orders are increased until the prediction error is a white noise.
The starting triplet (n, m, d), selected in the first step, is calculated for each subject before starting the monitoring phase of the risky surgical procedures.
The first 50 sweeps are recorded and stored, and the optimal model orders, among the ones that satisfy the whiteness criterion, are selected by minimizing the Akaike Information Criterion (AIC) function; though there are several methods for finding the optimal orders for model (Final Prediction Error, Minimum Description Length, etc.), for the specific application the best one is the AIC because where the model orders are low the AIC shows better performances [27].
The AIC function is expressed as:
$aic=\mathrm{ln}\left({\sigma}^{2}\right)+\frac{n+m}{N}\left(5\right)$
where σ^{2} is the variance of the prediction error err(k).
The d parameter is chosen as m/2 considering a previous investigation [18].
On the basis of the AIC function, evaluated on the first 50 recorded sweeps, the triplet (n, m, d) is selected for the specific patient.
The reference signal
In many applications of the ARX filter for the extraction of the single trial EP, the reference u(k) is the average of a sufficient number of trials, representing a standard EP shape. The choice of the reference signal is fundamental for the ARX filter capability in tracking the intertrials EP variations: in fact the filter tries to identify a waveform, similar to the input reference, in the recorded signal trial y(t). Thus, a fixed reference could be misleading in the case of a dynamic process. For this reason a fast adapting reference could be a better choice in order to be able to track also sudden changes in the waveform.
However, the ability of this reference to adapt to morphological variations in the EP wave depends on the number of sweeps included: the moving average introduces a delay in the updating of the reference that depends on the number of averaged trials. In order to enhance the adaptation ability of the reference signal, a moving average with an exponential decay is here proposed and tested.
The reference is computed as follows:
u _{ i }= μ·u _{ i  1}+ (1  μ)·y _{ i } (6)
The μ coefficient (forgetting coefficient) provides an exponential windowing on the averaged sweeps, so that only the most recent trials contribute significantly to the reference u _{ i }, while the oldest ones are progressively forgotten, with an exponential decay. Different values of μ have been tested in order to have a reference signal that adaptively varies according to the possible dynamical changes in the evoked potentials. Offline tests led to the value μ = 0.95 [28].
In order to compare the ARX performances with the reference built on a static average and with reference built on a moving average driven by a forgetting coefficient, as defined by the eq. (6), signal processing simulations have been performed.
Three known input signals, representing the x _{ i }(k) target signals, have been built from one generic EP that has been multiplied for 3 different coefficients: k _{ 1 }= 20, k _{ 2 }= 15, k _{ 3 }= 1 in order to simulate different amplitudes.
An AR model has been identified from an EEG trial of 125 ms not related to any stimulus. A white Gaussian noise has been utilized as the AR model input in order to generate EEG signals n _{ i }(k). With this method it was possible to collect, as the sum of random n _{ i }(k) and known targets signals k_{j}*x _{ i }(k), 50 trials using the k _{ 1 }coefficient, 1 using the k _{ 2 }and 1 using the k _{ 3 }. These input signals have been utilized for simulating a possible variation of the EP due to a stressful maneuver applied at the 51^{st} trial that produces a strong EP decrease at the 52^{nd} trial.
The ARX filter using the two different reference signals, i.e. the moving average and the moving average with forgetting coefficient, have been compared according to their ability to extract the input signals k_{j}*x _{ i }(k).
Correlation with the Hreflex amplitude
The degree of correlation between the amplitude of SEP and Hreflex has been computed through the Pearson's correlation coefficient in the 8 surgery sessions studied.
The Pearson's coefficient is defined as:
$r=\frac{\frac{{\displaystyle \sum (AsA\tilde{s})\cdot (AhA\tilde{h})}}{n1}}{Ss\cdot Sh}\left(7\right)$
where A_{s} and A_{h} are the amplitudes of SEP and Hreflex, A $\tilde{s}$ and A $\tilde{h}$ are the averages of the respective amplitudes, S_{s} and S_{h} are the standard deviations and n is the number of the considered trials. Pearson's coefficient value indicates the correlation level for a specific freedom degree given by the number of trials analyzed; a coefficient value of 1 indicates a perfect correlation whereas a value of 0 indicates no correlation.
Data analysis protocol
The sequence of the offline data analysis steps have been structured in order to comply with the requirements for a data processing during a surgery session.

Cleaning corrupted recorded trials
 1.
presence of saturated samples: if there is at least one sample with a value of +/1.25 V (the maximum/minimum acquisition board voltage input), the trial is rejected;
 2.
singletrial amplitude range: given the maximum sample value (Max) and the minimum sample value (Min) of the ithtrial, if MaxMin>0.7 V, the trial is rejected.

Building the initial reference signal

Choosing the optimal order for the model
The trial set, used to build the initial reference signal, is also used to evaluate the best model order. The following range of model orders have been used in order to calculate the AIC function:

Single sweep analysis.
The trials acquired during the surgical procedure are evaluated in combination with the Hreflex signal.
Results
Simulations
However the difference due to the reference is clearly visible in the extraction of the 51^{st} trial. The 51^{st} raw trial, as described above, is built from an EP reduced by 20% in respect to the EP of the 50^{th} trial. Utilizing the moving average with a forgetting coefficient, it is possible to notice a variation of extracted EP of about 50% on P30 peak amplitude whereas using as reference the average of the previous 50 trails the peak shows a decrease of only 5%.
Degree of correlation between SEP and Hreflex obtained computing the Pearson's coefficient.
age  sex  level  r  df  p = 0.05  P = 0.01 

31  F  D5L4  0.501  31  0.351  0.447 
13  F  D5L1  0.551  30  0.355  0.456 
17  M  D3L3  0.746  36  0.321  0.413 
20  F  D5L3  0.611  47  0.282  0.365 
14  F  D4L4  0.580  30  0.355  0.456 
17  F  D3L1  0.617  40  0.304  0.393 
16  F  D3L4  0.489  41  0.299  0.386 
53  F  D4L5  0.385  42  0.294  0.379 
A surgery case
The procedure followed in the signal processing method is here outlined using data obtained from one of the eight cases presented in Table 1.
Surgery history
Time  Maneuvers 

8:44  Incision 
10:31  Hammered 
10:50  Hammered 
11:13  Derotation 
11:18  Decortications 
11:58  Suture 
12:12  End 
In step 1 and step 2 there is no detectable changes in the monitored waveforms as they are recorded during the initial surgery phase and no stressful maneuvers have been performed yet.
At step 3 the vertebral column begins to be hammered and both signal amplitudes decrease considerably. The two signals, however, show different behaviors. The Hreflex recovers very quickly, whereas the SEP shows a slower recovery dynamic. The same behavior can be observed after other stressful maneuvers, as decortications and tractions.
Conclusions and discussion
The application target was to implement a signal processing method able to monitor single trial SEPs during surgery and to be very sensitive to the evoked potential variations. The ARX filter, using a reference generated by a recursive exponential averaging, fulfills the requirements and provides better performances than the classical moving average reference. In order to test the method we used a lower high pass filter (1 Hz) than the one utilized in usual neuromonitoring applications (30 Hz) in order to test the signal processing method in less favorable conditions and demonstrating its capability to detect also small single potential variations.
The paper presents offline results, but the method is designed for online application. During surgery, the time available to complete a single trial extraction is 10 sec, i.e the interstimuli interval. The proposed algorithm completes the analysis in less than 1 sec (running on Laptop with Intel P4 512MB) so it is suitable for online application. Another critical point is the choice of the model orders (n, m, d). The intersubject variability imposes the calculation of the optimal orders for each new patient as described in the method section and this procedure has to be completed between the beginning of the surgery, when the patient falls asleep, and the starting of the surgical risky maneuvers (about 15 minutes). Actually this procedure is completed in less than 2 minute because, as described in the method section, it needs only 50 trials obtained at higher stimulation frequency (5 Hz). There is also some intraindividual variability, even if less critical than intersubject variability. For this reason the test of the model is carried out each trial through the whiteness test of the prediction error err(k). In the event of a not whiteness of the err(k) the model orders are increased.
Other signal processing methods [29, 30] have been proposed in literature for SEP monitoring during surgery, but the presented application compared the single trial SEP with the Hreflex using the same stimulus delivered at very low frequency. For this reason we develop a signal processing method, based on previous works [22–26], which gives information at each single trial recorded. The method improves the sensibility to the SEP variation as presented in the result section. Though 8 surgery sessions are not enough to have clinical and medical results, some comments are still possible, even at this early stage of research. The analysis of these 8 clinical cases recorded engenders a general consideration about the correlation between SEP and Hreflex amplitudes that are statistically correlated as shown by the Pearson's coefficient. As expected, the Pearson's coefficient shows a fair correlation between the amplitude of the single trial SEP and the Hreflex. As described previously, the behavior of the two signals is quite different regarding the response of the spinal cord to a stressful maneuver. However, a correlation exists between the EP and Hreflex. The example above, in fact, shows that the same maneuver affects both signals. The presence of correlation between the changes that occur in both signals means that the changes are related to the same cause: in this case, this could mean that modification of the spinal cord functionality is induced by the same surgical maneuvers. Thus, the two signals, anatomically and functionally independent, generated by a single transcutaneous stimulus, are both related to the functionality of the spinal cord and are both affected by the same cause. Regarding the neuromonitoring technique where the method has been applied, the combined neuromonitoring has shown the expected behavior: even if SEPs and Hreflex monitor anatomically different pathways, when they are affected by the same surgical maneuvers, they show significantly correlated amplitude changes.
A larger set of monitored surgeries will be collected in order to obtain a strong and significant clinical result that confirms a higher reliability of this new combined monitoring technique.
Declarations
Authors’ Affiliations
References
 Guérit JM: Neuromonitoring in the operating room: why, when, and how to monitor? Electroencephalography & Clinical Neurophysiology 1998, 106: 1–21. 10.1016/S00134694(97)000771View ArticleGoogle Scholar
 Laureau E, Marciniak B, Hebrad A, Herbaux B, Guieu JD: Neuromonitoring and anesthesia in surgery of the spine. Clin Neurophysiol 1998, 28: 299–320. 10.1016/S09877053(98)800026View ArticleGoogle Scholar
 Minahan RE: Intraoperative neuromonitoring. Neurologist 2002, 8: 209–226. 10.1097/0012789320020700000001View ArticleGoogle Scholar
 Guérit JM, Dion RA: Stateofart of neuromonitoring for prevention of immediate and delayed paraplegia in thoracic and thoracoabdominal aorta surgery. Ann Thorac Surg 2002, 74: 1867–1869. 10.1016/S00034975(02)041309View ArticleGoogle Scholar
 Lin J, Zhang J, Ye Q, Sheng J, Qiu G: Failure and complication following surgical treatment of scoliosis – an analysis of 101 cases. Chin Med Sci J 1999, 14: 174–179.Google Scholar
 Cervellati S, Bettini N, Bianco T, Parisini P: Neurological complications in segmental spinal instrumentation: analysis of 750 patients. Eur Spine J 1996, 5: 161–166. 10.1007/BF00395507View ArticleGoogle Scholar
 Wiedemayer H, Fauser B, Sandalcioglu IE, Schafer H, Stolke D: The impact of neurophysiological intraoperative monitoring on surgical decisions: a critical analysis of 423 cases. J Neurosurg Spine 2000, 96: 255–262.View ArticleGoogle Scholar
 Bracchi F, Rossi L, Merzagora AC: Apparatus and Methods for Monitoring the Functionality of Spinal Cord. International Patent Pending: PCT/IT2004/000440 Google Scholar
 Nuwer MR: Spinal Cord Monitoring With Somatosensory Techniques. J Clin Neurophysiol 1998, 15: 183–193. 10.1097/0000469119980500000002View ArticleGoogle Scholar
 Lille F, Petit B, Margules S, Mazel C, RoyCamille R: Somatosensory evoked potentials during spinal orthopedic surgery in adult patients. Neurophysiol Clin 1993, 23: 179–92. 10.1016/S09877053(05)802291View ArticleGoogle Scholar
 Bracchi F, Grossi P, Trovati L, Viganò P: HReflex spinal cord monitoring during vertebral column stabilization surgery. In Handbook of Spinal Cord Monitoring. Edited by: Boyd J. London: Kluwer; 1992:253–258.Google Scholar
 Bracchi F, Guerrieri M, Frassi P, Vigano P, Robinson MA: HReflex intraoperative spinal cord monitoring. In Neurological Complications of Spinal Surgery Edited by: G.I.C.D. Paris. 1994, 113–119.Google Scholar
 Pelosi L, Lamb J, Grevitt M, Mehdian SM, Web JK, Blumhardt LD: Combined monitoring of motor and somatosensory evoked potential in orthopaedic spinal surgery. Clin Neurophysiol 2002, 113: 1028–1091. 10.1016/S13882457(02)000275View ArticleGoogle Scholar
 Slimp JC: Electrophysiologic intraoperative monitoring for spine procedures. Phys Med Rehabil Clin N Am 2004, 15: 85–105.View ArticleGoogle Scholar
 Balzer JR, Rose RD, Welch WC, Sclabassi RJ, Robert JMD: Simultaneous somatosensory evoked potential and electromyographic recording during lumbosacral decompression and instrumentation. Neurosurgery 1998, 42: 1318–1324. 10.1097/0000612319980600000074View ArticleGoogle Scholar
 Schieppati M: The Hoffmann Reflex: a means of assessing spinal reflex excitability and its descending control in man. Prog Neurobiol 1987, 28: 345–376. 10.1016/03010082(87)900074View ArticleGoogle Scholar
 Fung KSM, Chan FHY, Lam FK, Poon PWF: A tracing evoked potential estimator. Med Biol Eng Comput 1999, 37: 218–227. 10.1007/BF02513290View ArticleGoogle Scholar
 Cerutti S, Baselli G, Liberati D, Pavesi G: Single sweep analysis of visual evoked potentials through a model of parametric identification. Biol Cybern 1987, 56: 111–120. 10.1007/BF00317986View ArticleGoogle Scholar
 Laguna P, Meste O, Poon PW, Caminal P, Rix H, Thakor NT: Adaptive filter for eventrelated bioelectric signals using an impulse correlated reference input: comparison with signal averaging technique. IEEE Trans Biomed Eng 1992, 39: 1032–1044. 10.1109/10.161335View ArticleGoogle Scholar
 Walter DO: A posteriori "Wiener filtering" of average evoked response. Electroencephalography & Clinical Neurophysiology 1968, (Suppl 27):61–70.Google Scholar
 Heinze HJ, Kunkel H: ARMAfiltering of evoked potentials. Methods Inf Med 1984, 33: 29–36.Google Scholar
 Spreckelsen MV, Bromm B: Estimation of singleevoked cerebral potentials by means of parametric modeling and Kalman filtering. IEEE Trans Biomed Eng 1988, 35: 691–700. 10.1109/10.7270View ArticleGoogle Scholar
 Mainardi LT: Single sweep analysis of event related auditory potential for the monitoring of sedation in cardiac surgery patients. Comput Methods Programs Biomed 2000, 63: 219–227. 10.1016/S01692607(00)001127View ArticleGoogle Scholar
 Magni R, Giunti S, Bianchi A, Reni G, Bandello F, Durante A, Cerutti S, Brancato R: Singlesweep analysis using an autoregressive with exogenous input (ARX) model. Doc Ophthalmol 1994, 86: 95–104. 10.1007/BF01224631View ArticleGoogle Scholar
 Liberati D, Dicorrado S, Mandelli S: Topographic mapping of single sweep evoked potentials in the brain. IEEE Trans Biomed Eng 1992, 39: 943–51. 10.1109/10.256428View ArticleGoogle Scholar
 Cerutti S, Liberati D, Mascellani P: Autoregressive modelling and filtering of EEG signal generating mechanism in patients under surgical interventions. Model Simulat Contr 1985, 2: 11–23.Google Scholar
 Haykin S: Adaptive filter theory. Prentice Hall; 2002.Google Scholar
 Bracchi F, Perale G, Rossi L, Gaggiani A, Bianchi A: A PCBased System for HReflex and Single Sweep SEP Coupled Monitoring of Spinal Cord Function in Vertebral Column Surgery. Proceedings of the 25th Annual International Conference of IEEE Engineering in Medicine and Biology Society: 17–21 September 2003; Cancun Google Scholar
 Lam BSC, Hu Y, Lu WW, Luk KDK, Chang CQ, Qui W, Chan FHY: Multiadaptive filtering technique for surface somatosensory evoked potentials processing. Med Eng Phys 2007, 27: 257–266. 10.1016/j.medengphy.2004.09.007View ArticleGoogle Scholar
 Qui W, Fung KS, Chan FHY, Lam FK, Poon PW, Hamernik RP: Adaptive filtering of evoked potential with radialbasisfunction neural network prefilter. IEEE Trans Biomed Eng 2002, 49: 255–232.Google Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.