From flagellar undulations to collective motion: predicting the dynamics of sperm suspensions

Swimming cells and microorganisms are as diverse in their collective dynamics as they are in their individual shapes and propulsion mechanisms. Even for sperm cells, which have a stereotyped shape consisting of a cell body connected to a flexible flagellum, a wide range of collective dynamics is observed spanning from the formation of tightly packed groups to the display of larger-scale, turbulence-like motion. Using a detailed mathematical model that resolves flagellum dynamics, we perform simulations of sperm suspensions containing up to 1000 cells and explore the connection between individual and collective dynamics. We find that depending on the level of variation in individual dynamics from one swimmer to another, the sperm exhibit either a strong tendency to aggregate, or the suspension exhibits large-scale swirling. Hydrodynamic interactions govern the formation and evolution of both states. In addition, a quantitative analysis of the states reveals that the flows generated at the time scale of flagellum undulations contribute significantly to the overall energy in the surrounding fluid, highlighting the importance of resolving these flows.

where Λ is the tension and M = K Bt × dt/ds is the bending moment. The external force per unit length, f , and torque per unit length, τ , arise due to the viscous stresses along the flagellum and, in the case of multiple swimmers, steric interactions between neighboring swimmers. The torques τ D per unit length arise due to the time-dependent preferred curvature. At the free end (s = l), the boundary conditions are such that the tension and moment are zero, while at the attachment point to the cell head, we require a clamped-end condition be satisfied.
To solve these equations numerically, we first discretize the flagellum into N flag segments of length ∆L. The segments have positions Y n and orientationst n for n = 1, . . . , N flag . The orientations are the discrete representation of the centerline tangents at the segment positions. Considering the tension and bending moment at the midpoints between adjacent segments and applying central differencing, we obtain the discretized force and torque balances where M n+1/2 = (K B /∆L)t n ×t n+1 . Multiplying through by ∆L and introducing F n = f n ∆L and T n = τ n ∆L as the total applied force and torque, respectively, on segment n, we may write Eqs. [1] and [2] as where F C n = Λ n+1/2 − Λ n−1/2 , T B n = M n+1/2 − M n−1/2 , and T C n = (∆L/2)t n × (Λ n+1/2 + Λ n−1/2 ). The driving torques, T D n , arising from the preferred curvature are given by where s n = (n − 1/2)∆L for n = 1, . . . , N flag , and The tension Λ n+1/2 enforces flagellum inextensibility at the level of each segment. In the discrete setting, they are the Lagrange multipliers associated with the constraints, In the discretized system, the boundary conditions at the free end are satisfied by taking Λ N flag +1/2 = M N flag +1/2 = 0. The clamped-end condition on the cell head is recovered by treating the head as another segment of size 2a with the bending moment and inextensibility constraint computed with respect to the attachment point rather than the head center.
Until now, we have discussed the model in the context of a single swimmer. To allow for N swimmers, we must have N force and moment balances, one pair for each of the flagella. As the driving torques, tension, and bending moments are internal to each flagellum, the force and moment balances are coupled only through the external forces and torques. In our simulations, the coupling arises due to steric repulsion between the segments and heads and hydrodynamic interactions. Accordingly, the total external force on segment n is given by F n = F S n + F H n where F S n is the steric force and F H n is the hydrodynamic force on the segment. The only external torque on the segment is the viscous torque, T n = T H n The steric force F S n on segment n is given by the generic repulsive barrier force where the sum runs over I(n), the set of all segments/heads, m, with |Y n − Y m | < χR nm and χ = 1.1. The parameter R nm specifies the center-to-center distance for particles n and m at contact and χ sets the range over which the barrier force acts. If n and m are both segments, the contact distance is R nm = 2b, where as if n and m are both heads, we have R nm = 2a. Lastly, if n and m are a head-segment pair then R nm = a + b. The parameter F S sets the strength of the repulsion at contact, which in our simulations is

