Last modified on 10 February 2015, at 19:32

Sensory Systems/Computer Models/Vestibular Simulation

Computer Simulation of the Vestibular SystemEdit

Semicircular CanalsEdit

Model without CupulaEdit

Simplified semicircular canal, without cupula.

Let us consider the mechanical description of the semi-circular canals (SCC). We will make very strong and reductive assumptions in the following description. The goal here is merely to understand the very basic mechanical principles underlying the semicircular canals.

The first strong simplification we make is that a semicircular canal can be modeled as a circular tube of “outer” radius R and “inner” radius r. (For proper hydro mechanical derivations see (Damiano and Rabbitt 1996) and Obrist (2005)). This tube is filled with endolymph.

The orientation of the semicircular canal can be described, in a given coordinate system, by a vector  \vec n that is perpendicular to the plane of the canal. We will also use the following notations:

 \theta Rotation angle of tube [rad]
 \dot{\theta} \equiv \frac{d \theta}{dt} Angular velocity of the tube [rad/s]
 \ddot{\theta} \equiv \frac{d^2 \theta}{dt^2} Angular acceleration of the tube [rad/s^2]
 \phi Rotation angle of the endolymph inside the tube [rad], and similar notation for the time derivatives
 \delta = \theta - \phi movement between the tube and the endolymph [rad].

Note that all these variables are scalar quantities. We use the fact that the angular velocity of the tube can be viewed as the projection of the actual angular velocity vector of the head  \vec \omega onto the plane of the semicircular canal described by  \vec n to go from the 3D environment of the head to our scalar description. That is,

 \dot{\theta} = \vec \omega \cdot \vec n

where the standard scalar product is meant with the dot.

To characterize the endolymph movement, consider a free floating piston, with the same density as the endolymph. Two forces are acting on the system:

  1. The inertial moment  I \ddot{\phi} , where I characterizes the inertia of the endolymph.
  2. The viscous moment  B \dot{\delta} , caused by the friction of the endolymph on the walls of the tube.

This gives the equation of motion

 I \ddot{\phi} = B \dot{\delta}

Substituting  \phi = \theta - \delta and integrating gives

 \dot{\theta} = \dot{\delta} + \frac{B}{I} \delta .

Let us now consider the example of a velocity step  \dot{\theta}(t) of constant amplitude  \omega . In this case, we obtain a displacement

 \delta = \frac{I}{B} \omega \cdot (1-e^{-\frac{B}{I}t})

and for  t \gg \frac{I}{B} , we obtain the constant displacement

 \delta \approx \frac{I}{B} \omega .

Now, let us derive the time constant  T_1 \equiv \frac{I}{B} . Fora thin tube,  r \ll R , the inertia is approximately given by

 I = m l^2 \approx 2 \rho \pi^2 r^2 R^3 .

From the Poiseuille-Hagen Equation, the force F from a laminar flow with velocity v in a thin tube is

 F = \frac{8 \bar{V} \eta l}{r^2}

where  \bar{V} = r^2 \pi v is the volume flow per second,  \eta the viscosity and  l = 2 \pi R the length of the tube.

With the torque  M = F \cdot R and the relative angular velocity  \Omega = \frac{v}{R} , substitution provides

 B = \frac{M}{\Omega} = 16 \eta \pi ^2 R^3

Finally, this gives the time constant  T_1

 T_1 = \frac{I}{B} = \frac{\delta r^2}{8 \eta}

For the human balance system, replacing the variables with experimentally obtained parameters yields a time constant  T_1 of about 0.01 s. This is brief enough that in equation (10.5) the  \approx can be replaced by " = ". This gives a system gain of

 G \equiv \frac{\delta}{\omega} = \frac{I}{B} = T_1

Model with CupulaEdit

Effect of the cupula.

Our discussion until this point has not included the role of the cupula in the SCC: The cupula acts as an elastic membrane that gets displaced by angular accelerations. Through its elasticity the cupula returns the system to its resting position. The elasticity of the cupula adds an additional elastic term to the equation of movement. If it is taken into account, this equation becomes

 \ddot{\theta} = \ddot{\delta} + \frac{B}{I} \dot{\delta} + \frac{K}{I} \delta

