Meshfree and efficient modelling of swimming cells

Locomotion in Stokes flow is an intensively-studied problem because it describes important biological phenomena such as the motility of many species' sperm, bacteria, algae and protozoa. Numerical computations can be challenging, particularly in three dimensions, due to the presence of moving boundaries and complex geometries; methods which combine ease-of-implementation and computational efficiency are therefore needed. A recently-proposed method to discretise the regularised Stokeslet boundary integral equation without the need for a connected 'mesh' is applied to the inertialess locomotion problem in Stokes flow. The mathematical formulation and key aspects of the computational implementation in MATLAB/GNU Octave are described, followed by numerical experiments with biflagellate algae and multiple uniflagellate sperm swimming between no-slip surfaces, for which both swimming trajectories and flow fields are calculated. These computational experiments required minutes of time on modest hardware; an extensible implementation is provided in a github repository. The nearest neighbour discretisation dramatically improves convergence and robustness, a key challenge in extending the regularised Stokeslet method to complicated, three dimensional, biological fluid problems.


I. INTRODUCTION
Inertialess locomotion in Stokes flow describes the motility of many types of sperm, bacteria, algae and protozoa. This topic has received extensive attention from mathematical modellers, starting with the classic work of Taylor [1] and continuing to the present day [2][3][4].
From the early work into the swimming of sea-urchin spermatozoa [5], to investigations into the orientation of biflagellates in shear flows [6], there has been a lot of interest into modelling biological swimmers. This interest has been extended recently towards understanding and developing novel microswimmers. Topical examples of these involve studies into the microscale flow dynamics of ribbons and sheets [7], and the modelling of self-propelling toroidal swimmers based on the hypotheses of Taylor and Purcell [8], as well as the study of phoretic toroidal swimmers [9]. Such works have the potential to enable the use of targeted drug delivery, amongst other things, through being able to guide microswimmers through complex biological environments [10], and improve diagnostics and management of male infertility by analysis of imaging data.
Of particular recent interest is the collective behaviours of microswimmers. The differences in these behaviours appear to have significant biological implications, an example of which is the collective swimming of bovine sperm in the presence of viscoelasticity, behaviour which is not apparent in a purely viscous fluid [11]. Other species of sperm exhibit collective behaviours which impact both swimming and the ability to effectively fertilise the egg, some species of opossum sperm are often seen swimming as a cooperative pair [12]. In addition to collective behaviours, the effects of interactions with other particles and/or boundaries have been recently shown to create interesting dynamics [13][14][15][16].
While each of the models presented above are in some sense idealised, the ability to further reduce detailed swimmer models to simplified representations provides the opportunity for extracting significant scientific information which may not be accessible otherwise. Such models allow for creation of a coarse-grained representation of a swimmer [4] reducing complex behaviour into a set of swimming 'modes' and their associated limit cycles. Detailed fluid dynamic modelling can also allow for calculation of parameters for continuum models [17] and give understanding of hidden aspects of swimmers' characteristics such as energy transport along a flagellum [18] or internal moment generation [19], as well as providing insight into the exact mechanisms for the collective swimming behaviours mentioned above.
Numerical methods are generally required to model finite amplitude motions, wall effects and swimmer-swimmer interactions. A range of numerical approaches exist, with perhaps the most extensively-studied being those based on singular, or regularized singular solutions of the Stokes flow equations, specifically resistive force theory [5], slender body theory [20], boundary integral methods [21], and regularized Stokeslet methods [22][23][24][25]. These techniques remove the need to mesh the volume of fluid, requiring only the solution of integral equations formulated on the surface of the swimming body/bodies and lines such as cilia and flagella, reducing both the cost of meshing/remeshing a continually moving domain, and the number of degrees of freedom of the resulting linear system. Other techniques can be used to perform computational analysis of swimmers, such as the use of the force coupling method to investigate the dynamics of suspensions of up to 1000 swimmers [26], and the immersed boundary method [28] for understanding the role of fluid elastic stress on flagellar swimming [29].
As reviewed recently [30], regularized Stokeslet methods have the further major advantage of removing the need to evaluate weakly-singular surface integrals, and enabling slender bodies such as cilia and flage To improve on the computational efficiency of the regularized Stokeslet method while retaining most of its simplicity of implementation, a method was proposed by Smith [30], involving taking a coarser discretisation for the unknown traction than that used for numerical quadrature of the kernel, enabled by the use of nearest-neighbour discretisation.
The method proved significantly more accurate for significantly lower computational cost, potentially enabling more complex and realistic problems to be investigated with given computational resources.
In this article we generalise the nearest-neighbour discretisation three dimensional regularized Stokeslet method to inertialess locomotion, in particular focusing on uniflagellate pushers modelling human sperm and a model biflagellate. In section II we will briefly review the mathematical definition of the inertialess free-swimming problem in the boundary integral formulation. In section III we implement the nearest-neighbour discretisation of the free-swimming problem with a single swimmer in an unbounded fluid. We then formulate the task of tracking the trajectory of the cell as an initial-value problem. The discretisation is then generalised in section IV to incorporate rigid boundaries and multiple swimming cells. Finally, in section V we present the results of numerical experiments with uniflagellate and biflagellate swimmers, and in section VI we discuss the method and further practical applications. Key aspects of the implementation in Matlab R /GNU Octave are given, and a github repository is provided with the full code necessary to generate the results in the report, as well as templates for applying the method to novel problems in very low Reynolds number locomotion.

