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.
Human sperm propel themselves by propagating a travelling wave along a single, active flagellum; this motility is essential for migration through the female reproductive tract and natural fertilization. Recent work with microfluidic devices [1,2] has suggested the ability to direct and sort cells through their own motility, a potentially valuable advance in assisted reproduction therapy and in the livestock industry. Cell scattering at simple geometric features, such as the outside of a corner, appear to be dependent on viscosity and temperature; developing mechanical models to understand, interpret and optimize these effects for their exploitation is therefore of considerable interest. We will develop a mathematical model of a cell interacting with its environment, and its computational implementation, and study the dynamics of a realistic model sperm swimming over a backstep feature to study the effect of elastic, viscous and geometric parameters. The model will combine geometric nonlinearity of the elastic flagellum with non-local hydrodynamic interactions, and will be solved numerically via an implicit finite difference method for the elastohydrodynamic equations, combined with a hybrid slender body theory/boundary integral method for the low Reynolds number fluid dynamics.
The motor apparatus driving the flagellar waveform is a remarkably phylogenetically conserved structure known as the axoneme. The axoneme in human sperm comprises nine doublet microtubules, linked to each other and a central pair by passive elastic structures, with additional stiffening from outer dense fibres and the fibrous sheath (for a recent review focused on mechanically relevant features, see ). Motor proteins bound to the microtubules exert forces on adjacent doublets in a coordinated manner to induce bending moments along the length of the flagellum, causing bending, which is in turn resisted by the surrounding fluid. The fluid mediates interactions with surrounding surfaces and other cells; the flagellar waveform emerges from this nonlinear coupling.
Machin  showed that in order to generate experimentally observed waveforms the flagellum must actively bend along its length, and developed a linearized theory that has formed the basis of many subsequent studies. The theory that bending is produced by relative sliding of internal microtubules was subsequently proposed by Satir , and the sliding mechanism was modelled in early studies by Brokaw [6,7], using the formalism of an active internal moment per unit length in an elastic filament. The regulation of the active motor proteins that cause this sliding, and their oscillatory behaviour, is however a subject of continuing enquiry [8–10], with modelling playing an important role in comparing regulatory theories . A number of studies since the 1970s have provided significant insights into how potential mechanisms of dynein regulation can produce the types of bending waves observed in nature (e.g. [8,12–15]).
The importance of large-amplitude elastohydrodynamic flagellar modelling was established by Gadêlha et al. , who delineated the range of validity of small-amplitude elastic theory and showed that for sufficiently high viscosity relative to flagellar stiffness, a buckling instability can give rise to waveform asymmetry without domain boundaries or asymmetric internal actuation. The numerical implementation of Gadêlha et al.'s study built on a model of passive flexible fibres in shear flow , although replacing the non-local hydrodynamics of the latter with a local drag-velocity law. The combination of three-dimensional, time-dependent flow, with the hydrodynamic interactions arising from fixed and moving boundaries, with active filament mechanics is computationally demanding; the majority of sperm models until the last decade made similar approximations for the fluid dynamics, or small-amplitude linearization of the flagellar wave.
Liron, Gueron and colleagues (e.g. [18,19]) modelled cilia arrays, taking both non-local fluid dynamics and geometric nonlinearity into account, building on earlier work by for example Lighthill  and Hines & Blum . However this formalism, expressed in terms of bending angles rather than flagellar position, does not appear to have been generalized to a free-swimming cell with the associated boundary condition resulting from the presence of a head. More recent work using the finite-element and finite-volume methods and cluster computing has also been focused on cilia ; another successful recent approach is the regularized stokeslet method combined with a generalized immersed boundary method .
While the fluid dynamic interaction of sperm with plane boundaries has received significant attention since the work of Rothschild over 50 years ago , motivating a number of experimental and theoretical studies [24–27], the interaction of sperm with ‘non-trivial’ geometric obstacles involving angles and curves or complex interfaces is a subject of growing recent interest [28–31].
Denissenko et al.  showed how sperm scatter at a range of angles when encountering the outside of a corner in an artificial microchannel maze, and that the scattering angle is modulated by viscosity; Kantsler et al.  studied the effect of very close interactions of sperm and the biflagellate algae Chlamydomonas with these features. The geometric nature of the female reproductive tract is also highly convoluted, further motivating the need for models which can accommodate complex wall shapes. These studies suggest tantalizing opportunities to direct and sort motile sperm on passive 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.  with the Lighthill–Gueron–Liron (LGL) theorem  for non-local slender body theory and the boundary element  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.
3. 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.
At microscopic scales, fluid dynamics can be modelled by the incompressible Stokes flow equations 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)=Xt 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’), 3.2 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 (2.2). It will also be convenient to make use of the regularized stokeslet Sϵ of Cortez , which corresponds to a spatially smoothed force; a frequently used implementation in three dimensions  takes the form 3.3 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 , an extension of the work of Lighthill , 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 , 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, 3.4 Here and in what follows, 0≤s≤L is an arclength parametrization for the flagellum, and fvis is the viscous force per unit length exerted by the fluid on the flagellum. The coefficients ξ∥ and ξ⊥ are parallel and perpendicular resistance coefficients similar to those of Gray & Hancock  and take the form 3.5 the parameter q being a length scale chosen intermediate in magnitude between b and L. The symbols , and 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 (2.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. , we will rewrite these terms in another commonly used form, , 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  enable the representation of these surfaces by boundary integrals; rigidity of the surfaces enables the use of single layer boundary integral representations , 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 , 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 3.6 Similar equations, but with the first two terms replaced by a single slender body integral , 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.
The elastohydrodynamic formulation we will work with was derived by Tornberg & Shelley , and extended to an internally driven flagellum by Gadêlha et al. ; 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 . The internal elastic contact force Fint and moment Mint exerted on the proximal flagellum [0,s0) by the distal flagellum (s0,L), respectively, are given by 3.7 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 (s0,s0+δs) and taking the limit as yields 3.8 Noting that , the local term of equation (2.6) can then be written as 3.9 For brevity, we will write the non-local (integral) velocities from equation (2.6) as V (written out explicitly in the appendix, equation (.2)). Non-dimensionalizing with scales L for position, 1/ω for time, ωL for velocity and E/L2 for tension and moment density yields the following dimensionless elastohydrodynamic equation: 3.10 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. ) 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. ) by the inclusion of the V term on the left-hand side.
Similar to Gadêlha et al. , the inextensibility constraint ∂t(Xs⋅Xs)=0 can be used with the elastohydrodynamic equation (2.10) to deduce an ordinary differential equation which must be satisfied by the line tension T, 3.11 The above equation is derived via the identity 3Xss⋅Xsss+Xs⋅Xssss=0 and its derivative with respect to s. As previously [17,16], we introduce the term λSp4(1−Xs⋅Xs) to the left-hand side of equation (2.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 (2.10) and (2.11). The assumption of zero contact force and moment at the distal (s=1) tip of the flagellum combined with the elasticity equations (2.7) yield (in dimensionless variables) 3.12 Taking the dot product of the first equation with Xs, using the identity Xs⋅Xsss=−Xss⋅Xss 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 FH and MH and non-dimensionalize them with the elastic scalings E/L2 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 FH and MH, respectively. With the appropriate scalings, the proximal boundary conditions are then 3.13 where . From these equations, we also derive the tension condition at the proximal end, . The calculation of the quantities FH and MH 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 UH and ΩH of the head; while UH 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.  used travelling waves of internal moment, which calculations from experiment  confirm are a good model. We therefore specify in dimensionless units, .
3.3 Numerical implementation
The elastohydrodynamic equation (2.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 . 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 3.14 variables with tildes denoting that values from the previous iteration are taken.
The non-local hydrodynamic term V in equation (2.10) is approximated by forming the slender body/boundary integral problem of determining fvis, ϕH and ϕW using the most recent approximations to and 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 3.15 Each iteration requires the solution of a linear system for the unknown discrete values of X(sl,tn+1), T(sl,tn+1), UH and ΩH, where l=0,…,Ns denotes the spatial grid coordinate and n=0,1,… the time step. We found that Ns=160 and 200 time steps per beat were sufficient to yield accurate results. The discrete form of equations (2.10) and (2.11) provide 4(Ns+1)+6=650 linear equations, the additional six equations arising from the translational and angular velocity of the cell head. The nonlinear correction is then a system of 3(Ns+Nh+Nb) linear equations, where Nh and Nb are the number of elements on the head and domain boundary, respectively.
To implement the boundary conditions (2.12) and (2.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 3.16 where ΔFH, ΔMH are corrections for the effect of the flagellum. The calculations of 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 non-local hydrodynamic correction V and the non-local corrections to the force and moment balance ΔFH, ΔΩH require the solution of a slender body theory–boundary integral hydrodynamic problem. Calculation of the grand resistance matrix 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 . 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.
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. , we consider planar waveform actuation, which is appropriate for cells swimming through high viscosity fluids such as cervical mucus . The semi-axes of the ellipsoidal head, modelled with the boundary element method, are ax=0.05L, ay=0.03L, az=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 Nm2, beat frequency 10–20 Hz giving an angular frequency ω≈100 rad s−1 . Taking a flagellar radius of 0.5 μm, viscosity μ≈0.14 Pa⋅s (similar to mucus analogue ) 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 m0=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 , 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 ); 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 (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 ). 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 ) 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 the effect of varying sperm number over finer increments for backstep height zero (a,b) and h=0.2L (c,d), with results summarized in figure 4a. Simulations were performed for Sp=13,13.5,…,17 over both a strip geometry and a backstep of height h=0.2L, so that the sperm cells initially start 0.2L above the surface, and then increase to around 0.4L after the backstep. In figure 3a–d, colour is matched to increasing sperm number, so that light green corresponds to Sp=13 and red to Sp=17. 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 ) 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  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.
Figure 5b plots asymmetry versus arclength for sperm numbers in the range 13–17, the effect being largest at higher sperm number. The asymmetry at the tip of the flagellum for a strip versus no boundary is shown in figure 5c as a function of sperm number.
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 non-local 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 . 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. . 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–31,45]. Viscous interactions of course become increasingly important in high viscosity fluids such as mucus and laboratory analogues. Kantsler et al.  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 .
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  and are relevant to energy-producing bioreactors . 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.
Trajectory and waveform data supporting the findings in this paper have been submitted to Dryad; see reference .
T.D.M.J. designed the research, implemented algorithms, analysed results and co-wrote the manuscript. H.G. designed the research and contributed to the writing of the manuscript. D.J.S. designed the research, helped with algorithmic implementation, co-wrote the manuscript and supervised the project.
This work was supported by EPSRC grant EP/K007637/1 to D.J.S., supporting T.D.M.J. T.D.M.J. is supported by a Royal Commission for the Exhibition of 1851 Research Fellowship. H.G. acknowledges a Hooke Fellowship held at the University of Oxford. The computations described in this paper were performed using the University of Birmingham's BlueBEAR HPC service, which was purchased through SRIF-3 funds. See http://www.bear.bham.ac.uk for more details.
Conflict of interests
The author declare no competing interests.
The authors thank Dr Eamonn Gaffney (University of Oxford), Prof. Ray Goldstein (University of Cambridge), Dr Vasily Kantsler, Dr Petr Denissenko (University of Warwick), Prof. John Blake (University of Birmingham) and Dr Jackson Kirkman-Brown (University of Birmingham/Birmingham Women's NHS Foundation Trust) for continuing valuable discussions.
Appendix A. Calculation of hydrodynamic terms
Following non-dimensionalization, the hydrodynamic model yields the following equation for the dimensionless fluid velocity away from the flagellum: A 1 The non-local contribution to the velocity V on the slender body is similarly given by A 2
At each step of the iterative solution to the nonlinear problem, the collocation code solves the integral equation A 3 for the unknown hydrodynamic force per unit length fvis and unknown stresses ϕH, ϕW.
The collocation code discretizes the flagellum with 160 elements, with the non-local contribution to the LGL slender body theory computed by the midpoint rule with constant force per unit length over each element. The force per unit area on the ellipsoidal head of 32 mesh elements is calculated using routines from BEMLIB  with 20 point Gauss–Legendre quadrature as described in detail in the appendix. The wall boundary is discretized into elements of width 0.075L, using regularized stokeslets with ϵ=0.01L. Integration is performed with repeated Gauss–Legendre quadrature with 4×4 points per element for the near-singular wall integrals, and a 2×2 point rule elsewhere.
To implement the boundary conditions (2.12) and (2.13), Gadêlha et al.  approximated the force and moment on the head by a grand resistance matrix  multiplying the velocity and angular velocity. In dimensional variables, the grand resistance matrix expresses the force and moment on a moving rigid body as The blocks and are 3×6 matrices yielding the force and moment terms, respectively. For example, a sphere of radius a in the absence of hydrodynamic interactions would have dimensionless grand resistance matrix given by This approach is convenient because the linearity of the relationship means that the head velocity UH and angular velocity ΩH can be dealt within the implicit formulation as unknowns in the linear algebra problem. To generalize to a non-local hydrodynamic model taking into account the effect of the flagellum and nearby boundary, the force and moment will be decomposed as consisting of part which is linear in velocity and angular velocity via resistance matrices and a remaining contribution from the flagellum. The matrices and are determined via the boundary integral method, taking into account the potentially highly significant effect of the wall feature, but not the subleading effect of the flagellum, which is accounted for as a correction, as described below.
Elastic scalings are used to non-dimensionalize all forces and moments, i.e. E/L2, E/L for FH and MH, respectively, with E/L3 for force per unit length fvis and E/L4 for stress ϕH, ϕW. The additional corrections ΔFH and ΔMH referred to in equation (2.16) are determined as part of the iterative process by performing a slender body/boundary integral calculation of , and with the most recent approximation to available, yielding in dimensionless variables A 5 where is the head centroid. Using also the most recent iterates for and , the corrections are then given by A 6 These corrections appear on the right-hand side of the linear system.
- Received November 25, 2014.
- Accepted February 17, 2015.
© 2015 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.