Spermatozoa scattering by a microchannel feature: an elastohydrodynamic model

Sperm traverse their microenvironment through viscous fluid by propagating flagellar waves; the waveform emerges as a consequence of elastic structure, internal active moments and low Reynolds number fluid dynamics. Engineered microchannels have recently been proposed as a method of sorting and manipulating motile cells; the interaction of cells with these artificial environments therefore warrants investigation. A numerical method is presented for large-amplitude elastohydrodynamic interaction of active swimmers with domain features. This method is employed to examine hydrodynamic scattering by a model microchannel backstep feature. Scattering is shown to depend on backstep height and the relative strength of viscous and elastic forces in the flagellum. In a ‘high viscosity’ parameter regime corresponding to human sperm in cervical mucus analogue, this hydrodynamic contribution to scattering is comparable in magnitude to recent data on contact effects, being of the order of 5°–10°. Scattering can be positive or negative depending on the relative strength of viscous and elastic effects, emphasizing the importance of viscosity on the interaction of sperm with their microenvironment. The modulation of scattering angle by viscosity is associated with variations in flagellar asymmetry induced by the elastohydrodynamic interaction with the boundary feature.


Summary
Sperm traverse their microenvironment through viscous fluid by propagating flagellar waves; the waveform emerges as a consequence of elastic structure, internal active moments and low Reynolds number fluid dynamics. Engineered microchannels have recently been proposed as a method of sorting and manipulating motile cells; the interaction of cells with these artificial environments therefore warrants investigation. A numerical method is presented for largeamplitude elastohydrodynamic interaction of active swimmers with domain features. This method is employed to examine hydrodynamic scattering by a model microchannel backstep feature. Scattering is shown to depend on backstep height and the relative strength of viscous and elastic forces in the flagellum. In a 'high viscosity' parameter regime corresponding to human sperm in cervical mucus analogue, this hydrodynamic contribution to scattering is comparable in magnitude to recent data on contact effects, being of the order of 5 • -10 • . Scattering can be positive or negative depending on the relative strength of viscous and elastic effects, emphasizing the importance of viscosity on the interaction of sperm with their microenvironment. The modulation of scattering angle by viscosity is associated with variations in flagellar asymmetry induced by the elastohydrodynamic interaction with the boundary feature. microdevices, however a better understanding of the subtle nonlinear physics of how flagellated swimmers interact with geometric features must be developed; to aid with this understanding we will develop a mathematical and computational approach which accounts for elasticity, viscosity and their interaction, without the need for large-scale computational resources. To this end, we will bring together the active elastic formulation of Gadêlha et al. [16] with the Lighthill-Gueron-Liron (LGL) theorem [18] for non-local slender body theory and the boundary element [32] and regularized stokeslet methods [33,34] to capture the influence of a non-trivial nearby surface. We will use this approach to explore how sperm scatter near geometric features due to elastohydrodynamic interaction over hundreds of flagellar beats with a single computer core, and quantify how the balance of viscosity and elasticity modulates this effect via changes to the flagellar waveform.

Mathematical model
The mathematical model of a sperm interacting with a geometric feature will be derived from (i) the Stokes flow equations, with a non-local hydrodynamic model, and (ii) geometrically nonlinear elasticity for an internally actuated flagellum. We will first derive the equations for the two parts of problem, before describing (iii) the numerical implementation.

