Deconvolution of the vestibular evoked myogenic potential

Journal of Theoretical Biology 294 (2012) 87–97

Contents lists available at SciVerse ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

Deconvolution of the vestibular evoked myogenic potential n ¨ ¨ ¨ Bernd Lutkenh oner , Turker Basel

ENT Clinic, M¨ unster University Hospital, M¨ unster, Germany

a r t i c l e i n f o

abstract

Article history: Received 6 September 2011 Accepted 27 October 2011 Available online 6 November 2011

The vestibular evoked myogenic potential (VEMP) and the associated variance modulation can be understood by a convolution model. Two functions of time are incorporated into the model: the motor unit action potential (MUAP) of an average motor unit, and the temporal modulation of the MUAP rate of all contributing motor units, briefly called rate modulation. The latter is the function of interest, whereas the MUAP acts as a filter that distorts the information contained in the measured data. Here, it is shown how to recover the rate modulation by undoing the filtering using a deconvolution approach. The key aspects of our deconvolution algorithm are as follows: (1) the rate modulation is described in terms of just a few parameters; (2) the MUAP is calculated by Wiener deconvolution of the VEMP with the rate modulation; (3) the model parameters are optimized using a figure-of-merit function where the most important term quantifies the difference between measured and model-predicted variance modulation. The effectiveness of the algorithm is demonstrated with simulated data. An analysis of real data confirms the view that there are basically two components, which roughly correspond to the waves p13–n23 and n34–p44 of the VEMP. The rate modulation corresponding to the first, inhibitory component is much stronger than that corresponding to the second, excitatory component. But the latter is more extended so that the two modulations have almost the same equivalent rectangular duration. & 2011 Elsevier Ltd. All rights reserved.

Keywords: VEMP EMG Electromyogram Variance modulation Inverse filtering

1. Introduction The primary function of the otolith organs in the inner ear is to register changes in linear motion: while the saccule is sensitive to vertical accelerations (perceived, for example, in an elevator), the utricule responds to horizontal accelerations. Remarkably, the otolith organs also respond to high-level sounds (reviewed, e.g., by Halmagyi et al., 2005). This latter property is gaining increasing importance for functional testing in humans. What is exploited here is that loud sound pulses can elicit reflexes, which alter the tone of certain muscles. This opens up the possibility to study the function of an otolith organ by analyzing the electromyogram (EMG) of a suitably selected muscle (see, e.g., the recent review article by Rosengren et al., 2010). The data analysis is commonly confined to the stimulus-triggered average of the EMG (Colebatch et al., 1994), for which the term vestibular evoked myogenic potential (VEMP) has been coined. ¨ ¨ A recent modeling study (Lutkenh oner et al., 2010) showed that the VEMP is inevitably associated with a temporal modulation of the EMG variance. This prediction was subsequently confirmed by a reanalysis of archived data from 672 clinical VEMP investigations

n ¨ Correspondence to: ENT Clinic, Kardinal-von-Galen-Ring 10, 48129 Munster, Germany. 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.10.033

¨ ¨ (Lutkenh oner et al., 2011). The study also led to the conclusion that the variance modulation conveys information not being available in the VEMP itself. The data turned out to be more complex than anticipated in the simulations, though. Our interpretation was that VEMP and associated variance modulation comprise at least two components: a main component of inhibitory nature and at least one minor component of excitatory nature. But this interpretation remains speculative as long as the line of reasoning is purely qualitative. So the question arises how a quantitative understanding of the data can be accomplished. ¨ ¨ Our general model (Lutkenh oner et al., 2010) can readily handle any number of VEMP components, and so it should be suitable to explain the measured data. However, it is not easy to verify this supposition. Given detailed knowledge of the constituents of the model, calculation of the VEMP and the associated variance modulation is straightforward, but the problem to be solved is the other way around: what is given are the data, and the task is to come to conclusions as to the constituents of the model. It is typical for such an inverse problem to have many (maybe an infinite number of) solutions that fit the data about equally well. The dilemma can only be resolved by means of constraints and assumptions, which cause simple and plausible solutions to be favored over complicated and implausible ones. This aspect turned out to be critical also for the problem considered here. In what follows, we first recapitulate the model. Then we develop an algorithm for solving the inverse problem,

88

B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 294 (2012) 87–97

paying particular attention to constraints. The algorithm is finally tested using both simulated and measured data.

2. Theory 2.1. Model The model considered here is essentially identical with our ¨ ¨ previously developed general model (Lutkenh oner et al., 2010), which itself is based on the model by Wit and Kingma (2006). The notations that are used represent a minor modification1 of those used in our most recent VEMP-related article, in which we pre¨ ¨ sented a detailed analysis of a simple analytical model (Lutkenh oner and Basel, 2011). The model explains the VEMP and the associated variance modulation in terms of two time-dependent functions. The first function, h(t), reflects the time course of a motor unit action potential (MUAP). The model ignores the fact that, in the real world, the MUAP time course differs from motor unit to motor unit. Thus, h(t) has to be interpreted as representing a mean MUAP. In contrast to our previous work, the MUAP function is normalized such that Z þ1 2 h ðtÞ dt ¼ 1: ð1Þ 1

Thus, h(t) is measured in the units of s  1/2. The second function is the rate modulation r(t), which is dimensionless. This function describes the temporal variation of the aggregate discharge rate of the ensemble of motor units. To be more explicit, the aggregate discharge rate is assumed to be r0(1þr(t)), where r0 is the discharge rate without modulation. ¨ ¨ According to the model (Lutkenh oner et al., 2010), the VEMP is proportional to the convolution of the functions r(t) and h(t), i.e. Z þ1 vðtÞ ¼ Z rðtÞhðtt 0 Þdt0 , ð2Þ 1

whereas the associated variance modulation corresponds to the convolution of r(t) and h2(t), i.e., Z þ1 2 rðt 0 Þh ðtt 0 Þdt 0 : ð3Þ mðtÞ ¼ 1

