# The 3D coarse-graining formulation of interacting elastohydrodynamic
filaments and multi-body microhydrodynamics

Paul Fuchter  Hermes Bloomfield-Gadêlha<EMAIL_ADDRESS>Department of Engineering Mathematics and Bristol Robotics Laboratory,
University of Bristol, Bristol, UK.

###### Abstract

Elastic filaments are vital to biological, physical and engineering systems,
from cilia driving fluid in the lungs to artificial swimmers and micro-
robotics. Simulating slender structures requires intricate balance of elastic,
body, active, and hydrodynamic moments, all in three-dimensions. Here, we
present a generalised 3D coarse-graining formulation that is efficient,
simple-to-implement, readily extendable and usable for a wide array of
applications. Our method allows for simulation of collections of 3D elastic
filaments, capable of full flexural and torsional deformations, coupled non-
locally via hydrodynamic interactions, and including multi-body
microhydrodynamics of structures with arbitrary geometry. The method exploits
the exponential mapping of quaternions for tracking three-dimensional
rotations of each interacting element in the system, allowing for computation
times up to 150 times faster than a direct quaternion implementation. Spheres
are used as a ‘building block’ of both filaments and solid micro-structures
for straightforward and intuitive construction of arbitrary three-dimensional
geometries present in the environment. We highlight the strengths of the
method in a series of non-trivial applications including bi-flagellated
swimming, sperm-egg scattering, and particle transport by cilia arrays.
Applications to lab-on-a-chip devices, multi-filaments, mono-to-multi
flagellated microorganisms, Brownian polymers, and micro-robotics are
straightforward. A Matlab code is provided for further customization and
generalizations.

##  I Introduction

Elastic passive and active filaments are the building blocks of numerous
biological, physical, engineering and robotic systems. These include, but are
not limited to, polymer and fibre dynamics, cilia driving fluid in the lungs,
microtubules regulating deformations in cancerous cells, flagella propelling
spermatozoa and algae (Fig. 1(a,b)), and robotic micro-swimmers in
microfluidic devices, among many others [1, 2, 3, 4, 5]. Simulating these
structures fully in three-dimensions is challenging due to the convoluted
balance of several interacting components. Moments arise from elasticity,
internal activity, contact and the inertialess fluid that the filaments are
embedded in. This complexity is augmented by the ‘non-local’ presence of other
structures embedded in the fluid, such as walls and solid bodies. Numerous
computational architectures have been developed to date attempting to simulate
these systems [6, 7, 8, 9, 10, 11], a testimony of the growing importance of
filament fluid-structure interactions to the scientific community at large,
and we direct the reader to the recent review on the topic by [1] and
references therein.

A hierarchy of fluid and elastic theories exist from continuous to discrete,
with varying levels of complexity depending upon the precision required.
Methods for resolving inertialess hydrodynamics include local theories, such
as resistive force theory (RFT) [12], and more complex non-local theories,
such as slender-body theories (SBT), regularised Stokeslets [13] and the
Rotne-Prager-Yamakawa (RPY) approximation [14], and boundary-element methods
(BEM) [15]. Similarly, elastic deformation of filaments may be described using
Cosserat rod theory, which accommodates bending, twist, shear and
extensibility in three-dimensions [16], Kirchoff rod approximation [17], where
filaments can only bend and twist, or even constraining inextensible filaments
to planar deformations (with zero twist) [18]. Immersed boundary methods
exploit direct numerical solution of Navier-Stokes equations with a Lagrangian
mesh for the filament [19], whilst in bead models the filament is discretised
into segments, and the elastic coupling is described via discrete ‘elastic
bonds’ [20, 21], among many others.

Figure 1: (a) Uni-flagellated sea urchin sperm [22]. (b) Multiflagellated
green algae Chlamydomonas [5]. (c) Representation of interactions that the
method is capable of simulating: multiple filaments, filaments connected to
arbitrary solid bodies, multi-filament-body structures and micro-structures
with arbitrary shapes, free or fixed.

Classical formulations of solving filament elastohydrodynamics revolves around
the continuum limit of moment balance [23, 24]. They are critical for
analytical progress and direct model interpretability from the resulting
partial differential equations (PDEs) [18]. However, the PDEs have a high-
order, are highly nonlinear, and are coupled with a boundary value problems
(BVP), which require numerous boundary conditions. Unsurprisingly, these
equations are challenging to solve even in lower dimensions [24, 25]. They are
also numerically stiff and thus computationally expensive/time-consuming,
often requiring penalization strategies to regularise numerical errors [23].
They also require specific numerical architecture that is not easily
transferable among different applications, especially in 3D, nor easily
generalizable to systems involving non-trivial geometry and/or possessing
many, and distinct, interacting units, such as elastic filaments interacting
with solid structures (Fig. 1(c)). As we shall see below, the formulation
derived here will attempt to resolve these challenges.

Recent developments in coarse-graining formulations offer intuitive model
interpretability, straightforward implementation, and numerical advantages
over previous methods. By asymptotically integrating the momentum density over
‘coarse’ segments along the filament, the coarse-graining method arises from
first principles, and removes the need to solve for unknown contact forces,
Lagrange multipliers, and associated boundary conditions to enforce filament
inextensibility [11, 26]. The elastohydrodynamic PDE is recast into a simpler
set of ordinary differential equations (ODEs) with trivial and intuitive
matricial form, only depending on the filament parametrization, as we explore
further in this paper. Moreau et al. [11] highlighted the high efficiency of
the formalism in 2D, and showed excellent agreement with classical
formulations, even for relatively large coarse-graining segments. Following
this, the coarse-graining method has gained momentum, and was adopted and
generalised in several ways: (i) non-local hydrodynamics have been accounted
for passive and active filaments in 2D, including wall interaction [27, 28,
29], (ii) 3D filament deformations were accounted for in [30, 31], though
limited to a local hydrodynamic resistive theory (RFT), as in [11], (iii) it
has been applied to study sperm flagellum propulsion in 2D [32], (iv) expanded
into a hybrid immersed boundary [33], as well as, stochastic method [34], and
(v) used for analytical and numerical study of the well-posedness of Newtonian
and viscoelastic elastohydrodynamic propulsion [35, 36].

Despite the above momentum in the literature, there is currently no coarse-
graining formulation that allows simulation of collections of interacting
elastohydrodynamic filaments with arbitrary solid structures in 3D. Here, we
capitalise on the strengths of coarse-graining formulation and develop a fast
simulation tool to resolve multiple interacting elastohydrodynamic filaments
and solid body microhydrodynamics in 3D. Our generalised coarse-graining
formulation is efficient, intuitive and easy-to-implement and customise. Our
method allows for simulation of 3D elastic filaments with both flexural and
torsional deformations, coupled hydrodynamically to free or fixed solid bodies
and microstructures present in the environment (see Fig. 1(c)). We use spheres
as the hydrodynamic ‘building block’ of filaments and solid structures for
straightforward application of the method to novel and arbitrary three-
dimensional geometries, though we note that any hydrodynamic theory may be
employed instead. The method is validated against previous experiments and
simulations in the literature.

We show that despite the substantial advantage of employing the coarse-
graining formulation, runtimes can still be hindered by the choice of
coordinate-system parametrization. We exploit the benefits of the exponential
mapping of quaternions [37, 38] for high-precision tracking of basis-rotations
in three-dimensions. This allows fast and efficient computation, up to 150
times faster than a direct quaternion parametrisation. The exponential map
also avoids undesired repeated change of basis used to circumvent coordinate-
singularities introduced by Euler-angles [17, 30, 31], critical whilst
tracking multiple filaments simultaneously. Finally, we highlight the
strengths of the new method in a series of non-trivial applications and report
novel results. These include the beating of multi-flagellated swimmers, sperm-
egg elastohydrodynamic scattering, and particle transport by cilia arrays.

##  II Methods

We consider collections of 3D inextensible and unshearable elastic filaments
undergoing deformations of bend and twist (Kirchoff rods) that are embedded in
a viscous inertialess fluid, interacting hydrodynamically with each other and
with solid bodies and micro-structures with arbitrary shapes, free or fixed,
present in the micro-environment. We begin by outlining the coarse-graining
formulation solving the dynamics of a single 3D elastohydrodynamic filament
connected to a solid body, before generalising to an arbitrary number of
filaments and solid bodies.

Figure 2: (a) An elastohydrodynamic structure comprised of solid body and
filament. The filament is described by centreline $\mathbf{x}(s)$ with
director basis $[\mathbf{d}_{1}(s),\mathbf{d}_{2}(s),\mathbf{d}_{3}(s)]$ shown
at two arbitrary points $s=s_{1}$ and $s=s_{2}$. (b) The structure is made up
of spheres. Body constructed using $N_{body}=133$ spheres and filament using
$N=13$ rigid segments with $n=3$ spheres per segment. (c) Body position
$\mathbf{x}_{b}$ and body director basis
$[\mathbf{d}_{1}^{b},\mathbf{d}_{2}^{b},\mathbf{d}_{3}^{b}]$ labelled. Segment
position $\mathbf{x}_{i}$ labelled for as well the the director basis for
segments $1$ and $3$.

###  II.1 The coarse-graining formalism in 3D

We consider first a single Kirchoff filament attached to an arbitrary solid
body (Fig. 2(a)). Let $\mathbf{x}(s)$ be the position of the centreline of the
filament, where $s$ is an arclength between $s=0$ and $s=L$, where $L$ is the
length of the filament. To describe deformations of the filament we will use
the orthonormal local director basis
$[\mathbf{d}_{1}(s),\mathbf{d}_{2}(s),\mathbf{d}_{3}(s)]$, defined at every
point along the centreline. As such, $\mathbf{d}_{3}(s)$ is tangent to the
curve $\mathbf{x}(s)$ and given by
$\mathbf{d}_{3}(s)=\partial\mathbf{x}/\partial s$, and $\mathbf{d}_{1}(s)$ and
$\mathbf{d}_{2}(s)$ point normal to the curve $\mathbf{x}(s)$. Derivatives of
the director basis with respect to arclength are related to the twist vector
$\boldsymbol{\kappa}(s)$ by

$\frac{\partial\mathbf{d}_{i}}{\partial
s}=\boldsymbol{\kappa}\times\mathbf{d}_{i},$ (1)

where
$\boldsymbol{\kappa}=\kappa_{1}\mathbf{d}_{1}+\kappa_{2}\mathbf{d}_{2}+\kappa_{3}\mathbf{d}_{3}$.
The constituents of the twist vector $\kappa_{1}$, $\kappa_{2}$ and
$\kappa_{3}$ are known as the filament curvatures [17]. The third curvature in
the $\mathbf{d}_{3}$ direction can also be referred to as the twist density.
We direct the reader to Ref. [39] for a thorough discussion on the interplay
between the three components of the twist vector and the differences between
$\kappa_{3}$, torsion and twist.

The balance of forces and torques along the filament at low Reynolds number
reads [40]

