C. elegans locomotion in three dimensions

The muscles

doi:10.3908/wormatlas.1.7 | From Omer Yuval

Framed curve: Kirchhoff-Love rod

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\).

Model ingredients

• Total moment \[\vec{M} = \vec{\tau} \times \vec{y} + z \vec{\tau},\]
• Bending moment \[\vec{y} = A ((\alpha-\alpha^0) \vec{e}_1 + (\beta - \beta^0) \vec{e}_2) + B (\alpha_t \vec{e}_1 + \beta_t \vec{e}_2).\]
• Twisting moment \[z = C (\gamma - \gamma^0) + D \gamma_t.\]
• Environment forces: (resistive forces, not fundamental derivation) \[ \vec{D}_{\text{lin}} = K_\tau (\vec\tau \otimes \vec\tau) \vec{x}_t + K_\nu ( {\mathbb{I}}- \vec\tau \otimes \vec\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 equation for anti-derivative of \(\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 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}_{h,u}^{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} |_h^2 + C( \gamma_h^n - \gamma^0 )^2 \Bigr) | \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

Some 3d results - head twitching

Some 3d results - head twitching

Optimisation structure

Given sequence of positions \(\vec{x}_d\) find controls \(\alpha^0, \beta^0, \gamma^0\) and initial frame \(\vec{e}^{1}_0\), \(\vec{e}^2_0\) to minimise \[ \begin{align} & \int_0^T \norm{\vec{x} - \vec{x}_d}^2 \, \mathrm{d}t \\ & + \sum_{c \in \{\alpha^0, \beta^0, \gamma^0\}} \int_0^T (w_c \norm{c}^2 + w_c^t \norm{c_t}^2 + w_c^u \norm{\frac{c_u}{| \vec{x}_u |}}^2 ) \, \mathrm{d}t \\ & + \norm{ \frac{1}{|\vec{x}_u|} (\vec{e}^1_{0,u} + \vec{e}^2_{0,u})}^2 \end{align} \] subject to pde, i.e. \((\vec{e}^1_0, \vec{e}^2_0, \alpha^0, \beta^0, \gamma^0) \mapsto (\vec{x}, \vec{e}^1, \vec{e}^2)\).

Preliminary results (i)

Preliminary results (ii)