Alexei B Ryabov
email: a.ryabov@icbm.de

Lars Rudolf
email: rudolf@mpipks-dresden.mpg.de

Bernd Blasius
email: blasius@icbm.de

Vertical distribution and composition of phytoplankton under the influence of an upper mixed layer

Keywords: competition, coexistence, deep chlorophyll maximum, DCM PACS: 87.23.Cc, 92.20.Jt 2000 MSC: 35K57, 92D25

The vertical distribution of phytoplankton is of fundamental importance for the dynamics and structure of aquatic communities. Here, using an advection-reaction-diffusion model, we investigate the distribution and competition of phytoplankton species in a water column, in which inverse resource gradients of light and a nutrient can limit growth of the biomass. This problem poses a challenge for ecologists, as the location of a production layer is not fixed, but rather depends on many internal parameters and environmental factors. In particular, we study the influence of an upper mixed layer (UML) in this system and show that it leads to a variety of dynamic effects: (i) Our model predicts alternative density profiles with a maximum of biomass either within or below the UML, thereby the system may be bistable or the relaxation from an unstable state may require a long-lasting transition.

(ii) Reduced mixing in the deep layer can induce oscillations of the biomass; we show that a UML can sustain these oscillations even if the diffusivity is

Introduction

The survival and competition of species in an heterogeneous environment has fascinated ecologists for long times (see e.g. [START_REF] Holmes | Partial Differential Equations in Ecology: Spatial Interactions and Population Dynamics[END_REF]; Tilman and Kareiva 1997; Huisman et al. 1999;[START_REF] Neuhauser | Mathematical challenges in spatial ecology[END_REF]. In many systems the spatial diversity of natural populations originates mainly from some underlying abiotic heterogeneity of the environment. If growth conditions vary between different locations this spatial variation should be reflected in the density distribution of natural populations. After the seminal papers by [START_REF] Skellam | Random dispersal in theoretical populations[END_REF] and by [START_REF] Kierstead | The size of water masses containing plankton bloom[END_REF], the dynamics of such systems have often been analyzed in terms of favorable and unfavorable patches (see e.g. [START_REF] Okubo | Diffusion and Ecological Problems[END_REF][START_REF] Cantrell | Spatial heterogeneity and critical patch size: area effects via diffusion in closed environments[END_REF][START_REF] Birch | Bounding biomass in the Fisher equation[END_REF]). This approach assumes that space can be divided into regions

A c c e p t e d m a n u s c r i p t

of positive and negative net growth, between which the organisms are transported by diffusion and advection. These suggestions, although realistic in many situations, do not hold in some resource-consumer systems, in which the size, the form and even the location of the species' favorable patch may vary, reflecting external and internal perturbations.

An important example is provided by the dynamics of phytoplankton in an incompletely mixed water column (see e. ). The production layer in such systems is highly variable and dynamic. First it is species specific: less nutrient limited species will have a positive production rate close to the surface, whereas smaller light requirements will result in subsurface productivity (see e.g. [START_REF] Venrick | Phytoplankton seasonatity in the central North Pacific: The endless summer reconsidered[END_REF]). Second, as each species shades light and consumes nutrients it acts as an ecosystem-engineer, modifying its abiotic environment [START_REF] Jones | Organisms as ecosystem engineers[END_REF]). As a consequence, the position of the favorable layer for a phytoplankton species depends on the current abundance and distribution of its own biomass, as well as that of all other species.

A c c e p t e d m a n u s c r i p t

Theoretical models have been a useful approach to describe nutrient limited phytoplankton growth in constant and seasonally driven environments [START_REF] Huppert | A model of phytoplankton blooms[END_REF](Huppert et al. , 2004)). The dynamics, competition and structuring within a vertical water column have been investigated in a series of modeling investigations (Radach and However, in a few recent investigations either surface or deep maxima were observed under almost the same conditions [START_REF] Venrick | Phytoplankton seasonatity in the central North Pacific: The endless summer reconsidered[END_REF]; Holm-Hansen and Hewes 2004). These observations suggest the existence of bistability in the spatial phytoplankton configuration, however they cannot easily be reproduced in current models and may indicate that some important mechanisms, which are contributing to the spatial organization of aquatic communities, are still not well understood.

One major ingredient determining whether a species can establish close to the surface is the presence of an upper mixed layer (UML). UMLs commonly occur in oceans and lakes due to mechanical perturbation of the surface waters (e.g., by wind, waves, storms) and are characterized by strong turbulent mixing up to a depth of about 30 m to 200 m or more [START_REF] Deuser | Variability of hydrography and particle flux: Transient and long-term relationships[END_REF][START_REF] Venrick | Phytoplankton seasonatity in the central North Pacific: The endless summer reconsidered[END_REF][START_REF] Law | A mesoscale phytoplankton bloom in the polar Southern Ocean stimulated by iron fertilization[END_REF]. To our knowledge the first study to demonstrate

A c c e p t e d m a n u s c r i p t

the strong influence of a UML on the spatial configuration of phytoplankton was performed by Yoshiyama andNakajima (2002, 2006), who considered a single species model where the water column was divided into an infinitely mixed UML and poorly mixed lower layers, with a small diffusivity across the separating boundary layer (thermocline). The model outcome were vertical phytoplankton patterns with a sharp edge at the thermocline, with the possibility of bistability in the phytoplankton distribution, characterized by phase transitions between biomass maxima in the deep and upper layers. However, these articles did not analyze such important questions as the influence of finite mixing in the upper layer, the competition of species occupying different layers, or the dynamics of the system when deep phytoplankton maxima are not stationary, and thus many important questions about the role of the UML for structuring phytoplankton configurations still remain open.

In this study, we consider a mathematical model for the growth of a nutrient and light limited phytoplankton community in a vertical incompletely mixed water column. A UML is introduced by smoothly varying the strength of diffusivity with vertical depth, so that discontinuities are avoided. We find that a UML ameliorates the growth conditions close to the surface and can have drastic effects because organisms close to the surface occlude light and prevent growth in all deeper layers. Thus the presence or absence of a UML turns out to be a major factor controlling the vertical distribution and competition outcome in the whole water column.

In particular, the presence of a UML yields the following key model out- As a result, in the two species model, we observe bistability both in the spatial distribution and in the species composition. Fourth, we show that a light limited species can obtain a competitive advantage if the diffusivity in the deep layers is reduced below a critical value. In this case we identify a dynamic competitive exclusion effect, where the species composition is determined by the dynamic state such that oscillatory states in the deep layers are displaced by steady solutions in the UML. Finally, we present a graphical approach for analysing the mechanism of resource competition in a spatially extended system.

Model

Our model describes the dynamics of n phytoplankton species, competing for resources in a vertical water column of depth Z B (in the following we consider only the cases n = 1 and n = 2). We study a bottom-up control of phytoplankton, which means that we did not consider the influence of zooplankton and higher trophic levels. Let P i (z, t) denote the density of species i at time t and depth z. Assume that there are two limiting factors:

∂P i ∂t = growth -loss -sinking + mixing (1) = μ i P i -mP i -v ∂P i ∂z + ∂ ∂z D(z) ∂P i ∂z ∂N ∂t = -uptake + recycling + mixing (2) = -α n i=1 μ i P i + εαm n i=1 P i + ∂ ∂z D(z) ∂N ∂z
where μ i (N, I) describes the local growth rate of species i, m is the mortality, v is the phytoplankton sinking velocity, D(z) is the depth dependent turbulent diffusivity, α is the nutrient content of a phytoplankton cell and ε is the phytoplankton recycling coefficient.

