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

The 3D coarse-graining formulation of interacting elastohydrodynamic filaments and multi-body microhydrodynamics | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results The 3D coarse-graining formulation of interacting elastohydrodynamic filaments and multi-body microhydrodynamics View ORCID Profile Paul Fuchter , View ORCID Profile Hermes Bloomfield-Gadêlha doi: https://doi.org/10.1101/2023.01.17.524366 Paul Fuchter Department of Engineering Mathematics and Bristol Robotics Laboratory, University of Bristol , Bristol, UK Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Paul Fuchter Hermes Bloomfield-Gadêlha Department of Engineering Mathematics and Bristol Robotics Laboratory, University of Bristol , Bristol, UK Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Hermes Bloomfield-Gadêlha For correspondence: hermes.gadelha{at}bristol.ac.uk Abstract Full Text Info/History Metrics Preview PDF 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 – 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 – 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. Download figure Open in new tab FIG. 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. 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. Classical formulations of solving filament elastohydro-dynamics 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 elasto-hydrodynamic 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 coarsegraining 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 – 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 elastohy-drodynamic 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 coarsegraining formulation is efficient, intuitive and easy-to-implement and customise. Our method allows for simulation of 3D elastic filaments with both flexural and tor-sional 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 microenvironment. 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. A. The coarse-graining formalism in 3D We consider first a single Kirchoff filament attached to an arbitrary solid body ( Fig. 2(a) ). 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 where κ = κ 1 d 1 + κ 2 d 2 + κ 3 d 3 + The constituents of the twist vector κ 1 , κ 2 and κ 3 are known as the filament curvatures [ 17 ]. The third curvature in the 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 κ 3 torsion and twist. Download figure Open in new tab FIG. 2. (a) An elastohydrodynamic structure comprised of solid body and filament. The filament is described by centreline x ( s ) with director basis [ d 1 ( s ), d 2 ( s ), 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 x b and body director labelled. Segment position x i labelled for as well the the director basis for segments 1 and 3. The balance of forces and torques along the filament at low Reynolds number reads [ 40 ] 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 coarsegraining 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, x b , described by the local basis ( Fig. 2(c) ). Sphere positions are given by 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 ≤ N body and filament spheres by i > 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 . Invoking force- and torque-free conditions, the total momentum balance of the elastohydrodynamic filament-body system in Eqs. ( 2 , 3 ) reads where F i and 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 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 and 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 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) × 3 Nn . The matrix A is comprised of ( 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 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 where is analogous to the matrix in Eq. 8 , and contain the forces and torques, respectively, exerted on all spheres by the fluid. 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. B. 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 v i and angular velocities ω i with forces and torques acting on each sphere via the non-local hydrody-namic mobility matrix , giving where and Ω contain the velocities v i and angular velocities ω 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 η and mutual distance, please see Appendix C for details. Exploiting the invertibility of the non-local hydrodynamic mobility tensor above, the generalised coarsegraining formulation ( Eq. 9 ) for the non-local elastohy-drodynamic coupling simply reads However, the governing Eq. 11 lists, unnecessarily, all linear v i and angular ω 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, v i and ω i are a function of the velocity of the centre of mass of the body , its angular velocity ω b , and the angular velocity of each segment , part of the same body-filament structure. By rigidly attaching the first segment to the body ( Fig. 2 ), we set ω 1 = ω b . Al-together, this reduces the dimensionality of the system from 6( N body + Nn ) to 3 + 3 N in Eq. ( 11 ). The dimensionality reduction is achieved via where is the reduced velocity system for each sphere in terms of for i > 1. Substituting Eq. ( 12 ) into ( 11 ) gives the full non-local, coarse-graining elasto-hydodynamic system in 3D where contains the reduced unknown velocities of each interacting unit of the elastohydrodynamic system. Knowledge of 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 d d i /d t = ω × 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. C. 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. 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 q = q + q 1 i + q 2 j + q 2 k = [ q 0 q 1 q 2 q 3 ] ⊤ , where i 2 = j 2 = k 2 = ijk = −1. Unit quaternions can be used to rotate a vector basis v in space via v ′ = qvq −1 . 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 axis 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. Generalising above to account for the multiple angular velocities of the system , the state shape-vector evolve in time according to and is solved together with governing Eq. (13) , where the matrix is given in Appendix. E. The exponential mapping, on the other hand, reparameterises 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 | r | sinc | r | r ⊤ ] ⊤ , 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 [ 44 ], where the matrix D is given by The determinant | D | = sinc 2 | r | 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 coordinate-singularity [ 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. Eq. ( 17 ) can be generalised to account for the angular velocities of each sphere in the system so that the state shape-vector can be evolved in time using 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 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 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 ω 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 . where η is the fluid viscosity. The stiffness parameter is closely related to parameters, such as, the elastohydrodynamic number E h = 8 π S 4 , and the sperm number where ξ ┴ 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, x b = x 1 , q b = q 1 and r b = r 1 . Initial quaternions are given by and generators by , where θ i = ( i −1) π /( N −1) and , 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 × 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 . 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 . Hereafter, we proceed our numerical study employing exponential mapping for more efficient computation. Download figure Open in new tab FIG. 3. (a) Relaxation of filament with N = 40, n = 2 and . 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 . 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 ]. 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 δ 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 = α 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 , and the torque-free condition with the kinematic constraint, . 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 δ 0 = 0.1( L ) and δ 0 = 0.02( L ) with α 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 , with K 0 = 8.25, if s ≤ 0.5 L , and K 0 = 16.5( L − s )/ L , if s > 0.5 L , and a fixed stiffness parameter . 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 A. 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 – 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. 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 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 , we induce straight swimming by prescribing an intrinsic time-dependent curvature κ = ±4[1 + sin(2πs - t)] 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) , in-indicating that Stokes law for the spherical body could provide a good approximation when the body is large (or the flagella is small). dicating 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. Download figure Open in new tab FIG. 4. (a) Model Chlamydomonas setup. Angle 2 θ between two flagella attached to body of radius R . Here, 2 θ = 23 degrees and R = 0.35( L ). (b) Periodic bending wave given by time-dependent intrinsic curvature κ = ± (4 + 4sin(2πs - t)] d 2 shown over one period T . Induces straight swimming of Chlamydomonas indicated by arrow. The 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 κ 1 in the d 1 direction causes out of plane swimming. Trajectory of Chlamydomonas when both flagella have equal out of plane intrinsic curvature κ 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 κ 1 = 0.5 ± 0.3 (f) Trajectory of Chlamydomonas when both flagella have equal and opposite out of plane intrinsic curvature κ 1 = ±0.5 (g) Asymmetric opposite sign out of plane curvatures κ 1 = ± 0.5 + 0.3 turns straight swimming with rolling into a helical trajectory, with rolling still present. 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 κ 1 in the d 1 direction of the two flagella, so that the intrinsic time-dependent curvature is κ = κ 1 d 1 ( s ) ± 4[1 + sin(2πs - t)] d 2 (s) In the symmetric case, when κ 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 κ 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 ]. B. 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.2 L and the egg as a sphere of radius 1 L . 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 5 L , 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 set an active internal moment density along the flagella, m = 12 k cos ( ks − t ) 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 Fig. 5(b) 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 ( Fig. 5(c) ). Download figure Open in new tab FIG. 5. (a) Waveform of model sperm shown over one beat period T. (b) Trajectory of sperm past egg with different initial condition for θ , the angle between the flagellum and . (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 θ , the sperm reaches a minimum separation from the egg before swimming away and increasing in speed. 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. C. 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 – 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 × 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 κ = (1 + 5 sin(1.5 πs − t + ϕ )) d 1 , Fig. 6(b) . 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, Fig. 6(a) . 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 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. Download figure Open in new tab FIG. 6. (a) Initial condition for 5 × 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 π 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). 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 straight-forward 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, by-passing 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 , namely: momentum balance , hydrodynamic coupling , dimensional reduction , geometry of deformation and basis tracking , and constitutive relations, internal activity and other boundary constraints . 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) biflagellated 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 elastomicrohydrodynamics 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. 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 . OpenUrl [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 . OpenUrl [11]. ↵ Moreau C , Giraldi L , Gadêlha H . The asymptotic coarsegraining 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 . OpenUrl Abstract [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 . OpenUrl [29]. ↵ Walker BJ , Gaffney EA . Regularised non-uniform segments and efficient no-slip elastohydrodynamics . Journal of Fluid Mechanics . 2021 ; 915 : A51 . OpenUrl [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 . OpenUrl [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 differentsized 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 . OpenUrl [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: Spermtemplated 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 ). 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 a × b in matrix form is given by 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 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 , , and have identical structure to the single filament case and i = 1, 2 denotes the filament number. When writing in general form rows and columns in Eq. B2 are rearranged accordingly. Adding additional structures results adding extra matrices along the diagonal in Eq. B2 . For example, adding a single free sphere to the system would add 6 additional rows to , 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 where t and r refer to translational and rotational components, respectively. The submatrices have general form Components of μ tt , μ tr , μ rt and μ rr are given by [ 41 ] where a i is the radius of sphere i , η is the fluid viscosity, r i,j is the separation of sphere i and j and ∊ is the Levi-Civita symbol. When i = j , and reduce to standard force and torque relation 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, couples the force on sphere j with the velocity of sphere i ; couples the torque on sphere j with the velocity of sphere i and 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 is given by the velocity 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 , ⌈ k/n ⌉ rounds k/n up to the nearest integer and is the angular velocity of the i th segment in the filament. The quantity ⌈ k/n ⌉ gives the segment number that sphere i is attached to. The angular velocity of a sphere in the system is given by Equations (D1) and (D2) allow all v i and ω i in the system to be written in terms of velocity of the centre of mass of the body its angular velocity ω b and the angular velocity of each filament segment . 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 + 3 N velocities. Writing equations (D1) and (D2) 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 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 Eq. 15 ) By rigidly attaching the first segment of the filament to the body, , the dimension of the vector on the RHS of Eq. E2 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 in terms of the reduced velocity system via the matrix Back to top Previous Next Posted January 19, 2023. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following The 3D coarse-graining formulation of interacting elastohydrodynamic filaments and multi-body microhydrodynamics Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share The 3D coarse-graining formulation of interacting elastohydrodynamic filaments and multi-body microhydrodynamics Paul Fuchter , Hermes Bloomfield-Gadêlha bioRxiv 2023.01.17.524366; doi: https://doi.org/10.1101/2023.01.17.524366 Share This Article: Copy Citation Tools The 3D coarse-graining formulation of interacting elastohydrodynamic filaments and multi-body microhydrodynamics Paul Fuchter , Hermes Bloomfield-Gadêlha bioRxiv 2023.01.17.524366; doi: https://doi.org/10.1101/2023.01.17.524366 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Biophysics Subject Areas All Articles Animal Behavior and Cognition (8103) Biochemistry (18917) Bioengineering (15076) Bioinformatics (44918) Biophysics (22805) Cancer Biology (19936) Cell Biology (27143) Clinical Trials (138) Developmental Biology (14065) Ecology (21240) Epidemiology (2067) Evolutionary Biology (25638) Genetics (16270) Genomics (23674) Immunology (18877) Microbiology (42910) Molecular Biology (18221) Neuroscience (94388) Paleontology (708) Pathology (3012) Pharmacology and Toxicology (5147) Physiology (8206) Plant Biology (16151) Scientific Communication and Education (2106) Synthetic Biology (4604) Systems Biology (10326) Zoology (2405) window.__CF$cv$params={r:'a47b8ceb8b730ade',t:'MTc5MTUyOTY2OA==',u:'01a11f7d5b767028b74fa254962eeda2',ut:'tu_EkfZg.Zv.VT2Xi.DtFK0CE3LkinTxlPKcvBqXkt0-1791529671-1.2.1.1-Ac6neY0RVarLaeNtThzYEv51WgozikUTQYvfTLGJpsuANsYwILicGF6Vf1ABiMYH0ntltopGOktCTI3cwzg1PUQeGbRyw2hIScxUOCIqm.o',i:60};(function(){if(!document.body)return;var s=document.createElement('script');s.src='/cdn-cgi/challenge-platform/scripts/precursor/main.js';document.head.appendChild(s);})();