Efficient parameterization of MUAP signal for identification of MU and volume conductor characteristics using neural networks

Journal of Neuroscience Methods 164 (2007) 325–338

Efficient parameterization of MUAP signal for identification of MU and volume conductor characteristics using neural networks A. Reffad, R.E. Bekka ∗ , K. Mebarkia, D. Chikouche Electronics Department, Engineering Faculty, Ferhat ABBAS University of Setif, 19000 Setif, Algeria Received 14 February 2007; received in revised form 15 April 2007; accepted 22 April 2007

Abstract The motor unit action potential (MUAPs) shapes depend on the anatomy and the physiology of the contracted muscle. The aim of this work is the identification of some characteristics of the motor unit (MU) and the volume conductor, namely the MU depth, the innervation zone width and the thickness of fat and skin layers based on MUAP signal parameters. The relationship between these characteristics and MUAP parameters are non-linear and complex. Thus, the use of the neural networks approach becomes an efficient tool to put in evidence this relationship. We have used the similarity and the homogeneity of the parameter criterions to choose which parameters are appropriate for the extraction. Two identification systems are presented and compared, a global system and a separate one. In order to evaluate the performance of each system, we have tested them using several simulated MUAP signals corrupted with additive Gaussian noise at different signal to noise ratios (SNR). A new test is introduced in which the electrode radius, the bar electrode dimensions and inclination angles for the detection system, fixed during the training process, are changed. © 2007 Elsevier B.V. All rights reserved. Keywords: Electromyography; Identification; Motor unit action potential; Parameterization; Neural networks

1. Introduction In order to better understand the relationship between internal anatomical and physiological characteristics and external parameters extracted from the recorded signals, researchers developed several models which allow the use of different methods to estimate the parameters of surface EMG signal recorded non-invasively on the skin by means of appropriate electrodes. The EMG signal is the summation of all action potentials of active motor units (MUs) within the detection volume. A MU consists of a motoneuron and the muscle fibers it innervates. The sources of the electrical potential are the depolarized zones of the active muscle fiber. Two intracellular action potentials (IAPs) are generated at the neuromuscular junction (NMJ), then the corresponding depolarized zones travel in two opposite directions towards the tendons with a propagation velocity. The electric activity of the muscle fibers, belonging to MU, generates the motor unit action potential (MUAP) (Basmajian and De Luca,

∗

Corresponding author. Tel.: +213 36 92 84 18; fax: +213 36 92 84 18. E-mail addresses: [email protected] (A. Reffad), bekka [email protected] (R.E. Bekka). 0165-0270/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.jneumeth.2007.04.015

1985; Merletti et al., 2003). The anatomical and physiological features of the muscle (for example depth of the fibers, number of fibers, conduction velocity of the fiber etc.) can alter the MUAP shape and duration (Lindstrom and Magnusson, 1977; Stulen and De Luca, 1981; Merletti et al., 1992, 1999a,b). These variations have been estimated by several parameters (Stegeman and Linssen, 1992; Roeleveld et al., 1997a,b; Disselhorst-Klug et al., 1998). It had been shown that the surface MU potential depends on the MU depth, the electrode impedance, the position of the active MU under the skin surface and the thickness of the subcutaneous fat layer (Roeleveld et al., 1997c,d). Lateva et al. (1996) showed that changing the position of the recording electrode in the muscle fiber direction causes slight changes in the MUAP amplitude and they observed a rapid decrease in the recorded MUAP amplitude when the electrode position is moved over the skin surface in the perpendicular direction to the muscle fiber. The signal has been reduced to a set of discriminative parameters (Kelly et al., 1990) in order to facilitate the identification process of characteristics. To give more information about the relationship between the parameters of the MUAP signal and the MU characteristics, Merletti et al. (1999a) have estimated both spectral and amplitude variables from surface signals and displayed them

326

A. Reffad et al. / Journal of Neuroscience Methods 164 (2007) 325–338

as functions of the model parameters. The association of such SEMG variables to muscle characteristics is crucial because it is used in identification systems such as those involved in the diagnosis of neuromuscular disorders (Rau and Disselhorst-Klug, 1997; Hogrel, 2005). Bekka et al. (2002) have proposed an artificial neural network approach to estimate MU characteristics from a single surface MUAP. The developed system is dedicated to identify five muscle characteristics of the planar Merletti’s model: the neuromuscular junction (NMJ) width, the MU territory radius, the MU depth, the anisotropy ratio, and the number of fibers. Their identification system is composed of four modules. The characteristics estimation stage uses four multilayer neural networks trained on simulated MAUPs which are generated from a planar Merletti model. They have also shown that the artificial neural technique performs well when the signal to noise ratio is greater than 20 dB. However, their study has not considered a model with subcutaneous layers and the signal parameterization has not been considered. In order to find the optimal combinations of the identification systems input parameters that can provide the best performance of characteristics estimation, we propose a new parameterization of double differential MUAP signal. This method is performed in two steps: first, we eliminate the non-sensitive parameters in the individual behavior, then the remaining parameters are selected by using the similarity and homogeneity criterions to reject redundant parameters and select the representative ones. Two MU characteristics (the MU depth h and the innervation zone width wi) and two volume conductor characteristics (the thickness of fat layer tf and the thickness of skin layer ts) have been considered. The extraction method is based on the identification principle where an artificial neural network (ANN) approach is proposed in order to link these characteristics with the double differential MUAP signal parameters. For this purpose, two neural networks identification systems have been considered and compared, a global identification system (GIS) and a separate one (SIS). According to Bekka approach (Bekka et al., 2002), the first system globally extracts the characteristics of the MU and the volume conductor, i.e. these characteristics are extracted at once by the same neural module. The second one extracts these characteristics separately using four neural modules. The SIS system has been proposed to overcome training problems encountered in the GIS system. The performance of the identification systems has been evaluated on their generalization ability and their sensitivity to the noise. Thus, the paper proposes two identification systems, GIS and SIS, to extract information about the MU and the volume conductor that can be derived from surface MUAP signals. So, the key issue is the number and type of parameters to be considered as input for identification systems. 2. Simulation of the MUAP The MUAP is simulated using the analytical model with multilayer cylindrical description of the volume conductor proposed by Farina et al. (2004a). In this model, bone, muscle, fat, and skin layers have been considered. Except the external layer (air layer), all other layers are limited in radial direction. All the