A similar pair of integral equations was obtained by Schoonhoven et al. (1988), although in a different context. In this article, v(t) always refers to the normalized VEMP, which may be understood as the VEMP derived from a normalized EMG (with unit variance). As to our model this means that the proportionality constant in Eq. (2) is

Z¼

a pffiffiffiffiffiffi r0 , a

ð4Þ

where a=a is a factor depending on the distribution of the MUAP amplitudes (ratio of arithmetic and quadratic pffiffiffiffiffiffiffiffiffi mean). The factor ¨ ¨ may be expected to be of the order of 2=3  0:82 (Lutkenh oner et al., 2010) so that Z can, to a first approximation, be interpreted pffiffiffiffiffiffi as r0 . 2.2. Deconvolution algorithm The primary question addressed in this article is how to find the model that best fits given data. Our answer is schematically depicted in Fig. 1. The goal of the algorithm illustrated in this figure is to determine r(t) and h(t) in such a way that the functions v(t) and m(t), calculated using Eqs. (2) and (3), optimally explain their experimentally obtained counterparts, vðtÞ and mðtÞ. The 1 ¨ ¨ The model formulated in the p present paper can be derived from Lutkenh oner ffiffiffiffi and Basel (2011) by denoting hðtÞ= D as h0 (t). The prime is omitted in this article, for the sake of simplicity. Besides, vnarm(t) is simply called v(t).

function playing the pivotal role in this algorithm is r(t). In the present subsection we act as though r(t) were known, apart from an unknown scaling factor c. This means rðtÞ ¼ cr^ ðtÞ, where r^ ðtÞ is a given function. The iterative optimization of r^ ðtÞ will be considered in a later subsection. ~  c Z hðtÞ by Given r^ ðtÞ, Eq. (2) allows us to determine hðtÞ deconvolution. The deconvolution is conveniently done in the frequency domain using the equation ^ ÞHðf ~ Þ, Vðf Þ ¼ Rðf

ð5Þ

^ Þ, and Hðf ~ Þ are the Fourier transforms of vðtÞ, r^ ðtÞ, where Vðf Þ, Rðf ~ ^ Þ ¼ 0, the and hðtÞ, respectively. Except for frequencies with Rðf ~ Þ¼ straightforward solution to the deconvolution problem is Hðf ^ Vðf Þ=Rðf Þ. However, without taking special measures, noise would ^ Þ9 is relatively be enormously amplified at frequencies where 9Rðf small. Avoiding this problem requires some kind of regularization: ^ Þ, Vðf Þ is multiplied by a function instead of being divided by Rðf inv ^ Þ if 9Rðf ^ Þ9 is relatively large, R^ ðf Þ, the value of which is about 1=Rðf ^ and about zero if 9Rðf Þ9 is relatively small. A function fulfilling these requirements is inv 1 1 R^ ðf Þ ¼ : ^ Þ 1 þ e2 =9Rðf ^ Þ92 Rðf

ð6Þ

The inherent regularization may be understood as optimal (Wiener) filtering. The latter method (see, e.g., Press et al., 2007) is based on the idea that the measured data are composed of signal s(t) and noise n(t). The signal is assumed to be a blurred version of an uncorrupted signal u(t); the blurring is modeled as the convolution of u(t) with some response function r(t). Estimating the uncorrupted signal from the measured data evidently requires a deconvolution. The purpose of optimal filtering is to reduce the effect of noise so that the result of the deconvolution is as close as possible to the uncorrupted signal. If S(f) and N(f) are the Fourier transforms of signal and noise, respectively, the filter that is optimal in the least-squares sense is   Sðf Þ2 1 Fðf Þ ¼ : ð7Þ  2 ¼ 2 2 2   9Sðf Þ9 þ Nðf Þ 1 þ 9Nðf Þ9 =9Sðf Þ9 In the present context, the signal and the uncorrupted signal are represented by VEMP and MUAP, respectively, whereas the rate modulation is the response function that causes the blurring. Thus, S(f) corresponds to the right-hand side of Eq. (5). As to the noise, it appears reasonable to assume that background EMG activity is dominating, and since the latter is a superposition of MUAPs, the frequency spectrum of the EMG may be assumed to be proportional to the spectrum of the MUAP (Wit and Kingma, 2 ~ Þ92 , where e2 is some propor2006). This means, 9Nðf Þ9  e2 9Hðf tionality factor. So we finally obtain

Fðf Þ ¼

1 1 : ¼ ~ Þ92 =9Rðf ^ Þ92 9Hðf ~ Þ92 Þ ^ Þ92 1þ ðe2 9Hðf 1 þ e2 =9Rðf

ð8Þ

In conclusion, the second factor on the right-hand side of Eq. (6) represents the optimal (Wiener) filter. ~ ~ Þ, an inverse Fourier transform provides hðtÞ. Given Hðf Evalua~ tion of the integral in Eq. (1) with h(t) replaced by hðtÞ yields c2Z2. ~ To obtain h(t), the function pffiffiffiffiffiffiffiffiffiffi hðtÞ has to be divided by cZ, which corresponds to sgnðcÞ c2 Z2 (the factor Z, introduced in Eq. (4), is conceived as a positive value). Since the sign of parameter c, sgn(c), is unknown at this point, we tentatively assume a positive sign. This is uncritical as sgn(c) is irrelevant for now, and a wrong decision can be easily corrected later. Given r^ ðtÞ and the estimate obtained for h(t), the scaling parameter c can be determined from ^ the measured variance modulation, mðtÞ. If mðtÞ is the variance modulation calculated using Eq. (3) with r(t) replaced by r^ ðtÞ, the

B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 294 (2012) 87–97

89