Equations ( 1) and ( 2) are coupled by means of the growth rate μ i (N, I) which depends on the local resource availability and also controls the nutrient uptake. Assuming that the limitation of growth follows the Monod kinetics (e.g. [START_REF] Turpin | Physiological mechanisms in phytoplankton resource competition[END_REF]) and both resources are essential (von Liebig's law of minimum), we obtain

μ i (N, I) = μ max min N H (i) N + N , I H (i) I + I , (3)
where μ max is the maximum growth rate, and

H (i)
N and H

I(z) = I in exp -K bg z -k z 0 n i=1 P i (ξ, t)dξ , (4)
where I in is the incident light intensity, K bg is the turbidity of the water without biomass, and k is the phytoplankton light-absorption coefficient.

To describe the water column stratification, we assume that the diffusivity D(z) takes a large value D U in the upper layer and a much smaller value D D in the deep layers. A gradual transition from one area to another can be written in terms of a generalized Fermi function

D(z) = D D + D U -D D 1 + e (z-Z mix )/w , (5)
where Z mix describes the depth of the UML and the parameter w characterizes the width of the transient layer. In all numerical experiments we chose D U = 50 cm 2 /s, modeling well-mixed waters in the UML, and 1.

D D = 0.1...
We simulate the model up to a depth of z = Z B and assume impenetrable

A c c e p t e d m a n u s c r i p t

boundaries at the surface and at the bottom for the phytoplankton

v i P i -D z ∂P i ∂z z=0,Z B = 0 . (6)
Since at the bottom of a sufficiently deep water column (we assume Z B = 300 m) the phytoplankton biomass vanishes together with flux, one can also use the Dirichlet boundary condition at the bottom P i (Z B ) = 0 and this will not influence our results. For the nutrient distribution we assume an impenetrable surface and a constant concentration at the bottom

∂N ∂z z=0 = 0, N(Z B ) = N B . (7)
As initial condition we assumed a small uniformly distributed concentration of phytoplankton (P (z, 0) < 1 cell m -3 ), whereas for the nutrient we used two different initial profiles, describing a nutrient saturated water column (N(z, 0) = N B ) and a nutrient depleted upper layer (N(z, 0) = 0 if z ≤ Z mix ). We have also explored the influence of other initial conditions, however we did not find any new solutions beside the ones described. Checking the stability of solutions, in each case we simulated dynamics for 50,000 system days (approximately 130 years).

The model was integrated using a backward difference method, based on the finite volume scheme (Pham Thi et al. 2005). For the numerical solution we have discretized all variables on a grid which consisted of 600 points.

Diffusion terms were approximated by a second order central discretization scheme, the advection term was represented by a third-order upwind biased formula, integration was made via the trapezoidal rule. The resulting system of ordinary differential equations was solved by the CVODE package (http://www.netlib.org/ode). For model validation we have compared our

A c c e p t e d m a n u s c r i p t

simulation results with already published results [START_REF] Huisman | Reduced mixing generates oscillations and chaos in the oceanic deep chlorophyll maximum[END_REF]) and further verified that the results remain unchanged if we double the number of points in the grid. Furthermore, in some limiting cases it was possible to compare our simulation results with analytic solutions.

The single species model

We first concentrate on the dynamics of a single species population and describe the formation of a DCM in a water column without a UML. Suppose that we start with an initially nutrient rich system (N(z, t) = N B ). Thereby the nutrient limitation is negligible and we observe a rapid formation of a transient phytoplankton maximum close to the surface (Fig. 1A). This phytoplankton profile P (z, t) is, however, not stable. With the depletion of the nutrient in the surface layer the production layer, i.e. the layer where μ(N, I) ≥ m, shifts downwards (see also Fig. 2A), until the system reaches a stable DCM configuration (Fig. 1B). In this equilibrium the upward flux of nutrient compensates the nutrient consumption and any further sinking of the production layer is balanced by light limitation [START_REF] Klausmeier | Algal games: the vertical distribution of phytoplankton in poorly mixed water columns[END_REF].

Note that the spatio-temporal evolution of the concentration profile P (z, t) depends on the growth conditions at all vertical positions. Due to the water turbidity and phytoplankton self-shading the light intensity I(z, t) is reduced in deeper layers, thereby increasing the limitation by light (dashed lines in Fig. 1). In contrast, the nutrient concentration N(z, t) is close to zero within the bulk of phytoplankton biomass and above it, and increases almost linearly with depth below the phytoplankton peak (see dot-dashed lines in Fig. 1 for

A c c e p t e d m a n u s c r i p t

the nutrient limitation of growth). This shaping of the spatial dependence of the growth limiting factors self-consistently depends on the full phytoplankton density profile P (z, t), a fact which makes the problem very hard to understand without mathematical simulation.

A rough insight into the time evolution of P (z, t) can be gained by considering the centers of biomass Z m = zP dz/ P dz and of phytoplankton net production Z g = zgP dz/ gP dz (black and gray arrows in Fig. 1). Here g(z) describes the net phytoplankton production which includes growth, loss, sinking and mixing, i.e., the product gP equals the right-hand-side of Eq.

(1). In an incompletely mixed water column without a UML the position of the center of mass, Z m , follows that of Z g . Thereby the phytoplankton production around Z g shifts the mass center Z m , changing the local nutrient consumption and the light absorption, which in turn has an influence again on the position of the growth center Z g . In this feedback-loop, the system eventually reaches a stable equilibrium configuration where both centers coincide, Z g = Z m , giving rise to a DCM (Fig. 1B).

Now suppose that there is strong mixing in the upper layer. If the bulk biomass is located sufficiently deep, then the mixing in the upper layer has practically no effect and an identical DCM can persist (Fig. 1C) independent of whether or not a UML is present. Note that we always assume that the depth of the UML is smaller than the compensation depth (i.e., a depth at which μ(I) = m in the absence of biomass, see e.g. [START_REF] Sverdrup | On conditions for the vernal blooming of phytoplankton[END_REF], Huisman et al. 1999b).

By contrast, if the phytoplankton biomass is initially located close to the surface, it will be almost uniformly distributed within the UML. The position

A c c e p t e d m a n u s c r i p t

of Z m is then fixed approximately in the middle of the UML and is almost uncoupled from Z g (Fig. 1D). Therefore a gradual shift of the center of mass into the deep layers is no longer possible and the transition to a DCM can only take place if the light intensity below the UML is sufficiently large to provide positive net growth in the deep layers -otherwise the phytoplankton remains trapped in the UML. We denote this stable configuration of a nearly uniform phytoplankton profile in the UML as an upper chlorophyll maximum (UCM, note that all acronyms are listed in Table 2). Figs. 1C and1D show that in the presence of a UML, depending on the initial conditions, the system can undergo two very different spatial configurations of either a deep or an upper phytoplankton maximum [START_REF] Yoshiyama | Catastrophic transition in vertical distributions of phytoplankton: alternative equilibria in a water column[END_REF]. Moreover, if the diffusivity is lower than the minimal persistence threshold

D min = v 2 4(μ(N B , I in ) -m) (8)
the sinking phytoplankton population cannot survive and goes extinct [START_REF] Riley | Quantitative ecology of the plankton of the Western North Atlantic[END_REF][START_REF] Shigesada | Effects of competition and shading in planktonic communities[END_REF]. Our numerical simulations show that in a system without a UML these three model outcomes (i.e., stationary DCM, oscillating DCM or extinction) are the only possible long-term solutions.

Further, these solutions are globally attracting, which means that they are

A c c e p t e d m a n u s c r i p t

reached from any initial condition.

In the presence of a UML these solutions can be found as well, however the dynamics may be more complicated. Most notably, as already mentioned, the system may be bistable: under the same parameters as in Fig. 2A, an initially nutrient saturated water column gives rise to a UCM (Fig. 2C), whereas an initially nutrient depleted water column leads to a DCM (Fig. 2D). Note that we observe bistability only in a certain parameter region in which a DCM is not affected by the upper layer and a UCM contains enough biomass to limit 3A). Note that the lower border of the bistability range is quite well described by the analytical criterion Eq. ( 16), which is derived from the condition that the phytoplankton net production rate below the mixed layer is not positive (dash-dotted lines in Fig. 3). Interestingly, outside the region of bistability close to the critical point smooth transitions from deep to surface maxima are possible. In this case intermediate density profiles with no clear separation between DCM and UCM appear and the phytoplankton biomass is located in both parts of the water column. To visualize the distribution of biomass in this regime, parameter combinations, for which the median of the biomass distribution is located at Z mix (i.e., half of the biomass is distributed within and half below the UML), are indicated as thick solid lines in Fig. 3. In the (N B , v) parameter plane the bifurcation lines have almost vertical structure (Fig. 3C). Most notably, the minimal value of N B which can lead to upper biomass maxima practically is independent from v. Thus sinking has no strong influence on the stationary solutions in a well mixed layer. By contrast, oscillatory solutions arise only for a large enough sinking velocity.

Again, the population can persist beyond the persistence threshold v > v max , where v max is the maximal sinking velocity which allows the survival of a population in a system without a UML. v max can be obtained from Eq. ( 8)

for a given diffusivity.

The two species model

We now turn to the influence of a UML on the competition of two species which are differently adapted to the conditions at the surface or in the deep water. For this, we extend the model to contain two phytoplankton species which differ in their respective half saturation constants H I and H N in Eq. ( 3). More specifically, we consider the competition between a most light-limited species (I-species) characterized by a low H N and a large H I value and so is forced to live close to the light source, and a most nutrientlimited species (N-species) with low H I and large H N that will usually do better in deep water. Note that the results of the previous section hold for

A c c e p t e d m a n u s c r i p t

the N-species.

To avoid the influence of different growth, mortality and consumption rates, we study the simplest possible case and keep the other species' parameters identical in both species (see Table 1). Thereby, in a well mixed uniform environment (e.g., a chemostat) these species would have parallel consumption vectors, so that the success of one species over the other depends only on the resource concentrations in the absence of biomass; with the consequence that in a chemostat such species cannot coexist and the outcome of their competition cannot be bistable [START_REF] Tilman | Resources: a graphical-mechanistic approach to competition and predation[END_REF][START_REF] Tilman | Resource competition and community structure[END_REF]. In a well mixed water column, where the light intensity reduces with depth, the competition of these species is more complicated, however the same results still hold (Huisman and Weissing 1995); whereas in an incompletely mixed water column, assuming only the light gradient, both species can coexist in a narrow parameter region (Huisman et al. 1999). In the following we show that the outcome of competition is drastically altered in the presence of two opposing resource gradients combined with a space dependent diffusivity, resulting in new regions of competitive exclusion, bistability and coexistence. Finally reducing D D to very small values leads to a surprising result: the I-species always wins the competition and the steady UCM formed by the I-species becomes the single possible attractor in the system, independent of the initial conditions (region I, Fig. 4J). Thus, a strongly light limited species is able to establish a steady population in the UML, outcompeting a less light limited species which, in the absence of the former, would form an oscillatory maximum in the deep weakly mixed layers. Note that the I-species survives even if D D < D min , provided that the nutrient concentration in the upper

A c c e p t e d m a n u s c r i p t

layer is sufficiently high (Fig. 4K).

Fig. 6 shows the outcome of competition between two species in the (N B , D D ) coordinate plane. In the system without a UML (Fig. 6A) the N-species either wins or dominates (more than 75% of the biomass). 6B), when the N-species solely would still occupy the UML, but becomes a weaker competitor under stronger nutrient limitation. We also note that change of the vertical profiles causes non-monotonous shifts in the species composition.

Graphical approach to analyze the competition outcome

In a spatial system the outcome of competition is hard to deduce from first principles. First, the production layers of differently adapted species typ-

A c c e p t e d m a n u s c r i p t

ically will be located at different depths, which tends to reduce the strength of interspecific competition. Second, the competition becomes indirect, since the occurrence of phytoplankton at a certain depth can completely reshape the resource distributions in other parts of the water column. To obtain an insight into the mechanisms of competition in such an environment, we suggest a graphical approach, which is based on a well known method developed by Tilman for mixed systems (1980,1982).

If all resources and the biomass are uniformly distributed, the state of the system may be represented as a point in a multidimensional resource space (Fig. 7A). As the biomass grows and takes-up resources this system state point moves in the resource plane (along the consumption vector) until it hits the zero net growth isocline, characterized by a balance between growth and loss processes. From the condition m = μ i (N, I) in equation ( 3) we can find such values N and I so that the specific growth rate equals the mortality

N * = m μ max -m H N , I * = m μ max -m H I . (9)
These values determine the location of the zero net growth isoclines (Fig. 7).

The isoclines divide the resource plane into areas of positive and negative local population net growth and mark possible resource combinations in equilibrium.

In a spatially-extended system, where the resource values change with depth z, the state of the system can be represented as a curve in the resource space, which is parametrically determined by the resource values (N (z), I(z)). In this sense, the notion of a system state point in a homogeneous system naturally extends to a system state curve (SSC) in a spa-

A c c e p t e d m a n u s c r i p t

tial explicit system. A special case is given by a well mixed water column (Fig. 7A): while the nutrients are uniformly distributed, the light intensity decreases exponentially with depth and the SSC reduces to a line segment {N = const, I out ≤ I ≤ I in }, where I out is the light intensity at the bottom of the water column (Huisman and Weissing 1995).