Table 1 Parameter values used in simulations Description External radius of the model (mm) Radius of the bone layer (mm) Angle of inclination of the MU with respect to the detection system (◦ ) Distance between the MU and the geometrical center of the detection system in angular direction (◦ ) Conduction velocity (m/s) Sampling frequency Radius of the recording circular electrodes (mm) Inter-electrode distance (mm) Number of fibers Radius of MU territory (mm) Position of the first detection point (mm) Lower semi-fiber length Upper semi-fiber length (mm) Spread of the lower tendon region (mm) Spread of the upper tendon region (mm)

Value 50 20 0 0

4 4096 1 10 50 1 50 80 130 10 10

layers are isotropic, except the muscle in which conductivities are different in both the radial (σ θ = σ ρ ) and axial (σ z ) directions. The conductivity of the muscle tissue is higher in the axial direction than in the radial one. According to Roeleveld et al. work (1997c), we have assumed that (a) the air layer has zero conductivity, (b) σ θ = σ ρ = 0.1 ( m)−1 and σ z = 0.5 ( m)−1 for the muscle layer, (c) σ θ = σ ρ = σ z = 1 ( m)−1 for the skin layer and (d) σ θ = σ ρ = σ z = 1 ( m)−1 for the bone layer. The simulated potentials are given by 1D inverse Fourier transform of a function in the temporal frequency domain obtained by spatial and temporal filtering operations on the current source density. The MU is modeled as a group of parallel fibers normally distributed inside a cylinder of radius 1 mm whose axis is parallel to the longitudinal double differential (LDD) detection system. The MU is centered at a depth h below the fat layer. The NMJ’s of the individual fibers are normally distributed within an innervation zone of width wi. Terminations are also normally distributed within tendon regions. The model parameters used in our simulations are given in Table 1. 3. Principle of identification Many factors affect the MUAP morphology, essentially the characteristics of the volume conductor and the muscle, the electrodes form, size and position, and the acquisition or recording conditions (Gydikov et al., 1986; Farina et al., 2001a, 2002, 2004b; Helal and Bouissou, 1992). The extraction procedure of the MU and volume conductor characteristics from the parameters of the modeled MUAP signal is called identification by which the parameters used to simulate the MUAP signal will be identified (Fig. 1). We know that MUAP morphology changes with MU and volume conductor characteristics but we do not know how. This change in the morphology can be observed by some signal parameters as amplitude, peak to peak value, root mean square, area, duration and by several other ones (Stegeman and

A. Reffad et al. / Journal of Neuroscience Methods 164 (2007) 325–338

327

nal parameter, the case leading to the multi-solution problem. As an example, different motor units having the same depth provide different values of the same parameter. It is clear that this dependent nature between the characteristics and suggested parameters complicates the identification task. Fig. 2 shows the variation of four parameters of the MUAP signal according to the MU depth h, the innervation zone width Fig. 1. Identification system principle.

Linssen, 1992; Merletti et al., 1992, 1999a). The objective of this section is the parameterization of the MUAP signal and the presentation of the selected parameters that can be considered as the most efficient entries to the identification systems and the best representative of the desired information contained in the MUAP signal. 3.1. MUAP parameterization The EMG signal parameterization is exploited in different applications, such as classification, diagnosis, prostheses control and recognition pattern (Farina et al., 2001b; Kohn, 1989; Hudgins et al., 1994; Pattichis et al., 1995; Englehart and Hudgins, 2003; Parker et al., 2006). The ordinarily used parameters are the average rectified value (ARV), the root mean square value (RMS), the peak to peak value (PP), the median and mean frequency (MDF, MNF) of the myoelectric signal of power spectral density function (De Luca, 1979; Stegeman and Linssen, 1992; Merletti et al., 1992; Merletti and Parker, 1998; Staudenmann et al., 2005; Holtermann et al., 2005; Mesin et al., 2006). Roeleveld et al. (1997a,b) have tried to investigate the relationship between the monopolar and bipolar MUAP recorded signal parameters and the MU location. They proposed a parametric equation to link the MUAP magnitude and the radial position of the recording electrode based on experimental results. A simulation study was done by Zhou and Rymer (2003) to estimate the number of MUAP in the surface EMG using a range of parameters. Researchers showed that some parameters are not suitable because they are not sensitive to the number of MUAPs. The above-mentioned studies reveal the difficulty to find suitable parameters since most of them are affected by several muscular characteristics at once. Therefore, it is very difficult to find a function which links the MU and volume conductor characteristics to the MUAP parameters. The relationship between these characteristics and MUAP parameters are non-linear because of the non-linearities of the electrical properties of the multiple elements that intervene between the muscle and the recording electrode. The multidimensionality of the studied problem is the main subject behind this difficulty. In fact, different values of a given characteristic may provide the same value of a given MUAP signal parameter, which rises up the so-called compensation problem. As an example, two motor units with different depths and sizes, a superficial one with a small size and a deep one with a big size, can both provide same value of an amplitude MUAP parameter. In the other hand, the same value of a given characteristic can provide different values of a given MUAP sig-

Fig. 2. The relationships between indicated characteristics in mm (h: depth of the MU; tf: thickness of fat layer, ts: thickness of skin layer, wi: width of innervation zone) and the parameters of the MUAP signal. (a) Peak to peak PP value. (b) Root mean square value RMS. (c) Median frequency MDF. (d) Rectified area SR.

328

A. Reffad et al. / Journal of Neuroscience Methods 164 (2007) 325–338