Fig. 1. Scheme illustrating key aspects of the deconvolution algorithm. Time-domain operations are shown on the left, frequency-domain operations on the right. The algorithm is started by making an initial guess for parameters defining an arbitrarily scaled rate modulation function, r^ ðtÞ. The corresponding MUAP function, h(t), is then obtained by a deconvolution of the VEMP (in the frequency domain). Given both the rate modulation and the MUAP function, a model prediction for the variance ^ modulation, mðtÞ, can be calculated. A comparison with the variance modulation actually measured, mðtÞ, allows us to rescale the rate modulation using an optimized scaling factor, c^ . Data and modeling results are eventually used to evaluate a figure-of-merit function. The result is fed to a nonlinear optimization algorithm, which provides new parameters for the function r^ ðtÞ. The process is iteratively repeated until certain exit criteria are met. For more details see text.

2.3. A deconvolution example

least-squares estimate of c is c^ ¼

Z

þ1

^ mðtÞmðtÞdt 1

Z

þ1

^ 2 ðtÞdt: m

ð9Þ

1

Given this estimate and the value previously obtained for the normalization factor c2Z2, the parameter Z can be determined. Moreover, rescaling of r^ ðtÞ yields r(t). If c^ turns out to be negative, h(t) is multiplied by  1 now, correcting the wrong decision above. The model predictions for VEMP and associated variance modulation can finally be calculated using Eqs. (2) and (3).

The above algorithm is now illustrated with an example. The rate modulation shown in Fig. 2a reflects the hypothesis that the VEMP arises from two modulations of opposite polarity: an inhibitory modulation causing the VEMP peaks p13 and n23, and an excitatory modulation causing the VEMP peaks n34 and p44. The time of the inhibitory modulation was arbitrarily defined as zero; the excitatory modulation was assumed to occur 21 ms later, corresponding to the latency difference of the VEMP peaks p13 and n34 (or n23 and p44). The disparity in peak amplitude ( 0.9 and 0.15) reflects the fact

90

B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 294 (2012) 87–97

10-3 1 / Hz

Rate Modulation

0

10-4

-0.5

00 80 1 / 000 4 0 1/ 00 / 2 00 0 /1

103

0

1

RM convolved with inverse RM (1/s)

50000

1

1

Hz

Inverse Rate Modulation (1/s2)

10-5

-50000

1

102

/2

/5

50

1

-100000

00

/1

25

101

1 200

10-1

100

0

10-2 0

50

0

Time (ms)

100 200 Frequency (Hz)

300

Fig. 2. Illustration of the concept of an inverse rate modulation. (a) Example of a rate modulation. A strong inhibitory modulation is followed by a much weaker excitatory modulation. (b) Absolute value of the Fourier transformed rate modulation. (c) Inverse rate modulation. The value of the regularization parameter, e, was 1/1000 Hz  1. (d) Absolute value of the Fourier transform of the inverse rate modulation. The solid curve corresponds to the time domain function shown on the left; the dotted curves were obtained for other values of the regularization parameter e. (e) The rate modulation convolved with the inverse rate modulation. Without regularization, the curve would correspond to a Dirac delta function. (f) The curves in (d) multiplied by the curve in (b). Without regularization, the multiplication would yield the value 1, for all frequencies.

that the first two VEMP peaks are generally the dominating ones. Both modulations were assumed to have a Gaussian shape. Thus, the rate modulation was defined as   ! K X 1 tt k 2 r k exp  r^ ðtÞ ¼ ð10Þ 2 tk k¼1

As the Fourier integral is a time integral and r^ ðtÞ is dimensionless, ^ Þ is measured in units of time, which corresponds to the Rðf ^ Þ9 is shown in Fig. 2b. reciprocal of frequency. The graph of 9Rðf While the spectral profile of the rate modulation is fairly flat at low frequencies, there is a rapid decrease at higher frequencies.

(with K¼2 in the present example). The parameter tk determines the peak width of the kth component. The values assumed in the present example were 2 ms and 3 ms. Analytical evaluation of the Fourier transform of r^ ðtÞ yields2

inv ization parameter e. Fig. 2d shows 9R^ ðf Þ9 for 7 different values of 1 e, ranging from 1/125 Hz to 1/8000 Hz  1. Without regularization (i.e., e ¼0), the function

K pffiffiffiffiffiffi X ^ Þ ¼ 2p tk rk exp ð2p2 t2k f 2 i2pt k f Þ: Rðf

ð11Þ

k¼1

2 ¨ ¨ The case K¼1 was considered by Lutkenh oner and Basel (2011). If R(f) is the Fourier transform of r(t), the Fourier transform of the time-shifted function rðtt 0 Þ is exp ði2pt 0 f ÞRðf Þ (see, e.g., Brigham, 1974).

inv The function R^ ðf Þ, defined in Eq. (6), depends on the regular-

^ ÞR^ inv ðf Þ Uðf Þ ¼ Rðf

ð12Þ

would have a value of unity for all frequencies. Thus, the graph of 9U(f)9 allows us to assess deviations from this ideal condition (Fig. 2f). With e ¼1/1000 Hz  1 (solid curve), the deviations appear to be acceptable up to an upper limit of about 100 Hz. The limit can be increased by decreasing e, but the consequence is a potentially problematic boost effect at frequencies of about 200 Hz (cf. Fig. 2d).

B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 294 (2012) 87–97

