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

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

Simon F. Schoeller and Eric E. Keaveny1Department of Mathematics
Imperial College London
South Kensington Campus
London SW7 2AZ, UK

Abstract

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.

Swimming cells and microorganisms encompass the entire range of cell types and exhibit a great variety of cell geometries and swimming strategies used to move in viscosity dominated environments [1, 2]. Some organisms, such as Volvox colonies, can be nearly spherical and utilize for propulsion thousands of relatively short flexible flagella distributed along their surface. Others, like sperm cells, have a single, long flagellum that propagates bending waves, leading to a large time-periodic shape change of the cell. Though stereotyped, the propagated waveforms, frequencies, and head shapes of sperm cells vary appreciably across different species [3].

There is not only great variation amongst individual cells, but also in the patterns that they form and in how populations are organized. This can range from algal plume formation in bioconvection [4, 5] to turbulence-like swirling in bacterial baths [6, 7, 8, 9]. Great variation is exhibited by sperm suspensions of different species that have been observed to form vortices near planar boundaries [10], aggregate into coherent sperm trains [11], or exhibit turbulence-like swirling as first documented nearly \num70 years ago by Lord Rothschild [12, 13, 14, 15]. What is not clear, however, is how the individual differences in sperm cells, such as geometry, flagellum length, waveform and flexibility, give rise to the differences in the way sperm populations organize themselves. Mathematical modeling and simulation provide a route to explore the connection between individual and collective dynamics where experiments might otherwise be very challenging.

In the past \numrange1015 years, there has been much work on the development of mathematical models [8, 16, 17, 18, 19, 20, 21, 22, 23, 24] of microorganism suspensions, especially to understand the turbulence-like state found in bacterial baths. To deal with large numbers of cells, the models rely on a reduced description of the swimmers, often treating them as simple, rigid objects (e.g., rods, dumbbells, spheres, ellipsoids) that interact through steady, dipolar flows and steric repulsion. These models have been effective in reproducing the large-scale motion of the population and relating its formation to the sign of the coefficient of dipolar flow induced by individuals.

While such models now broadly shape our thinking about suspensions of swimming microorganisms, they reduce the diversity of the microscopic world into a handful of parameters and a limited number of degrees of freedom. They ignore the complexity and time-dependence of cell shapes, as well as the time-dependent and beyond-dipolar features of the flows generated by the shape changes. Indeed, many of these details have been included in models of individual or small collections of swimmers. In the case of sperm cells, the addition of flagellar motion leads to hydrodynamically-induced attraction, synchronization, and phase-locking of planar, flagellar waveforms [25, 26, 27]. In simulations of \numrange50100 swimmers [28, 29] these effects were found to induce aggregation and clustering of cells. This tendency to aggregate may aid sperm in forming coherent groups such as sperm trains, however, it remains unclear how a turbulence-like state could then be reached.

In this paper, we perform detailed simulations of up to \num1000 interacting sperm, resolving the coupled flagellar dynamics along with the flows generated at the undulation time-scale and at sub-flagellum length scales. We report that variation in the individual dynamics across the population, here controlled by the undulation frequency, can suppress the tendency to aggregate and instead lead to a density dependent turbulence-like state with hydrodynamics still playing a key, but now different role. Additionally, an analysis of the clustered and turbulence-like states reveals the strong influence of flagellar undulation on the quantities used to explore large-scale motion in swimmer suspensions. These results suggest that only minor variations in sperm behavior across species are needed to produce very distinct collective dynamics.

Mathematical model

Our swimmer model is closely related to several models [30, 31, 32, 33, 34] employed to describe the dynamics of sperm cells in a viscous fluid and extends directly from a model previously used to capture undulatory locomotion through structured media [35]. We briefly summarize it here in the context of a single swimmer and provide a more detailed description in the SI.

Our swimmers consist of two elements, the flagellum and the cell head. The flagellum is treated as an inextensible, yet flexible beam of length l and bending modulus KB, while an oblate spheroid of semi-major and minor axes a and b, respectively, represents the cell head. The overall swimmer length is taken to be d=2.1a+l, where the distance 0.1a accounts for a linkage between the flagellum and the head. The flagellum is parametrized by arclength s∈[0,l] such that the position of a point along the flagellum is Y(s) and the unit tangent at that point is ^t(s)=dY/ds. The flagellum is driven internally by the moments per unit length, τD, that arise due to the preferred curvature,

κ(s,t)

