Med. Biol. Eng. & Comput., 1977, 15, 648-655

Effect of motor-unit firing time statistics on e.m.g. spectra P. L a g o

N.B.

Jones

Graduate Division of Biomedical Engineering, Centre for Medical Research,Universityof Sussex ,Falmer, Brighton, Sussex BN1 9QT, England

Abstract--The generation of an e.m.g, can be modelled by linear filtering of N impulse trains, regarded as the neural inputs to N moto[ units, where the e.m.g, is the sum of the outputs of the N filters. Starting from this model, it is she wn that not only are the shapes of the N action potential images as viewed by the electrodes an important factor in defining the spectrum, but so is the interspike interval distribution of the impulse trains. Changes in the statistics of the pulse trains are shown to affect the low-frequency end of the power spectrum even i f the motor-unit action potentials do not change shape. The effect is illustrated by use of data already published on deltoid muscle and brachii biceps. A means of including motor-unit synchronisation is introduced and it is shown that its effect can be contrary to that due to the firing statistics. K e y w o r d s - - E . M . G . , Spectral analysis, Muscle, Spike trains

1 Introduction IN RECENT years, frequency analysis has been used as a means of data reduction and characterisation o f the e.m.g. Evidence has been found that modifications of the power spectral density of the e.m.g. signal are related to abnormal physiological states. This means that power-spectral-density analysis may eventually be used in the early detection of neuropathies or myopathies or as a means of observing the progress of a disease. Other modifications of power spectral density related to fatigue and level of contraction have been reported in the literature (KADEFORS et al., 1968; JOHANSSON et al., 1970; O'DONNELL et al., 1973; FAY et al., 1976). Although the comparison of the results is difficult mainly due to differences in the data reduction and experimental techniques, at least a clear pattern emerges: a shift of energies from the higher to the lower frequencies with fatigue. With respect to changes with load, however, no generally accepted pattern exists. F o r example, FAY et al. (1976), working on the pelvic floor muscles, reported a power increase in the band 0-50 Hz of up to one-third of the maximum contraction and a substantial increase in the band 0-100 Hz for two-thirds of the maximum contraction, with a corresponding decrease of energies in the high-frequency band. Using surface electrodes on the biceps, O'DONNELL et al. (1973) found the percentage of power at low frequencies (26-49 Hz) was higher for a moderate contraction than for a nearly maximum contraction and a continuous First received20th Septemberandin finalform 25th November1976

648

increase in the high-frequency band ( > 6 0 H z ) with increasing contraction. Because of the advantages which would result from a consistent theory to aid in interpreting e.m.g, spectra, it is essential to relate the elemental phenomena affecting e.m.g, generation to the observed spectral changes. So far, spectral modifications have been related to changes in the motor-unit action potential, fibre type mixture, conduction velocity, synchronisation phenomena, percentage of polyphasic action potentials and dropping-out of small motor units. (KADEFORS et al., 1968; LINDSTROM et al., 1970; LARSSON, 1969; PERSON and MISmN, 1964; O'DONNELL et al., 1973.) In this paper, another cause is established: the alteration of the statistical properties of the interspike interval, which have been found to depend both on fatigue and the level oI contraction (D~ LUCA and FOgREST, 1973). At the same time evidence is presented to show that ignoring these statistical properties affects the correct interpretation of the power spectral density, particularly in the lowfrequency band. 2 Models of the generation of myoelectric signals Let the muscle he constituted by a set of N motor units, each of them excited, if recruited, by a train of spikes, u~(t); i = 1, N. Each time a spike occurs, the associated motor unit contracts and an electrical signal is generated. Let h~(t); i = 1, N, be the observed motor-unit action potential (m.u.a.p.). The summing property of motor-unit action potentials has been established experimentally (Bin6 and PARTRIDGE, 1971). It follows, as shown in Fig. 1,

Medical & Biological Engineering & Computing

November 1977

that the summation of these motor-unit action potentials trains (m.u.a.p.t.) constitutes the observed myoelectric signal. For each spike train Hi(t) define a train of Dirac pulses,

~,(t)

Z ~(t- tki).

=

the following mathematical model: the time of occurrence of the spikes in the train u~(t) constitutes a sample from a stochastic point process of the renewal type. To characterise the generation mechanism of the train of pulses, it is sufficient to define the probability density of the interspike interval,

f(t).

k

The same m.u.a.p.t, will be generated if 5i(t) is the input signal to a linear dynamic system with impulse response hi(t). It follows that the schematic model considered before is equivalent to the linear dynamic system represented in Fig. 2. Now, if h~(t) values are assumed to be equal, hi(t) = h(t); i = 1, N, then the equivalent model represented in Fig. 3 is obtained (BRODY et al., 1974). It is a direct consequence of the linearity assumption.

3 Statistical properties of the spike train For constant-force isometric contractions, CLAMANN (1969) found the spike trains Hi(t); i = 1, N are not periodic, but random: the interspike interval behaves as a random variable and no evidence has been found of significant correlation between interspike intervals (CLAMANN,1969; DE LUCA and FORREST, 1973; SHIAVI and I'qEGIN, 1975a). These properties can be incorporated into

The simplicity of model-3 over model-2 is certainly a very strong argument for its use in the analysis and modelling of the e.m.g., but particular care should be taken when describing the pooled spike train. This is because the description of the statistical properties of the pooled train N

"("~(0

=

Z

i=l

u,O)

is not simple: the pooled train ut~>(t) can no longer be considered as generated by a stochastic point process of the renewal type, and all the mathematical simplicity is lost. Assuming the probability density of the interspike interval is the same for all motor units, it has been shown by Cox and SMITH 0954) that the probability density of the interspike interval of the pooled train is given by +oo

(1)

u,(t)

I

I I

'

tI

t K

k+ I

MUAPT$

t'

k~a

+

Fig.

1 Generation of the e.m.g.

.UN(t) Motor unit N

II

tU J

I

%

/

MUAPT N

A

_A

.

~j~J

9

~)i)

x(t)

Myoelectric Signol

Mod~l(2)

Fig. 2 Equivalent linear dynamical model

6=(P) [=r~

5L(t? Model (31

Fig. 3 Single-input linear dynamical model Medical & Biological Engineering & Computing

where F(t) is the probability distribution function of the interspike interval, /t the mean value of the interspike interval and n the number of iecruited motor units. It follows from eqn. 1 that ftp~(t) tends quite rapidly to the exponential distribution of mean #/n when n tends to infinity regardless of F(t). This property has been wrongly used (BRODY et al., 1974) to validate the Poisson model* for the pooled train for n sufficiently large, but in fact the stated property is just a necessary, in no way a sufficient, condition. The Poisson characteristics of * A renewal stochastic Poisson poi nt process

November 1977

649

the limiting process as n tends to infinity should be considered only as 'local' (Cox and SMITH, 1954). Only when it tends to infinity at the same time as n can the local Poisson characteristics of the pooled train then become global. This is the so-called Palm-Khintchine limit theorem (KmNTCHINE, 1960). In the particular situation considered, this theorem states that, under fairly general conditions, the limiting process is of the Poisson type if the following conditions are met: (i) n -+ o0

(ii) it ~ ~

but

where H(jw) is the Fourier transform of h(t) and Ooa(w) the power spectral density of 6f(t); i = 1, n, the following relation is obtained:

Oxx(W) = n. IHCjw)I 2.~0~(w) + 2re. n ( n - 1) . H2(0 ).fi(w)

It remains to relate tb~a(w) to the probability density of the interspike interval f(t). The autocorrelation function of

n/it = constant

dl(t); i = 1, n,

and obviously the second condition is not satisfied.

R~(r),

can be expressed as follows (ZADEH, 1957):

R~(r) =

4 P o w e r spectral analysis

For simplicity it is assumed that

{

(m:(r)+ m : ( - r)); r r 0

~

Let X(t) be the myoelectric signal, Rxx(r) its autocorrelation function, and Cxx(r) the covariance function. Let &(t) be the motor-unit action potentim train m.u.a.p.t., Without any loss in generality, it can be assumed the first n motor units are recruited. It follows from the relation

where m:(r) is the renewal density or intensity function of the renewal process as defined by Cox (1962). It follows again from the theory of stochastic point processes (Cox, 1962) that for a renewal process the Laplace transform of m:(r) and f(t), respectively, M:(s) and F(s) are related by

x(t) = ~ x,(t) i=1

M :(s)

F(s) 1 -

and the experimental evidence of no significant interaction between motor units (SmAw and NEG~N, 1975b), that

Cxx(r) = n. C:,x(r) . . . . . . .

(2)

where C~(T) is the covariance function of x~(t); i = 1, n. From eqn. 2 and the relationship E(5~(t)) = 1//2, it follows that

. . . . . . .

(6)

. . . . . . .

(7)

F(s)

and lim m:(r) = t+o~_

1 it

Using relationships 4, 5,6 and 7, it is found that

FOw) it lH(jw)12"

it~

.

h(u)du

then

1 + 1--F(jw)

F(-jw) t" + 1-F(-jw)J' w#O

Oxx(W) =

0

(8)

n2

-~- .2~.H~(o).a(w) It-

Oxx(W) = .n%.(w) on

+2re

n(n-1)

It2

9

{l

"hOd.du

~0

1

.6(w)

(3)

where ~xx(W) and ~xx(w) are the power spectra density of X(t) and xi(t); i = l , n , respectively" Finally as

9 xx(W) = H(jw).H(-jw).~na(w) -- ]n(jw)l=~o~(w) 650

(5)

2--d(r); r = 0 /2

hi(t) = h(t); i = 1, N.

Rxx(Z) = n. Rxx(r) +

(4)

its

; w = o

Let the normalised power spectral density Oxx(W), be defined by

F(jw) f(-jw) Oxx(W) = 1 + 1 - F(jw--~)+ 1 - F ( - j w ) ' w # 0

(9)

Oxx(W) describes the interrelation between the statistical properties of the spike train and the power

Medical & Biological Engineering & Computing

November 1977

spectral density of the e.m.g. The Poisson model for the pooled train of pulses (BRoDY et al., 1974) corresponds to the approximation

(c) Gaussian distribution (not truncated)

l

f ( t ) -- a/(2rc)a .exp

Oxx(W) = 1 for n sufficiently large but Oxx(W) does not depend on n, and so the approximation is not always justified. The power spectral analysis for a single train of delta functions has been developed before. According to COGGSHALL (1973), the power spectrum of a train of delta functions can be expressed as

Sx(w)=fl.

[ 1 + , ="~' 1 (F,(jw)+F,(-jw))]

(10)

where fl is the mean firing rate and F.(jw) the Fourier transform of f~(t), the probability density function associated with the sum of n consecutive intervals. It is possible to show that for a renewal process

~-tr2 ]

F(jw) = e - ~"W.e . . . . 2/2 I - - e - G ' 2 W2

O~(w) =

1 + e -~ w2- 2. cos(/tw) .e-~*2w~ (11)

Power spectral analysis of e.m.g, modelled by a renewal gaussian process with triangular motorunit action potentials has already been developed (LIBKIND, 1968; 1969). In this particular situation eqn. 8 gives

~xx(W)

n =

-

-

It

256. b 2 . sin4(Cw/4), sin2(Cw/2). (1 - e - ~ 2~,~) X

Sx(W) = 1 . 0 x x ( W )

( --

C z . w~ . (1 + e - ' 2 w2 _ 2 cos(pw), e-~o~ w2),

; w ~ 0

It

w~ 0

(12)

In general, the normalised power spectral density Oxx(W) as given by eqn. 9 has to be evaluated numerically using a n f.f.t, algorithm. Some important exceptions are:

where b a n d c characterise the motor-unit action potential signal, Fig. 4, as compared with Libkind's result.

(a) G a m m a distribution

~bxx(W) =

n. 256. b 2 . sin4( Cw/4) . sin2( Cw/2) C 2 . (1~ + 2 C ) . w 4

; w ~ 0

(13)

f(t) - B.F(p) " exp

Let tr tend to zero. i n the limit, the spike train will become periodic with period /z and ~ x x ( j w ) will -

;

If(t) = 0 ; F(jw) = e -jw". (1 + j w B ) -p

t/> a

t