Force-coupling method
The force and torque balances, Eqs. (3) and (4) respectively, establish a low Reynolds number mobility problem whose solution is the coupled motion of the particles through the fluid. To solve this mobility problem, we utilize the force-coupling method (FCM) [7,6,5], which, for the sake of clarity, we describe here in the case of spherical particles of equal radius, however, it readily extends to spheroidal particles [5], such as our cell heads, and allows for polydispersity in particle sizes.
In FCM forces and torques on the particles are projected onto the fluid using a truncated and regularized force multipole expansion. This leads to the incompressible Stokes flow for pressure field, p, and fluid velocity field, u. The Stokes flow is generated by the forces and torques each particle n exerts on the fluid, which, in the case of our swimmers, will be F H n = −F C n − F S n , and The force and torques are projected to the fluid via the Gaussian functions, which for particle n are given by ] .
The motion of the particles is obtained by volume averaging the resulting fluid velocity using the same Gaussian functions. Specifically, the translational and angular velocities are given by If the particles have radius A, FCM reproduces the correct Stokes drag and viscous torque with σ ∆ = A/ √ π and σ Θ = A/(6 √ π) 1/3 . In our simulations, we preserve these ratios between the particle and the Gaussian envelope sizes. The segments are taken to have radius b = ∆L/2.2, while the semi major-axes of the spheroidal heads are a = 3b and the semi minor-axes are b. All simulation parameters are summarized in table 1 and the datasets provided have the same units.

Updating positions and orientations
As our simulations are quasi-2D in that the motion of the swimmers is restricted to a plane, we have Ω n = Ω nẑ for all n and can introduce angle θ n such thatt n = (cos θ n , sin θ n ). Therefore, after computing the V n and Ω n for each segment and head, the positions and orientations are obtained by integrating in time subject to the inextensibility constraints given by [5]. To do this numerically, we rely on a second-order implicit BDF scheme [1] and Broyden's method [2] to find the updated positions and orientations, as well as the Lagrange multipliers Λ n+1/2 .

Swimmer center of mass and director
We define the center of mass of the m-th swimmer as Here, Y (m) 0 is the position of the m-th swimmer's head, and Y (m) k is that of its k-the segment. The center of mass velocity is then given by The director of swimmer m is given byn

Force multipole expansion for a single swimmer
The singular force density for a single swimmer is given by where the asymmetric dipole is related to the torque through G ij = ϵ ijk T k /2. We have suppressed the superscript index, since we are considering only one swimmer. Expanding this force density about the swimmer's center of mass, X cm , and defining the position of the center of mass relative to the n-th segment/head as The first term is zero because there is no external force on the swimmer. The remaining two terms are the force dipole and quadrupole, respectively, which drive the flow observed in the limit |x − X cm | ≫ d. We compute the symmetric, traceless and anti-symmetric parts of these tensors, given by as functions of time. These are shown in Fig. 2 in the main text, with the exception of G A (t) which is zero for all time as there is no net-torque acting on the swimmers.

Stochastic variation of the undulation frequency
We introduce stochastic changes of the undulation frequencies to explore how variability within a population, as well as the fluctuations of individual frequencies over time can impact collective dynamics. After every fixed average period T = 2π/ω, a new undulation frequency, ω B , is randomly assigned to a swimmer to replace its previous value, ω A . New frequencies are drawn from a log-normal distribution with parameters (µ ln , σ ln ) given by where the mean frequency is ω and the standard deviation of the distribution is σ ω = δ ω ω. To enforce continuity in time for each swimmer's preferred curvature κ(s, t), the phase must also be updated from φ A to φ B . Specifically, at the time t AB when the frequency changes, we require This can be fulfilled by taking for each individual swimmer. We note that an alternative approach to introducing stochasticity is to assign a random frequency to each swimmer at the beginning of the simulation and keep it fixed for all time. This was adopted previously in [8]. We found that on time-scales that are relatively short compared to the evolution of the suspension this leads to qualitatively similar results as those obtained by allowing the frequency to vary over time. Over longer times, however, we expect, based on our monodisperse frequency simulations, that keeping the frequencies fixed could potentially lead to preferential clustering of swimmers with similar frequencies and hence structurally different dynamics.

Identifying and counting swimmer clusters
To effectively identify swimmer clusters, we must account for relative swimmer orientations in addition to relative distances. Based on this, swimmer i is taken to belong to a cluster if the modified distance ) with at least one other swimmer j in the cluster is d β (i, j) < d/3. The parameter β controls the influence of particle alignment and is chosen to be β/d = 1 6 1 tan(π/8) = ( √ 2 − 1)/6. With this choice of β, swimmers with orientations differing by an angle of π/4 must be half as far apart as aligned swimmers to still be considered clustered. Using this notion of distance, we employ the hierarchical cluster algorithm in Matlab for cluster identification. The number of clusters over time for the two simulations with 1000 individual swimmers are shown in Fig. 1. Individual swimmers are counted as clusters of size one. We see that the number of clusters for the fixed frequency simulation is much lower than that for the simulation where the undulation frequency is allowed to vary stochastically. Additionally, we see that when stochasticity is introduced, cluster growth is halted very early in the simulation.