Hydrodynamics
At microscopic scales, fluid dynamics can be modelled by the incompressible Stokes flow equations 0 = −Vp + μ∇ 2 u, V · u = 0, (3.1) where u is velocity, p is pressure and μ is dynamic viscosity. For our problem, these equations will be augmented with the no-slip, no-penetration condition u(X) = X t for points X on the solid boundary, where subscript t denotes time derivative.
The linearity of the Stokes flow equations enables the construction of solutions to satisfy boundary conditions via discrete and/or continuous sums of suitably-weighted fundamental solutions. These techniques replace solid surfaces, such as the sperm flagellum, head and its surrounding microenvironment, by line or surface distributions of immersed forces. A concentrated point force located at y with strength F produces a velocity field (the 'stokeslet'), the symbol δ jk being the Kronecker delta tensor and the summation convention being used. The symbol S(x, y) will be used to denote the second rank tensor in equation (3.2). It will also be convenient to make use of the regularized stokeslet S of Cortez [35], which corresponds to a spatially smoothed force; a frequently used implementation in three dimensions [33] takes the form The parameter > 0 defines the length scale over which the point force is smoothed; this smoothness property is particularly convenient for the formulation of boundary integral methods. The LGL theorem [19], an extension of the work of Lighthill [20], derives from a line distribution of singular stokeslets and source dipoles: an approximate expression for the flow field at the surface of a moving slender body, accurate to O( b/L), where b is the radius and L is the flagellar length. Ignoring image systems, which are not required in our formulation, and using the properties of the stokeslet to reorder the source and field points, we have the expression for the approximate velocity field produced by the slender body v, S(X(s 0 , t), X(s, t)) · f vis (s, t) ds. (3.4) Here and in what follows, 0 ≤ s ≤ L is an arclength parametrization for the flagellum, and f vis is the viscous force per unit length exerted by the fluid on the flagellum.  [36] and take the form the parameter q being a length scale chosen intermediate in magnitude between b and L. The symbolŝ s,n andb are unit tangent, normal and binormal. Whereas Gueron and Liron [18,19] considered the dynamics of a cilium projecting from a plane boundary, and hence the associated image systems, in this study we will not require these terms because surfaces will be represented via boundary integrals. Equation (3.4) can be considered a non-local extension of resistive force theories, which retain only the first three terms. To couple LGL to the elastohydrodynamic model of Gadêlha et al. [16], we will rewrite these terms in another commonly used form, −(1/ξ ⊥ )(I + (γ − 1)ŝŝ) · f vis , with γ = ξ ⊥ /ξ playing a similar role to the drag anisotropy ratio of resistive force theory, but depending on the choice of q. The precise value of q is not critical provided that b q L because changes to the resistance coefficients are accompanied by changes to the integrals; for our study with b = 0.01L, we choose q = 0.1L, leading to γ ≈ 1.4.
To model a sperm, we will consider a cell with a rigid head as well as a flagellum, swimming near a rigid step-like surface. The linearity of the Stokes flow equations means that a solution satisfying the additional no-slip boundary conditions associated with the head and the wall may be constructed by linear superposition. Moreover, the Lorentz reciprocal relation and its regularized analogue [33] enable the representation of these surfaces by boundary integrals; rigidity of the surfaces enables the use of single layer boundary integral representations [37, p. 32]. In this study, we will use a hybrid approach, representing the head via a surface distribution of singular stokeslets with stress φ H , discretized via BEMLIB [32], and the wall by regularized stokeslets and boundary elements, with stress φ W [33,34]. The full fluid dynamic model for the velocity field on the surface of the flagellum is therefore Similar equations, but with the first two terms replaced by a single slender body integral − L 0 S · f vis ds, hold on the surface of the head and the wall. In the next section, we will discuss the equations of an internally driven elastic flagellum, and their coupling to the fluid mechanics.