II. THE FREE-SWIMMING PROBLEM
The dynamics of a Newtonian fluid at very low Reynolds numbers, associated with locomotion of cells, is described by the Stokes flow equations. The dimensionless form of the equations is, augmented with the no-slip, no-penetration boundary condition u(X) =Ẋ for boundary points X, where overdot denotes time-derivative. We note here that, for the kinematicdriven problems in the present paper, the viscosity term has been non-dimensionalised out of the PDE; for a force-driven problem the viscosity term would appear in the dimensionless group (the sperm number). Initially we will consider a single swimmer in a three dimensional unbounded fluid which is stationary at infinity. Two classical problems in Stokes flow are the resistance problem -which involves calculating the force and moment on a rigid body made to translate and rotate in stationary fluid, and the mobility problem -which involves calculating the rigid body motion due to an imposed force and moment.
The free-swimming problem in Stokes flow is a variant of the mobility problem. Rather than -or perhaps in addition to -the body being driven by imposed forces, it translates and rotates as a result of changing its shape. In this section we will briefly review this problem, which has been solved numerically in many previous studies, and introduce our notation.
As usual for the regularized Stokeslet method, the fluid velocity u j at location x (suppressing time-dependence) is approximated by a surface integral over the surface ∂D of the swimmer, The regularisation error associated with equation (2) has been discussed previously [22] and will not be reviewed here. In this paper we will treat the approximation as exact. The surface of the body will undergo motions that may be described by a model formulated in a body frame -for example a frame in which the head of the cell does not move. If the body frame coordinates are ξ, and the body frame is described by the matrix of basis vectors (equivalently a rotation matrix) B = (b (1) |b (2) |b (3) ) and origin x 0 then the laboratory frame coordinates and velocities are, x =ẋ 0 +Ḃ · ξ + B ·ξ.
Denoting the rigid body velocity and angular velocity of the frame by U and Ω respective, we then have, Applying the condition u(x) =ẋ on ∂D in equation (2) yields the regularized Stokeslet boundary integral equation, where it is understood that repeated indices (such as j in the above) are summed over, and unrepeated indices (such as i in the above) range over {1, 2, 3}.
If at time t, the body frame origin x 0 and orientation B are known, and a model is given for the swimmer shape ξ and motionξ in the body frame, then then unknowns of the problem are the surface traction f (X) for X ∈ ∂D, the translational velocity U and angular velocity Ω. The problem is closed by augmenting equation (7) with the force and moment balance equations; here we assume that the inertia and moment of inertia of the swimmer are negligible. The full problem is then given by, where ijk is the Levi-Civita symbol.
Numerical discretisation of the problem (8) will in general involve N vector degrees of freedom for the traction f , three unknowns for the components of the translational velocity U and three unknowns for the components of the angular velocity Ω, totalling 3N + 6 scalar unknowns in total. Through numerical collocation, problem (8) can be formulated as 3N + 6 linear equations. In the next section we will describe a nearest-neighbour regularized Stokeslet discretisation of this problem.

A. A single swimmer in an unbounded fluid
The discretisation of the regularized Stokeslet method is discussed in detail in [30]; in brief we suggest that a good balance of ease-of-implementation and numerical efficiency can be achieved by discretising the integrals via a quadrature rule, with the key modification of using a finer discretisation for the rapidly-varying regularised Stokeslet and a coarser discretisation for the more slowly-varying traction. A simple way to achieve this is through nearest-neighbour interpolation of the traction. The resulting method contains the original and extensively-used method of Cortez and colleagues [22] as the limiting case in which the discretisations are equal.
Replacing the integrals in problem (8) with numerical quadrature yields the discrete problem, where dS(X[q]) denotes the quadrature weight associated with the local surface metric.
A subtlety here concerns the calculation of the nearest-neighbour matrix when dealing with time-evolving geometries, and in particular the case when different bodies approach closely. As ). Applying this interpolation to problem (9) yields, Computationally, problem (11) corresponds to 3N + 3 + 3 linear equations in 3N + 3 + 3 scalar unknowns (F j [n] := g j [n]dS(x[n]) for n = 1, . . . , N , followed by U j and Ω j ). These equations can be expressed in block form as, . . .
where the blocks have entries given by, for m, n = 1, . . . , N, and the velocity U and angular velocity Ω are expressed as 3 × 1 column vectors.

