Computer Simulation of the Vestibular SystemEdit
Model without CupulaEdit
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 that is perpendicular to the plane of the canal. We will also use the following notations:
- Rotation angle of tube [rad]
- Angular velocity of the tube [rad/s]
- Angular acceleration of the tube [rad/s^2]
- Rotation angle of the endolymph inside the tube [rad], and similar notation for the time derivatives
- 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 onto the plane of the semicircular canal described by to go from the 3D environment of the head to our scalar description. That is,
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:
- The inertial moment , where I characterizes the inertia of the endolymph.
- The viscous moment , caused by the friction of the endolymph on the walls of the tube.
This gives the equation of motion
Substituting and integrating gives
Let us now consider the example of a velocity step of constant amplitude . In this case, we obtain a displacement
and for , we obtain the constant displacement
Now, let us derive the time constant . Fora thin tube, , the inertia is approximately given by
From the Poiseuille-Hagen Equation, the force F from a laminar flow with velocity v in a thin tube is
where is the volume flow per second, the viscosity and the length of the tube.
With the torque and the relative angular velocity , substitution provides
Finally, this gives the time constant
For the human balance system, replacing the variables with experimentally obtained parameters yields a time constant of about 0.01 s. This is brief enough that in equation (10.5) the can be replaced by " = ". This gives a system gain of
Model with CupulaEdit
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
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
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
where "~" indicates the Laplace transformed variable. With from above, and defined by
we get the
For humans, typical values for 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:
Since , and since
Typically we are interested in the cupula displacement as a function of head velocity :
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!
For Linear, Time-Invariant systems (LTI systems), the input and output have a simple relationship in the frequency domain :
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:
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 = ; 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 =  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) mp.show()
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.
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:
where is the effective density of the material, the displacements along the i-axis, the i-component of the volume force, and the components of the Cauchy’s strain tensor. are the coordinates.
For linear elastic, isotropic material, Cauchy’s strain tensor is given by
where and are the Lamé constants; is identical with the shear modulus. , and is the stress tensor
This leads to Navier’s Equations of motion
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 the simultaneous decrease in diameter. The Lamé constants and are related to E and by:
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.
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.
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, . Based on this internal model, the expected firing rate of the internal estimate of the angular velocity, , 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 .
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.
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.
- Thomas Haslwanter (2013). "Vestibular Processing: Simulation of the Velocity Storage [Python"]. http://work.thaslwanter.at/CSS/Code/vestibular_feedback.py.