Elastohydrodynamics
The elastohydrodynamic formulation we will work with was derived by Tornberg & Shelley [17], and extended to an internally driven flagellum by Gadêlha et al. [16]; the central feature of this approach is to formulate the problem in terms of the flagellar position X(s, t) and line tension T(s, t). Alternative approaches based on bending angles and curvatures [38,39] have also been pursued, as has complex curvature [40]. The internal elastic contact force F int and moment M int exerted on the proximal flagellum [0, s 0 ) by the distal flagellum (s 0 , L), respectively, are given by where E is constant elastic modulus and m(s, t) is a prescribed active moment density representing the internal flagellar motors. Balancing elastic and viscous forces acting on a segment of flagellum (s 0 , s 0 + δs) and taking the limit as δs → 0 yields Noting thatŝ = X s , the local term of equation (3.6) can then be written as For brevity, we will write the non-local (integral) velocities from equation ( time, ωL for velocity and E/L 2 for tension and moment density yields the following dimensionless elastohydrodynamic equation: The parameter Sp = L(ξ ⊥ ω/E) 1/4 is the sperm number, which quantifies the relative importance of viscous and elastic effects. This model can be seen as an extension of linear models (such as Camalet et al. [14]) by the inclusion of the nonlinear terms on the right-hand side, and an extension of hydrodynamically local models (such as Gadêlha et al. [16]) by the inclusion of the V term on the left-hand side. Similar to Gadêlha et al. [16], the inextensibility constraint ∂ t (X s · X s ) = 0 can be used with the elastohydrodynamic equation (3.10) to deduce an ordinary differential equation which must be satisfied by the line tension T, The above equation is derived via the identity 3X ss · X sss + X s · X ssss = 0 and its derivative with respect to s. As previously [16,17], we introduce the term λSp 4 (1 − X s · X s ) to the left-hand side of equation (3.11) to dampen numerical errors in flagellar length. The value used in this study is λ = 80, though as found by Gadêlha et al. the solution is insensitive to the precise value of λ. The final part of the mathematical model is the specification of the boundary conditions for equations (3.10) and (3.11). The assumption of zero contact force and moment at the distal (s = 1) tip of the flagellum combined with the elasticity equations (3.7) yield (in dimensionless variables) 0 = −X sss + mn + TX s and 0 = X ss at s = 1. (3.12) Taking the dot product of the first equation with X s , using the identity X s · X sss = −X ss · X ss and the second equation yields the distal tension boundary condition, T = 0. At the proximal end of the flagellum, the boundary conditions are given by considering the force and moment exerted by the fluid on the head. We denote these quantities F H and M H and non-dimensionalize them with the elastic scalings E/L 2 and E/L, respectively. In the inertialess Stokes flow regime, the total force and moment acting on the head are zero, so by Newton's third law, the force and moment on the flagellum at s = 0 are also given by F H and M H , respectively. With the appropriate scalings, the proximal boundary conditions are then F H = X sss − mn − TX s and M H ∧ X s = −X ss + Mn, at s = 0, (3.13) where M = 1 0 m ds. From these equations, we also derive the tension condition at the proximal end, F H · X s = −X ss · X ss − T. The calculation of the quantities F H and M H with non-local hydrodynamic interaction is described in more detail in the next section and the appendix. Finally, we introduce the translational and angular velocity U H and Ω H of the head; while U H and two components of the angular velocity are constrained by knowledge of the function X, there is an independent rotational component of the motion that defines the principal bending plane of the flagellum. These quantities will be determined by kinematic considerations and the implementation of the boundary conditions.
To complete the mathematical model, it is necessary to specify the internal active moment m(s, t). Gadêlha et al. [16] used travelling waves of internal moment, which calculations from experiment [3] confirm are a good model. We therefore specify in dimensionless units, m(s, t) = m 0 cos(ks − t).

Numerical implementation
The elastohydrodynamic equation (3.10) is treated with a Crank-Nicolson-type finite difference discretization, with the second-order central differences in the interior, and third-order one-sided difference for the boundary conditions, using coefficients taken from Fornberg [41]. The higher order boundary stencil produced comparable errors to the central stencil on polynomial test functions. Both linear and nonlinear terms are treated implicitly; nonlinearity of these equations is dealt with by performing an iterative process on every time step, with the operator on the left-hand side at t + dt being linearized as variables with tildes denoting that values from the previous iteration are taken. The non-local hydrodynamic term V in equation (3.10) is approximated by forming the slender body/boundary integral problem of determining f vis , φ H and φ W using the most recent approximations toX andX t available; details are given in the appendix.
At the first iteration of each time step, the converged values from the previous time step are used as starting guesses for all variables, except for X which is approximated via linear extrapolation. The nonlinear iteration is terminated when the maximum difference in position between successive iterations relative to the distance travelled by the flagellum over the time step falls below 0.5%. Similarly, the auxiliary equation for the tension at t + dt is linearized as  .13), the force and moment on the head are a priori unknown and need to be determined as part of the coupled problem. The force and moment are decomposed into a linear part, given by the grand resistance matrix associated with rigid body motion in the vicinity of the wall, and an additional subleading correction resulting from the influence of the flagellum. Following non-dimensionalization with the elasticity scalings, the force and moment on the head may then be expressed as where F H , M H are corrections for the effect of the flagellum. The calculations of R and the corrections are described in the appendix. In summary, each time step requires a number of iterations to solve the nonlinear problem and each iteration involves the solution of a sparse linear system arising from the finite difference discretization of the elastohydrodynamic equations. The 'right-hand side' terms arising from the nonlocal hydrodynamic correction V and the non-local corrections to the force and moment balance F H , Ω H require the solution of a slender body theory-boundary integral hydrodynamic problem. Calculation of the grand resistance matrix R requires the separate solution of a boundary integral problem with multiple right-hand sides to determine the force and moment resistances associated with the rigid body modes of the head and the wall interaction. The code is implemented in Fortran 90 (gfortran, GNU Compiler Collection); linear systems are equilibriated and solved by LU factorization with the LAPACK routines dgeequ and dgesv, respectively, and the boundary integrals over the sperm head are calculated with routines from BEMLIB [32]. A typical run of 200 beats with 500 boundary elements required approximately 24 h walltime on a single core of a 2.2 GHz Intel Sandy Bridge E5-2660 node.