B. Computing swimmer trajectories via an initial-value problem
The position and orientation of a swimmer can be expressed as a position vector and a frame of basis vectors b (j) . (2) so it is sufficient to formulate the problem in terms of two basis vectors only, or six scalar degrees of freedom. Of course, this formulation still contains redundant information -three Euler angles constrain precisely the body frame, however the basis vector approach is very straightforward to implement.
Noting thatẋ we may then formulate the calculation of trajectories as a system of 9 ordinary differential equations, where evaluation of the functions U ( solving the swimming problem (8), for example via the discretisation (11). The 'outer' problem (14) can be solved using built-in functions such as ode45 in Matlab R or lsode in GNU Octave.
For practical purposes, when using a built-in initial value problem solver such as ode45, the tractions f i (X), required to compute the rate of energy dissipation and the flow field, may not be automatically available. To record this information, we may introduce the Augmenting the swimming problem (14) with equations (15) then yields an approximation to the force distribution available external to ode45 by numerically differentiating H i (x, t) with respect to time.

A. Boundaries and fixed obstacles
Mammalian sperm usually migrate and fertilise within a thin film of viscous fluid between opposed surfaces, and are typically imaged between a microscope slide and coverslip.
Indeed, the major effect of boundaries on microswimmer flow fields has long been recognised [31]. Therefore it is important to take boundary effects into account in fluid dynamic simulations. The 'Blakelet' and its regularized counterpart found by Ainley and colleagues [32] (see also recent work by Cortez) is an elegant and efficient way to model a single infinite plane boundary; certain other geometrically simple situations possess similar fundamental solutions. However, it is important for full generality to take into account more complex boundary, and perhaps also fixed obstacles, present in the flow.
Representing the boundary by B, the swimming problem becomes, The blocks have entries given by, n] for m, n = 1, . . . , N, where the total number of force unknowns is , and the right hand sides are given by,

B. Multiple swimmers
The last situation we will consider is where there are multiple swimmers -which are not necessarily discretised by equal size sets -as well as a boundary. The numerical discretisation is somewhat more complicated, and so we modify our notation in an attempt to make the implementation more interpretable. Suppose that we now have N sw swimmers, described by collocation points with ith components x . The discretisation will follow the ordering convention, which is inherited by the right hand side velocities and the force discretisation. If the number of force points associated with swimmer r is N s (r), and the number of force points associated with the boundary is N b , then the number of vector force unknowns is N f = Nsw r=1 N s (r)+N b . The size of x is then 3N f and the total number of scalar degrees of freedom in the system is 3N f + 6N sw . We will define the index ι(r) to be the location of the rth swimmer in the x i vector, with ι(1) = 1 and ι(r) = r−1 γ=1 N s (γ) for 1 < r N sw .
The quadrature points may be denoted X [1], . . . , X[Q] as previously; the Stokeslet matrix is then constructed as, To construct the remaining blocks, we introduce the notation 1 (n) to be the column vector of length n with every entry equal to 1 and 0 (m×n) to be the m × n matrix of zeros. We also define the N f × N sw matrices, Then, . . .
with ⊗ denoting the Kronecker product.
Recalling that ν[·, ·] denotes the nearest-neighbour matrix, we define the N s (r)×1 column vectors, and the N sw × N f matrices, Then the 3N sw × 3N f blocks A F and A M are, Finally, denoting the orientation matrix of the rth swimmer by B

(r)
ij and its body frame waveform as ξ (r) j , the terms of the right hand side take the form, where V (r) Now that we have defined A S ij , A U , A Ω , A F , A M and V i , the 3(N f + 2) × 3(N f + 2) linear system is of the form given by equation (17).

V. RESULTS AND ANALYSIS
We now turn our attention to the application of this method to two model problems: (1) a single biflagellate swimming in an infinite fluid, and (2) multiple sperm cells swimming between two boundaries. The implementation for both these model problems is provided in the associated github repository. After presenting the results for these swimming problems we will discuss the convergence of the method for the two types of swimmer provided and compare with the results obtained through the classic Nyström discretisation (when the force and quadrature discretisations are the same).