=K0sin(ks−ωt+φ)⋅{2(l−s)/l,s>l/21,s≤l/2,

(1)

where ω=2π/T is the undulation frequency, k is the wavenumber, φ is the phase, and K0 is the amplitude. The linear decay in the preferred curvature amplitude near the free-end is introduced to better reproduce observed flagellar waveforms. Allowing the flagellum to also be subject to external applied forces, f, and torques, τ, per unit length due to viscous stresses, the force and moment balances along the flagellum are

dΛds+f

=0

(2)

dMds+τD+^t×Λ+τ

=0.

(3)

where Λ is the tension and M=KB^t×d^t/ds is the bending moment.
At one end (s=0), the flagellum is attached via a clamped-end condition to the cell head, while the other end (s=l) remains free.

The continuous beam equations, (2) and (3), are discretized to obtain force and torque balances on Nflag segments of the flagellum, as well as the cell head (see SI for details). For the n-th segment of the flagellum, the balances are

FCn+FHn

=0,

(4)

TBn+TCn+TDn+THn

=0,

(5)

where TBn are torques arising due to bending, FCn and TCn are constraint forces and torques associated with tension, TDn are the torques due to the preferred curvature, and FHn and THn are the hydrodynamic forces and torques on the segment. The force and torque balances yield a low Reynolds number mobility problem for the translational and angular motion of the segments that is identical to that for a collection of rigid particles. We solve this mobility problem using a regularized multipole approach for Stokes flows known as the force-coupling method (FCM) [36, 37, 38, 35, 39]. In the limit of large Nflag, FCM provides an approximation of the hydrodynamics consistent with a regularized version of slender body theory [40, 41, 42, 43, 44, 31, 33]. After solving the mobility problem for the velocities and angular velocities of the flagellum segments and cell head, the differential equations for their positions and orientations are integrated in time using a second-order backwards differentiation scheme. Broyden’s method is then used to solve the resulting system of equations for the updated positions and orientations, and the Lagrange multipliers associated with the forces and torques arising from the inextensibility constraint.

We choose the parameters in our simulations to reproduce a cell geometry, flagellar waveform, and swimming speed close to that reported for ram sperm [45, 14] and in the range of other mammalian sperm [46, 47]. The head size a is such that a/l=0.0470 and b=a/3. The wavenumber is k=2π/l, and the curvature amplitude is K0=12.76/l. Finally, we set the dimensionless parameter known as the Sperm number, which provides a measure of the ratio of the viscous to elastic forces [48, 49, 2, 35, 50], to Sp=(4πωη/KB)1/4l=12.0, where η is the dynamic viscosity. With these parameters, 23.7 undulation periods are required for the sperm to swim its flagellum length.

We perform three dimensional simulations in domains of size L×L×Lz corresponding to 2D periodic thin films with finite thickness Lz=0.277d. The corresponding fluid flow is resolved on a grid of size Nx×Ny×Nz. Free surface conditions at z=0 and z=Lz are established by restricting the motion of the swimmers to the mid-plane z=Lz/2 and applying periodic boundary conditions in all three directions. This particular choice in boundary conditions and film thickness corresponds to those commonly employed in experiments on swimming microorganisms [6, 16, 51, 52].

When simulating multiple swimmers, a short-ranged repulsive force between nearly touching segments is included to capture steric repulsion between swimmers (see SI). The strength of this force at contact is determined by the parameter FS. This and all simulation parameters, including their values in simulation units, are summarized in Table 1.

Parameter

Description

Value in simulation units

η

Viscosity

1

Sp

Sperm number

12

K0

Curvature amplitude

0.2

KB

Bending modulus

1800

b

Segment radius and head height

1

a

Head half axis

3

d

Swimmer length

70.1

FS/d

Steric reference force

15π

Nx,y

Grid dimensions (in-plane)

{2048,3072}

Nz

Grid dimension (normal)

64

L

Linear domain size

931.3

Lz

Domain height

19.40

Nflag

Number of flagellum segments per swimmer

29

N

Number of swimmers

{100,…,1000}

Nt/T

Time-steps per undulation

300

σω/ω

Relative std. dev. of frequencies

{0,0.2}

Table 1: Simulation parameters \includegraphics

[width=]panel_streamlines_lowres.png

Figure 1: Fluid flow induced by an individual swimmer. A Instantaneous streamlines (grey) and fluid velocity field (red arrows) at z=Lz/2 produced by a single swimmer at one instant in time. B Period-averaged streamlines and fluid velocity field in a co-moving frame at z=Lz/2. Two snapshots of the swimmer are also shown.

Time-dependent and higher-order features of the swimmer flow field

Examining the flow field produced by a single swimmer and the force-moments associated with this flow, we see clear differences between their instantaneous and period-averaged values. The flow at the mid-plane of the thin film and the swimmer at an instant in time are shown in Fig. 1A, while Fig. 1B shows the flow averaged over one undulation period. These flow fields are comparable to those given by other computational studies of individual swimmers propelled by an elastic filament [30, 28, 29, 53, 27, 54]. While the period-averaged flow closely resembles the dipolar flow generated by a so-called pusher, we find that the instantaneous flow field at any point in time is markedly different, containing features of higher-order forces singularities.

\includegraphics

[width=]panel_moments.pdf

Figure 2: The symmetric, trace-less dipole, Gij, and quadrupole, KS,Aijk, coefficients as a function of time. The values are reported in the swimmer frame where the origin is the swimmer’s center of mass and ^e1 is the average swimming direction. Components that are zero at all times are excluded. Dipole entries are in units of F0d, while the quadrupole entries have units of F0d2.

To analyze this further, we compute (see SI and [49]) the in-plane dipole and quadrupole force moments for the swimmer as functions of time. The symmetric, trace-less dipole tensor G, the symmetric (in the first two indices), trace-less quadrupole tensor KS, and the anti-symmetric quadrupole tensor KA, are plotted over two undulation periods in Fig. 2. The time periodic nature of the force moments, especially their change in sign over the course of a period, are consistent with experimental measurements [51, 54] on flagellate microswimmers. The anti-symmetric dipole moment is always zero, as there is no net-torque on the swimmer. From these time series, we can extract the non-vanishing time-averages, the largest of which are the dipole coefficients ⟨(GS11,GS22)⟩=(−1.050,0.432)×10−3F0d, and the quadrupole coefficients ⟨(KS111,KS221)⟩=(−2.321,0.822)×10−3F0d2, where F0 is the in-plane drag on the cell head moving at the average swimming velocity. We note that for both the dipole and quadrupole these average values are much smaller than the maximum values attained during an undulation period. Quadrupole components (KS111 and KS221) are indeed expected to contribute substantially to the fluid flow around spermatozoa [55, 33, 33] and have also been linked to swimmer velocity correlations in bacterial suspensions [56]. The maximum value of the dipole coefficient in our model is more than a factor of 5.5 greater than its mean value, while for the quadrupole this increases to a factor of 7. Additionally, the relatively large values of the quadrupole coefficients produce strong flows that, although they decay faster than the dipolar flow, are non-negligible for significant distances from the swimmer. We estimate that the instantaneous quadrupolar flow dominates the dipolar flow up to separations of approximately 3d.

From this, a picture emerges that not only are the time-dependent aspects of the dipolar flow much stronger than the average flow, but that higher-order contributions are non-negligible when time dependence is included. Both of these features can impact the swimmer-swimmer hydrodynamic interactions, especially in semi-dilute to dense suspensions. As we will see in the proceeding sections, these features, linked to the dynamics at the level of individual cells, do indeed propagate to larger scales and can significantly affect the observed coherent structures and the processes by which they form.

\includegraphics

[width=]SpermSuspensionFig.png

Figure 3: Evolution of suspensions of 1000 swimmers. A Snapshots of a system with fixed undulation frequency taken after 30 (left) and 198 (right) periods. B Snapshots after 340 undulation periods from two stochastically-varying frequency RFT simulations starting from isotropic (left) and polar (right) initial configurations. C Snapshots of a simulation with randomly varying frequencies starting from a polar configuration at 60 (left) and 198 (right) undulations.

Suspensions with a monodisperse distribution of undulation frequencies form clusters

Using the methods described in the previous section, we simulate the dynamics of a suspension of N=1000 swimmers with Nflag=29 for \num200 undulation periods. In this simulation, all swimmers are driven by a prescribed preferred curvature, Eq. 1, with the same undulation frequency, ω=KB(Sp/l)4/(4πη), and wave-number k=2π/l. Each swimmer has a different phase, φ, drawn randomly from a uniform distribution. The swimmers are initially distributed uniformly in the domain and are aligned to swim in the −x-direction. Swimmer-swimmer hydrodynamic interactions are incorporated by solving the full mobility problem that couples the motion of all swimmer segments and heads, while steric interactions are included through short-ranged, pairwise, repulsive forces between all flagellum segments and cell heads. The effective area fraction is ν=Nd2/(4L2)=1.42. Based on the film thickness, the effective volume fraction, as used in [22], is approximately \num2.56. For suspensions of pusher dipoles at this effective volume fraction, one would expect rapid decay of the initial polar order followed by the onset of large-scale, fluctuating motion [22, 19].

Fig. 3A shows the suspension after \num30 and \num198 undulation. periods. The corresponding video showing the complete evolution is provided as ESM. While the suspension quickly evolves from the initial polar state, the dynamics are very different from the long-wavelength bending instabilities predicted by both simulations of rigid, steady pushers [18, 22] and kinetic or continuum models [17, 19, 20, 21, 24]. Instead, we find that the swimmers have a strong tendency to aggregate and form clusters within which flagella are aligned. A quantitative analysis of this cluster growth is included in the SI.
Once small clusters have formed, they remain intact and can attract other clusters in their vicinity. As two attracting clusters approach each other, they either rotate to swim in the same direction and merge to form a bigger, synchronized cluster, or instead swim in opposing directions and move apart at an enhanced relative velocity. Earlier simulations of sperm [28, 29] showed similar aggregation dynamics in two dimensional suspensions, but here we see that aggregation also dominates in larger, denser suspensions where long-range hydrodynamics could be expected to produce large-scale swirls and jets.

\includegraphics

[width=]panel_attraction.pdf

Figure 4: Relative center of mass displacements (in units of swimmer length) over four periods of undulation for two parallel swimmers with relative phase shifts Δφ=0,π/2 and π. For clarity, the cell heads attached to the left ends of the flagella are omitted.

The aggregation process is driven by the hydrodynamic attraction and phase-locking behavior that has been found previously with models of undulatory swimmers [26, 28, 29, 57, 34, 27]. We measure this for our system by examining in detail the pairwise interactions between two swimmers moving in the −x-direction with phase differences Δφ=0,π/2, and π. The vector fields in Fig. 4 show the displacement over four undulation periods of one swimmer centered at (x,y) relative to another placed at the origin. When the swimmers are far apart, all three vector fields resemble the displacement field associated with pusher dipoles. At smaller separations, however, the fields differ and the influence of the phase difference can be readily seen. We see that there is clear attraction, if the swimmers are in phase. When the phase difference is π/2, the swimmers attract and simultaneously move relative to each other as to align the crests of the flagellar waves. If they are out of phase, the swimmers again move relative to one another attempting to align wave crests, however, the relative shift is much larger. The large displacements observed near the swimmer’s ends are due to the steric forces between the head of one swimmer and the flagellum of the other.

Suspensions with stochastically varying frequencies reach a turbulence-like state

In the previous section, we saw that flagellum undulations can produce flows that lead to hydrodynamically-induced aggregation and cluster formation. While this provides a mechanism for the formation of coherent groups, such as wood mouse sperm trains [11], the patterns differ markedly from the large-scale swirling found in experiments on ram sperm suspensions [14, 15]. In real samples, sperm cells have distributions of waveforms, undulation frequencies or geometries. Accordingly, the parameters appearing in our model should not necessarily be the same for all individuals and should be chosen to accurately describe variations from individual to individual across the population.

We explore how variations in individual dynamics across the population impact large-scale suspension dynamics by introducing stochastic modulations of the undulation frequency into the model, allowing flagellar waves to differ from individual to individual and to vary over time [46, 58]. As we solve for the flagellum dynamics, changing the frequency also leads to changes in the waveform. The individual frequencies are drawn from a log-normal distribution after time intervals of length T=2π/ω as described in the SI. We set the mean frequency to ω and, in view of matching experimental data [45, 58], set the standard deviation σω to be ω/5.

We repeat our simulations of \num1000 swimmers, but now allowing for stochastic variations in the frequencies. Snapshots of the suspension evolution are shown in Fig. 3C and corresponding videos can be found as part of the ESM. The stochastic variations in the waveforms suppress the ability of neighboring flagella to synchronize. Rather than forming clusters, the swimmers now form waves, swirls and vortices, as seen in concentrated, quasi-two-dimensional ram sperm suspensions [14]. By running an equivalent simulation starting from an isotropic initial configuration, it was verified that the qualitative features of the fully non-linear state are independent of the initial conditions. We also note that a sufficiently large σω is needed to successfully suppress aggregation and observe large-scale motion of the suspension. Simulations of smaller ensembles (see SI) show that lower values of σω still exhibit cluster formation (e.g., σω=ω/20 or σω=ω/50). In these cases, stochastic frequency changes are rarely large enough to break the synchronization of adjacent flagella.

While hydrodynamic interactions no longer lead to aggregation, they still play a strong role in the evolution of the suspension. We perform simulations where we have removed the hydrodynamic interactions by solving the mobility problem using resistive force-theory (RFT) instead of FCM. As describe in the SI, the drag coefficients appearing in the RFT are chosen as to closely match the swimming speeds and flagellar waves obtained with the full FCM model. Fig. 3B shows snapshots from two different RFT simulations with ν=1.42 and N=1000 where the initial swimmer directions are either all aligned or distributed isotropically. The corresponding videos are provided as part of the ESM. The time-evolution of these suspensions is very different from those with hydrodynamic interactions. Most strikingly, for the initially polar suspension, long bending-waves are entirely absent and at long times we observe âlaningâ states similar to those found in simulations of self-propelled rods without hydrodynamic interactions [29, 59, 9]. When the initial state is isotropic, laning is is not observed (Fig. 3B), indicating further differences between the RFT and full simulations.
Additionally, in thin films, the hydrodynamic interactions between the swimmers are longer ranged than they would be in bulk. As a result, we find that as we increase the film thickness, though the instability of the polar state takes longer to develop, we still observe (see SI) the eventual onset of large-scale motion qualitatively similar to that found in the thinner film cases.

Order parameters exhibit density dependence

The polar and nematic order of the suspension is dependent on the effective area fraction, ν. We perform a series of smaller stochastic simulations (σω=ω/5) ranging from \numrange170500 swimmers run to t=400T and compare the order parameters

S1(t)

=⟨−^x⋅^nn(t)⟩n,

(6)

S2(t)

=⟨2(^x⋅^nn(t))2−1⟩n,

(7)

which measure the global alignment with the initial swimming direction (S1) and the global nematic order with respect to the initial direction (S2) [22]. Here, ^nn is the n-th swimmer’s mean orientation and ⟨⋅⟩n denotes the average over all swimmers. Videos of the simulations are provided as ESM.

\includegraphics

[width=]panel_S12_v4.pdf

Figure 5: Evolution of the global order parameters S1 (A) and S2 (B) over 400 periods for initially aligned suspensions with σω=ω/5, L=8.86d, and ν=0.54−1.59. Simulations are performed with full hydrodynamic interactions (HI, green) and with steric interactions only (RFT, blue), which show the opposite trend with ν.

We find that the directional order, given by S1, decreases with time and the decay rate increases with swimmer density (Fig. 5A). The time-scale of the decay is comparable to that found for rod models [22, 23]. The nematic order, given by S2, initially decreases with time. However, for higher densities, we see large fluctuations emerge after the initial decay (Fig. 5B). The appearance of these fluctuations during the relatively gentle decay of S1 reflects periodic folding and rotation of large patches of aligned swimmers. Fig. 5 also shows S1 and S2 for initially polar RFT simulations with the same effective area fractions. With hydrodynamic interactions removed, we find that higher density suspensions retain their alignment for longer times due to enhanced local caging.

Flagellar undulations contribute significantly to the energy and velocity spectra

\includegraphics

[width = ]panel_2_spectra.png

Figure 6: Fluid energy and center of mass velocity spectra for suspensions with effective area fractions ν=0.54−1.59 in domains of size L=8.86d and L=13.29d with both fixed (σω=0) and stochastically (σω=ω/5) varying frequencies. For the cases where L=8.86d the averaging is performed over times between 200T and 400T, while for L=13.29d, the interval is 150T and 200T.
A Fluid energy spectrum S(k)/N. The error bars indicate the sample standard deviation. B Fluid energy spectrum ¯S(k)/N of the filtered velocity field from which the short time-scales have been removed. C Center of mass velocity spectrum for all stochastic simulations. The data are normalized by the total energy, E=∑m⟨|U(m)cm|2⟩t, for each case. The error bars indicate one standard error of the mean. D Center of mass velocity spectrum of the filtered center of mass velocity fields for the lowest and highest ν. The spectra for the respective unfiltered fields from C are included for comparison.

To further quantify the large-scale motion, we compute the fluid energy spectrum,

S(k)≡⟨k⟨|^u(k,t)|2⟩∥k∥=k⟩t,

(8)

where ^u(k,t) is the spatial Fourier transform of the fluid velocity and ⟨⋅⟩t denotes time-averaging while ⟨⋅⟩∥k∥=k signifies averaging over all wave vectors of magnitude k. Fig. 6A shows S(k)/N from simulations with both fixed and stochastically varying frequencies and for different values of ν. With the exception of the pronounced peak near the undulation wavelength, the spectra are similar to those found in experiments on concentrated sperm suspensions [14]. In particular, we also find that the spectra decay like k−3 for large k (Fig. 6A). Given that this power-law decay appears over sub-swimmer length scales, we associate it with random forcing of the fluid by the flagella. A calculation presented in the SI demonstrates that an identical decay is obtained when the fluid is driven by spatially uncorrelated forcing. This is different from the k−4 decay observed at small scales in suspensions of self-propelled rods [23]. This decay was attributed to the sharp jump in the forcing along the length of each swimmer. Thus, the details regarding how the swimmers propel themselves can lead to changes in the spectrum at length scales comparable to the swimmer size.

Further, we find that the fluid energy spectrum is dominated by contributions at the time-scale of flagellar undulations, even at long length scales. Despite qualitative differences in their dynamics and the quantitative differences in the evolution of the order parameters, S(k)/N is nearly identical for all stochastic simulations, regardless of ν. The collapse illustrates that S(k) scales with the number of swimmers. This is in contrast with large increases with ν seen at low k in simulations of self-propelled rod suspensions [23]. Surprisingly, the S(k)/N curves for the fixed and stochastically varying frequency simulations also do not differ drastically, despite the striking differences in their suspension evolution. We expect that this is due to the energy injected at the short time-scales associated with flagellar undulations being large compared with that due to the large-scale motion of the suspension. To limit the contribution of short time-scales to S(k), we apply a low-pass filter to ^u(k,t) by computing a running time-average of u(x,t) with a window of 8T before performing the transform. The spectra, ¯S(k)/N, computed using the filtered velocity fields are shown in Fig. 6B. With the short time-scales suppressed, the dependence on ν becomes visible. We now observe larger values of ¯S(k)/N at low k when the sperm density is high, revealing quantitatively the large-scale fluid motion at the time scale of suspension evolution. Additionally, the spectrum for the fixed-frequency simulations is now almost flat at small k indicative of the lack of long range correlations.

Following from previous work on bacterial suspensions [9, 59], we also examine the spectrum of the center of mass velocity field, Ucm(x,t)=∑Nm=1U(m)cm(t)δ(x−X(m)cm(t)), where X(m)cm and U(m)cm are the center of mass position and velocity, respectively, as defined in the SI. Fig. 6C shows the swimmer velocity spectrum, ⟨|^Ucm(k)|2⟩t, for different values of ν. For each ν, the lowest spatial modes (lowest k) provide the largest contribution to the total spectrum, with the contribution increasing as ν increases. From these maximum values, the spectra then quickly decrease with k to a flat profile at scales below the flagellum length, corresponding to the level of noise generated in constructing the swimmer velocity field.

Due to flagellar undulations, the swimmers’ center of mass velocities oscillate about their mean values at the time-scale of the undulation frequency. We can again remove contributions of the short time-scale by filtering the swimmers’ center of mass velocities using a running time-average with window size 8T. Fig. 6D shows the spectra ⟨|ˆ¯Ucm(k)|2⟩t associated with the filtered swimmer velocity field. Filtering reduces values of the swimmer velocity spectrum for all k, indicating contributions at short time-scales at all length scales. For lower values of ν, the spectrum is reduced by an approximately constant factor across all k, while for higher values of ν, the reduction is more pronounced at shorter scales.

Discussion and Conclusions

In this study, we have used numerical simulation to show that variability in sperm flagellum dynamics across a suspension can lead to large changes in collective dynamics. We have incorporated these variations through stochastic changes in the actuation frequency that, in turn, lead to changes in waveform and amplitude. In actual sperm suspensions, variations are also likely to be present in flagellum length and cell geometry, and further, the flagellar waveform is likely to be more complex and fully three-dimensional [47, 60]. We expect these features to have a similar, aggregation-limiting effect as the stochastic variation in frequency. As such, measurements quantifying the distributions of individual sperm properties could not only provide a better classification of the individual cells, but also a better understanding of their collective dynamics.

Our fixed frequency simulations show that flagellar synchronization leads to aggregation and clustering. Certain properties associated with sperm that are found to self-organize into coherent groups, such as wood mouse sperm, might promote synchronization. It is known that these sperm typically have flagella that are longer than those of sperm from other species. By allowing for larger amplitudes, longer flagella may enhance the higher-order, time-dependent flows that we find are responsible for aggregation. Longer, thinner flagella would also bend and flex with greater ease in response to flows generated by neighboring sperm, further reinforcing the tendency to synchronize and aggregate. Additionally, microtubule sliding driving the undulations may depend on the external load experienced by the flagellum. This could allow for a non-trivial coupling between the actuation mechanism and the surrounding flow field, which, in turn, may also promote synchronisation of neighbouring flagella and cluster formation.

The turbulence-like state achieved when there is sufficient variation in the undulation frequency is both quantitatively and qualitatively similar to that observed in suspensions of ram and bull semen. While longer time simulations would be needed to explore this state in more detail, simulations starting from either polar or isotropic states eventually exhibit similar qualitative dynamics and have the same fluid energy spectrum at the final simulation times. Both this and the many previous results obtained using reduced models suggest that the observed jets and swirls are features of a final turbulence-like state. The presence of this state has been empirically correlated with sample fertility [61]. We have shown that its onset depends on sperm density, and its observation could therefore be serving as an indicator of sperm count, a well-known measure of sample fertility.

While in the simulations presented here the fluid domain is fully three-dimensional, the motion of the sperm cells is restricted to a plane. Even though layering in sperm suspensions has been observed in experiments [14] and sperm are attracted to move along planar surfaces [62], the imposed two dimensional motion in our simulations may enhance the role of steric interactions. For example, in a real sample when two sperm cells collide, it is likely that one will pass over the other and they will both experience some change in their swimming directions. In our simulations, the sperm cannot pass by each other as easily and the collisions can dramatically affect their swimming directions, leading to near alignment or anti-alignment. The interactions of self-driven filaments [57] and the trajectories of interacting sperm pairs [34] have recently been explored in simulation. These studies have shown that perturbations from planarity can lead to more complex trajectories, as well as a dependence on the elasticity of the flagella.
With further advances in numerical methods for simulating flexible filaments and more accurate models for three-dimensional flagellar waveforms, it will be possible to explore these effects through direct simulation of sperm suspensions.
Our theoretical model for the elastic flagellum and the FCM approach to the mobility problem carry over to fully three dimensional simulations. While there will be additional computational costs associated with the increased size of the fluid domain, more limiting is the technical challenge involved in keeping track of the fully three-dimensional deformation of the flagellum and, at the same time, using implicit time integration to handle numerical stiffness. We are currently developing the numerical methods to address this challenge.

Additionally, it would be of interest to ascertain suspension dynamics at even larger length scales and for longer timescales than those that may be accessible through direct simulation. Continuum models incorporating the time dependence of the flows generated by the swimmers have recently been developed [63, 64, 65] and could possibly be adapted to capture the effects of flagellar undulations seen in our simulations and used to explore nonlinear suspension dynamics at these scales.

Along with variations amongst the cells, collective dynamics are likely to be influenced by more complex features of the environment, such as nearby boundaries that could curve sperm trajectories [58], or particles or other micro-scale structures in the surrounding fluid that are known to enhance swimming speeds of undulatory swimmers [35]. Non-Newtonian features of the surrounding fluid, such as viscoelasticity, typically associated with biological fluids have also been been found to promote clustering [66]. These, and other important, perhaps more complex effects, such as chemotaxis, could further contribute to the richness and diversity of sperm collective dynamics.

Competing interests

The authors declare no competing interests.

Authors’ Contributions

S.F.S. and E.E.K. designed research, performed research, analyzed
data, and wrote the paper.

Acknowledgements

The authors would like to thank Pierre Degond and Demetrios Papageorgiou for helpful discussions.

Funding

Supporting Information

see below

Supporting Information

Swimmer model

As each swimmer is identical, for clarity, we describe the model for a single swimmer and subsequently describe how it generalizes when considering a suspension. Recall that a swimmer is composed of two elements, its cell head and its flagellum, where the cell head is treated as a rigid oblate spheroid with semi-major and minor axes a and b, respectively. Using the arc-length parametrization where ^t(s)=dY/ds, the force and moment balance equations are

dΛds+f

=0

dMds+τD+^t×Λ+τ

=0.

where Λ is the tension and M=KB^t×d^t/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 Nflag segments of length ΔL. The segments have positions Yn and orientations ^tn for n=1,…,Nflag. 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

Λn+1/2−Λn−1/2ΔL+fn=0

(9)

Mn+1/2−Mn−1/2ΔL+

12^tn×(Λn+1/2+Λn−1/2)+τDn+τn=0,

(10)

where Mn+1/2=(KB/ΔL)^tn×^tn+1. Multiplying through by ΔL and introducing Fn=fnΔL and Tn=τnΔL as the total applied force and torque, respectively, on segment n, we may write Eqs. [9] and [10] as

FCn+Fn

=0,

(11)

TBn+TCn+TDn+Tn

=0,

(12)

where FCn=Λn+1/2−Λn−1/2, TBn=Mn+1/2−Mn−1/2, and TCn=(ΔL/2)^tn×(Λn+1/2+Λn−1/2). The driving torques, TDn, arising from the preferred curvature are given by

TDn=KB(κ(sn,t)−κ(sn+1,t))^z

where sn=(n−1/2)ΔL for n=1,…,Nflag, and

κ(s,t)

=K0sin(ks−ωt+φ)⋅{2(l−s)/l,s>l/21,s≤l/2,

as well as, κ(s1,t)=κ(sNflag+1,t)=0.

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,

Yn+1−Yn+ΔL2(^tn+^tn+1)=0.

(13)

In the discretized system, the boundary conditions at the free end are satisfied by taking ΛNflag+1/2=MNflag+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 Fn=FSn+FHn where FSn is the steric force and FHn is the hydrodynamic force on the segment. The only external torque on the segment is the viscous torque, Tn=THn

The steric force FSn on segment n is given by the generic repulsive barrier force

FSn

=−FS∑m∈I(n)((χRnm)2−|Yn−Ym|2(χRnm)2−R2nm)4(Yn−Ym)d,

where the sum runs over I(n), the set of all segments/heads, m, with |Yn−Ym|<χRnm and χ=1.1. The parameter Rnm 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 Rnm=2b, where as if n and m are both heads, we have Rnm=2a. Lastly, if n and m are a head-segment pair then Rnm=a+b. The parameter FS sets the strength of the repulsion at contact, which in our simulations is FS=9018KB/d2.

Force-coupling method

The force and torque balances, Eqs. (11) and (12) 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) [36, 37, 38], 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 [38], 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

η∇2u−∇p=

∑nFHnΔn(x)−12∑nTHn×∇Θn(x)

∇⋅u

=0,

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 FHn=−FCn−FSn, and THn=−TBn−TCn−TDn. The force and torques are projected to the fluid via the Gaussian functions, which for particle n are given by

Δn(x)

=(2πσ2Δ)−3/2exp[−|x−Yn|2/(2σ2Δ)],

Θn(x)

=(2πσ2Θ)−3/2exp[−|x−Yn|2/(2σ2Θ)].

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

Vn

=∫uΔn(x)d3x

Ωn

=12∫(∇×u)Θn(x)d3x.

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^z for all n and can introduce angle θn such that ^tn=(cosθn,sinθn). Therefore, after computing the Vn and Ωn for each segment and head, the positions and orientations are obtained by integrating in time

dYndt

=Vn

(14)

dθndt

=Ωn.

(15)

subject to the inextensibility constraints given by [13]. To do this numerically, we rely on a second-order implicit BDF scheme [67] and Broyden’s method [68] 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

X(m)cm

=1a2+Nflagb2⎛⎝a2Y(m)0+b2Nflag∑n=1Y(m)n⎞⎠.

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

U(m)cm

=1a2+Nflagb2⎛⎝a2V(m)0+b2Nflag∑n=1V(m)n⎞⎠.

The director of swimmer m is given by

^n(m)

=−1Nflag+1Nflag∑n=0^t(m)n.

Force multipole expansion for a single swimmer

The singular force density for a single swimmer is given by

fi(x)=∑nFniδ(x−Yn)+∑nGnij∂jδ(x−Yn),

where the asymmetric dipole is related to the torque through Gij=ϵijkTk/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, Xcm, and defining the position of the center of mass relative to the n-th segment/head as Rn=Xcm−Yn, we have

fi(x)

=

=(∑nFni)=0δ(x−Xcm)+

+(∑nGnij+FniRnj)∂jδ(x−Xcm)+

+(∑nGnijRnk+12FniRnjRnk)∂k∂jδ(x−Xcm)+

+ higher order terms.

The first term is zero because there is no external force on the swimmer. The remaining two terms are the force dipole and quadrupole,

Gij

≡∑nGnij+FniRnj

Kijk

≡∑nGnijRnk+12FniRnjRnk,

respectively, which drive the flow observed in the limit ∣∣x−Xcm∣∣≫d. We compute the symmetric, trace-less and anti-symmetric parts of these tensors, given by

GSij

=12(Gij+Gji)−13Gkkδij

GAij

=12(Gij−Gji)

KSijk

=12(Kijk+Kjik)−13Kmmkδij

KAijk

=12(Kijk−Kjik)

as functions of time. These are shown in Fig. 2 in the main text, with the exception of GA(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

μln

=log(ω2√ω2+σ2ω)

σln

=√log(δ2ω+1),

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 tAB when the frequency changes, we require

ks−ωAtAB+φA−(ks−ωBtAB+φB)∈2πZ.

This can be fulfilled by taking

φB=2π((ωB−ωA)tAB+φA2π−⌊(ωB−ωA)tAB+φA2π⌋),

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 [28]. 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

\includegraphics

[width=.5]no_clusters_1000_SI.pdf

Figure 7: Number of clusters over time in simulations of 1000 swimmers starting from a polar configuration. Time is normalized by the average swimming velocity U0 and the swimmer length d.

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

dβ(i,j)

≡∣∣X(i)cm−X(j)cm∣∣+βtan(12arccos(^n(i)⋅^n(j)))

=∣∣X(i)cm−X(j)cm∣∣+β⎛⎜
⎜⎝√1−^n(i)⋅^n(j))√1+^n(i)⋅^n(j)⎞⎟
⎟⎠

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=161tan(π/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. 7. 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 Lz→0 and, as a result, the instability to exhibit a faster growth rate for thinner films. Figs. 8 and 9 show results from simulations containing 85 swimmers in domains of linear size L=4.43d and for thicknesses Lz=0.277d,0.554d, and 1.107d. The effective area fraction is ν=1.08 and the swimmers are initially in a polar state.

\includegraphics

[width=.75]LzSweepPanel.png

Figure 8: Snapshots at t=100T of simulations with 85 swimmers all starting from a polar initial state for film thicknesses Lz=0.277d,0.554d,1.107d (left to right).\includegraphics

[width=.45]LzSweepS_1x.pdf
\includegraphics[width=.45]LzSweepS_2x.pdf

Figure 9: Order parameters as a function of time for the simulations of 85 swimmers in films with thicknesses Lz=0.277d,0.554d, and 1.107d. The swimmers are initially aligned with the −x direction.

In accordance with the slower formation of the bending waves as shown in Fig. 8, the slower decay (see Fig. 9) of the order parameters S1 and S2 (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 Lz used in Fig. 9 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.

\includegraphics

[width=1.]SigSweepPanel.png

Figure 10: Snapshots at times t=50T and t=200T of simulations with 85 initially aligned swimmer with different frequency distributions. Each simulation has domain of size L=4.43d and film thickness Lz=0.277d. From left to right, the frequency standard deviations are σω=ω/50,ω/20,ω/10,ω/5, while the mean frequency, ω=2π/T, and all other parameters are the same.

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. 10). 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) [2, 41] 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

Vn

=(α∥^tn^t⊤n+α⊥(I−^tn^t⊤n))Fn,

Ωn

=βTn,β=(8πb3η)−1,

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 Fig. 11.
For the cell heads, we have for all swimmers m

V(0)m

=γF(0)m,

Ω(0)m

=λT(0)m,

where γ and λ are mobility coefficients taken from FCM simulations of a single spheroid. The numerical values used in the simulations are

γ=0.0280(ηb)−1,

λ=0.002430η−1b−3,

α⊥=0.1100(ηb)−1,

α∥=0.1700(ηb)−1.

\includegraphics

[width=.5]waveform_HI_RFT.pdf

Figure 11: Waveform of an individual swimmer with full hydrodynamics (top) and RFT only (bottom) in the center of mass frame. The swimming direction is to the left and approximately half an undulation period is depicted. The lighter the color of the center-line, the further the snapshot lies in the past.

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 as

^g(k)

=∫R2g(x)e−ik⋅xd2x,

g(x)

=14π2∫R2^g(k)eik⋅xd2k.

We begin by considering a force density, f(x), restricted to a square domain,

fL(x)

=f(x)\mathbbm1Ω(x),

where Ω=[−L/2,L/2]×[−L/2,L/2] and \mathbbm1Ω is the indicator function defined on the set Ω. The Fourier transform of the restricted force density is then

^fL(k)

=∫R2fL(x)e−ik⋅xd2k.

The Fourier transform of the fluid velocity resulting from the restricted force density is given by

^uL(k)=1ηk2(I−^k^k⊤)^fL(k)

where (1/ηk2)(I−^k^k⊤) is the Fourier representation of the inverse Stokes operator, k=|k|, and ^k=k/k.

The fluid energy density spectrum associated with ^uL(k) is

SL(k)=k2πL2∫2π0⟨|^uL(k)|2⟩dθ,

where θ is related to the two-dimensional wave-vector through k=(kcosθ,ksinθ) and ⟨⋅⟩ denotes the ensemble average. The energy density spectrum of the infinite system will then be given by S(k)=limL→∞SL(k).

As

|^uL(k)|2

=^uL(k)⋅^uL(−k)

=1η2k4(I−^k^k⊤):(^fL(k)^f⊤L(−k)),

the fluid energy density spectrum is directly related to the forcing through

SL(k)=12πη2k3L2∫2π0(I−^k^k⊤):⟨^fL(k)^f⊤L(−k)⟩dθ.

For SL(k)∝k−3 to hold, we must have ⟨^fL(k)^f⊤L(−k)⟩ independent of k. If, in addition, we require statistically independent and identical force components that do not depend on θ, we have that

⟨^fL(k)^f⊤L(−k)⟩=AL2I.

where A is a constant independent of L. The L2 dependence on the system size guarantees that SL(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

⟨fL(x+y)f⊤L(x)⟩

=1L2∫Ω⟨fL(x+y)f⊤L(x)⟩d2x

=1L2⟨∫ΩfL(x+y)f⊤L(x)d2x⟩.

As the force density is strictly restricted to Ω, we may replace the integral over Ω with one over all of R2. Thus,

⟨fL(x+y)f⊤L(x)⟩

=1L2⟨∫R2fL(x+y)f⊤L(x)d2x⟩,

which, from the convolution theorem, may be expressed as

⟨fL(x+y)f⊤L(x)⟩

=14π2L2∫R2⟨^fL(k)^f⊤L(−k)⟩eik⋅yd2k.

Finally, using the expression for ⟨^fL(k)^f⊤L(−k)⟩ given above, we see that

⟨fL(x+y)f⊤L(x)⟩

=AL2I4π2L2∫R2e% ik⋅yd2k

=AIδ(y)

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.