$\displaystyle\mathbf{n}_{s}+\mathbf{f}$ $\displaystyle=0,$ (2)
$\displaystyle\mathbf{m}_{s}+\mathbf{x}_{s}\times\mathbf{n}+\boldsymbol{\tau}$
$\displaystyle=0,$ (3)

where $\mathbf{n}(s)$ and $\mathbf{m}(s)$ are the contact force and contact
moment, $\mathbf{f}(s)$ and $\boldsymbol{\tau}(s)$ are the force and torque
densities exerted on the filament by the fluid, and subscript $s$ denotes
differentiation with respect to $s$. The elastic moment is known
constitutively for a Kirchoff rod, whereas the contact force remains unknown.
The coarse-graining formalism is derived from first principles by exploiting
the integral form of the momentum balance above. This circumvents the need for
solving a high-order PDE with a second-order boundary value problem using
Lagrange multipliers [11]. For this, we discretise the body into $N_{body}$
spheres and the filament into $N$ rigid segments each made of $n$ spheres, as
shown in Fig. 2(b) for $N_{body}=133$, $N=13$ and $n=3$. Positions of spheres
in the body are defined relative to the frame of reference that moves with the
entire structure, $\mathbf{x}_{b}$, described by the local basis
$[\mathbf{d}_{1}^{b},\mathbf{d}_{2}^{b},\mathbf{d}_{3}^{b}]$ (Fig. 2(c)).
Sphere positions are given by $\mathbf{y}_{i}$ for $i=1,2,...,M$, where
$M=N_{body}+Nn$ is the total number of spheres in the system. Body spheres are
labelled by $i\leq N_{body}$ and filament spheres by $i>N_{body}$. Segment
endpoints are located at $\mathbf{x}_{i}$ for $i=1,2,...,N$, and
$\mathbf{x}_{1}$ is defined relative to $\mathbf{x}_{b}$ in the local basis
$[\mathbf{d}_{1}^{b},\mathbf{d}_{2}^{b},\mathbf{d}_{3}^{b}]$.

Invoking force- and torque-free conditions, the total momentum balance of the
elastohydrodynamic filament-body system in Eqs. (2,3) reads

$\displaystyle\sum_{i=1}^{N_{body}+Nn}\mathbf{F}_{i}$ $\displaystyle=0,$ (4)
$\displaystyle\sum_{i=1}^{N_{body}+Nn}\left((\mathbf{y}_{i}-\mathbf{x}_{b})\times\mathbf{F}_{i}+\mathbf{T}_{i}\right)$
$\displaystyle=0,$ (5)

where $\mathbf{F}_{i}$ and $\mathbf{T}_{i}$ are the total force and torque,
respectively, exerted on sphere $i$ by the fluid. Integrating Eq. 3 from
$s=s_{j}$ to $s=L$ yields

$\sum_{i=N_{body}+(j-1)n+1}^{N_{body}+Nn}\left((\mathbf{y}_{i}-\mathbf{x}_{j})\times\mathbf{F}_{i}+\mathbf{T}_{i}\right)=-\mathbf{m}_{j},$
(6)

for $j=2,...,N$, where $\mathbf{m}_{j}$ is the bending moment at $s_{j}$
defined constitutively by

$\mathbf{m}_{j}=E_{b}\kappa_{1}^{j}\mathbf{d}_{1}^{j}+E_{b}\kappa_{2}^{j}\mathbf{d}_{2}^{j}+E_{t}\kappa_{3}^{j}\mathbf{d}_{3}^{j},$
(7)

where $E_{b}$ is the bending stiffness, $E_{t}$ is the torsional stiffness and
$\kappa_{i}^{j}$ are the curvatures at $s_{j}$ for $i=1,2,3$. The coarse-
grained governing equations (4, 5, 6) can be written in a simpler and compact
matrix form, given by