If e is significantly greater than 1/1000 Hz  1, even low-frequency components are notably underestimated, which is clearly not acceptable. All in all, assigning a value of about 1/1000 Hz  1 to e seems to be a good compromise here. An inverse Fourier transform of U(f) yields the function u(t), which is shown in Fig. 2e for the case e ¼1/1000 Hz  1. As this step requires the calculation of a frequency integral and U(f) is dimensionless, u(t) is measured in units of frequency, which corresponds to the reciprocal of time. The more pulse-like this function is, the higher, in theory, the fidelity of the deconvolution result is. Thus, u(t) can be used for assessing the performance of the deconvolution procedure. The two measures considered in this article are the area and the width of the peak centered at time zero (filled area in Fig. 2e). In the present example, the peak width (at the base) is 7.2 ms and the area is 1.048, which is close to the ideal value of unity. The above chain of thoughts began in the time domain (Fig. 2a), continued in the frequency domain (Fig. 2b, d, and f), and finally returned to the time domain (Fig. 2e). The inverse Fourier transform inv inv of R^ ðf Þ, i.e., r^ ðtÞ, was left out of account so far. To complete this inv example, r^ ðtÞ is shown in Fig. 2c for e ¼1/1000 Hz  1; the function shall be called inverse rate modulation (although the term pseudoinverse would be more appropriate). The curve in Fig. 2e could now be calculated also in the time domain: by the convolution of the rate modulation (Fig. 2a) and the inverse rate modulation (Fig. 2c).

So far, the function r^ ðtÞ was assumed to be known. But as a matter of fact, only a rough guess is available, at best. Thus, the deconvolution problem considered above is part of a more comprehensive problem, which also comprises the optimization of r^ ðtÞ. The fundamental question arises as to how to define the figure-of-merit function that scores the match between model and data. A figure-of-merit function consisting of several terms is needed here. Most important is the term that measures the difference between the variance modulation predicted by the model, m(t), and the variance modulation actually measured, mðtÞ. A suitable measure of this difference is Z þ 1 Z þ1 ðmðtÞmðtÞÞ2 dt m2 ðtÞdt: ð13Þ Qm ¼ 1

The difference between model-predicted VEMP, v(t), and measured VEMP, vðtÞ, is accordingly quantified as Z þ 1 Z þ1 ðvðtÞvðtÞÞ2 dt v2 ðtÞdt: ð14Þ Qv ¼ 1

1

The latter term is generally of marginal importance, though. The point is that the MUAP is estimated by deconvolution of the measured VEMP. Thus, convolution of the estimated MUAP and the rate modulation is, to a first approximation, nothing else than the reversal of the deconvolution. This explains why modelpredicted and measured VEMP are usually highly concordant. The error term Qv is nonetheless indispensable, because it may become relevant when the regularization is problematic or inadequate. A mathematical model explaining measured data is not necessarily reconcilable with what we know about the nature of the investigated phenomenon. To counteract this potential problem, penalty terms are added to the two error terms. The first penalty term is Q1 ¼

Z

þ1

exp 1



 9rðtÞ91

lr

dt:

The exponential function serves as a barrier function: with lr ¼ 0.005, the value of the integrand rapidly increases when 9r(t)9 grows beyond one. Thus, the penalty term Q1 can be used to impose a ‘‘soft’’ upper limit on the value of 9r(t)9. A second penalty term, defined as Z þ 1 Z tb 2 2 h ðtÞdt h ðtÞdt, ð16Þ Q 2 ¼ 1 ta

1

aims to achieve that the MUAP function is temporally compact, which means that the function is basically confined to a limited time interval (the interval [ta,tb] in terms of the above equation). If the goal is roughly accomplished, Q2 is about zero; any activity outside the specified interval causes Q2 to increase. In this study, [ta,tb] corresponded to the time range 5 ms to 35 ms. Additional penalty terms may be required to ensure other desirable properties of r^ ðtÞ. If these terms are combined into a single term Q3, the complete figure-of-merit function may be written as Q ¼ wv Q v þwm Q m þ w1 Q 1 þw2 Q 2 þw3 Q 3

ð17Þ

where wv, wm, w1, w2, and w3 are the weighting factors. The latter were set to one in this study, for the sake of simplicity. Although not explicitly stated, Q depends on the parameters of r^ ðtÞ. The latter have to be adjusted so that Q is minimized. Numerous methods for solving this kind of problem are readily available. The method used here will be described further below. 2.5. Mathematical representation of the rate modulation

2.4. Iterative parameter optimization

1

91

ð15Þ

The rate modulation function could be set up in numerous ways. Here, the function defined in Eq. (10) was used. If the number of components, K, is small, there may be no need to constrain the parameters of this function by means of a specific penalty term. However, the situation changes as K increases: it may happen then that the estimated r^ ðtÞ shows a much higher complexity than is supported by the data, which means that some components basically reflect noise. To counteract this problem, a smoothness constraint may be helpful. For the function considered here, smoothness can be enforced by requesting that the values of the parameters tk (1rkrK) do not fall below a certain limit, tmin. Another useful constraint is that the times tk (1rkrK) have a certain lower limit. Both constraints can be implemented using barrier functions similar to Eq. (15). This leads to the penalty term:   X   K K X t t t t Q3 ¼ exp min k þ exp min k ð18Þ k¼2

lt

k¼1

lt

The parameter values used in this study were t min ¼ tmin ¼ 1 ms, and lt ¼0.004 ms. Note that the penalty term Q3 does not impose a constraint on t1. The reason is that t1 is one of the two parameters that are predefined. The optimization problem is, as a matter of fact, overdetermined, because even under ideal conditions there are two parameters which cannot be determined from the data. One aspect requiring special attention is that it is impossible to determine the timing of the functions r(t) and h(t) in a unique way: supposing that r(t) is a solution to the problem, r(t  t0) is a solution as well, because any time shift of r(t) can be compensated for by shifting h(t) in the opposite direction. Thus, the shift parameter t0 has to be defined more or less arbitrarily. A second aspect to be considered is that r^ ðtÞ is, conceptually, a function with arbitrary scaling (the estimation of an optimal scaling factor is considered in Eq. (9)). Thus, one of the parameters rk (1rkrK) can be arbitrarily set to any value except zero. This indeterminacy is resolved by arbitrarily setting t1 ¼0 and r1 ¼ 1; the sign of r1 expresses the expectation that the first (main) VEMP component is inhibitory.

