On the Solution of the Cupula Equation using Fractional Calculus.

Moustafa ElShahed

The response of the semicircular canal is examined theoretically. The mathematical modeling of cupula using fractional calculus has been shown to be very effective in representing the response of the semicircular canal. With the help of Laplace transformation and MATHEMATICA an exact solution is obtained for the equation of motion.

The semicircular canals are the primary organs for the sensing of angular or rotary motion. They are part of the organs of equilibrium contained in the non-auditory portion of the inner ear. Within the swollen region of each semicircular canal in the skull, there is a flap valve called the cupula,which detects motion of the endolymph-the viscous fluid in the canal. When our head moves,the endolymph deflects the cupula.This movement of the cupula is recognized by hair cells and then transmitted to the brain. For this reason, understanding the mechanics of the semicircular canal system with respect to head motion would be extremely useful.

The three semicircular ducts in each vestibular apparatus, known respectively as the anterior, posterior, and horizontal semicircular ducts, are arranged at right angles to each other so that they represent all three planes in space. Each semicircular duct has an enlargement at one of its ends called the ampulla and the ducts are filled with a viscous fluid called endolymph.The ampulla contains the cupula, a gelatinous dividing partition with the same density as endolymph. [Guyton 1992 P. 414].

When the head suddenly begins to rotate in any direction, the endolymph in the membranous semicircular ducts, because of its inertia, tends to remain stationary while the semicircular ducts themselves turn. This causes relative fluid flow in the ducts in a direction opposite to the rotation of the head; this relative flow lasts for a few seconds, thus activating the semicircular duct receptors; and the person perceives this onset of rotation, which is called angular acceleration. When the rotation suddenly stops, exactly the opposite effects take place: That is, the endolymph continues to rotate while the semicircular duct itself stops, and the duct receptors are activated in the opposite direction; the person perceives the termination of rotation [Cuyton 1992 P. 416].
The main idea of  this paper is the mathematical study of these phenomena and the obtaining of analytical solution for the equation of motion.

2-Mathematical Formulation

The most popular model for the response of the semicircular canals to the angular acceleration is due to Wilhelm Steinhausen,who in 1933 proposed that the canals respond to angular acceleration in the same manner as would heavily damped torsion pendulum. The "torsion pendulum equation" which he proposed may be written in the form [Elshehawey et. al. 2000]:

(x^″) (t) + (B/J) (x^′) (t) + (K/J) x (t)=f (t) .        (1)

Where x(t), (x^′) (t), and (x^″) (t) represent respectively the mean angular displacement, mean velocity and mean angular acceleration of the endolymph relative to the duct and f(t) is the component of the angular acceleration of the head perpendicular to the plane of the canal. The coefficient J is the inertia term defined by the mass and distribution of the endolymph. The damping coefficient B denotes the ratio of the torque resulting from viscous forces to the mean angular velocity of the endolymph. The stiffness of the cupula is represented by K.

There may also be viscous forces present in the cupula gelatinous layer, Due to this potential viscoelastic nature of the cupula material. The fractional derivative approach provides a powerful tool for modeling systems with damping materials. If the fractional derivative model is used to represent the damping characteristics, the equation of motion as discussed by [Suarez and Shokooh 1995]can be written in the form:

(Subscript[D, t]^2) x (t) + (B/J) (Subscript[D, t]^α) 
          x (t) + (K/J) x (t) = 
      f (t)        (2)
Note that the operator (Subscript[D, t]) of the old model was replaced by the fractional derivative operator  (Subscript[D, t]^α)  where α can now take on any value.