$\bNiceArray{wl{1cm}cccwr{1cm}wl{1cm}cccwr{1cm}}[first-col,code-for-first-
col=\scriptscriptstyle,]\text{Eq.~{}(\ref{total_force})}&\mathbf{I}\Cdots\mathbf{I}\mathbf{0}\Cdots\mathbf{0}\\\
\text{Eq.~{}(\ref{total_torque})}\left[\boldsymbol{\Delta}_{1}\right]_{\times}\Cdots\left[\boldsymbol{\Delta}_{M}\right]_{\times}\mathbf{I}\Cdots\mathbf{I}\\\
\text{Eq.~{}(\ref{sub_torque})
}j=2\mathbf{0}\Cdots\mathbf{0}\Block{3-2}{\mathbf{A}}\mathbf{0}\Cdots\mathbf{0}\Block{3-2}{\mathbf{B}}\\\
\Vdots\Vdots\Ddots\Vdots\Vdots\Ddots\Vdots\\\ \text{Eq.~{}(\ref{sub_torque})
}j=N\mathbf{0}\Cdots\mathbf{0}\mathbf{0}\Cdots\mathbf{0}\\\
\CodeAfter\leavevmode\hbox to11.23pt{\vbox
to7.15pt{\pgfpicture\makeatletter\hbox{\hskip 5.61441pt\lower 5.49055pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}}{{}}{{}}}{}{{}}{}{{{}}{{}}}{}
{}{}{}{}{{{}{}}}{}{{}}{}{}{}{{{}{}}}
{}{{}{}}{}{}{}{{}}{{}}{{}{}}{{}{}}{{{{}{}{{}} }}{{}}}
{}\pgfsys@moveto{0.0pt}{5.69055pt}\pgfsys@stroke\pgfsys@invoke{
}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-3.9479pt}{7.55705pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$3M$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\leavevmode\hbox
to11.23pt{\vbox to7.15pt{\pgfpicture\makeatletter\hbox{\hskip 5.61441pt\lower
5.49055pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{{}}{{}}{{}}}{}{{}}{}{{{}}{{}}}{}
{}{}{}{}{{{}{}}}{}{{}}{}{}{}{{{}{}}}
{}{{}{}}{}{}{}{{}}{{}}{{}{}}{{}{}}{{{{}{}{{}} }}{{}}}
{}\pgfsys@moveto{0.0pt}{5.69055pt}\pgfsys@stroke\pgfsys@invoke{
}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-3.9479pt}{7.55705pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$3M$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\begin{bmatrix}\mathbf{F}_{1}\\\
\vdots\\\ \mathbf{F}_{M}\\\ \mathbf{T}_{1}\\\ \vdots\\\
\mathbf{T}_{M}\end{bmatrix}=-\begin{bmatrix}\mathbf{0}\\\ \mathbf{0}\\\
\mathbf{m}_{2}\\\ \vdots\\\ \mathbf{m}_{N}\end{bmatrix},$ (8)

where we have used the matrix form of the cross product, denoted by
$[\mathbf{a}]_{\times}$ (see Appendix A), and defined
$\boldsymbol{\Delta}_{i}=\mathbf{y}_{i}-\mathbf{x}_{b}$. The blocks of zeros
multiply forces and torques experienced by spheres in the body, as they do not
contribute to the bending moment relations on the filament. The matrices
$\mathbf{A}$ and $\mathbf{B}$ multiply forces and torques experienced by the
spheres in the filament, and are size $3(N-1)\times 3Nn$. The matrix
$\mathbf{A}$ is comprised of $(N-1)Nn$ sub-matrices given by
$A_{i,j}=[\mathbf{y}_{N_{body}+j}-\mathbf{x}_{i+1}]_{\times},\text{ if }j>in$,
and zero otherwise. The matrix $\mathbf{B}$ has the same structure as
$\mathbf{A}$ with sub-matrices given by $B_{i,j}=\mathbf{I},\text{ if }j>in$,
and zero otherwise.

The coarse-grained matrix formulation in Eq. 8 encodes the force and torque
balance for a single filament attached to an arbitrary solid body (Fig. 2(b)).
By setting $N_{body}=0$ or $N=0$, we model a free filament and a free solid
body, respectively, further highlighting the simplicity of the formulation.
The general form of the coarse-grained formulation in Eq. 8 for an arbitrary
number of solid bodies each with an arbitrary number of attached filaments is
simply given by

$\boldsymbol{\mathcal{M}}_{F}\begin{bmatrix}\boldsymbol{\mathcal{F}}\\\
\boldsymbol{\mathcal{T}}\end{bmatrix}=\boldsymbol{\mathcal{K}},$ (9)

where $\boldsymbol{\mathcal{M}}_{F}$ is analogous to the matrix in Eq. 8,
$\boldsymbol{\mathcal{F}}$ and $\boldsymbol{\mathcal{T}}$ contain the forces
and torques, respectively, exerted on all spheres by the fluid.
$\boldsymbol{\mathcal{K}}$ is analogous to the right-hand-side in Eq. 8 and
defines the constitutive relation for each 3D Kirchoff filament in the system.
In the Appendix B, we provide the explicit form of Eq. 9 and include a
practical example of the coarse-grained matrix in the case of two free
filaments. It is worth noting that the coarse-graining formulation in Eq. 9 is
general and does not depend on the specific choice of hydrodynamic coupling.
For this reason, Eq. 9 can be used in conjunction with a variety of
computational and analytic methods for solving the fluid mechanics around
slender bodies. In this paper, we exploit the mathematical capabilities of the
Rotne-Prager-Yamakawa (RPY) non-local hydrodynamic approximation [41] for low
Reynolds number fluids.

###  II.2 Non-local hydrodynamic coupling and dimensionality reduction

We approximate the non-local hydrodynamic forces and torques on each
interacting sphere using the Rotne-Prager-Yamakawa (RPY) tensor [41]. This
relates linear velocities $\mathbf{v}_{i}$ and angular velocities
$\boldsymbol{\omega}_{i}$ with forces and torques acting on each sphere via
the non-local hydrodynamic mobility matrix $\boldsymbol{\mathcal{M}}_{H}$,
giving

$\begin{bmatrix}\boldsymbol{\mathcal{V}}\\\
\boldsymbol{{\Omega}}\end{bmatrix}=\boldsymbol{\mathcal{M}}_{H}\begin{bmatrix}\boldsymbol{\mathcal{F}}\\\
\boldsymbol{\mathcal{T}}\end{bmatrix},$ (10)

where $\boldsymbol{\mathcal{V}}$ and $\boldsymbol{{\Omega}}$ contain the
velocities $\mathbf{v}_{i}$ and angular velocities $\boldsymbol{\omega}_{i}$
of all spheres. The mobility tensor encodes the translational and rotational
movement of each sphere, as well as the non-local coupling between spheres,
depending on the sphere radius $a_{i}$, fluid viscosity $\eta$ and mutual
distance, please see Appendix C for details.

Exploiting the invertibility of the non-local hydrodynamic mobility tensor
above, the generalised coarse-graining formulation (Eq. 9) for the _non-local_
elastohydrodynamic coupling simply reads

$\boldsymbol{\mathcal{M}}_{F}\boldsymbol{\mathcal{M}}_{H}^{-1}\begin{bmatrix}\boldsymbol{\mathcal{V}}\\\
\boldsymbol{{\Omega}}\end{bmatrix}=\boldsymbol{\mathcal{K}}.$ (11)

However, the governing Eq. 11 lists, unnecessarily, all linear
$\mathbf{v}_{i}$ and angular $\boldsymbol{\omega}_{i}$ velocities for each
filament’s sphere as unknowns. Within the same body-filament structure,
translation and rotations of each is sphere are coupled, as they are part of
the same assembly. Hence, linear velocities can be written in terms of the
angular velocities to reduce the dimensionality of a given body-filament
structure (Fig. 2), see Appendix D for details. As a result, $\mathbf{v}_{i}$
and $\boldsymbol{\omega}_{i}$ are a function of the velocity of the centre of
mass of the body $\dot{\mathbf{x}}_{b}$, its angular velocity
$\boldsymbol{\omega}_{b}$, and the angular velocity of each segment
$\boldsymbol{\omega}_{i}^{seg}$, part of the same body-filament structure. By
rigidly attaching the first segment to the body (Fig. 2), we set
$\boldsymbol{\omega}_{1}=\boldsymbol{\omega}_{b}$. Altogether, this reduces
the dimensionality of the system from $6(N_{body}+Nn)$ to $3+3N$ in Eq. (11).
The dimensionality reduction is achieved via

$\begin{bmatrix}\boldsymbol{\mathcal{V}}\\\
\boldsymbol{\Omega}\end{bmatrix}=\boldsymbol{\mathcal{Q}}\boldsymbol{\mathcal{W}},$
(12)

where $\boldsymbol{\mathcal{W}}$ is the reduced velocity system for each
sphere in terms of
${\mathbf{x}}_{b},{\boldsymbol{\omega}_{b}},{\boldsymbol{\omega}_{i}^{seg}}$,
for $i>1$. Substituting Eq. (12) into (11) gives the full non-local, coarse-
graining elastohydodynamic system in 3D

$\boldsymbol{\mathcal{M}}_{F}\boldsymbol{\mathcal{M}}_{H}^{-1}\boldsymbol{\mathcal{Q}}\boldsymbol{\mathcal{W}}=\boldsymbol{\mathcal{K}},$
(13)

where $\boldsymbol{\mathcal{W}}$ contains the reduced unknown velocities of
each interacting unit of the elastohydrodynamic system. Knowledge of
$\boldsymbol{\mathcal{W}}$ provide linear and angular velocities of each
sphere, but without the filament shape information. The filament shape,
position and orientation can be found by numerically integrating these
velocities relative to the laboratory fixed frame of reference, for example,
using the rate of change of the directors base
$\mathrm{d}\mathbf{d}_{i}/{\mathrm{d}t}=\boldsymbol{\omega}\times\mathbf{d}_{i}$.
The latter however introduces severe accuracy issues due to preserving
orthogonality and unit length of the directors [37]. We circumvent these
difficulties in the next section by exploiting the exponential mapping of
quaternions to accurately track rotations of the basis vectors in 3D.

###  II.3 Tracking director rotations in 3D using exponential mapping of
quaternions

Euler angles are commonly used to track 3D rotations of vector basis as
function of three parameters. The challenge of any three-parameter
parametrisation is the existence of coordinate-singularities. With Euler
angles this corresponds to performing two successive rotations about the same
axis, effectively losing a degree of freedom, known as gimbal lock [42].
Coordinate-singularities can be circumvented by continuously changing the
laboratory frame of reference, though this method is inappropriate for
tracking multiple interacting filaments at once as we require here.
Quaternions are canonically used to avoid coordinate-singularities by
introducing four-parameter rotations, although, this can increase considerably
the computational time, as we detail in the results section. Instead, we focus
our study on the exponential mapping of quaternions to track rotations in 3D
(referred here as exponential mapping). The exponential map effectively
reparametrises quaternions into a three-parameter rotation scheme, thus re-
introducing singularities. However, unlike Euler angles, the coordinate-
singularities of the exponential map are trivial to avoid, requiring only a
re-scaling rather than redefining entirely the laboratory frame of reference
[42]. This enables faster and more efficient simulations of multiple
interacting filaments and solid bodies via coarse-graining formulation
introduced above.

Figure 3: (a) Relaxation of filament with $N=40$, $n=2$ and $\mathcal{S}=3$.
Initial condition dark blue and shown at subsequent times $t=0.5,1,1.5,2,2.5$,
ending in a straight configuration. (b) Runtime in hours to solve for
relaxation of single filament using quaternion and exponential map
parametrisation, for various number of segments $N$. Number of spheres per
segment is $n=2$ for all cases. (c) Distance that the end tip of helically
actuated filament makes with the axis of rotation, averaged over time during
steady state. Sperm number
$\mathrm{Sp}^{4}=((4\pi)(\log{L/a}+0.5))\mathcal{S}^{4}$. Results agree well
with simulations from Schoeller et al. [8] and experiments from Coq. et al
[43]. (d) Swimming speed of two swimmer system normalised by swimming speed of
a single swimmer. Inset shows in-phase swimmers at several times over one
period, with the initial distance between swimmers $l$ labelled. Results show
good agreement with data from Delmotte et al. [21].

Next, we will briefly outline quaternions and exponential mapping, and their
coupling with coarse-graining formulation. Quaternions are an extension to the
complex plane into a complex four-dimensional space. A quaternion is given by
$\mathbf{q}=q_{0}+q_{1}i+q_{2}j+q_{2}k=\begin{bmatrix}q_{0}&q_{1}&q_{2}&q_{3}\end{bmatrix}^{\top}$,
where $i^{2}=j^{2}=k^{2}=ijk=-1$. Unit quaternions can be used to rotate a
vector basis $\mathbf{v}$ in space via
$\mathbf{v}\prime=\mathbf{q}\mathbf{v}\mathbf{q}^{-1}$. The axis of rotation
is given by the imaginary part of $\mathbf{q}$ and the angle of rotation by
the real part of $\mathbf{q}$. Specifically,

$\mathbf{q}=\begin{bmatrix}\cos({\theta/2})\\\
\sin{(\theta/2)}\hat{\mathbf{n}}\end{bmatrix}$ (14)

gives the quaternion corresponding to a rotation about the axis
$\hat{\mathbf{n}}$ by $\theta$. The angular velocity of a set of basis vectors
is related to the quaternion by

$\dot{\mathbf{q}}=\frac{1}{2}\begin{bmatrix}-q_{1}&-q_{2}&-q_{3}\\\
q_{0}&q_{3}&-q_{2}\\\ -q_{3}&q_{0}&q_{1}\\\
q_{2}&-q_{1}&q_{0}\end{bmatrix}\boldsymbol{\omega},$ (15)

providing a simple framework for numerical integration of quaternions, though
this increases the dimensionality of the system. Generalising above to account
for the multiple angular velocities of the system $\boldsymbol{\mathcal{W}}$,
the state shape-vector $\boldsymbol{\mathcal{X}}_{quat}$ evolve in time
according to

$\dot{\boldsymbol{\mathcal{X}}}_{quat}=\boldsymbol{\mathcal{C}}\boldsymbol{\mathcal{W}},$
(16)

and is solved together with governing Eq. (13), where the matrix
$\boldsymbol{\mathcal{C}}$ is given in Appendix. E.

The exponential mapping, on the other hand, re-parameterises quaternions back
to a three-parameter parametrization [37]. We define a pure vector known as
the generator $\mathbf{r}$, which is related to a quaternion via the
exponential map,
$\mathbf{q}=\exp({\mathbf{r}})=\begin{bmatrix}\cos{|\mathbf{r}|}&\operatorname{sinc}{|\mathbf{r}|}\mathbf{r}^{\top}\end{bmatrix}^{\top}$,
where $\operatorname{sinc}{x}=\sin{x}/x$. The angular velocity of a set of
basis vectors defined by generator $\mathbf{r}$ can be written conveniently in
matrix form [44],

$\boldsymbol{\omega}=2\mathbf{D}\dot{\mathbf{r}},$ (17)

where the matrix $\mathbf{D}$ is given by

$\begin{aligned} \mathbf{D}=&\left\\{\left[\begin{array}[]{lll}1&0&0\\\
0&1&0\\\
0&0&1\end{array}\right]+\frac{\sin^{2}|\mathbf{r}|}{|\mathbf{r}|^{2}}\left[\begin{array}[]{ccc}0&-r_{3}&r_{2}\\\
r_{3}&0&-r_{1}\\\ -r_{2}&r_{1}&0\end{array}\right]\right.\\\
&\left.-\frac{|\mathbf{r}|-\sin|\mathbf{r}|\cos|\mathbf{r}|}{|\mathbf{r}|^{3}}\left[\begin{array}[]{ccc}r_{2}^{2}+r_{3}^{2}&-r_{1}r_{2}&-r_{1}r_{3}\\\
-r_{1}r_{2}&r_{1}^{2}+r_{3}^{2}&-r_{2}r_{3}\\\
-r_{1}r_{3}&-r_{2}r_{3}&r_{1}^{2}+r_{2}^{2}\end{array}\right]\right\\}\end{aligned}.$
(18)

The determinant $|\mathbf{D}|=\operatorname{sinc}^{2}|\mathbf{r}|$ goes to
zero as $|\mathbf{r}|\to n\pi$, for integer $n\geq 1$. As such, a rescaling is
necessary whenever the magnitude of a generator approaches $n\pi$ to avoid
coordinate-singularity [37]. This can be achieved by rescaling

$\mathbf{r}\to\mathbf{r}-\pi\frac{\mathbf{r}}{|\mathbf{r}|},$ (19)

whenever $|\mathbf{r}|\geq\pi/2$. This rescaling does not change the rotation
performed by the generator $\mathbf{r}$ but keeps it far away from any
coordinate-singularity. Eq. (17) can be generalised to account for the angular
velocities of each sphere in the system $\boldsymbol{\mathcal{W}}$ so that the
state shape-vector $\boldsymbol{\mathcal{X}}_{gen}$ can be evolved in time
using

$\boldsymbol{\mathcal{W}}=\boldsymbol{\mathcal{D}}\dot{\boldsymbol{\mathcal{X}}}_{gen},$
(20)

which can be substituted, directly, into the governing system of equations
(13). This provides a convenient matrix form to evolve in time the non-local,
multi-filament-body elastohystodynamic ODE system

$\boldsymbol{\mathcal{M}}_{F}\boldsymbol{\mathcal{M}}_{H}^{-1}\boldsymbol{\mathcal{Q}}\boldsymbol{\mathcal{D}}\dot{\boldsymbol{\mathcal{X}}}_{gen}=\boldsymbol{\mathcal{K}}.$
(21)

After prescribing an initial state vector and boundary conditions, the
governing system of equations, Eqs. (13,16) for quaternions, and Eq. (21) for
exponential mapping, can be solved directly for the time evolution of the
system’s configuration state. We approximate curvatures in
$\boldsymbol{\mathcal{K}}$ using finite differences of the director basis
(using Eq. 1). We use MATLAB solver ode15s to handle numerical integration,
but note that any ODE solver can be used, highlighting the simplicity of the
matrix representation. A Matlab code is provided to serve as a basis for rapid
customization and further generalizations.

##  III Exponential mapping of quaternions and comparison with PDE
formulations

We proceed with a comparison between the computational efficiency of
quaternions and exponential mapping within the coarse-graining framework. We
also compare our simulations with previous methods with similar levels of
elastohydrodynamic accuracy [8, 21], which have been equally validated against
experiments, though employing distinct fluid-structure interaction approaches.
Finally, in the next section, we provide exemplars of the capability of the
method to a series of complex elastohydrodynamic systems in three-dimensions
that are straightforward to implement using the CG formulation.

We non-dimensionalise the system using length-scale $L$, timescale $T$ given
either by the characteristic frequency $\omega$ of the system or the
relaxation time of the fibre, and force density $E_{b}/L^{3}$. For simplicity,
we set $E_{b}=E_{t}$. Finally, we consider any filament attachment to the body
as a rigid connection, i.e. the filaments are clamped at body. The
dimensionless system of equations is thus multiplied by the non-dimensional
stiffness parameter $\mathcal{S}^{4}=L^{4}\frac{\eta\omega}{E_{b}}$, where
$\eta$ is the fluid viscosity. The stiffness parameter is closely related to
parameters, such as, the elastohydrodynamic number
$\mathrm{E_{h}}=8\pi\mathcal{S}^{4}$, and the sperm number
$\mathrm{Sp}^{4}=(\frac{\xi_{\perp}}{\eta}\mathcal{S}^{4})$, where
$\xi_{\perp}$ denotes the perpendicular resistive force drag coefficient from
RFT used elsewhere [11, 31], and characterises the ratio of viscous and
elastic forces in the system.

We first compare the computational performance between quaternion and
exponential mapping for the simple case of a filament free from external
forces and torques, initially bent and undergoing relaxation dynamics, by
recording the computational time taken to solve for the dynamics. Given no
body exists, $\mathbf{x}_{b}=\mathbf{x}_{1}$, $\mathbf{q}_{b}=\mathbf{q}_{1}$
and $\mathbf{r}_{b}=\mathbf{r}_{1}$. Initial quaternions are given by
$\mathbf{q}_{i}=\begin{bmatrix}\cos({\theta_{i}/2})&\sin{(\theta_{i}/2)}\hat{\mathbf{n}}^{\top}\end{bmatrix}^{\top}$
and generators by $\mathbf{r}_{i}=\theta_{i}/2\hat{\mathbf{n}}$, where
$\theta_{i}=(i-1)\pi/(N-1)$ and
$\hat{\mathbf{n}}=\begin{bmatrix}0&1&0\end{bmatrix}$, for $i=1,...,N$. This
corresponds to a filament bent into a semi-circle, as shown in Fig. 3(a) for
$N=40$, $n=2$. The initial shape unbends towards the unstressed straight
configuration, with the centre of mass position conserved to within $2\times
10^{-2}$ $(L)$. We set $n=2$ and vary the number of segments $N$. Figure. 3(b)
shows the runtime in hours to solve from $t=0$ to $t=20$, for various values
of $N$ for both the quaternion and exponential mapping, with the stiffness
parameter $\mathcal{S}=3$. For $N=10$ segments, runtimes are comparable at 3.8
s and 5.6 s for the exponential map and quaternions, respectively. As $N$
increases the exponential map results in significantly reduced runtimes. For
example, at $N=30$, runtimes are 13.7 s and 3 min for the exponential map and
quaternions, respectively. This difference in runtime increases exponentially,
and for $N=90$ segments the exponential map parametrisation runs over 150
times faster. This increase in runtime is true across all ranges of
$\mathcal{S}$. Hereafter, we proceed our numerical study employing exponential
mapping for more efficient computation.

In order to further test that deformations and displacements of the filament
are being simulated accurately, we compare our results with experiments and
simulations carried out by previous works. Schoeller et al. [8] simulate a
single helically actuated filament, where one end of the filament is held a
constant distance $\delta_{0}$ and angle $\alpha_{0}$ from an axis of
rotation. The base of the filament is rotated at an angular velocity $\omega$
about the axis of rotation. We implement this using a filament with $N=20$ and
$n=1$. Initial generators are all set to
$\mathbf{r}_{i}=\alpha_{0}/2\begin{bmatrix}0&1&0\end{bmatrix}^{\top}$. To
enforce the motion of the base of the filament, we replace the force-free
condition with the kinematic constraint
$\dot{\mathbf{x}}_{1}=-\delta_{0}\begin{bmatrix}\sin{t}&\cos{t}&0\end{bmatrix}^{\top}$,
and the torque-free condition with the kinematic constraint,
$\dot{\mathbf{r}}_{1}=\frac{\alpha_{0}}{2}\begin{bmatrix}\cos{t}&-\sin{t}&0\end{bmatrix}^{\top}.$
We solve the system between $t=0$ and $t=10000$, ensuring the system reaches a
steady state. Once a steady periodic state is achieved, we measure the
distance that the tip of the filament makes with the axis of rotation,
averaged over one period. Figure. 3(c) shows the tip distance compared with
results from Ref. [8] for $\delta_{0}=0.1(L)$ and $\delta_{0}=0.02(L)$ with
$\alpha_{0}=0.2618$. Results agree excellently with both Ref. [8] and
experimental data from Ref. [43]. It is worth noting that Ref. [8] has also
been validated and contrasted with experiments independently. The agreement
with canonical elastohydrodynamic methods in 3D, such as in [8], using the
moment balance PDEs directly, together with Lagrange multipliers to enforce
the inextensibility constraint (not required here), indicates equivalence
between canonical PDE methods and the simpler coarse-grained method employed
here in 3D, as also demonstrated in 2D by [11].

We conduct one further test that places significance on the non-local
hydrodynamics between spheres in the system. Following Delmotte et al. [21],
we simulate the case of a single swimmer and then two identical swimmers
coupled via non-local hydrodynamics. We prescribe a time dependent intrinsic
curvature to a filament made of $N=20$ segments with $n=1$ spheres per segment
so that our hydrodynamic description matches the beads model employed in Ref.
[21]. We use the curvature $\boldsymbol{\kappa}^{0}(s,t)=-K_{0}\sin{(ks-\omega
t)}\mathbf{d}_{1}(s,t)$, with $K_{0}=8.25$, if $s\leq 0.5L$, and
$K_{0}=16.5(L-s)/L$, if $s>0.5L$, and a fixed stiffness parameter
$\mathcal{S}=22.6^{1/4}$. Numerical parameter values are chosen to match Ref.
[21], which itself validated this swimming case with experiments [45] and
simulations [46]. For the case of a single swimmer, the curvature wave causes
planar beating that propels the filament forwards, characterized by a distance
travelled per stroke of $V=0.0671(L/T)$. This swimming speed compares well
with $V=0.066(L/T)$ from simulations in Ref. [21] and $V=0.07(L/T)$ from
experiments in Ref. [45].

We also measure the swimming speed of two active filaments swimming together,
Fig 3(d), both with the same time-dependent curvature, as the single swimmer
above, but set to either in-phase or anti-phase (inset of Fig 3(d) shows in-
phase swimmers). The filaments are set an initial distance $l$ apart and their
initial configuration is taken from the steady-state of the single swimmer
case. The normalised swimming speed of the two-swimmer system is given in Fig.
3(d), and shows a reduction in swimming speed for in-phase swimmers as they
approach each other, and an increase in swimming speed for anti-phase swimmers
as they approach each other. Figure. 3(d) shows excellent agreement with Ref.
[21] for in-phase and anti-phase swimmers. This is despite major differences
in both modelling and computational formulations. In particular, Ref. [21]
employs a novel gears model which solves the moment balance PDE with contact-
contact forces between interlocked spheres making up the elastic filament
which rotate relative to each other similarly to a gear, in such a way that
the unknown contact forces play the role of Lagrange multipliers in the
system. For anti-phase swimmers at small $l$, we used a modified version of
the RPY tensors allowing for overlapping spheres given in Ref. [41], as the
beats of the two swimmers begin to overlap slightly [21].

##  IV 3D non-local Elastohydrodynamic applications

Figure 4: (a) Model Chlamydomonas setup. Angle $2\theta$ between two flagella
attached to body of radius $R$. Here, $2\theta=23$ degrees and $R=0.35(L)$.
(b) Periodic bending wave given by time-dependent intrinsic curvature
$\boldsymbol{\kappa}=\pm(4+4\sin(2\pi s-t)\mathbf{d}_{2}$ shown over one
period $T$. Induces straight swimming of Chlamydomonas indicated by arrow. The
$\mathbf{d}_{2}$ director basis along the entire length of both filaments
points directly into the page, meaning the bending wave produced by this time-
dependent curvature is perfectly in the plane of the page. (c) Swimming speed
of Chlamydomonas dependent upon body radius and non-local hydrodynamic
coupling. Reduced coupling indicates that terms in the RPY tensors coupling
spheres in separate parts of the Chlamydomonas are set to zero (hydrodynamic
coupling between flagellum-flagellum and flagellum-body are turned off). (d)
Adding a constant out of plane curvature of magnitude $\kappa_{1}$ in the
$\mathbf{d}_{1}$ direction causes out of plane swimming. Trajectory of
Chlamydomonas when both flagella have equal out of plane intrinsic curvature
$\kappa_{1}=0.5$. (e) Adding asymmetry between flagella turns the circular
trajectory into a helical trajectory. Here, the flagella have intrinsic out of
plane curvature $\kappa_{1}=0.5\pm 0.3$. (f) Trajectory of Chlamydomonas when
both flagella have equal and opposite out of plane intrinsic curvature
$\kappa_{1}=\pm 0.5$. (g) Asymmetric opposite sign out of plane curvatures
$\kappa_{1}=\pm 0.5+0.3$ turns straight swimming with rolling into a helical
trajectory, with rolling still present.

###  IV.1 Bi-flagellate Chlamydomonas elastohydrodynamic-swimming in 3D

Chlamydomonas is a widespread model organism whose flagella are structured
similarly to those found in mammals, making it an important microorganism in
experiments and modelling of flagella motion [47]. However, modelling
Chlamydomonas is not an easy task due to the multiple elastic and solid
interacting parts. It is composed of a solid body and two elastic tail-like
appendages coupled by contact-contact and hydrodynamic interactions, whilst
the whole system is free from external forces and torques during free-swimming
motion. As such Chlamydomonas modelling is often abstracted with minimal
models [48, 49, 50]. Fauci et al. [51, 52] modelled Chlamydomonas using the
immersed boundary method, but since these studies we are unaware of any
simulation work of Chlamydomonas that includes full elastohydrodynamic
interactions. Here, we demonstrate that simulating this model organism is
straightforward with the coarse-graining formalism presented here.

Figure 5: (a) Waveform of model sperm shown over one beat period T. (b)
Trajectory of sperm past egg with different initial condition for $\theta$,
the angle between the flagellum and $\hat{{z}}$. (c) Swimming speed dependent
on time. Sperm slows down as it gets closer to the egg, and then begins to
speed up as the sperm swims away from the egg. (d) Swimming speed dependent on
sperm-egg separation, shown from $t=0(T)$ to $t=1200(T)$. Arrows indicate the
flow of time. Depending on $\theta$, the sperm reaches a minimum separation
from the egg before swimming away and increasing in speed.

The coarse-grained Chlamydomonas swimmer is made of two active elastic
filaments with an angle $2\theta$ between them, attached to single spherical
body of radius $R$, as depicted in Fig. 4(a). The body is made out of
$N_{body}=184$ spheres and the filaments contain $N=15$ segments with $n=1$
sphere per segment. With $\mathcal{S}=3$, we induce straight swimming by
prescribing an intrinsic time-dependent curvature $\boldsymbol{\kappa}=\pm
4[1+\sin(2\pi s-t)]\mathbf{d}_{2}(s)$ in each filament. The constant term in
the prescribed curvature induces a non-zero average curvature along the
flagella (see Fig. 4(b)), allowing the flagella to progressively ‘pull’ the
spherical body over time, known as a puller swimmer. We vary the body radius
(keeping $N_{body}$ = 184) and record the swimming speed of the Chlamydomonas
swimmer. Fig. 4(c) shows the swmimming speed of each Chlamydomonas as body
size is increased with full non-local hydrodynamic coupling and reduced
coupling (where hydrodynamic coupling between flagellum-flagellum and
flagellum-body are turned off). Increasing the body size linearly decreases
the swimming speed of Chlamydomonas, as expected from the linearity of Stokes
flow. Reduced hydrodynamic coupling results in a slower swimming speed (up to
60% slower for the smallest body radius used here), highlighting the
importance of non-local coupling for Chlamydomonas swimming (not possible with
minimal models), and how this microorganism may exploit this effect to achieve
higher swimming speeds. The anti-phase Chlamydomonas beating is thus critical
to unlock higher propulsion, similarly to hydrodynamic coupled single-swimmers
in Fig. 3(d), that achieve higher speeds when moving in anti-phase near to
each other. Interestingly, the difference in swimming speeds between the full
and reduced hydrodynamic coupling decrease with the body size in Figure. 4(c),
indicating that Stokes law for the spherical body could provide a good
approximation when the body is large (or the flagella is small). This because
the non-local hydrodynamic effect becomes increasingly small, when compared
with the zeroth order contribution from Stokes law, as the body size
decreases.

Motivated by recent empirical studies observing three-dimensional beating of
Chlamydomonas flagella [50, 53], we show here that it is straightforward to
implement 3D beating of Chlamydomonas using the CG framework. For the
simplicity of this demonstration, we introduce a constant out-of-plane
curvature $\kappa_{1}$ in the $\mathbf{d}_{1}$ direction of the two flagella,
so that the intrinsic time-dependent curvature is
$\boldsymbol{\kappa}=\kappa_{1}\mathbf{d}_{1}(s)\pm 4[1+\sin(2\pi
s-t)]\mathbf{d}_{2}(s)$. In the symmetric case, when $\kappa_{1}$ is the same
for the two flagella, the out-of-plane beating causes the swimming path of
Chlamydomonas to be no longer straight. Instead, the Chlamydomonas is trapped
into perfect circular paths with radius $r$, linearly proportional to the mean
out-of-plane curvature $\kappa_{1}$, shown in Fig. 4(d). Interestingly, any
incongruousness in the mean intrinsic curvature between the two flagella is
sufficient to allow Chlamydomonas to swim progressively via an array of
helical paths in 3D, depending on differences in the orientation of their
beating-planes, as observed experimentally [50, 53].

By adding an asymmetry between the two flagella’s out of plane curvatures, the
circular paths seen in Fig. 4(d) turn into helical paths shown in Fig. 4(e).
If instead the two out of plane curvatures are equal and opposite, the
Chlamydomonas swims along a straight trajectory and rolls about its own axis,
shown in Fig. 4(f). Again, breaking the symmetry between the two flagella
produces more complex swimming. In the case of opposite sign and asymmetric
out of plane curvatures, the trajectory of Chlamydomonas becomes helical
(distinct from the previous helical case) (Fig. 4(g)), in further agreement
with 3D observations [50, 53].

###  IV.2 Sperm-egg elastohydrodynamic scattering

The journey of sperm cells to the egg is one of the most important biological
processes found in nature. The fertilization process involves the interaction
of multiple bodies with varying size, active and passive components, and solid
and elastic structures, providing endless opportunities biolgoical soft-matter
and fluid-structure research [2]. Most of sperm swimming modelling focuses on
hydrodynamic simulations using prescribed swimming beat patterns (kinematic
constrained) to resolve swimming trajectories [54]. However, this bypasses the
potential for elastohydrodynamic modulations in the system. In this example,
we model for the first time the elastohydrodynamic scattering between a sperm
and and egg freely floating in the fluid. The CG formalism allows
investigation of elastohydrodynamic swimming modulation and simple
construction of elastic and passive solid bodies, that are fixed or free to
move in the fluid. The latter is particularly challenging given that free
force/torque balance governs the movement of the body as a whole. Here, we
briefly test the hypothesis that spherical egg alone is able to
hydrodynamically attract the sperm to the egg via non-local hydrodynamic
interactions, and flagella bending re-orientation, as have been observed for
bacteria swimming near spherical obstacles [55, 56].

For simplicity, we model the sperm head with a spherical shape of radius
$0.2L$ and the egg as a sphere of radius $1L$. We note that eggs vary in size
depending on the species [57]. The sperm flagella is made of $N=16$ segments
with $n=1$ sphere per segment. Using an initial sperm-egg separation (measured
from centre of sperm head to centre of egg) of $5L$, we vary $\theta$, the
angle that the sperm cell makes with the $z$-axis at $t=0$. Rather than
prescribing a time-dependent intrinsic curvature, we set an active internal
moment density along the flagella, $\mathbf{m}=12k\cos{(ks-t)}\mathbf{d}_{2}$,
that is integrated from $s_{j}$ to $L$ and added to bending moment at $s_{j}$.
This causes the waveform shown in Fig. 5(a). Simulations are run with
$\mathcal{S}=6$. Fig. 5(b) shows the initial condition for a single value of
$\theta$, and the subsequent trajectories of sperm cells towards the egg for
four values of $\theta$. The sperm cell that approaches closest to the egg
shows the most noticeable curve in trajectory, as well as the largest change
in speed as it swims past the egg (Fig. 5(c)).

Fig. 5(d) shows the sperm swimming speed as a function of sperm-egg
separation; initially separated by 5 $(L)$, the sperms approach the egg and
their swimming speed slows. This reduction in swimming speed is up to 28 % for
the sperm with the closest approach. Once the sperm reaches the closest
separation from the egg, its swimming speed abruptly increases before
gradually returning to its free-space swimming speed (i.e. the swimming speed
in absence of solid bodies), as the sperm-egg separation increases. These
simulations indicates that hydrodynamic effects alone are unable to assist
sperm to find eggs with similar sizes as the sperm, even when swimming in
close proximity. This may explain need for other form of guidance mechanisms
through evolution, such as chemotaxis in freshwater fish species [58]. Most
interestingly, Fig. 5(d) seems to indicate the possibility of ‘accelerating’
microorganisms by solid body scattering, though this depends on the proximity
of the swimmer to the free body. The curves in Fig. 5(d) are skewed, and above
the average speed for certain sperm-egg separations are permitted, so that the
accumulated distance travelled is higher for those cases (see purple and
yellow curves in Fig. 5(d)). Trajectories are qualitatively similar to
trajectories presented in Ref. [54]. In all, the exemplary results in Fig. 5
motivate further studies on the design and optimization of micro-environments
to promote, or suppress, cell transport and progression using microfluidic
systems.

Figure 6: (a) Initial condition for $5\times 5$ cilia array with wall made of
spheres and tracer particles positioned above. Tracer particle radius $r=0.01\
(L)$ (plotted larger for visual purposes). (b) Beating of a single cilium over
one period with time progressing from black to white. Cilia beat in the $y-z$
plane and induce net fluid flow in the positive $y$ direction. The phase of
cilia beating varies by $2\pi$ across the array, with cilia at the same $y$
having identical phase. Tracer particles pushed by fluid flow created by cilia
array. (a) Birds-eye view of trajectory of tracer particles. Their direction
of travel is from negative $y$ to positive $y$. (b) Side-view of trajectory of
tracer particles. (c) PDF of tracer particle velocity in the $y$ direction,
dependent on initial height. Higher particles are further from cilia array and
as a result have lower average velocity in the $y$ direction (inset).

###  IV.3 Elastohydrodynamic cilia array and particle transport

Arrays of cilia are found in many places in our body and are responsible a
wide array of biological processes, from driving fluid flow during embryonic
growth to clearing mucus in the lungs [4]. Elastohydrodynamic simulations of
cilia array have been developed for decades now to investigate several multi-
scale aspects of this system, including studies on synchronization phenomena
and metachronal waves [59, 60, 61, 62, 63]. Here, we consider this canonical
system using the CG formalism for the first time. The method outlined here is
capable of simulating arrays of cilia and tracking fluid flow, which we
demonstrate with a simple example. We arrange $25$ cilia into a $5\times 5$
grid arranged on top of a wall made of spheres, as shown in Fig. 6(a). Each
cilia is made of $N=7$ segments with $n=1$ spheres per segment and governed by
a prescribed curvature $\boldsymbol{\kappa}=(1+5\sin(1.5\pi
s-t+\phi))\mathbf{d}_{1}$, Fig. 6(b). Each row (cilia with same location in
$y$) in the cilia array has a different value for the phase $\phi$ such that a
metachronal wave is travelling in the positive $y$ direction, Fig. 6(a).
Specifically, $\phi=(i-1)/5\times 2\pi$ where $i$ corresponds to the row
number of the cilium, for $i=1,...,5$. The beating of each individual cilium
is shown in Fig. 6(b) over one period of oscillation. We note that one may use
modified mobility tensors accounting for the hydrodynamics of a wall, but we
choose to use spheres to model the wall for generality. In our example, the
wall is fixed, but is easily made free, and can also be transformed into a
non-planar ciliated surface (e.g. a curved wall) to investigate the effect of
different wall geometries [61], something that is not immediately possible
using modified hydrodynamic tensors.

The metachronal wave induces a net fluid flow above the cilia array which
transports tracer particles positioned above the array. The resultant
trajectory of tracer particles with an initial height of 1.4 $(L)$ can be seen
in Fig. 6(c,d) from a birds-eye and side-view, respectively. The inset in Fig.
6(c) shows the oscillations of a single tracer particle over the beating
period, and correlates with ciliary beating driving them. Over long time, the
tracer particles are successfully transported to the positive side of the $y$
axis, in accordance with the biased waveform in Fig. 6(b), with particles near
the centre of the array moving faster than those at the sides due to the
finite size of the array. Other boundary effects of the patch of cilia can be
seen in Fig. 6(c), in which particles tend to move upward (downward) for
negative (positive) $y$. The PDF of particle velocity in the $y$ direction for
different initial heights in Fig. 6(e) shows that the particle transport is
not effective for this ciliary beating. The transport velocity is only weakly
shifted to positive values due to small component of the beating in the
direction of transport, with the average velocity of transport decaying with
particle height, shown in Fig. 6(f). Distributions of accessible velocities by
the tracer particles vary noticeably with particle height due to the fast
spacial characteristic decay of the velocity field in low Reynolds number.

##  V Discussion and conclusions

We offer a simple and intuitive method for simulating fully three-dimensional
interacting elastic filaments coupled via non-local multi-body
microhydrodynamics (Matlab code provided). By using asymptotic integration of
momentum balance along each coarse-grained segment [11], we avoid challenging
numerical treatment of high-order fluid-structure interaction PDEs and the
numerically stiff Lagrange multipliers required to enforce the inextensibility
constraint. Our method allows straightforward construction of complicated
multi-filament-body systems involving moving, fixed and coupled-filament
conditions, allowing easier treatment of cumbersome force-free and torque-free
constraints. Our CG formulation features the construction of complex solid
body geometries using spheres as building blocks, and the exponential mapping
of quaternions to overcome coordinate-singularities in 3D with a trivial
rescaling.

The CG framework recasts the 3D non-local elastohydrodynamic PDEs into an ODE
system of equations, bypassing altogether the need of numerically meshing
partial derivatives. This allows the use of generic numerical solvers. As
such, no experience is needed in numerical methods, and only basic knowledge
of linear algebra and matrix manipulation are required to use this framework.
The full elastohydrodynamic system is simply condensed into block-operators,

$\boldsymbol{\mathcal{M}}_{F}\boldsymbol{\mathcal{M}}_{H}^{-1}\boldsymbol{\mathcal{Q}}\boldsymbol{\mathcal{D}}\dot{\boldsymbol{\mathcal{X}}}_{gen}=\boldsymbol{\mathcal{K}},$
(22)

each encoding distinct model interactions of the shape field
$\dot{\boldsymbol{\mathcal{X}}}_{gen}$, namely: momentum balance
$\boldsymbol{\mathcal{M}}_{F}$, hydrodynamic coupling
$\boldsymbol{\mathcal{M}}_{H}^{-1}$, dimensional reduction
$\boldsymbol{\mathcal{Q}}$, geometry of deformation and basis tracking
$\boldsymbol{\mathcal{D}}$, and constitutive relations, internal activity and
other boundary constraints $\boldsymbol{\mathcal{K}}$. Eq. (22) is succinct
and facilitates model explainability, generalisations and/or simplifications,
as well as analytical and numerical analysis, better treatment of boundary
conditions, and easier implementation in different platforms. The formalism
allows model customization without the need of algorithmic or numerical re-
design. For example, the reduction in dimensionality operator may be removed
if needed, the non-local hydrodynamic block may be replaced by other local or
non-local drag theories, or by numerical solutions from a Navier-Stokes solver
(e.g. immersed boundary methods). The exponential mapping of quaternions may
be replaced by other coordinate parametrizations of choice and different
material properties of the filament may be equally invoked. Other types of
forces arising from, for example, gravity, magnetism, electrostatic and
stochastic interactions, are straightforward to implement, offering further
flexibility whilst employing this methodology.

The method has been validated and contrasted against previous experimental and
computational studies showing excellent agreement with the literature (Fig.
3). We show that the exponential mapping parametrization provides a simpler
implementation, and generally faster running times, than a direct quaternion
implementation, but both parametrizations are robust while tracking multiple
interacting filaments. Finally, we showcase the method in three complex
elastohydrodynamic examples: (1) bi-flagellated Chlamydomonas swimming in 3D
(Fig. 4), (2) sperm-egg elastohydrodynamic coupling (Fig. 5), and (3)
canonical particle transport by cilia array (Fig. 6). The elastohydrodynamic
models and results presented for examples (1) and (2) are novel. Indeed, the
computational realization of these systems using the classical PDE formulation
is still challenging today, specifically regarding the implementation of the
boundary conditions and other filament constraints.

Our simulations revealed that the non-local hydrodynamic coupling increases
the elastohydrodynamic propulsion of the Chlamydomonas swimming. Most
importantly, the 3D elastohydrodynamic model is able to predict a bewildering
array of complex 3D swimming trajectories that may arise via simple symmetry-
breaking features of the flagella beat in 3D, from circular to helical
trajectories (Fig. 4), in agreement with [50, 53]. This offers numerous
modelling opportunities motivated by recent 3D observations of this important
model microorganism [50, 53]. We also report novel sperm-scattering results in
which sperm swimming can either be enhanced or reduced by the presence of a
free-body (Fig. 5). Finally, we observe canonical small-scale oscillations of
tracer particles in synchrony with ciliary beating, and show that the non-
local particle transport decays with height (Fig. 6), thus mimicking generic
properties of ciliary systems, for the first time using the CG framework.

Numerous open questions exist in elasto-microhydrodynamics systems in which
this method could be readily employed: self-organisation of flagellar beat
[64, 65], elastohydrodynamic modulation of active bundles [66], non-trivial
body geometries [67], physicochemical active guidance [58], ciliated
microorganisms [68], cilia carpets [69], artificial swimmers [70, 71], soft-
robotics [72, 73], locomotive robotics [74] and the role of torsional
instabilities in filament systems [75], to name a few. Extending the framework
to include inertial forces or the full Cosserat rod Theory would widen the
scope for the method’s use. In all, we hope that the ease-of-implementation
advantages of the method will facilitate the study of non-trivial fluid-
structure interaction systems, and assist researchers to fast-track
implementation and accelerate the modelling cycle, serving communities within,
and away from, mathematical and physical sciences that are equally interested
in simulating these systems.

##

Data accessibility. The MATLAB codes used to run the simulations presented in
this paper are available on github via the link:
Funding. Funded by DTP-EPSRC.
Competing interests. We declare we have no competing interests.
Acknowledgements. This work was carried out using the computational facilities
of the Advanced Computing Research Centre, University of Bristol -
http://www.bristol.ac.uk/acrc/.

##  Appendix A Matrix cross product

The cross product $\mathbf{a}\times\mathbf{b}$ in matrix form is given by

$\mathbf{a}\times\mathbf{b}=[\mathbf{a}]_{\times}\mathbf{b}=\begin{bmatrix}0&-a_{3}&a_{2}\\\
a_{3}&0&-a_{1}\\\ -a_{2}&a_{1}&0\end{bmatrix}\begin{bmatrix}b_{1}\\\ b_{2}\\\
b_{3}\end{bmatrix}.$ (23)

##  Appendix B Force and torque balance of two free filaments

For a single filament with no body, the block of zeros in the large matrix in
Eq. 8 reduces to zero size, and the number of spheres in the system is given
by $M=Nn$. If $\mathbf{M}_{F}$ encodes the force and torque balance of a
single filament, e.g.

$\mathbf{M}_{F}\begin{bmatrix}\boldsymbol{\mathcal{F}}\\\
\boldsymbol{\mathcal{T}}\end{bmatrix}=\boldsymbol{\mathcal{K}},$ (24)

then the force and torque balance of two free filaments is given by

$\begin{bmatrix}(\mathbf{M}_{F})_{1}&\mathbf{0}\\\
\mathbf{0}&(\mathbf{M}_{F})_{2}\end{bmatrix}\begin{bmatrix}\boldsymbol{\mathcal{F}}_{1}\\\
\boldsymbol{\mathcal{T}}_{1}\\\ \boldsymbol{\mathcal{F}}_{2}\\\
\boldsymbol{\mathcal{T}}_{2}\end{bmatrix}=\begin{bmatrix}\boldsymbol{\mathcal{K}}_{1}\\\
\boldsymbol{\mathcal{K}}_{2}\end{bmatrix},$ (25)

where $(\mathbf{M}_{F})_{i}$, $\boldsymbol{\mathcal{F}}_{i}$,
$\boldsymbol{\mathcal{T}}_{i}$ and $\boldsymbol{\mathcal{K}}_{i}$ have
identical structure to the single filament case and $i=1,2$ denotes the
filament number. When writing in general form

$\boldsymbol{\mathcal{M}}_{F}\begin{bmatrix}\boldsymbol{\mathcal{F}}\\\
\boldsymbol{\mathcal{T}}\end{bmatrix}=\boldsymbol{\mathcal{K}},$ (26)

rows and columns in Eq. 25 are rearranged accordingly. Adding additional
structures results adding extra matrices along the diagonal in Eq. 25. For
example, adding a single free sphere to the system would add 6 additional rows
to $\boldsymbol{\mathcal{M}}_{F}$, corresponding to one identity matrix
multiplying the force experienced by the sphere, and one identity matrix
multiplying the torque experienced by the sphere. Six extra zeros would be
added to the right-hand-side corresponding to the force- and torque-free
conditions.

##  Appendix C Hydrodynamic mobility matrix

Using the Rotne-Prager-Yamakawa approximation gives the hydrodynamic mobility
matrix

$\boldsymbol{\mathcal{M}}_{H}=\begin{bmatrix}\boldsymbol{\mu}^{tt}&\boldsymbol{\mu}^{tr}\\\
\boldsymbol{\mu}^{rt}&\boldsymbol{\mu}^{rr}\end{bmatrix},$ (27)

where $t$ and $r$ refer to translational and rotational components,
respectively. The submatrices have general form

$\displaystyle\boldsymbol{\mu}=\begin{bmatrix}\mu_{11}&\mu_{12}&\ldots&\mu_{1N_{sph}}\\\
\mu_{21}&\mu_{22}&\ldots&\mu_{2N_{sph}}\\\ \vdots&\vdots&\ddots&\vdots\\\
\mu_{N_{sph}1}&\mu_{N_{sph}2}&\ldots&\mu_{N_{sph}N_{sph}}\end{bmatrix}.$ (28)

Components of
$\boldsymbol{\mu}^{tt},\boldsymbol{\mu}^{tr},\boldsymbol{\mu}^{rt}$ and
$\boldsymbol{\mu}^{rr}$ are given by [41]

$\displaystyle\mu_{ij}^{tt}=\begin{cases}\frac{1}{6\pi\eta
a_{i}}\mathbf{I}&i=j,\\\ \frac{1}{8\pi\eta
r_{ij}}\left(\left(1+\frac{a_{i}^{2}+a_{j}^{2}}{3r_{ij}^{2}}\right)\mathbf{I}+\right.\\\
\left.\left(1-\frac{a_{i}^{2}+a_{j}^{2}}{r_{ij}^{2}}\right)\hat{\boldsymbol{r}}_{ij}\hat{\mathbf{r}}_{ij}\right)&i\neq
j,\end{cases}$ (29) $\mu_{ij}^{rr}=\begin{cases}\frac{1}{8\pi\eta
a_{i}^{3}}\mathbf{I}&i=j,\\\ \frac{1}{16\pi\eta
r_{ij}^{3}}\left(3\hat{\boldsymbol{r}}_{ij}\hat{\boldsymbol{r}}_{ij}-\mathbf{I}\right)&i\neq
j,\end{cases}$ (30) $\mu_{ij}^{tr}=\mu_{ij}^{rt}=\begin{cases}0&i=j,\\\
\frac{1}{8\pi\eta
r_{ij}^{2}}\boldsymbol{\varepsilon}\cdot\hat{\boldsymbol{r}}_{ij}&i\neq
j,\end{cases}$ (31)

where $a_{i}$ is the radius of sphere $i$, $\eta$ is the fluid viscosity,
$\mathbf{r}_{i,j}$ is the separation of sphere $i$ and $j$ and
$\boldsymbol{\varepsilon}$ is the Levi-Civita symbol. When $i=j$,
$\mu^{tt}_{ij}$ and $\mu^{rr}_{ij}$ reduce to standard force and torque
relation for a sphere in the Stokes regime, and there is no translation-
rotation component. For $i\neq j$, the tensors dictate how sphere $j$ affects
the movement of sphere $i$. For example, $\mu^{tt}_{ij}$ couples the force on
sphere $j$ with the velocity of sphere $i$; $\mu^{tr}_{ij}$ couples the torque
on sphere $j$ with the velocity of sphere $i$ and $\mu^{rr}_{ij}$ couples the
torque on sphere $j$ with the angular velocity of sphere $i$. The tensors used
here allow for spheres of different size, allowing different structures to be
constructed out of differently sized spheres, depending upon the precision
required.

##  Appendix D Reducing dimensionality of system

Consider a single filament attached to a solid body. Spheres are rigidly
attached within each segment, and the velocity at $\mathbf{x}_{i}$ is given by
the velocity at $\mathbf{x}_{i-1}$ plus a contribution from the angular
velocity of segment $i-1$. Additionally, the velocity of spheres in the body
is given by the velocity $\dot{\mathbf{x}}_{b}$ plus a contribution from the
angular velocity of the body $\boldsymbol{\omega}_{b}$ As such, the velocity
of a sphere in the system is given by

$\mathbf{v}_{i}=\begin{cases}\dot{\mathbf{x}}_{b}-\boldsymbol{\Delta}_{i}\times\boldsymbol{\omega}_{b}&i\leq
N_{body},\\\
\dot{\mathbf{x}}_{b}-(\mathbf{x}_{1}-\mathbf{x}_{b})\times\boldsymbol{\omega}_{b}-\Delta
s\sum_{j=1}^{\lceil
k/n\rceil-1}\mathbf{d}_{3}^{j}\times\boldsymbol{\omega}_{j}^{seg}-\mathbf{d}^{\lceil
k/n\rceil}_{3}\times\boldsymbol{\omega}_{\lceil
k/n\rceil}^{seg}&i>N_{body},\end{cases}\\\ $ (32)

where $\Delta s$ is the length of each segment, $k=i-N_{body}$, $\lceil
k/n\rceil$ rounds $k/n$ up to the nearest integer and
$\boldsymbol{\omega}_{i}^{seg}$ is the angular velocity of the
$i\textsuperscript{th}$ segment in the filament. The quantity $\lceil
k/n\rceil$ gives the segment number that sphere $i$ is attached to. The
angular velocity of a sphere in the system is given by

$\boldsymbol{\omega}_{i}=\begin{cases}\boldsymbol{\omega}_{b}&i\leq
N_{body},\\\ \boldsymbol{\omega}_{\lceil
k/n\rceil}^{seg}&i>N_{body}.\end{cases}$ (33)

Equations (32) and (33) allow all $\mathbf{v}_{i}$ and
$\boldsymbol{\omega}_{i}$ in the system to be written in terms of velocity of
the centre of mass of the body $\dot{\mathbf{x}}_{b}$, its angular velocity
$\boldsymbol{\omega}_{b}$ and the angular velocity of each filament segment
$\boldsymbol{\omega}_{i}^{seg}$. By additionally rigidly attaching the first
segment of the filament to the body, the dimensionality of the system is
reduced from $6(N_{body}+Nn)$ to $3+3N$ velocities. Writing equations (32) and
(33) in matrix form gives Eq. (12) in the main text.

##  Appendix E Implementing quaternion parametrisation

Consider again the single filament attached to a body. Assign a quaternion
$\mathbf{q}_{b}$ to the body frame vectors and $\mathbf{q}_{i}$ to filament
segment frame vectors. The state vector for the system is

$\mathbf{X}_{quat}=\begin{bmatrix}\mathbf{x}_{b}^{\top}&\mathbf{q}_{b}^{\top}&\mathbf{q}_{1}^{\top}&\ldots&\mathbf{q}_{N}^{\top}\end{bmatrix}^{\top}.$
(34)

The rate of change of the state vector is given by

$\dot{\mathbf{X}}_{quat}=\bNiceMatrix\mathbf{I}&\mathbf{0}\Cdots\mathbf{0}\\\
\mathbf{0}\mathbf{P}(\mathbf{q}_{b})\mathbf{0}\Cdots\mathbf{0}\\\
\mathbf{0}\mathbf{0}\mathbf{P}(\mathbf{q}_{1})\mathbf{0}\Cdots\mathbf{0}\\\
\Vdots\Vdots\\\ \mathbf{0}\Cdots\mathbf{0}\mathbf{P}(\mathbf{q}_{N})\\\
\CodeAfter\line missing{3.0}{5-6}\bNiceMatrix\dot{\mathbf{x}}_{b}\\\
\boldsymbol{\omega}_{b}\\\ \boldsymbol{\omega}_{1}^{seg}\\\ \vdots\\\
\boldsymbol{\omega}_{N}^{seg}$ (35)

where $\mathbf{P}(\mathbf{q})$ is given by (see Eq. 15)

$\mathbf{P}(\mathbf{q})=\frac{1}{2}\begin{bmatrix}-q_{1}&-q_{2}&-q_{3}\\\
q_{0}&q_{3}&-q_{2}\\\ -q_{3}&q_{0}&q_{1}\\\ q_{2}&-q_{1}&q_{0}\end{bmatrix}.$
(36)

By rigidly attaching the first segment of the filament to the body,
$\boldsymbol{\omega}_{1}^{seg}=\boldsymbol{\omega}_{b}$, the dimension of the
vector on the RHS of Eq. 35 can be reduced by 3, in addition to removing one
column from the matrix pre-multiplying the vector. Generalising to arbitrary
systems gives the dynamics of the state vector
$\boldsymbol{\mathcal{X}}_{quat}$ in terms of the reduced velocity system
$\boldsymbol{\mathcal{W}}$ via the matrix $\boldsymbol{\mathcal{C}}$,

$\dot{\boldsymbol{\mathcal{X}}}_{quat}=\boldsymbol{\mathcal{C}}\boldsymbol{\mathcal{W}}.$
(37)

## References

  * [1] du Roure O, Lindner A, Nazockdast EN, Shelley MJ.  Dynamics of Flexible Fibers in Viscous Flows and Fluids.  Annual Review of Fluid Mechanics. 2019;51(1).
  * [2] Gaffney EA, Gadêlha H, Smith DJ, Blake JR, Kirkman-Brown JC.  Mammalian Sperm Motility: Observation and Theory.  Annual Review of Fluid Mechanics. 2011 Jan;43(1).
  * [3] Lauga E, Powers TR.  The hydrodynamics of swimming microorganisms.  Reports on Progress in Physics. 2009 Sep;72(9).
  * [4] Cartwright JHE, Piro O, Tuval I.  Fluid-dynamical basis of the embryonic development of left-right asymmetry in vertebrates.  Proceedings of the National Academy of Sciences. 2004 May;101(19).
  * [5] Elgeti J, Winkler RG, Gompper G.  Physics of microswimmers—single particle motion and collective behavior: a review.  Reports on Progress in Physics. 2015 Apr;78(5).
  * [6] Ryan PM, Wolgemuth CW.  A finite volume algorithm for the dynamics of filaments, rods, and beams.  Journal of Computational Physics. 2022 Oct;466.
  * [7] Koshakji A, Chomette G, Turner J, Jablonski J, Haynes A, Carlucci D, et al.  A robust computational framework for simulating the dynamics of large assemblies of highly-flexible fibers immersed in viscous flow.  Journal of Computational Physics. 2023 Feb;474.
  * [8] Schoeller SF, Townsend AK, Westwood TA, Keaveny EE.  Methods for suspensions of passive and active filaments.  Journal of Computational Physics. 2021 Jan;424.
  * [9] Maxian O, Mogilner A, Donev A.  Integral-based spectral method for inextensible slender fibers in Stokes flow.  Phys Rev Fluids. 2021 Jan;6:014102.
  * [10] Maxian O, Sprinkle B, Peskin CS, Donev A.  Hydrodynamics of a twisting, bending, inextensible fiber in Stokes flow.  Phys Rev Fluids. 2022 Jul;7:074101.
  * [11] Moreau C, Giraldi L, Gadêlha H.  The asymptotic coarse-graining formulation of slender-rods, bio-filaments and flagella.  Journal of The Royal Society Interface. 2018 Jul;15(144).
  * [12] Gray J, Hancock GJ.  The Propulsion of Sea-Urchin Spermatozoa.  Journal of Experimental Biology. 1955 12;32(4):802-14.
  * [13] Cortez R.  Regularized Stokeslet segments.  Journal of Computational Physics. 2018 Dec;375.
  * [14] Wajnryb E, Mizerski KA, Zuk PJ, Szymczak P.  Generalization of the Rotne–Prager–Yamakawa mobility and shear disturbance tensors.  Journal of Fluid Mechanics. 2013 Sep;731.
  * [15] Pozrikidis C.  Shear flow over cylindrical rods attached to a substrate.  Journal of Fluids and Structures. 2010 Apr;26(3).
  * [16] Cao DQ, Tucker RW.  Nonlinear dynamics of elastic rods using the Cosserat theory: Modelling and simulation.  International Journal of Solids and Structures. 2008 Jan;45(2).
  * [17] Nizette M, Goriely A.  Towards a classification of Euler–Kirchhoff filaments.  J Math Phys. 1999;40(6).
  * [18] Wiggins CH, Goldstein RE.  Flexive and Propulsive Dynamics of Elastica at Low Reynolds Number.  Physical Review Letters. 1998;80(17).
  * [19] Fauci LJ, Peskin CS.  A computational model of aquatic animal locomotion.  Journal of Computational Physics. 1988 Jul;77(1).
  * [20] Elgeti J, Kaupp UB, Gompper G.  Hydrodynamics of Sperm Cells near Surfaces.  Biophysical Journal. 2010 Aug;99(4).
  * [21] Delmotte B, Climent E, Plouraboué F.  A general formulation of Bead Models applied to flexible fibers and active filaments at low Reynolds number.  Journal of Computational Physics. 2015 Apr;286.
  * [22] Gong A, Rode S, Kaupp UB, Gompper G, Elgeti J, Friedrich BM, et al.  The steering gaits of sperm.  Philosophical Transactions of the Royal Society B: Biological Sciences. 2020 Feb;375(1792).
  * [23] Tornberg AK, Shelley MJ.  Simulating the dynamics and interactions of flexible fibers in Stokes flows.  Journal of Computational Physics. 2004 May;196(1).
  * [24] Gadelha H, Gaffney EA, Smith DJ, Kirkman-Brown JC.  Nonlinear instability in flagellar dynamics: a novel modulation mechanism in sperm migration?  Journal of The Royal Society Interface. 2010.
  * [25] Gadêlha H, Gaffney EA.  Flagellar ultrastructure suppresses buckling instabilities and enables mammalian sperm navigation in high-viscosity media.  Journal of The Royal Society Interface. 2019 Mar;16(152).
  * [26] Jabbarzadeh M, Fu HC.  A numerical method for inextensible elastic filaments in viscous fluids.  Journal of Computational Physics. 2020 Oct;418.
  * [27] Hall-McNair AL, Montenegro-Johnson TD, Gadêlha H, Smith DJ, Gallagher MT.  Efficient implementation of elastohydrodynamics via integral operators.  Physical Review Fluids. 2019 Nov;4(11).
  * [28] Walker BJ, Ishimoto K, Gadêlha H, Gaffney EA.  Filament mechanics in a half-space via regularised Stokeslet segments.  Journal of Fluid Mechanics. 2019;879:808–833.
  * [29] Walker BJ, Gaffney EA.  Regularised non-uniform segments and efficient no-slip elastohydrodynamics.  Journal of Fluid Mechanics. 2021;915:A51.
  * [30] El Alaoui-Faris Y, Pomet JB, Régnier S, Giraldi L.  Optimal actuation of flagellar magnetic microswimmers.  Physical Review E. 2020 Apr;101(4).
  * [31] Walker BJ, Ishimoto K, Gaffney EA.  Efficient simulation of filament elastohydrodynamics in three dimensions.  Physical Review Fluids. 2020 Dec;5(12).
  * [32] Neal CV, Hall-McNair AL, Kirkman-Brown J, Smith DJ, Gallagher MT.  Doing more with less: The flagellar end piece enhances the propulsive effectiveness of human spermatozoa.  Physical Review Fluids. 2020 Jul;5(7).
  * [33] Ntetsika M, Papadopoulos P.  Numerical simulation and predictive modeling of an inextensible filament in two-dimensional viscous shear flow using the Immersed Boundary/Coarse-Graining Method and Artificial Neural Networks.  Computer Methods in Applied Mechanics and Engineering. 2022 Nov;401.
  * [34] Ntetsika M, Papadopoulos P.  A Hybrid Immersed Boundary/Coarse-Graining Method for Modeling Inextensible Semi-Flexible Filaments in Thermally Fluctuating Fluids.  CMES - Computer Modelling in Engineering Sciences. 2021;129(3):1243-58.
  * [35] Mori Y, Ohm L.  Well-posedness and applications of classical elastohydrodynamics for a swimming filament.  arXiv preprint arXiv:220803350. 2022.
  * [36] Ohm L.  Well-posedness of a viscoelastic resistive force theory and applications to swimming.  arXiv preprint arXiv:220914160. 2022.
  * [37] Boyle M.  The Integration of Angular Velocity.  Advances in Applied Clifford Algebras. 2017 Sep;27(3).
  * [38] Rucker C.  Integrating Rotations Using Nonunit Quaternions.  IEEE Robotics and Automation Letters. 2018 Oct;3(4).
  * [39] Goriely A.  The Mathematics and Mechanics of Biological Growth. vol. 45 of Interdisciplinary Applied Mathematics.  New York, NY: Springer New York; 2017.
  * [40] Antman S.  Nonlinear Problems of Elasticity; 2nd ed.  Dordrecht: Springer; 2005.
  * [41] Zuk PJ, Wajnryb E, Mizerski KA, Szymczak P.  Rotne–Prager–Yamakawa approximation for different-sized particles in application to macromolecular bead models.  Journal of Fluid Mechanics. 2014 Feb;741.
  * [42] Grassia FS.  Practical Parameterization of Rotations Using the Exponential Map.  Journal of Graphics Tools. 1998 Jan;3(3).
  * [43] Coq N, Roure Od, Fermigier M, Bartolo D.  Helical beating of an actuated elastic filament.  Journal of Physics: Condensed Matter. 2009 May;21(20).
  * [44] Boyle M.  Angular velocity of gravitational radiation from precessing binaries and the corotating frame.  Physical Review D. 2013 May;87(10).
  * [45] Bilbao A, Padmanabhan V, Rumbaugh K, Vanapalli S, Blawzdziewicz J.  Navigation and chemotaxis of nematodes in bulk and confined fluids.  In: APS March Meeting Abstracts. vol. 2013; 2013. p. Y39-002.
  * [46] Majmudar T, Keaveny EE, Zhang J, Shelley MJ.  Experiments and theory of undulatory locomotion in a simple structured medium.  Journal of The Royal Society Interface. 2012 Aug;9(73).
  * [47] Sasso S, Stibor H, Mittag M, Grossman AR.  From molecular manipulation of domesticated Chlamydomonas reinhardtii to survival in nature.  eLife. 2018 Nov;7.
  * [48] Friedrich BM, Jülicher F.  Flagellar Synchronization Independent of Hydrodynamic Interactions.  Physical Review Letters. 2012 Sep;109(13).
  * [49] Bennett RR, Golestanian R.  Emergent Run-and-Tumble Behavior in a Simple Model of Chlamydomonas with Intrinsic Noise.  Physical Review Letters. 2013 Apr;110(14).
  * [50] Cortese D, Wan KY.  Control of Helical Navigation by Three-Dimensional Flagellar Beating.  Physical Review Letters. 2021 Feb;126(8).
  * [51] Fauci LJ.  Computational modeling of the swimming of biflagellated algal cells.  Contemporary Mathematics. 1993;141.
  * [52] Fauci LJ.  A computational model of the fluid dynamics of undulatory and flagellar swimming.  American Zoologist. 1996;36(6).
  * [53] Wilson LG, Bees MA. The chiral beat of algal flagella: force and torque via imaging. arXiv; 2022.  Available from: https://arxiv.org/abs/2210.08669.
  * [54] Ishimoto K, Cosson J, Gaffney EA.  A simulation study of sperm motility hydrodynamics near fish eggs and spheres.  Journal of Theoretical Biology. 2016 Jan;389.
  * [55] Takagi D, Palacci J, B Braunschweig A, J Shelley M, Zhang J.  Hydrodynamic capture of microswimmers into sphere-bound orbits.  Soft Matter. 2014;10(11).
  * [56] Spagnolie SE, Moreno-Flores GR, Bartolo D, Lauga E.  Geometric capture and escape of a microswimmer colliding with an obstacle.  Soft Matter. 2015 Apr;11(17).
  * [57] Ahlstrom EH, Moser HG.  Characters useful in identification of pelagic marine fish eggs.  CalCOFI rep. 1980;21.
  * [58] Kholodnyy V, Gadêlha H, Cosson J, Boryshpolets S.  How do freshwater fish sperm find the egg? The physicochemical factors guiding the gamete encounters of externally fertilizing freshwater fish.  Reviews in Aquaculture. 2020;12(2).
  * [59] Gueron S, Levit-Gurevich K, Liron N, Blum JJ.  Cilia internal mechanism and metachronal coordination as the result of hydrodynamical coupling.  Proceedings of the National Academy of Sciences. 1997 Jun;94(12).
  * [60] Meng F, Bennett RR, Uchida N, Golestanian R.  Conditions for metachronal coordination in arrays of model cilia.  Proceedings of the National Academy of Sciences. 2021 Aug;118(32).
  * [61] Westwood TA, Keaveny EE. Coordinated motion of active filaments on spherical surfaces. arXiv; 2021.  Available from: https://arxiv.org/abs/2109.13578.
  * [62] Elgeti J, Gompper G.  Emergence of metachronal waves in cilia arrays.  Proceedings of the National Academy of Sciences. 2013 Mar;110(12).
  * [63] Gilpin W, Bull MS, Prakash M.  The multiscale physics of cilia and flagella.  Nature Reviews Physics. 2020 Feb;2(2).
  * [64] Oriola D, Gadêlha H, Casademunt J.  Nonlinear amplitude dynamics in flagellar beating.  Royal Society Open Science. 2017;4(3).
  * [65] Chakrabarti B, Saintillan D.  Spontaneous oscillations, beating patterns, and hydrodynamics of active microfilaments.  Physical Review Fluids. 2019 Apr;4(4).
  * [66] Pochitaloff M, Miranda M, Richard M, Chaiyasitdhi A, Takagi Y, Cao W, et al.  Flagella-like beating of actin bundles driven by self-organized myosin waves.  Nature Physics. 2022 Oct;18(10).
  * [67] Wang P, Azad MARA, Yang X, Martelli PR, Cheung KY, Shi J, et al.  Self-adaptive and efficient propulsion of Ray sperms at different viscosities enabled by heterogeneous dual helixes.  Proceedings of the National Academy of Sciences. 2021 Jun;118(23).
  * [68] Drescher K, Leptos KC, Tuval I, Ishikawa T, Pedley TJ, Goldstein RE.  Dancing Volvox: Hydrodynamic Bound States of Swimming Algae.  Physical Review Letters. 2009 Apr;102(16).
  * [69] Wang W, Liu Q, Tanasijevic I, Reynolds MF, Cortese AJ, Miskin MZ, et al.  Cilia metasurfaces for electronically programmable microfluidic manipulation.  Nature. 2022 May;605(7911).
  * [70] Magdanz V, Khalil ISM, Simmchen J, Furtado GP, Mohanty S, Gebauer J, et al.  IRONSperm: Sperm-templated soft magnetic microrobots.  Science Advances. 2020 Jul;6(28).
  * [71] Magdanz V, Sanchez S, Schmidt OG.  Development of a Sperm-Flagella Driven Micro-Bio-Robot.  Advanced Materials. 2013;25(45).
  * [72] Milana E, Zhang R, Vetrano MR, Peerlinck S, De Volder M, Onck PR, et al.  Metachronal patterns in artificial cilia for low Reynolds number fluid propulsion.  Science Advances. 2020 Dec;6(49).
  * [73] Li C, Iscen A, Sai H, Sato K, Sather NA, Chin SM, et al.  Supramolecular–covalent hybrid polymers for light-activated mechanical actuation.  Nature Materials. 2020 Aug;19(8).
  * [74] Renda F, Seneviratne L.  A Geometric and Unified Approach for Modeling Soft-Rigid Multi-Body Systems with Lumped and Distributed Degrees of Freedom.  In: 2018 IEEE International Conference on Robotics and Automation (ICRA); 2018. .
  * [75] Chakrabarti B, Liu Y, LaGrone J, Cortez R, Fauci L, du Roure O, et al.  Flexible filaments buckle into helicoidal shapes in strong compressional flows.  Nature Physics. 2020 Jun;16(6).