Results
The numerical scheme is applied to predict the trajectory of a sperm-like cell swimming in an unbounded fluid at varying Sp, over a 'backstep' (the latter being shown in figure 1a), the limiting case of zero backstep height being referred to as a 'strip'. As in Gadêlha et al. [16], we consider planar waveform actuation, which is appropriate for cells swimming through high viscosity fluids such as cervical mucus [42]. The semi-axes of the ellipsoidal head, modelled with the boundary element method, are a x = 0.05L, a y = 0.03L, a z = 0.04L, which correspond to a 5 × 3 × 4 µm head for a flagellum of length L = 50 µm. The swimmer is initially at rest, with a straight flagellum, and a 'soft start' is applied whereby the internal shear moment is initially low and smoothly increases to its maximum, reaching 99% after five beats. The sperm number of a human gamete can be approximated by using bending stiffness E ≈ 5 × 10 −21 Nm 2 , beat frequency 10-20 Hz giving an angular frequency ω ≈ 100 rad s −1 [3]. Taking  0.5 µm, viscosity μ ≈ 0.14 Pa · s (similar to mucus analogue [42]) yields the normal resistance coefficient ξ ⊥ ≈ 0.503 and sperm number Sp ≈ 15.8. Therefore, we will consider a range of sperm numbers between 13 and 17, fixing the magnitude of the internally generated shear-force m 0 = 240 and wavenumber k = 6π . The resulting waveforms are shown in figure 1a. As sperm number increases, beat amplitude is suppressed, as is observed for sperm in high viscosity medium [42], leading to a reduction in side-to-side yaw. All simulations in infinite fluid, i.e. with no nearby boundaries, produced trajectories which were straight overall, once the within-beat yaw was accounted for (data submitted to Dryad repository [43]); flagellar waveforms for Sp = 13, 15 and 17 are shown in figure 1b,c. Figure 2 shows a planar projection of the trajectories (X(0, t), Y(0, t)) and the tangent angle θ := arctan(dY/dX(s = 0)) (in degrees) of those trajectories, of cells swimming over backsteps of varying height. The derivative dY/dX is calculated numerically by sampling the trajectory at the temporal midpoint of each beat-cycle and taking centred differences. Colour indicates the trajectory over the backstep of the height in figure 2a,c,e with green denoting h = 0 and red denoting h = 0.5. Simulations were performed over backsteps of height h = 0.05, 0.1, . . . , 0.5 and are displayed up to the time at which X(0, t) ≥ 1.
The results in figure 2a,c,e suggest that the backstep affects swimmers at different sperm numbers differently, producing a range of scattering angles. However, it is important in these results to factor out the effects of the strip from the backstep. Taking the (lightest) green trajectory, representing a strip, as a baseline comparison, it is evident that for all sperm numbers the hydrodynamic effect of the backstep is to deflect the swimmer downwards relative to a strip trajectory. Figure 2b,d,f reveals that this downward deflection is not smooth, rather there is a sharp bump at x = 0 where the head initially passes over the backstep, and a further bump at around x = 0.3 where the effect of the step itself becomes subleading relative to boundary interactions between the head and the lower wall.
Simulations were also performed comparing the effect of the backstep to a 'cliff' geometry, with the lower portion of the backstep missing (data submitted [43]). After passing the backstep, cells swam straight as though in an infinite fluid, suggesting that the majority of the angular deflection occurs due to interaction with the lower boundary; boundary forces change suddenly over a step jump, and the cell acts as though it were above a higher boundary. Additionally, simulations over a strip at Sp = 13 for different starting heights (data submitted [43]) showed that attraction to the surface initially increased and then decreased as height above the surface increased, which suggests that hydrodynamic boundary attraction is responsible for the behaviour in figure 2a,b. Figure 3 shows  Figure 3a,b shows for a sperm swimming over a strip, the boundary repels the swimmer more at this close distance as sperm number is increased. This effect is to be expected because increasing the sperm number increases the relative strength of viscous to elastic forces, thus the effect of the boundary is likely to be enhanced as Sp increases. The initial dip in figure 3b is an artefact of the numerical soft start of our system, as the waveform emerges from a straight initial state. Figure 3c,d shows a larger range of scattering angles than for fixed sperm number over various backstep heights, of the order of 10 • . Furthermore, additional simulations (data submitted [43]) showed that this hydrodynamic deflection was not sensitive to the phase of the waveform as it passed over the backstep, in contrast to scattering due to contact forces (R. Goldstein 2014, personal communication). Figure 4a shows the effects of changing sperm number, giving the deflection for a strip, a backstep, and  their difference. A slight increase in the magnitude of this difference is observed as sperm number is increased, owing to increased hydrodynamic interaction mediated by viscosity. Figure 4b summarizes the effect of varying both backstep height and sperm number simultaneously, quantified by the 'final deflection angle' θ d , i.e. the value of θ for which X = L. At Sp = 13 deflection is always negative, whereas for Sp = 15, 17 deflection is always positive. The relationship between θ d and h is non-monotone at the lower sperm number but is monotonic in the higher range. At Sp = 13, the  deflection angle initially increases in magnitude, then decreases after the maximum at around h = 0.15L. This riser height corresponds to a distance of 0.35L between the cell and the boundary, which is where boundary attraction is strongest at this sperm number. For Sp = 15, 17, the deflection angle decreases monotonically with backstep height in the range we have considered. This effect probably occurs because at these sperm numbers the strip causes the cell to pitch away. However in all cases, increasing backstep height to 0.5L results in a plateau.
The effects of the backstep on the waveform are summarized in figure 5, which show the waveform shape with and without the boundary, and quantitative measures of the asymmetry of the waveform. Recall that the flagellar actuation is symmetric; waveform asymmetry is produced due to increased hydrodynamic drag arising from proximity to the wall [44] affecting closer portions of the flagellum more than further portions. Figure 5a shows waveforms at sperm number Sp = 13, 17 in infinite fluid as well as over a strip. In infinite fluid, the waveform is symmetrical for all sperm numbers considered, while the presence of a boundary gives rise to a waveform asymmetry that increases with Sp.
'Asymmetry' is quantified by sampling the flagellar wave every 41 numerical time steps (relative to a beat cycle of 200 time steps), projecting into the body frame and calculating the average lateral position relative to the body frame centreline over a fixed period, in this case beats 82-90. This quantity is plotted as a function of arclength in figure 5b; its distal (s = 1) value is plotted in figure 5c.

