Sensory Systems/Computer Models/Vestibular Simulation

Auditory_System_Simulation · Somatosensory_System_Simulation

Computer Simulation of the Vestibular System edit

Semicircular Canals edit

Model without Cupula edit

 
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   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:

  1. The inertial moment  , where I characterizes the inertia of the endolymph.
  2. 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 Cupula edit

 
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

 

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

 

we obtain

 

Typically we are interested in the cupula displacement   as a function of head velocity  :

 

For typical head movements (0.2 Hz < f < 20 Hz), 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 Systems edit

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:

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

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)
mp.show()

Otoliths edit

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:

 

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:

 

and

 

Central Vestibular Processing edit

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 Storage edit

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,  . 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 .[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.

 

References edit

  1. Thomas Haslwanter (2013). "Vestibular Processing: Simulation of the Velocity Storage [Python]".