Filament mechanics in a half-space via regularised Stokeslet segments

We present a generalisation of efficient numerical frameworks for modelling fluid-filament interactions via the discretisation of a recently-developed, non-local integral equation formulation to incorporate regularised Stokeslets with half-space boundary conditions, as motivated by the importance of confining geometries in many applications. We proceed to utilise this framework to examine the drag on slender inextensible filaments moving near a boundary, firstly with a relatively-simple example, evaluating the accuracy of resistive force theories near boundaries using regularised Stokeslet segments. This highlights that resistive force theories do not accurately quantify filament dynamics in a range of circumstances, even with analytical corrections for the boundary. However, there is the notable and important exception of movement in a plane parallel to the boundary, where accuracy is maintained. In particular, this justifies the judicious use of resistive force theories in examining the mechanics of filaments and monoflagellate microswimmers with planar flagellar patterns moving parallel to boundaries. We proceed to apply the numerical framework developed here to consider how filament elastohydrodynamics can impact drag near a boundary, analysing in detail the complex responses of a passive cantilevered filament to an oscillatory flow. In particular, we document the emergence of an asymmetric periodic beating in passive filaments in particular parameter regimes, which are remarkably similar to the power and reverse strokes exhibited by motile 9+2 cilia. Furthermore, these changes in the morphology of the filament beating, arising from the fluid-structure interactions, also induce a significant increase in the hydrodynamic drag of the filament.