wi, the thickness of fat layer tf, and the thickness of skin layer ts. Each of these parameters decreases with the increases of the characteristic. If we take the peak to peak (PP) value as a discriminate parameter of the MU depth h and try to establish the linked function between them, we can suggest h = A/PP; A is a constant to be determined. But PP value also depends on the other characteristics, i.e. it is the result of the contribution of several characteristics and we are confronted to the compensation problem. In this case, the peak to peak value is viewed as an unknown or not discriminated parameter for the MUAP signal to extract the MU depth. Neural networks can be then an alternative approach to provide a solution to the mentioned problems. Bekka et al. (2002) have proposed an identification system based on neural networks to extract some muscle properties from the single differential MUAP. They have been confronted to the compensation problem due to the nature of the studied problem and the choice of the used parameters. Our analysis of the compensation problem leads us to assume that the solution of this problem when using the neural networks as a mathematical model depends essentially on the choice of the set of parameters which might be adequate to provide a better estimation. We have to find among the suggested parameters the sensitive ones and eliminate those which are redundant keeping a best representative vector parameter. 3.2. Parameters definition Additionally to the RMS and MDF parameters, several new parameters have been defined to enrich the set of double differential MUAP signal parameters; these parameters are defined in Table 2 and illustrated in Fig. 3. Sorting of these parameters permits to select the adequate parameters that check a certain criteria. Fig. 3a depicts the amplitude parameters. The surface and distance parameters are represented in Fig. 3b, where M1 and M2 are zeros crossing points limiting the MUAP signal that is going to be processed. Fig. 3c illustrates the curve parameters. Normally, the curvature degree at a point M of a curve represents the radius r of the tangent circle to the curve at this point. For the sake of simplicity, the calculation of curvature is represented by the calculation of the absolute X2 coefficient in the interpolated second-order polynomial function around this point. In our case, we have considered the interpolation around the peak for each lobe. Of course, these parameters will neither be used randomly nor completely. In this section, we have defined twenty three parameters, which depend mainly on the double differential MUAP signal morphology in order to increase the possibility to find good parameters. The next section reveals the new proposed approach which allows the selection of these parameters. 3.3. Parameters selection 3.3.1. Individual behavior of the parameters A sensitive parameter is the one that varies with the variation of each of the characteristics of the MU and volume conductor

Fig. 3. Double differential MUAP parameters. (a) Amplitude parameters. (b) Surface and distance parameters. (c) Curve parameters.

A. Reffad et al. / Journal of Neuroscience Methods 164 (2007) 325–338

329

Table 2 Parameters definition Parameters

Description

Amplitude parameters Fig. 3a P Amplitude of the principal peak P1 Amplitude of the first peak P2 Amplitude of the second peak PP1 Peak to peak value between the first and the principal peak PP2 Peak to peak value between the second and the principal peak DDP Absolute value of the difference between P2 and P1 Rap1 First amplitude ratio P1/P Rap2 Second amplitude ratio P2/P Surface parameters Fig. 3b S1 Area under the first lobe S2 Area under the second lobe SP Area over the principal lobe SR Rectified area between the points M1 and M2 Raps1 Area ratio S1/SP Curve, length and distance parameters Fig. 3b and c CR1 First lobe curve degree CR2 Second lobe curve degree CRp Principal lobe curve degree L Curve length signal between the points M1 and M2 Lp Curve length of the principal lobe d1 Distance between extreme points of the first and the principal lobes d2 Distance between extreme points of the second and the principal lobes Dac Sum of the distances d1 and d2 Other parameters DDD RLpD RMS MDF

Duration signal between points M1 and M2 Length-duration ratio Lp/DDD Root mean square Median frequency

whereas a non-sensitive is the one whose value remains constant. We have examined the individual behavior of the suggested parameters according to each characteristic separately. Thus, non-sensitive parameters have been eliminated. Fig. 4 shows this individual behavior for some sensitive parameters. Fig. 5 shows individual behavior of parameters that do not satisfy the sensitivity criteria. The RLpD is sensitive to the MU depth h (Fig. 5a), but it is neither sensitive to the thickness of skin layer ts nor to the innervation zone width wi. RLpD remains constant when ts and wi change from 2.5 to 3.5 mm and from 6 to 8 mm, respectively (Fig. 5c and d). Parameters DDD, Lp and L are also not sensitive. They remain constant (Fig. 5b–d). Therefore, parameters DDD, Lp, L and RLpD will be omitted. The rest of suggested parameters are sensitive and have not been represented because they have similar individual behaviors, therefore, we consider them as redundant parameters. After this pre-selection, we obtain the following sensitive parameters: P, P1, P2, PP1, PP2, S1, S2, SP, SR, D1, D2, DAC, DDP, CRp, CR1, CR2,Rap1, Rap2, Raps1, RMS and MDF. For the reason that the redundant parameters provide the same information, it is necessary to choose the most representative ones. In the section below, the global behaviors of

Fig. 4. Individual behavior of sensitive parameters normalized with respect to the corresponding first value according to the (a) MU depth h. (b) The thickness of fat layer tf. (c) The thickness of skin layer ts. (d) The innervation zone width wi.

the sensitive parameters is examined to ovoid the redundancy effect. 3.3.2. Global behavior of the parameters It was possible to eliminate redundant parameters on the basis of their individual behavior. However, this elimination

330

A. Reffad et al. / Journal of Neuroscience Methods 164 (2007) 325–338

