- Email: [email protected]
An analytical model of the vestibular evoked myogenic potential
Journal of Theoretical Biology 286 (2011) 41–49
Contents lists available at ScienceDirect
Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi
An analytical model of the vestibular evoked myogenic potential n ¨ ¨ ¨ Bernd Lutkenh oner , Turker Basel
ENT Clinic, M¨ unster University Hospital, Kardinal-von-Galen-Ring 10, 48129 M¨ unster, Germany
a r t i c l e i n f o
abstract
Article history: Received 10 January 2011 Received in revised form 27 April 2011 Accepted 6 July 2011 Available online 14 July 2011
The vestibular evoked myogenic potential (VEMP) can be modeled (scaling factors aside) as a convolution of the motor unit action potential (MUAP) of a representative motor unit, h(t), with the temporal modulation of the MUAP rate of all contributing motor units, r(t). Accordingly, the variance modulation associated with the VEMP can be modeled as a convolution of r(t) with the square of h(t). To get a deeper theoretical understanding of the VEMP phenomenon, a specific realization of this general model is investigated here. Both r(t) and h(t) were derived from a Gaussian probability density function (in the latter case taking the first derivative). The resulting model turned out to be simple enough to be evaluated analytically in the time and in the frequency domain, while still being realistic enough to account for the basic aspects of the VEMP generation. Perhaps the most significant conclusion of this study is that, in the case of noisy data, it may be difficult to falsify the hypothesis of a rate modulation of infinitesimal duration. Thus, certain aspects of the data (particularly the peak amplitudes) can be interpreted using a short-modulation approximation rather than the general model. The importance of this realization arises from the fact that the approximation offers an exceptionally simple and convenient way for a model-based interpretation of experimental data, whereas any attempt to use the general model for that purpose would result in an ill-posed inverse problem that is far from easy to solve. & 2011 Elsevier Ltd. All rights reserved.
Keywords: VEMP modeling Canonical VEMP Variance modulation Motor unit action potential Electromyogram
1. Introduction Reflexes elicited by stimulating the otoliths in the inner ear with high-level sound pulses modulate the tone of various muscles. The accompanying short-term perturbation of the electromyogram (EMG) can be measured as the vestibular evoked myogenic potential (VEMP). The main component of the VEMP, classically recorded over the sternocleidomastoid muscles (Colebatch et al., 1994), arises from a brief inhibition of the muscle tone. Thus, a sine qua non for a successful measurement is a tonic muscle contraction. Recording of the VEMP has become a standard clinical test of otolith function, with documented diagnostic utility in varied vestibular and central nervous system disorders (Rosengren et al., 2010). VEMPs with reduced amplitude on the affected side were found, for example, in patients with vestibular schwannoma (Ushio et al., 2009) and unilateral Menie re’s disease (Kingma and Wit, 2011). Patients with Tullio phenomenon (hypersensitivity of the vestibular system to sound), by contrast, show an augmented VEMP amplitude on the affected side, associated with an abnormally low VEMP threshold (Colebatch et al., 1998; Brantberg and Verrecchia, 2009).
n
Corresponding author. Tel.: þ49 251 83 56864; fax: þ49 251 83 56882. ¨ ¨ E-mail address: [email protected] (B. Lutkenh oner).
0022-5193/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2011.07.002
Wit and Kingma (2006) developed a simple model for the generation of the VEMP. A key constituent of the model is the motor unit action potential (MUAP). All motor units are assumed to have the same MUAP, apart from a unit-dependent factor that is thought to account for varying distances between motor unit and recording electrodes. The EMG, being the algebraic sum of the MUAP trains of all motor units (Day and Hulliger, 2001), can be simulated by adding MUAPs with random occurrence times and amplitudes. As the MUAPs are zero-mean signals, their random occurrence causes partial cancellation (positive and negative waves tend to extinguish each other). What remains is stochastic activity with an expected value of zero. A simple modification transforms this EMG model into a model for the VEMP: Instead of being constant in time, the probability of generating a MUAP drops to zero for a short period, corresponding to complete inhibition of the MUAP generation. The inhomogeneity in the MUAP rate disturbs the partial cancellation of MUAPs so that the expected value of the EMG temporarily deviates from zero. This is what is known as the VEMP. Building upon the pioneering work of Wit and Kingma (2006), we developed a somewhat more general and, with respect to ¨ ¨ computation time, more efficient model (Lutkenh oner et al., 2010). Moreover, we showed that the model can be greatly simplified by assuming that the duration of the inhibition is short (as compared with the time constants that determine the MUAP).
42
B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 286 (2011) 41–49
A comparison between model and data basically confirmed the ¨ ¨ theoretical predictions (Lutkenh oner et al., 2011). In particular, we were able to verify that the sound-induced inhibition affects not only the mean of the EMG (causing the VEMP), but also the variance. Some findings were unanticipated, though. For the most part, the observed complications can be attributed to the VEMP being composed of multiple components. But there is another point deserving attention: the assumed short duration of the inhibition. In a study of the firing of single motor units, Colebatch and Rothwell (2004) found inhibitions with durations between 2 and 8 ms (mean 3.6 ms). This may be too long to be compatible with the key assumption underlying the short-inhibition approximation, namely that the MUAP function is roughly constant for the whole duration of the inhibition. The main purpose of the present study was to develop an understanding of how the duration of the inhibition impacts on the VEMP and the associated variance modulation. The model presented here was chosen so that it is, on the one hand, realistic enough to account for the basic aspects of the VEMP generation, and, on the other hand, simple enough to allow analytical evaluation.1 Special attention is given to the question whether the short-inhibition approximation is applicable even when the underlying assumption is not fully met. This question is important because the approximation offers an exceptionally simple and convenient way for a model-based interpretation of measured data. The general model, by contrast, leads to a non-linear deconvolution problem that is far from easy to solve. With deconvolution in mind, our analytical model is considered in both the time- and the frequency domain so that it may serve as a reference model also for future work on VEMP deconvolution.
2. General theory To begin with, we briefly recall the theory developed pre¨ ¨ viously (Lutkenh oner et al., 2010). The notations were partially changed, to better meet the requirements of the current article. Moreover, the theory was slightly generalized so that both inhibition and excitation can be considered. Besides, the model is also considered in the frequency domain. 2.1. General model in the time domain The model is based on the simplifying assumption that the MUAPs are identical for all motor units, apart from an amplitude factor that may differ from unit to unit. The MUAP time course is described by a dimensionless zero-mean function h(t). The normalization of h(t) is arbitrary. Here, we normalize so that the maximum absolute value is one. The quantity Z þ1 D¼ h2 ðtÞdt ð1Þ
2
s20 ¼ a r0 D
ð2Þ
Note that r0D is, loosely speaking, the average number of MUAPs significantly overlapping in time. If the rate exhibits a temporal modulation r(t) such that
rðtÞ ¼ r0 ð1 þrðtÞÞ
ð3Þ
then the mean of the EMG, i.e. the VEMP, is ~ vðtÞ ¼ ar0 DvðtÞ
ð4Þ
with ~ ¼ D1 vðtÞ
Z
þ1
rðt 0 Þhðtt 0 Þdt 0
ð5Þ
1
For the variance we get
s2 ðtÞ ¼ s20 ð1 þ mðtÞÞ
ð6Þ
with mðtÞ ¼ D1
Z
þ1
rðt 0 Þh2 ðtt 0 Þdt 0
ð7Þ
1
The functions r(t) and m(t) will be referred to as the rate modulation and the variance modulation, respectively. Thepstanffiffiffiffiffiffiffiffiffi dard deviation of the EMG is, according to Eq. (2), s0 ¼ a r0 D. Dividing the VEMP signal v(t) by s0 yields the normalized VEMP qffiffiffiffiffiffiffiffiffi ~ vnorm ðtÞ ¼ r0 Da=avðtÞ ð8Þ The normalized VEMP is a dimensionless measure, as is the variance modulation. Yet there is an important difference. Unlike the variance modulation, the normalized VEMP depends on the base rate r0, i.e. on the muscle tone. In a modeling study, this ~ complication can be avoided by considering vðtÞ instead of v(t) or vnorm(t). If needed, the latter two functions can be obtained by ~ simple rescaling. With regard to the fact that the function vðtÞ represents the VEMP reduced to the simplest form possible, it will be referred to as the canonical VEMP. A comparison of Eqs. (5) and (7) shows that canonical VEMP and variance modulation make a perfectly matched pair. The ratio a=a is generally close to 1, and a pffiffiffiffiffiffiffiffiffi r0 D would be about 4 (corresponding to typical value for D ¼0.01 s and r0 ¼1600 s 1). Thus, normalized and canonical VEMP may be expected to differ by a factor of the order of 4. 2.2. Approximation for short rate modulations
1
may then be interpreted as the effective MUAP duration, because it measures how long h(t) has an absolute value close to maximum.2 The MUAPs of individual motor units may be 1 More tedious calculations were assisted by the computational software Mathematica 7 (Wolfram Research, Inc.). 2 Eq. (1) is not meant as a definition of effective duration for arbitrary signals. However, for a temporally compact, normalized signal such as the MUAP function h(t) considered here, the integral on the right-hand side may be interpreted that way. For the MUAP function introduced in Eq. (17), calculation of an effective duration according to Gabor (1946) yields
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Z Z þ1
Dt ¼
obtained by multiplying h(t) by a unit-dependent amplitude factor. Detailed knowledge about the distribution of this factor is not required. Only the arithmetic mean, a, and the quadratic mean, a, are needed. The above assumptions allow us to integrate the MUAP trains of the individual motor units into a single train. In concrete terms, we assume that the motor units as a whole fire at a rate r(t). In the case of a constant rate, which means r(t) r0, the recorded EMG has zero mean and variance
þ1
t 2 h2 ðtÞdt
2p
1
h2 ðtÞdt ¼
1
pffiffiffiffiffiffi 3py
If the rate modulation r(t) corresponds to a pulse that is so short that h(t) remains approximately constant for the entire pulse duration, Eqs. (5) and (7) can be greatly simplified. The modulation effect turns out to be proportional to Z þ1 d¼ rðtÞdt ð9Þ 1
(footnote continued) The D resulting differ by pffiffiffiffiffiffi pffiffiffi from Eq. (1) is provided in Eq. (18). The two measures the factor 2 3=e 1:27. Regarding the fact that the factor 2p in Gabor’s definition is relatively arbitrary – other authors omit this factor (see, e.g., Eq. (53) of Boashash (1992)) – our usage of the term effective duration appears to be justified.
B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 286 (2011) 41–49
where 9d9 is the equivalent rectangular duration of the rate modulation. The sign of d indicates the type of modulation (positive sign for excitation and negative sign for inhibition). The equations for canonical VEMP and variance modulation reduce to ~ ¼ dD1 hðtÞ vðtÞ
ð10Þ
mðtÞ ¼ dD1 h2 ðtÞ
ð11Þ
43
0
While these two equations are strictly valid only for a rate modulation of infinitesimal duration, they may be expected to be at least approximately valid when the duration of the rate modulation is short compared to the time constants that determine the MUAP function h(t). This is why they will be referred to as the short-modulation approximation. What ‘short’ actually means and whether the approximation is applicable to typical experimental conditions will be examined in the model simulations below.
r(t) -1 1 h2(t)
0
2.3. General model in the frequency domain Canonical VEMP and variance modulation may be considered as filtered versions of the rate modulation r(t), because they are obtained by a convolution with h(t) and h2(t), respectively (cf. Eqs. (5) and (7)). This filtering point of view is best analyzed in the frequency domain. The Fourier transform of a function x(t) is Z 1 xðtÞexpð2piftÞdt ð12Þ Xðf Þ ¼
h(t) -1 0.5 ~ v(t)
1
~ (e.g. Brigham, 1974). If the Fourier transforms of vðtÞ, m(t), r(t), h(t), and h2(t) are denoted as V~ ðf Þ, M(f), R(f), H(f), and H2(f), respectively, the convolution theorem (see e.g. Brigham, 1974) allows us to rewrite Eqs. (5) and (7) as V~ ðf Þ ¼ D1 Rðf ÞHðf Þ
ð13Þ
Mðf Þ ¼ D1 Rðf ÞH2 ðf Þ
ð14Þ
0
m(t) -0.5 -10
3. Simple Gaussian model 3.1. Time domain To put the general theory outlined above into concrete terms, specific assumptions about the rate modulation function r(t) and the normalized MUAP function h(t) will be made now. The rate modulation function was derived from the probability density function of a Gaussian distribution t2 ð15Þ rðtÞ ¼ r0 exp 2 2t The parameter t controls the width of the function, 9r09 determines the maximum modulation effect, and the sign of r0 distinguishes between inhibition (negative sign) and excitation. Since complete inhibition (i.e. r(t)¼ 1) represents a limiting condition, the requirement r0 Z 1 has to be fulfilled. Solving the integral in Eq. (9) yields pffiffiffiffiffiffi d ¼ 2pr0 t ð16Þ Fig. 1a shows the graph of r(t) for r0 ¼ 1 and t ¼3 ms. For h(t) we choose the first derivative of a Gaussian, as suggested by Wit and Kingma (2006): ! t ðt=yÞ2 1 hðy; tÞ ¼ exp ð17Þ y 2
0 time (ms)
10
Fig. 1. Model in the time-domain. (a) Rate modulation function r(t). (b) Normalized MUAP function h(t) (black curve) and its square (gray curve). ~ (c) Canonical VEMP vðtÞ (black curve) and variance modulation m(t) (gray curve).
The black curve in Fig. 1b shows this function for y ¼4 ms. There is a minimum at y and a maximum at y (dotted vertical lines). The function is normalized so that the values at minimum and maximum are 1 and 1, respectively. The square of h(t) is displayed as a gray curve. Solving the integral in Eq. (1) yields
D¼
pffiffiffiffi e p y 2:41y 2
ð18Þ
Thus, in the present example, the MUAP function has an effective duration of about 9.64 ms. The canonical VEMP is obtained by solving the integral in Eq. (5) ~ ¼ vðtÞ
23=2 e1=2 ytr0 2
ðy þ t2 Þ3=2
t exp
t2 2
!
2ðy þ t2 Þ
ð19Þ
The black curve in Fig. 1c shows this function for r0 ¼ 1, t ¼3 ms ~ has and y ¼4 ms. A closer inspection of Eq. (19) reveals that vðtÞ basically the same mathematical structure as h(y;t). Indeed, a
44
B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 286 (2011) 41–49
simple transformation yields
0
~ ¼ v~ 0 hðyt ; tÞ vðtÞ
ð20Þ
a
with
yt ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi y2 þ t2
ð21Þ
R(f)
and 2 v~ 0 ¼ 23=2 e1 r0 yt=yt
ð22Þ
~ are at t^ v ¼ yt and t^ v (dashed vertical lines The extrema of vðtÞ in Fig. 1c). Eq. (22) may be rewritten as )1 pffiffiffi ( ðytÞ2 2 1þ r0 ð23Þ v~ 0 ¼ e 2yt
-0.008 0.01
v~ 0 ¼ d=D
b
H2(f)
which shows that the maximum condition for v~ 0 =r0 is y ¼ t (for any other combination of y and t, the term in curly brackets is greater than one). Thus, the global maximum of v~ 0 is pffiffiffi 2=e9r0 9 0:529r0 9. In the limit t-0, Eq. (22) may be rewritten as
0
ð24Þ
consistent with Eq. (10). The variance modulation is obtained by solving the integral in Eq. (7) ! t3 ty2 t2 t2 ð25Þ exp 2 þ mðtÞ ¼ 23=2 r0 2 2 3=2 5=2 ðy þ2t2 Þ ðy þ 2t2 Þ y þ 2t2
Im[H(f)] -0.01 0.008
0
The latter times are marked by short vertical lines in Fig. 1c (almost coinciding with the dotted vertical lines). For t Z y, by contrast, 9m(t)9 has is a single minimum at t¼0.
M(f) -0.008
3.2. Frequency domain
-200
The Fourier transform of r(t) is pffiffiffiffiffiffi Rðf Þ ¼ 2p t r0 expð2p2 t2 f 2 Þ
ð27Þ
and the Fourier transforms of h(t) and h2(t) are 2
c
~ Im[V(f)]
The function is displayed as a gray curve in Fig. 1c (same parameters as used for the canonical VEMP). For t o y, 9m(t)9 has a minimum at t ¼0 and maxima at 7 t^ m with qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi t^m ¼ y 1 þ ðt=yÞ2 2ðt=yÞ4 ð26Þ
2
Hðf Þ ¼ ið2pÞ3=2 e1=2 y f expð2p2 y f 2 Þ
ð28Þ
pffiffiffiffi 2 2 H2 ðf Þ ¼ e py=2ð12p2 y f 2 Þexpðp2 y f 2 Þ
ð29Þ
The Fourier transforms of canonical VEMP and variance modulation are easily obtained using Eqs. (13) and (14). Since the functions r(t), h2(t), and m(t) are even functions, their Fourier ~ transforms are real-valued. The functions h(t) and vðtÞ, by contrast, are odd so that their Fourier transforms are imaginary (see e.g. Brigham, 1974). The Fourier transforms of the functions in Fig. 1 are displayed in Fig. 2 (showing the imaginary part in the case of H(f) and V~ ðf Þ). To facilitate the comparison of the various Fourier transforms, only the absolute value is considered in Fig. 3. Moreover, logarithmic scales are used on both axes. Owing to the fact that 9R(f)9 is a monotonically decreasing function with a relatively flat profile at lower frequencies, 9H(f)9 and 9H2(f)9 are similar to 9V~ ðf Þ9 and 9M(f)9, respectively. In quantitative terms, however, there are significant differences. The functions 9H(f)9 and 9V~ ðf Þ9 have a maximum at 1/(2py) and f^ v ¼ 1=ð2pyt Þ, respectively. The functions 9H2(f)9 and 9M(f)9 share a zero at pffiffiffi ð30Þ fz ¼ 2=ð2pyÞ
-100 0 100 frequency (Hz)
200
Fig. 2. Model in the frequency domain. (a) The Fourier transform of the rate modulation function, R(f), is real valued. (b) The Fourier transform of the normalized MUAP function, H(f), is imaginary (black curve), whereas the Fourier transform of its square, H2(f), is real valued (gray curve). (c) The Fourier transform of the canonical VEMP, V~ ðf Þ, is imaginary (black curve). The Fourier transform of the variance modulation, M(f), is real valued (gray curve).
from which they rise to a maximum at pffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4=3ðt=yÞ2 6 f^ m ¼ 1 2py 1 þ 2ðt=yÞ2
pffiffiffi 6=ð2pyÞ and ð31Þ
respectively. Maxima and zeros are indicated by dotted vertical lines in Fig. 3.
4. Results 4.1. Varying the effective duration of the modulation The above model is now used to study how canonical VEMP and variance modulation are affected by the duration of the modulation, which is controlled by the parameter t. This para~ meter affects not only the waveforms of m(t) and vðtÞ, but also the amplitudes. To facilitate the comparison of different conditions, the model parameter r0 was used to compensate for the t-dependence of the amplitude factor v~ 0 . To be specific, v~ 0 was
B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 286 (2011) 41–49
|R(f)| 10-3
10-4
10-5 |H(f)|
10-2
|H2(f)|
45
minimum at t ¼0. If t is significantly smaller than y, the variance modulation resembles h2(t), whereas it resembles r(t) if t is significantly larger than y. The latter fact is illustrated by the dashed curve, which shows the graph of r(t) for t ¼8 ms (properly rescaled). The t-dependence of the model is examined more closely in Fig. 5. The solid curve in Fig. 5a is based on Eq. (23) and shows the value r0 that is required to keep v~ 0 at the constant level 0.1. The value of r0 cannot be less than 1, because this would be inconsistent with the physiological interpretation of the rate modulation function r(t), as explained above. This constraint imposes a lower limit on t, namely rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! pffiffiffi e2 2y ð32Þ 1 1 v~ 20 tmin ¼ 2 e9v~ 0 9 In the present example, with
v~ 0 ¼ 0:1 and y ¼4 ms, we get
tmin 0:39 ms (indicated by a dotted vertical line). For compar-
10-3
ison, also the short-modulation approximation is considered (dashed curve): plugging Eqs. (16) and (18) into Eq. (24) yields r0 ¼ 23=2 ev~ 0 y=t. The curve in Fig. 5b shows the t-dependence of d, calculated using Eq. (16). For comparison, the short-modulation approximation leads to pffiffiffiffi dt-0 ¼ v~ 0 D ¼ e p=2v~ 0 y 2:41v~ 0 y ð33Þ
10-4
(cf. Eqs. (18) and (24)). With v~ 0 ¼ 0:1 and y ¼4 ms we get dt-0 E 0.96 ms (dashed horizontal line). The figure demon-
~ |V(f)| 10-3 |M(f)| 10-4
10-5 101
102 frequency (Hz)
Fig. 3. Model in the frequency domain. Absolute values displayed in a doublelogarithmic plot. (a) Fourier transform of the rate modulation function. (b) Fourier transform of the normalized MUAP function (black curve) and its square (gray curve). (c) Fourier transform of the canonical VEMP (black curve) and the associated variance modulation (gray curve).
kept constant at the level 1/10. In Fig. 4, time-domain curves are shown on the left and frequency-domain curves on the right. The canonical VEMP is considered in the upper panels, the variance modulation in the bottom panels. The time constant of the normalized MUAP function, y, was kept constant at 4 ms. The solid curves were obtained by varying the parameter t between 1 and 8 ms, in steps of 1 ms. The short-modulation approximation is shown as a dotted curve. The results for the canonical VEMP confirm the analytical conclusion that the parameter t has, in essence, a time scaling effect: an increase in t increases, in the time domain, the span between maximum and minimum and shifts, in the frequency-domain, the position of the maximum absolute value to lower frequencies. Entirely different is the situation for the variance modulation, where it is of crucial importance whether the parameter t is greater or less than y (the special case t ¼ y is shown as a gray curve). For t o y, the time-domain curves show two minima and a maximum in between, as in Fig. 1. For t Z y, by contrast, there is only a single
strates that the discrepancy between exact model and approximation rapidly increases with t. The times of the extrema of VEMP and variance modulation are examined in Fig. 5c (considering only the extrema at nonnegative times). The black solid curve shows the time of the VEMP extremum, t^v . The curve monotonically increases starting from the value y ¼4 ms. The gray curve represents t^ m , the time of the extremum in the case of the variance modulation. Provided that t is less than about 1.5 ms, the two curves basically coincide, but then they rapidly diverge, owing to the fact that t^ m begins to decrease, eventually becoming zero at t ¼ y (indicated by a dotted vertical line). The corresponding amplitudes are shown in Fig. 5d. The ~ t^ v Þ, corresponds canonical VEMP at the time of the extremum, vð to v~ 0 , which was kept at the constant level 0.1 (black line). The variance modulation at the time of the extremum, mðt^m Þ, is ~ t^v Þ, provided that t is not much larger than comparable to vð 4 ms; beyond that point it is roughly proportional to t (gray curve). When analyzing experimental data, t^v can generally be determined more reliably than t^ m . Since for sufficiently small values of t these two times do not differ very much, it could be useful to evaluate m(t) at the time t^v rather than t^m . The consequence of doing so can be assessed from the dashed curve, which shows the t-dependence of mðt^v Þ. The error that is made by ~ evaluating m(t) at the extremum of vðtÞ is negligible for t o2 ms, and it may be acceptable even for t being twice as large, owing to the fact that the extremum of m(t) is rather flat for t E y. 4.2. The diversity of variance modulations for a given VEMP The typical VEMP has a stereotyped time course, dominated by the peaks p13 and n23. This fact forms the basis for the following considerations, in which the canonical VEMP is invariant and has a time course that corresponds to realistic experimental conditions. According to Eqs. (20) and (21), different combinations of the parameters y and t may result in exactly the same canonical VEMP. We now examine to what extent the associated variance modulation may vary.
46
B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 286 (2011) 41–49
0.1
~ |V(f)|
~ v(t)
10-3
0
10-4
-0.1 14
8
8
4 1
0 1
10-3
3
-0.1
|M(f)|
m(t)
2
4
1 2
5
3
6 7
-0.2 -20
-10
10-4
4
8
0
10
20
1
time (ms)
10
100
frequency (Hz)
Fig. 4. Varying the duration of the rate modulation while keeping the normalized MUAP function invariant. Canonical VEMP considered in the upper half, variance modulation considered in the bottom half. Time domain considered on the left, frequency domain on the right. The numbers attached to the curves indicate the value of the duration parameter t.
The curve in Fig. 6a shows the canonical VEMP corresponding to yt ¼5 ms and v~ 0 ¼ 0:1 (cf. Eqs. (20) and (21)). The positive and the negative peak are 10 ms apart, as are the peaks p13 and n23 in a typical VEMP experiment. The solid curves in Fig. 6b show the variance modulation for values of t between 1pand ffiffiffi 4 ms (the gray curve represents the special case y ¼ t ¼ 5= 2 3:54 ms). The short-modulation approximation is represented by the dotted curve. The figure demonstrates that the waveform of m(t) decisively depends on t. Conversely, this means that the time course of the variance modulation can provide a clue as to the duration of the modulation.
varied while the canonical VEMP was kept invariant (corresponding to yt ¼5 ms and v~ 0 ¼ 0:1). In Fig. 7a, the true d is shown as a black solid curve, whereas dest is shown as a dashed curve. For t-0, the two curves approach the same limiting value, which, according to Eq. (33), is about 1.20 ms (dotted line). But as t increases, the two curves diverge. The relative estimation error, 9(dest d)/d9, is shown as a dashed curve in Fig. 7b. It is about 30% at the worst (for t E3.1 ms).
4.3. Usability of the short-modulation approximation
5.1. The VEMP as a filtered rate modulation
In the short-modulation approximation, the modulation is fully characterized by the quantity d, with the absolute value corresponding to the equivalent rectangular duration of the modulation and the sign distinguishing between excitation and inhibition. Since VEMP, v(t), and MUAP function, h(t), differ only by a constant factor (cf. Eq. (10)), there is an easy way to estimate d from experimental data. The basis for doing so is Eq. (11), where h(t) may be replaced by a properly rescaled v(t). If the rescaling is done so that the resulting h(t) has a maximum absolute value of one, where t^v is the time corresponding to that maximum, d may be estimated as
VEMPs are generally measured with the objective of evaluating the function of the vestibular system. Thus, in terms of the present model, it is the rate modulation r(t) that represents the effect of interest. But a VEMP measurement allows us to observe r(t) only indirectly. Figuratively speaking, recording the VEMP means that the rate modulation is observed through a filter with an impulse response corresponding to the MUAP function h(t) (see Eqs. (5) and (13)). The VEMP-associated variance modulation may be considered, accordingly, as the rate modulation observed through a filter with the impulse response h2(t) (see Eqs. (7) and (14)). A pair of equations that closely resembles our Eqs. (5) and (7) was obtained by Schoonhoven et al. (1988), although in a different context: they developed a method for the analysis of the compound action potential of a peripheral nerve (cf. their Eqs. (9) and (10)). It would certainly be desirable to be able to undo the filtering by applying some kind of inverse filter. If the function h(t) were known, r(t) could indeed be largely recovered by deconvolution, although the ill-posed nature of the inverse problem to be solved would require special provisions, such as optimal (Wiener)
dest ¼ Dmðt^v Þ
ð34Þ
The short-modulation approximation implies that a strong modulation of short duration may have exactly the same effect as a weaker modulation of longer duration. However, the results of this modeling study partially contradict that simple view, and so the question arises as to how closely dest matches the true value of d as calculated using Eq. (16). This question is investigated in Fig. 7. As in the previous figure, the parameter t was
5. Discussion
B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 286 (2011) 41–49
47
0.1
0
a
~ v(t)
Parameter r0
-0.2 -0.4
0
-0.6 -0.1
-0.8
0
-1
m(t)
2
δ (ms)
-2
3
-0.1
-20
-5 Time of Extremum (ms)
3.54 4
-3 -4
-10
0 time (ms)
10
20
Fig. 6. Varying the duration of the rate modulation while keeping the canonical VEMP invariant. The canonical VEMP shown in (a) may be associated with each of the variance modulations shown in (b). The numbers attached to the curves indicate the value of the duration parameter t.
8 t
6 -1
4
a
t
2 0
-2 δ (ms)
v(t )
Amplitude of Extremum
b
1
-1
-0.1
-3
m(t )
m(t )
-4
-0.2 2
4 τ (ms)
6
8
Fig. 5. Effect of the duration parameter t (same model as in Fig. 4). (a) The value r0 that is required to keep v~ 0 at the constant level 0.1 (cf. Eq. (22)). (b) Parameter d calculated using Eq. (16). (c) Times of the extrema of VEMP (black curve) and variance modulation (gray curve). (d) Amplitudes of the extrema of VEMP (black line) and variance modulation (gray curve). The dashed curve shows the variance modulation at the time of the VEMP extremum.
filtering (see e.g. Press et al., 2007). But h(t) is generally unknown, which vastly complicates the situation. At least partial information about h(t) could be easily extracted from the data. On the assumption that all MUAPs have the same waveform, the modulus of a Fourier-transformed MUAP, which would be 9H(f)9 in the notation of the present article, can be expected to have the same shape as the EMG spectrum (Wit and Kingma, 2006). Fig. 3 may help to understand what possibilities this opens up: given 9V~ ðf Þ9 (derived from the measured VEMP, v(t)) and an estimate of 9H(f)9, it is possible to estimate 9R(f)9, scaling issues left aside. While there is, of course, no unique way to reconstruct a time-domain function from the modulus of its Fourier transform, the above approach might provide valuable constraints. A different approach
30
relative estimation error (percent)
0
b
20
10
0 0
1
2
3
4
5
τ (ms) Fig. 7. Estimation of d using the short-modulation approximation (simulation for yt ¼5 ms and v~ 0 ¼ 0:1). (a) Comparison between true d (black curve) and estimate dest (dashed curve). (b) Relative estimation error for dest. The worst-case error is below 30%.
48
B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 286 (2011) 41–49
is suggested by the frequency-domain representation of our general model. Combining Eqs. (13) and (14) yields H2 ðf Þ V~ ðf Þ ¼ Hðf ÞMðf Þ
ð35Þ
which in the time-domain reads Z þ1 Z þ1 0 ~ h2 ðt 0 Þvðtt Þdt0 ¼ hðt 0 Þmðtt 0 Þdt0 1
ð36Þ
1
Thus, h(t) is formally the solution of the non-linear integral equation (36). However, this inverse problem can be assumed to be seriously ill-posed, which means that it is a challenging task to find a solution which not only explains the data, but is also reasonable from a physiological point of view. With this in mind, it is remarkable that the short-modulation approximation offers an amazingly simple way of escaping these difficulties. 5.2. The duration of the rate modulation matters The main impetus for this study was the concern that, in typical experiments, the duration of the rate modulation might be too long to be reconcilable with the assumption of a short rate modulation. Stimulating with clicks of 0.1 ms duration, Colebatch and Rothwell (2004) found a mean duration of 3.6 ms for the resulting inhibition of motor unit activity. If tone pulses are used rather than clicks, the inhibition might last even somewhat longer. The present analyses show that such durations are indeed not short in the sense of the model. In the simple Gaussian model considered here, the duration of the rate modulation is controlled by the parameter t. As to the VEMP, t has, in essence, a time-scaling effect: the VEMP has exactly the same waveform as the MUAP function h(t), apart from a t-dependent time dilation (see Fig. 4). Thus, supposing that h(t) is unknown, the VEMP gives no clue as to the value of t. Quite different is the situation in the case of the variance modulation, as best seen in Fig. 6: one and the same VEMP may be associated with fundamentally different variance modulations. This is clearly contradictory to a basic feature of the short-modulation approximation, namely that the variance modulation corresponds to the VEMP squared (apart from a scaling factor). As a matter of fact, a clear deviation from the short-modulation approximation (represented by the dotted curve in Fig. 6b) is found already for t ¼ 1 ms. 5.3. Making a case for the short-modulation approximation The fact that the duration of the rate modulation may well have relevant effects on both the VEMP and the associated variance modulation seems to render the short-modulation approximation invalid. However, such a conclusion would be premature. Without knowledge of the MUAP function h(t), most effects are not easily recognized as such. Thus, under realistic experimental conditions (i.e. in the case of noisy data), it may be difficult, if not impossible, to falsify the short-modulation approximation. This is particularly true for the VEMP. For the model considered in this study, the VEMP is basically a time-dilated version of h(t) so that, without using additional information or making assumptions, it is absolutely impossible to reach conclusions about the duration of the rate modulation. The assumption of an infinitesimal duration is, in that situation, just as arbitrary as any other assumption, but it has the inestimable advantage that the resulting model (here called ‘short-modulation approximation’) is extraordinarily simple from a mathematical point of view: getting the MUAP function h(t) merely requires rescaling of the VEMP. Thus, as to the latter, the short-modulation approximation is, by definition, in perfect agreement with the experimental data. The crucial question, consequently, is to what extent
the approximation allows us to explain the variance modulation. Fig. 6 provides an answer. The duration parameter t was systematically varied between 0 and 4 ms while the VEMP curve was kept constant (by concurrently changing other model parameters). At the times of the two VEMP peaks (dotted vertical lines in Fig. 6) and beyond (towards 7N), the effect of t on the variance modulation is relatively minor so that the differences between the curves may be considered to be insignificant (unless data with an outstanding signal-to-noise ratio are available). The short-modulation approximation consequently allows us to make reasonable predictions regarding the peak amplitudes, and for this reason it can be used to estimate the equivalent rectangular duration of the rate modulation, 9d9. The model calculations presented in Fig. 7b resulted in a worst-case estimation error of less than 30% (dashed curve). This appears to be quite acceptable, considering that much greater uncertainties are presumably introduced by the simplifying assumptions on which the model is based. 5.4. Limitations of the simple Gaussian model The simple Gaussian model may help to interpret experimental findings qualitatively, but it has not been devised to match measured data quantitatively. Apart from scaling factors, the model has only two parameters (the time constants y and t) so that the possibilities for simultaneously fitting the VEMP and the associated variance modulation are limited. Moreover, the choice of the functions h(t) and r(t) was motivated more by mathematical convenience than by physiological considerations. Thus, these functions are most likely suboptimal with regard to real conditions. On the other hand, given our current state knowledge, some arbitrariness in the definition of the two functions is unavoidable, and we are not aware of another simple model that would give us an obvious advantage over the one we investigated here.
6. Conclusions The simple Gaussian model introduced in this study is arguably one of the most straightforward VEMP models conceivable. It stands out in that all pertinent aspects can be examined analytically, often on the basis of surprisingly simple formulas. Despite the simplicity of the model, the predicted VEMP curves closely resemble those experimentally observed. Thus, all in all, the model appears to be ideally suited for getting a basic theoretical understanding of the VEMP phenomenon. This does not only apply to questions arising from previous work. Considering the fact that the model has a simple representation also in the frequency-domain, it might provide guidance for future work on VEMP deconvolution. According to the underlying general model, both the VEMP and the associated variance modulation may be considered as filtered versions of the effect of interest, which is represented by the temporal modulation of the MUAP rate, r(t). The main component of the VEMP experimentally observed, i.e. the p13–n23 complex, has roughly the same appearance as the VEMP predicted by our simple Gaussian model. According to this model, the time course of the VEMP mainly reflects the impulse response of the filter, which is the MUAP function h(t), whereas the rate modulation r(t) has basically a smoothing and broadening effect that depends on the duration of the rate modulation. The same applies to the variance modulation, where the filter corresponds to h2(t). But owing to the fact that this filter has a non-negative impulse response, a peculiarity arises: with increasing duration of the rate modulation, the variance modulation gradually loses its similarity
B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 286 (2011) 41–49
with the filter function h2(t) and finally looks more like the rate modulation r(t). The effect would, in principle, allow us to estimate the duration of the rate modulation. However, our ¨ ¨ experimental data (Lutkenh oner et al., 2011) point out complications, presumably due to the multi-component structure of the VEMP, which was disregarded in this study. The duration of the rate modulation does not affect all aspects of the VEMP to the same degree. A moderate duration increase affects mainly the amplitudes of VEMP and variance modulation, whereas the effect on the waveforms may be insignificant compared to the inevitable noise in the data. This study suggests that the effect on the peak amplitudes is well explained by the short-modulation approximation of the model, which offers an exceptionally simple and convenient way for a model-based interpretation of experimental results. References Boashash, B., 1992. Estimating and interpreting the instantaneous frequency of a signal. 1. Fundamentals. Proc. IEEE 80, 520–538. Brantberg, K., Verrecchia, L., 2009. Testing vestibular-evoked myogenic potentials with 90-dB clicks is effective in the diagnosis of superior canal dehiscence syndrome. Audiol. Neurootol. 14, 54–58. Brigham, E.O., 1974. The Fast Fourier Transform. Prentice-Hall, Englewood Cliffs, NJ. Colebatch, J.G., Halmagyi, G.M., Skuse, N.F., 1994. Myogenic potentials generated by a click-evoked vestibulocollic reflex. J. Neurol. Neurosurg. Psychiatry 57, 190–197.
49
Colebatch, J.G., Day, B.L., Bronstein, A.M., Davies, R.A., Gresty, M.A., Luxon, L.M., Rothwell, J.C., 1998. Vestibular hypersensitivity to clicks is characteristic of the Tullio phenomenon. J. Neurol. Neurosurg. Psychiatry 65, 670–678. Colebatch, J.G., Rothwell, J.C., 2004. Motor unit excitability changes mediating vestibulocollic reflexes in the sternocleidomastoid muscle. Clin. Neurophysiol. 115, 2567–2573. Day, S.J., Hulliger, M., 2001. Experimental simulation of cat electromyogram: evidence for algebraic summation of motor-unit action-potential trains. J. Neurophysiol. 86, 2144–2158. Gabor, D., 1946. Theory of communication. J. Inst. Electr. Eng. (London) 93, 429–457. Kingma, C.M., Wit, H.P., 2011. Asymmetric vestibular evoked myogenic potentials in unilateral Meniere patients. Eur. Arch. Otorhinolaryngol. 268, 57–61. ¨ ¨ Lutkenh oner, B., Stoll, W., Basel, T., 2010. Modeling the vestibular evoked myogenic potential. J. Theor. Biol. 263, 70–78. ¨ ¨ Lutkenh oner, B., Rudack, C., Basel, T., 2011. The variance modulation associated with the vestibular evoked myogenic potential. Clin. Neurophysiol 122, 1448–1456. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 2007. Numerical Recipes. The Art of Scientific Computing third ed. Cambridge University Press, Cambridge. Rosengren, S.M., Welgampola, M.S., Colebatch, J.G., 2010. Vestibular evoked myogenic potentials: past, present and future. Clin. Neurophysiol. 121, 636–651. Schoonhoven, R., Stegeman, D.F., van Oosterom, A., Dautzenberg, G.F., 1988. The inverse problem in electroneurography—I: conceptual basis and mathematical formulation. IEEE Trans. Biomed. Eng. 35, 769–777. Ushio, M., Iwasaki, S., Murofushi, T., Sugasawa, K., Chihara, Y., Fujimoto, C., Nakamura, M., Yamaguchi, T., Yamasoba, T., 2009. The diagnostic value of vestibular-evoked myogenic potential in patients with vestibular schwannoma. Clin. Neurophysiol. 120, 1149–1153. Wit, H.P., Kingma, C.M., 2006. A simple model for the generation of the vestibular evoked myogenic potential (VEMP). Clin. Neurophysiol. 117, 1354–1358.
Contents lists available at ScienceDirect
Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi
An analytical model of the vestibular evoked myogenic potential n ¨ ¨ ¨ Bernd Lutkenh oner , Turker Basel
ENT Clinic, M¨ unster University Hospital, Kardinal-von-Galen-Ring 10, 48129 M¨ unster, Germany
a r t i c l e i n f o
abstract
Article history: Received 10 January 2011 Received in revised form 27 April 2011 Accepted 6 July 2011 Available online 14 July 2011
The vestibular evoked myogenic potential (VEMP) can be modeled (scaling factors aside) as a convolution of the motor unit action potential (MUAP) of a representative motor unit, h(t), with the temporal modulation of the MUAP rate of all contributing motor units, r(t). Accordingly, the variance modulation associated with the VEMP can be modeled as a convolution of r(t) with the square of h(t). To get a deeper theoretical understanding of the VEMP phenomenon, a specific realization of this general model is investigated here. Both r(t) and h(t) were derived from a Gaussian probability density function (in the latter case taking the first derivative). The resulting model turned out to be simple enough to be evaluated analytically in the time and in the frequency domain, while still being realistic enough to account for the basic aspects of the VEMP generation. Perhaps the most significant conclusion of this study is that, in the case of noisy data, it may be difficult to falsify the hypothesis of a rate modulation of infinitesimal duration. Thus, certain aspects of the data (particularly the peak amplitudes) can be interpreted using a short-modulation approximation rather than the general model. The importance of this realization arises from the fact that the approximation offers an exceptionally simple and convenient way for a model-based interpretation of experimental data, whereas any attempt to use the general model for that purpose would result in an ill-posed inverse problem that is far from easy to solve. & 2011 Elsevier Ltd. All rights reserved.
Keywords: VEMP modeling Canonical VEMP Variance modulation Motor unit action potential Electromyogram
1. Introduction Reflexes elicited by stimulating the otoliths in the inner ear with high-level sound pulses modulate the tone of various muscles. The accompanying short-term perturbation of the electromyogram (EMG) can be measured as the vestibular evoked myogenic potential (VEMP). The main component of the VEMP, classically recorded over the sternocleidomastoid muscles (Colebatch et al., 1994), arises from a brief inhibition of the muscle tone. Thus, a sine qua non for a successful measurement is a tonic muscle contraction. Recording of the VEMP has become a standard clinical test of otolith function, with documented diagnostic utility in varied vestibular and central nervous system disorders (Rosengren et al., 2010). VEMPs with reduced amplitude on the affected side were found, for example, in patients with vestibular schwannoma (Ushio et al., 2009) and unilateral Menie re’s disease (Kingma and Wit, 2011). Patients with Tullio phenomenon (hypersensitivity of the vestibular system to sound), by contrast, show an augmented VEMP amplitude on the affected side, associated with an abnormally low VEMP threshold (Colebatch et al., 1998; Brantberg and Verrecchia, 2009).
n
Corresponding author. Tel.: þ49 251 83 56864; fax: þ49 251 83 56882. ¨ ¨ E-mail address: [email protected] (B. Lutkenh oner).
0022-5193/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2011.07.002
Wit and Kingma (2006) developed a simple model for the generation of the VEMP. A key constituent of the model is the motor unit action potential (MUAP). All motor units are assumed to have the same MUAP, apart from a unit-dependent factor that is thought to account for varying distances between motor unit and recording electrodes. The EMG, being the algebraic sum of the MUAP trains of all motor units (Day and Hulliger, 2001), can be simulated by adding MUAPs with random occurrence times and amplitudes. As the MUAPs are zero-mean signals, their random occurrence causes partial cancellation (positive and negative waves tend to extinguish each other). What remains is stochastic activity with an expected value of zero. A simple modification transforms this EMG model into a model for the VEMP: Instead of being constant in time, the probability of generating a MUAP drops to zero for a short period, corresponding to complete inhibition of the MUAP generation. The inhomogeneity in the MUAP rate disturbs the partial cancellation of MUAPs so that the expected value of the EMG temporarily deviates from zero. This is what is known as the VEMP. Building upon the pioneering work of Wit and Kingma (2006), we developed a somewhat more general and, with respect to ¨ ¨ computation time, more efficient model (Lutkenh oner et al., 2010). Moreover, we showed that the model can be greatly simplified by assuming that the duration of the inhibition is short (as compared with the time constants that determine the MUAP).
42
B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 286 (2011) 41–49
A comparison between model and data basically confirmed the ¨ ¨ theoretical predictions (Lutkenh oner et al., 2011). In particular, we were able to verify that the sound-induced inhibition affects not only the mean of the EMG (causing the VEMP), but also the variance. Some findings were unanticipated, though. For the most part, the observed complications can be attributed to the VEMP being composed of multiple components. But there is another point deserving attention: the assumed short duration of the inhibition. In a study of the firing of single motor units, Colebatch and Rothwell (2004) found inhibitions with durations between 2 and 8 ms (mean 3.6 ms). This may be too long to be compatible with the key assumption underlying the short-inhibition approximation, namely that the MUAP function is roughly constant for the whole duration of the inhibition. The main purpose of the present study was to develop an understanding of how the duration of the inhibition impacts on the VEMP and the associated variance modulation. The model presented here was chosen so that it is, on the one hand, realistic enough to account for the basic aspects of the VEMP generation, and, on the other hand, simple enough to allow analytical evaluation.1 Special attention is given to the question whether the short-inhibition approximation is applicable even when the underlying assumption is not fully met. This question is important because the approximation offers an exceptionally simple and convenient way for a model-based interpretation of measured data. The general model, by contrast, leads to a non-linear deconvolution problem that is far from easy to solve. With deconvolution in mind, our analytical model is considered in both the time- and the frequency domain so that it may serve as a reference model also for future work on VEMP deconvolution.
2. General theory To begin with, we briefly recall the theory developed pre¨ ¨ viously (Lutkenh oner et al., 2010). The notations were partially changed, to better meet the requirements of the current article. Moreover, the theory was slightly generalized so that both inhibition and excitation can be considered. Besides, the model is also considered in the frequency domain. 2.1. General model in the time domain The model is based on the simplifying assumption that the MUAPs are identical for all motor units, apart from an amplitude factor that may differ from unit to unit. The MUAP time course is described by a dimensionless zero-mean function h(t). The normalization of h(t) is arbitrary. Here, we normalize so that the maximum absolute value is one. The quantity Z þ1 D¼ h2 ðtÞdt ð1Þ
2
s20 ¼ a r0 D
ð2Þ
Note that r0D is, loosely speaking, the average number of MUAPs significantly overlapping in time. If the rate exhibits a temporal modulation r(t) such that
rðtÞ ¼ r0 ð1 þrðtÞÞ
ð3Þ
then the mean of the EMG, i.e. the VEMP, is ~ vðtÞ ¼ ar0 DvðtÞ
ð4Þ
with ~ ¼ D1 vðtÞ
Z
þ1
rðt 0 Þhðtt 0 Þdt 0
ð5Þ
1
For the variance we get
s2 ðtÞ ¼ s20 ð1 þ mðtÞÞ
ð6Þ
with mðtÞ ¼ D1
Z
þ1
rðt 0 Þh2 ðtt 0 Þdt 0
ð7Þ
1
The functions r(t) and m(t) will be referred to as the rate modulation and the variance modulation, respectively. Thepstanffiffiffiffiffiffiffiffiffi dard deviation of the EMG is, according to Eq. (2), s0 ¼ a r0 D. Dividing the VEMP signal v(t) by s0 yields the normalized VEMP qffiffiffiffiffiffiffiffiffi ~ vnorm ðtÞ ¼ r0 Da=avðtÞ ð8Þ The normalized VEMP is a dimensionless measure, as is the variance modulation. Yet there is an important difference. Unlike the variance modulation, the normalized VEMP depends on the base rate r0, i.e. on the muscle tone. In a modeling study, this ~ complication can be avoided by considering vðtÞ instead of v(t) or vnorm(t). If needed, the latter two functions can be obtained by ~ simple rescaling. With regard to the fact that the function vðtÞ represents the VEMP reduced to the simplest form possible, it will be referred to as the canonical VEMP. A comparison of Eqs. (5) and (7) shows that canonical VEMP and variance modulation make a perfectly matched pair. The ratio a=a is generally close to 1, and a pffiffiffiffiffiffiffiffiffi r0 D would be about 4 (corresponding to typical value for D ¼0.01 s and r0 ¼1600 s 1). Thus, normalized and canonical VEMP may be expected to differ by a factor of the order of 4. 2.2. Approximation for short rate modulations
1
may then be interpreted as the effective MUAP duration, because it measures how long h(t) has an absolute value close to maximum.2 The MUAPs of individual motor units may be 1 More tedious calculations were assisted by the computational software Mathematica 7 (Wolfram Research, Inc.). 2 Eq. (1) is not meant as a definition of effective duration for arbitrary signals. However, for a temporally compact, normalized signal such as the MUAP function h(t) considered here, the integral on the right-hand side may be interpreted that way. For the MUAP function introduced in Eq. (17), calculation of an effective duration according to Gabor (1946) yields
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Z Z þ1
Dt ¼
obtained by multiplying h(t) by a unit-dependent amplitude factor. Detailed knowledge about the distribution of this factor is not required. Only the arithmetic mean, a, and the quadratic mean, a, are needed. The above assumptions allow us to integrate the MUAP trains of the individual motor units into a single train. In concrete terms, we assume that the motor units as a whole fire at a rate r(t). In the case of a constant rate, which means r(t) r0, the recorded EMG has zero mean and variance
þ1
t 2 h2 ðtÞdt
2p
1
h2 ðtÞdt ¼
1
pffiffiffiffiffiffi 3py
If the rate modulation r(t) corresponds to a pulse that is so short that h(t) remains approximately constant for the entire pulse duration, Eqs. (5) and (7) can be greatly simplified. The modulation effect turns out to be proportional to Z þ1 d¼ rðtÞdt ð9Þ 1
(footnote continued) The D resulting differ by pffiffiffiffiffiffi pffiffiffi from Eq. (1) is provided in Eq. (18). The two measures the factor 2 3=e 1:27. Regarding the fact that the factor 2p in Gabor’s definition is relatively arbitrary – other authors omit this factor (see, e.g., Eq. (53) of Boashash (1992)) – our usage of the term effective duration appears to be justified.
B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 286 (2011) 41–49
where 9d9 is the equivalent rectangular duration of the rate modulation. The sign of d indicates the type of modulation (positive sign for excitation and negative sign for inhibition). The equations for canonical VEMP and variance modulation reduce to ~ ¼ dD1 hðtÞ vðtÞ
ð10Þ
mðtÞ ¼ dD1 h2 ðtÞ
ð11Þ
43
0
While these two equations are strictly valid only for a rate modulation of infinitesimal duration, they may be expected to be at least approximately valid when the duration of the rate modulation is short compared to the time constants that determine the MUAP function h(t). This is why they will be referred to as the short-modulation approximation. What ‘short’ actually means and whether the approximation is applicable to typical experimental conditions will be examined in the model simulations below.
r(t) -1 1 h2(t)
0
2.3. General model in the frequency domain Canonical VEMP and variance modulation may be considered as filtered versions of the rate modulation r(t), because they are obtained by a convolution with h(t) and h2(t), respectively (cf. Eqs. (5) and (7)). This filtering point of view is best analyzed in the frequency domain. The Fourier transform of a function x(t) is Z 1 xðtÞexpð2piftÞdt ð12Þ Xðf Þ ¼
h(t) -1 0.5 ~ v(t)
1
~ (e.g. Brigham, 1974). If the Fourier transforms of vðtÞ, m(t), r(t), h(t), and h2(t) are denoted as V~ ðf Þ, M(f), R(f), H(f), and H2(f), respectively, the convolution theorem (see e.g. Brigham, 1974) allows us to rewrite Eqs. (5) and (7) as V~ ðf Þ ¼ D1 Rðf ÞHðf Þ
ð13Þ
Mðf Þ ¼ D1 Rðf ÞH2 ðf Þ
ð14Þ
0
m(t) -0.5 -10
3. Simple Gaussian model 3.1. Time domain To put the general theory outlined above into concrete terms, specific assumptions about the rate modulation function r(t) and the normalized MUAP function h(t) will be made now. The rate modulation function was derived from the probability density function of a Gaussian distribution t2 ð15Þ rðtÞ ¼ r0 exp 2 2t The parameter t controls the width of the function, 9r09 determines the maximum modulation effect, and the sign of r0 distinguishes between inhibition (negative sign) and excitation. Since complete inhibition (i.e. r(t)¼ 1) represents a limiting condition, the requirement r0 Z 1 has to be fulfilled. Solving the integral in Eq. (9) yields pffiffiffiffiffiffi d ¼ 2pr0 t ð16Þ Fig. 1a shows the graph of r(t) for r0 ¼ 1 and t ¼3 ms. For h(t) we choose the first derivative of a Gaussian, as suggested by Wit and Kingma (2006): ! t ðt=yÞ2 1 hðy; tÞ ¼ exp ð17Þ y 2
0 time (ms)
10
Fig. 1. Model in the time-domain. (a) Rate modulation function r(t). (b) Normalized MUAP function h(t) (black curve) and its square (gray curve). ~ (c) Canonical VEMP vðtÞ (black curve) and variance modulation m(t) (gray curve).
The black curve in Fig. 1b shows this function for y ¼4 ms. There is a minimum at y and a maximum at y (dotted vertical lines). The function is normalized so that the values at minimum and maximum are 1 and 1, respectively. The square of h(t) is displayed as a gray curve. Solving the integral in Eq. (1) yields
D¼
pffiffiffiffi e p y 2:41y 2
ð18Þ
Thus, in the present example, the MUAP function has an effective duration of about 9.64 ms. The canonical VEMP is obtained by solving the integral in Eq. (5) ~ ¼ vðtÞ
23=2 e1=2 ytr0 2
ðy þ t2 Þ3=2
t exp
t2 2
!
2ðy þ t2 Þ
ð19Þ
The black curve in Fig. 1c shows this function for r0 ¼ 1, t ¼3 ms ~ has and y ¼4 ms. A closer inspection of Eq. (19) reveals that vðtÞ basically the same mathematical structure as h(y;t). Indeed, a
44
B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 286 (2011) 41–49
simple transformation yields
0
~ ¼ v~ 0 hðyt ; tÞ vðtÞ
ð20Þ
a
with
yt ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi y2 þ t2
ð21Þ
R(f)
and 2 v~ 0 ¼ 23=2 e1 r0 yt=yt
ð22Þ
~ are at t^ v ¼ yt and t^ v (dashed vertical lines The extrema of vðtÞ in Fig. 1c). Eq. (22) may be rewritten as )1 pffiffiffi ( ðytÞ2 2 1þ r0 ð23Þ v~ 0 ¼ e 2yt
-0.008 0.01
v~ 0 ¼ d=D
b
H2(f)
which shows that the maximum condition for v~ 0 =r0 is y ¼ t (for any other combination of y and t, the term in curly brackets is greater than one). Thus, the global maximum of v~ 0 is pffiffiffi 2=e9r0 9 0:529r0 9. In the limit t-0, Eq. (22) may be rewritten as
0
ð24Þ
consistent with Eq. (10). The variance modulation is obtained by solving the integral in Eq. (7) ! t3 ty2 t2 t2 ð25Þ exp 2 þ mðtÞ ¼ 23=2 r0 2 2 3=2 5=2 ðy þ2t2 Þ ðy þ 2t2 Þ y þ 2t2
Im[H(f)] -0.01 0.008
0
The latter times are marked by short vertical lines in Fig. 1c (almost coinciding with the dotted vertical lines). For t Z y, by contrast, 9m(t)9 has is a single minimum at t¼0.
M(f) -0.008
3.2. Frequency domain
-200
The Fourier transform of r(t) is pffiffiffiffiffiffi Rðf Þ ¼ 2p t r0 expð2p2 t2 f 2 Þ
ð27Þ
and the Fourier transforms of h(t) and h2(t) are 2
c
~ Im[V(f)]
The function is displayed as a gray curve in Fig. 1c (same parameters as used for the canonical VEMP). For t o y, 9m(t)9 has a minimum at t ¼0 and maxima at 7 t^ m with qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi t^m ¼ y 1 þ ðt=yÞ2 2ðt=yÞ4 ð26Þ
2
Hðf Þ ¼ ið2pÞ3=2 e1=2 y f expð2p2 y f 2 Þ
ð28Þ
pffiffiffiffi 2 2 H2 ðf Þ ¼ e py=2ð12p2 y f 2 Þexpðp2 y f 2 Þ
ð29Þ
The Fourier transforms of canonical VEMP and variance modulation are easily obtained using Eqs. (13) and (14). Since the functions r(t), h2(t), and m(t) are even functions, their Fourier ~ transforms are real-valued. The functions h(t) and vðtÞ, by contrast, are odd so that their Fourier transforms are imaginary (see e.g. Brigham, 1974). The Fourier transforms of the functions in Fig. 1 are displayed in Fig. 2 (showing the imaginary part in the case of H(f) and V~ ðf Þ). To facilitate the comparison of the various Fourier transforms, only the absolute value is considered in Fig. 3. Moreover, logarithmic scales are used on both axes. Owing to the fact that 9R(f)9 is a monotonically decreasing function with a relatively flat profile at lower frequencies, 9H(f)9 and 9H2(f)9 are similar to 9V~ ðf Þ9 and 9M(f)9, respectively. In quantitative terms, however, there are significant differences. The functions 9H(f)9 and 9V~ ðf Þ9 have a maximum at 1/(2py) and f^ v ¼ 1=ð2pyt Þ, respectively. The functions 9H2(f)9 and 9M(f)9 share a zero at pffiffiffi ð30Þ fz ¼ 2=ð2pyÞ
-100 0 100 frequency (Hz)
200
Fig. 2. Model in the frequency domain. (a) The Fourier transform of the rate modulation function, R(f), is real valued. (b) The Fourier transform of the normalized MUAP function, H(f), is imaginary (black curve), whereas the Fourier transform of its square, H2(f), is real valued (gray curve). (c) The Fourier transform of the canonical VEMP, V~ ðf Þ, is imaginary (black curve). The Fourier transform of the variance modulation, M(f), is real valued (gray curve).
from which they rise to a maximum at pffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4=3ðt=yÞ2 6 f^ m ¼ 1 2py 1 þ 2ðt=yÞ2
pffiffiffi 6=ð2pyÞ and ð31Þ
respectively. Maxima and zeros are indicated by dotted vertical lines in Fig. 3.
4. Results 4.1. Varying the effective duration of the modulation The above model is now used to study how canonical VEMP and variance modulation are affected by the duration of the modulation, which is controlled by the parameter t. This para~ meter affects not only the waveforms of m(t) and vðtÞ, but also the amplitudes. To facilitate the comparison of different conditions, the model parameter r0 was used to compensate for the t-dependence of the amplitude factor v~ 0 . To be specific, v~ 0 was
B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 286 (2011) 41–49
|R(f)| 10-3
10-4
10-5 |H(f)|
10-2
|H2(f)|
45
minimum at t ¼0. If t is significantly smaller than y, the variance modulation resembles h2(t), whereas it resembles r(t) if t is significantly larger than y. The latter fact is illustrated by the dashed curve, which shows the graph of r(t) for t ¼8 ms (properly rescaled). The t-dependence of the model is examined more closely in Fig. 5. The solid curve in Fig. 5a is based on Eq. (23) and shows the value r0 that is required to keep v~ 0 at the constant level 0.1. The value of r0 cannot be less than 1, because this would be inconsistent with the physiological interpretation of the rate modulation function r(t), as explained above. This constraint imposes a lower limit on t, namely rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! pffiffiffi e2 2y ð32Þ 1 1 v~ 20 tmin ¼ 2 e9v~ 0 9 In the present example, with
v~ 0 ¼ 0:1 and y ¼4 ms, we get
tmin 0:39 ms (indicated by a dotted vertical line). For compar-
10-3
ison, also the short-modulation approximation is considered (dashed curve): plugging Eqs. (16) and (18) into Eq. (24) yields r0 ¼ 23=2 ev~ 0 y=t. The curve in Fig. 5b shows the t-dependence of d, calculated using Eq. (16). For comparison, the short-modulation approximation leads to pffiffiffiffi dt-0 ¼ v~ 0 D ¼ e p=2v~ 0 y 2:41v~ 0 y ð33Þ
10-4
(cf. Eqs. (18) and (24)). With v~ 0 ¼ 0:1 and y ¼4 ms we get dt-0 E 0.96 ms (dashed horizontal line). The figure demon-
~ |V(f)| 10-3 |M(f)| 10-4
10-5 101
102 frequency (Hz)
Fig. 3. Model in the frequency domain. Absolute values displayed in a doublelogarithmic plot. (a) Fourier transform of the rate modulation function. (b) Fourier transform of the normalized MUAP function (black curve) and its square (gray curve). (c) Fourier transform of the canonical VEMP (black curve) and the associated variance modulation (gray curve).
kept constant at the level 1/10. In Fig. 4, time-domain curves are shown on the left and frequency-domain curves on the right. The canonical VEMP is considered in the upper panels, the variance modulation in the bottom panels. The time constant of the normalized MUAP function, y, was kept constant at 4 ms. The solid curves were obtained by varying the parameter t between 1 and 8 ms, in steps of 1 ms. The short-modulation approximation is shown as a dotted curve. The results for the canonical VEMP confirm the analytical conclusion that the parameter t has, in essence, a time scaling effect: an increase in t increases, in the time domain, the span between maximum and minimum and shifts, in the frequency-domain, the position of the maximum absolute value to lower frequencies. Entirely different is the situation for the variance modulation, where it is of crucial importance whether the parameter t is greater or less than y (the special case t ¼ y is shown as a gray curve). For t o y, the time-domain curves show two minima and a maximum in between, as in Fig. 1. For t Z y, by contrast, there is only a single
strates that the discrepancy between exact model and approximation rapidly increases with t. The times of the extrema of VEMP and variance modulation are examined in Fig. 5c (considering only the extrema at nonnegative times). The black solid curve shows the time of the VEMP extremum, t^v . The curve monotonically increases starting from the value y ¼4 ms. The gray curve represents t^ m , the time of the extremum in the case of the variance modulation. Provided that t is less than about 1.5 ms, the two curves basically coincide, but then they rapidly diverge, owing to the fact that t^ m begins to decrease, eventually becoming zero at t ¼ y (indicated by a dotted vertical line). The corresponding amplitudes are shown in Fig. 5d. The ~ t^ v Þ, corresponds canonical VEMP at the time of the extremum, vð to v~ 0 , which was kept at the constant level 0.1 (black line). The variance modulation at the time of the extremum, mðt^m Þ, is ~ t^v Þ, provided that t is not much larger than comparable to vð 4 ms; beyond that point it is roughly proportional to t (gray curve). When analyzing experimental data, t^v can generally be determined more reliably than t^ m . Since for sufficiently small values of t these two times do not differ very much, it could be useful to evaluate m(t) at the time t^v rather than t^m . The consequence of doing so can be assessed from the dashed curve, which shows the t-dependence of mðt^v Þ. The error that is made by ~ evaluating m(t) at the extremum of vðtÞ is negligible for t o2 ms, and it may be acceptable even for t being twice as large, owing to the fact that the extremum of m(t) is rather flat for t E y. 4.2. The diversity of variance modulations for a given VEMP The typical VEMP has a stereotyped time course, dominated by the peaks p13 and n23. This fact forms the basis for the following considerations, in which the canonical VEMP is invariant and has a time course that corresponds to realistic experimental conditions. According to Eqs. (20) and (21), different combinations of the parameters y and t may result in exactly the same canonical VEMP. We now examine to what extent the associated variance modulation may vary.
46
B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 286 (2011) 41–49
0.1
~ |V(f)|
~ v(t)
10-3
0
10-4
-0.1 14
8
8
4 1
0 1
10-3
3
-0.1
|M(f)|
m(t)
2
4
1 2
5
3
6 7
-0.2 -20
-10
10-4
4
8
0
10
20
1
time (ms)
10
100
frequency (Hz)
Fig. 4. Varying the duration of the rate modulation while keeping the normalized MUAP function invariant. Canonical VEMP considered in the upper half, variance modulation considered in the bottom half. Time domain considered on the left, frequency domain on the right. The numbers attached to the curves indicate the value of the duration parameter t.
The curve in Fig. 6a shows the canonical VEMP corresponding to yt ¼5 ms and v~ 0 ¼ 0:1 (cf. Eqs. (20) and (21)). The positive and the negative peak are 10 ms apart, as are the peaks p13 and n23 in a typical VEMP experiment. The solid curves in Fig. 6b show the variance modulation for values of t between 1pand ffiffiffi 4 ms (the gray curve represents the special case y ¼ t ¼ 5= 2 3:54 ms). The short-modulation approximation is represented by the dotted curve. The figure demonstrates that the waveform of m(t) decisively depends on t. Conversely, this means that the time course of the variance modulation can provide a clue as to the duration of the modulation.
varied while the canonical VEMP was kept invariant (corresponding to yt ¼5 ms and v~ 0 ¼ 0:1). In Fig. 7a, the true d is shown as a black solid curve, whereas dest is shown as a dashed curve. For t-0, the two curves approach the same limiting value, which, according to Eq. (33), is about 1.20 ms (dotted line). But as t increases, the two curves diverge. The relative estimation error, 9(dest d)/d9, is shown as a dashed curve in Fig. 7b. It is about 30% at the worst (for t E3.1 ms).
4.3. Usability of the short-modulation approximation
5.1. The VEMP as a filtered rate modulation
In the short-modulation approximation, the modulation is fully characterized by the quantity d, with the absolute value corresponding to the equivalent rectangular duration of the modulation and the sign distinguishing between excitation and inhibition. Since VEMP, v(t), and MUAP function, h(t), differ only by a constant factor (cf. Eq. (10)), there is an easy way to estimate d from experimental data. The basis for doing so is Eq. (11), where h(t) may be replaced by a properly rescaled v(t). If the rescaling is done so that the resulting h(t) has a maximum absolute value of one, where t^v is the time corresponding to that maximum, d may be estimated as
VEMPs are generally measured with the objective of evaluating the function of the vestibular system. Thus, in terms of the present model, it is the rate modulation r(t) that represents the effect of interest. But a VEMP measurement allows us to observe r(t) only indirectly. Figuratively speaking, recording the VEMP means that the rate modulation is observed through a filter with an impulse response corresponding to the MUAP function h(t) (see Eqs. (5) and (13)). The VEMP-associated variance modulation may be considered, accordingly, as the rate modulation observed through a filter with the impulse response h2(t) (see Eqs. (7) and (14)). A pair of equations that closely resembles our Eqs. (5) and (7) was obtained by Schoonhoven et al. (1988), although in a different context: they developed a method for the analysis of the compound action potential of a peripheral nerve (cf. their Eqs. (9) and (10)). It would certainly be desirable to be able to undo the filtering by applying some kind of inverse filter. If the function h(t) were known, r(t) could indeed be largely recovered by deconvolution, although the ill-posed nature of the inverse problem to be solved would require special provisions, such as optimal (Wiener)
dest ¼ Dmðt^v Þ
ð34Þ
The short-modulation approximation implies that a strong modulation of short duration may have exactly the same effect as a weaker modulation of longer duration. However, the results of this modeling study partially contradict that simple view, and so the question arises as to how closely dest matches the true value of d as calculated using Eq. (16). This question is investigated in Fig. 7. As in the previous figure, the parameter t was
5. Discussion
B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 286 (2011) 41–49
47
0.1
0
a
~ v(t)
Parameter r0
-0.2 -0.4
0
-0.6 -0.1
-0.8
0
-1
m(t)
2
δ (ms)
-2
3
-0.1
-20
-5 Time of Extremum (ms)
3.54 4
-3 -4
-10
0 time (ms)
10
20
Fig. 6. Varying the duration of the rate modulation while keeping the canonical VEMP invariant. The canonical VEMP shown in (a) may be associated with each of the variance modulations shown in (b). The numbers attached to the curves indicate the value of the duration parameter t.
8 t
6 -1
4
a
t
2 0
-2 δ (ms)
v(t )
Amplitude of Extremum
b
1
-1
-0.1
-3
m(t )
m(t )
-4
-0.2 2
4 τ (ms)
6
8
Fig. 5. Effect of the duration parameter t (same model as in Fig. 4). (a) The value r0 that is required to keep v~ 0 at the constant level 0.1 (cf. Eq. (22)). (b) Parameter d calculated using Eq. (16). (c) Times of the extrema of VEMP (black curve) and variance modulation (gray curve). (d) Amplitudes of the extrema of VEMP (black line) and variance modulation (gray curve). The dashed curve shows the variance modulation at the time of the VEMP extremum.
filtering (see e.g. Press et al., 2007). But h(t) is generally unknown, which vastly complicates the situation. At least partial information about h(t) could be easily extracted from the data. On the assumption that all MUAPs have the same waveform, the modulus of a Fourier-transformed MUAP, which would be 9H(f)9 in the notation of the present article, can be expected to have the same shape as the EMG spectrum (Wit and Kingma, 2006). Fig. 3 may help to understand what possibilities this opens up: given 9V~ ðf Þ9 (derived from the measured VEMP, v(t)) and an estimate of 9H(f)9, it is possible to estimate 9R(f)9, scaling issues left aside. While there is, of course, no unique way to reconstruct a time-domain function from the modulus of its Fourier transform, the above approach might provide valuable constraints. A different approach
30
relative estimation error (percent)
0
b
20
10
0 0
1
2
3
4
5
τ (ms) Fig. 7. Estimation of d using the short-modulation approximation (simulation for yt ¼5 ms and v~ 0 ¼ 0:1). (a) Comparison between true d (black curve) and estimate dest (dashed curve). (b) Relative estimation error for dest. The worst-case error is below 30%.
48
B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 286 (2011) 41–49
is suggested by the frequency-domain representation of our general model. Combining Eqs. (13) and (14) yields H2 ðf Þ V~ ðf Þ ¼ Hðf ÞMðf Þ
ð35Þ
which in the time-domain reads Z þ1 Z þ1 0 ~ h2 ðt 0 Þvðtt Þdt0 ¼ hðt 0 Þmðtt 0 Þdt0 1
ð36Þ
1
Thus, h(t) is formally the solution of the non-linear integral equation (36). However, this inverse problem can be assumed to be seriously ill-posed, which means that it is a challenging task to find a solution which not only explains the data, but is also reasonable from a physiological point of view. With this in mind, it is remarkable that the short-modulation approximation offers an amazingly simple way of escaping these difficulties. 5.2. The duration of the rate modulation matters The main impetus for this study was the concern that, in typical experiments, the duration of the rate modulation might be too long to be reconcilable with the assumption of a short rate modulation. Stimulating with clicks of 0.1 ms duration, Colebatch and Rothwell (2004) found a mean duration of 3.6 ms for the resulting inhibition of motor unit activity. If tone pulses are used rather than clicks, the inhibition might last even somewhat longer. The present analyses show that such durations are indeed not short in the sense of the model. In the simple Gaussian model considered here, the duration of the rate modulation is controlled by the parameter t. As to the VEMP, t has, in essence, a time-scaling effect: the VEMP has exactly the same waveform as the MUAP function h(t), apart from a t-dependent time dilation (see Fig. 4). Thus, supposing that h(t) is unknown, the VEMP gives no clue as to the value of t. Quite different is the situation in the case of the variance modulation, as best seen in Fig. 6: one and the same VEMP may be associated with fundamentally different variance modulations. This is clearly contradictory to a basic feature of the short-modulation approximation, namely that the variance modulation corresponds to the VEMP squared (apart from a scaling factor). As a matter of fact, a clear deviation from the short-modulation approximation (represented by the dotted curve in Fig. 6b) is found already for t ¼ 1 ms. 5.3. Making a case for the short-modulation approximation The fact that the duration of the rate modulation may well have relevant effects on both the VEMP and the associated variance modulation seems to render the short-modulation approximation invalid. However, such a conclusion would be premature. Without knowledge of the MUAP function h(t), most effects are not easily recognized as such. Thus, under realistic experimental conditions (i.e. in the case of noisy data), it may be difficult, if not impossible, to falsify the short-modulation approximation. This is particularly true for the VEMP. For the model considered in this study, the VEMP is basically a time-dilated version of h(t) so that, without using additional information or making assumptions, it is absolutely impossible to reach conclusions about the duration of the rate modulation. The assumption of an infinitesimal duration is, in that situation, just as arbitrary as any other assumption, but it has the inestimable advantage that the resulting model (here called ‘short-modulation approximation’) is extraordinarily simple from a mathematical point of view: getting the MUAP function h(t) merely requires rescaling of the VEMP. Thus, as to the latter, the short-modulation approximation is, by definition, in perfect agreement with the experimental data. The crucial question, consequently, is to what extent
the approximation allows us to explain the variance modulation. Fig. 6 provides an answer. The duration parameter t was systematically varied between 0 and 4 ms while the VEMP curve was kept constant (by concurrently changing other model parameters). At the times of the two VEMP peaks (dotted vertical lines in Fig. 6) and beyond (towards 7N), the effect of t on the variance modulation is relatively minor so that the differences between the curves may be considered to be insignificant (unless data with an outstanding signal-to-noise ratio are available). The short-modulation approximation consequently allows us to make reasonable predictions regarding the peak amplitudes, and for this reason it can be used to estimate the equivalent rectangular duration of the rate modulation, 9d9. The model calculations presented in Fig. 7b resulted in a worst-case estimation error of less than 30% (dashed curve). This appears to be quite acceptable, considering that much greater uncertainties are presumably introduced by the simplifying assumptions on which the model is based. 5.4. Limitations of the simple Gaussian model The simple Gaussian model may help to interpret experimental findings qualitatively, but it has not been devised to match measured data quantitatively. Apart from scaling factors, the model has only two parameters (the time constants y and t) so that the possibilities for simultaneously fitting the VEMP and the associated variance modulation are limited. Moreover, the choice of the functions h(t) and r(t) was motivated more by mathematical convenience than by physiological considerations. Thus, these functions are most likely suboptimal with regard to real conditions. On the other hand, given our current state knowledge, some arbitrariness in the definition of the two functions is unavoidable, and we are not aware of another simple model that would give us an obvious advantage over the one we investigated here.
6. Conclusions The simple Gaussian model introduced in this study is arguably one of the most straightforward VEMP models conceivable. It stands out in that all pertinent aspects can be examined analytically, often on the basis of surprisingly simple formulas. Despite the simplicity of the model, the predicted VEMP curves closely resemble those experimentally observed. Thus, all in all, the model appears to be ideally suited for getting a basic theoretical understanding of the VEMP phenomenon. This does not only apply to questions arising from previous work. Considering the fact that the model has a simple representation also in the frequency-domain, it might provide guidance for future work on VEMP deconvolution. According to the underlying general model, both the VEMP and the associated variance modulation may be considered as filtered versions of the effect of interest, which is represented by the temporal modulation of the MUAP rate, r(t). The main component of the VEMP experimentally observed, i.e. the p13–n23 complex, has roughly the same appearance as the VEMP predicted by our simple Gaussian model. According to this model, the time course of the VEMP mainly reflects the impulse response of the filter, which is the MUAP function h(t), whereas the rate modulation r(t) has basically a smoothing and broadening effect that depends on the duration of the rate modulation. The same applies to the variance modulation, where the filter corresponds to h2(t). But owing to the fact that this filter has a non-negative impulse response, a peculiarity arises: with increasing duration of the rate modulation, the variance modulation gradually loses its similarity
B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 286 (2011) 41–49
with the filter function h2(t) and finally looks more like the rate modulation r(t). The effect would, in principle, allow us to estimate the duration of the rate modulation. However, our ¨ ¨ experimental data (Lutkenh oner et al., 2011) point out complications, presumably due to the multi-component structure of the VEMP, which was disregarded in this study. The duration of the rate modulation does not affect all aspects of the VEMP to the same degree. A moderate duration increase affects mainly the amplitudes of VEMP and variance modulation, whereas the effect on the waveforms may be insignificant compared to the inevitable noise in the data. This study suggests that the effect on the peak amplitudes is well explained by the short-modulation approximation of the model, which offers an exceptionally simple and convenient way for a model-based interpretation of experimental results. References Boashash, B., 1992. Estimating and interpreting the instantaneous frequency of a signal. 1. Fundamentals. Proc. IEEE 80, 520–538. Brantberg, K., Verrecchia, L., 2009. Testing vestibular-evoked myogenic potentials with 90-dB clicks is effective in the diagnosis of superior canal dehiscence syndrome. Audiol. Neurootol. 14, 54–58. Brigham, E.O., 1974. The Fast Fourier Transform. Prentice-Hall, Englewood Cliffs, NJ. Colebatch, J.G., Halmagyi, G.M., Skuse, N.F., 1994. Myogenic potentials generated by a click-evoked vestibulocollic reflex. J. Neurol. Neurosurg. Psychiatry 57, 190–197.
49
Colebatch, J.G., Day, B.L., Bronstein, A.M., Davies, R.A., Gresty, M.A., Luxon, L.M., Rothwell, J.C., 1998. Vestibular hypersensitivity to clicks is characteristic of the Tullio phenomenon. J. Neurol. Neurosurg. Psychiatry 65, 670–678. Colebatch, J.G., Rothwell, J.C., 2004. Motor unit excitability changes mediating vestibulocollic reflexes in the sternocleidomastoid muscle. Clin. Neurophysiol. 115, 2567–2573. Day, S.J., Hulliger, M., 2001. Experimental simulation of cat electromyogram: evidence for algebraic summation of motor-unit action-potential trains. J. Neurophysiol. 86, 2144–2158. Gabor, D., 1946. Theory of communication. J. Inst. Electr. Eng. (London) 93, 429–457. Kingma, C.M., Wit, H.P., 2011. Asymmetric vestibular evoked myogenic potentials in unilateral Meniere patients. Eur. Arch. Otorhinolaryngol. 268, 57–61. ¨ ¨ Lutkenh oner, B., Stoll, W., Basel, T., 2010. Modeling the vestibular evoked myogenic potential. J. Theor. Biol. 263, 70–78. ¨ ¨ Lutkenh oner, B., Rudack, C., Basel, T., 2011. The variance modulation associated with the vestibular evoked myogenic potential. Clin. Neurophysiol 122, 1448–1456. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 2007. Numerical Recipes. The Art of Scientific Computing third ed. Cambridge University Press, Cambridge. Rosengren, S.M., Welgampola, M.S., Colebatch, J.G., 2010. Vestibular evoked myogenic potentials: past, present and future. Clin. Neurophysiol. 121, 636–651. Schoonhoven, R., Stegeman, D.F., van Oosterom, A., Dautzenberg, G.F., 1988. The inverse problem in electroneurography—I: conceptual basis and mathematical formulation. IEEE Trans. Biomed. Eng. 35, 769–777. Ushio, M., Iwasaki, S., Murofushi, T., Sugasawa, K., Chihara, Y., Fujimoto, C., Nakamura, M., Yamaguchi, T., Yamasoba, T., 2009. The diagnostic value of vestibular-evoked myogenic potential in patients with vestibular schwannoma. Clin. Neurophysiol. 120, 1149–1153. Wit, H.P., Kingma, C.M., 2006. A simple model for the generation of the vestibular evoked myogenic potential (VEMP). Clin. Neurophysiol. 117, 1354–1358.