To simplify the discussion, in the following we always focus on the system state curve (SSC) in equilibrium. Fig. 7A shows a typical simulation outcome for the SSC of a single-species population in an incompletely mixed water column. Note that we use logarithmic scaling to magnify the location of the SSC close to the zero net growth isoclines. Due to the diffusive mixing the SSC does not settle at the zero net growth isoclines, as would be the case for a uniform distribution of the resources. Instead the SSC extends into the area of positive ("favorable range"), as well as into the area of negative population net growth ("unfavorable range"). To give a crude insight into the density variation along a SSC we indicate the central range of the curve which contains 90% of the biomass as a thick line.

For the analysis of a multi-species system, we find it convenient to calculate the SSC independently for each species in the absence of all other species. So each species attains its own SSC which would indicate its singlespecies resource configuration in equilibrium. Fig. 7B shows an apparent example, illustrating the case of competitive exclusion. In this example the SSC of the N-species lies below the null growth isocline of I-species. Thus, the N-species reduces the resource levels at all vertical positions to values which are smaller than the resource requirements of the I-species and so do not permit positive net-growth of the invading I-species. By contrast, the I-

A c c e p t e d m a n u s c r i p t

species' SSC is located mostly above the null growth isocline of the N-species and so allows a positive net-growth of the N-species. As a consequence, in this example the I-species will always be excluded as it has larger resource requirements.

Fig. 7C demonstrates an example of the coexistence of two species. In this figure both curves intersect in such a way, that the SSC of one species is below that of its competitor in an essential part of its favorable range, and vice versa. In these conditions both species can coexist because both are superior competitors at different depths. In general, however, if the SSC of one species lies above the zero null growth isocline of its competitor, there is no unambiguous answer to the question "whether or not the latter species can invade?" The possibility of invasion will depend on the principal a UML, now the SSC of the I-species lies below the main part of the Nspecies' SSC. In this example, the N-species can exist in the presence of the I-species, however its fraction is small so that its contribution to the total resource distribution is negligible. As a result, in the presence of a UML the outcome of competition is bistable and depends on initial conditions.

Our investigations show that the analysis of a species' SSC can further our understanding about resource competition in heterogeneous multi-species systems. What is more, one can obtain rigorous results for biologically similar species. A more detailed analysis will be presented in a forthcoming article (Ryabov and Blasius, in preparation).

Discussion

In this article we investigated the influence of an upper mixed layer on the distribution and competition of phytoplankton species in a water column, in which inverse resource gradients (of light and a nutrient) can limit growth of the biomass. In this system the location of a production layer is not fixed, rather it depends on the initial and boundary conditions, on the species requirements, and on the stage of the process. Together with the presence of differently mixed areas (i.e., the UML and lower layers), this leads to a plethora of phenomena, including bistability of phytoplankton profiles, changes in the competition outcome, and new critical conditions for the survival of a phytoplankton population. Thus our study not only proves the UML to be an important factor with the potential to shape the spatial distribution and species composition of phytoplankton, but also reveals insights of general ecological importance.

A c c e p t e d m a n u s c r i p t

While previous theoretical investigations have usually focused on either a fully mixed or an incompletely mixed system, the presence of a UML requires a combination of these approaches. Including both factors, by dividing the water column into two separate compartments, Yoshiyama andNakajima (2002, 2006) showed the existence of bistability in the spatial distribution of a phytoplankton monoculture. Our work confirms and extends these findings on the following key points. First, assuming a gradual change of diffusivity with depth, our modeling approach integrates the whole water column into a single framework. This allows us, for example, to investigate the influence of mixing in both the upper and the lower layers. Secondly, we analyze the competition of species which are differently adapted to the availability of nutrients and light. As we show, this has drastic effects because the species composition strongly correlates with the spatial patterning. And finally, our analysis includes the case of a stratified lower layer, when a stable DCM cannot persist and oscillatory or chaotic solutions appear, which again have a strong influence on the species competition in the system.

In comparison with lower layers three factors make a UML more favorable for phytoplankton species. Firstly, a UML has only one border with an unfavorable environment below the euphotic zone. In contrast, a deep production layer has two such boundaries and diffusion of cells upward and downward from it leads to additional losses due to either light or nutrient limitations. Secondly, the strong mixing in the upper layer allows for the survival of a sinking phytoplankton population. This additionally supports a population of phytoplankton even if the diffusivity in the lower layers becomes very small. Higher diffusivity also reduces the effective mortality rate

A c c e p t e d m a n u s c r i p t

of sinking phytoplankton species (see Appendix A.3). Thirdly, a UML promotes a nearly uniform distribution of nutrients, which makes the nutrient consumption more efficient and gives an additional competitive advantage for a species inhabiting a UML. All these factors decrease the total loss rate and, hence, reduce the resource requirements. However, we note that a deep and turbid UML can also play a negative role, as it can lead to the extinction of species due to light limitation (Huisman et al. 1999b).

Inclusion of a UML in the single species model allows for the existence of two alternative density configurations with either a deep or an upper phytoplankton maximum. Note that in a small parameter range we observe a third intermediate profile with an essential part of the biomass above and below the thermocline. Moreover, there is a range of parameters for which the system is bistable and the appearence of either deep or surface biomass maximum is determined by the initial distribution of nutrients. Both the consumption of nutrient and the self-shading of light are necessary conditions for this behavior. Since the biomass obstructs the upward nutrient flow, it makes the upper layer unfavorable, provided that a deep maximum of phytoplankton has established. The shading of light is an opposite mechanism, which prevents population growth in the lower layers. The third requirement for bistability is the presence of strong mixing in the upper layer. This decouples the locations of the production layer from that of the bulk biomass and prevents a drift of the population toward a DCM, as it would occur in a system without a UML. Similar bistable behaviour has also been suggested to occur in field data. For example, in a recent study of Antarctic waters (Holm-Hansen and Hewes 2004) it was found that adjacent spots could ran-

A c c e p t e d m a n u s c r i p t

domly exhibit surface maxima or DCM, even though temperature, salinity and water density profiles of the water columns were practically identical.

The occurrence of bistability may have important ecological consequences and usually goes together with catastrophic shifts. Imagine that we slowly change a parameter, i.e., the incident light intensity, so that its value passes through the bistability range. Then there will be a point, where the current state suddenly loses its stability, thereby inducing a rapid transition from one phytoplankton profile to another. This is also important for a border between oligotrophic and eutrophic waters, where this transition might be induced by seasonal changes and leads to a shift of the interface between a deep and a surface biomass. maximum. In this scenario bistability will lead to a lag in the transition from one state to another and back.

Outside the bistability range the system possesses only one attractor.

However, the transition from an unstable to a stable solution may take a long time (especially close to the bistability range) and includes two stages.

During the first stage, the nutrient concentration slowly evolves to a level which allows for the formation of a stable biomass profile. This stage occurs on the slow diffusive time scale and may typically last 5-10 years in deep waters. Once the nutrient values have crossed a critical value, a rapid transition to the stable profile is triggered. This second stage develops on the biological time scale and lasts approximately 10-50 days.

As a consequence, in this, and in the bistable, regime the system is very sensitive to disturbance. If, due to some perturbation, a large quantity of nutrients is transported into the UML, the system can rapidly switch from a DCM to a UCM and remain in this state for a long time. Similar, a The mechanisms ensuring persistence ultimately are related to the fact that locally, within the well-mixed upper layer, diffusivity is much larger than the persistence threshold. Thus obviously, a stationary distribution of a strongly light limited species within the UML will not be affected by low values of diffusivity in deeper layers. However more astonishingly, that the persistence threshold is absent for the oscillatory solutions for which the bulk of the biomass remains located in the weakly mixed lower layers for most of the time in the cycle. Again, this effect is related to the mobility of the production layer, which for an oscillatory DCM undergoes cycles in depth. These vertical cycles of the production layer can act as an "oscillatory pump", so that during some, possibly very small, phase of the cycle the production layer is located within the UML. These short cyclic visits of Possible interesting further extensions of our model would be the inclusion of zooplankton or higher trophic levels, the inclusion of an Allee effect or the extension of our study into a three-dimensional system.

A c c e p t e d m a n u s c r i p t

Acknowledgments

We thank Thilo Gross and Wolfgang Ebenhöh for advise and useful discussions. This study was supported by German DFG (SFB 555) and German VW-Stiftung.

A. Analytical derivations

A.1. Total equilibrium phytoplankton biomass in the single-species model At equilibrium the left hand sides of Eqs. ( 1) and ( 2) equal zero

μ(N, I)P -mP -v ∂P ∂z + ∂ ∂z D(z) ∂P ∂z = 0 (10)
-αμ(N,

I)P + εαmP + ∂ ∂z D(z) ∂N ∂z = 0 (11)
After integration of these equations with respect to the boundary conditions, Eqs. ( 6) and [START_REF] Cantrell | Spatial heterogeneity and critical patch size: area effects via diffusion in closed environments[END_REF],

Z B 0 μ(N, I)P (z)dz -m Z B 0 P (z)dz = 0 -α Z B 0 μ(N, I)P (z)dz + εαm Z B 0 P (z)dz = -D(Z B ) ∂N(z) ∂z z=Z B ,
we obtain the total biomass

W = Z B 0 P (z)dz = D D αm(1 -ε) ∂N(z) ∂z z=Z B . (12)
This formula relates the phytoplankton biomass and the nutrient flux and can be interpreted as a conservation law in the system.

A c c e p t e d m a n u s c r i p t

To describe the nutrient flow we first consider the region, which is free of the phytoplankton biomass. Using Eq. ( 11) we obtain

D(z) ∂N(z) ∂z = const . (13)
This equation corresponds to the stationary diffusion flow and goes together with a linear gradient of nutrients, below the phytoplankton maximum. If the total phytoplankton biomass is located within the UML, the nutrient flow can be estimated as

D(z) ∂N(z) ∂z = D D N B -N * Z B -Z mix , (14)
where N * from Eq. ( 9) approximates the concentration of the nutrient within the UML. Substituting this into Eq. ( 12), we obtain the total phytoplankton biomass in the UML (see Fig. 8).

W = D D αm(1 -ε) N B -N * Z B -Z mix . (15)

A.2. Border of the stability of a UCM

A sufficient condition for stability of a UCM is the light limitation of growth below the UML, i.e., the light intensity below the UML should be smaller than the critical light intensity I * , see [START_REF] Deuser | Variability of hydrography and particle flux: Transient and long-term relationships[END_REF]. Using Eq. ( 4), we obtain

