The three-dimensional coarse-graining formulation of interacting elastohydrodynamic filaments and multi-body microhydrodynamics

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 generalized three-dimensional (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 3D 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 microstructures 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.


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 ( figure 1a,b), and robotic micro-swimmers in microfluidic devices, among many others [2][3][4][5][6]. Simulating these structures fully in three dimensions (3D) 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 [7][8][9][10][11][12], 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 du-Roure et al. [2] 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) [13], and more complex non-local theories, such as slender-body theories (SBT), regularized Stokeslets [14], the Rotne-Prager-Yamakawa (RPY) approximation [15] and boundary-element methods (BEM) [16]. Similarly, elastic deformation of filaments may be described using Cosserat rod theory, which accommodates bending, twist, shear and extensibility in 3D [17], the Kirchoff rod approximation [18], where filaments can only bend and twist, or even constraining inextensible filaments to planar deformations (with zero twist) [19]. Immersed boundary methods exploit direct numerical solution of Navier-Stokes equations with a Lagrangian mesh for the filament [20], while in bead models the filament is discretized into segments, and the elastic coupling is described via discrete elastic bonds [21,22].
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) [19]. However, the PDEs have a high order, are highly nonlinear and are coupled with 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 regularize numerical errors [23]. They also require specific numerical architecture that is not easily transferable to 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 (figure 1c). As we shall see below, the formulation derived here will attempt to resolve these challenges.
Recent developments in coarse-graining (CG) formulations [12] 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 CG 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 [12,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. [12] highlighted the high efficiency of the formalism in two dimensions (2D), and showed excellent agreement with classical formulations, even for relatively large CG segments. Following this, the CG method has gained momentum, and was adopted and generalized in several ways: (i) non-local hydrodynamics have been accounted for passive and active filaments in 2D, including wall interactions [27][28][29], (ii) 3D filament deformations were accounted for in [30,31], though limited to a local hydrodynamic resistive force theory (RFT), as in [12], (iii) it has been applied to study sperm flagellum propulsion in 2D [32], (iv) expanded into a hybrid immersed boundary method [33], as well as a 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 CG formulation that allows simulation of collections of interacting elastohydrodynamic filaments with arbitrary solid structures in 3D. Here, we capitalize on the strengths of CG formulation and develop a fast simulation tool to resolve multiple interacting elastohydrodynamic filaments and solid body microhydrodynamics in 3D. Our generalized CG formulation is efficient, intuitive and easy to implement and customize. 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 (figure 1c). We use spheres as the hydrodynamic 'building block' of filaments and solid structures for straightforward application of the method to novel and arbitrary 3D 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 CG 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 parametrization. The exponential map also avoids undesired repeated change of basis used to circumvent coordinate-singularities introduced by Euler angles [18,30,31], critical while 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. Representation of interactions that the method is capable of simulating: multiple filaments, filaments connected to arbitrary solid bodies, multi-filament-body structures and microstructures with arbitrary shapes, free or fixed.

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 microstructures with arbitrary shapes, free or fixed, present in the micro-environment. We begin by outlining the CG formulation solving the dynamics of a single 3D elastohydrodynamic filament connected to a solid body, before generalizing to an arbitrary number of filaments and solid bodies.

The coarse-graining formalism in three dimensions
We consider first a single Kirchoff filament attached to an arbitrary solid body (figure 2a). Let 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 [d 1 (s), d 2 (s), d 3 (s)], defined at every point along the centreline. As such, d 3 (s) is tangent to the curve x(s) and given by d 3 (s) = ∂x/∂s, and d 1 (s) and d 2 (s) point normal to the curve x(s). Derivatives of the director basis with respect to arclength are related to the twist vector κ(s) by The constituents of the twist vector κ 1 , κ 2 and κ 3 are known as the filament curvatures [18]. The third curvature in the d 3 direction can also be referred to as the twist density. We direct the reader to Goriely [39] for a thorough discussion on the interplay between the three components of the twist vector and the differences between κ 3 , torsion and twist. The balance of forces and torques along the filament at low Reynolds number reads [40] n s þ f ¼ 0 ð2:2Þ and where n s and m s are the contact force and contact moment, f(s) and τ(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 CG 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 BVP using Lagrange multipliers [12]. For this, we discretize the body into N body spheres and the filament into N rigid segments each made of n spheres, as shown in figure 2b 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, described by the local basis ½d b 1 , d b 2 , d b 3 and located at x b (figure 2c). Sphere positions are given by y k for k = 1, 2, …, M, where M = N body + Nn is the total number of spheres in the system. Body spheres are labelled by k ≤ N body and filament spheres by k > N body . Segment endpoints are located at x i for i = 1, 2, …, N, and x 1 is defined relative to x b in the local basis Each segment is oriented in space according to its local director basis ½d i 1 , d i 2 , d i 3 , for i = 1, …, N. As such, subsequent segment positions are given by for i = 2, 3, ..., N, where Δs = L/N, ensuring the filament remains exactly inextensible. Invoking force-free and torque-free conditions, the total momentum balance of the elastohydrodynamic filament-body system in equations (2.2) and (2.3) reads where F i and T i are the total force and torque, respectively, exerted on sphere i by the fluid. Integrating equation (2.3) from s = s j to s = L yields for j = 2, …, N, where m j is the bending moment at s j defined constitutively by where E b is the bending stiffness, E t is the torsional stiffness, k j a are the curvatures at s j for α = 1, 2, 3, and k 0 a ðs j Þ is the preferred curvature at s j . If k 0 a ðs j Þ ¼ 0 for α = 1, 2, 3, the filament will relax to a straight configuration. By making k 0 a ðs j Þ non-zero, the filament will relax to non-straight configurations, and by making k 0 a ðs j Þ a travelling wave, the filament can be made to swim. The CG governing equations (2.5)-(2.7) can be written in a simpler and compact matrix form, given by (2.9) where we have used the matrix form of the cross product, denoted by ½a Â (see appendix A), and defined Δ i = y i − 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 A and B multiply forces and torques experienced by the spheres in the filament, and are size 3(N − 1) × 3Nn. The matrix A comprises (N − 1)Nn sub-matrices given by A i,j ¼ ½y N body þj À x iþ1 Â , if j . in, and zero otherwise. The matrix B has the same structure as A with sub-matrices given by B i,j = I, if j > in, and zero otherwise.
The CG matrix formulation in equation (2.9) encodes the force and torque balance for a single filament attached to an arbitrary solid body (figure 2b). By setting N body = 0 or N = 0, we model a free filament or a free solid body, respectively, further highlighting the simplicity of the formulation. The general form of the CG formulation in equation (2.9) for an arbitrary number of solid bodies each with an arbitrary number of attached filaments is simply given by where M F is analogous to the matrix in equation (2.9), F and T contain the forces and torques, respectively, exerted on all spheres by the fluid. The vector K is analogous to the right-hand side in equation (2.9) and defines the constitutive relations for every 3D Kirchoff filament in the system. In appendix B, we provide the explicit form of equation (2.10) and include a practical example of the CG matrix in the case of two free filaments. It is worth noting that the CG formulation in equation (2.10) is general and does not depend on the specific choice of hydrodynamic coupling. For this reason, equation (2.10) 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 RPY non-local hydrodynamic approximation [41] for low Reynolds number fluids.

Non-local hydrodynamic coupling and dimensionality reduction
We approximate the non-local hydrodynamic forces and torques on each interacting sphere using the RPY tensor [41]. This relates linear velocities v i and angular velocities ω i with forces and torques acting on each sphere via the non-local hydrodynamic mobility matrix M H , giving where V and V contain the velocities v i and angular velocities ω i of all spheres, respectively. 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 η and mutual distance, see appendix C for details. We note that as spheres begin to approach each other, the RPY approximation begins to break down, and thus additional lubrication corrections may need to be included. Zuk et al. [41] discuss the RPY approximation and its use in macromolecular bead models. The formulation presented here allows for any hydrodynamic method to be used instead of the RPY approximation, so depending on the specific problem a different hydrodynamic method may be employed. Exploiting the invertibility of the non-local hydrodynamic mobility tensor above, the generalized CG formulation (equation (2.10)) for the non-local elastohydrodynamic coupling simply reads However, the governing equation (2.12) lists, unnecessarily, all linear v i and angular ω i velocities for each filament's spheres as unknowns. Within the same body-filament structure, translation and rotations of each 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 (figure 2), see appendix D for details. As a result, v i and ω i are a function of the velocity of the centre of mass of the body _ x b , its angular velocity ω b , and the angular velocity of each segment v seg i , from the same body-filament structure. By rigidly attaching the first segment to the body (figure 2), we set ω 1 = ω b . Altogether, this reduces the dimensionality of the system from 6(N body + Nn) to 3 + 3N in equation (2.12). The dimensionality reduction is achieved via matrix Q, where W is the reduced velocity system for all spheres in terms of i , for i > 1 (see appendix D for details). Substituting equation (2.13) into (2.12) gives the full non-local CG elastohydrodynamic system in 3D, where W contains the reduced unknown velocities of all interacting units of the elastohydrodynamic system. Knowledge of W provides 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 relative to the laboratory fixed frame of reference, for example, using the rate of change of the directors base dd i /dt = ω × 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.

Tracking director rotations in three dimensions using exponential mapping of quaternions
Euler angles are commonly used to track 3D rotations of a vector basis as a function of three parameters. The challenge of any three-parameter parametrization is the existence of coordinatesingularities. 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]. Coordinatesingularities can be circumvented by continuously changing the laboratory frame of reference [31], 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, and have a history of use in the literature [9,[43][44][45]. Using a naive royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20230021 quaternion implementation 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 re-parametrizes 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 rescaling 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 the CG formulation introduced above. Next, we will briefly outline quaternions and exponential mapping, and their coupling with CG formulation. Quaternions are an extension of the complex plane into a complex fourdimensional space. A quaternion is given by The axis of rotation is given by the imaginary part of q and the angle of rotation by the real part of q. Specifically, gives the quaternion corresponding to a rotation about the axisn by θ. The angular velocity of a set of basis vectors is related to the quaternion by providing a simple framework for numerical integration of quaternions, though this increases the dimensionality of the system. Generalizing above to account for the multiple angular velocities of the system W, the state shape-vector X quat evolves in time according to where the matrix C is given in appendix E and is solved together with governing equation (2.14). The exponential mapping, on the other hand, re-parametrizes quaternions back to a three-parameter parametrization [37]. We define a pure vector known as the generator r, which is related to a quaternion via the exponential map, q ¼ expðrÞ ¼ ½cos jrj sincjrjr``, where sinc x ¼ sin x=x. The angular velocity of a set of basis vectors defined by generator r can be written conveniently in matrix form [46], where the matrix D is given by The determinant jDj ¼ sinc 2 jrj goes to zero as |r| → nπ, for integer n ≥ 1. As such, a rescaling is necessary whenever the magnitude of a generator approaches nπ to avoid coordinatesingularity [37]. This can be achieved by rescaling whenever |r| ≥ π/2. This rescaling does not change the rotation performed by the generator r but keeps it far away from any coordinate-singularity. Equation (2.18) can be generalized to account for the angular velocities of each sphere in the system W so that the state shape-vector X gen can be evolved in time using which can be substituted, directly, into the governing system of equations (2.14). This provides a convenient matrix form to evolve in time the non-local, multi-filament--body elastohydrodynamic ODE system After prescribing an initial state vector and boundary conditions, the governing system of equations, equations (2.14) and (2.17) for quaternions, and equation (2.22) for exponential mapping, can be solved directly for the time evolution of the system's state vector. We approximate curvatures in K using finite differences of the director basis (using equation (2.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 implementation is provided here [47] to serve as a basis for rapid customization and further generalizations. We hope that the simplicity of the matrix system allows researchers to use a numerical platform of their choice for quick implementation. For results taken using several parameter values, we make use of a high-performance computing system to run multiple simulations simultaneously.

Exponential mapping of quaternions and comparison with partial differentional equation formulations
We proceed with a comparison between the computational efficiency of quaternions and exponential mapping within the CG framework. We also compare our simulations with previous methods with similar levels of elastohydrodynamic accuracy [9,22], 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 3D that are straightforward to implement using the CG formulation. We non-dimensionalize the system using length-scale L, timescale T given either by the characteristic frequency ω 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, meaning that the body and first segment of filaments rotate with equal angular velocity, and may not bend or twist relative to each other. The dimensionless system of equations is thus multiplied by the dimensionless stiffness parameter S 4 ¼ L 4 ðhv=E b Þ, where η is the fluid viscosity. The stiffness parameter is closely related to parameters such as the elastohydrodynamic number E h ¼ 8pS 4 , and the sperm number Sp 4 ¼ ððj ? =hÞS 4 Þ, where ξ ⊥ denotes the perpendicular resistive force drag coefficient from RFT used elsewhere [12,31], and characterizes the ratio of viscous and elastic forces in the system. For most problems, the stiffness parameter generally lies between S ¼ 1 and S ¼ 20. As the stiffness parameter is lowered in elastohydrodynamic systems, especially below S ¼ 1, solving the system becomes more numerically stiff and will take more time to solve.
We first compare the computational performance between quaternion and exponential mapping for the simple case of a royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20230021 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, x b = x 1 , q b = q 1 and r b = r 1 . Initial quaternions are given by q i ¼ ½cosðu i =2Þ sin ðu i =2Þn``and generators by This corresponds to a filament bent into a semicircle, as shown in figure 3a for N = 40, n = 2. The initial shape unbends towards the unstressed straight configuration, with the centre of mass position conserved to within 2 × 10 −2 (L). We set n = 2 and vary the number of segments N. Figure 3b 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 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 parametrization runs over 150 times faster. This increase in runtime is true across all ranges of S. Hereafter, we proceed with 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. [9] simulate a single helically actuated filament, where one end of the filament is held a constant distance δ 0 and angle α 0 from an axis of rotation. The base of the filament is rotated at an angular velocity ω about the axis of rotation. We implement this using a filament with N = 20 and n = 1. Initial generators are all set to r i ¼ a 0 =2½0 1 0`. To enforce the motion of the base of the filament, we replace the force-free condition with the kinematic constraint _ x 1 ¼ Àd 0 ½sin t cos t 0`, and the torque-free condition with the kinematic constraint _ r 1 ¼ ða 0 =2Þ½cos t À sin t 0`. We solve the system between t = 0 and t = 10 000, 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 3c shows the tip distance compared with results from Schoeller et al. [9] for δ 0 = 0.1 (L) and δ 0 = 0.02 (L) with α 0 = 0.2618. Results agree excellently with both Schoeller et al. [9] and experimental data from Coq et al. [48]. It is worth noting that Schoeller et al. [9] has also been validated and contrasted with experiments independently. The agreement with other elastohydrodynamic methods in 3D, such as in Schoeller et al. [9], using the moment balance PDEs directly, together with Lagrange multipliers to enforce the inextensibility constraint, indicates equivalence between canonical PDE methods and the simpler CG method employed here in 3D, as also demonstrated in 2D by Moreau et al. [12].
We conduct one further test that places significance on the non-local hydrodynamics between spheres in the system. Following Delmotte et al. [22], we simulate the case of a single swimmer and then two identical swimmers coupled via nonlocal hydrodynamics. We prescribe a time-dependent intrinsic curvature to a filament made of N = 16 segments with n = 1 sphere per segment so that our hydrodynamic description matches the beads model employed in Delmotte et al. [22]. We use the curvature κ 0 (s, t) = −K 0 sin(ks − ωt)d 1 (s, t), with K 0 = 8.25, if s ≤ 0.5L, and K 0 = 16.5(L − s)/L, if s > 0.5L, and a fixed stiffness parameter S ¼ 22:6 1=4 . Numerical parameter values are chosen to match Delmotte et al. [22], which itself validated this swimming case with experiments [49] and simulations [50]. 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 Delmotte et al. [22] and V = 0.07 (L/T ) from experiments in Bilbao et al. [49].
We also measure the swimming speed of two active filaments swimming together, figure 3d, both with the same time-dependent curvature as the single swimmer above, but set to either in-phase or anti-phase (inset of figure 3d 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 normalized swimming speed of the two-swimmer system is given in figure 3d, 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 3d shows excellent agreement with Delmotte et al. [22] for in-phase and anti-phase swimmers. This is despite major differences in both modelling and computational formulations. In particular, Delmotte et al. [22] 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 Zuk et al. [41], as the beats of the two swimmers begin to overlap slightly. Using the modified version of the tensors ensures that the hydrodynamic mobility matrix remains positive-definite and non-singular if the separation between two spheres approaches zero, and does not affect runtime significantly. Additionally, as swimmers begin to touch, steric repulsive forces begin to play an important role, which may be added to the methodology, as we outline in the discussion.

Bi-flagellate Chlamydomonas elastohydrodynamicswimming in three dimensions
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 [51]. 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, while the whole system is free from external forces and torques during free-swimming motion. As such, Chlamydomonas modelling is often abstracted with minimal models [52][53][54]. There is a large amount of literature studying bi-flagellated swimmers using far from minimal models, with many prescribing the shape of the flagella in the body frame to solve for the swimming kinetics in the laboratory-fixed frame of reference, but not solving for the shape of the flagella [55][56][57]. There is also work where the full elastohydrodynamic system is solved; Fauci et al. [58,59] modelled Chlamydomonas using the immersed boundary method, and more recently Nourian & Shum [44] developed a numerical method to study bi-flagellated bacteria. The coarse-grained Chlamydomonas swimmer is made of two active elastic filaments with an angle 2θ between them, attached to single spherical body of radius R, as depicted in figure 4a. The body is made out of N body = 184 spheres and the filaments contain N = 15 segments with n = 1 sphere per segment. With S ¼ 3, we induce straight swimming by prescribing a time-dependent preferred curvature, k 0 ¼ +4½1 þ sinð2ps À tÞd 2 ðsÞ in each filament. We note that this curvature is the preferred curvature of the system that contributes to the bending moment, and the emergent waveform curvature depends on the total balance of moments on the flagella, as outlined in the Methods section. The constant term in the prescribed curvature induces a non-zero average curvature along the flagella (figure 4b), 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. Figure 4c shows the swimming 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 is 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 figure 3d, 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 decreases with the body size in figure 4c, indicating that Stokes law for the spherical body could provide a good approximation when the body is large (or the flagella are small). This is because the non-local hydrodynamic effect becomes increasingly small, when compared with the zerothorder contribution from Stokes law, as the body size decreases.
Motivated by recent empirical studies observing threedimensional beating of Chlamydomonas flagella [54,60], 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 k 0 1 in the d 1 direction of the two flagella, so that the time-dependent preferred curvature is k 0 ¼ k 0 1 d 1 ðsÞ+ 4½1 þ sinð2ps À tÞd 2 ðsÞ. In the symmetric case, when k 0 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-ofplane curvature k 0 1 , shown in figure 4d. Interestingly, any incongruity 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 [54,60].
By adding an asymmetry between the two flagella's outof-plane curvatures, the circular paths seen in figure 4d turn into helical paths shown in figure 4e. If instead the two outof-plane curvatures are equal and opposite, the Chlamydomonas swims along a straight trajectory and rolls about its own axis, shown in figure 4f. 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; figure 4g), in further agreement with 3D observations [54,60].

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 in royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20230021 biological soft-matter and fluid-structure research [3]. Most of sperm swimming modelling focuses on hydrodynamic simulations using prescribed swimming beat patterns (kinematic constraints) to resolve swimming trajectories [61]. 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 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 [62,63].
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 [64]. The sperm flagella are 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 θ, the angle that the sperm cell makes with the z-axis at t = 0. Rather than prescribing a time-dependent intrinsic curvature, we prescribe an active internal moment density along the flagella. For this, we modify equation (2.3) accordingly, adding an extra term m a which is the prescribed active moment density, m s + x s × n + τ + m a = 0. When integrated, this modifies equation (2.7), giving an extra term on the right-hand side corresponding to the active moment density integrated from s j to L, We set the active internal moment density m a = 12k cos (ks − t)d 2 . This causes the waveform shown in figure 5a. Simulations are run with S ¼ 6. Figure 5b shows the initial condition for a single value of θ, and the subsequent trajectories of sperm cells towards the egg for four values of θ. 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 (figure 5c). Figure 5d shows the sperm swimming speed as a function of sperm-egg separation; initially separated by 5L, 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 freespace swimming speed (i.e. the swimming speed in absence of solid bodies), as the sperm-egg separation increases. These simulations indicate that hydrodynamic effects alone are unable to assist sperm to find eggs with similar sizes as  royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20230021 the sperm, even when swimming in close proximity. This may explain the need for other form of guidance mechanisms through evolution, such as chemotaxis in freshwater fish species [65]. Most interestingly, figure 5d 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 figure 5d are skewed, and above the average speed for certain sperm-egg separations, so that the accumulated distance travelled is higher for those cases (see purple and yellow curves in figure 5d). We highlight that similar effects of swimmer-surface interactions have been discussed previously [61,[66][67][68][69], and specifically speed-modulation of swimmers due to boundaries in [70]. In all, the exemplary results in figure 5 motivate further studies on the design and optimization of micro-environments to promote, or suppress, cell transport and progression using microfluidic systems.

Elastohydrodynamic cilia array and particle transport
Arrays of cilia are found in many places in our body and are responsible for a wide array of biological processes, from driving fluid flow during embryonic growth to clearing mucus in the lungs [5]. Elastohydrodynamic simulations of cilia array have been developed for decades to investigate several multi-scale aspects of this system, including studies on synchronization phenomena and metachronal waves [71][72][73][74][75]. 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 × 5 grid arranged on top of a wall made of spheres, as shown in figure 6a. Each cilium is made of N = 7 segments with n = 1 sphere per segment and governed by a prescribed curvature κ = (1 + 5sin(1.5πs − t + ϕ)) d 1 , figure 6b. Each row (cilia with same location in y) in the cilia array has a different value for the phase ϕ such that a metachronal wave is travelling in the positive y direction, figure 6a. Specifically, ϕ = (i − 1)/5 × 2π where i corresponds to the row number of the cilium, for i = 1, …, 5.
The beating of each individual cilium is shown in figure 6b over one period of oscillation. We note that one may use modified mobility tensors accounting for the hydrodynamics of a wall, which would speed up computation, 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 [73], 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 figure 6c,d from a side view and birds-eye view, respectively. The inset in figure 6c shows the oscillations of a single tracer particle over the beating period, and correlates with the 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 figure 6b, 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 figure 6c, in which particles tend to move upward (downward) for negative ( positive) y. The probability density function (PDF) of particle velocity in the y-direction for different initial heights in figure 6e shows that the particle transport is not effective for this ciliary beating. The transport velocity is only weakly shifted to positive values due to the small component of the beating in the direction of transport, with the average velocity of transport decaying with particle height, shown in figure 6f. Distributions of accessible velocities by the tracer particles vary noticeably with particle height due

Discussion and conclusion
We offer a simple and intuitive method for simulating fully 3D interacting elastic filaments coupled via non-local multibody microhydrodynamics (Matlab code provided [47]). By using asymptotic integration of momentum balance along each coarse-grained segment [12], 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-filamentbody systems involving moving, fixed and coupled-filament conditions, allowing easier treatment of cumbersome forcefree 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, each encoding distinct model interactions of the shape field _ X gen , namely: momentum balance M F , hydrodynamic coupling M À1 H , dimensional reduction Q, geometry of deformation and basis tracking D, and constitutive relations, internal activity and other boundary constraints K. Equation (5.1) is succinct and facilitates model explainability, generalizations and/or simplifications, as well as analytical and numerical analysis, better treatment of boundary conditions and easier implementation on different platforms. The formalism allows model customization without the need of algorithmic or numerical redesign. 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.
Adding non-trivial biologically inspired features of microswimmers is also straightforward with this formalism. For example, the elastic bending moment in equation (2.8) can be modified to allow for regions of different bending stiffness; the bending stiffness E b may be a function of arclength or depend on the plane of bending. Regions of inactivity may be included by, for example, making the active moment density m a in equation (4.1) zero for some range in arclength. If using corrections to the RPY tensors that allow for overlapping spheres, sphere radii can be readily changed along the length of the flagellum to account for varying filament radius. Segments could comprise spheres that are placed off the centreline, to incorporate unusual flagella geometries. Including additional forces and torque is also fairly straightforward. For example, they can be added to the right-hand side of the system of equations, giving where F Ã ðX Þ and T Ã ðX Þ contain additional forces and torques on spheres that are functions of the current state of the system. These may include repulsive steric forces that stop structures from overlapping, or attractive forces between spheres on different structures that bind filaments to other filaments or solid bodies. These forces can be switched on or off as a function of sphere separation or force magnitude, to incorporate the effect of new elastic bonds forming or breaking beyond a force critical value. Other types of forces arising from, for example, gravity, magnetism, electrostatic and royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20230021 stochastic interactions are straightforward to implement, offering further flexibility while employing this methodology. The method has been validated and contrasted against previous experimental and computational studies showing excellent agreement with the literature (figure 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. There are several reasons why the direct quaternion implementation runs slower than the exponential mapping in our framework. First and foremost, the quaternion parametrization requires solving for angular velocities first, and then calculating for the rate of change of quaternions. In the exponential mapping implementation, we solve directly for the rate of change of the exponential map. A large contribution to the runtime of the quaternion implementation is due to the variable timestep solver in Matlab, which adjusts the timestep in order to achieve the desired tolerance. With the quaternion parametrization, equation (2.16) implicitly contains the constraint that the quaternion remains unitary, whereas the exponential map gives a unit quaternion by definition, allowing larger timesteps in the exponential map implementation with the same tolerances.
We showcase the method in three complex elastohydrodynamic examples: (i) bi-flagellated Chlamydomonas swimming in 3D (figure 4), (ii) sperm-egg elastohydrodynamic coupling (figure 5), and (iii) canonical particle transport by cilia array ( figure 6). The elastohydrodynamic models and results presented for examples (i) and (ii) 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 (figure 4), in agreement with [54,60]. This offers numerous modelling opportunities motivated by recent 3D observations of this important model microorganism [54,60]. We also report novel sperm-scattering results in which sperm swimming can either be enhanced or reduced by the presence of a free body ( figure 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 (figure 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-organization of flagellar beat [76,77], elastohydrodynamic modulation of active bundles [78], non-trivial body geometries [79], physicochemical active guidance [65], ciliated microorganisms [80], cilia carpets [81], artificial swimmers [82,83], soft-robotics [84,85], locomotive robotics [86] and the role of torsional instabilities in filament systems [87], 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. As the size of the system becomes very large, methodologies which are optimized for large systems [9] will probably be more computationally efficient, but we note that significant runtime improvements could be achieved by tailoring the numerical scheme to specific applications, for example, using adaptive coarse-graining similar to the adaptive mesh used in [88], reducing hydrodynamics to RFT in single swimmer studies as in [27], in combination with using specific integration schemes [89] rather than built-in solvers used here, as well as using features such as Matlab MEX functionality to calculate the hydrodynamic mobility matrix. 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.
facilities of the Advanced Computing Research Centre, University of Bristol: http://www.bristol.ac.uk/acrc/. We thank the anonymous referees for suggested improvements in this work.

Appendix A. Matrix cross product
The cross product a × b in matrix form is given by Àa 3 a 2 a 3 0 Àa 1 Àa 2 a 1 0 For a single filament with no body, the block of zeros in the large matrix in equation (2.9) reduces to zero size, and the number of spheres in the system is given by M = Nn. If M F encodes the force and torque balance of a single filament, e.g.
then the force and torque balance of two free filaments is given by where (M F ) i , F i , T i and K i have identical structure to the single filament case and i = 1, 2 denotes the filament number. When writing in general form royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20230021 rows and columns in equation (B 2) are rearranged accordingly. Adding additional structures results in adding additional matrices along the diagonal in equation (B 2). For example, adding a single free sphere to the system would add six additional rows to 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 RPY approximation gives the hydrodynamic mobility matrix where t and r refer to translational and rotational components, respectively. Components of μ tt , μ tr , μ rt and μ rr are given by [41] m tt ij ¼ where a i is the radius of sphere i, η is the fluid viscosity, r ij is the separation of sphere i and j, r ij ¼ r ij ,r ij ¼ r ij =r ij and 1 is the Levi-Civita symbol. When i = j, m tt ij and m rr ij reduce to standard force and torque relations for a sphere in the Stokes regime, and there is no translation-rotation component. For i ≠ j, the tensors dictate how sphere j affects the movement of sphere i. For example, m tt ij couples the force on sphere j with the velocity of sphere i; m tr ij couples the torque on sphere j with the velocity of sphere i and m 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 x i is given by the velocity at x i−1 plus a contribution from the angular velocity of segment i − 1. Additionally, the velocity of spheres in the body are given by the velocity _ x b plus a contribution from the angular velocity of the body ω b . As such, the velocity of a sphere in the system is given by where Δs is the length of each segment, k = i − N body , dk=ne rounds k/n up to the nearest integer, l ¼ k À nðdk=ne À 1Þ and v seg i is the angular velocity of the ith segment in the filament. The quantity dk=ne gives the segment number that sphere i is attached to and l gives the sphere number on that segment (e.g. l ≤ n). The angular velocity of a sphere in the system is given by Equations (D 1) and (D 2) allow all v i and ω i in the system to be written in terms of velocity of the centre of mass of the body _ x b , its angular velocity ω b and the angular velocity of each filament segment v seg i . 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 (D 1) and (D 2) in matrix form gives equation (2.13) in the main text.

Appendix E. Implementing quaternion parametrization
Consider again the single filament attached to a body. Assign a quaternion q b to the body frame vectors and q i to filament segment frame vectors. The state vector for the system is The rate of change of the state vector is given by where P(q) is given by (see equation (2.16)) PðqÞ ¼ 1 2 Àq 1 Àq 2 Àq 3 q 0 q 3 Àq 2 Àq 3 q 0 q 1 q 2 Àq 1 q 0 2 6 6 4 3 7 7 5 : ðE 3Þ By rigidly attaching the first segment of the filament to the body, v seg 1 ¼ v b , the dimension of the vector on the righthand side of equation (E 2) can be reduced by 3, in addition to removing one column from the matrix pre-multiplying the vector. Generalizing to arbitrary systems gives the dynamics of the state vector X quat in terms of the reduced velocity system W via the matrix C, _ X quat ¼ CW: ðE 4Þ