92

B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 294 (2012) 87–97

2.6. Numerical calculations ¨ ¨ To be consistent with our previous publications (Lutkenh oner ¨ ¨ et al., 2010, 2011; Lutkenh oner and Basel, 2011), the above theory was formulated in continuous time. But the numerical calculations presented here were, of course, done in discrete time (using Matlab 6.5, The MathWorks Inc., Natick, MA, USA). This means, in particular, that a fast Fourier transform (FFT) was used to switch from the time to the frequency domain, and an inverse FFT was used for the reverse switch. The results of these operations were always rescaled in such a way that they conform to continuous-time theory. In the case of the Matlab routines FFT and IFFT, this requires multiplication with the sampling interval (0.1 ms) and the sampling rate (10,000 Hz), respectively. Prior to Fourier transformation, the data (time range 40 ms to 80 ms) were windowed using a cosine-tapered window (Matlab routine TUKEYWIN); the ratio of taper to constant sections was 0.5. The resulting data vector was padded with trailing zeros to length 213. As to the integrals in Eq. (9), the calculations accounted for the times between 10 ms and 60 ms. The minimization of the figure-of-merit function Q was accomplished using Matlab routine FMINSEARCH, which is an implementation of the simplex search method of Lagarias et al. (1998). Starting at an initial guess, optimal parameters are searched without using analytical or numerical gradients. The algorithm does not necessarily find the global minimum so that the choice of the initial guess requires careful attention. With this in mind, the figure-of-merit function, Q, was evaluated for thousands of trial parameter sets and the parameter set yielding the smallest Q became the initial guess. If the number of components was small (Kr3), the trial parameter sets were obtained using a course discretization of the parameter space. If there were more components, random points of the parameter space were chosen. If the rate modulation comprises just a few components, its specific mathematical representation (see above) generally ensures that the main peak of the estimated function is centered at time zero. However, the built-in standardization may fail if there are more components than necessary to explain the data (an analysis with K ¼20 will be presented below). Time shifting the estimated rate modulation can, ex post facto, fix the problem.

both the VEMP and the variance modulation. The noise reduction at the beginning and at the end of the curves is caused by the tapering that preceded the deconvolution. The question now is how well rate modulation and MUAP function can be reconstructed from these noisy data. The answer is provided by the black curves in Fig. 3c and d. A comparison with the gray curves shows that the deconvolution algorithm worked out well. Apart from noise, the estimated curves distinctly differ only in a single respect from the model underlying the synthetic data: the estimated rate modulation (Fig. 3c) shows a peak at time zero that is smaller in amplitude and somewhat broader than assumed in the model (the estimated r1. is  0.74 instead of 0.9; the estimated t1 is 2.49 instead of 2). But this is not a serious discrepancy. On the contrary, such an estimation error was to be expected. Theoretical considerations ¨ ¨ (Lutkenh oner et al., 2010) led to the conclusion that the effect of a brief inhibition can be characterized by the equivalent rectangular duration (ERD) of the inhibition, regardless of the time course of the inhibition. This implies that measured data with typical background noise do not allow us to distinguish a strong inhibition of short duration from a weaker inhibition of longer duration. For a Gaussian peak as in Eq. (10), with amplitude rk and a width characterized by the parameter tk, the ERD is pffiffiffiffiffiffi dk ¼ 2p9r k 9tk ð19Þ ¨ (cf. Lutkenh¨ oner and Basel, 2011). Applying this formula to both the model underlying the simulated data and the estimate yields ERDs of 4.51 ms and 4.59 ms, respectively. Thus, model and estimate agree almost perfectly in this respect. It remains to be said that the regularization parameter e had the value 1/1000 Hz  1, as was the case for the solid curve in Fig. 2d. The deconvolution algorithm essentially operates with Fourier transformed data, and so it appears appropriate to compare model and estimate also in the frequency domain. This is done in the right half of Fig. 3. The differences between exact model and estimate seem to be more pronounced than in the time domain, but this concerns mainly higher frequencies, where the signal-to-noise ratio is bad. Besides, the use of a logarithmic scale emphasizes errors that are small in absolute terms. 3.2. Deconvolution under more realistic test conditions

3. Results 3.1. Deconvolution of simulated data The performance of the algorithm developed above is investigated now with a simple example. The rate modulation function shown as a gray curve in Fig. 3c is identical with the one considered in Fig. 2a. The MUAP function, represented by the gray curve in Fig. 3d, corresponds to the first derivative of a Gaussian, as suggested by Wit and Kingma (2006) and also assumed in our ¨ ¨ analytical model (Lutkenh oner and Basel, 2011).3 Adding white noise to the convolution of rate modulation and MUAP function yielded the VEMP that is shown as a gray curve in Fig. 3a. The associated variance modulation, again with white noise added, is represented by the gray curve in Fig. 3b. The standard deviation of the noise amounts to 10% of the maximal signal amplitude, for 3 The parameters of the rate modulation function (cf. Eq. (10)) are r1 ¼  0.9, r2 ¼0.15, t1 ¼0 ms, t2 ¼ 21 ms, t1 ¼2 ms, and t2 ¼3 ms. The time constant y of the pffiffiffiffiffiffi ¨ ¨ MUAP function (cf. Eq. (17) of Lutkenh oner and Basel, 2011) is 21  4:58 ms. The time interval between the two main peaks of the simulated VEMP is then exactly 10 ms, supposing that there is no noise in the data (the quantity yt in Eq. (21) of ¨ ¨ Lutkenh oner and Basel (2011) refers to one half of that interval). The value of the parameter Z was 25 s  1/2.