We envisage identifying the characteristics whose values vary from 1 to 11 mm for the MU depth and the innervation zone width and from 1 to 5 mm for the thickness of fat and skin layers. We have sampled the characteristic space by step of 2.5 mm for the MU depth and the innervation zone width and step of 1 mm for the thickness of fat and skin layers. Therefore, we have obtained 625 combinations of the considered characteristics. For each combination, a MUAP signal has been simulated and parameterized with sensitive parameters. Series of a given parameter values according to each signal constitute a vector parameter which represents its global behavior. Fig. 6 shows the global behavior of the indicated parameters according to the MU and volume conductor characteristics combination. This figure depicts two examples of redundant parameters (PP2, DAC and S1) (Fig. 6a), (CR1, CR2 and CRp) (Fig. 6b) and no redundant ones (Rap1, Rap2, Raps1, and DDp) (Fig. 6c). To overcome the redundancy effect, we first proposed a procedure of non-supervised classification of the sensitive parameters. The procedure is based on the similarity, between the sensitive parameters, that use the correlation coefficient (CC). The principle of this procedure is to regroup the parameters that have the highest correlation coefficient. It is described by the following steps: (1) Take an arbitrary parameter from the set of parameters to be classified. (Each parameter is a 625 data vector.) (2) Calculate the CC between that parameter and all the rest of parameters (taken one by one). (3) Consider that parameters having a CC of 0.99 with the first parameter as a class. (4) Actualize the set of parameters omitting the parameters of the new defined class. (5) Repeat phases 1 to 4 till the classification of all parameters. After the classification phase, we have found the following classes: • • • • • • • • Fig. 5. Individual behavior of non-sensitive parameters normalized with respect to the corresponding first value according to the (a) MU depth h. (b) The thickness of fat layer tf. (c) The thickness of skin layer ts. (d) The innervation zone width wi.

is restricted (the parameters are computed according to one characteristic while the other ones are fixed). Thus, to complete the pre-selection made previously, we should examine the global behavior of each parameter according to a specific four dimension space defined by the MU and volume conductor characteristics combination (h, tf, ts, wi).

Class 1 = {P,P1,P2,PP1,PP2,S1,S2,SP,SR,D1,D2,DAC} Class 2 = {DDP} Class 3 = {Rap1} Class 4 = {Rap2} Class 5 = {Raps1} Class 6 = {CRp, CR1, CR2} Class 7 = {RMS} Class 8 = {MDF}

We notice that the classes 1 and 6 contain several elements which represent the same information described by the same manner; these elements are therefore redundant. The elimination of these redundant parameters leads to choose a representative parameter of each class. We have stated that the representative parameter of each class is the most homogeneous one. A parameter is mathematically considered as most homogeneous when it possesses the smallest coefficient of variation Cv (Baillargeon, 1989) corresponding to the entire set of redundant parameters.

A. Reffad et al. / Journal of Neuroscience Methods 164 (2007) 325–338

331

The expression of the coefficient of variation is given by the following expression: σ Cv = μ where σ is the standard deviation and μ is the arithmetic mean. We have arranged parameters in a descending order of homogeneity. The redundant parameters of the two classes are arranged as follows: • Class 1 = {P1,P2,PP1,D1,DAC,PP2,D2,P,S2,S1,SR,SP} • Class 6 = {CR1, CRP, CR2} We have chosen P1 from class 1 and CR1 from class 6, because they are the most homogeneous. Obviously there is no selection problem regarding the mono-parameter classes. Consequently, they will be directly selected. Finally the set of selected parameters are: P1, DDP, Rap1, Rap2, Raps1, CR1, MDF and RMS. 4. MU and volume conductor characteristics identification

Fig. 6. Global behavior of normalized parameters with respect to the corresponding first value according to the characteristic combination. (a and b) Examples of redundant parameters. (c) No redundant parameters.

It is believed that artificial neural networks (ANN) are appropriate to analyse natural systems like spinal networks and motor units potentials because of the non-linear relationship which characterize biological phenomena. Feed forward networks like the multilayer perceptron (MLP) are examples of supervised networks which easily model nonlinear systems and possess a high degree of generalization. A typical feed forward network consists of an input layer, a hidden neuron layer and an output layer of neurons. Input layer simply transmits inputs through weighted links to hidden neurons where weighted inputs are accumulated and processed by a transfer function to generate an output to be sent to the output layer. A similar process takes place in the neurons in the output layer where outputs are generated (Luo and Unbehauen, 1998). The MLP is considered as a powerful tool that has been used extensively for classification, non-linear regression, speech recognition, and many other applications (Schulte et al., 2005; Feldkamp et al., 1998; Kelly et al., 1990). The MLP uses the backpropagation algorithm which reduces the mean squared error between the actual output and the desired one for a set of input/output couples. This error evaluates the MLP network performance (Willis et al., 1991). We have used the MLP to investigate the relationship between the selected MUAP signal parameters (inputs) and the characteristics mentioned previously (outputs). Both identification systems described hereafter, use the neural networks approach. The identification module presented in Fig. 1 is substituted by a neural network module. Unfortunately, there is no theory for searching the optimum network configuration; therefore all architectures of the neural networks presented below were determined from a series of repeated experiments by a trial and error principle taking the best configurations according to the resulting performance. The limits of the identification systems have been analysed via a new test by varying model parameters, such as the electrode

332

A. Reffad et al. / Journal of Neuroscience Methods 164 (2007) 325–338

Fig. 7. Neural network module of the global identification system.

radius, the bar electrode dimensions and the orientation angles of the detection system that were fixed during the training process. 4.1. Global identification system GIS Once the parameter set is chosen, this module estimates the characteristics at the same time. Fig. 7 shows the neural network module of the global identification system. The inputs of the neural network are the eight double differential MUAP signal parameters selected previously. In order for the net to learn, one needs to present a number of training examples (the training set). In our case the training set comprises 625 different pairs of the set of selected parameters and their corresponding characteristic combinations. Several configurations were tested until an acceptable performance was obtained. For each experiment, ten runs of the network were used; each run had different weight initialization. The configuration was varied by changing the number of neurons in the first and the second hidden layers. The number of epochs denotes the number of times that the whole training set was presented at random at the input of the module during the training phase. The networks are trained on an average training value of 1000. A neural network configuration of (8:15:15:4), i.e. two hidden layers with 15 neurons each, reached its peak performance of 1.4 × 10−4 after 2500 epochs of training. The network suffered from a decrease in the performance when more hidden neurons have been added. This result confirms the theory that having too many or too few neurons in a hidden layer can have a negative effect on the network performance (Karayiannis and Venetsanopoulos, 1993). Fig. 8 shows the desired outputs (targets) with solid line and the learned outputs in doted line of the MU depth h, the innervation zone width wi, the thickness of fat layer tf and the thickness of skin layer ts according to the learned examples. For reason of clarity, we have presented the first 100 and 300 examples for the innervation zone width and the skin layer, respectively. The main learning problem in the global identification system is the compensation effect of the characteristic, i.e. for different characteristic combinations we have the same value of a given parameter. Because of the global nature of identification, each characteristic requires a particular adaptation of network weights so that it will be better identified. These particular adaptations of the global network cannot be exactly satisfied for all charac-