Discussion
A numerical method for simulating the swimming of monoflagellate cells over geometric features was presented and applied to model sperm interacting with microchannel backstep feature. The scheme incorporates non-local hydrodynamics with large-amplitude active filament mechanics. We believe this method to be the simplest generalization of previous work that is capable of taking into account nonlocal hydrodynamic interaction geometrical features. The linearity of the Stokes flow equations entails that the largest error in our method arises from the LGL slender body theory, which is at worst on the order of the square root of the slenderness ratio. Accuracy of the method of regularized stokeslets is on the order of the regularization parameter near the boundary, and its square far from the boundary where the swimmer is located. Future work may consider boundary integral modelling of the flagellum also; however, we do not expect that this would qualitatively change swimmer trajectories.
The interaction between the cell and the lower boundary involves the competing effects of asymmetric hydrodynamic forces leading to waveform asymmetry and boundary repulsion, and the pitching behaviour associated with swimmer/boundary interaction [26]. At lower sperm number and at greater distances from the boundary, waveform asymmetry is smaller, and the cell pitches towards the boundary. At higher sperm number and closer distances from the boundary, waveform asymmetry is larger and the cell pitches away. The effect of the backstep is a sudden drop in the lower boundary, which changes the relative importance of these effects; waveform asymmetry is reduced relative to hydrodynamic attraction, and the net result is a deflection towards the lower boundary after the backstep relative to the expected trajectory over a strip (figure 2).
Analysing sperm scattering over a backstep, we found that hydrodynamic effects may be comparable in magnitude in the relatively high viscosity range considered to the contact interactions found experimentally by Kantsler et al. [2]. A transition is predicted from scattering towards the backstep at lower viscosity to scattering away from the backstep at higher viscosity. Qualitatively this behaviour is similar to the temperature-related transition in Kantsler et al.'s observations (with lower temperature corresponding to higher viscosity); the correspondence is not exact however, with Kantsler et al.'s observations being carried out with bull sperm in low viscosity buffer, and with cells exhibiting very close interaction with the boundary, compared with our longer range interactions and sperm number representative of human cells in mucus analogue that we chose to focus on in this study. Clearly integrating both surface interactions and hydrodynamics will be necessary to develop a comprehensive model, particularly at higher sperm number/viscosity.
The role of hydrodynamic interactions in determining surface attraction and more complex effects associated with boundary features continues to receive significant theoretical attention and is stimulating novel mathematical approaches [28][29][30][31]45]. Viscous interactions of course become increasingly important in high viscosity fluids such as mucus and laboratory analogues. Kantsler et al. [2] noted the need to take both elastic and steric interactions into account; modelling very short length scale or contact interactions, with either glass, epithelium, cumulus or even ciliated surfaces, and their effect on the flagellar wave, is a topic of importance, though numerical simulation requires taking account of the rapidly varying hydrodynamic force and electrostatic interactions as the swimmer approaches these boundaries. We hope that the numerically implicit method, potentially also combined with adaptive refinement of the boundary element meshes, will enable accurately resolved simulation of sperm-like swimmers in very near surface-contact in future work. Other valuable methods for modelling three-dimensional sperm motility and elastic-fluid interaction include models based exclusively on regularized stokeslets [46,47] and techniques such as stochastic rotation dynamics [48].
While we have used our model to examine a swimmer representative of human sperm, the approach is applicable to a much wider range of eukaryotic cells, including the sperm of other species and, with a slight reworking of the head boundary condition, biflagellate organisms such as the green alga Chlamydomonas. These species are of particular interest as they have been used as models for flagellar synchronization [49] and are relevant to energy-producing bioreactors [50]. For these systems, the model may also be extended to include a non-local hydrodynamic contribution from other swimmers. Larger swimming organisms, such as Caenorhabditis elegans, have also been shown to be significantly affected by interactions with a structured microenvironment [51,52]. Another application is the design and optimization of biomimetic artificial microswimmers (e.g. [53,54]). Because the model includes internal periodic actuation via prescribed bending moments, it might be used to optimize actuation for various purposes such as forward progress, subject to constraints such as fixed mechanical energy. Furthermore, the inclusion of geometrical boundary features and the use of sperm number allow such optimization to be tailored to specific environments. The elastohydrodynamic model can additionally be used to solve the inverse problem of estimating internal moments from observed flagellar data, potentially allowing us to examine how nature has optimized swimming in various environments and informing truly biomimetic design.
Despite the linearity of the Stokes flow equations, the interaction of sperm with their microenvironment presents a subtle nonlinear mechanics problem. Sperm scattering depends nonlinearly on the ratio between viscous and elastic forces, with even a simple backstep feature producing attractive or repulsive scattering of cells depending on parameter values. These scattering effects may be valuable in sorting cells in microdevices, in addition to giving insight into the complexity of how sperm interact with their microenvironment. The combination of mechanical models and experiment will provide the best way to understand and exploit these effects for biomedical applications.
Data accessibility. Trajectory and waveform data supporting the findings in this paper have been submitted to Dryad; see reference [43].