I in exp -K bg Z mix -k Z mix 0 P (z)dz < I * . (16)
If the total phytoplankton biomass is located in the UML, equations [START_REF] Holmes | Partial Differential Equations in Ecology: Spatial Interactions and Population Dynamics[END_REF] and ( 16) give the following criterion for stability of the upper maximum

ln (I in /I * ) -K bg Z mix < D D (N B -N * )k αm(1 -ε)(Z B -Z mix ) . ( 17
)
This line is shown in Figs. 3A -3C as the lower boundary of the bistability range and is in a good agreement with the numerical simulations.

A c c e p t e d m a n u s c r i p t A.3. Losses in the UML

The upper mixed layer is more favorable for sinking phytoplankton species.

To show this, consider a water column where the diffusivity does not depend on depth. Its is easy to show that washout from the water column can be interpreted as an additional mortality term. Substituting P = P exp (vz/2D) into Eq. ( 1), we obtain

∂ P ∂t = μ P -m + v 2 4D P + D ∂ 2 P ∂z 2 .
Note that ∂ t P has the same sign as ∂ t P , thereby both functions grow and decline simultaneously. Introduce the new mortality

m = m + v 2 4D , (18)
and substitute it in the expression for the limiting resource values [START_REF] Deuser | Variability of hydrography and particle flux: Transient and long-term relationships[END_REF]. Since m > m, the new limiting values I * and N * in the presence of sedimentation should be larger, see [START_REF] Deuser | Variability of hydrography and particle flux: Transient and long-term relationships[END_REF]. This results in a shift of the zero net growth isoclines towards higher values of resources. However, in a well mixed layer the term v 2 /4D vanishes, leading to lower resource requirements. Thus a UML creates more favorable conditions for the sinking phytoplankton biomass.

  g. Jamart et al. 1977; Klausmeier and Litchman 2001; Beckmann and Hense 2007). In many regions of the world's ocean the form of vertical phytoplankton profiles is determined mainly by two factors. First, the reduction of the light intensity with depth makes deep layers unfavorable for photosynthetic phytoplankton species. Second, an opposing gradient of nutrients can often maintain a positive net production only in deep subsurface layers. These deep states, known as deep chlorophyll maxima (DCM), constitute one of the most striking characteristics of nutrient poor waters in ocean ecosystems and freshwater lakes (Steele and Yentsch 1960; Anderson 1969; Abbott et al. 1984; Cullen 1982; Tittel et al. 2003

A c c e p t e d m a n u s c r i p t range dynamics are

  comes: First, the spatial density profile can become bistable with vertical maxima either close to the surface or in deep layers. Close to the bistability characterized by time scale separation and the relaxation from an unstable state may require a long-lasting transition. Second, in the presence of a UML the persistence threshold of sinking phytoplankton disappears. Thus, our model provides a new mechanism for the persistence of a population in a flow[START_REF] Hershey | Stable isotopes resolve the drift paradox for Baetis mayflies in an arctic river[END_REF][START_REF] Huisman | How do sinking phytoplankton species manage to persist?[END_REF][START_REF] Ryabov | Population growth and persistence in a heterogeneous environment: the role of diffusion and advection[END_REF]. Third, a UML can strongly modify the competition outcome between different phytoplankton species by providing a vertical niche for species which are better adapted to the conditions close to the surface.

A c c e p t e d m a n u s c r i p t

  the concentration of a nutrient, N(z, t), and the light intensity I(z, t), both of which are a function of the vertical position z. Coupling of the nutrient and the phytoplankton dynamics leads to the following system of reactionadvection-diffusion equations[START_REF] Radach | The vertical structure of phy[END_REF][START_REF] Jamart | A theoretical study of phytoplankton growth and nutrient distribution in the Pacific Ocean off the northwestern U.S. coast[END_REF]; Klausmeier and Litchman 2001;[START_REF] Huisman | Reduced mixing generates oscillations and chaos in the oceanic deep chlorophyll maximum[END_REF]

I

  are the corresponding half saturation constants. Varying H N and H I allows to model, for A c c e p t e d m a n u s c r i p t instance, a species which is better adapted for light (a smaller H I ) or for the nutrient (a smaller H N ). Light dissipates with depth as it is absorbed by the biomass, water, clay particles and many other absorbing substances. Assume that the light intensity decreases exponentially according to Lamber-Beer's law (see e.g. Shigesada and Okubo 1981, Kirk 1994)

Fig. 2

 2 Fig.2depicts the typical spatio-temporal evolution of the phytoplankton density. Without a UML (Fig.2A), an initial phytoplankton maximum at the surface slowly moves downward until the distribution converges to a stable DCM equilibrium. However, as is shown in Fig.2B, DCM's are not necessarily steady states. For small values of diffusivity D D the sinking of biomass may destabilize the density profile, yielding sustained regular or chaotic oscillations of biomass in the deep layer[START_REF] Huisman | Reduced mixing generates oscillations and chaos in the oceanic deep chlorophyll maximum[END_REF]).

A c c e p t e d m a n u s c r i p t 1 0

 1 Fig. 3A, large values of I in in general lead to a DCM and large N B to a UCM, while for intermediate resource levels we observe a region with bistable behavior. The bistability range is reduced for smaller values of I in or N B and disappears at a critical point (I in ≈ 350 μmol photons m -2 s -1 and N B ≈ 25 mmol/m 3 in Fig.3A). Note that the lower border of the bistability

Fig. 3B demonstrates

  Fig. 3B demonstrates the bistability range in the (N B , D D ) parameter plane, where the bifurcation lines have an almost hyperbolic form. A UCM

Figs. 3D -A c c e p t e d m a n u s c r i p t is not affectedA c c e p t e d m a n u s c r i p ta

 3D Figs. 3D -3F present the two-species stability diagrams in the different parameter planes, listing all possible model outcomes. Basically, the overall shape of the transition lines between bistable and UCM/DCM states remain identical to those of the one-species model. The most notable difference is that the bistability range (regime V) has become much wider and is extended toward smaller values of N B (compare top and bottom panels in Fig.3). By contrast, the transition line to the UCM remains largely unchanged, since it is mainly determined by the minimal depth of the deep maximum at which it

  Furthermore, we observe no phytoplankton for diffusivity below D min . Comparison of Fig. 6A with Figs. 6B and 6C show that a UML plays a crucial role improving the competitive abilities of the I-species, which occupies a shallower depth because of its higher light requirements. Shifts in the competition outcome can be caused by a change in the resource availability or by a bifurcation in the vertical biomass distribution. The thick lines in Fig. 6B and 6C indicate those transitions in the species composition which are caused by a change of vertical biomass profiles. These lines coincide with the bifurcation lines in Fig. 3E. In the presence of a UML the I-species can win under the following conditions: (i) when the diffusivity in the lower layers is reduced (D D < 0.1cm 2 /s) and the N-species alone would exhibit oscillating behavior; (ii) in the range of bistability, where either a UCM of I-species or a DCM with a large fraction of N-species appears; (iii) in the UML for sufficiently strong D D and low N B ( 15 mmol n/m 3 in Fig.

A c c e p t e d m a n u s c r i p t

  eigenvalue of a reaction-diffusion equation, characterizing the growth of the biomass[START_REF] Cantrell | Spatial heterogeneity and critical patch size: area effects via diffusion in closed environments[END_REF][START_REF] Ryabov | Population growth and persistence in a heterogeneous environment: the role of diffusion and advection[END_REF], a problem which unfortunately can only be solved in some simple cases.As is shown in Fig.7Dthe analysis of SSCs can also help to gain insight into the role of a UML. Here we use the same set of parameters as in Fig.7B, but in the presence of a UML. In this case the main part of the I-species biomass is located close to the surface, while the N-species still forms a DCM.Hence the presence of the UML leaves the SSC of the N-species unaffected while the SSC of the I-species becomes steeper and is shifted in direction of smaller resource values. This change in the form of the I-species' SSC results in new intersections of both system state curves with the result to improve the competing abilities of the I-species. While the presence of the N-species still completely prevents an invasion of the I-species, as in the case without

A c c e p t e d m a n u s c r i p t prolonged lack ofA c c e p t e d m a n u s c r i p t

  light can induce a UCM, when the reduced phytoplankton density can not prevent diffusion of nutrients to the surface. While outside of the bistability range the reverse transition in principle is possible without further external perturbation, such a flip back to the original state will usually take a long time. Thus, the effective size of the bistability regime, i.e. the parameter regime in which two dynamic states are preserved for long times and disturbance can induce long-lasting shifts in the system, is much larger than the true bistability range and comprises a large part of the model's parameter space. Since the relaxation time is long, in a seasonal environment such a system might exhibit only an unstable state (most probable an upper chlorophyll maximum).The effects of a UML are even more pronounced in a system that includes two competing species which are differently adapted to light and nutrient limitation. To analyze this situation we presented a graphical approach which helps to understand the mechanisms of competition in spatially extended systems. Moreover, we found that in the range of parameters where the two species can independently form an upper or a deep biomass maximum, the two species model demonstrates bistability both in the spatial distribution and in the competition outcome. Compared to the single-species model, the bistability range is considerably enlarged. Note that the bistability in the competition outcome is induced by the UML, whereas in a homogeneously mixed water column these species either coexist or only one species wins. Another remarkable finding is the survival of a sinking phytoplankton population even if the diffusivity in the deep layers cannot prevent population washout (Shigesada and Okubo 1981; Speirs and Gurney 2001; Straube and Pikovsky 2007). While the retention of phytoplankton in the upper layer can be provided by many factors such as the regulation of algal buoyancy or the formation of eddies, etc (Raymont, 1980), our model shows another dynamic mechanism: in the absence of phytoplankton in the deep layers, the nutrient can diffuse upward into the UML, where the population can start to grow as the nutrient concentration reaches a sufficient level. However if the light limitation is not strong enough the biomass shifts into the deeper layers where it cannot outgrow the sinking and ultimately declines being limited by light. Then the nutrient again can diffuse upwards and the cycle repeats.Thus we identify a novel mechanism of how a population can overcome the drift paradox, as the persistence threshold disappears in the presence of a UML.

A c c e p t e d m a n u s c r i p t

  the UML are sufficient to refuel biomass and to prevent extinction of the population. Thereby, however, the biomass can be temporarily reduced to very small values so that different model outcomes might be expected by including an Allee effect.Finally, we have identified an interesting competitive exclusion effect, where a species can be outcompeted in dependence of its dynamical state. This occurs for low values of deep diffusivity, when the DCM of the N-species becomes oscillatory[START_REF] Huisman | Reduced mixing generates oscillations and chaos in the oceanic deep chlorophyll maximum[END_REF]) and a further reduction of diffusivity favors to the domination of the strongly light limited I-species, which still occupies the UML. At the first glance this result is counterintuitive, because usually a decrease of nutrient transport results in an increased depth of the biomass maximum and here we suddenly obtain a UCM solution instead of a DCM. The replacement of the oscillatory solutions occurs because the I-species forms a non-oscillatory UCM and shades light, inhibiting the rapid outburst of the biomass of the N-species in the lower layers. As a consequence the (non-oscillatory) UCM configuration in the UML is able to replace the oscillatory deep maximum. As some climate models predict higher water stratification with an increase of temperature[START_REF] Bopp | Potential impact of climate change on marine export production[END_REF][START_REF] Sarmiento | Response of ocean ecosystems to climate warming[END_REF]), our findings may help to predict and understand future changes in phytoplankton patterns ongoing with global climate change.

