C. elegans locomotion in three dimensions

The mechanical model

• viscoelastic rod
• fixed length (inextensibility)
• active bending moment
• environment drag

The muscles

Framed curve

Frame equations

Any (smooth) frame \((\vec{\tau}, \vec{e}_1, \vec{e}_2)\) defines three scalar fields \(\alpha, \beta\) and \(\gamma\) which satisfy:

\[ \frac{1}{|\vec{x}_u|} \vec{\tau}_u = \alpha \vec{e}^1 + \beta \vec{e}^2, \\ \frac{1}{|\vec{x}_u|} \vec{e}^1_u = -\alpha \vec{\tau} + \gamma \vec{e}^2, \\ \frac{1}{|\vec{x}_u|} \vec{e}^2_u = -\beta \tau - \gamma \vec{e}^1. \]

This corresponds to (vector) curvature \(\vec{\kappa} = \alpha \vec{e}^1 + \beta \vec{e}^2\) and twist \(\gamma\).

Darboux vector

• We denote by \(\vec{\Omega} = \vec{\tau} \times \vec{\kappa} + \gamma \vec{\tau}\) (Darboux vector). Then \[ \frac{1}{|\vec{x}_u|} \vec{e}^j_u = \vec{\Omega} \times \vec{e}^j. \]

• We can do similar for time derivatives. Let \(\vec{\omega} = \vec{\tau} \times \vec{\tau}_t + m \vec{\tau}\) then we have \[ \vec{e}^j_t = \vec{\omega} \times \vec{e}^j. \]

Darboux vector (ii)

• Since the \(u\) and \(t\) time derivatives commute (assuming fixed length), we can deduce that \[ \frac{1}{|\vec{x}_u|} \vec{\omega}_u - \vec{\Omega}_t = \vec{\omega} \times \Omega. \] In particular, \[ \frac{m_u}{| \vec{x}_u |} - \gamma_t + \frac{\vec{x}_{ut}}{| \vec{x}_u |} \cdot ( \vec{\tau} \times \vec{\kappa} ) = 0. \]

Model ingredients

• Passive moment: (linear viscoelastic) \[ \vec{M}^p = \vec{\tau} \times ( A ( \alpha \vec{e}_1 + \beta \vec{e}_2 ) + B ( \alpha_t \vec{e}_1 + \beta_t \vec{e}_2 ) + C \gamma \vec\tau + D \gamma_t \vec\tau ). \]
• Active moment: (from muscles - in application \(\gamma_0 = 0\)) \[ \vec{M}^a = - \vec{\tau} \times ( A \alpha_0 \vec{e}_1 + B \beta_0 \vec{e}_2 ) - D \gamma_0 \vec\tau, \]
• Environment forces: (resistive forces, not fundamental derivation) \[ \vec{D}_{\text{lin}} = \K \vec{x}_t + K_\nu ( {\id}- \tau \otimes \tau ) \vec{x}_t, \qquad \vec{D}_\text{rot} = -K_\text{rot} m \vec{\tau} \]

Summary of equations and unknowns

Equations look like fourth order parabolic equation for \(x\), with local length constraint, coupled to second order parabolic for \(\gamma\).

Discretization ideas

• We split the domain \((0,1)\) into equal elements of length \(h\) and choose a time step \(\Delta t\)

• We use a combination of piecewise linear finite elements \(V_h\) and piecewise constant functions \(Q_h\).

• Use mass lumping for vertex-wise updates.

Two tangent vectors

For a discrete parametrisation \(\vec{x}_h \in V_h^3\), we introduce two different tangent vector fields.

• \(\vec{\tau}_h \in Q_h^3\) as the piecewise constant normalised derivative of \(\vec{x}_h\): \(\vec{\tau}_h = {\vec{x}_{h,u}} / { |\vec{x}_{h,u}| }\).

• \(\tilde{\vec{\tau}}_h \in V_h^3\) with vertex values given by \[ \tilde{\vec{\tau}}_h( u_j, \cdot ) = \frac{ \vec{\tau}_{h}( u_i^-, \cdot) + \vec{\tau}_h( u_i^+, \cdot ) }{| \vec{\tau}_{h}( u_i^-, \cdot) + \vec{\tau}_h( u_i^+, \cdot ) |} \quad \mbox{ for } i = 1, \ldots, N, \]

where \(\vec\tau_h( u_i^\pm, \cdot )\) is \(\vec{\tau}_h\) evaluated on the left (or right) element to the vertex \(u_i\).

Discretization choices

Semi-discrete scheme

Semi-discrete stability