The estimation errors in the above example were caused, on the one hand, by noise in the simulated data and, on the other hand, by the regularization incorporated into the deconvolution. When analyzing real data, another source of error must be taken into account: a model is at best a rough approximation of the reality. As to that, the conditions considered in Fig. 3 were unrealistic, because the synthetic data were generated with a rate modulation that was a specific realization of the rate modulation function implemented into the deconvolution algorithm (Eq. (10)). In Fig. 4 it is investigated how the deconvolution algorithm performs under less ideal conditions. While the MUAP function remained unchanged (gray curve in Fig. 4d), the rate modulation is a staircase function now (gray curve in Fig. 4c). This means that the rate modulation is constant for certain periods of time and changes instantaneously from one moment to the next. An inhibitory modulation is followed by two non-overlapping excitatory modulations of smaller amplitude.4 The inhibitory component has 4 The inhibitory modulation was centered at time zero and had the amplitude  0.9. The two excitatory modulations were centered at 6 ms and 21 ms, and the amplitudes 0.4 and 0.2, respectively. The durations of the three modulations pffiffiffiffiffiffiwere pffiffiffiffiffiffi pffiffiffiffiffiffi were 2 2p, 4 2p, and 5 2p ms. As to the ERD, the modulations correspond to Gaussian peaks with time constants of 2, 4, and 5 ms, given that the peak amplitudes are the same (cf. Eq. (19)).

B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 294 (2012) 87–97

10-2

1 / Hz

normalized VEMP

1

93

0

10-3

10-4

10-3

0 1 / Hz

Variance Modulation

-1

10-4

-0.2 10-5

10-3 1 / Hz

Rate Modulation

0

-0.5

10-4

10-5

10-1

10-2 Hz-1/2

MUAP (s-1/2)

5 0

10-3

-5 -10 0

50 Time (ms)

0

100 200 Frequency (Hz)

300

Fig. 3. Deconvolution of simulated data (first example). (a) and (b) VEMP and associated variance modulation, respectively. The gray curves represent the simulated data; the black curves represent the model prediction that was made on the basis of the deconvolution result. (c) and (d) Rate modulation and MUAP, respectively. The black curves represent the deconvolution result; the gray curves represent the model that was used for generating the data. (e) to (h) Frequency-domain representation (absolute value of the Fourier transform) of the curves on the left.

exactly the same ERD as in the previous example. VEMP and associated variance modulation, synthesized as in the previous example, are shown as gray curves in Fig. 4a and b. A first analysis is based on the correct assumption that the rate modulation comprises three components. The results of this analysis are represented by the solid black curves in Fig. 4. Both

the VEMP (Fig. 4a) and the associated variance modulation (Fig. 4b) are explained well, even though the estimated rate modulation (Fig. 4c) is only a coarse approximation of the staircase function underlying the synthesized data. The estimated MUAP function (Fig. 4d) agrees almost perfectly with the function underlying the synthesized data, which is evidence in support of

B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 294 (2012) 87–97

normalized VEMP

94

basically equivalent to their rectangular counterparts in the rate modulation underlying the synthesized data. In general, it is not a priori known how many components the rate modulation has. The dotted curves in Fig. 4 illustrate the consequences of falsely assuming only two components. While this wrong assumption does not compromise the ability of the model to explain the VEMP (Fig. 4a), the optimal fit to the variance modulation is clearly deficient (Fig. 4b). This suggests that the answer to the question as to how many components are needed (supposing that there are just a few of them) may be found in the data themselves.

0

-1

Variance Modulation

3.3. Deconvolution of measured data 0.2

0

-0.2

Rate Modulation

0.5

0

-0.5

MUAP (s-1/2)

10 5 0 -5 -10 0 Time (ms)

50

Fig. 4. Deconvolution of simulated data (second example). The figure is organized in the same way as the left half of Fig. 3. The main difference between the two examples is that the rate modulation is a staircase function now, comprising three rather than two components (gray curve in (c)). The deconvolution algorithm cannot faithfully recover this rate modulation from the data, because the build-in rate modulation is composed of Gaussian peaks. But the estimated function shows at least the correct structure (black curve in (c)). Moreover, the estimated MUAP (black curve in (d)) is basically a noisy replica of the true MUAP (gray curve in (d)). The model prediction made on the basis of the deconvolution result (black curves in (a) and (b)) explains both the VEMP and the variance modulation (gray curves in (a) and (b)) fairly well. Thus, considering the noisy nature of the simulated data, the goodness of fit gives no reason to question the assumption of Gaussian shaped rate modulations. This example demonstrates that noise in the data severely limits the possibility to recover temporal details of the rate modulation. Nevertheless, grossly wrong assumptions as to its general structure do not go unnoticed: the dotted curves show a solution based on the false assumption that the rate modulation consists of only two components; the fit of the variance modulation is clearly deficient now.

¨ ¨ the idea (Lutkenh oner et al., 2010) that the effect of a brief inhibition (or excitation) can be characterized by its ERD. From that point of view, the Gaussian peaks in the estimated rate modulation are

