Sensory Systems/Computer Models/Vestibular Simulation

Sensory Systems

Vision
Hearing
Balance
Feeling
Smell
Taste

Technological Aspects

Implants
Models

In Animals

Birds
Fish
Marine Animals
Arthropods
Other Animals

Computer Simulation of the Vestibular System

Semicircular Canals

Model without Cupula

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=ml^{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 Cupula

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 {LaplaceTransform}}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}}forx\ll 1$

we obtain

$s_{1}\approx -{\frac {1}{T_{1}}},ands_{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 < 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

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

$G(s)={\frac {num(s)}{den(s)}}={\frac {n(0)*{{s}^{0}}+n(1)*{{s}^{1}}+n(2)*{{s}^{2}}+...}{d(0)*{{s}^{0}}+d(1)*{{s}^{1}}+d(2)*{{s}^{2}}+...}}$

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

$G(s)={\frac {1}{7s+1}}$

and can be simulated as follows:

With Simulink

Step-response simulation of a lowpass filter with Simulink.

Commandline

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

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 }}$

and

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

Central Vestibular Processing

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

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 .

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.

${\begin{aligned}&{\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\end{aligned}}$