Effect of film thickness on suspension dynamics
As the flows induced by singular force distributions decay more slowly in 2D than in 3D, we expect the hydrodynamic interactions to play a stronger role as L z → 0 and, as a result, the instability to exhibit a faster growth rate for thinner films. Figs. 2 and 3 show results from simulations containing 85 swimmers in domains of linear size L = 4.43d and for thicknesses L z = 0.277d, 0.554d, and 1.107d. The effective area fraction is ν = 1.08 and the swimmers are initially in a polar state. In accordance with the slower formation of  the bending waves as shown in Fig. 2, the slower decay (see Fig. 3) of the order parameters S 1 and S 2 (as defined in the main text) for thicker films indicate that the shorter range of the hydrodynamic interactions lead to slower growth rates of the instability. The smallest value of L z used in Fig. 3 matches that for the larger simulations presented in the main text. The order parameter decay rate in those simulations, however, is faster still due to the longer modes afforded by the larger system size.

Effect of undulation frequency variability on aggregation
As mentioned in the main text, low variability in the swimmers' undulation frequencies does not entirely suppress cluster formation and still allows for neighboring swimmers to synchronize and align their waveforms. Snapshots from simulations of a system of size L = 4.43d with 85 initially aligned swimmers and different values of σ ω reveal that for small σ ω , cluster formation persists (see Fig. 4). As σ ω increases, the number of visually identifiable clusters is reduced and we see instead the emergence of a bending wave.

Tuning the resistive force theory
To remove the hydrodynamic interactions between the swimmers yet still allow for propulsion via undulation, we solve the mobility problem using a drag-based resistive force theory (RFT) [4,3] rather than FCM. With the drag-based model, the velocity and angular velocity of segment n of a flagellum are related to the force and torque on the segment through where α ∥ , α ⊥ , and β are mobility coefficients for the segments. These mobility coefficients are initially estimated based on FCM simulations of straight filaments and subsequently adjusted to ensure that the swimmer waveform and free swimming velocity are comparable to those given by the simulations with FCM, see

Energy density spectrum of a randomly forced Stokesian fluid
In this section, we show that a k −3 behavior of the fluid energy spectrum can be related to spatially uncorrelated forcing of a 2D fluid. In the following, the Fourier transform and its inverse are defined aŝ We begin by considering a force density, f (x), restricted to a square domain, where Ω = [−L/2, L/2] × [−L/2, L/2] and Ω is the indicator function defined on the set Ω. The Fourier transform of the restricted force density is then The Fourier transform of the fluid velocity resulting from the restricted force density is given bŷ is the Fourier representation of the inverse Stokes operator, k = |k|, andk = k/k.
The fluid energy density spectrum associated withû L (k) is where θ is related to the two-dimensional wave-vector through k = (k cos θ, k sin θ) and ⟨·⟩ denotes the ensemble average. The energy density spectrum of the infinite system will then be given by S(k) = lim L→∞ S L (k). As ) , the fluid energy density spectrum is directly related to the forcing through For S L (k) ∝ k −3 to hold, we must have ⟩ independent of k. If, in addition, we require statistically independent and identical force components that do not depend on θ, we have that ⟨f where A is a constant independent of L. The L 2 dependence on the system size guarantees that S L (k) is well-defined and non-zero as L → ∞ and that the total energy will scale with the system size.
Assuming that the forcing is a statistically stationary spatial process, we may express the correlations of the force density as As the force density is strictly restricted to Ω, we may replace the integral over Ω with one over all of R 2 . Thus, which, from the convolution theorem, may be expressed as Finally, using the expression for ⟨f L (k)f ⊤ L (−k) ⟩ given above, we see that for any L. Allowing L → ∞ then yields ⟨f (x + y)f ⊤ (x)

⟩ = AIδ(y)
and we see that a k −3 decay in the fluid energy density spectrum is given by uncorrelated spatial forcing of the fluid.