Fig. 8. The targets (—) and their corresponding learned outputs (· · ·) of the learning sets, of the characteristics after the learning process for the global estimation. (a) The MU depth. (b) The thickness of fat layer. (c) The thickness of skin layer. (d) The innervation zone width.

A. Reffad et al. / Journal of Neuroscience Methods 164 (2007) 325–338

333

Fig. 9. Neural network modules of the separate identification. (a) System for the MU depth. (b) System for the thickness fat layer. (c) System for the thickness skin layer. (d) System for the innervation zone width.

teristics. Therefore, all outputs of the neural network cannot be estimated with the same performance. These results have led us to think about separate identification system in which each characteristic is identified independently of the other characteristics to improve the performance. 4.2. Separate identification systems SIS We have proposed four separate identification modules to estimate separately the mentioned characteristics. By the trial and error principle, we have minimized the number of the net inputs for each module. The object of that minimization is to make the identification task easier by finding the parsimonious parameters. Fig. 9a shows the neural network module which identifies the MU depth. We note that this module has a configuration of (4:8:8:1) and only four input parameters P1, MDF, DDP and RMS. After 5000 epochs of training, the module has reached a performance of 6.78 × 10−11 , the targets and outputs are entirely confused (Fig. 10a). The neural network module of the separate system identifying the thickness of fat layer tf is depicted in Fig. 9b. It has a configuration of (7:14:14:1) with seven input parameters. The obtained performance after 12000 epochs of learning is 7 × 10−5 . Fig. 9c shows the module identifying the thickness of skin layer. It has five input parameters MDF, DDP, RMS, Rap1 and CR1, with a configuration of (5:13:14:1). A performance of 1.64 × 10−5 is reached after a training value of 10000 epochs. The module identifying the innervation zone width (Fig. 9d) was trained with (6:15:16:1) as configuration with six input parameters. After 12000 epochs a performance of 6.3 × 10−5 was obtained. Fig. 10b–d shows targets and learned outputs which are almost confused. Several series of trials were executed to find the parsimonious parameters for each identification module. We

Fig. 10. The targets (—) and their corresponding learned outputs (· · ·) of the learning sets, of the characteristics after the learning process for the separate estimation. (a) The MU depth. (b) The thickness of fat layer. (c) The thickness of skin layer. (d) The innervation zone width.

334

A. Reffad et al. / Journal of Neuroscience Methods 164 (2007) 325–338

Table 3 Estimation of the characteristics with the GIS and the SIS systems Characteristics

Characteristics identified globally (GIS)

Characteristics identified separately (SIS)

h (mm)

tf (mm)

ts (mm)

wi (mm)

hg (mm)

tfg (mm)

tsg (mm)

wig (mm)

hs (mm)

tfs (mm)

tss (mm)

wis (mm)

1 1 3.7 3.5 4.0 4.0 5 5 6.2 6.5 7.4 7.0 8.0 8.0 9.0 9.0

3 4.5 2.2 4.2 2.5 2.5 2.5 4.5 3.2 1.5 3.5 2.5 1.5 1.5 1.5 1.5

3 3.25 2.2 5 3.5 2.5 3.5 4.5 3.1 2.5 3.5 3.5 1.5 4.5 1.5 4.5

5.5 8.5 4.3 11 4 4 10 10 6.8 5.5 7.4 9 7 7 4 10

1.00 0.49 3.47 3.48 3.48 3.48 5.98 6.00 5.99 5.96 6.04 6.16 8.45 8.50 8.26 8.49

2.90 5.23 2.75 4.15 3.65 3.60 1.15 2.49 3.68 2.13 5.64 4.16 0.70 1.03 2.98 2.40

3.27 3.14 2.24 4.74 3.16 2.23 3.79 5.26 3.01 2.28 3.18 3.27 1.54 4.27 1.39 3.85

5.95 9.33 4.98 11 5.80 5.60 9.21 7.87 7.42 7.13 9.56 10.08 5.53 5.17 7.20 11

1 1 3.5 3.5 3.5 3.5 6 6 6 6 8.5 6 7.88 8.49 8.5 8.5

3 4.95 2.50 4.26 3.55 3.52 1.18 2.15 3.63 2.11 3.30 2.04 1.09 1.26 1.48 1.57

3.18 2.91 2.07 5 3.19 2.10 3.71 5 3.07 2.44 3.18 3.36 1.43 4.24 1.64 4.12

5.46 9.15 4.45 10.99 5.36 5.44 10.24 10.03 6.74 7.13 7.79 8.72 7.74 6.97 6.66 10.99

noticed that the MU depth estimation, the identification system converges very quickly with the four input parameters. It appears that the parameters MDF, RMS and DDP are solicited by all the separate identification modules. Therefore, they can be considered as the most significant parameters reflecting MU geometry.

Table 4 Average relative error between the estimated characteristics and the real characteristics values

5. Results and performance evaluation

Average relative errors of these estimated characteristics given in Table 3 are summarized in Table 4. From these results, we can consider that separate identified values are more accurate than global identified ones. To test the noise sensitivity performance of the identification systems, a modeled MUAP signal with the characteristics (h, tf, ts, wi) equal to (3.5, 3, 3, 8.5 mm) assumed to be an original signal is corrupted with additive Gaussian noise for different values of signal to noise ratio (SNR). The values of the SNR vary from 10 to 40 dB with a step of 10 dB. The reconstructed MUAP signals (the signals modeled with the identified characteristics) by the GIS and SIS systems are depicted in Fig. 11a and b, respectively. It is clear from Fig. 11b that there is a good fit between the MUAP simulated by the fixed model parameters (original signal) and those reconstructed from sep-