The term fractional calculus is by no means new. It is a generalization of the ordinary differentiation by non-integer derivatives. The subject is as old as calculus itself and goes back to times when Leibniz, Gauss, and Newton invented differentiation as a new tool in a calculation. In a letter to L`Hospital in 1695 Leibniz raised the following question: 

There are several definitions of derivative of fractional order.The fractional  operator (Subscript[D, t]^α)  based on Riemann- Liouville integral  is defined as (Miller and Ross, 1993):

\(D\_t\%α\ x \((t)\) = \(1\/\(Γ \((n - α)\)\)\) \(d\^n\/dt\^n\) \((∫\_0\%t\(\( 
                      x \((τ)\)\)\/\((t - τ)\)\^\(α - n + 1\)\) dτ)\), \ 
      n - 1 < α < n . \ \ \ \ \ \ \ \ \ \ \ \((3)\)\)
where Γ is the gamma function. Also the Caputo's definition can be written as :

\(D\_t\%α\ x \((t)\) = \(1\/\(Γ \((n - α)\)\)\) \((∫\_0\%t\(\((t - τ)\)\^\(n - α - 
                      1\)\) \(\(d\^n\) x \((τ)\)\)\/dt\^n)\) dτ, 
      n - 1 < α < n . \ \ \ \ \ \ \ \ \ \ \ \ \((4)\)\)

It is easy to observe that
                \(Lim\+\(α -> n\)\ D\_t\%α\ x \((t)\) = \(\(d\^n\) x \((t)\)\)\/dt\^n, \ \ \ n = 1, 2, \(\(...\) \(...\)\)\)

Constitutive equations, giving the relations between the uni-directional stress σ  and the respective strain ϵ, for an elastic material and a viscous material are as follows:
          σ (t)=k  .   ϵ (t)=k  .  (D^0)ϵ(t)               (5)
                 \(σ \((t)\) = \(k . \(\(d\ ϵ \((t)\)\)\/\(d\ t\)\) = k\ \  . \ \ \ D\^1  ϵ \((t)\)\ \ \ \ \ \ \ \ \ \ \ \((6)\)\)\)

Equations (5) and (6) are the well known Hook's law of elasticity and Newton's law of viscosity, respectively. The symbol (D^α) appearing in the right -hand side of the above equations denotes differentiation of arbitrary order α with respect to time.

Since the cases α = 0 and α = 1 correspond to purely elastic and purely viscous materials, respectively, it seems reasonable to assume that the equation: 
           σ (t)=k  .  (D^α)ϵ(t)               (7)           
 which is a straightforward generalization of equation (5) and (6), can represent the constitutive law for viscoelastic materials [Stiassnie 1979]. According to the definition of the fractional derivative of order α, equation (7) is rewritten as:

σ (t)=   (k/(Γ (1 - α))) (d/dt) ((∫_0^t)(( ϵ (τ))/(((t - τ))^α)) \[DifferentialD]τ),0<α<1.

3-Analysis

One of the most efficient and elegant methods to solve equation (2) is by means of the Laplace transform.

The formula for Laplace Transform of Riemann- Liouville fractional derivative is

\(∫\_0\%∞\( e\^\(-st\)\) {D\_t\%α\ x \((t)\)} dt = \(s\^α\) 
              x\ \((s)\) - ∑\+\(k = 0\)\%\(n - 1\)\(s\^k\) 
                D\_t\%\(α - k - 1\)\ x \((0)\)\),    (n - 1 < α < n)

The formula for Laplace Transform of Caputo fractional derivative is

\(∫\_0\%∞\( e\^\(-st\)\) {D\_t\%α\ x \((t)\)} dt = \(s\^α\) 
              x\ \((s)\) - ∑\+\(k = 0\)\%\(n - 1\)s\^\(α - k - 1\)\ \(x\^k\) \((0)\)\),     (n - 1 < α < n)

The definition of the fractional differentiation of  Riemann- Liouville played an important role in the development of the theory of fractional derivatives and integrals and for its applications  in pure mathematics. But applied problems require definitions of fractional derivatives allowing the utilization of physically interpretable initial conditions, which contain x(0), x'(0),... etc .The Laplace transform of a fractional derivative assumes the following simple form[Igor 1999]:  
                        L [(Subscript[D, t]^α) x (t)] = (s^α) x (s).                             (8)
   \(assuming\ zero\ initial\ displacement\ and\ \(\(velocity\)\(.\)\)\)

\(Employing\ the\ Laplace\ transform\ to\ equation \((2)\)\ in\ the\ light\ of\ \((8)\)\ we\ obtain\)
        X(s)=((f (s))/((s^2) + (B/J) (s^α) + (K/J))) .                                        (9)
We will examine the response of the canals to the unit impulse excitation δ(t) .  Then eq (9) takes the form:
 X(s)=(1/((s^2) + (B/J) (s^α) + (K/J))) .                                        (10)

Equation (5) can be written in the form [Igor 1999, p.155]:

X (s)=(∑_(k = 1)^∞)\(\(\(\((\(-1\))\)\^k\) \(\((K\/J)\)\^k\) 
                  s\^\(\(-α\)\ k - α\)\)\/\((s\^\(2 - α\)\  + B\/J)\)\^\(k + 1\)\) .              (11)

The inversion of  equation (10) takes the form:

\(\(\ \)\(x \((t)\) = \ ∑\+\(i = 1\)\%∞\(\((\(-1\))\)\^k\/\(k!\)\) \(\((K\/J)\)\^k\) 
            t\^\(2\ k + 1\)\ \ \ \(E\_\(α, β\)\%k\) \((\(-\(B\/J\)\) t\^\(2 - α\))\)\ \ \ \ \ \ \ \ \ \ \ \ \ \((12)\)\)\)

Where (Subscript[E, α,β ])(z) is the two parameter Mittage-Leffler function defined by:

\(\(E\_\(α, β\)\) \((z)\) = \ ∑\+\(i = 1\)\%∞ z\^k\/\(Γ \((αk + β)\)\), \ \ \ \ \ \ \(\((α > 0, \ \ \ \ \ \ β > 0)\)\(.\)\(\ \ \ \ \ \ \)\)\)

It has been shown that for most damping materials α =(1/2) is quite satisfactory Begley[1986].Then equation (10) takes the form:

X (s)=(1/((s^2) + (B/J) (Sqrt[s]) + (K/J))) .                                (13)

In order to obtain the inversion of Eq (13) we first write the denominator as: 
                    D (s) = (s^2) + (B/J) (Sqrt[s]) + (K/J), 
   we also set  s=(λ^2)   to obtain
                          P (λ) = (λ^4) + (B/J) λ + 
            (K/J). 
  Denoting the roots of  P(λ) by (Subscript[λ, 1]),  \(TraditionalForm\`\(\(λ\_2\)\(,\)\(\ \ \)\(λ\_3\)\(,\)\(λ\_4\)\(\ \)\)\). 
    Here, it may be pointed that [Suarez and Shokooh (1995-1997)] , [Sakakbira (1997)] and [Rossikhin and Shitikova 1999-2001] used a difficult methods to obtain the roots of the characteristic polynomial P(λ). However, in the present paper we have been able to use an easy method. We use the MATHEMATICA program to obtain the roots and we get::

\(λ\_1 = \(1\/2\) Ψ - \(1\/2\) \@\(\(-Ψ\^2\) - \(2  B\)\/\(J\ Ψ\)\)\),
                              \(λ\_2 = \(1\/2\) Ψ + \(1\/2\) \@\(\(-Ψ\^2\) - \(2  B\)\/\(J\ Ψ\)\)\),
                               \(λ\_3 = \(1\/2\) Ψ - \(1\/2\) \@\(\(-Ψ\^2\) + \(2  B\)\/\(J\ Ψ\)\)\),
                                (Subscript[λ, 4])=(1/2) Ψ+(1/2)\(\@\(\(-Ψ\^2\) + \(2  B\)\/\(J\ Ψ\)\)\).
Where    N=4 ((((2/3)))^((1/3))) (K/(J^((1/3)))),  M = (2^((1/3))) (3^((2/3))) (J^((2/3))),
    H=\(\((9 
                           B\^2 + \(\@\(\(3\)\(\ \)\)\) \@\(27  B\^4 - 256  J\ K\^3\))\)\^\(1\/3\)\)and    Ψ=(Sqrt[(N/H) + (H/M)]),

Then  equation (13) takes the form:
       X (s)=(∑_(i = 1)^4)(J/(((Sqrt[s]) - (Subscript[λ, i])) (B + 
                      4  J (Subscript[λ, i]^3)))).           (14)

Then the inversion of equation (14) gives:
 x (t)=(∑_(i = 1)^4)\(\(J\/\((B + 4  J\ λ\_i\%3)\)\)[
              1\/\@\(t\ π\) + \(λ\_i\) e\^\(t\ λ\_i\%2\)\ Erfc \((\(-λ\_i\)\ \@t)\)]\).      (15)

Secondly,we examine the response of the canals to a constant angular acceleration .  If  f(t)=A,a constant,then equation(9) takes the form:
X (s) =(J/(K s))+(1/K)(∑_(i = 1)^4)\(\(\(-B\)\ J - \(J\^2\) λ\_i\%3\)\/\(\((\@s - λ\_i)\) \((B + 4  J\ λ\_i\%3)\)\)\) .                (16)

Then the inversion of equation (16) takes the form:
x (t)=
A((J/K)+(1/K)  \(∑\+\(i = 1\)\%4\(\(\(-B\)\ J - \(J\^2\) λ\_i\%3\)\/\((B + 4  J\ λ\_i\%3)\)\)[
                      1\/\@\(t\ π\) + \(λ\_i\) e\^\(t\ λ\_i\%2\)\ Erfc \((\(-λ\_i\)\ \@t)\)]\) ).(17)

Now we will examine the response of the semicircular canals to a sinusoidal input of the form 
                                                   f(t) =A cos ω t , 
A is  constant. Then equation (9) takes the form:

\(X\ \((s)\) = \(A\/R\) \((BJω\^2\/\(s\^2 + ω\^2\) + \(s \((J\ K + \(J\^2\) ω\^4)\)\)\/\(s\^2 + ω\^2\))\) + 
∑\+\(i = 1\)\%4\( A \((\(J\^3\) \(ω\^6\) λ\_i - \(J\^3\) \(ω\^4\) λ\_i\%3 + \(J\^2\) \(ω\^2\) \(λ\_i\) \((K - B\ λ\_i)\) - \ J\ K \((B + J\ λ\_i\%3)\))\)\)\/\(\((B + 4  J\ λ\_i\%3)\) \((\@s - λ\_i)\)\)\)
                                                                                                                         (18)
where                             R=((K^2) + (B^2) (ω^2) + 
              2  J K(ω^4) + (J^2) (ω^8)).

Then the inversion of equation (18) takes the form:
  x (t)=\(\(A\/R\) \((BJω\ \ sin\ ω\ t + \((J\ K + \(J\^2\) ω\^4)\)\ cos\ ω\ t)\)\ \ \  + \ 
∑\+\(i = 1\)\%4\(\( 
                      A \((\(J\^3\) \(ω\^6\) λ\_i - \(J\^3\) \(ω\^4\) λ\_i\%3 + \(J\^2\) \(ω\^2\) \(λ\_i\) \((K - B\ λ\_i)\) - \ J\ K \((B + J\ λ\_i\%3)\))\)\)\/\((B + 
                        4  J\ λ\_i\%3)\)\)
\ \ \ [
                1\/\@\(t\ π\) + \(λ\_i\) e\^\(t\ λ\_i\%2\)\ Erfc \((\(-λ\_i\)\ \@t)\)] . \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \((19)\)\)

4. Conclusions

The response of the semicircular canal is examined theoretically. The mathematical modeling of cupula using fractional calculus has been shown to be very effective in representing the response of the semicircular canal. With the help of Laplace transformation and MATHEMATICA an exact solution is obtained for the equation of motion.

5. References

[1] Bagley,R. L., and Torvic, P. J.,  A theoritical Basis for the Application of Fractional Calculus to viscoelasticity. Journal of Rheology. 27,201-210,(1983).

[2]Igor Podlubny., Fractional Differential Equations,Academic press, (1999)

[3] Gerd Baumann,Knocking on Leibniz's Door,Mathematica in Education and Research, 8, 58-65,(1999).

[4] Guyton, A. C.,Human Physiology and Mechanisms of Disease.,W.B.Saunders Company, Philadelphia (1992)P.414-419.

[5] K.S. Miller and B. Ross, An Introduction to the Fractional Calculus and Fractional Differential Equations, John Wiley & Sons, Inc., New York, (1993).

[6] K.B. Oldham and J. Spanier, The Fractional Calculus, Academic Press, New York, (1974).

[7] Suarez. L and Shokooh. L,  Response of systems with damping materials modeled using fractional calculus, Appl Mech Rev, 48, S118-S126, (1995).