A. Biflagellate in an infinite fluid
We will first apply the algorithm in section III A to model a biflagellate, superficially similar to various marine algae, swimming in an unbounded fluid. We model the beat pattern of the cell (figure 1a) following Sartori et al [33], writing the flagellar tangent angle ψ in the form where s and t are dimensionless arclength along the flagellum and time respectively. We find that choosing ψ 0 (s) = −2.5s, ψ 1 = 0.7 + 0.15 sin (2πs) , φ (s) = −2πs, 0 ≤ s ≤ 1, provides a sufficiently representative test case for the computational algorithm. Of course a more realistic beat for a genuine biflagellate species such as Chlamydomonas reinhardtii could be appended as required.  C. Convergence of the method with discretisation refinement A practical refinement heuristic for assessing the convergence (with increased discretisation) of the nearest neighbour method is given by Smith [30]. For testing the convergence of the present swimming problems we denote the maximum discretisation spacings from [30]  3. Generate a more refined head quadrature discretisation by halving h H q and repeat step 2.
4. Generate a more refined head force discretisation by halving h H f and repeat step 2.

Repeat steps 3 and 4 until a suitable level of convergence is reached.
mer, and indeed for each component of a single swimmer (the head and flagellum may be discretised differently for example). To this end we apply the existing convergence heuristic in stages as outlined in table I. To measure the convergence we compare the straight line distance travelled over one full beat of the swimmer's flagellum. In contrast to the classical (Nyström) discretisation [22], there is no tight coupling between the regularisation parameter and the discretisation length scales [30]. As a consequence of this we allow the choice of regularisation parameter to be guided by the geometry of the swimmer (chosen here to be related to the dimensions of the flagellum).
We have analysed the convergence of the results for the following cases: a single swimming biflagellate (as described in §V A), a single swimming sperm (as in §V B) with no boundary, and a single swimming sperm with boundary. We also assess the effect of the boundary through fixing the sperm discretisation and applying the heuristic of table I to the     requires many more degrees of freedom to reach the same levels of convergence. This convergence rate could be improved in the Nyström case by varying (as discussed in [22]), however as previously discussed the nearest neighbour discretisation is much more robust to this parameter.

VI. DISCUSSION
This report has described an extension of the nearest-neighbour regularized Stokeslet method [30] to enable the simulation of multiple force-and moment-free cells swimming in a bounded domain. Cell trajectory calculations were achieved by casting the task as an initialvalue problem; by integrating the force at each step it was additionally possible to store the evolving force distribution to enable post-calculation of the velocity field. The method was assessed on two problems of a type which may be of interest in the biological fluid mechanics community: swimming of a biflagellate in an unbounded domain, and motility of multiple human sperm between two no-slip surfaces.
Numerical experiments provide evidence that the method is relatively efficient and converges well, requiring minutes to solve the problems described above, without specialist computational hardware, and we note with interest the significantly improved convergence of this method when compared to the classic Nyström discretisation. While the construction of the matrices is somewhat tedious, the underlying concept of the method -a coarse/fine discretisation of the boundary integral equations to address the fact that the force distribu-tion varies more slowly than the kernel -should ensure that the method is comprehensible and extensible by non-specialists. Crucially, no true 'mesh' generation (i.e. with connectivity tables) is required to simulate a new swimmer of interest. We hope that these properties of ease-of-use, extensibility and efficiency make the method appealing to potential users, and in support of this aim we provide all Matlab R code used to generate this report in the repository github.com/djsmithbham/nearestStokesletSwimmers. Within this repository, a template file nnSwimmerTemplate.m is provided which sets out how new swimmers can be added to the existing codebase.
There are many potential extensions for this work spanning the whole field of locomotion at low Reynolds number. The convergence properties of this method mean that it may be valuable for high-throughput analysis of experimental data, or (perhaps with adaptations to deal efficiently with long-range interactions) suspensions of relatively large numbers of swimmers. It would be interesting to see if the modification of the method to take into account viscoelastic effects would allow for the collective swimming behaviour of sperm seen by Tung et al. [11] to be reproduced from an idealised model of swimming. There is potential for this method to be applied to the world of phoretic swimmers to examine the dynamics of many phoretic particles, or to the case of swimmers driven by magnetic fields. The computational efficiency of this method can also be exploited through modelling multiple swimmers in complex environments, for example ciliary flow. While such flows would previously have been simulated and then applied as a background flow to a swimmer, with this efficient method one would be able to model the ciliary beating patterns directly and could allow for a more realistic interaction between swimmers and their environment.