In order to evaluate the neural identification system performances, two kinds of tests are proposed in this study. The first is to test the generalization ability of the neural networks (common test). These networks should provide reasonable responses to inputs not belonging to the training set (Karayiannis and Venetsanopoulos, 1993). The second one is to test the noise sensitivity performance. With regard to the first one, several non-learned examples have been simulated and their identified values are obtained by the two kinds of systems (GIS and SIS). Table 3 reports non-learned examples of MU and volume conductor characteristics and their identified values by the global and separate systems.

GIS (%) SIS (%)

h

tf

ts

wi

11.98 8.074

41.51 21.30

7.57 6.70

22.83 13.22

Table 5 GIS and SIS identification of the characteristics when radius R changes from 1 to 5 mm by step of 0.5 mm, the signal is modeled with the characteristic combination (h, tf, ts, wi) equals to (3.5, 3, 3, 8.5 mm) R

1 1.5 2 2.5 3 3.5 4 4.5 5

Characteristics identified globally (GIS)

Characteristics identified separately (SIS)

h (mm)

tf (mm)

ts (mm)

wi (mm)

h (mm)

tf (mm)

ts (mm)

wi (mm)

3.48 3.48 3.49 3.49 3.5 3.5 3.51 3.51 3.51

3.11 3.16 3.24 3.33 3.44 3.55 3.66 3.82 3.98

2.98 2.94 2.88 2.85 2.8 2.77 2.71 2.7 2.59

8.4 8.41 8.41 8.45 8.5 8.64 8.89 9.37 10.3

3.5 3.5 3.5 3.5 3.5 3.5 3.5 3.5 3.5

2.96 3.02 3.12 3.23 3.28 3.34 3.43 3.56 3.76

3 2.97 2.9 2.84 2.79 2.79 2.74 2.64 2.57

8.48 8.72 9.08 9.87 10.07 10.83 10.99 10.99 11

A. Reffad et al. / Journal of Neuroscience Methods 164 (2007) 325–338

335

and SIS systems, a new test is introduced. Instead of identifying the characteristics with non-learned examples and fixing rest of model parameters like it is the case in the common test, the evaluation of the performance is accomplished indeed by changing one or several model parameters. In our study, we have tested the two identification systems regarding the change of the shape and

Fig. 11. The reconstructed MUAP signals from original signal corrupted by additive Gaussian noise with different values of SNR (Signals are normalized with respect to the peak value of original signal). (a) MUAPs reconstructed by the global system. (b) MUAPs reconstructed by the separate system. (c) Curves of the correlation coefficient CC between the original signal and the MUAPs reconstructed by the GIS and SIS systems.

arately identified characteristics when SNR takes the values: 40, 30 and 20 dB. However for SNR = 10 dB, there is no fitting between the original signal and the reconstructed one. It can be observed from Fig. 11a that global estimation is very affected for SNR = 10 and 20 dB. In these cases, there is no fitting between the original signal and those reconstructed from identified characteristics. It is also, interesting to observe that the estimate of the correlation coefficient between the original signal and those simulated from identified characteristics states that the noise sensitivity performance of the separate identification system is better than the global one (Fig. 11c). In order to show the credibility of GIS

Fig. 12. Relative error of the identified characteristics as a function of electrodes radius for the GIS and SIS systems. (a) The MU depth. (b) The thickness of fat layer. (c) The thickness of skin layer. (d) The innervation zone width.

336

A. Reffad et al. / Journal of Neuroscience Methods 164 (2007) 325–338

Fig. 13. Relative error estimation of the identified characteristics obtained as a function of the dimensions of bar electrodes (–, GIS), (䊉, SIS).

In this case, we have simulated different signals with different values of radius R. The neural network module for both systems has trained the examples at fixed electrodes radius of 1 mm. Table 5 summarizes the sensitivity of the GIS and SIS systems when the radius R takes values from 1 to 5 mm by step of 0.5 mm. Fig. 12 reports the relative error Err of the identified characteristics as a function of the electrode radius. When the electrode radius increases, we observe that the SIS system generally provides a better performance than the GIS system. Besides, the Err of the identified MU depth h is zero for the SIS system and remains less than 0.6% for the GIS system (Fig. 12a). Fig. 12b and c shows that the separate identification of layers thickness fat and skin) is relatively more accurate. However, in Fig. 12d, we notice that the results given by the GIS are better concerning the innervation zone width identification. We have also noticed that an increase of the electrode radius leads to a decrease of the identification performance.

separate identification of the MU depth h is independent of the electrode dimensions and the evaluated relative error expresses the accuracy of the identification quality (Fig. 13a). Nevertheless, the global system provides this estimation with Err smaller than 0.5%. For the fat thickness estimation (Fig. 13b), we observe that the performance of the GIS and SIS systems decreases when the bar electrodes dimensions increase and they affect similarly this performance. Nevertheless a better performance is obtained when a SIS system is used. The relative error Err computed in the case of the skin thickness depends on the variation of the longitudinal (along Z) and transversal (along θ) dimensions of bar electrodes. Indeed, it can be observed from Fig. 13c that the longitudinal bar electrodes give an improved performance in the case of the global identification. On the other hand, the separate identification performance is better when the transversal bar electrodes are used. Fig. 13d shows that the global identification system performance tends to decrease when the longitudinal dimension increases, i.e. this identification is not affected by the variation of the transversal dimension and gives better results when the longitudinal dimension electrodes is less than 5 mm (Err equals zero%). It should be noted that the GIS system for innervation zone identification remains more efficient than the SIS one.

5.2. Bar electrode change

5.3. Inclination angle change

The neural networks of the two identification systems described above have learned with circular electrodes with a constant value of radius. In this test, we will observe the effect of size variation of the rectangular electrodes on the systems performance. Fig. 13 depicts the relative error Err of the characteristics identified by the GIS and SIS systems as a function of electrode sizes in the two spatial directions θ and Z. The