The example considered in Fig. 4 was designed so that VEMP and variance modulation bear some similarities with our recently ¨ ¨ published experimental data (Lutkenh oner et al., 2011). But there are also conspicuous differences. For example, the two main VEMP peaks show a substantial asymmetry in the synthetic data, whereas they are roughly congruent in the measured data, apart from their polarity. All our attempts to achieve a significantly better match by manually adjusting the model parameters failed. Eventually, these failures gave the impetus for developing the deconvolution algorithm presented in this article. So the fundamental question now is how the algorithm copes with real data. ¨ ¨ To answer this question, the grand-averaged data of Lutkenh oner et al. (2011) are analyzed. VEMP and associated variance modulation are shown as gray curves in Fig. 5a and b, respectively. Solid black curves represent an analysis based on the assumption that the rate modulation function comprises two components. The agreement between model prediction and data can be considered fairly good. The rate modulation determined by the algorithm (solid curve in Fig. 5c) is reminiscent of the rate modulation assumed in Fig. 3, in that there is an inhibitory, main component followed by an excitatory, minor component. Table 1 provides not only the estimated model parameters, but also the ERDs of the two components, calculated using Eq. (19). While the second component is of much smaller amplitude, it has a longer duration than the first component so that the ERDs are similar. The estimated MUAP (solid curve in Fig. 5d) resembles the MUAP assumed in the above simulations (Figs. 3d and 4d). However, the amplitudes of the two dominant peaks clearly differ. It remains to be said that the value obtained for parameter Z was 26.3 s  1/2. According to Eq. (4), this to pcorresponds ffiffiffiffiffiffiffiffiffi a mean MUAP rate, r0, of 1035/s (assuming a=a ¼ 2=3)—in good ¨ ¨ agreement with a previous analysis of the same data (Lutkenh oner et al., 2011), which resulted in the value 1150 s  1. The assumption of a two-component rate modulation led to such a good agreement between data and model that there is little room for further improvement. In other words, the data give no clear indication of additional components. This does, of course, not exclude the possibility that the time course of the true rate modulation is more complex than suggested by the two-component solution. The dotted curve in Fig. 5c gives an idea of how the solution can look like if the model has more degrees of freedom. The rate modulation was allowed to consist of 20 components. Almost 60 parameters can be independently varied in this case so that there is a vast diversity of possible time courses. With that said, the optimal solution found by the deconvolution algorithm is surprisingly simple: it is, in essence, nothing but a modification of the two-component solution; the main difference is that the excitatory peak appears to be dichotomous. Considering the similarity with the two-component solution, it is not surprising that a conspicuous superiority of the more complex solution is found neither for the VEMP (Fig. 5a) nor the variance modulation (Fig. 5b). Fig. 5d shows a remarkable effect, though: while the MUAP related to the two-component solution (solid curve) shows small oscillations that extend beyond 50 ms, the MUAP related

0.2

95

10-3

0.1 1 / Hz

normalized VEMP

B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 294 (2012) 87–97

0

10-4

-0.1

10-5

-0.2

0

10-4

1 / Hz

Variance Modulation

10-3

10-5

-0.1

0 10-4

1 / Hz

Rate Modulation

10-3

-0.1

10-5 -0.2

10-2

0

Hz-1/2

MUAP (s-1/2)

5

-5 -10

10-3

10-4 0 Time (ms)

50

0

100 200 Frequency (Hz)

300

¨ ¨ Fig. 5. Deconvolution of grand-averaged experimental data (Lutkenh oner et al., 2011). The figure is organized in the same way as Fig. 3. The gray curves (only found in (a) and (b) and (e) and (f)) represent the data; the black curves represent a deconvolution based on the assumption that the rate modulation consists of two components (black curve in (c)). The model explains the data so well that adding more components to the rate modulation does not result in a major improvement. This is illustrated by the dotted curves, which represent a solution where the deconvolution algorithm was allowed to independently adjust 20 components. The data apparently do not require so many degrees of freedom: the estimated rate modulation (dotted curve in (c)) is, in essence, only a modification of the two-component solution, with a dichotomous excitatory component. Table 1 Parameters of the two-component solution considered in Figs. 5 and 6. The bottom row provides, furthermore, the equivalent rectangular duration (ERD) of each component. Component number, k

Parameter rk Parameter tk (ms) Parameter tk (ms) ERD (ms)

1

2

 0.225 0 1.56 0.877

0.048 17.23 6.17 0.749

to the more complex solution (dotted curve) basically fades away within about 30 ms. Thus, when using the more complex rate modulation, the deconvolution algorithm appears to use the abundance of

additional parameters mainly for minimizing the penalty term Q2 (cf. Eq. (16)), which serves to confine the MUAP to a pre-defined time interval. The right half of Fig. 5 is analogous to the right half of Fig. 3 and provides a frequency-domain representation of the curves in the left half. The good match between data and model becomes apparent also here. Again, there is nothing to suggest that the more complex solution (dotted curves) is at a major advantage over the two-component solution (solid curves). A comparison with the model simulation in Fig. 3 shows that the real data exhibit a less rapid decay at higher frequencies. To deepen the understanding of the two-component solution, the components will be considered separately now. The VEMP is studied in Fig. 6a, the associated variance modulation in Fig. 6b. The first component is represented by a solid black curve, the second component by a dotted curve; a gray curve shows the experimental

96

B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 294 (2012) 87–97

namely that certain aspects of the VEMP can be understood using a simple approximation. We came to the conclusion that this approximation – relying on the assumption that the VEMP is caused by a rate modulation of short duration – ‘‘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.’’ We took the second part of that statement as a challenge, and in this third article we now present a working solution to the inverse problem. Our deconvolution algorithm optimizes parameters that determine the rate modulation r(t). Thus, this function plays a pivotal role. The MUAP h(t), by contrast, is not described in terms of parameters, but is obtained by a deconvolution of the measured VEMP with the rate modulation. There are, of course, alternatives to this approach. So it would be perfectly feasible to let the MUAP h(t) play a pivotal role. At first glance, this approach may even appear to be more elegant, because Eqs. (2) and (3) can be combined into the single equation: Z þ1 Z þ1 2 h ðt 0 Þvðtt 0 Þdt 0 ¼ Z hðt 0 Þmðtt 0 Þdt 0 ð20Þ

normalized VEMP

0.2 0.1 0 -0.1

Variance Modulation

-0.2

0

1

-0.1 0

50 Time (ms)

Fig. 6. More detailed consideration of the two-component solution of the previous figure. The contribution of the first component is represented by a black solid curve, the contribution of the second component by a dotted curve; the experimental data are again represented by a gray curve: (a) VEMP and (b) variance modulation.