Introduction
The mechanics of flexible filaments on the microscale underpin much of biology, from the propulsive flagella of motile bacteria and spermatozoa to nodal cilia, the latter hypothesised to be responsible for the breaking of left-right symmetry in mammals (Berg & Anderson 1973;Gray 1928;Smith et al. 2019). Furthermore, the dynamics of elastic filaments are of intense interest in the physics of microdevices and surface flows, including near-wall dynamics. For instance, the consideration of soft deformable sensors has already motivated extensive studies of attached filaments (Guglielmini et al. 2012;Roper et al. 2006), as has the characterisation of attached filament forces for understanding the drag induced by slender appendages (Curtis et al. 2012;Pozrikidis 2011;Simons et al. 2014). Such appendages range from the primary cilium to carbon nanotube mats, with an extensive review of the field presented by du Roure et al. (2019), which notes that both theoretical and numerical developments are very much still required in this field. Indeed, with advances in microscopy enabling ever more detailed quantification of kinematics, often with confining geometry such as a cover slip or substrate, the development of validated, and ideally simple, methodologies would be beneficial in estimating mechanics from kinematics.
Furthermore, the complex mechanics of fluid-structure interaction is an important problem and has been well studied, as illustrated by the slender body theory of Tornberg & Shelley (2004) and Liu et al. (2018). Numerical simulations that attempt to move past slender body theory are frequently plagued by extensive numerical stiffness, such as the regularised Stokeslet simulations of Ishimoto & Gaffney (2018) and Olson et al. (2013). The recent advance of Moreau et al. (2018) and its subsequent extension by Hall-McNair et al. (2019) have sought to address this stiffness, with their methodologies significantly reducing the computational cost associated with filament-fluid interactions via the use of integrated force and moment balance equations, but have not considered even the simplest confined geometry that is an infinite planar wall. Hence our first objective will be to generalise these improved frameworks to dynamics in a half-space bounded by a wall, adapting the recent regularised Stokeslet segment approach of Cortez (2018) for drag calculations and elastohydrodynamics. We will validate in detail the computation of drag from kinematic data against the earlier work of Ramia et al. (1993), additionally validating our proposed elastohydrodynamic framework against the goldstandard boundary element method of Pozrikidis (2010Pozrikidis ( , 2011. These validations demonstrate the high accuracy of our approach in capturing the mechanics of filaments, and thus it is of extensive use for understanding cellular swimming. For instance, resistive force theories are very popular due to their ease of use (Gadêlha et al. 2010;Lauga et al. 2006;Moreau et al. 2018;Schulman et al. 2014;Sznitman et al. 2010;Utada et al. 2014), and they have been shown to be of reasonable accuracy in free space for small-bodied cells such as spermatozoa (Johnson & Brokaw 1979). With freedom to choose the resistive coefficients, such local drag theories have been shown to perform very well even near a boundary (Friedrich et al. 2010), as supported by the boundary element method validation studies for straight rods of Ramia et al. (1993). This suggests that the refinements of Katz et al. (1975) and Brenner (1962) to freespace resistive force theories, valid for straight filament motion parallel or perpendicular to a planar wall, may be useful for the analysis of microscopy data, and would be very popular if demonstrated to be accurate. In particular, in the common circumstance where subject cells are imaged swimming parallel to a coverslip, for example Friedrich et al. (2010); Riedel-Kruse & Hilfinger (2007), and thus have unchanging boundary separation, we hypothesise that these refined resistive force theories may be sufficient to accurately capture the hydrodynamic drag on a curved filament. Hence, as our first application, we use the regularised Stokeslet segment framework of Cortez (2018) to test this hypothesis in exemplar problems, both when a filament is moving parallel to a boundary and when it is attached to the surface, subsequently moving extensively perpendicular to the wall.
Further, the presented elastohydrodynamic framework is sufficiently flexible to enable sophisticated considerations of fluid-structure interactions, inherited from the approach of Moreau et al. (2018). Hence, as an application of non-local elastohydrodynamics in a half-space, we proceed to demonstrate the methodology's ease of use in the context of such complex physics by examining the mechanics of a flexible inextensible adhered filament in an oscillatory flow, generalising Pozrikidis' study of carbon nanomat surface drag to time dependent flows (Pozrikidis 2011). Such simulations are clear precursors for the application of the framework to simulations of oscillating systems in physiological Figure 1: The two-dimensional piecewise-linear filament. We consider an inextensible filament composed of N straight segments each of length vi, connected at material points xi with each segment making angle θi with the global ex axis. The endpoints x1 and x N +1 shall be referred to as the base and tip respectively, and we note that, given segment lengths vi, the entire filament may be described by the location of the base and the angles θi. and soft matter modelling, such as primary and motile ciliary systems respectively in the kidneys and lungs, in addition to numerical studies of active surfaces, for example Balazs et al. (2014); Shum et al. (2013).

Piecewise-linear filaments
We consider a planar slender inextensible filament represented by N piecewise-linear segments of constant length and described by arclength parameter s ∈ [0, L], where L is the filament length.The segment endpoints are taken to correspond to material points, with their positions being denoted x 1 , . . . , x N +1 and where x i and x i+1 correspond to the i th segment. We will refer to x 1 as the base of the filament and correspondingly x N +1 as the tip, and denote the constant arclength associated with the material point x i as s i . As the filament is assumed to be planar, without loss of generality we assume that it lies in a plane spanned by orthogonal unit vectors e x and e y , with e z completing the orthonormal right-handed triad e x e y e z , and we may write x i = x i e x + y i e y . Following Moreau et al. (2018) we note that, once the segment lengths v 1 , . . . , v N are prescribed, the filament may be completely described by the N + 2 scalars x 1 , y 1 , θ 1 , . . . , θ N , where θ i is defined as the angle between the i th segment and the unit vector e x as shown in figure 1. Explicitly, from these N + 2 variables we may recover the filament endpoints as and, using dots to denote derivatives with respect to time, the velocity of each material point is given byẋ again expressible using only the N + 2 variables due to the imposed geometrical constraints.

Force distributions via regularised Stokeslet segments
To describe the low-Reynolds number fluid dynamics pertinent to the filament we utilise the recent work of Cortez (2018), namely the method of regularised Stokeslet segments (RSS). Here we briefly recapitulate the formulation of this methodology as presented by Cortez, relating the force density applied on the fluid by the filament to the fluid velocity via non-local hydrodynamics.
As in the popular method of regularised Stokeslets (Cortez 2001), we begin by considering solutions of the three-dimensional smoothly-forced Stokes equations, which may be stated for viscosity µ and force f applied on the fluid at the origin as where u is the fluid velocity, p is the pressure and φ is a smooth approximation to the Dirac delta distribution dependent on the small parameter . Following section 2 of Cortez (2018) we take for which the solution to equation (2.3) is known and given by (2.5) By linearity of the Stokes equations the velocity at a pointx due to a distribution of regularised Stokeslets along the filament with force density f (s) is therefore given by Whilst the original method of regularised Stokeslets would entail approximating the force distribution by a finite number of smoothed point forces, the method of regularised Stokeslet segments instead considers a linear distribution of forces along the straight segments of the discretised filament. With the discretisation of the force distribution as a continuous piecewise-linear function along each of the segments, the fluid velocity is instead given by a sum of integrals over each of the segments, where each individual integral may be analytically evaluated. Parameterising the j th segment by α ∈ [0, 1] so that x = x j − αv for v = x j − x j+1 , and additionally writing the force density as f = f j + α(f j+1 − f j ) for force densities f j , f j+1 at x j , x j+1 respectively, the integral along the j th segment is given by where we have identified |v| = v j and R = |x − x j + αv| 2 + 2 . As noted by Cortez which follows from the definition of T m,p and consideration of d dα (α m R p ). The simple extension of this methodology to a half-space bounded by a no-slip planar wall via the image system of Ainley et al. (2008) (with the typographical error corrected (Smith 2009)) is presented in detail by Cortez (2018, sec. 3.3). Exemplified in figure 2, we will consider two configurations of half space, one in which the infinite planar boundary is parallel to the plane containing the filament, given by z = 0, and another where the boundary is situated perpendicular to this plane and given by y = 0. In both cases the hydrodynamics may be described by the same regularised singularity representation, each requiring computation of T m,p for the additional values (m, p) ∈ {(0, −5), (1, −5), (2, −5), (3, −5), (4, −5), (5, −5)} in order to include the additional necessary regularised singularities. Omitted from the previous work of Cortez, we compute The remaining values of T m,p may be computed via the recurrence relation equation (2.9). Simple computation of the coefficients of the T m,p results in a matricial form of the velocity contribution at x from the j th segment as a linear operator acting on f j and f j+1 , the details of which are cumbersome and given in the Supplementary Material. The overall velocity at the evaluation pointx is then given by the sum of these contributions, resulting in a matricial equation of the form where F = (f 1,x , f 1,y , . . . , f N +1,x , f N +1,y ) T herein denotes the composite vector of force densities applied on the fluid at the segment endpoints x j in the directions e x and e y , where we write f (x j ) = f j,x e x + f j,y e y . Taking the evaluation pointx as each x i in turn yields a square system of linear equations in 2(N + 1) variables, relating the velocities u(x i ) to the force distribution on the filament via non-local hydrodynamics. Application of the no-slip condition at the segment endpoints then gives a relation between the filament velocities and the forces applied on the fluid by the filament, which we write in brief asẊ = AF , (2.12) for the square matrix A composed row-wise of blocks M (x i ) for i = 1, . . . , N + 1, and whereẊ = (ẋ 1 ,ẏ 1 , . . . ,ẋ N +1 ,ẏ N +1 ) T . Given kinematic data, this linear system may be readily solved to give the force densities applied on the fluid by the filament.

Extension to efficient solution of elastohydrodynamic equations
We now proceed to combine the work of Moreau et al. (2018) and Cortez (2018) to give an efficient scheme for solving non-local planar elastohydrodynamics, similar in concept to the piecewise-constant force density approach of Hall-McNair et al. (2019) but with continuous piecewise-linear force discretisation and the additional inclusion of an infinite planar boundary. As formulated by Moreau et al. (2018), we integrate the pointwise conditions of force and moment balance on the filament, given by for contact force and couple denoted n, m respectively and where a subscript of s denotes differentiation with respect to arclength, noting that we have assumed slenderness of the filament. This yields the integrated equations (2.16) Here we have assumed conditions of zero contact force and couple at the tip of the filament, i.e. n(L) = m(L) = 0, retaining generality at the base for the time being. With the assumption of piecewise-linear distributions of force density f over each segment, we may write these integrals as linear operators on F , yielding the system where R = (n(s 1 )·e x , n(s 1 )·e y , m(s 1 ), . . . , m(s N )) T and we are writing m(s i ) = m(s i )e z by planarity of the filament. The matrix B is of dimension (N + 2) × 2(N + 1), which under the assumption of equispaced segment endpoints has rows B k given by where i = 1, . . . , N and ∆s = L/N is the length of each segment. Thus we may write the coupled elastohydrodynamic problem as where B represents integration over segments and A encodes the relationship between the forces on the surrounding fluid and the velocities of material points on the filament, here non-local and assumed to be invertible. As noted in section 2.1, the material velocitiesẋ i may be expressed in terms of the time derivatives of the base point x 1 and the segment angles θ 1 , . . . , θ N . Defining Q to be the 2(N + 1) × (N + 2) matrix such that (2.23) where the subscript P denotes that the i th row of Q is to be permuted to (2.24) This permutation of Q allows us to define the blocks Q 1 and Q 2 as being strictly lower triangular matrices of dimension (N + 1) × N with entries Given the matrix Q we may rewrite the full coupled elastohydrodynamic problem simply as with BA −1 Q being a square matrix of dimension (N + 2) × (N + 2). Finally, and here for example assuming force and moment-free conditions at the filament base, the standard linear constitutive relation m(s i ) = EIθ s (s i ) ≈ EI(θ i − θ i−1 )/∆s for bending stiffness EI, valid for i = 2, . . . , N , gives R explicitly in terms of θ, yielding a linear system of low dimension that may be readily solved forθ, with the temporal dynamics then computable via existing ODE methods. As noted by Moreau et al., this formulation is readily extensible to the imposition of a variety of boundary conditions, in particular that of a cantilevered filament base, achieved by replacing the overall zero-force and zero-moment equations withẋ 1 =ẏ 1 =θ 1 = 0. We non-dimensionalise this system by scaling spatial coordinates with filament length L, forces with EI/L 2 , and time with some characteristic timescale T , yielding the nondimensional system for elastohydrodynamic number E h , where the notation· denotes non-dimensional quan-tities. The dimensional and non-dimensional quantities are related by where we have multiplied the two force balance equations by ∆s = L/N and absorbed the dimensional scalings of x 1 , y 1 into Q for convenience, withθ = (x 1 ,ŷ 1 , θ 1 , . . . , θ N ) T . In order to pose the non-dimensional equations in terms of the commonly-used sperm number Sp ( for resistive force coefficient ξ and L/ = 10 3 , giving the approximate relation E h ≈ 15.2 · Sp 4 .

Implementation, verification and parameter choice
Both the calculation of force densities from kinematic data and the solution of full elastohydrodynamics were implemented in MATLAB, the latter utilising the inbuilt stiff ODE solver ode15s (Shampine & Reichelt 1997). Prior to the recent work of Hall-McNair et al. (2019) the solution of non-local elastohydrodynamics has required significant computational work and minimal timestep for the solution of this stiff problem (Ishimoto & Gaffney 2018;Olson et al. 2013), and we replicate in our implementation the low computational cost associated with the integrated elasticity equations of Moreau et al. (2018). In particular, typical simulations of a cantilevered filament in background flow, explored in detail in sections 3.4 and 3.5, have a typical runtime of 10s, where N = 40 and we simulate over 10 periods of oscillation of the background flow on modest hardware (Intel R Core TM i7-6920HQ CPU).
Verification was first performed on the reduction to an unbounded domain, with the full elastohydrodynamics being verified by comparison with the resistive force theory results of Moreau et al. (2018) for filament relaxation. Further drag calculations were compared against the work of Cortez (2018) for the case of uniform motion. The regularised image system in a half-space was then verified by explicit evaluation of the fluid velocity on the no-slip boundary, with the numerical result being zero to machine precision, and further by noting that far-field results were seen to converge to those of the free-space system. Further qualitative checks were performed, an example being the successful reproduction of the intuitive result that relaxation timescales increase for filaments with reduced separation from the no-slip stationary boundary, along with numerical comparison against the previous works of Pozrikidis regarding filaments in steady shear flow (Pozrikidis 2010(Pozrikidis , 2011. Results of additional verification against the boundary element method of Ramia et al. (1993) are shown in figure 3(a). Sufficient accuracy is typically achieved with N = 20 segments, similar to the discretisation used by Cortez (2018), though we typically take N = 40 to enable the capturing of high-curvature filaments.
As noted by Cortez (2018), the use of regularised Stokeslet segments allows the regularisation parameter to represent the radius of the filament being modelled, verified here by comparison with the boundary element method of Ramia et al. (1993) in figure 3, subject to the phenomenological condition that be less than the length of the segments, which appears necessary for convergence. We will take /L = 10 −3 unless otherwise stated, with typical slender filament aspect ratios yielding /L in the range (10 −2 , 10 −3 ) (Ishimoto & Gaffney 2016;Yonekura et al. 2003). Whilst the results that follow naturally depend on the size the of the regularisation parameter, it being used as a proxy for the filament radius, this dependence appears to be predominantly quantitative, with the same qualitative conclusions holding across the range of physically-relevant values of /L.

Endpoint effects
Inherent to the method of regularised Stokeslet segments as presented here is the presence of apparent large variations in computed force densities at the tips of filaments. This was initially observed by Cortez (2018, sec. 3.1, fig. 4), who remarked that these endpoint effects may be reduced when /L is small. Indeed, for /L in our range of physical interest these endpoint effects contribute minimally to overall drag calculations, in particular having little effect on the computed total drag exerted on a filament and the computed force densities away from the filament tip, whose presentation we focus on herein. Exploration of the effects of non-uniform segment lengths yields the observation that these oscillatory endpoint effects remain limited to approximately the 3 segments proximal to the ends of the filament, thus the integral contribution of these oscillations may be reduced by the clustering of segments near the filament tips. For our typical parameters of N = 40 and /L = 10 −3 we observe a difference in integrated drag along of filaments of less than 2% between linear, quadratic and Chebyshev segment endpoints, hence we proceed with calculations utilising segments of uniform length.

Boundary-corrected resistive force theory
The leading-order approximation of slender body hydrodynamics that is resistive force theory (RFT) was first introduced by Gray and Hancock (Gray & Hancock 1955;Hancock 1953), relating the local drag on a body to its local tangential and normal velocities u t , u n via (2.31) Here f t , f n are the local tangential and normal components of force applied on the fluid, with the constant resistive coefficients C t , C n typically being functions of filament aspect ratio. As RFT relates forces only to local velocities, calculations using this simple local drag theory are not able to account for the presence of a boundary. Brenner (1962) and Katz et al. (1975) posed corrections to resistive force theory for straight filaments in asymptotic regimes far-from or near-to an infinite no-slip boundary, with the latter having been verified against high-accuracy boundary element methods by Ramia et al. (1993). In the case of filament motion confined to a plane parallel to the boundary, as exemplified in figure 2(a), by this wall-corrected resistive force theory (W-RFT) a straight filament parallel to the boundary has resistive coefficients given by where h is the boundary separation of the filament and has been taken as the filament radius. Similarly, and as in figure 2(b), when filament motion is in a plane perpendicular to the boundary the coefficients are given by where all coefficients are as summarised by Brennen & Winet (1977) and equations (2.33) and (2.35) additionally require h. In free-space we will adopt the resistive force coefficients defined by the large-h limit of equations (2.32) and (2.34), and refer to this simpler theory as free-space resistive force theory (F-RFT). As given by Brenner and Katz et al., the leading order boundary corrections to C t and C n are O (L/h) 3 when L h, and O (h/L) when L h, though the overall error in the resistive force approximation remains logarithmic in the filament aspect ratio.

Evaluation of wall-corrected resistive force theory for straight filaments
For the case of a straight uniform filament aligned parallel to a planar boundary, utilising the approach of section 2.2 we compute the hydrodynamic drag on the slender body as it moves parallel to the wall along its tangent at unit non-dimensional velocity, comparing the solutions given by the methods of regularised Stokeslet segments, freespace RFT and wall-corrected RFT. Having normalised by the F-RFT solution, it being independent of the boundary separation h, we show the computed non-dimensional total drag in figure 3. Good agreement near the boundary can be seen between the W-RFT of Katz et al. and our implementation of the method of regularised Stokeslet segments, the former as previously validated in the limit h/L → 0 with high-accuracy boundary element methods by Ramia et al. (1993). Indeed, direct comparison of the RSS solution with Ramia et al. (1993, fig. 4c), where we note that here we are taking /L = 10 −2 to ensure the filament radius matches that used in the boundary element computations, which thus provide additional verification of our methodology. Further, far from the boundary the correction of Brenner (1962) lies within 10% of the RSS solution, evidencing good overall agreement between the two schemes as the normalised wall separation h/L increases. Analogous agreement between these methodologies can additionally be seen for total normal drag when moving normal to the boundary, thus we conclude that the corrections to resistive force theory of Brenner and Katz et al. agree closely with regularised Stokeslet segments for straight filaments in their respective asymptotic regimes. Hence, agreement in more generality may be reasonably hypothesised.

Computing hydrodynamic drag on filaments parallel to boundaries
Following the established good agreement of W-RFT and RSS for straight filaments we proceed to evaluate the agreement for curved filaments moving in a plane parallel to the no-slip boundary. Motivated by the previous use of free-space resistive force theory for drag determination from captured kinematic data (Friedrich et al. 2010;Gaffney et al. 2011;Ishijima 2011;Ooi et al. 2014), we compute the forces on a free filament from the smoothed kinematic data of Ishimoto et al. (2017) corresponding to the flagellum of a human spermatozoon We first consider the filament moving in very close proximity to the boundary, with normalised separation h/L = 0.01, within the typical range of slithering motion for spermatozoa (Nosrati et al. 2015). A summary of the resulting drag calculations is shown in figure 4, from which we observe that the wall-corrected RFT of Katz et al. (1975) can demonstrate strong pointwise agreement with regularised Stokeslet segments over a single frame ( figure 4(a)). However, calculations of the total drag force over the filament highlight a general trend of over-estimation of force densities by W-RFT, with the differences between W-RFT and RSS having a median of 37% over the 100frame range shown here, with differences measured relative to the RSS value. Free-space resistive force theory appears to do the opposite, systematically underestimating the magnitude of total drag, with increased median deviation of 45% from the RSS solution.
Most notable however is the tight clustering of the effective ratio of drag coefficients C n /C t , shown in figure 4(c), computed pointwise from the RSS solution and with an analogous distribution of the effective values of C n /µ and C t /µ shown in figure 4(d). This suggests that an appropriate choice of resistive coefficient would enable the accurate approximation of local filament drag using only a local theory, in concurrence with the findings of Friedrich et al. (2010) though seen here much closer to the boundary and at a local level, rather than Friedrich et al.'s comparison to observed sperm behaviour, which necessarily averages over the cell. We verify this result using additional waveform data, modifying the kinematic data of Ishimoto et al. (2017) to crudely approximate a pinned spermatozoa by fixing both the endpoint and local tangent in space, along with the idealised pinned waveforms of Curtis et al. (2012). In particular, the tight clustering of effective coefficient ratios at low boundary separations is retained, though we note reduced agreement compared to that present for free-swimming data. Surprisingly, these coefficients and ratios do not align with the corrected coefficients of Katz et al. (1975), derived for straight filaments parallel to a boundary. Additionally, for approximate separations 0.3 < h/L < 1 we observe very strong agreement between resistive force theories and the non-local solution, with median differences in total drag between methodologies consistently less than 15%, in many cases being less than 5%. This is coupled with additional agreement concerning the direction of the resultant filament drag, which is retained even at reduced separations ( figure 4(b), lower). However, lost at such intermediate separations is any clear choice of effective resistive coefficient ratio, with the distribution of effective ratios significantly broadening above h/L ≈ 0.05.
At a much greater distance from the boundary, with h/L = 10, as expected F-RFT and W-RFT give approximately equal estimates for the drag on the free-swimming filament, with the W-RFT solution having approached that of F-RFT. As in the case of nearboundary swimming, only small differences are present between the W-RFT and RSS solutions, with median differences between methods of around 6%, in this case with the magnitude of all computed drag forces having been reduced from their near-wall values by approximately a factor of two. Thus, in the medium and far-field of a boundary resistive force theories appear remarkably accurate for determining the total drag on even curved filaments moving parallel to the boundary, though surprisingly at this increased boundary separation there is little grouping of the effective resistive coefficient ratios.

Hyperactivation-induced tugging of tethered spermatozoa
Explored initially by Curtis et al. (2012) using both free-space and wall-corrected resistive force theory, and reconsidered by Simons et al. (2014) using regularised Stokeslets, we re-evaluate the observation that the hyperactivation of mammalian spermatozoa aids in surface escape via a beat-induced tugging effect on the tether point of boundary-attached spermatozoa, and consider in detail the drag on the filament. Adopting the idealised beat patterns used by Curtis et al., appropriately non-dimensionalised, we position the base of the filament 0.01L from the boundary and compute both the total and local drag from this kinematic data for both hyperactivated and normal beating, noting that the plane of filament beating is perpendicular to the boundary. Mirroring the setup of Curtis et al., we assume that the filament is clamped at its base, implemented by rotating the kinematic data to align the basal tangent in each instant at some angle θ 0 to the boundary. Figure 5 shows the results of the drag computation over a single beat period for both the hyperactivated and normal beat patterns for θ 0 = π/2. We see reaffirmed by the method of regularised Stokeslet segments the conclusion of Curtis et al., with the hyperactivated beat pattern exhibiting a change of sign in force component perpendicular to the boundary, whilst no such change is observed for the normally-beating flagellum. Sampling θ 0 ∈ [π/4, π/2], we note that this observation holds across a range of basal orientations for these beat patterns, in agreement with the RFT-established conclusions of Curtis et al.. Overall, figure 5(a) demonstrates fair agreement between the local free-space resistive force theory solution and that obtained using regularised Stokeslet segments, suggesting a surprising validity in using simple local theories in this circumstance, as concluded broadly by Simons et al.. However, the pointwise values of the drag density highlight a stark disagreement between local and non-local theories, with deviations of up to 43% for the tangential component of the drag for the frame of hyperactivated beating shown in figure 5(b). For reference, the integrated total drag along the filament on average gives differences of only 24% between RSS and F-RFT, suggesting that a serendipitous cancellation occurs when integrating over the flagellum.
In figure 5(c) we show the ratio of effective tangential and normal drag coefficients for the hyperactivated beating, along with contours of zero curvature for reference. Significant variation in this ratio can be seen, in stark contrast to those computed for nearwall parallel swimming above. Thus the conclusion of Friedrich et al. (2010), that resistive force theory is accurate to high precision, does not hold for a bound spermatozoon, and hence the use of simple resistive force theories for near-boundary drag calculations is not reliable in general. This precludes the general use of the near-boundary correction of Katz et al. (1975) in this circumstance, supported by the poor agreement seen in figure 5, whilst the far-field result of Brenner (1962) is inappropriate in such close boundary proximity.

Morphological bifurcation of cantilevered filaments in an oscillatory flow
In the absence of a resistive force theory capable of accurately quantifying the details of pinned filament motion perpendicular to boundaries, we extensively utilise the non-local theory of regularised Stokeslet segments to consider the fully-coupled elastohydrodynamics of a cantilevered filament in an oscillating background flow, as formulated in section 2.3. With a planar infinite no-slip boundary situated at y = 0, we consider the non-dimensional background flow with velocityû b =â sin 2πt −t 0 ŷe x for amplitudeâ and phaset 0 , giving rise to the modified non-dimensional system where the components ofÛ b are given by the background flow evaluated at the endpoints of filament segments, explicitlŷ Figure 5: Drag computation for pinned spermatozoa with a base clamped perpendicular to a boundary over a single period of 100 frames. (a) Normalised total force on the boundary over a single beat period as computed by regularised Stokeslet segments (dot-dashed), free-space resistive force theory (dotted), and wall-corrected resistive force theory (solid), in the cases of hyperactive and normal beating patterns (high and low amplitudes respectively). Normal beating gives rise to a force on the boundary of unchanging sign, whilst hyperactivated beating generates a force whose sign changes over the beat period, corresponding to a so-called tugging effect. (b) Tangential (upper panel) and normal components (lower panel) of force density as a function of arclength for a single frame, shown for the hyperactivated beat pattern. In contrast to the fair agreement between methods seen in (a), with mean error of 24% for the hyperactivated beat between F-RFT and RSS, tangential force density differs on average by 43% between methods for the frame shown. (c) Ratio between effective drag coefficients computed using regularised Stokeslet segments, with contours of zero filament curvature superimposed as black curves, showing significant variation in the effective ratio. Colour online.
Noting that n(s 1 ) and m(s 1 ) are a priori unknown for a cantilevered filament, we impose the appropriate constraints ofẋ 1 =ẏ 1 =θ 1 = 0 in place of the total force and momentbalance equations, again yielding a square linear system. Under the assumption of a straight initial configuration with θ i = π/2 for i = 1, . . . , N , which we make throughout, there is freedom in the three non-dimensional parametersâ,t 0 and E h , the latter being equivalent to the sperm number Sp. In appendix A we consider in detail the effects on the dynamics of varying the phaset 0 of the background flow, concluding that phase can effectively be neglected in long-time dynamics, and thus we will fixt 0 = 0, noting that this breaks the inherent left-right symmetry of the initial condition. For fixedâ we consider the beating of the filament as it is driven by the background flow with the sperm number being varied. Withâ = 2π, we showcase in figure 6 the diverse range of observed behaviours. For low sperm numbers, in this case below a threshold value of around Sp = 2.5, we observe symmetric, low amplitude beating of the filament. Increasing the sperm number beyond this point results in the remarkable emergence of an asymmetric beating, shown in figure 6(b), strongly resembling the characteristic beating of a cilium. In particular, we see reproduced in this passive, flow-driven filament the defining power and recovery strokes of beating cilia (Brennen & Winet 1977).
Increasing the sperm number further results in the filament buckling and its beating regaining approximate symmetry, as can be seen in figure 6(c) for Sp = 4.55. For very high sperm numbers we see the appearance of buckling modes of higher order, resolved in figure 7 using N = 150 segments to capture filaments with such exceptionally high curvatures. With such high sperm numbers, we observe the collapse of modes onto those corresponding to lower-order buckling as time progresses, in addition to a greatly increased relaxation time to periodic beating in comparison to stiffer filaments.
In order to classify the mode of periodic beating we introduce the symmetry measure S, defined for our piecewise-linear filament with segment endpoint coordinatesx i (t) as where the integral overt runs over half the period of the motion after convergence to periodic beating has occurred. We note that the left-right symmetric beats of figures 6(a) and 6(c) give S = 1, with motion that breaks left-right symmetry such as the distinct forward and reverse strokes of ciliary beating yielding reduced values of S. We find that values of S a small tolerance less than unity indeed correspond well to the ciliary-type beating of figure 6(b). We additionally distinguish between approximately-straight and highly buckled beating by considering the number of filament regions with curvature of unchanging sign, classifying filaments with three or more such regions as buckled. Shown in figure 7(c) are the regions of the Sp-â parameter space that approximately correspond to low-amplitude symmetric, ciliary, and buckled modes of filament beating.
3.5. Time-averaged total drag on cantilevered filaments in oscillatory flow corresponds to distinct beating morphologies Computation of the total hydrodynamic drag exerted on the filament by the background oscillatory flow over a single period reveals a bifurcation structure closely aligned with that of the morphology of beating. Shown in figure 8(a) and rescaled by flow amplitude, the normalised total drag in the direction parallel to the boundary is seen to be of constant sign, and minimal for the approximately symmetric beating present at low sperm numbers, with a region of non-trivial total drag appearing for Sp ≈ 1 where elastic effects further distance the beating from reciprocal motion. figure 8(c) demonstrates that such minimal total drag arises as the result of gross cancellation over the period of oscillation, consistent with the symmetry of the associated beating mode. Upon entering the region of parameter space consistent with ciliary-type beating there is a significant change in the magnitude of the experienced total drag, intuitively correlated The curvature of the filament over the latter part of the initial period of the background flow, from which we may identify the collapse of high order modes onto lower order buckling modes. The frame shown in (a) is indicated on the curvature plot by a dashed line, with colour scalings consistent between plots. (c) A bifurcation diagram highlighting the regions of parameter space corresponding to low-amplitude symmetric, ciliary, and buckled modes of beating. Buckled modes appear to the right of the dividing dashed line, as defined by the transition between less than three and three regions of constant curvature sign along the filament. Furthermore, ciliary beating occurs within contours delimiting where S is essentially unity, shown shaded. In the overlap region between ciliary and buckled beating we observe ciliary-type beats with buckled tips. Colour online.
with the emergence of the distinct forward and backward strokes of ciliary beating. With reference to the bifurcation diagram of figure 7(c), highly-buckled beating at high values of Sp results in little total drag in the direction parallel to the boundary, consistent with the truly time-reversible motion obtained in the limit of equation (3.1) as E h → ∞.
Considering similarly the total force applied on the filament base in the direction normal to the boundary over a single beat period, we again observe that the total force is of constant sign, corresponding to the filament pushing towards the boundary. Shown in figure 8(b), the maximal force occurs around Sp = 1 for the majority of sampled flow amplitudes, aligning with the local maximum of total parallel drag noted previously and pertaining to the non-reciprocal symmetric beating of stiff filaments. Comparison with figure 8(d) reveals in this case that significant total drag cancellation occurs for the asymmetric ciliary-type beating. In contrast, the total drag on stiffer filaments in this direction sees limited cancellation, exemplifying the significance of beating morphology in filament drag calculations.

Discussion
In this work we have examined the accuracy of resistive force theories in quantifying the mechanics of planar filament motion in a half space via regularised Stokeslet segments, further exploring the coupled non-local elastohydrodynamics of cantilevered filaments in an oscillatory flow. We have seen that the corrections to free-space resistive force theory of Katz et al. (1975) and Brenner (1962) perform well both in the near and far-field of planar boundaries when applied to straight filaments, in agreement with the verification of Ramia et al. (1993) and serving as additional validation of our methodology.
By considering the motion of curved filaments in planes parallel to a boundary, in line with our hypothesis we have found that these corrections capture quantitative filament dynamics to moderate accuracy in all but the most extreme boundary proximities. Remarkably however, when very close to the boundary we have nonetheless observed a tightly-distributed ratio of effective drag coefficients, and indeed tightly-distributed coefficients, suggesting that a resistive force theory with such coefficients could yield accurate estimates for hydrodynamic drag in this circumstance. With subjects imaged Comparison with (a) highlights gross total drag cancellation occurring for low Sp, with reduced cancellation for the asymmetric ciliary-type beating. (d) Maximum absolute instantaneous force applied on the base of the filament in the direction perpendicular to the boundary. By comparison with (b) we observe limited cancellation of total drag in this direction in stiff filaments, with significant cancellation occurring for ciliary-type beating. Each panel is normalised by the maximum instantaneous total parallel drag and rescaled by flow amplitude to enable meaningful comparison. Select contours of total drag are shown and the region of ciliary-type beating is outlined in black, the latter determined only from the kinematic symmetry measure S. Colour online.
in the relative near-field of a coverslip, the similar conclusion of Friedrich et al. (2010) is thus in part supported by our non-local calculations, albeit here at greatly-reduced boundary separations, and suggests the future development of an accurate, empirical resistive force theory in the near-field of boundaries. This has potential application to the study of a slithering mode of swimming, as reported to be prevalent in spermatozoa by Nosrati et al. (2015) though this mode is far from ubiquitous. Analytical exploration of this result would require consideration of the asymptotic limit where the separation of singularities from their images is much less than the lengthscale of the filament radius of curvature, with such detailed calculations expected to be a subject of significant future study.
We have seen that the agreement between methodologies for parallel-swimming filaments does not carry through to those moving in a plane perpendicular to the boundary, with no clear consensus on effective coefficient ratio for pinned filaments moving perpen-dicular to boundaries. Thus we conclude that constant-coefficient resistive force theories cannot be expected to give accuracy comparable to non-local methods for filaments moving perpendicular to boundaries. Despite this, here we have reverified the conclusion of Curtis et al. (2012), based originally on resistive force theory, though this appears to rely on serendipitous cancellations on integrating along the filament.
In particular, we have noted a potential lack of reliability in local drag theories when applied to pinned filaments, as evidenced by the loss of coherence in effective drag coefficient ratio when using artificially-pinned kinematic data in section 3.2, supported further by the findings of section 3.3 for tethered spermatozoa. Accordingly, we have utilised the non-local hydrodynamics of Cortez (2018) to examine the response of a cantilevered elastic filament to an oscillatory background flow, with this application highlighting the flexibility in our treatment of the coupled elastohydrodynamics. Having established convergence to a single limiting periodic behaviour from a range of initial conditions for fixed filament and flow parameters, as detailed in the appendix A, we have evidenced in passive flow-driven filaments the existence of a remarkable mode of beating typically characteristic of actively-beating cilia. The fact that the filament movement appears qualitatively similar to that of actively-beating cilia further highlights that pumping fluxes are sensitive to the details of a ciliary beat with, here, zero total flux of fluid above the filaments, in contrast to the induced fluxes of ciliary pumping.
The asymmetric ciliary-type mode of beating was found to be accompanied by two distinct highly-symmetric beating modes. Most prevalent in the Sp-â parameter space was that of low-amplitude beating, with the pinned filament retaining low curvature throughout and generating the maximal time-averaged total normal force into the boundary amongst beating morphologies. Highly-buckled modes corresponding to less-rigid filaments produced minimal time-averaged total drag over each period of the background flow, with the high-order buckling modes still resolved by our piecewise-linear formulation of the governing elastohydrodynamics. The near-symmetry of each of these modes resulted in little time-averaged total drag in the direction parallel to the boundary. In contrast, entering the region of parameter space corresponding to ciliary beating is seen to strongly correlate with a sharp increase in time-averaged total parallel drag, intuitively the result of kinematic asymmetry arising from the distinct power and recovery strokes of ciliary beating. Hence, elastohydrodynamically-induced morphological changes in general have a dramatic effect on filament drag mechanics. Furthermore, the efficiency of the presented methodology would facilitate in-depth mechanical studies of primary cilia with regards to transmitted force and mechanotransduction, noting that the interaction of bending modulus variations, effective boundary conditions at the ciliary base, changes in cross section of the cilium, and cilium length are reported to be under-explored even in the restricted context of renal physiology (Nag & Resnick 2017).
In summary, we have considered in detail the validity of traditional and corrected resistive force theories for filaments in the presence of a planar boundary, concluding that, whilst in general these local drag theories may not be relied upon to perform well against non-local solutions for curved filaments, for filament motion parallel to a boundary the use of a resistive force theory may be remarkably accurate with appropriately-chosen coefficients, and viable over a range of boundary separations. Verified against previous highaccuracy methods, we have additionally presented an efficient non-local formulation of the governing elastohydrodynamics in a half-space, exemplifying its flexibility by application to a cantilevered filament in oscillatory flow. This study of passive filaments revealed a surprising asymmetric beating morphology found typically in active cilia, highlighting a complex relationship between passive elastic fibers, internally-forced filaments, the flows that they drive and the forces that they induce. , here for Sp = 3.13, a = 2π. In the upper panel we see that all choices oft0 lead to convergence to the same behaviour, with even the approximate initial symmetry of thet0 = π/2 instance (dashed) collapsing onto the common limiting behaviour. A posteriori validation of convergence to periodic motion from all phases is provided by consideration of the per-cycle change inxmax, as shown in the lower panel and highlighting eventual periodicity. Symmetric counterparts fort0 ∈ [π, 2π] have been omitted. (c) Signed curvature of the filament during the first periods of the background flow, sampled at 100 frames per period, fort0 = 3π/8, showing the transition to ciliary beating from a straight initial configuration. Colour online.
Appendix A. Influence of flow phase on cantilevered filament dynamics One might expect the phaset 0 of the background oscillatory flow of section 3.4 to have significant impact on long-time filament dynamics, the most notable of such effects being the selection of left-right polarity in resulting behaviours. We examine the longtime dynamics for a range of phases, flow amplitudes and elastohydrodynamic numbers, tracking in particular the maximally-attained displacement of the filament tip from the centrelinex = 0 throughout each period of the background flow, denotedx max and exemplified by figure 9(a). In particular, convergence ofx max to a constant would serve to validate, a posteriori, the assumption of periodic motion. Sampling Sp andâ each from the range [0.5, 12], as expected we observe a duality in solutions inherited from the left-right symmetry of the problem, and we additionally note the convergence of all solutions to a single periodic motion, where the period is aligned with that of the background flow to working precision. This is demonstrated by figure 9(b), which shows the evolution ofx max over a number of cycles of the background flow for Sp = 3.13 andâ = 2π. We see convergence of all solutions to a common limit cycle of periodic behaviour, with even the approximately-symmetric transient motion oft 0 = π/2 (shown dashed) eventually collapsing onto the same mode of behaviour. Hence we resolve that exploration of the Sp-â parameter space, holdingt 0 constant, is sufficient to capture the long-time behaviours of the cantilevered filament system.