To investigate the effect of the inclination angle for the LDD detection system on the characteristics estimation, several MUAPs have been simulated with different inclination angles values. Fig. 14a and b shows the relative error Err of the identified characteristics obtained by GIS and SIS systems, respectively, as a function of the inclination angle alpha. Examination of the obtained curves reveals that both systems provide a negligible

the size of the electrodes and the inclination angle of the detection system. In all operation tests, we have used the modeled signal with the characteristic combination (h, tf, ts, wi) equals to (3.5, 3, 3, 8.5 mm). 5.1. Electrode radius change

A. Reffad et al. / Journal of Neuroscience Methods 164 (2007) 325–338

Fig. 14. Relative error Err estimation of the identified characteristics obtained as a function of the inclination angle. (a) Estimations given by the GIS system. (b) Estimations given by the SIS system.

Err when the MU depth is estimated (0% for the SIS estimation and lower than 1% for the GIS estimation). The other characteristics are better identified by the GIS system. We note that both systems give reasonable results when the inclination angle remains lower than 5◦ . 6. Conclusions In this work, we have focused our study on the extraction of two MU characteristics and two volume conductor characteristics. The extraction of these characteristics from the MUAP signal parameters is not an easy task. In fact, we were confronted to the compensation and multi-solution problems. Therefore, the identification cannot be accomplished by a unique parameter, that does not mean that several parameters can solve the problem. Then, the research of optimal parameters returns to look for a compromise between the nature of parameters and their number. Thus, the study of the parameters behavior is an unavoidable phase to select the suitable parameters. Our study shows that the use of the neuronal network approach is an efficient alternative to reveal the relationship between the MU and volume conductor characteristics and the MUAP parameters. The optimal parameters selection takes place in two steps. The first step consists to eliminate non-sensitive parameters obtaining thus a restricted set of sensitive parameters. During the second step, we choose the most representative parameter among the previous whole set of the sensitive but still redundant parameters. The selection of the representative parameter, during the second step, is based on the similarity and homogeneity

337

criterion. The similarity criterion permits us to overcome the redundancy problem using a procedure, proposed in this paper, which classifies the parameters according to their similarity degree indicated by the correlation coefficient. The homogeneity criterion allows us to choose a representative parameter among the set of redundant parameters. The parameters selection procedure based on the above criterions is a contribution approach used to select the most suitable parameters of the simulated signal among the amount of available ones. As a consequence, the neural networks trained with this selected parameters provided reasonable results. To avoid the compensation problem accompanying the GIS training, we have thought to extract the characteristics separately, therefore the SIS was proposed. In general, we consider that SIS is more accurate than the GIS; this result was expected and has been confirmed by the common test. But it remains to improve the SIS by searching other suitable parameters, like it was done for the MU depth where the error was zero. The performance of the two systems was tested using three kinds of tests, the common test, the noise test and the change test. The last test has allowed observing the limits of the used systems when the model parameters fixed during the training process have been changed. The test of the global and separate neuronal network systems with additive Gaussian noise has shown that the identification process performs well when the signal to noise ratio is greater than 10 dB in the case of the separate estimation and greater than 20 dB in the case of the global estimation. The identification performance was also illustrated by computing the correlation coefficient between the original MUAP signal and reconstructed MUAPs signals. The electrode radius effect on the identification systems was studied. For all characteristics, except MU depth, the performance of both systems decreases when the electrode radius increases. The identification of the MU depth using the SIS system was very accurate with zero error whatever the variation of the electrode radius. When GIS is used, this identification remains very acceptable with a relative error lower than 0.6%. In general, the performance for both GIS and SIS systems is very altered if the longitudinal electrodes are used and it appears that this alteration is more important at longitudinal dimensions greater than 5 mm. Concerning the inclination angle test, both systems provide an acceptable estimation of the characteristics when the detection system is deviated from angles lower than 5◦ . The method of the optimal parameters selection, presented in this work, simplifies partially the identification problem while avoiding the use of several defined parameters involving no additional information for a better training of the neural network systems. Acknowledgements Authors sincerely thank Prof. R. Merletti and Prof. D. Farina for having welcomed Miss A. Reffad in the Laboratory for Engineering of the Neuromuscular System (Politecnico di Torino, Italy) in 2003 for 1 month practicum. We are grateful for the opportunities provided to her in this laboratory and for the pos-

338

A. Reffad et al. / Journal of Neuroscience Methods 164 (2007) 325–338