data. Only the first component is relevant in the time range of VEMP peak p13, while the second component dominates in the time range of VEMP peak p44. In between, the data interpretation is more complicated because both components have to be taken into account. The most striking discrepancy between the two components – apart from the differences in timing and amplitude – is found in the time course of the variance modulation: while the first component shows two separate peaks that basically coincide with the VEMP peaks p13 and n23, the second component shows a single broad peak that roughly coincides with a zero crossing in the VEMP curve (between the peaks n34 and p44). According to our analytical ¨ ¨ model (Lutkenh oner and Basel, 2011), the discrepancy can be attributed to the durations of the two rate modulations: for a component of relatively short duration (first component), a variance modulation resembling the squared VEMP is to be expected, whereas for a component of longer duration (second component) the predicted variance modulation is more similar to the time course of the rate modulation.

4. Discussion 4.1. Methodological aspects This is the third of a series of theoretical articles devoted to the ¨ ¨ VEMP. The first article (Lutkenh oner et al., 2010) was inspired by the pioneering work of Wit and Kingma (2006) and provided the basis for the convolution model represented by Eqs. (2) and (3). The ¨ ¨ second article (Lutkenh oner and Basel, 2011) was devoted to a specific realization of this model—simple enough to be evaluated analytically in the time and in the frequency domain and yet realistic enough to capture the basic aspects of the VEMP phenomenon. The study corroborated a basic conjecture of the first article,

1

¨ ¨ (cf. Lutkenh oner and Basel, 2011), which contains only one unknown function: the MUAP h(t). The problem would then be to determine h(t) (and the constant of proportionality, Z) such that separate evaluations of the left and the right sides of the equation yield consistent results. The function r(t), which was eliminated from the problem, could finally be calculated by a deconvolution based on Eq. (2) or (3). Under ideal circumstances, it would not matter which of the two equations is used, but differences are to be expected if there are modeling errors and noise. Moreover, differing assumptions as to the indispensable regularization would also lead to deviating deconvolution results. Matters are complicated further by the fact that a deconvolution result that is optimal in a mathematical sense does not necessarily have a reasonable physiological interpretation. This brief consideration indicates that optimizing h(t) on the basis of Eq. (20) without simultaneously keeping an eye on r(t) is problematic. Last but not least, finding a suitable function that can be characterized by just a few parameters appears to be more difficult for h(t) than for r(t). For these reasons, we refrained from investigating the outlined alternative in more detail. 4.2. Interpretation of data The deconvolution algorithm worked well not only on the simulated data, but also on the real data. This success was not assured from the outset, because our convolution model is simplistic in that it builds on the concept of a mean MUAP rather than allowing each motor unit to have its individual MUAP. The good agreement between data and model prediction suggests that the simplicity of the model does not seriously hamper its ability to quantitatively explain measured data. The analysis presented in Fig. 5 suggests that the assumption of a brief inhibition followed by a somewhat longer excitation is completely sufficient to explain the data: the inhibitory component is essentially responsible for the p13–n23 complex of the VEMP, the excitatory component for the n34–p44 complex. This confirms the qualitative interpretation given in our recent experi¨ ¨ mental paper (Lutkenh oner et al., 2011). The analysis also confirms the presumption that, at the time of the n23, contributions of an inhibitory and an excitatory component overlap. As to the VEMP, the two contributions are of the same sign and so there is an additive effect; as to the variance modulation, the contributions are of opposite sign so that partial cancellation is observed (cf. Fig. 6). But contrary to our former view, these effects are clearly not the main reason for the huge discrepancy of the p13

B. L¨ utkenh¨ oner, T. Basel / Journal of Theoretical Biology 294 (2012) 87–97

and n23 counterparts in the variance modulation. More important seems to be the fact that the two dominant MUAP peaks considerably differ in amplitude (cf. Fig. 5d). This disparity affects, in particular, the variance modulation, which depends on the MUAP squared (amplitude differences are more pronounced after squaring). Although the data could be explained quite well, it cannot be excluded that the assumed two-component rate modulation provides only a crude approximation of the real situation. According to Eqs. (2) and (3), the VEMP and the associated variance modulation may be interpreted as filtered versions of the rate modulation, with h(t) and h2(t) being the impulse responses of the filters. As these filters have a low-pass characteristic, some high-frequency features of the rate modulation might have basically no impact on the data. From a theoretical point of view, this possibility cannot be denied. Yet, the concern is probably unfounded: peristimulus time histograms from single motor units (Colebatch and Rothwell, 2004) do not suggest a much higher complexity than was assumed in our study.

5. Conclusions Our convolution model is capable of explaining VEMP and associated variance modulation quantitatively. Thus, the simplifying assumption underlying this model – namely that the time course of the MUAP is the same for all motor units – does not seriously hamper a model-based interpretation of experimental data. A deconvolution sheds new light on the VEMP phenomenon in that it provides an undistorted view of the phenomenon of

97

interest: the temporal modulation of the MUAP rate. That way, it is possible to identify and characterize individual components of the VEMP and to distinguish between inhibition and excitation.

References 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. Colebatch, J.G., Rothwell, J.C., 2004. Motor unit excitability changes mediating vestibulocollic reflexes in the sternocleidomastoid muscle. Clin. Neurophysiol. 115, 2567–2573. Halmagyi, G.M., Curthoys, I.S., Colebatch, J.G., Aw, S.T., 2005. Vestibular responses to sound. Ann. N. Y. Acad. Sci. 1039, 54–67. Lagarias, J.C., Reeds, J.A., Wright, M.H., Wright, P.E., 1998. Convergence properties of the Nelder–Mead simplex method in low dimensions. SIAM J. Optim. 9, 112–147. ¨ ¨ Lutkenh oner, B., Stoll, W., Basel, T., 2010. Modeling the vestibular evoked myogenic potential. J. Theor. Biol. 263, 70–78. ¨ ¨ Lutkenh oner, B., Basel, T., 2011. An analytical model of the vestibular evoked myogenic potential. J. Theor. Biol. 286, 41–49. ¨ ¨ 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. 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.