An elegant way to solve such differential equations is the Laplace-Transformation. The Laplace transform turns differential equations into algebraic equations: if the Laplace transform of a signal x(t) is denoted by X(s), the Laplace transform of the time derivative is

 \frac{dx(t)}{dt} \xrightarrow{Laplace Transform} s \cdot X(s) - x(0)

The term x(0) details the starting condition, and can often be set to zero by an appropriate choice of the reference position. Thus, the Laplace transform is

 s^2 \tilde{\theta} = s^2 \tilde{\delta} + \frac{B}{I} s \tilde{\delta} + \frac{K}{I} \tilde{\delta}

where "~" indicates the Laplace transformed variable. With  T_1 from above, and  T_2 defined by

 T_2 = \frac{B}{K}

we get the

 \frac{ \tilde{\delta} }{ \tilde{\theta} } = \frac{T_1 s^2}{T_1 s^2 + s + \frac{1}{T_2}}

For humans, typical values for  T_2 = B/K are about 5 sec.

To find the poles of this transfer function, we have to determine for which values of s the denominator equals 0:

 s_{1,2} = \frac{1}{T_1} \Big(-1 \pm \sqrt{1-4\frac{T_1}{T_2}} \Big)

Since  T_2 \gg T_1 , and since

 \sqrt{1-x} \approx 1 - \frac{x}{2} for x \ll 1

we obtain

 s_1 \approx - \frac{1}{T_1}, and s_2 \approx - \frac{1}{T_2}

Typically we are interested in the cupula displacement  \delta as a function of head velocity  \dot{\theta} \equiv s \tilde{\theta} :

 \frac{\tilde{\delta}}{s \tilde{\theta}}(s) = \frac{T_1 T_2 s}{(T_1 s +1)(T_2 s + 1)}

For typical head movements (0.2 Hz < f < 20Hz), the system gain is approximately constant. In other words, for typical head movements the cupula displacement is proportional to the angular head velocity!

Bode plot of the cupula displacement of a function of head velocity, with T1 = 0.01 sec, T2 = 5 sec, and an amplification factor of (T1+ T2)/ (T1* T2) to obtain a gain of approximately 0 for the central frequencies.

Control SystemsEdit

For Linear, Time-Invariant systems (LTI systems), the input and output have a simple relationship in the frequency domain :

 Out(s) = G(s)*In(s)

where the transfer function G(s) can be expressed by the algebraic function


In other words, specifying the coefficients of the numerator (n) and denominator (d) uniquely characterizes the transfer function. This notation is used by some computational tools to simulate the response of such a system to a given input.

Different tools can be used to simulate such a system. For example, the response of a low-pass filter with a time-constant of 7 sec to an input step at 1 sec has the following transfer function


and can be simulated as follows:

With SimulinkEdit
Step-response simulation of a lowpass filter with Simulink.

If you work on the command line, you can use the Control System Toolbox of MATLAB or the module signal of the Python package SciPy:

MATLAB Control System Toolbox:

% Define the transfer function
num = [1];
tau = 7;
den = [tau, 1];
mySystem = tf(num,den)
% Generate an input step
t = 0:0.1:30;
inSignal = zeros(size(t));
inSignal(t>=1) = 1;
% Simulate and show the output
[outSignal, tSim] = lsim(mySystem, inSignal, t);
plot(t, inSignal, tSim, outSignal);

Python - SciPy:

# Import required packages
import numpy as np
import scipy.signal as ss
import matplotlib.pylab as mp
# Define transfer function
num = [1]
tau = 7
den = [tau, 1]
mySystem = ss.lti(num, den)
# Generate inSignal
t = np.arange(0,30,0.1)
inSignal = np.zeros(t.size)
inSignal[t>=1] = 1
# Simulate and plot outSignal
tout, outSignal, xout = ss.lsim(mySystem, inSignal, t)
mp.plot(t, inSignal, tout, outSignal)


Consider now the mechanics of the otolith organs. Since they are made up by complex, visco-elastic materials with a curved shape, their mechanics cannot be described with analytical tools. However, their movement can be simulated numerically with the finite element technique. Thereby the volume under consideration is divided into many small volume elements, and for each element the physical equations are approximated by analytical functions.

FE-Simulations: Small, finite elements are used to construct a mechanical model; here for example the saccule.