itive criticism provided by Prof. R. Merletti in reviewing the manuscript. References Basmajian JV, De Luca CJ. Muscle alive. 5th ed. Williams and Wilkins; 1985. Baillargeon G. Probabilit´es statistique et techniques de r´egression. Les editions SMG; 1989. Bekka RE, Boudaoud S, Chicouche D. The use of a neural network system in the identification of motor unit characteristics from surface detected action potentials: a simulation study. J Neurosci Methods 2002;116:89–98. Disselhorst-Klug C, Silny J, Rau G. Estimation of the relationship between the non-invasively detected activity of single motor units and their characteristic pathological changes by modelling. J Electromyogr Kinesiol 1998;8:323–35. De Luca CJ. Physiology and mathematics of myoelectric signals. IEEE Trans Biomed Eng 1979;26(6):313–25. Englehart K, Hudgins B. A robust, real-time control scheme for multifunction myoelectric control. IEEE Trans Biomed Eng 2003;50(7):848–54. Farina D, Mesin L, Martina S, Merletti R. A surface EMG generation model with multilayer cylindrical description of the volume conductor. IEEE Trans Biomed Eng 2004a;53(3):415–26. Farina D, Merletti R, Enoka RM. The extraction of neural strategies from the surface EMG. J Appl Physiol 2004b;96:1486–95. Farina D, Cescon C, Merletti R. Influence of anatomical, physical, and detectionsystem parameters on surface EMG. Biol Cybern 2002;86:445–56. Farina D, Fortunato E, Merletti R. Effect of electrode shape on spectral features of surface detected motor unit action potentials. Acta Physiol Phamacol Bulg 2001a;21:487–96. Farina D, Merletti R, Nazzaro M, Caruso I. Effect of joint angle on EMG variables in leg and thigh muscles. IEEE Eng Med Biol Mag 2001b;20:62–71. Feldkamp LA, Puskorius GV, Haykins S, Kosko B. A signal processing framework based on dynamic neural networks with application to problems in adaptation, filtering and classification. Proc IEEE 1998;86(11):2259– 77. Gydikov A, Gerilovsky L, Radicheva N, Trayanova N. Influence of the muscle fiber end geometry on the extracellular potentials. Biol Cybern 1986;54:1–8. Helal JN, Bouissou P. The spatial integration effect of surface electrode detecting myoelectric signal. IEEE Trans Biomed Eng 1992;39:1161–7. Holtermann A, Roeleveld K, Karlsson JS. Inhomogeneities in muscle activation reveal motor unit recruitment. J Electromyogr Kinesiol 2005;15:131–7. Hudgins B, Parker PA, Scott RN. Control of artificial limbs using myoelectric pattern recognition. Med Life Sci Eng 1994;13:21–38. Hogrel JV. Clinical applications of surface electromyography in neuromuscular disorders. Clin Neurophysiol 2005;35:59–71. Kelly MF, Parker PA, Scott RN. The application of neural networks to myoelectric signal analysis: a preliminary study. IEEE Trans Biomed Eng 1990;37(3):221–30. Kohn AF. Fractal number and spectral skewness: two features for the pattern classification of motor unit action potentials. 11th Annu Int IEEE Conf Eng Med Biol Mag 1989:1885–6. Karayiannis NB, Venetsanopoulos AN. Artificial neural networks: learning algorithms, performance evaluation and applications. Kluwer academic publishers; 1993. Lindstrom L, Magnusson R. Interpretation of myoelectric power spectra: a model and its applications. IEEE Proc 1977;65:653–62.

Luo FL, Unbehauen R. Applied neural networks for signal processing. Cambridge university press; 1998. Lateva ZC, McGill KC, Burgar CG. Anatomical and electrophysiological determinants of the human thenar compound muscle action potential. Muscle Nerve 1996;19:1457–68. Merletti R, Farina D, Gazzoni M. The linear electrode array: a useful tool with many applications. J Electromyogr Kinesiol 2003;13:37–47. Merletti R, Lo Conte L, Avignone E, Guglielminotti P. Modeling of surface myoelectric signals Part I: model implementation. IEEE Trans Biomed Eng 1999a;46(7):810–20. Merletti R, Roy SH, Kupa E, Roatta S, Granata A. Modeling of surface myoelectric signals. Part II. Model-based signal interpretation. IEEE Trans Biomed Eng 1999b;46(7):821–9. Merletti R, Parker P. Electromyography. In: John G, Webster R, editors. Wiley encyclopaedia of electrical and electronics engineering, vol. 6. New York: Wiley; 1998. p. 524–40. Merletti R, Knaflitz M, Deluca CJ. Electrically evoked myoelectric signals. CRC Crit Rev Biomed Eng 1992;19:293–340. Mesin L, Joubert M, Hanekom T, Merletti R, Farina D. A finite element model for describing the effect of muscle shortening on surface EMG. IEEE Trans Biomed Eng 2006;53(4):593–600. Parker P, Englehart K, Hudgins B. Myoelectric signal processing for control of powered limb prostheses. J Electromyogr Kinesiol 2006;16:541–8. Pattichis CS, Schizas CN, Middleton LT. Neural network models in EMG diagnosis. IEEE Trans Biomed Eng 1995;42(5):486–96. Rau G, Disselhorst-Klug C. Principles of high-spatial-resolution surface EMG (HSR-EMG): single motor unit detection and application in the diagnosis of neuromuscular disorders. J Electromyogr Kinesiol 1997;7(4): 233–9. Roeleveld K, Stegeman DF, Vingerhoets HM, Van Oosterom A. The motor unit potential distribution over the skin surface and its use in estimating the motor unit location. Acta Physiol Scand 1997a;161(4):465–72. Roeleveld K, Stegeman DF, Vingerhoets HM, Van Oosterom A. Motor unit potential contribution to surface electromyography. Acta Physiol Scand 1997b;160:175–83. Roeleveld K, Blok JH, Stegeman DF, Van Oosterom A. Volume conduction models for surface EMG: confrontation with measurements. J Electromyogr Kinesiol 1997c;7(4):221–32. Roeleveld K, Stegeman DF, Falck B, Stalberg E. Motor unit size estimation: confrontation of macro EMG and surface EMG. Electroenceph Clin Neurophysiol 1997d;105:181–8. Stegeman DF, Linssen W. Muscle fiber action potential changes and surface EMG: a simulation study. J Electromyogr Kinesiol 1992;2:130–40. Stulen FB, De Luca CJ. Frequency parameters of the myoelectric signal as a measure of muscle conduction velocity. IEEE Trans Biomed Eng 1981;28:515–23. Staudenmann D, Kingma I, Stegeman DF, van Die¨en JH. Towards optimal multichannel EMG electrode configurations in muscle force estimation: a high density EMG study. J Electromyogr Kinesiol 2005;15:1–11. Schulte E, Dimitrova NA, Dimitrov GV, Rau G, Disselhorst-Klug C. Estimation of the muscle fibre semi-length under varying joint positions on the basis of non-invasively extracted motor unit action potentials. J Electromyogr Kinesiol 2005;15:290–9. Willis MJ, Di Massimo C, Montague GA, Tam MT, Morris AJ. Artificial neural networks in process engineering. IEEE Proc D 1991;138(3):256–66. Zhou P, Rymer WZ. Estimation of the number of motor unit action potentials in the surface electromyogram. Proc IEEE EMBS Eng 2003:372–4.