Lemma If \(\alpha^0, \beta^0, \gamma^0\) are independent of time, any solution to the above problem satisfies: \[ \begin{multline*} \int_0^1 ( \K \vec{x}_{h,t} \cdot \vec{x}_{h,t} + K^\rot_h m_h^2 ) | \vec{x}_{h,u} | \\ + \frac{1}{2} \ddt \int_0^1 \bigl( A \bigl( ( \alpha_h - \alpha^0 )^2 + ( \beta_h - \beta^0 )^2 \bigr)_h + C ( \gamma_h - \gamma^0 )^2 \bigr) | \vec{x}_{h,u} | \\ + \int_0^1 \bigl( ( B ( \alpha_{h,t}^2 + \beta_{h,t}^2 ) )_h + D \gamma_{h,t}^2 \bigr) | \vec{x}_{h,u} | = 0. \end{multline*} \]

Full discrete scheme

Full discrete scheme: frame update

Using the abbreviations: \[ \vec{k}_i^n = \tilde{\vec\tau}_h^{n-1}( u_i ) \times \tilde{\vec\tau}_h^n( u_i ), \vec{l}_i^n = \tilde{\vec{\tau}}^n_h( u_i ), \varphi_i^n = \Delta t \, m_h^n( u_i ), \] we apply the Rodrigues formula twice (for \(j=1,2\)): \[ \begin{multline*} \tilde{\vec{e}}_{j,h}^n( u_i ) = \vec{e}_{j,h}^{n-1}( u_i ) ( \tilde{\vec\tau}_h^{n-1}( u_i ) \cdot \tilde{\vec\tau}_h^{n}( u_i )) + \vec{k}_i^n \times \vec{e}_{j,h}^{n-1}( u_i ) \\ + \vec{e}_{j,h}^{n-1}( u_i ) \cdot \vec{k}_i^n \vec{k}_i^n \frac{1}{1 + \tilde{\vec\tau}_h^{n-1}( u_i ) \cdot \tilde{\vec\tau}_h^{n}( u_i )} \end{multline*} \] \[ \begin{multline*} \vec{e}_{j,h}^n( u_i ) = \tilde{\vec{e}}_{j,h}^n( u_i ) \cos( \varphi_i ) + \vec{l}_i^n \times \tilde{\vec{e}}_{j,h}^n( u_i ) \sin( \varphi_i ) \\ \nonumber \qquad + ( \tilde{\vec{e}}_{j,h}^n( u_i ) \cdot \vec{l}_i^n )\vec{l}_i^n ( 1 - \cos( \varphi_i ) ). \end{multline*} \]

Frame update equations maintain frame orthonormality

• The first rotation maps \(\tilde{\vec{\tau}}_h^{n-1}\) onto \(\tilde{\vec{\tau}}_h^n\)
• The second rotates the frame about the new \(\tilde{\vec\tau}_h^n\) (leaving \(\tilde{\vec\tau}_h^n\) unaffected).
• Since we apply the same rotations to the two frame vectors, and these rotations map \(\tilde{\vec\tau}_h^{n-1}\) to \(\tilde{\vec\tau}_h^n\), this update procedure results preserves the orthogonality of the frame vectors at each vertex.
• In practice floating point errors accumulate. The frames could be periodically re-orthonormalized.

Constraint equation gives control of length element

Lemma. If we approximate the constraint equation by \[ \int_0^1 \frac{ \vec{x}_{h,u}^n }{ | \vec{x}_{h,u}^n | } \cdot \vec{x}_{hu}^{n+1} q_h = \int_0^1 L q_h, \] then \[ | \vec{x}_{h,u}^{n+1} | = \frac{ L }{ 1 - \frac{1}{2} | \vec{\tau}_h^{n+1} - \vec{\tau}_h^{n} |^2 } \]

Relaxation

\[ \alpha^0 = 2 \sin( 3 \pi u / 2 ), \quad \beta^0 = 3 \cos( 3 \pi u / 2 ), \quad \gamma^0 = 5 \cos( 2 \pi u ). \]

Relaxation - energy (refinement in space and time)

\[ \mathcal{E}(t^n) := \int_0^1 \Bigl( A | \vec{\kappa}_h - \alpha^0 \vec{e}_{1,h}^n - \beta^0 \vec{e}_{2,h} |^2 + C( \gamma_h^n - \gamma^0 )^2 \Bigr)_h | \vec{x}_{h,u}^n | \]

Relaxation - accuracy tests

• Error in length element: \[ \mathcal{F}_1(t^n) := \left| \int_0^1 | \vec{x}_{h,u}^n | - L \right|. \]

• Error in frame orthonormality: \[ \mathcal{F}_2(t^n) := \left( \sum_{0 \le j_1 \le j_2 \le 2} \int_0^1 | \vec{e}_{j_1,h}^n \cdot \vec{e}_{j_2,h}^n - \delta_{j_1, j_2} |^2 | \vec{x}_{h,u}^n | \right)^{{1}/{2}}. \]

Relaxation - accuracy tests

Planar locomotion

Some 3d results - head twitching

Some 3d results - head twitching

Open questions

• Identify appropriate material parameters!
• Can we show full-discrete stability or error bounds for our numerical scheme (robustness)?
• Can we verify simplifying assumptions (RFT, effects of nonlinear (visco)elasticity, incompressibility)?
• Can we do more to relate this to the biology?