A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t

  toplankton growth dynamics; a mathematical model. Mem. Soc. Roy. Sci. de Liege, 6e serie: 113-146

A c c e p t e d m a n u s c r i p tFigure 1 :

 1 Figure 1: Typical vertical phytoplankton profiles in the single species model for a system without a UML (top) and with a UML (bottom, depth Z mix indicated by the horizontal dashed line). Plotted are the density of phytoplankton as function of depth, P (z) (black solid line), and the growth limiting terms with respect to light, I/(I + H I ) (gray dashed line), and to the nutrient, N/(N + H N ) (gray dot-dashed line); positive growth is possiblewhere both regulating terms are larger than m/μ max (level of zero net growth indicated by the vertical dotted line). Black and gray arrows show the centers of biomass Z m and net production Z g , respectively. (A) without a UML, a non-stable phytoplankton maximum close to the surface is driven downwards (indicated by the arrows) and evolves to the stable state (B). Under the same conditions in a system with a UML, we observe two alternative stable configurations: (C) a phytoplankton profile with a maximum in the deep layers (DCM), or (D) a profile with a maximum in the upper layer (UCM).

A c c e p t e d m a n u s c r i p tFigure 2 :Figure 3 :

 23 Figure 2: Spatio-temporal evolution of typical phytoplankton profiles. Plotted is the phytoplankton density P (z, t) in color-coding (×10 7 cells m -3 ) for different of D D and initial conditions, in a system without a UML (A-B) and in the presence of a UML (C-H). Black lines track the evolution of the center of biomass, Z m , white dashed lines indicate the depth of the UML. Dynamics without a UML: (A) gradual evolution of a DCM, D D = 0.3 cm 2 /s, (B) oscillations of the biomass for small mixing, D D = 0.12 cm 2 /s. Bistability (with a UML): depending on the initial conditions either (C) a stable UCM or (D) a stable DCM is formed (value of D D as in A). Transition from an unstable state: for any initial condition (E) only a DCM is stable, D D = 0.2 cm 2 /s, or (F) only a UCM is stable, D D = 0.4 cm 2 /s, however the transient process may last a long time. Oscillations of the biomass: (G) oscillations are not affected by a UML, D D = 0.12 cm 2 /s (value of D D as in B), or (H) oscillations are induced by a UML, D D = 0.04 cm 2 /s (value of D D < D min = 0.0408 cm 2 /s).

A c c e p t e d m a n u s c r i p tFigure 5 :

 5 Figure 5: Vertical cross section for N B = 20 mmol/m 3 through the (N B , D D ) parameter plane of Figs. 3B and E, showing the center of mass, Z m , of the I-species (triangles) and the N -species (circles) as a function of D D . (A) Monoculture of the I-species, (B) monoculture of the N -species, (C) two-species system. In (B) and (C) several dynamic regimes (separated by vertical dotted lines) can be distinguished (see text). The competition outcome in (C) : I -I-species wins, II -bistability: either I-species wins or the species coexist due to oscillations, III -coexistence due to oscillations, IV -N -species wins, V -bistability, VI -coexistence. In the regimes of coexistence, the fraction of biomass is indicated by the intensity of gray color.

  Maier-Reimer 1975; Shigesada and Okubo 1981; Varela et al. 1992; Huisman and Weissing 1995; Huisman et al. 1999; Klausmeier and Litchman 2001; Diehl 2002; Hodges and Rudnick 2004; Huisman et

al. 2006; Beckmann and Hense 2007). It was shown that a given set of parameters may lead to either a surface or a deep chlorophyll maximum, whereby the location of the maximum is entirely determined by the environmental conditions. These model results are in agreement with many field studies (e.g., Aristegui et al. 2003; Matondkar et al. 2005; Weston et al. 2005).

Table 1 :

 1 Parameters values and their meaning

	Symbol Interpretation	Units	Value
	Independent variables
	t	Time	h	-
	z	Depth	m	-
	Dependent variables
	P (z, t) Population density	cells m -3
	I(z, t)	Light intensity	μmol photons m -2 s -1
	N(z, t) Nutrient concentration	mmol nutrient m -3
	Parameters
	I in	Incident light intensity	μmol photons m -2 s -1 600 (100 -600)
	K bg	Background turbidity	m -1	0.045
	k	Absorption coefficient of phytoplankton m 2 cell -1	6×10 -10
	Z B	Depth of the water column	m	300
	Z mix	Depth of the upper mixed layer	m	50
	w	Characteristic width of the thermocline m	1
	D D	Turbulent diffusivity in the deep layers cm 2 s -1	0.3 (0.04 -1)
	D U	Turbulent diffusivity in the UML	cm 2 s -1	50
	μ max	Maximum specific growth rate	h -1	0.04
	H I	Half saturation constant of light	μmol photons m -2 s -1 20; 98
		limited growth for N-and I-species
	H N	Half saturation constant of nutrient	mmol nutrient m -3	0.0425; 0.015
		limited growth for N-and I-species
	m	Specific loss rate	h -1	0.01
	α	Nutrient content of phytoplankton	mmol nutrient cell -1	1 ×10 -9
	ε	Nutrient recycling coefficient	dimensionless	0.5
	v	Sinking velocity	m h -1	0.042
	N B	Nutrient concentration at Z B	mmol nutrient m -3	5-100

Table 2 :

 2 Acronyms

	Symbol Interpretation
	DCM	Deep chlorophyll maximum
	ODCM Oscillatory or chaotic
		deep chlorophyll maximum
	UCM	Upper chlorophyll maximum
	UML	Upper mixed layer
	SSC	System state curve

A c c e p t e d m a n u s c r i p t


Bernd Lütkenhöner

Wolfgang Stoll

Türker Basel

Modeling the vestibular evoked myogenic potential

Keywords: VEMP, vestibular function testing, sacculus, sonomotor responses, electromyogram, EMG

Measuring the vestibular evoked myogenic potential (VEMP) promises to become a routine method for assessing vestibular function, although the technique is not yet standardized. To overcome the problem that the VEMP amplitude depends not only on the inhibition triggered by the acoustic stimulation of the vestibular end organs in the inner ear, but also on the tone of the muscle from which the potential is recorded, the VEMP is often normalized by dividing through a measure of the electromyogram (EMG) activity. The underlying idea is that VEMP amplitude and EMG activity are proportional. But this would imply that the muscle tone is irrelevant for a successful VEMP recording, contradicting experimental evidence. Here, an analytical model is presented that allows to resolve the contradiction. The EMG is modeled as the sum of motor unit action potentials (MUAPs). A brief inhibition can be characterized by its equivalent rectangular duration (ERD), irrespective of the actual time course of the inhibition. The VEMP resembles a polarity-inverted MUAP under such circumstances. Its amplitude is proportional to both the ERD and the MUAP rate. The EMG activity, by contrast, is proportional to the square root of the MUAP rate. Thus, the normalized VEMP still depends on the muscle tone. To avoid confounding effects of the muscle tone, the standard deviation of the EMG could be considered. But the inhibition effect on the standard deviation is small so that the measuring time would have to be much longer than usual today.

Introduction

Although several structural specializations have evolved within the vertebrate labyrinth that reduce the effect of ambient pressure changes on the balance organs, sufficiently loud sounds may nevertheless stimulate the vestibular end organs in the inner ear [START_REF] Carey | Evolutionary changes in the cochlea and labyrinth: Solving the problem of sound transmission to the balance organs of the inner ear[END_REF]. In a report about vestibular side effects of extremely loud sounds, von [START_REF] Von Békésy | Über akustische Reizung des Vestibularapparates[END_REF] speculated that the movement of the stapes might cause eddies that spread to the vestibular organs. In humans, the vestibular organ with the highest sound sensitivity seems to be the sacculus [START_REF] Townsend | The averaged inion response evoked by acoustic stimulation: its relation to the saccule[END_REF][START_REF] Todd | A saccular origin of frequency tuning in myogenic vestibular evoked potentials? Implications for human responses to loud sounds[END_REF][START_REF] Sheykholeslami | The otolithic organ as a receptor of vestibular hearing revealed by vestibular-evoked myogenic potentials in patients with inner ear anomalies[END_REF], presumably due to its close proximity to the stapes. Other vestibular end organs appear to be acoustically responsive as well, though (see e.g. Young et al., 1977;Carey et al., 2004;Xu et al., 2009). The sacculus has bilateral excitatory connections to the neck extensor muscle motoneurons and bilateral inhibitory connections to the neck flexor muscle motoneurons [START_REF] Uchino | Sacculocollic reflex arcs in cats[END_REF]. Thus, stimulation of the sacculus can modulate the tonic electromyogram (EMG) activity of the respective muscles. In the averaged EMG, such modulations emerge as sonomotor responses [START_REF] A C C E P T E D ; Jacobson | The vestibular evoked myogenic potential and other sonomotor evoked potentials[END_REF].