Here we will only show the physical equations for the visco-elastic otolith materials. The movement of each elastic material has to obey Cauchy’s equations of motion:

 \rho \frac{\partial^2 u_i}{\partial t^2} = \rho B_i + \sum_{j} \frac{\partial T_{ij}}{\partial x_j}

where  \rho is the effective density of the material,  u_i the displacements along the i-axis,  B_i the i-component of the volume force, and  T_{ij} the components of the Cauchy’s strain tensor.  x_j are the coordinates.

For linear elastic, isotropic material, Cauchy’s strain tensor is given by

 T_{ij} = \lambda e \delta_{ij} + 2 \mu E_{ij}

where  \lambda and  \mu are the Lamé constants;  \mu is identical with the shear modulus.  e = div(\vec u) , and  E_{ij} is the stress tensor

 E_{ij} = \frac{1}{2} \Big( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \Big).

This leads to Navier’s Equations of motion

 \rho \frac{\partial ^2 u_i}{\partial t^2} = \rho B_i + (\lambda + \mu) \frac{\partial e}{\partial x_i} + \mu \sum_{j} \frac{\partial ^2 u_i}{\partial x_j^2}

This equation holds for purely elastic, isotropic materials, and can be solved with the finite element technique. A typical procedure to find the mechanical parameters that appear in this equation is the following: when a cylindrical sample of the material is put under strain, the Young coefficient E characterizes the change in length, and the Poisson’s ratio  \nu the simultaneous decrease in diameter. The Lamé constants  \lambda and  \mu are related to E and  \nu by:

 E = \frac{\mu (3 \lambda + 2 \mu)}{\lambda + \mu}


 \nu = \frac{\lambda}{2(\lambda + \mu)}

Central Vestibular ProcessingEdit

Central processing of vestibular information significantly affects the perceived orientation and movement in space. The corresponding information processing in the brainstem can often be modeled efficiently with control-system tools. As a specific example, we show how to model the effect of velocity storage.

Velocity StorageEdit

The concept of velocity storage is based on the following experimental finding: when we abruptly stop from a sustained rotation about an earth-vertical axis, the cupula is deflected by the deceleration, but returns to its resting state with a time-constant of about 5 sec. However, the perceived rotation continues much longer, and decreases with a much longer time constant, typically somewhere between 15 and 20 sec.

Vestibular Modeling: The blue curve describes the deflection of the cupula as a response to a velocity step, modeled as a high-pass filter with a time-constant of 5 sec. The green curve represents the internal estimate of the angular velocity, obtained with an internal model of the cupula-response in a negative feedback look, and a feed-forward gain-factor of 2.

In the attached figure, the response of the canals to an angular velocity stimulus ω is modeled by the transfer function C, here a simple high-pass filter with a time constant of 5 sec. (The canal response is determined by the deflection of the cupula, and is approximately proportional to the neural firing rate.) To model the increase in time constant, we assume that the central vestibular system has an internal model of the transfer function of the canals, \hat{C}. Based on this internal model, the expected firing rate of the internal estimate of the angular velocity, \hat{\omega}, is compared to the actual firing rate. With a the gain-factor k set to 2, the output of the model nicely reproduces the increase in the time constant. The corresponding Python code can be found at [1].

It is worth noting that this feedback loop can be justified physiologically: we know that there are strong connections between the left and right vestibular nuclei. If those connections are severed, the time constant of the perceived rotation decreases to the peripheral time-constant of the semicircular canals.

Central Vestibular Processing can often be described with control-system models. Here "omega" is the head velocity, "C" the transfer function of the semicircular canals, and "k" a simple gain factor. The "hat"-ed variables indicate internal estimates.

Mathematically, negative feedback with a high gain has the interesting property that it can practically invert the transfer function in the negative feedback loop: if k>>1, and if the internal model of the canal transfer function is similar to the actual transfer function, the estimated angular velocity corresponds to the actual angular velocity.

  & \hat{\omega }=(\omega C-\hat{\omega }\hat{C})k \\ 
 & \hat{\omega }(1+\hat{C}k)=\omega Ck \\ 
 & \frac{{\hat{\omega }}}{\omega }=\frac{C}{1/k+\hat{C}}\,\,\xrightarrow[if\,C\approx \hat{C}]{k>>1}1  


Auditory_System · Somatosensory_System