The sonomotor response currently receiving the most attention is the vestibular evoked myogenic potential (VEMP) recorded from surface electrodes over the sternocleidomastoid muscle. [START_REF] Colebatch | Myogenic potentials generated by a clickevoked vestibulocollic reflex[END_REF] were the first to record a response from this recording site. Clicks of 140 or 145 dB SPL were presented at a rate of 3/s while the subjects activated their neck flexors, and 512 stimulus-related epochs of surface EMG (sEMG) were averaged. The averaged sEMG, i.e. the VEMP, showed early peaks with mean latencies of 13 and 23 ms, respectively. Considering their polarity, they were denoted as p13 and n23. This pioneering work triggered not only numerous research studies, but also lead to continuously increasing clinical interest in the VEMP [see e.g. the

A c c e p t e d m a n u s c r i p t

-4 -reviews by [START_REF] Ferber-Viart | Vestibular evoked myogenic potentials in humans: a review[END_REF], [START_REF] Colebatch | Vestibular evoked potentials[END_REF], [START_REF] Welgampola | Characteristics and clinical applications of vestibular-evoked myogenic potentials[END_REF], [START_REF] Hamann | Vestibulär evozierte myogene Potentiale[END_REF], [START_REF] Rauch | Vestibular evoked myogenic potentials[END_REF][START_REF] A C C E P T E D ; Jacobson | The vestibular evoked myogenic potential and other sonomotor evoked potentials[END_REF][START_REF] A C C E P T E D ; Jacobson | The vestibular evoked myogenic potential and other sonomotor evoked potentials[END_REF], or the recent monograph by [START_REF] Murofushi | Vestibular Evoked Myogenic Potential -Its Basics and Clinical Applications[END_REF]]. The clinical interest results from the opinion that recording the VEMP is the only readily available method for unilaterally assessing the sacculus [START_REF] Clarke | Unilateral examination of utricle and saccule function[END_REF][START_REF] Wuyts | Vestibular function testing[END_REF].

The latencies of the early VEMP peaks are basically invariant, apart from a possible prolongation in the case of a retrolabyrinthine lesion [START_REF] Murofushi | Diagnostic value of prolonged latencies in the vestibular evoked myogenic potential[END_REF]. Thus, parameters of interest are mainly the amplitudes of the VEMP peaks and measures that depend on these amplitudes such as the VEMP threshold. However, the VEMP amplitude depends not only on the inhibition triggered by the acoustic stimulation of the vestibular end organs in the inner ear, but also on the level of the tonic muscle activation. [START_REF] Colebatch | Myogenic potentials generated by a clickevoked vestibulocollic reflex[END_REF] quantified the latter in terms of the level of the averaged rectified sEMG that was observed just before the stimulus presentation, and they found a highly linear relationship to the VEMP amplitude. Alternatively, the muscle tone may be quantified in terms of the pre-stimulus root-mean-square (RMS) value of the sEMG, which was reported to be linearly related to the VEMP amplitude as well [START_REF] Lim | The influence of voluntary EMG activity and click intensity on the vestibular click evoked myogenic potential[END_REF]. Subsequent studies corroborated these findings [START_REF] Ochi | Variance of vestibular-evoked myogenic potentials[END_REF][START_REF] Welgampola | Vestibular activation by bone conducted sound[END_REF][START_REF] Akin | The influence of voluntary tonic EMG level on the vestibular-evoked myogenic potential[END_REF]. Nowadays, the prevailing idea seems to be that the VEMP amplitude is proportional to the sEMG level. Many authors consequently normalize the VEMP amplitude by dividing it by a measure of the muscle tone (e.g. [START_REF] Karino | Independent effects of simultaneous inputs from the saccule and lateral semicircular canal. Evaluation using VEMPs[END_REF][START_REF] Welgampola | Characteristics and clinical applications of vestibular-evoked myogenic potentials[END_REF][START_REF] Miyamoto | Preliminary study on vestibular-evoked myogenic potential induced by bone-conducted stimuli[END_REF][START_REF] Brantberg | Age-related changes in vestibular evoked myogenic potentials[END_REF][START_REF] Sandhu | Effects of eye position on the vestibular evoked myogenic potential[END_REF][START_REF] Seo | Vestibular evoked myogenic potential induced by bone-conducted stimuli in patients with conductive hearing loss[END_REF][START_REF] Brantberg | Testing vestibular-evoked myogenic potentials with 90-dB clicks is effective in the diagnosis of superior canal dehiscence syndrome[END_REF]. But there is a problem. Proportionality between VEMP amplitude and sEMG level would imply that the signalto-noise ratio of the VEMP is independent of the muscle tone, which would be contradictory to the Answering this question requires a solid theoretical basis. The theory that will be developed here was inspired by the model of [START_REF] Wit | A simple model for the generation of the vestibular evoked myogenic potential (VEMP)[END_REF]. They simulated the VEMP generation by adding a large number of motor unit action potentials which randomly occurred in time with equal probability, apart from a narrow time window for which complete inhibition was assumed. The theory presented below is based on a similar model, but the mathematical framework is different:

Analytical formulas are derived that allow calculating the quantities of interest (mean and standard deviation of the sEMG) directly from the assumptions. The general theory requires the evaluation of integrals, but regarding a typical VEMP experiment the situation is less complicated:

Approximations can be used, owing to the fact that the inhibition time window may be assumed to be short (see e.g. [START_REF] Colebatch | Motor unit excitability changes mediating vestibulocollic reflexes in the sternocleidomastoid muscle[END_REF]. Thus, in the end, the question raised above can be answered on the basis of surprisingly simple formulas in an intuitively appealing way.

Model

Motor unit action potential and electromyogram

The elementary unit of function in the motor system is the motor unit, which consists of the cell body of the motor neuron, the axon of the motor neuron that runs in the peripheral nerve, the neuromuscular junction, and the muscle fibers innervated by that neuron [START_REF] Rowland | Diseases of the motor unit[END_REF]. When the motor neuron fires and an action potential reaches the neuromuscular junction, the muscle fibers depolarize and contract, and the currents associated with the depolarization spread in the surrounding volume conductor. If there were only one active motor unit, an electrode placed on the skin above the muscle would record a motor unit action potential (MUAP) whenever the unit fires.

A c c e p t e d m a n u s c r i p t

-6 -However, under realistic conditions many motor units are simultaneously active, and the contributions of the individual motor units can barely be recognized in the recorded signal, but are intermingled in an interference pattern [START_REF] Stegeman | Surface EMG models: properties and applications[END_REF]. This signal represents the surface electromyogram (sEMG).

The sEMG is the algebraic sum of the MUAP trains of all active motor units [START_REF] Day | Experimental simulation of cat electromyogram: evidence for algebraic summation of motor-unit action-potential trains[END_REF]. In the present study, the motor units are assumed to have the same MUAP time course, except for an amplitude factor that may differ from unit to unit. This simplifying assumption allows to integrate the MUAP trains of all motor units into a single MUAP train, where each MUAP is defined by a time of occurrence and an amplitude factor. Thus, the sEMG that is observed in a certain (finite) time range of interest may be calculated using the formula

) (t v T = ∑ = - T N i i i t t a 1 ) ( μ , (1)
where ) (t μ is the time course of a normalized MUAP occurring at time zero, and t i and a i are time of occurrence and amplitude, respectively, of the i-th MUAP. The parameter T N denotes the number of MUAPs that occur between -T/2 and +T/2, where T is assumed to be large enough to ensure that MUAPs occurring outside this range have no measurable effect on the time range of interest. The normalized MUAP, ) (t μ , is dimensionless, which means that the unit of measure (Volt) is attached to the amplitude factor i a . The function

) (t μ is assumed to fulfill the requirement 0 ) ( = ∫ +∞ ∞ - dt t μ (2)
and to have the value one at its maximum. The sample MUAP function shown in the inset of Fig. 1 (same curve also in the inset of Fig. 2A) corresponds to a single sine wave, defined as where Δ ˆ is the time span between negative and positive peak. In our simulations, Δ ˆ was 10 ms, corresponding to the interval between the peaks p13 and n23 of the VEMP. The zero crossing between the two peaks is considered as the occurrence time of the MUAP (marked by a vertical line). This special MUAP function was chosen because of its simplicity [a good alternative would have been the first derivative of a Gaussian probability density function, as in the model of [START_REF] Wit | A simple model for the generation of the vestibular evoked myogenic potential (VEMP)[END_REF]]. From a physiological point of view, such a MUAP function is, of course, relatively artificial. Therefore, it is emphasized that the theory developed here does in no way depend on Eq. ( 3) and could easily be applied to more realistic model functions or even MUAPs actually measured. As to the latter, it must be taken into account that a MUAP is usually measured under conditions that are not comparable to a VEMP measurement. Features such as the MUAP duration strongly depend not only on the recording conditions (e.g. near-versus far-field), but also on aspects such as muscle fiber length and conduction velocity [START_REF] Dumitru | Determinants of motor unit action potential duration[END_REF]. Moreover, Eq. ( 1) is based on the simplifying assumption that the time course of the MUAP is identical for all motor units. Thus, ) (t μ has to be interpreted as a function describing the time course of an average MUAP, which may be expected to resemble a compound muscle action potential (CMAP). The relationship between MUAP and CMAP was investigated, for example, by [START_REF] Dumitru | Physiologic basis of potentials recorded in electromyography[END_REF] and

McGill et al. (2001).

To calculate an sEMG using Eq. 1, random numbers have to be generated for T N , and for t i and a i (

T N i ≤ ≤ 1
). For the times of occurrence t i , a Poisson process with rate 0 ρ is assumed. This means that the random variable T N has the expected value T 0 ρ . To simulate a real experiment, where T N would vary from trial to trial, a random number has to be drawn from a Poisson distribution.

Regarding the amplitudes a i , [START_REF] Wit | A simple model for the generation of the vestibular evoked myogenic potential (VEMP)[END_REF] suggested to use a Gamma distribution of

A c c e p t e d m a n u s c r i p t

-8 -order 2. 1 The theory developed in the text that follows is more general, and specific assumptions about the distribution of the amplitudes are not required at this point.

Vestibular Evoked Myogenic Potential (VEMP)

With a simple extension, the above model can be transformed into a model for the VEMP. Instead of unconditionally summing over all the T N MUAPs in Eq. 1, some of the MUAPs are considered inhibited. The inhibition is controlled by an inhibition function ) (t ξ that specifies the probability of inhibition. 2 Equivalently, it can be postulated that the occurrence times of the MUAPs are drawn from a time-dependent Poisson process (see e.g. [START_REF] Cox | The theory of stochastic processes Chapman and Hall[END_REF]. This means that the probability of a MUAP between t and t+dt is assumed to be dt t) ( ρ , with a rate function (

)

) ( 1 ) ( 0 t t ξ ρ ρ - = . (4)
To prevent misunderstanding it is emphasized that, if not explicitly stated otherwise, the term 'rate' refers to the ensemble of all motor units.

1 In the simulations based on Eq. 1, Poisson-distributed random numbers N T and Gamma-distributed random numbers a i were drawn using the routines POISSRND and GAMRND, respectively, of the Matlab Statistical Toolbox (The MathWorks Inc., Natick, MA, USA). The MUAP occurrence times t i were obtained by mapping random numbers generated by the Matlab routine RAND (uniformly distributed on the interval (0,1)) onto the time range (-T/2, T/2).

2 Whether or not a specific MUAP in Eq. 1, say the i-th one, is considered inhibited depends on the value of

) ( i t ξ . A
value of zero means 'not inhibited', a value of one means 'inhibited'. For values in between, a random number is drawn from a uniform distribution on the interval (0,1). In our simulations, this was done using the Matlab routine -10 -curves in Fig. 2B show deflections similar to those in the theoretical curves (gray background). In Fig. 2A, however, the signal-to-noise ratio is insufficient.

Analytical considerations

Simulations based on Eq. 1 can be performed with little mathematical effort, but they are numerically expensive, because a huge number of trials have to be simulated to get estimates that approximate the expected values. Apart from that, simulations can only be done for a finite number of parameter constellations so that drawing general conclusions is problematic. This is why the model is considered analytically now. The consistency between the simulations based on Eq. 1 and the analytical model is demonstrated in Fig. 1, where the simulation results are represented by the white curves and the analytical model by the thick black curves in the background.

General theory

To prepare the transition from Eq. 1 to the analytical model, it shall be assumed that the time span are assumed to be drawn independently of each other and of the MUAP occurrence times t i . The probability that a MUAP amplitude is between a and a+da is assumed to be f(a)da. With these assumptions, the following equation is obtained for the expected value of ) (t v T in the limit T→∞:

) (t v = ∑ ∫ ∫ = - ∞ ∞ → - T N i i i i i T i i T T T da dt t t a N t a f 1 2 / 2 / 0 ) ( ) ( ) ( lim μ ρ . ( 5
)
Simplifying this equation yields

) (t v = ∫ +∞ ∞ - - ' ) ' ( ) ' ( dt t t t a μ ρ , ( 6
)
where

∫ ∞ = 0 ) ( da a f a a (7)
is the arithmetic mean of the MUAP amplitudes. A similar formula can be derived for the variance of the sEMG (see Appendix):

[ ]

) ( var t v = ∫ +∞ ∞ - - ' ) ' ( ) ' ( 2 2 dt t t t a μ ρ (8) with ∫ ∞ = 0 2 ) ( da a f a a . ( 9
)
The quantity a represents the quadratic mean of the MUAP amplitudes.

The equations 6 and 8 are more general than needed in this study. Taking into account that the rate function ) (t ρ has the structure specified in Eq. 4, the two equations may be rewritten as

) (t v = ∫ +∞ ∞ - - - ' ) ' ( ) ' ( 0 dt t t t a μ ξ ρ (10) and [ ] ) ( var t v = ( ) ) ( 1 2 0 t a Ξ - ⋅ ⋅ Δ ρ (
∫ +∞ ∞ - - Δ = Ξ ' ) ' ( ) ' ( 1 ) ( 2 dt t t t t μ ξ . ( 12
)
The quantity Δ in the latter two equations is defined as

Δ = ∫ +∞ ∞ - dt t) ( 2 μ . ( 13
)
It may be interpreted as a measure of the effective MUAP duration, because the integral basically measures how long the normalized MUAP has an absolute value close to the maximum value, i.e. one.

Approximations for short inhibition windows

Auditory stimuli that are used to elicit a VEMP are usually short, and the inhibition time window may be assumed to be short as well. Supposing that it is short enough to ensure that the MUAP function ) (t μ is roughly constant for the duration of the inhibition, Eq. 10 may be approximated as

) (t v ) ( 0 t aμ δ ρ ⋅ - ≈ (14) with δ = ∫ +∞ ∞ - ' ) ' ( dt t ξ . ( 15
)
The latter quantity shall be called the equivalent rectangular duration (ERD) of the inhibition, because it indicates how long a rectangular window with complete inhibition must be to produce the same inhibition effect as the actual inhibition function ) (t ξ . According to Eq. 14, the VEMP is a polarity-inverted MUAP, except for scaling issues. Its amplitude is proportional to the term δ ρ 0 ,

A c c e p t e d m a n u s c r i p t

-13 -which may be interpreted as the expected number of inhibited MUAPs. The approximation of ) (t Ξ for short inhibition windows is

) ( ) ( 2 t t μ δ Δ ≈ Ξ . ( 16
)
Thus, the variance of the sEMG shows a modulation which corresponds to the MUAP function squared. The ratio Δ / δ measures the ERD in units of the effective MUAP duration.

Approximations for a slowly changing inhibition function

Very small VEMP amplitudes (or no VEMPs at all) are to be expected if the inhibition function ) (t ξ changes much slower than the MUAP function ) (t μ . Equation 10 can be approximated in that case as 0

) ( ≈ t v
, due to Eq. 2. The corresponding approximation for Eq. 12 is

) ( ) ( t t ξ ≈ Ξ
. Thus, the variance of the sEMG shows a modulation which follows the time course of the inhibition.

Normalization of the VEMP

To normalize the VEMP, the recorded potential is divided by a measure of the muscle tone. Within the framework of the present theory, the most natural measure of the muscle tone is the standard deviation (or RMS value) of the uninhibited sEMG, which is a ⋅ Δ 0 ρ (see Eq. 11). Normalization of the approximation for short inhibition windows (Eq. 14) yields

) (t v norm ) ( 0 0 t a a μ ρ δ ρ ⋅ Δ - ≈ . ( 17
)
The first fraction on the right of this equation was intentionally not simplified (by canceling the factor 0 ρ ), because the terms in numerator and denominator have instructive interpretations: δ ρ 0 is the number of inhibited MUAPs, whereas Δ 0 ρ indicates with how many other MUAPs In what follows, the various factors will be considered one by one.

The dependence of the normalized VEMP on the effective MUAP duration, Δ , is not a major issue.

In theory, Δ could even be estimated from the data: If it is assumed that the MUAP function ) (t μ and the VEMP have a similar shape (see Eq. 14), the integral in Eq. 13 could be evaluated with a normalized MUAP function derived from the measured VEMP. But such a procedure may seriously fail if the assumption is violated. Thus, a model-based evaluation of Δ appears more appropriate. For the sine wave displayed in the inset of Fig. 1, the result is Δ = Δ ˆ = 10 ms, corresponding to the typical latency difference of the VEMP peaks p13 and n23.

The ratio a a / in Eq. 17 does not represent a serious problem either, because the lack of knowledge regarding the amplitude distribution of the MUAPs can be largely compensated for by making reasonable assumptions. If the MUAP amplitudes are drawn from a Gamma distribution of order k, the ratio a a / has the value ) 1 /( + k k

. For k = 2 (as in Fig. 1B) this value is 0.82, and with increasing k it approaches 1. Similar values can be expected for a real experiment, because the assumption of a Gamma distribution is considered appropriate for both normal and pathological conditions (e.g. [START_REF] Slawnych | Motor unit number estimation: sample size considerations[END_REF]. All in all, the amplitude distribution of the MUAPs seems to be of secondary importance with regard to VEMP modeling. Thus, as long as the goal is to achieve a qualitative understanding of the VEMP, it appears acceptable to assume MUAPs of unit amplitude, for the sake of simplicity. This yields 1 / = a a .

A c c e p t e d m a n u s c r i p t

-15 -

The actual problem is the dependence of the normalized VEMP on the MUAP rate, 0 ρ . The normalization reduced this problem, at least. Without normalization, the amplitude of the VEMP is proportional to 0 ρ , whereas after normalization it is proportional to 0 ρ . Supposing, for example, that the uncertainty regarding 0 ρ corresponds to a factor of 16, the uncertainty regarding and11). An increase of the MUAP rate from 400/s to 6400/s corresponds to a factor of 16. In the curves on the left, the VEMP amplitude increases by about the same factor, whereas in the normalized curves on the right there is only an increase by about a factor of 4. As to the standard deviation of the uninhibited sEMG, increasing the MUAP rate by a factor of 16 led to an increase by a factor of 4, whereas the normalized counterpart is always 1, by definition. Most interesting is the inhibition-induced modulation of the normalized standard deviation, which appears to be independent of the MUAP rate. To explain this observation, the variance of the sEMG (Eq. 11) has to be normalized. The normalized variance, i.e. the square of the normalized standard deviation, is

) ( 1 t Ξ -
. This term is indeed independent of the MUAP rate 0 ρ (see Eqs 12 and 16).

A c c e p t e d m a n u s c r i p t

-16 -

Relationship between VEMP amplitude and sEMG level

Both the VEMP amplitude and the standard deviation of the sEMG depend on the MUAP rate 0 ρ .

However, the latter is generally unknown when analyzing experimental data, and therefore it is useful to establish a more direct relationship between the two measures. The following consideration is based on the approximation for short inhibition windows, according to which the peak amplitude of the VEMP has an absolute value of

a s ⋅ = δ ρ 0 . ( 18
)
The standard deviation of the uninhibited sEMG is

a ⋅ Δ = 0 ρ σ . ( 19
)
The latter formula may be rewritten as

( )⋅ Δ = 2 2 0 / a σ ρ
Inserting this result into the equation for s yields:

2 2 a a s σ δ ⋅ Δ = (20)
Thus, the VEMP amplitude is proportional to the variance of the sEMG. If both quantities are appropriately normalized (division by a and 2 a , respectively), the factor of proportionality is the ratio of ERD and effective MUAP duration.

A c c e p t e d m a n u s c r i p t

-17 -

Discussion

Resolving a seeming contradiction

Starting point for this theoretical study was the realization that the common view of a proportional relationship between VEMP amplitude and sEMG level (in this study quantified in terms of the standard deviation of the sEMG, σ ) is contradictory to the experience that a successful VEMP recording requires a certain muscle tone. Equation 20 allows to resolve the contradiction, although this is not obvious at first glance. The importance of the muscle tone follows from the fact that, according to Eq. 20, σ / s is proportional to σ . This means the higher the muscle tone, the better is the signal-to-noise ratio of the VEMP. However, the equation appears to conflict with the view that VEMP amplitude and sEMG level are proportional to each other. Instead, it suggests that the VEMP amplitude is proportional to the sEMG level squared. This raises the question as to why contradictory conclusions were drawn from experimental data.

The pivotal element in the following consideration is the ERD, δ . So far, there was no reason to address the question as to how δ might depend on experimental conditions. But this question becomes crucially import now. By reanalyzing data published in the classical article by [START_REF] Colebatch | Myogenic potentials generated by a clickevoked vestibulocollic reflex[END_REF] it will be exemplified how the parameter δ can be used to reconcile model and experiment. Data points read from Fig. 4B of that article are displayed as filled circles in Fig. 4A.

A linear fit (thick gray line in the background) is quite satisfactory so that the data could be seen as counterevidence to the quadratic law suggested by Eq. 20. But such a conclusion would be premature. The dashed curve shows a quadratic function fitted to all data points with a mean sEMG level <70 μV (4 data points excluded). The function explains these data at least as well as the linear function.

A c c e p t e d m a n u s c r i p t

-18 -To remove the discrepancy between model and data at higher sEMG levels, the inhibition function ) (t ξ is assumed to depend on the MUAP rate 0 ρ . The inhibition function determines the proportion of MUAPs that are inhibited, and it appears reasonable to assume that the relative inhibition effect decreases with increasing excitation. In Eq. 20, the inhibition effect is not expressed in terms of the inhibition function ) (t ξ , but in terms of the ERD δ , which is an integral of the inhibition function (Eq. 15). In order to make the model consistent with the data, δ is assumed to be a function of σ :

( )

2 2 0 / 1 ref σ σ δ σ δ + = (21)
With this empirical extension, the model explains the data also at higher sEMG levels (black solid curve). The estimated reference sEMG level was ref σ = 37.5 μV. At low sEMG levels, Eq. 21 may be approximated as ( ) 0 δ σ δ ≈ , and the VEMP amplitude is proportional to 2 σ (dotted curve). At high sEMG levels, Eq. 21 may be approximated as ( )

σ σ δ σ δ / 0 ref ≈
so that the VEMP amplitude becomes proportional to σ . The latter approximation explains why the curve representing the model basically coincides, over a wide range of levels, with the linear fit published by [START_REF] Colebatch | Myogenic potentials generated by a clickevoked vestibulocollic reflex[END_REF]. The curve in Fig. 4B

Estimation of the ERD

The above considerations exemplified how to compare different experimental conditions by calculating relative ERDs. At least in theory, the data of a VEMP experiment allow to determine the ERD in absolute terms, too. This information is not contained in the classical VEMP curve (mean of the sEMG), but requires inspecting the standard deviation of the sEMG (as explained in the context of Fig. 3). Getting meaningful data for such an analysis may be a challenging task, though. The problem is that the standard deviation is affected to a much lesser degree by the inhibition than the mean (see Fig 2A). Thus, to make an analysis of the standard deviation possible, the measuring time has to be largely increased. This appears feasible in scientific investigations, where data could also be accumulated over multiple sessions and subjects, if necessary. Whether estimating the ERD could become an option for clinical investigations as well is difficult to assess at this point. Figure 2B suggests that not more than 200 sEMG trials are required to obtain clear evidence of an inhibition effect in the standard deviation. In a real experiment, however, an ERD of 5 ms, as assumed in that example, might be difficult to achieve.

Estimation of the MUAP rate

Given the ERD, another quantity of potential interest could be determined as well: the MUAP rate noise ratio of 6 (assuming that averaging n trials enhances the signal-to-noise ratio by n ), which seems to be a realistic order of magnitude for a typical VEMP experiment. Considering the inhibition effects found in peristimulus time histograms of motor unit discharges [START_REF] Colebatch | Motor unit excitability changes mediating vestibulocollic reflexes in the sternocleidomastoid muscle[END_REF], the assumption δ = 2 ms appears to be reasonable. With these assumptions, Eq.

22 predicts a MUAP rate of 1350/s. This value is consistent with an estimate derived by [START_REF] Wit | A simple model for the generation of the vestibular evoked myogenic potential (VEMP)[END_REF]. Referring to [START_REF] Keenan | Influence of amplitude cancellation on the simulated surface electromyogram[END_REF] they assumed that the number of motor units contributing to the sEMG recorded from the sternocleidomastoid muscle is between 50 and 100.

Moreover, referring to a compilation of MUAP rates in a review article by [START_REF] Enoka | Motor unit physiology: some unresolved issues[END_REF] they assumed a MUAP rate of 25/s per motor unit. For the ensemble of all contributing motor units they consequently estimated a MUAP rate between 1250/s and 2500/s.

Conclusions

Building upon the pioneering work of [START_REF] Wit | A simple model for the generation of the vestibular evoked myogenic potential (VEMP)[END_REF], an analytical model has been developed that explains the VEMP by the inhibition of MUAPs. An inhibition that is short compared to the duration of the VEMP (a reasonable assumption, considering e.g. the study of [START_REF] Colebatch | Motor unit excitability changes mediating vestibulocollic reflexes in the sternocleidomastoid muscle[END_REF]) can be characterized in terms of its ERD, irrespective of its actual time course. In theory, the ERD can be estimated from measured data by considering the standard deviation of the sEMG. But in practice, this is problematic because the inhibition effect on the standard deviation is small. Thus, the measuring time would have to be much longer than usual today.

For an ERD being independent of the sEMG level, the VEMP amplitude predicted by the model is proportional to the variance of the sEMG. The normalized VEMP amplitude, defined as the ratio of VEMP amplitude and standard deviation of the sEMG, is consequently proportional to the standard For that purpose, it has to be assumed that the ERD decreases with increasing sEMG level, which means that the relative inhibition effect (proportion of inhibited MUAPs) is assumed to diminish with increasing excitation. All in all, the model suggests that the optimal condition for measuring the VEMP is an intermediate muscle tone.

Increasing the muscle tone to the highest possible level would cause discomfort without significantly enhancing the signal-to-noise ratio.

Appendix: Variance of the sEMG

The variance of ) (t v T in the limit T→∞ may be calculated as [ ] ) where E[…] denotes the expected value. To calculate the first term on the right of this equation,

) ( var t v = [ ] ( ) 2 2 ) ( ) ( lim t v t v E T T - ∞ → , ( A1
) ( 2 t v T is expanded, yielding [ ] ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ - - ⋅ + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ - = ∑ ∑ ∑ - = + = = 1 1 1 1 2 2 2 ) ( ) ( 2 ) ( ) ( T T T T N i j N i j i j i N i i i t t t t a a E t t a E t v E μ μ μ . ( A2
∑ ∑ ∫ ∫ ∫ ∫ - = + = - - ∞ ∞ - - 1 1 1 2 / 2 / 2 / 2 / 0 0 ) ( ) ( ) ( ) ( ) ( ) ( 2 T T N i N i j T T j i j i j i j i j i T j T i T T da da dt dt t t t t a a a f a f N t N t μ μ ρ ρ , which can be rewritten as 2 2 / 2 / 2 2 ) ' ( ) ' ( ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - ∫ - T T T T T dt t t t a N N N μ ρ .
In the limit T→∞, T N is negligible compared to 2 T N , and, with Eq. 6, the above term reduces to

( ) 2 2 ) ( ' ) ' ( ) ' ( t v dt t t t a = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - ∫ +∞ ∞ - μ ρ .
Using this result, Eq. A1 may be rewritten as

)] ( var[ t v = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ - ∑ = ∞ → T N i i i T t t a E 1 2 2 ) ( lim μ . ( A3
)
The remaining calculation is analogous to the calculation of ) (t v . The result is Eq. 8.

  experience that a certain muscle tone is an indispensable requirement for a successful VEMP recording. This contradiction indicates that the relationship between VEMP amplitude and sEMG level is still insufficiently understood. The question arises as to what extent the proposed VEMP amplitude really eliminates the influence of the muscle tone, thus providing an unequivocal signature of the inhibition arising from the vestibular end organs.

  Figure 1A illustrates how the sEMG arises from an interference of MUAPs. The thin curves show

  interferes. The normalized VEMP is nothing else than a signal-to-noise ratio. The signal is proportional to the number of inhibited neurons, whereas the noise is proportional to the square root of the number of interfering MUAPs. Equation17shows that the normalization did not completely eliminate the dependence of the VEMP on aspects of the muscle activity.

  in Fig.3, which compares a situation before normalization (on the left) with a situation after normalization (on the right). Three different MUAP rates are considered: 400/s (top), 1600/s (middle), and 6400/s (bottom). Complete inhibition within a time window of 1.25 ms, centered at t = 0, was assumed. As in the simulations before, the MUAP function corresponded to the sine wave shown in the inset of Fig.1. The MUAPs had unit amplitude, for the sake of simplicity. The lower curve in each panel shows the mean of the sEMG, i.e. the VEMP; the upper curve shows the standard deviation of the sEMG. The thin curves are based on 10,000 sEMG epochs (simulation using Eq. 1), whereas the thick gray curves in the background are based on the analytical solutions (Eqs 10

  shows the function ( ) 0 /δ σ δ , which shall be called the relative ERD. If the muscle tone is high, the relative ERD is about five times smaller than if the muscle tone is low. The normalized VEMP is considered in Fig. 4C. Data points and curves were derived from the counterparts in Fig. 4A by dividing the p13-n23 amplitude by the mean sEMG activity. The normalized VEMP linearly increases with the sEMG activity if the latter is low (dotted line). At intermediate levels, where the majority of the data points are found, this increase slows down, and beyond a level of 60 μV the normalized VEMP is roughly independent of the sEMG activity.

  be used now to estimate the order of magnitude of the MUAP rate in a typical VEMP experiment. The assumptions Δ = 10 ms and ( ) that the average of 100 trials has a signal-to-

  deviation of the sEMG. Since the normalized VEMP amplitude may be interpreted as a signal-tonoise ratio, the model explains why a certain muscle tone is an indispensable requirement for a successful VEMP recording. The model is also able to explain why experimental data often suggest a more or less proportional relationship between VEMP amplitude and sEMG level (quantified, for example, by the RMS value of the sEMG).

)

  Evaluation of the second term in Eq. A2 yields

Fig. 1

 1 Fig. 1 Model explaining how surface electromyogram (sEMG) and vestibular evoked myogenic

Fig. 2 .

 2 Fig. 2. VEMPs for inhibition windows of 1.25 ms (A) and 5 ms (B), respectively. The thin curves

Fig. 3 :

 3 Fig. 3: VEMP (lower curve) and standard deviation of the sEMG (upper curve) before and after

Fig. 4 :

 4 Fig. 4: Reanalysis of the data in Colebatch et al.'s (1994) figure 4B. A VEMP amplitude

  Fig. 1

Fig. 2 A

 2 Fig. 2

Fig. 3 A

 3 Fig. 3

Fig. 4

 4 Fig. 4