1. Introduction
Flagella or cilia are ubiquitous organelles in eukaryotic cells. They play a crucial role in physiological processes in animals, such as cerebrospinal fluid flow (Olstad et al. Reference Olstad, Ringers, Hansen, Wens, Brandt, Wachten, Yaksi and Jurisch-Yaksi2019), maintenance of the circadian clock (Tu et al. Reference Tu2023), and mucociliary clearance in the respiratory system (Nawroth et al. Reference Nawroth, Van Der Does, Firth and Kanso2019). In the microscopic world, unicellular organisms use their flagella/cilia to efficiently forage and travel, employing diverse locomotion strategies (Tam & Hosoi Reference Tam and Hosoi2011; Omori et al. Reference Omori, Ito and Ishikawa2020). For instance, Paramecium is known to regulate its ciliary beating to swim forwards, swim backwards and reorient (Machemer & Eckert Reference Machemer and Eckert1973; Okamoto & Nakaoka Reference Okamoto and Nakaoka1994; Ishikawa & Hota Reference Ishikawa and Hota2006). Escherichia coli exhibits a run-and-tumble locomotion by bundling and unbundling its flagella (Berg & Brown Reference Berg and Brown1972). The biflagellate alga Chlamydomonas transitions between two swimming gaits by modulating the synchronisation of its two flagella (Polin et al. Reference Polin, Tuval, Drescher, Gollub and Goldstein2009). Spermatozoa, however, adopt a distinct locomotion strategy by assembling into bundles with their heads attached. This cooperative behaviour has been found to increase their swimming speed (Woolley et al. Reference Woolley, Crockett, Groom and Revell2009; Fisher & Hoekstra Reference Fisher and Hoekstra2010), providing potential advantages in sperm competition (Moore et al. Reference Moore, Komrskova, Jenkins and Breed2002; Immler Reference Immler2008).
The increased velocity along the average path of sperm bundles during one beat cycle has been attributed to either their straighter swimming trajectory (Fisher et al. Reference Fisher, Giomi, Hoekstra and Mahadevan2014; Pearce et al. Reference Pearce, Hoogerbrugge, Hook, Fisher and Giomi2018) or their synchronised flagellar beating when the difference between the flagellar beating phase
$\Delta \phi$
is a constant for beat cycles (Woolley et al. Reference Woolley, Crockett, Groom and Revell2009; Zhang et al. Reference Zhang, Klingner, Le Gars, Misra, Magdanz and Khalil2023). While it has been observed experimentally that in-phase flagellar synchronisation (
$\Delta \phi =0$
) can increase the swimming speed of paired sperm cells (Woolley et al. Reference Woolley, Crockett, Groom and Revell2009; Zhang et al. Reference Zhang, Klingner, Le Gars, Misra, Magdanz and Khalil2023), the computational investigation indicates that increased speed results only from anti-phase flagellar synchronisation (
$\Delta\phi=\pi$
) or large flagellar phase lags (
$\Delta \phi \gt \pi /4$
) (Cripe et al. Reference Cripe, Richfield and Simons2016). This discrepancy highlights the complexity of the collective dynamics of sperm bundles. Rich behaviours are also discovered in a simple system comprising two adjacent but separate microswimmers, e.g. hydrodynamic attraction/repulsion (Yang et al. Reference Yang, Elgeti and Gompper2008; Carichino et al. Reference Carichino, Drumm and Olson2021), alignment (Olson & Fauci Reference Olson and Fauci2015; Taketoshi et al. Reference Taketoshi, Omori and Ishikawa2020), oscillation (Pooley et al. Reference Pooley, Alexander and Yeomans2007; Carichino et al. Reference Carichino, Drumm and Olson2021) and synchronisation (Di Leonardo et al. Reference Di Leonardo, Búzás, Kelemen, Vizsnyiczai, Oroszi and Ormos2012; Tătulea-Codrean & Lauga Reference Tătulea-Codrean and Lauga2022; Samatas & Lintuvuori Reference Samatas and Lintuvuori2023). These behaviours depend on their waveforms and relative displacement, phase and orientation (Pooley et al. Reference Pooley, Alexander and Yeomans2007; Elfring & Lauga Reference Elfring and Lauga2011a). Furthermore, current models infer that the hydrodynamic synchronisation of co-swimming cells requires geometrically asymmetric waveforms or the presence of a viscoelastic fluid environment (Elfring & Lauga Reference Elfring and Lauga2009, Reference Elfring and Lauga2011b; Elfring et al. Reference Elfring, Pak and Lauga2010). Their swimming speed and efficiency would increase drastically if co-swimming cells mechanically adhere into pairs (Simons & Rosenberger Reference Simons and Rosenberger2021). These factors also contribute to the rich dynamics observed in other microbial and artificial swimming systems (Drescher et al. Reference Drescher, Dunkel, Cisneros, Ganguly and Goldstein2011; Elgeti et al. Reference Elgeti, Winkler and Gompper2015; Pramanik et al. Reference Pramanik, Verstappen and Onck2024).
Despite many studies on the system of multiple separate microorganisms, the influence of flagellar beating on the locomotion of paired spermatozoa remains largely unexplored experimentally. A few challenges may account for the insufficiency of investigation. First, only a tiny proportion of sperm cells form pairs, restricting the sample size for experimental observation. Second, sperm locomotion is influenced by complex mechanical and hydrodynamic cell–cell and cell–environment interactions, which depend on their flagellar beat patterns and external factors, e.g. the geometry of surrounding environments (Raveshi et al. Reference Raveshi, Abdul Halim, Agnihotri, O’Bryan, Neild and Nosrati2021), fluid viscoelasticity (Tung et al. Reference Tung, Lin, Harvey, Fiore, Ardón, Wu and Suarez2017; Zaferani et al. Reference Zaferani, Javi, Mokhtare, Li and Abbaspourrad2021) and chemoattractants (Friedrich & Jülicher Reference Friedrich and Jülicher2007; Li et al. Reference Li, Chakrabarti, Castilla, Mahajan and Saintillan2023; Zaferani & Abbaspourrad Reference Zaferani and Abbaspourrad2023). As a result, it is difficult to control the flagellar beat pattern and experimentally investigate its influence on the locomotion of sperm pairs. Experimental studies of paired spermatozoa have been limited to comparing out-of-phase and in-phase flagellar beat patterns due to these challenges (Woolley et al. Reference Woolley, Crockett, Groom and Revell2009; Zhang et al. Reference Zhang, Klingner, Le Gars, Misra, Magdanz and Khalil2023).
(a) Schematic of a pair of bovine spermatozoa swimming in a 20
$\unicode {x03BC}\textrm {m}$
deep chamber with their heads attached. The paired sperm cells at two instants are overlaid at the same frame of reference. Any point on the ith cell can be described by the position vector
$\boldsymbol {r}_{i}$
with respect to the laboratory frame. The ith cell experiences a hydrodynamic force
$\boldsymbol {f}_{i}$
at an arbitrary point
$\boldsymbol {r}_i$
. The heads of the sperm pair are attached, with no relative translational motion allowed. But they can oscillate relatively about the pivot point (head tip) with instantaneous angular velocities
$\boldsymbol { {\Omega }}_i,\ i = 1, 2$
. The instantaneous translational velocity of the sperm pair
$\boldsymbol {U}$
is represented by that of the pivot point. Inset: the comoving frame is spanned by the orthonormal unit vectors
$\boldsymbol {e}_1$
and
$\boldsymbol {e}_2$
, and its origin overlaps with the pivot point. The flagellar shape in the comoving frame can be described by the local tangent angle
$\psi$
. Forces
$F_{ij}$
exist between the ith head and the jth head, and are specified differently based on the simplification of the head–head attachment. Three models are developed to investigate the locomotion of paired spermatozoa. (b) In model 1, the head orientation difference
$\Delta \alpha$
is unconstrained, such that
$\Delta \alpha$
can be positive, zero or negative. (c) In model 2,
$\Delta \alpha$
is constrained such that
$\Delta \alpha \gt 0$
. (d) Time-varying head angular velocities
$\boldsymbol {{\Omega }}_i,\ i = 1, 2$
, are extracted from experiments and prescribed in model 3. Scale bar: 10
$\unicode {x03BC}\textrm {m}$
.
In our experiments, bovine sperm pairs of two cells with their heads attached were observed in a chamber with a half-depth
$h=10\,\unicode {x03BC}\textrm {m}$
(figure 1
a). They swam in a plane parallel to the boundary surface of the chamber with primarily planar flagellar beats, similar to the previous experimental observations (Winet et al. Reference Winet, Bernstein and Head1984; Woolley Reference Woolley2003; Woolley et al. Reference Woolley, Crockett, Groom and Revell2009). The paired sperm cells were experimentally observed to transition between different modes of flagellar synchronisation: in-phase (
$\Delta \phi =0$
), anti-phase (
$\Delta \phi =\pi$
) and lagged synchronisation (
$\Delta \phi$
is a constant not equal to 0 or
$\pi$
). To investigate the influence of flagellar phase lag
$\Delta \phi$
on the swimming of sperm pairs, we develop a three-dimensional hydrodynamic model, referred to as model 1. The attachment between the heads involves adhesion proteins on the head surface and is influenced by factors such as cyclic adenosine monophosphate (cAMP) and sperm antagglutin on the head surface (Lindahl & Sjöblom Reference Lindahl and Sjöblom1981; Flaherty et al. Reference Flaherty, Swann, Primakoff and Myles1993). Considering the as yet not fully understood mechanical coupling due to the head–head attachment, in model 1, we simplify the mechanical head–head coupling as an adhesive force
$\boldsymbol {F}^{\mathrm {a}}_{ij}$
between the ith head and jth head. In model 1, the relative oscillation between the heads in a sperm pair is not constrained, such that their orientation difference
$\Delta \alpha =\alpha _2-\alpha _1$
can be positive, zero or negative (figure 1
b). In addition, we experimentally observed that the relative oscillation between the heads of paired sperm cells started and paused seemingly randomly, consistent with previous experimental observations (Woolley et al. Reference Woolley, Crockett, Groom and Revell2009; Zhang et al. Zhang et al. Reference Zhang, Klingner, Le Gars, Misra, Magdanz and Khalil2023). Their orientation difference
$\Delta \alpha$
was always a positive value during swimming. Based on these experimental observations, we speculate that a steric force
$\boldsymbol {F}^{{\rm s}}_{ij}$
may exist between the heads, besides the adhesive force. Therefore, based on model 1, we develop a second model, referred to as model 2, to explore this potential case where both adhesive and steric forces exist between the sperm heads. In model 2, the relative oscillation between the heads is constrained, such that their orientation difference is
$\Delta \alpha \gt 0$
during the whole course of swimming (figure 1
c). Both models 1 and 2 are used to investigate the phase-lag dependence of the swimming performances of paired spermatozoa regarding three parameters: average swimming speed, average power consumption, and swimming efficiency. For each specific flagellar phase lag
$\Delta \phi$
, its value is prescribed and fixed during the simulation. In models 1 and 2, the head angular velocities
$\boldsymbol {{\Omega }}_i(t)$
(
$i=1,2$
) are determined based on force-balance and torque-balance conditions for sperm cells at the low Reynolds number. To further explore the influence of head oscillation on the swimming trajectory of sperm pairs, a third model (referred to as model 3) is developed by providing model 1 with experimentally measured head angular velocities
$\boldsymbol {{\Omega }}_i(t)$
(
$i=1,2$
) (figure 1
d). For more realistic representations, flagellar waveforms are reconstructed from our experimental measurements and prescribed in all our models.
2. Tracking of paired spermatozoa in a chamber
Cryopreserved bovine spermatozoa were obtained from Semex Inc (Guelph) and stored in liquid nitrogen. Semen straws were thawed in a
$37\,^\circ\mathrm {C}$
water bath for 2 min, before suspending the cells in 2 ml high glucose Dulbecco’s Modified Eagle’s Medium (DMEM, D6546 Sigma Aldrich). Sperm cells were washed twice by centrifugation at 300 g for 5 min, and resuspended in a 2 ml clean medium. Then 0.3 % Methyl cellulose (M0512, Sigma Aldrich) was added to increase the viscosity of the medium. Five microlitres of sperm suspension were pipetted into slides with chamber depth 20
$\unicode {x03BC}\textrm {m}$
for immediate videomicroscopy.
Videomicroscopy was performed in an inverted Nikon microscope with a FastCam SA1.1 high-speed camera and a 40
$\times$
objective in phase contrast mode, obtaining video sequences with 500 frames per second. In our experiments, paired spermatozoa with their heads attached were observed (see supplementary movie 1). We track the flagella using the customised script in Matlab. The algorithm detects first the head tip and then the junction between the head and flagellum, from which the orientation of each cell is derived. Subsequently, the flagellum of each cell is tracked using the method reported by Geyer et al. (Reference Geyer, Jülicher, Howard and Friedrich2013) and Riedel-Kruse et al. (Reference Riedel-Kruse, Hilfinger, Howard and Jülicher2007). The tracked images need to be examined manually, and modified when the detection sometimes fails, due to e.g. dirt particles and overlapped flagella. Along each flagellum, 45 points are tracked. These flagellum points are off to both sides of the flagellum’s centreline and are not equally spaced. A Savitzky–Golay filter with degree 3 and a span of 5 sequential flagellum points is used to filter the flagellum points. These filtered flagellum points are then interpolated with splines. The arc length of the flagellum is determined by summing the lengths of the splines, and points at equal distances of 0.25
$\unicode {x03BC}\textrm {m}$
along the flagellum are then determined. These equidistant points along the ith flagellum of a sperm pair are time-varying and used to determine their velocity
$\boldsymbol {v}_i$
relative to the ith head.
3. Characterisation of the locomotion of paired spermatozoa
Bovine spermatozoa are approximately 60
$\unicode {x03BC}\textrm {m}$
in length. Approximately
$0.1{-}3\,\%$
of them, depending on conditions, formed bundles, most of which were sperm pairs (Morcillo i Soler et al. Reference Morcillo i Soler, Hidalgo, Fekete, Zalanyi, Khalil, Yeste and Magdanz2022). Considering their approximately planar kinematics, we can describe their projected locomotion on the two-dimensional plane where they swim, and neglect the out-of-plane component, as shown in figure 1. The ith sperm head of a sperm pair can be described by its orientation
$\alpha _i(t)$
and the position of its head tip
$\boldsymbol{r}_{\rm p}(t)$
with respect to the laboratory frame. The flagellar shape can be described by the tangent angle
$\psi$
with respect to the comoving frame spanned by the orthonormal unit vectors
$\boldsymbol {e}_1$
and
$\boldsymbol {e}_2$
. Here,
$\boldsymbol {e}_1$
and
$\boldsymbol {e}_2$
are oriented along the long and short axes of the projection of the ellipsoidal head on the swimming plane, respectively, as illustrated in figure 1 (inset). The tangent angle
$\psi (l,t),\ 0\le {l}\le {2L}$
, is enclosed between
$\boldsymbol {e}_1$
and the local tangent vector to the flagellar centreline, where
$L$
is the half-length of the flagellum. The flagellar shape with respect to the laboratory frame at a time t can be characterised by
\begin{equation} \boldsymbol {r}_{\rm f}(l,t)=\boldsymbol{r}_{\rm p}(t)-{a}\boldsymbol {e}_{1}-\int _0^{l}\mathrm {d}l^{\prime}\left [\cos \psi (l^{\prime},t)\,\boldsymbol {e}_{1}+\sin \psi (l^{\prime},t)\,\boldsymbol {e}_{2}\right ], \end{equation}
where
$a$
is the major axis of the projection of the head on the swimming plane (Friedrich et al. Reference Friedrich, Riedel-Kruse, Howard and Jülicher2010). Here,
$\boldsymbol {r}_{\rm f}(0,t)$
corresponds to the head–flagellum junction, and
$\boldsymbol {r}_{\rm f}(2L,t)$
corresponds to the distal end of the flagellum. This description can be applied to the two flagella of a sperm pair. Any point on the ith cell with respect to the laboratory frame is represented by
$\boldsymbol {r}_i(s,t)$
, where
$s$
is the coordinate of this point.
4. Mechanical and hydrodynamic cell–cell interactions
In our experiments, paired bovine sperm cells, one head on top of the other, oscillated their heads during swimming. From the top view, a portion of the heads, e.g. the head tips, remained overlapping during the whole course of swimming (figure 1). The head tips were the pivot point about which the heads oscillated. The attachment between the heads allowed their relative oscillation within the swimming plane, but constrained their relative translational motion. Similar experimental observations have been reported previously (Woolley et al. Reference Woolley, Crockett, Groom and Revell2009; Zhang et al. Reference Zhang, Klingner, Le Gars, Misra, Magdanz and Khalil2023). The head of a bovine spermatozoon resembles an approximately flattened ellipsoid, with detailed cellular morphologies (Pesch & Bergmann Reference Pesch and Bergmann2006; Carvalho et al. Reference Carvalho, Silva, Sartori and Dode2013). In our models, we ignore the detailed morphologies of the bovine sperm head and only consider its three principal dimensions – length, width and height. Thus the sperm head is simplified as an ellipsoid with dimensions
$9\times5\times0.4$
$\unicode {x03BC}\textrm {m}$
(length
$\times$
width
$\times$
height) in our models, based on previous experimental measurements (Pesch & Bergmann Reference Pesch and Bergmann2006; Carvalho et al. Reference Carvalho, Silva, Sartori and Dode2013). The flagellum is approximately a tube with half-length
$L=25\,\unicode {x03BC}\textrm {m}$
and radius
$\rho=0.25$
$\unicode {x03BC}\textrm {m}$
(Pesch & Bergmann Reference Pesch and Bergmann2006; Carvalho et al. Reference Carvalho, Silva, Sartori and Dode2013). These dimensions for the flagellum are used in our models.
Previous computational studies have shown that the hydrodynamic interaction between two adjacent but separate flagella may lead them to swim away from each other in the three-dimensional space (Simons et al. Reference Simons, Fauci and Cortez2015; Carichino et al. Reference Carichino, Drumm and Olson2021). However, compared with hydrodynamic forces, the mechanical head–head coupling is strong, so that paired sperm cells continue to swim together (Woolley et al. Reference Woolley, Crockett, Groom and Revell2009; Zhang et al. Reference Zhang, Klingner, Le Gars, Misra, Magdanz and Khalil2023). Therefore, the mechanical head–head coupling and its influence on the swimming of paired spermatozoa cannot be ignored. We simplify the as yet not fully understood head–head coupling as a pair of adhesive forces on the pivot point
$\boldsymbol {F}^{{\rm a}}_{ij},\ i=1, 2$
, elaborated on in § 4.1. Here, the force on the ith head
$\boldsymbol {F}^{{\rm a}}_{ij}$
results from the interaction with the jth head, such that
$F^{{\rm a}}_{12}=-F^{{\rm a}}_{21}$
. For the second potential case where both adhesive and steric forces exist between the heads, we develop model 2 based on model 1, detailed in § 4.2.
4.1. Model of paired spermatozoa with adhesive forces between their heads
Let
$\boldsymbol {{\Omega }}_i(t),\ i=1, 2$
, denote the instantaneous angular velocity of the ith head about the pivot point, and let
$\boldsymbol {U}(t)$
denote the instantaneous translational velocity of the pivot point, with respect to the laboratory frame. The velocity
$\boldsymbol {v}_i(s,t)$
at an arbitrary point on the ith flagellum with respect to the comoving frame is obtained from our experimentally observed time-varying flagellar shapes. Given
$\boldsymbol {v}_i(s,t)$
, the velocity
$\boldsymbol {u}_i(s,t)$
at this point in the laboratory frame is
\begin{equation} \boldsymbol {u}_i(s,t)=\boldsymbol {v}_i(s,t)+\boldsymbol {U}(t)+\boldsymbol {{\Omega }}_i(t)\times \boldsymbol {r}^{\prime}_i(s,t), \end{equation}
where
$\boldsymbol {r}^{\prime}_i(s,t)$
is the position vector of this point in the comoving frame (Bayly et al. Reference Bayly, Lewis, Ranz, Okamoto, Pless and Dutcher2011). For the velocity of any point on the head, (4.1) is also applicable, where the first term on the right-hand side vanishes. The sperm cells are spatially discretised, as detailed in Appendix A. For the planar locomotion, the velocity
$^{k}\boldsymbol {u}_i(t)$
at the kth point on the ith cell in the laboratory frame can be represented by
\begin{equation} ^{k}\boldsymbol {u}_i(t)={^{k}\boldsymbol {v}_i(t)}+{^{k}{\boldsymbol{B}}_i(t)}\,\mathcal {U}(t), \end{equation}
where
\begin{equation} \begin{aligned} \mathcal {U}=\left [\begin{matrix}U_{x} \\ U_{y} \\ {\Omega }_1 \\ {\Omega }_2 \\ F^{{\rm a}}_{21x} \\[2pt] F^{{\rm a}}_{21y}\\ \end{matrix}\right ], \end{aligned} \qquad \begin{aligned} ^{k}{\boldsymbol{B}}_i = \left \{ \begin{array}{ll} \left [\begin{matrix} 1 & \quad 0 & \quad -\,{^{k}r^{\prime}_{1y}} & \quad 0 & \quad 0 & \quad 0\\[1pt] 0 & \quad 1 & \quad ^{k}r^{\prime}_{1x} & \quad 0 & \quad 0 & \quad 0\\[1pt] 0 & \quad 0 & \quad 0 & \quad 0 & \quad 0 & \quad 0\\ \end{matrix}\right ], & \quad i=1, \\[20pt] \left [\begin{matrix} 1 & \quad 0 & \quad 0 & \quad -\,{^{k}r^{\prime}_{2y}} & \quad 0 & \quad 0\\[1pt] 0 & \quad 1 & \quad 0 & \quad ^{k}r^{\prime}_{2x} & \quad 0 & \quad 0 \\[1pt] 0 & \quad 0 & \quad 0 & \quad 0 & \quad 0 & \quad 0\\ \end{matrix}\right ], & \quad i=2. \end{array} \right . \end{aligned} \end{equation}
Here, the subscripts
$x$
and
$y$
in the variables represent their components along the
$x$
- and
$y$
-axes, respectively. The angular speeds
${\Omega }_1$
and
${\Omega }_2$
are defined such that
$\boldsymbol {{\Omega }}_1={\Omega }_1\boldsymbol {e}_3$
and
$\boldsymbol {{\Omega }}_2={\Omega }_2\boldsymbol {e}_3$
, where
$\boldsymbol {e}_3=\boldsymbol {e}_1\times \boldsymbol {e}_2$
.
The microscale geometry and swimming speed of spermatozoa render them low-Reynolds-number swimmers. They are generally deemed neutrally buoyant, as their swimming speed dominates their sedimentation speed (Gong et al. Reference Gong, Rode, Kaupp, Gompper, Elgeti, Friedrich and Alvarez2020). Consequently, the swimming paired sperm cells are force-free and torque-free (Lauga & Powers Reference Lauga and Powers2009). The ith cell in a sperm pair satisfies
\begin{equation} \int _s\boldsymbol {f}_i(s,t)\,\mathrm {d}s + \boldsymbol {F}^{{\rm a}}_{ij}(t) = \boldsymbol {0}, \quad \int _s\boldsymbol {r}_i(s,t)\times \boldsymbol {f}_i(s,t)\,\mathrm {d}s + \boldsymbol {T}^{{\rm a}}_{ij}(t) = \boldsymbol {0}, \end{equation}
where
$\boldsymbol {f}_i(s,t)$
is the hydrodynamic force at point
$\boldsymbol {r}_i(s,t)$
, as illustrated in figure 1. The adhesive torque on the ith head about the origin of the laboratory frame is
$\boldsymbol {T}^{{\rm a}}_{ij} = \boldsymbol{r}_{\rm p}\times \boldsymbol {F}^{{\rm a}}_{ij}$
. The total torque balances about any point.
The loss modulus of the fluid in our experiments dominates the storage modulus (Morcillo i Soler et al. Reference Morcillo i Soler, Hidalgo, Fekete, Zalanyi, Khalil, Yeste and Magdanz2022; Zhang et al. Reference Zhang, Klingner, Le Gars, Misra, Magdanz and Khalil2023). Consequently, we can neglect the elasticity of the fluid and calculate the hydrodynamic force
$\boldsymbol {f}_i(s,t)$
using Stokes equations. To incorporate the wall effect due to the top and bottom slides of the chamber, we introduce two walls in our models, which are at
$z=0\,\unicode {x03BC}\textrm {m}$
and
$z=20\,\unicode {x03BC}\textrm {m}$
, respectively, parallel with the
$x$
–
$y$
plane. According to previous studies, most spermatozoa are within 0.2 times their body length from the bottom surface (Winet et al. Reference Winet, Bernstein and Head1984; Elgeti et al. Reference Elgeti, Kaupp and Gompper2011). We therefore assume that the sperm cells are in the middle between the walls in our models. In the experiments, the sperm heads are physically attached. In the models, they need to be prevented from concurrently occupying the same space. Therefore, a pair of sperm cells is positioned on two parallel two-dimensional planes, with minimum distance 0.25
$\unicode {x03BC}\textrm {m}$
between their heads. Given that the height of the head is 0.4
$\unicode {x03BC}\textrm {m}$
, the two swimming planes are set at
$z=9.675\,\unicode {x03BC}\textrm {m}$
and
$z=10.325\,\unicode {x03BC}\textrm {m}$
, respectively. This arrangement ensures that the two cells remain in close proximity without overlapping in space.
To determine the hydrodynamic force
$\boldsymbol {f}_i$
on the ith cell, we use the regularised Stokeslets method, which is an effective approximation of the low-Reynolds-number flow. We discretise the ith cell into
$n_i$
points, and the walls into
$n_{{\rm w}}$
points, as detailed in Appendix A. Thus a total of
$n=n_1+n_2+n_{{\rm w}}$
points are included in our models. The relationship between the force
$\boldsymbol {f}_i$
and the velocity
$\boldsymbol {u}_i$
can be expressed in a vector form
\begin{equation} \left [\begin{matrix}\boldsymbol {u}_1^{\rm T} & \quad \boldsymbol {u}_2^{\rm T} & \quad \boldsymbol {u}_{{w}}^{\rm T}\end{matrix}\right ]^{\rm T}=-{\boldsymbol{A}}\left [\begin{matrix}\boldsymbol {f}_1^{\rm T} & \quad \boldsymbol {f}_2^{\rm T} & \quad \boldsymbol {f}_{{\rm w}}^{\rm T}\end{matrix}\right ]^{\rm T}, \end{equation}
where
$\boldsymbol {f}_{{\rm w}}$
is the hydrodynamic force on the stationary wall surfaces, and
${\boldsymbol{A}}$
is a
$3n\times 3n$
matrix (Appendix B). The velocity of the
$n_{{\rm w}}$
points on the wall surfaces
$\boldsymbol {u}_{{\rm w}}$
can be specified by a
$3n_{{\rm w}}\times 1$
zero matrix,
$\boldsymbol {u}_{{\rm w}}={\boldsymbol{0}}_{3n_{\textrm{w}}\times {1}}$
. Combining (4.2)–(4.5), we derive a linear system to describe the dynamics of paired sperm cells,
\begin{equation} \left [\begin{matrix}\boldsymbol {v} \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}\right ]=-\left [\begin{matrix}{\boldsymbol{A}} & \quad {\boldsymbol{B}}\\ {\boldsymbol{C}} & \quad {\boldsymbol{D}} \\ {\boldsymbol{E}} & \quad {\boldsymbol{F}} \\ {\boldsymbol{G}} & \quad {\boldsymbol{H}}\\ {\boldsymbol{I}} & \quad {\boldsymbol{J}}\\ {\boldsymbol{K}} & \quad {\boldsymbol{L}}\\ {\boldsymbol{M}} & \quad {\boldsymbol{N}}\\\end{matrix}\right ]\left [\begin{matrix}\boldsymbol {f}\\ \mathcal {U}\end{matrix}\right ], \end{equation}
where
${\boldsymbol{B}}= [\begin{matrix} ^{1}{\boldsymbol{B}}_1^{\mathrm{T}} & \quad ^{2}{\boldsymbol{B}}_1^{\mathrm{T}} & \quad \cdots & \quad ^{n_1}{\boldsymbol{B}}_1^{\mathrm{T}} & \quad ^{1}{\boldsymbol{B}}_2^{\mathrm{T}} & \quad ^{2}{\boldsymbol{B}}_2^{\mathrm{T}} & \cdots & \quad ^{n_2}{\boldsymbol{B}}_2^{\mathrm{T}} & \quad {\boldsymbol{0}}_{6\times {3n_{{\rm w}}}} \end{matrix} ]^{\mathrm{T}}$
, the velocity is
$\boldsymbol {v}= [\begin{matrix} \boldsymbol {v}_1^{\mathrm{T}} & \quad \boldsymbol {v}_2^{\mathrm{T}} & \quad {\boldsymbol{0}}_{1\times {3n_{{\rm w}}}} \end{matrix} ]^{\mathrm{T}}$
, and the force is
$\boldsymbol {f}= [\begin{matrix} \boldsymbol {f}_1^{\mathrm{T}} & \quad \boldsymbol {f}_2^{\mathrm{T}} & \quad \boldsymbol {f}_{{\rm w}}^{\mathrm{T}} \end{matrix} ]^{\mathrm{T}}$
. The last six equations in (4.6) represent the force balance along the x- and y-axes, and the torque balance along the z-axis, respectively, for each pair of sperm cells. According to the force-balance and torque-balance conditions, the blocks
${\boldsymbol{C}},{\boldsymbol{D}},\ldots ,{\boldsymbol{N}}$
of this linear system can be derived (Appendix C). In (4.6), velocity
$\boldsymbol {v}$
at any time is known, which is provided from experimental observations, but
$\boldsymbol {f}$
and all the variables in
$\mathcal {U}$
are unknown and need to be determined from (4.6). The velocity
$\boldsymbol {u}_i$
at any point on the ith sperm pair with respect to the laboratory frame can be further calculated from (4.2) after
$\mathcal {U}$
is determined.
4.2. Model of paired spermatozoa with adhesive and steric forces between their heads
To prevent sperm cells from unrealistically passing through each other, a steric force has been introduced and modelled as an elastic spring with a specific form (Cripe et al. Reference Cripe, Richfield and Simons2016; Simons & Rosenberger Reference Simons and Rosenberger2021). The steric force in their model acts only over a very short range to avoid significantly influencing sperm swimming. Likewise, instead of focusing on the specific value of steric force, we expect that the steric force in our model 2 can effectively repel the heads when they approach a total overlap, while having a minimum influence on their relative oscillation when their orientation difference is
$\Delta \alpha \gt 0$
. To this end, we develop model 2 based on model 1. We first consider a scenario in which both heads have the same angular velocity, i.e.
$\boldsymbol {{\Omega }}_1=\boldsymbol {{\Omega }}_2=\boldsymbol {{\Omega }}$
. The velocity
$^{k}\boldsymbol {u}_i(t)$
at the kth point on the ith cell in the laboratory frame can continue to be represented by (4.2), but the matrices
$\mathcal {U}$
and
$^{k}{\boldsymbol{B}}_i$
need to be modified to
\begin{equation} \begin{aligned} \mathcal {U}^{\prime}=\left [\begin{matrix}U_{x} \\ U_{y} \\ {\Omega } \end{matrix}\right ], \end{aligned} \quad \begin{aligned} ^{k}{\boldsymbol{B}}^{\prime}_i = \left [\begin{matrix} 1 & \quad 0 & \quad -\,{^{k}r^{\prime}_{iy}} \\0 & \quad 1 & \quad ^{k}r^{\prime}_{ix} \\0 & \quad 0 & \quad 0 \\\end{matrix}\right ], \quad i=1, 2, \end{aligned} \end{equation}
respectively. As
$\boldsymbol {F}^{\mathrm {a}}_{12}=-\boldsymbol {F}^{\mathrm {a}}_{21}$
and
$\boldsymbol {F}^{\mathrm {s}}_{12}=-\boldsymbol {F}^{\mathrm {s}}_{21}$
, the force-balance and torque-balance conditions on the whole sperm pair are
\begin{equation} \sum _{i=1}^2\int _s\boldsymbol {f}_i(s,t)\,\mathrm {d}s = \boldsymbol {0}, \quad \sum _{i=1}^2\int _s\boldsymbol {r}_i(s,t)\times \boldsymbol {f}_i(s,t)\,\mathrm {d}s = \boldsymbol {0}. \end{equation}
Combining (4.2), (4.5), (4.7) and (4.8), we derive a linear system to describe the dynamics of paired sperm cells with the same head angular velocity
$\boldsymbol {{\Omega }}$
,
\begin{equation} \left [\begin{matrix}\boldsymbol {v} \\ 0 \\ 0 \\ 0 \end{matrix}\right ]=-\left [\begin{matrix}{\boldsymbol{A}} & {\boldsymbol{B}}^{\prime} \\ {\boldsymbol{C}}^{\prime} & {\boldsymbol{D}}^{\prime} \\ {\boldsymbol{E}}^{\prime} & {\boldsymbol{F}}^{\prime} \\ {\boldsymbol{G}}^{\prime} & {\boldsymbol{H}}^{\prime} \end{matrix}\right ]\left [\begin{matrix}\boldsymbol {f}\\ \mathcal {U}^{\prime}\end{matrix}\right ], \end{equation}
where
${\boldsymbol{B}}^{\prime}= [\begin{matrix} ^{1}{{\boldsymbol{B}}^{\prime}_1}^{\mathrm{T}} & ^{2}{{\boldsymbol{B}}^{\prime}_1}^{\mathrm{T}} & \cdots & ^{n_1}{{\boldsymbol{B}}^{\prime}_1}^{\mathrm{T}} & ^{1}{{\boldsymbol{B}}^{\prime}_2}^{\mathrm{T}} & ^{2}{{\boldsymbol{B}}^{\prime}_2}^{\mathrm{T}} & \cdots & ^{n_2}{{\boldsymbol{B}}^{\prime}_2}^{\mathrm{T}} & {\boldsymbol{0}}_{3\times {3n_{\textrm{w}}}} \end{matrix} ]^{\mathrm{T}}$
, the velocity is
$\boldsymbol {v}= [\begin{matrix} \boldsymbol {v}_1^{\mathrm{T}} & \boldsymbol {v}_2^{\mathrm{T}} & {\boldsymbol{0}}_{1\times {3n_{\textrm{w}}}} \end{matrix} ]^{\mathrm{T}}$
, and the force is
$\boldsymbol {f}= [\begin{matrix} \boldsymbol {f}_1^{\mathrm{T}} & \boldsymbol {f}_2^{\mathrm{T}} & \boldsymbol {f}_{\textrm{w}}^{\mathrm{T}} \end{matrix} ]^{\mathrm{T}}$
. The last three equations in (4.9) represent the force balance along the x- and y-axes, and the torque balance along the z-axis, for the sperm pair, respectively. According to the force-balance and torque-balance conditions, the blocks
${\boldsymbol{C}}^{\prime},{\boldsymbol{D}}^{\prime},\ldots ,{\boldsymbol{H}}^{\prime}$
of this linear system can be derived (Appendix D). In (4.9), velocity
$\boldsymbol {v}$
at any time is known, which is provided from experimental observations, but
$\boldsymbol {f}$
and all the variables in
$\mathcal {U}^{\prime}$
are unknown and need to be determined from (4.9). The velocity
$\boldsymbol {u}_i$
at any point on the ith sperm pair with respect to the laboratory frame can be further calculated from (4.2) after
$\mathcal {U}^{\prime}$
is determined.
Our simulation consists of
$n_{{\rm t}}$
discrete time steps. If at the beginning of the jth time step, the position
$\boldsymbol {r}(t_j)$
of any point on the surface of the sperm cells and the velocity
$\boldsymbol {v}(t_{j})$
are known, then
$\boldsymbol {r}(t_{j+1})$
,
$\Delta \alpha (t_{j+1})$
,
$\boldsymbol {f}(t_{j+1})$
,
$\boldsymbol {U}(t_{j+1})$
,
$\boldsymbol {{\Omega }}_1(t_{j+1})$
and
$\boldsymbol {{\Omega }}_2(t_{j+1})$
at the beginning of the
$(j+1)$
th time step can be determined by solving (4.6) or (4.9). In our simulations, the initial position
$\boldsymbol {r}(t_1)$
is given, and the velocity
$\boldsymbol {v}$
at any time step is provided from experimental observations. Assuming that the initial orientation difference between the heads is
$\Delta \alpha (t_1)\geq 0$
, now we can formulate model 2 by combining (4.6) and (4.9) following the method described in Algorithm1.
Model 2.
In model 2 (i.e. Algorithm1), sperm heads are fused with zero relative angular velocity,
$\Delta \boldsymbol {{\Omega }}=\boldsymbol {{\Omega }}_2-\boldsymbol {{\Omega }}_1=\boldsymbol {0}$
, when
$\Delta \alpha$
is within a very small range close to zero. This range depends on the variation of
$\Delta \alpha$
over one discrete time step. Note that both adhesive and steric forces between the heads are intrinsically included in model 2, constraining the relative oscillation of the heads such that
$\Delta \boldsymbol {{\Omega }}=\boldsymbol {0}$
when
$\Delta \alpha$
is within the range, despite the specific values of the forces and where they are applied being unknown. Only an adhesive force between the heads is included in model 2 when
$\Delta \alpha$
is out of the range. Thereby, we finish developing model 2, in which sperm heads are effectively prevented from totally overlapping, while having a minimum influence on their relative oscillation when
$\Delta \alpha \gt 0$
.
4.3. Validation of the models
Model 3 can be developed based on model 1 (i.e. (4.6)) by reducing the vector
$\mathcal {U}$
in model 1 to
$\mathcal {U}= [\begin{matrix}U_{x} & U_{y} & F^{\mathrm {a}}_{21x} & F^{\mathrm {a}}_{21y}\\\end{matrix} ]^{\mathrm{T}}$
and removing the two equations associated with the torque-balance condition from (4.6). Our models build on regularised Stokeslets method, involving a choice of regularised parameter
$\epsilon$
and spatial discretisation of the sperm cells and walls. To validate our discretisation of the head and flagellum, following the approach of Cortez et al. (Reference Cortez, Fauci and Medovikov2005) and Gillies et al. (Reference Gillies, Cannon, Green and Pacey2009), we compare our predicted resistive force on an isolated head and an isolated flagellum translating in an unbounded fluid with their exact solution. The non-dimensional force error is
$e_{{\rm f}}\lt 2.4\times 10^{-4}$
(Appendix E). In addition, we compute the cumulative distribution of velocity error for a sphere translating in an unbounded fluid with velocity
$U_z=1\,\rm m\,s^-{^1}$
. All the points on the surface of the sphere have a velocity error
$e_{{\rm v}}=|U_{z}-1|\lt 0.018\,\rm m\,s^-{^1}$
(Appendix E). To regularise the walls in our models, the regularised parameter
$\epsilon$
is chosen in the form
$\epsilon =\xi \rho ^{m}$
(Ainley et al. Reference Ainley, Durkin, Embid, Boindala and Cortez2008; Gillies et al. Reference Gillies, Cannon, Green and Pacey2009). In our test cases,
$\xi =0.5$
and
$m=0.9$
are chosen, which yields predictions that closely match the published results (Appendix F). Therefore, these values of
$\xi$
and
$m$
are used in our models for all subsequent simulations. The walls in our models are square. The simulations in the test cases are insensitive to the wall size when the side length of the walls ranges from 60 to 160
$\unicode {x03BC}\textrm {m}$
, as shown in Appendix F. Therefore, the side length of each wall is set to 80
$\unicode {x03BC}\textrm {m}$
in our models, which is fixed in all our subsequent simulations.
4.3.1. Force and torque errors
We use our models to simulate the locomotion of a sperm pair. In the simulations, the viscosity is chosen as
$\mu =1\,\mathrm {Pa\,s}$
, and the flagellar beat patterns are prescribed based on the experimental measurements. For model 3, prescribed head angular velocities
$\boldsymbol {{\Omega }}_i(t),$
$i=1, 2$
, are further required.
Regarding models 1 and 2, force-balance conditions along the x- and y-axes and torque-balance conditions along the z-axis are applied. As the mechanical coupling between the heads is a pair of internal forces, the total hydrodynamic force on the whole sperm pair along the x- and y-axes should satisfy
$F_x=\sum _{i=1}^2F_{ix}=0$
and
$F_y=\sum _{i=1}^2F_{iy}=0$
, where
$F_{ix}$
and
$F_{iy}$
are the x- and y-components of the total hydrodynamic force on the ith cell, respectively. Likewise, the total hydrodynamic torque on the whole sperm pair along the z-axis should satisfy
$T_z=\sum _{i=1}^2T_{iz}=0$
, where
$T_{iz}$
is the z-component of the total hydrodynamic torque on the ith cell. As shown in figure 2, the absolute values
$|F_x|$
,
$|F_y|$
and
$|T_z|$
are minimal. The minimal values result from the numerical error. The absolute value of the total hydrodynamic force on the sperm pair along the z-axis,
$|F_z|$
, is also very small. We take the total hydrodynamic force on cell 1,
$|F_1|$
, as the benchmark, and
$|F_z|$
is three orders of magnitude lower than
$|F_1|$
, despite no force-balance condition along the z-axis applied to the sperm cells. Likewise, taking
$|T_1|$
as the benchmark, the absolute values of the total hydrodynamic torque on the sperm pair along the x- and y-axes,
$|T_x|$
and
$|T_y|$
, are three orders of magnitude lower than
$|T_1|$
, despite no torque-balance conditions along the x- and y-axes applied to the sperm cells.
Regarding model 3, force-balance conditions along the x- and y-axes are applied to the sperm cells. The absolute values
$|F_x|$
and
$|F_y|$
are minimal due to the numerical error, as shown in figure 2(a,iii). In addition,
$|F_z|$
is three orders of magnitude lower than
$|F_{1}|$
, and
$|T_x|$
,
$|T_y|$
and
$|T_z|$
are three orders of magnitude lower than
$|T_{1}|$
, despite no force-balance condition along the z-axis and no torque-balance conditions along the x-, y- and z-axes applied to the sperm pair.
(a) The absolute values of the x-, y- and z-components of the total hydrodynamic force on a sperm pair in models 1 (i), 2 (ii), and 3 (iii) are minimal compared to the absolute value of the total hydrodynamic force on cell 1 of this sperm pair,
$|F_1|$
. (b) The absolute values of the x-, y- and z-components of the total hydrodynamic torque on the sperm pair in models 1 (i), 2 (ii), and 3 (iii) are minimal compared to the absolute value of the total hydrodynamic torque on cell 1 of this sperm pair,
$|T_1|$
. The force and torque components during a shorter time range, 0–0.1 s, are shown in the insets.
4.3.2. Comparison of the predicted and experimentally observed trajectories of sperm pairs
Three sperm pairs, referred to as (a) sperm pair 1, (b) sperm pair 2, and (c) sperm pair 3, were tracked using high-speed videomicroscopy. In each sperm pair, the trajectory of the pivot point is predicted using models 1, 2 and 3, respectively. (d) The head orientations of cell 1 (
$\alpha _1$
) and cell 2 (
$\alpha _2$
) in sperm pair 3 are time-varying. Orientations in model 3 overlap with the experimentally measured ones. (e) Regarding the pivot point in sperm pair 3, the predictions and experimental observation have similar tendencies in the displacement
$d$
. The predictions and experimental observation have similar tendencies in the travelling distance
$S$
.
Following previous studies (Yundt et al. Reference Yundt, Shack and Lardner1975; Gillies et al. Reference Gillies, Cannon, Green and Pacey2009), we assess our models by comparing the predicted and experimentally observed swimming trajectories. Three sperm pairs are simulated using our models. The trajectory of each sperm pair is the time-varying position of their pivot point, as shown in figure 3. We calculate the displacement
$d$
and travelling distance
$S$
of the pivot point during one flagellar beat cycle
$T$
. The displacement during the jth flagellar beat cycle
$d$
is defined as
\begin{equation} d = \frac {1}{n_{\mathrm {t}}}\sum _{i=1}^{n_{\mathrm {t}}}\left |\boldsymbol{r}_{\rm p}(t_i+T)-\boldsymbol{r}_{\rm p}(t_i)\right |, \end{equation}
where
$n_{\mathrm {t}}$
is the number of the time steps during one flagellar beat cycle
$T$
. The displacement
$|\boldsymbol{r}_{\rm p}(t_i+T)-\boldsymbol{r}_{\rm p}(t_i)|$
depends on the initial time
$t_i$
during one beat cycle. We thus define
$d$
by averaging
$|\boldsymbol{r}_{\rm p}(t_i+T)-\boldsymbol{r}_{\rm p}(t_i)|$
over one beat cycle using (4.10). The travelling distance of the pivot point during the jth flagellar beat cycle
$S$
is defined as
$S=\sum _{i=1}^{n_{\mathrm {t}}}\left |\boldsymbol{r}_{\rm p}(t_i)-\boldsymbol{r}_{\rm p}(t_{i-1})\right |$
. Taking sperm pair 3, for example (figure 3
e), there are 40 flagellar beat cycles. Compared with the experimentally observed displacement averaged over 40 flagellar beat cycles
$\langle d_{\mathrm {obse}}\rangle$
, the predicted one,
$\langle d_{\mathrm {theo}}\rangle$
, in models 1, 2 and 3, is
$0.2\,\%$
shorter,
$7.8\,\%$
longer, and
$11.7\,\%$
longer, respectively. The predicted travelling distances of the pivot point averaged over the 40 beat cycles,
$\langle S_{\mathrm {theo}}\rangle$
, in models 1, 2 and 3 are
$32.1\,\%$
,
$29.4\,\%$
and
$25.1\,\%$
, respectively, shorter than the experimentally observed one
$\langle S_{\mathrm {obse}}\rangle$
.
In our experiments, the trajectory of the pivot point wiggles around its average path, resulting in the time-varying longitudinal and lateral displacement within one flagellar beat cycle. Our models predict the fine structure of the locomotion of the sperm pairs such as the wiggling. However, the agreement between the predicted and experimentally observed swimming trajectories varies significantly between the models. We quantify the agreement between the predicted and experimentally observed trajectories using Fréchet distance
$d_{\mathrm {Fr}}$
, a lower value of which indicates a higher similarity (Eiter & Mannila Reference Eiter and Mannila1994). Taking sperm pair 3, for example, compared with the experimentally observed swimming trajectory,
$d_{\mathrm {Fr}}$
is 45.6, 56.6 and 21.1
$\unicode {x03BC}\textrm {m}$
for the predicted trajectory in models 1, 2 and 3, respectively. As shown in figure 3
c, compared with models 1 and 2, model 3 predicts a trajectory with higher similarity to the experimentally observed one, revealing the significant role of head oscillation in the swimming trajectory of sperm pairs. The improved prediction accuracy results from the available correct information, i.e. the head angular velocities.
5. Swimming performances of paired spermatozoa
A swimming sperm pair was observed experimentally to exhibit three different modes of flagellar synchronisation: in-phase, anti-phase and lagged synchronisation. (a) Interrupted by phase slips, the sperm pair transitioned between the modes of flagellar synchronisation. Scale bar: 10
$\unicode {x03BC}\textrm {m}$
. (b) To characterise the flagellar beating, their shapes are mapped into a two-dimensional space spanned by the first two shape scores
$\beta _1$
and
$\beta _2$
. The phase of the flagellar beating is obtained through binning tracked flagellar shapes according to shape similarity. (c) The phase lag between the flagella of the sperm pair
$\Delta \phi$
is unwrapped and clearly shows the phase lags, slips and synchronisation. Inset: the flagellar phases
$\phi _1$
and
$\phi _2$
during a time interval when phase slips occur.
In our experiments, paired sperm cells beat their two flagella with the same period, exhibiting different modes of flagellar synchronisation, as shown in figure 4(a). To formalise the locomotion of the paired spermatozoa, the high-dimensional description of flagellar beating by (3.1) can be further mapped to a low-dimensional space spanned by the first two shape scores
$\beta _1$
and
$\beta _2$
using principal component analysis (Geyer et al. Reference Geyer, Jülicher, Howard and Friedrich2013; Werner et al. Reference Werner, Rink, Riedel-Kruse and Friedrich2014). The points in the
$(\beta _1,\beta _2)$
space represent a series of flagellar shapes and form a closed loop (figure 4
b). We define a limit cycle of the beating of the ith flagellum by fitting the closed loop parametrised by a phase
$\phi _i$
(Appendix G). The flagellar phase
$\phi _i$
is defined according to flagellar shape similarity such that it is independent of its time derivative, whereby a unique phase is assigned for each tracked flagellar shape (Geyer et al. Reference Geyer, Jülicher, Howard and Friedrich2013). Accordingly, a swimming sperm pair can be completely characterised with the position of the pivot point
$\boldsymbol{r}_{\rm p}$
, the orientation of the heads
$\alpha _i$
, and the phase of flagellar beating
$\phi _i$
. The flagellar synchronisation of paired cells is achieved when the phase lag between their flagellar beating
$\Delta \phi =\phi _2-\phi _1$
remains unchanged over one beat cycle. The sperm pairs were observed experimentally to maintain a constant flagellar phase lag
$\Delta \phi$
for several beat cycles, followed by phase slips, and then establish another constant flagellar phase lag, as shown in figure 4(c). This repeated process allowed the flagellar synchronisation to transition between in-phase, anti-phase and lagged synchronisation.
Three distinct sperm pairs respond differently to the flagellar phase lag
$\Delta \phi$
in (a) average swimming speed
$\langle U\rangle$
, (b) swimming power
$\langle P\rangle$
, and (c) normalised swimming efficiency
$\eta /\eta _0$
. Here,
$\eta _0$
is the efficiency when
$\Delta \phi =0$
. Models 1 and 2 are used to calculate
$\langle U\rangle$
,
$\langle P\rangle$
and
$\eta /\eta _0$
for each sperm pair. Each simulation includes 14 different values of flagellar phase lag
$\Delta \phi$
:
$0$
,
$\pi /6$
,
$\pi /3$
,
$\pi /2$
,
$2\pi /3$
,
$5\pi /6$
,
$\pi$
,
$7\pi /6$
,
$4\pi /3$
,
$3\pi /2$
,
$5\pi /3$
,
$11\pi /6$
,
$23\pi /12$
and
$95\pi /48$
.
Inspired by the experimental observations, we further investigate numerically the influence of the flagellar phase lag
$\Delta \phi$
on the swimming of paired cells. Both models 1 and 2 are used in our simulations. In our simulations, the flagellar beat patterns are reconstructed from their limit cycles (Appendix G). The phase-lag dependence of flagellar waveforms is neglected. We assess the swimming performances of sperm pairs regarding three parameters: average swimming speed
$\langle U\rangle$
, average power consumption
$\langle P\rangle$
, and swimming efficiency
$\eta$
. In our simulations, the three parameters are determined after the sperm pair has established regular beating, when its displacement during one flagellar beat cycle,
$d$
, stops varying between beat cycles. The average swimming speed over one beat cycle is determined by
$\langle U\rangle =d/T$
, representing the swimming speed along the average path. Here,
$d$
is calculated using (4.10). The average power consumption per cell
$\langle P\rangle$
for a sperm pair is determined by averaging the instantaneous power of cells 1 and 2 over one beat cycle (Appendix H). Compared to speed, swimming efficiency
$\eta$
may be more critical for spermatozoa due to their limited energy reserves. We define
$\eta =\langle U\rangle ^2/\langle P\rangle$
. This means that a longer swimming displacement per unit of energy consumption leads to higher efficiency. Given other conditions (e.g. the flagellar waveform due to the viscosity) unchanged,
$\langle U\rangle$
is independent of fluid viscosity
$\mu$
, whereas the average power
$\langle P\rangle$
depends linearly on it, and
$\eta$
has an inverse dependence on it. To compare
$\langle P\rangle$
and
$\eta$
at different fluid viscosities and various flagellar phase lags, we use
$\mu =1\,\mathrm {Pa\,s}$
in all our simulations for swimming performances, and normalise
$\eta$
with its value at
$\Delta \phi =0$
, denoted by
$\eta _0$
.
In all our simulations for swimming performances, the period of flagellar beating is the same for the cells within a sperm pair, but distinct among the pairs. The flagellar waveform varies between the pairs. Therefore, sperm pairs with the same flagellar phase lag exhibit distinct swimming performances, as shown in figure 5. In addition, a small variation in
$\Delta \phi$
may lead the swimming performances to change significantly, especially for
$\eta /\eta _0$
when the value of
$\langle P\rangle$
is small. For instance, the value of
$\eta /\eta _0$
for sperm pair 2 in model 1 significantly increases from 0.40 to 0.95 as
$\Delta \phi$
increases from
$23\pi /12$
to
$95\pi /48$
(figure 5
c). In all our simulations for swimming performances, compared with the small flagellar phase lags (approximately 0 radians), large flagellar phase lags (approximately
$\pi$
radians) lead to a higher
$\langle P\rangle$
, which aligns with previous studies (Elfring & Lauga Reference Elfring and Lauga2011a; Cripe et al. Reference Cripe, Richfield and Simons2016). However, a high swimming speed or efficiency does not necessarily result from flagellar phase lags close to
$\pi$
radians. For instance, in model 2, sperm pair 2 achieves its highest
$\langle U\rangle$
and
$\eta /\eta _0$
when its flagellar phase lag
$\Delta \phi$
is close to
$\pi /2$
radians. In addition, the swimming performances are influenced significantly by the mechanical head–head coupling. For the same sperm pair with the same
$\Delta \phi$
but different head–head coupling, it may have a large difference in
$\langle U\rangle$
,
$\langle P\rangle$
or
$\eta /\eta _0$
, as shown in figure 5. Taking sperm pair 2, for example, when
$\Delta \phi =\pi /3$
, the value of
$\langle U\rangle$
in model 1 is equal to that in model 2. But when
$\Delta \phi =\pi /2$
, the value of
$\langle U\rangle$
for sperm pair 2 in model 1 diverges significantly from that in model 2. One reason is that the orientation difference between the heads
$\Delta \alpha (t)$
is constrained in model 2 such that
$\Delta \alpha (t)\ge 0$
for the whole course of swimming, whereas no such constraint exists in model 1, as shown in figure 6(a). This constraint in the relative head oscillation leads to different swimming trajectories, as shown in figure 6(b). However, in our experiments, we did not observe the sperm pairs with their head orientation difference at
$\Delta \alpha (t)\lt 0$
during swimming. Therefore, steric force probably exists between the heads of the sperm pairs in the experiments, repelling each other.
For sperm pair 2, (a) its head orientation difference
$\Delta \alpha (t)$
and (b) its swimming trajectories during five flagellar beat cycles are shown. When
$\Delta \phi =\pi /3$
, sperm pair 2 in models 1 and 2 has the same time-varying
$\Delta \alpha (t)$
. When
$\Delta \phi =\pi /2$
, the head orientation difference
$\Delta \alpha (t)$
for sperm pair 2 in model 1 can be a negative value within a time interval, during which the relative oscillation of the heads of sperm pair 2 in model 2 are fused such that their relative angular speed is
$\Delta {\Omega }=0$
.
Furthermore, the flagellar waveform also varies between the cells in the same sperm pair, resulting in a difference in the power
$\langle P_1\rangle$
and
$\langle P_2\rangle$
, as shown in figure 7. In our simulations,
$\langle P\rangle$
for each sperm pair is always a positive value. But for sperm pair 2 in model 1, the average power for cell 1 is
$\langle P_1\rangle =-5.5$
pW when
$\Delta \phi =23\pi /12$
, as shown in figure 7(b).
Sperm pairs (a) 1, (b) 2 and (c) 3 are shown. Models 1 and 2 are used to calculate the average power consumption for cell 1 (
$\langle P_1\rangle$
) and cell 2 (
$\langle P_2\rangle$
) in each sperm pair. Each simulation includes 14 different values of flagellar phase lag
$\Delta \phi$
:
$0$
,
$\pi /6$
,
$\pi /3$
,
$\pi /2$
,
$2\pi /3$
,
$5\pi /6$
,
$\pi$
,
$7\pi /6$
,
$4\pi /3$
,
$3\pi /2$
,
$5\pi /3$
,
$11\pi /6$
,
$23\pi /12$
and
$95\pi /48$
.
6. Discussion
The collective behaviours of mechanically and hydrodynamically coupled flagellated cells are complex and not fully understood. Here, using paired bovine spermatozoa as an example, we develop hydrodynamic models to unveil more nuanced aspects of their locomotion. Compared to current theoretical models, our models offer some advantages for accurate predictions. First, our models include sperm heads, which are often missing in current models for attached or detached sperm cells (Cripe et al. Reference Cripe, Richfield and Simons2016; Simons & Rosenberger Reference Simons and Rosenberger2021; Carichino et al. Reference Carichino, Drumm and Olson2021). In addition, we consider both the mechanical and hydrodynamic cell–cell interactions, whereas only one type of interaction has typically been incorporated in current theoretical models (Pearce et al. Reference Pearce, Hoogerbrugge, Hook, Fisher and Giomi2018; Carichino et al. Reference Carichino, Drumm and Olson2021).
Our simulations show that the mechanical head–head coupling significantly affects the swimming trajectory of sperm pairs (figures 3 and 6). The head orientations
$\alpha _1$
and
$\alpha _2$
are history-dependent, while the displacement
$d$
and travelling distance
$S$
are not. The time-dependent trajectory of paired spermatozoa is sensitive to the initial values of their head orientation, angular velocity and translational velocity. A small variation in the initial values would lead to a significant cumulative error in the predicted trajectory. Comparing the predictions in the displacement during one beat cycle
$d$
with the experimentally observed ones, we find that both underpredictions and overpredictions appear in our simulations, depending on models and sperm pairs (figure 3). In a previous computational study for single sperm cells, only underpredictions in
$d$
were found (Gillies et al. Reference Gillies, Cannon, Green and Pacey2009). In our simulations, the discrepancy in the displacement
$d$
between the predictions and experimental observations likely results from the missing information in the environment and flagella. For instance, we postulate that the sperm pairs swim in the middle plane between the walls, which may not align with reality. Closer proximity of sperm cells to the wall surface results in higher swimming speed (Smith et al. Reference Smith, Gaffney, Blake and Kirkman-Brown2009). Additionally, although our mathematical models treat flagellar beating as planar, it is not strictly planar in reality, particularly at the distal end of the flagellum. A slight tilt of the planar motion out of the swimming plane may contribute to the discrepancy. Another potential contributor is the missing information on whether sperm pairs come into contact with the chamber walls transiently.
Correct information, e.g. the time-varying angular velocities of heads, improves the accuracy in predicting the trajectory of sperm pairs. However, the underlying mechanism of the mechanical head–head coupling is not fully understood, limiting our prediction accuracy. The trajectory of each sperm pair during one beat cycle differs from that during another beat cycle. Such seemingly stochastic swimming has also been observed experimentally in individual sperm cells and is attributed to the active phase and amplitude fluctuations of flagellar beating (Ma et al. Reference Ma, Klindt, Riedel-Kruse, Jülicher and Friedrich2014). In the case of paired spermatozoa, cell–cell interactions are another contributor. The combined effects of flagellar fluctuations and cell–cell interactions are also reflected in the radial distribution of the points around the closed loop (figure 4 b).
Three different modes of flagellar synchronisation were observed in our experiments, reminiscent of the similar flagellar beat patterns found in Chlamydomonas (Leptos et al. Reference Leptos, Wan, Polin, Tuval, Pesci and Goldstein2013; Wan et al. Reference Wan, Leptos and Goldstein2014; Quaranta et al. Reference Quaranta, Aubin-Tam and Tam2015; Wan & Goldstein Reference Wan and Goldstein2016) and model microfilaments (Guo et al. Reference Guo, Fauci, Shelley and Kanso2018,Reference Guo, Man, Wan and Kanso2021). Current models focus primarily on flagellar synchronisation at a rigid/elastic base, using either low-order representations of flagella as oscillators (Uchida & Golestanian Reference Uchida and Golestanian2011; Golestanian et al. Reference Golestanian, Yeomans and Uchida2011; Klindt et al. Reference Klindt, Ruloff, Wagner and Friedrich2017; Hickey et al. Reference Hickey, Golestanian and Vilfan2023) or more realistic representations as elastic beams (Riedel-Kruse et al. Reference Riedel-Kruse, Hilfinger, Howard and Jülicher2007; Osterman & Vilfan Reference Osterman and Vilfan2011; Tam & Hosoi Reference Tam and Hosoi2011; Goldstein et al. Reference Goldstein, Lauga, Pesci and Proctor2016; Guo et al. Reference Guo, Man, Wan and Kanso2021). In theoretical models for swimming flagella, their waveforms were typically prescribed as sinusoidal waves (Elfring et al. Reference Elfring, Pak and Lauga2010; Cripe et al. Reference Cripe, Richfield and Simons2016). In our models, flagellar waveforms are prescribed based on experimental observations, which is a more realistic representation but brings in some drawbacks, such as the exclusion of the phase-lag dependence of flagellar waveforms and the necessity for a priori information in flagellar beating.
The view prevails since Taylor (Reference Taylor1951) that co-swimming sperm cells tend to form in-phase synchronisation to minimise the energy dissipation in the fluid. However, our results do not fully support this perspective. While in-phase synchronisation leads to relatively low power consumption, the minimum power can occur at a small flagellar phase lag due to the mechanical and hydrodynamic interactions (figure 5). Similarly, a recent computational study on two-dimensional co-swimming sheets has demonstrated that minimum power can occur at any flagellar phase lag, depending on the specific kinematics of the sheets (Liao & Lauga Reference Liao and Lauga2021). Moreover, a recent study on swimming paired flagella has indicated that the anti-phase synchronisation can induce the highest swimming speed and efficiency (Cripe et al. Reference Cripe, Richfield and Simons2016), which is also found in our calculations. However, our simulations show that the highest swimming speed and efficiency of paired spermatozoa can result from flagellar phase lags besides
$\pi$
radians. To further confirm the existence of a statistical relationship between the swimming performances of paired spermatozoa and their flagellar beat patterns, a larger sample size is required.
Supplementary movie.
A supplementary movie is available at https://doi.org/10.1017/jfm.2025.79.
Acknowledgements.
K.Z. acknowledges financial support from the China Scholarship Council (CSC).
Funding.
This work was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Programme under grant 866494 project-MAESTRO.
Declaration of interests.
The authors declare no conflict of interest.
Author contributions.
K.Z. and I.S.M.K. designed the research; K.Z. developed the models, analysed the data, and drafted the manuscript; A.L. and V.M. conducted the experiments; K.Z., A.K., V.M., S.M. and I.S.M.K. revised the manuscript.
Appendix A. Spatial discretisation of a sperm cell
In our models, the flagellum is discretised with cross-sections equally spaced along its length with a distance equal to its radius
$\rho$
. The cross-section is a circle. Six regularised Stokeslets are equally spaced along each circular cross-section, as illustrated in figure 8(a). The unit vectors
$\boldsymbol {e}_1$
and
$\boldsymbol {e}_2$
respectively align with the two principal axes of the ellipsoidal head, as illustrated in figure 1. Let
$\boldsymbol {e}_3$
denote the unit vector aligned with the third principal axis of the ellipsoidal head, where
$\boldsymbol {e}_3=\boldsymbol {e}_1\times \boldsymbol {e}_2$
, as illustrated in figure 8(b). The sperm head is discretised with
$n_{\mathrm {h}}$
cross-sections along its long axis (i.e. the axis aligned with vector
$\boldsymbol {e}_1$
). These
$n_{\mathrm {h}}$
cross-sections are perpendicular to
$\boldsymbol {e}_1$
and divide the elliptical cross-section of the ellipsoid at the
$(\boldsymbol {e}_1,\boldsymbol {e}_3)$
plane into arcs. These arcs have the same arc length of approximate
$\rho$
by adjusting the distance between the
$n_{\mathrm {h}}$
cross-sections. In each of the
$n_{\mathrm {h}}$
cross-sections, regularised Stokeslets are equally spaced at distance approximately
$\rho$
. To discretise the walls, regularised Stokeslets are distributed on each wall in a square pattern with side length
$10\rho$
, as illustrated in figure 8(c).
Appendix B. Calculation of hydrodynamic forces on paired sperm cells
For the total
$n=n_1+n_2+n_{{\rm w}}$
regularised Stokeslets, the fluid response
$\boldsymbol {u}$
at
$^{m}\boldsymbol {r}$
to the hydrodynamic forces
$^k\boldsymbol {f}$
is calculated by
\begin{equation} \boldsymbol {u} (^{m}\boldsymbol {r} )=-\frac {1}{8\pi \mu }\sum _{k=1}^{n}\sum _{p=1}^{3}{S}^{\epsilon }_{pq} ({^k\boldsymbol {r}},{^{m}\boldsymbol {r}} )\,{^kT}\;{}^{k}_{p}{\boldsymbol {f}},\quad m=1, 2,\ldots ,n, \end{equation}
where
$^kT$
is the quadrature weight of the kth regularised Stokeslets,
$\epsilon$
is the regularisation parameter, and
$^{k}_{p}{\boldsymbol {f}}$
is the pth component of the force on the cells by the fluid (Cortez et al. Reference Cortez, Fauci and Medovikov2005).
Appendix C. Force balance and torque balance on each of the paired sperm cells
As the adhesive force and hydrodynamic forces on cell 1 are balanced, we have
\begin{equation} \begin{aligned} {\boldsymbol{C}} & = \left [\begin{matrix} ^{1}{\boldsymbol{C}}_1 & \quad ^{2}{\boldsymbol{C}}_1 & \quad \cdots & \quad ^{n_1}{\boldsymbol{C}}_1 & \quad {\boldsymbol{0}}_{1\times 3(n_2+n_{{\rm w}})}\end{matrix}\right ], \\ {\boldsymbol{E}} & = \left [\begin{matrix} ^{1}{\boldsymbol{E}}_1 & \quad ^{2}{\boldsymbol{E}}_1 & \quad \cdots & \quad ^{n_1}{\boldsymbol{E}}_1 & \quad {\boldsymbol{0}}_{1\times 3(n_2+n_{{\rm w}})} \end{matrix}\right ], \end{aligned} \quad \begin{aligned} {\boldsymbol{D}} & = \left [\begin{matrix} 0 & \quad 0 & \quad 0 & \quad 0 & \quad -1 & \quad 0 \end{matrix}\right ], \\ {\boldsymbol{F}} & \quad = \left [\begin{matrix} 0 & \quad 0 & \quad 0 & \quad 0 & \quad 0 & \quad -1 \end{matrix}\right ], \end{aligned} \end{equation}
where
$^{k}{\boldsymbol{C}}_1= [\begin{matrix} 1 & 0 & 0 \end{matrix} ]$
and
$^{k}{\boldsymbol{E}}_1 = [\begin{matrix} 0 & 1 & 0 \end{matrix} ]$
,
$k=1, 2,\ldots , n_1$
. The torques due to these forces on cell 1 are balanced about the origin of the laboratory frame, so we have
\begin{equation} {\boldsymbol{G}} =\left [\begin{matrix} ^{1}{\boldsymbol{G}}_1 & \quad ^{2}{\boldsymbol{G}}_1 & \quad \cdots & \quad ^{n_1}{\boldsymbol{G}}_1 & \quad {\boldsymbol{0}}_{1\times 3(n_2+n_{{\rm w}})}\end{matrix}\right ], \quad {\boldsymbol{H}} = \left [\begin{matrix} 0 & \quad 0 & \quad 0 & \quad 0 & \quad r_{{\rm p}y} & \quad -r_{{\rm p}x} \end{matrix}\right ], \end{equation}
where
$^{k}{\boldsymbol{G}}_1= [\begin{matrix} -{^kr_{1y}} & ^kr_{1x} & 0 \end{matrix} ]$
,
$k=1, 2,\ldots , n_1$
. Here,
$^kr_{1x}$
and
$^kr_{1y}$
are the components of the position vector of the kth point on cell 1 along the x- and y-axes of the laboratory frame, respectively. Likewise,
$r_{{\rm p}x}$
and
$r_{{\rm p}y}$
are the x- and y-components of
$\boldsymbol{r}_{\rm p}$
. Based on the force balance on cell 2, we have
\begin{equation} \begin{aligned} {\boldsymbol{I}} & = \left [\begin{matrix} {\boldsymbol{0}}_{1\times 3n_1} & \quad\!\! ^{1}{\boldsymbol{I}}_2 & \quad\!\! ^{2}{\boldsymbol{I}}_2 & \quad\!\! \cdots & \quad\!\! ^{n_2}{\boldsymbol{I}}_2 & \quad\!\! {\boldsymbol{0}}_{1\times 3n_{{\rm w}}}\end{matrix}\right ], \\ {\boldsymbol{K}} & = \left [\begin{matrix} {\boldsymbol{0}}_{1\times 3n_1} & \quad\!\! ^{1}{\boldsymbol{K}}_2 & \quad\!\! ^{2}{\boldsymbol{K}}_2 & \quad\!\! \cdots & \quad\!\! ^{n_2}{\boldsymbol{K}}_2 & \quad\!\! {\boldsymbol{0}}_{1\times 3n_{{\rm w}}}\end{matrix}\right ], \end{aligned} \quad\!\! \begin{aligned} {\boldsymbol{J}} & = \left [\begin{matrix} 0 & \quad\!\! 0 & \quad\!\! 0 & \quad\!\! 0 & \quad\!\! 1 & \quad\!\! 0 \end{matrix}\right ], \\ {\boldsymbol{L}} & = \left [\begin{matrix} 0 & \quad\!\! 0 & \quad\!\! 0 & \quad\!\! 0 & \quad\!\! 0 & \quad\!\! 1 \end{matrix}\right ], \end{aligned} \end{equation}
where
$^{k}{\boldsymbol{I}}_2= [\begin{matrix} 1 & 0 & 0 \end{matrix} ]$
and
$^{k}{\boldsymbol{K}}_2= [\begin{matrix} 0 & 1 & 0 \end{matrix} ]$
,
$k=1, 2,\ldots , n_2$
. Based on the torque balance on cell 2 about the origin of the laboratory frame, we have
\begin{equation} {{\boldsymbol{M}}=\left [\begin{matrix} {\boldsymbol{0}}_{1\times 3n_1} & \quad\!\! ^{1}{\boldsymbol{M}}_2 & \quad\!\! ^{2}{\boldsymbol{M}}_2 & \quad\!\! \cdots & \quad\!\! ^{n_2}{\boldsymbol{M}}_2 & \quad\!\! {\boldsymbol{0}}_{1\times 3n_{{\rm w}}}\end{matrix}\right ]},\quad\!\! {{\boldsymbol{N}} = \left [\begin{matrix} 0 & \quad\!\! 0 & \quad\!\! 0 & \quad\!\! 0 & \quad\!\! -r_{{\rm p}y} & \quad\!\! r_{{\rm p}x}\end{matrix}\right ]}, \end{equation}
where
$^{k}{\boldsymbol{M}}_2= [\begin{matrix} -{^kr_{2y}} & ^kr_{2x} & 0 \end{matrix} ]$
,
$k=1, 2,\ldots , n_2$
. Here,
$^kr_{2x}$
and
$^kr_{2y}$
are the x- and y-components of the position vector of the kth point on cell 2, respectively.
Spatial discretisation of the sperm (a) flagellum, (b) head and (c) wall. A portion of the flagellum and wall is shown.
Appendix D. Force balance and torque balance on a sperm pair
As the forces on a sperm pair are balanced, we have
\begin{equation} \begin{aligned} {\boldsymbol{C}}^{\prime} & = \left [\begin{matrix} ^{1}{\boldsymbol{C}}^{\prime}_1 & \quad ^{2}{\boldsymbol{C}}^{\prime}_1 & \quad \cdots & \quad ^{n_1}{\boldsymbol{C}}^{\prime}_1 & \quad ^{1}{\boldsymbol{C}}^{\prime}_2 & \quad ^{2}{\boldsymbol{C}}^{\prime}_2 & \quad \cdots & \quad ^{n_2}{\boldsymbol{C}}^{\prime}_2 & \quad {\boldsymbol{0}}_{1\times 3n_{{\rm w}}}\end{matrix}\right ], \\ {\boldsymbol{E}}^{\prime} & = \left [\begin{matrix} ^{1}{\boldsymbol{E}}^{\prime}_1 & \quad ^{2}{\boldsymbol{E}}^{\prime}_1 & \quad \cdots & \quad ^{n_1}{\boldsymbol{E}}^{\prime}_1 & \quad ^{1}{\boldsymbol{E}}^{\prime}_2 & \quad ^{2}{\boldsymbol{E}}^{\prime}_2 & \quad \cdots & \quad ^{n_2}{\boldsymbol{E}}^{\prime}_2 & \quad {\boldsymbol{0}}_{1\times 3n_{{\rm w}}} \end{matrix}\right ], \end{aligned} \quad \begin{aligned} {\boldsymbol{D}}^{\prime} & = {\boldsymbol{0}}_{1\times 3}, \\ {\boldsymbol{F}}^{\prime} & = {\boldsymbol{0}}_{1\times 3}, \end{aligned} \end{equation}
where
$^{k}{\boldsymbol{C}}^{\prime}_1 = [\begin{matrix} 1 & 0 & 0 \end{matrix} ]$
and
$^{k}{\boldsymbol{E}}^{\prime}_1 = [\begin{matrix} 0 & 1 & 0 \end{matrix} ]$
,
$k=1, 2,\ldots , n_1$
, with blocks
$^{k}{\boldsymbol{C}}^{\prime}_2 = [\begin{matrix} 1 & 0 & 0 \end{matrix} ]$
and
$^{k}{\boldsymbol{E}}^{\prime}_2 = [\begin{matrix} 0 & 1 & 0 \end{matrix} ]$
,
$k=1, 2,\ldots , n_2$
. The torques due to these forces on the sperm pair are balanced about the origin of the laboratory frame, so we have
\begin{equation} {\boldsymbol{G}}^{\prime} =\left [\begin{matrix} ^{1}{\boldsymbol{G}}^{\prime}_1 & \quad ^{2}{\boldsymbol{G}}^{\prime}_1 & \quad \cdots & \quad ^{n_1}{\boldsymbol{G}}^{\prime}_1 & \quad ^{1}{\boldsymbol{G}}^{\prime}_2 & \quad ^{2}{\boldsymbol{G}}^{\prime}_2 & \quad \cdots & \quad ^{n_2}{\boldsymbol{G}}^{\prime}_2 & \quad {\boldsymbol{0}}_{1\times 3n_{{\rm w}}}\end{matrix}\right ], \quad {\boldsymbol{H}}^{\prime} = {\boldsymbol{0}}_{1\times 3}, \end{equation}
where
$^{k}{\boldsymbol{G}}^{\prime}_1= [\begin{matrix} -{^kr_{1y}} & ^kr_{1x} & 0 \end{matrix} ]$
,
$k=1, 2,\ldots , n_1$
. Here,
$^kr_{1x}$
and
$^kr_{1y}$
are the components of the position vector of the kth point on cell 1 along the x- and y-axes of the laboratory frame, respectively, with blocks
$^{k}{\boldsymbol{G}}^{\prime}_2= [\begin{matrix} -{^kr_{2y}} & ^kr_{2x} & 0 \end{matrix} ]$
,
$k=1, 2,\ldots , n_2$
, where
$^kr_{2x}$
and
$^kr_{2y}$
are the components of the position vector of the kth point on cell 2 along the x- and y-axes of the laboratory frame, respectively.
Appendix E. Verification of the spatial discretisation of spermatozoa
Following the approach of Cortez et al. (Reference Cortez, Fauci and Medovikov2005) and Gillies et al. (Reference Gillies, Cannon, Green and Pacey2009), we choose an optimal regularised parameter
$\epsilon$
by comparing our predictions of the resistive force on two isolated ellipsoids translating in the x- and y-directions in an unbounded fluid with their exact solution. The first ellipsoid has the same dimensions as the sperm head (
$9\times5\times0.4$
$\unicode{x03BC}\textrm {m}$
), and is discretised using the same method as for the sperm head. The resistive force on a tube is often approximated using the slender body theory (Johnson Reference Johnson1980; Rodenborn et al. Reference Rodenborn, Chen, Swinney, Liu and Zhang2013). We represent the sperm flagellum with the second ellipsoid with an elongated needle-like geometry of
$50\times0.25\times0.25$
$\unicode {x03BC}\textrm {m}$
. is discretised using the same method as for the sperm flagellum. The non-dimensional resistive force error
$e_{{\rm f}}$
is defined as
\begin{equation} e_{{\rm f}} = \sqrt {\frac {1}{n_{{\rm c}}\,\zeta _{{\rm max}}^2}\sum _{i=1}^{n_{{\rm c}}}\zeta _i^2}, \end{equation}
with
$\zeta _i=\zeta _{\mathrm {exact}}-\zeta _{\mathrm {pred}}$
, where
$\zeta _{\mathrm {exact}}$
is the exact resistive force on the ellipsoids (Kim Reference Kim1986), and
$\zeta _{\mathrm {pred}}$
is our predicted resistive force on the ellipsoids. The number of the separate simulations is
$n_{\mathrm {c}}=4$
, i.e. the two ellipsoids translating in the x- and y-directions. Here,
$\zeta _{\mathrm {max}}$
is the maximum value of
$\zeta _i,\ i = 1,2,3,4$
. The force error is minimised when
$\epsilon =0.36\rho$
, as illustrated in figure 9(a). Using
$\epsilon =0.36\rho$
, we further calculate the velocity
$U_{z}$
of the points on an isolated sphere moving with a unit velocity along the z-axis in an unbounded fluid. This sphere has a volume equivalent to that of the sperm head, and is discretised into 682 points in the same manner as for the sperm head. Cumulative distribution of velocity error for the points
$e_{{\rm v}}=|U_{z}-1|$
shows that all the points have a velocity error less than
$0.018\,{\rm m\,s}^{-1}$
(figure 9
b). Note that the regularised parameter
$\epsilon =0.36\rho$
used here is to minimise the numerical error for the sperm head and flagellum moving in an unbounded fluid. Its value needs to be changed if walls are included.
(a) Resistive force error
$e_{{\rm f}}$
for two isolated ellipsoids is dependent on the regularised paratmeter
$\epsilon$
. (b) Cumulative distribution of the velocity error
$e_{{\rm v}}$
for the points on the sphere.
Appendix F. Verification of the regularised parameter
In our test cases,
$\xi =0.5$
and
$m=0.9$
are chosen, which work well. Due to the lack of established convergence analysis for the regularised Stokeslets method (Ainley et al. Reference Ainley, Durkin, Embid, Boindala and Cortez2008), we follow the approach of Ainley et al. (Reference Ainley, Durkin, Embid, Boindala and Cortez2008) and Gillies et al. (Reference Gillies, Cannon, Green and Pacey2009), and compare the predicted results in our test cases with published ones, as shown in figure 10. In our test cases, we calculate the resistive force
$\boldsymbol {F}_{\mathrm {s}}$
on a sphere translating in the x- and y-directions between the walls with velocity
$\boldsymbol {U}_{\mathrm {s}}$
. This sphere has a volume equivalent to that of the sperm head. The walls are perpendicular to the z-axis. The resistive force
$\boldsymbol {F}$
on the sphere translating in an unbounded fluid with the same velocity
$\boldsymbol {U}_{\mathrm {s}}$
is
$\boldsymbol {F}_{\mathrm {s}\infty }=6\pi \mu r_{\mathrm {s}}\boldsymbol {U}_{\mathrm {s}}$
, where
$\mu$
is the viscosity of the fluid. The ratio
$F_{\mathrm {s}}/F_{\mathrm {s}\infty }$
depends on the ratio
$h/r_{\mathrm {s}}$
. As shown in figure 10(a), our theoretical results agree with those of Ramia et al. (Reference Ramia, Tullock and Phan-Thien1993) and Gillies et al. (Reference Gillies, Cannon, Green and Pacey2009).
(a) Compared to the resistive force
$F_{\mathrm {s}\infty }$
on the sphere translating in an unbounded fluid, it experiences a larger resistive force
$F_{\mathrm {s}}$
when translating between the walls with the same velocity. (b) For the isolated tube parallel to the walls, its normal resistance coefficient
$C_{\mathrm {n}}$
is larger than its tangential resistance coefficient
$C_{\mathrm {t}}$
. The tube translating in an unbounded fluid has a smaller tangential resistance coefficient
$C_{\mathrm {t}\infty }$
than that when translating between the walls. Our predictions for (c)
$F_{\mathrm {s}}/F_{\mathrm {s}\infty }$
for the sphere, and (d)
$C_{\mathrm {n}}/C_{\mathrm {t}}$
and
$C_{\mathrm {t}}/C_{\mathrm {t}\infty }$
for the tube are insensitive to the wall size.
In addition, we calculate the normal resistance coefficient
$C_{\mathrm {n}}$
and tangential resistance coefficient
$C_{\mathrm {t}}$
of an isolated tube translating in the radial and longitudinal directions between the walls. The tube is parallel to the walls. The tube has dimensions identical to those of the sperm flagellum, and is discretised in the same manner as for the flagellum. The ratios
$C_{\mathrm {n}}/C_{\mathrm {t}}$
and
$C_{\mathrm {t}}/C_{\mathrm {t}\infty }$
depend on the ratio
$h/L$
, where
$C_{\mathrm {t}\infty }$
is the tangential resistance coefficient of an isolated tube translating in an unbounded fluid, as shown in figure 10(b). Our theoretical results agree with those of Ramia et al. (Reference Ramia, Tullock and Phan-Thien1993) and Gillies et al. (Reference Gillies, Cannon, Green and Pacey2009). Furthermore, as illustrated in figures 10(c,d), the wall size has only a slight impact on these ratios as long as the sphere and tube are covered by the walls. When the side length of the walls ranges from 60 to 160
$\unicode {x03BC}\textrm {m}$
, the maximum variation in
$F_{\mathrm {s}}/F_{\mathrm {s}\infty }$
of the sphere is
$0.6\,\%$
, and the maximum variations in
$C_{\mathrm {n}}/C_{\mathrm {t}}$
and
$C_{\mathrm {t}}/C_{\mathrm {t}\infty }$
of the tube are
$0.9\,\%$
and
$2.5\,\%$
, respectively. Note that several different values of the wall half-depth
$h$
are included in our test cases, but
$h=10\,\unicode {x03BC}\textrm {m}$
is kept fixed for all our further simulations, as the depth of the chamber that we used is 20
$\unicode {x03BC}\textrm {m}$
.
Appendix G. Reconstruction of flagellar shapes and their limit cycles
Using principal component analysis, we can approximate the tangent angle along each flagellum
$\psi (l,t)$
as
\begin{equation} \psi (l,t)\approx {^0\psi (l)+\beta _1(t)\,\boldsymbol {V}_1(l)+\beta _2(t)\,\boldsymbol {V}_2(l)}, \end{equation}
where
$^0\psi (l)$
is the mean tangent angle (Geyer et al. Reference Geyer, Jülicher, Howard and Friedrich2013; Werner et al. Reference Werner, Rink, Riedel-Kruse and Friedrich2014). Here,
$\boldsymbol {V}_1$
and
$\boldsymbol {V}_2$
are the first two flagellar shape modes, and
$\beta _1$
and
$\beta _2$
are their respective scores. This approximation accounts for
$95\,\%$
of the variance of the tangent angles. To reduce the noise of the tangent angles, we approximated the flagellar shapes using smoothing splines before calculating the shape modes.
The flagellar shapes at different times were thereby mapped to the
$(\beta _1,\beta _2)$
space, represented by the points in the space (figure 4
b). By fitting the points, we obtained the limit cycle. Flagellar phases were defined by binning flagellar shapes according to shape similarity. We then reconstructed flagellar tangent angles from the limit cycle. Finally, flagellar beat patterns with various phase lags were reconstructed from the tangent angles using (3.1).
Appendix H. Calculation of power consumption
After a sperm pair has established regular beating in our simulations, the instantaneous power dissipated in the fluid by the ith sperm cell at a time
$t_j$
is determined by
\begin{equation} P_i\left (t_j\right )=\sum _{k=1}^{n_i}\left ({^k\boldsymbol {f}_i\left (t_j\right )}\boldsymbol {\cdot }{^k\boldsymbol {u}_i\left (t_j\right )}\right ), \end{equation}
where
$n_i$
is the number of the regularised Stokeslets points on the ith sperm cell. The force
$^k\boldsymbol {f}_i$
and the velocity
$^k\boldsymbol {u}_i$
are determined from (4.6) or (4.9), depending on which model is used. The average power over one beat cycle for the ith cell is
$\langle P_i\rangle =\sum _{j=1}^{n_{\mathrm {t}}}P_i(t_j)/n_{\mathrm {t}}$
, where
$n_{\mathrm {t}}$
is the number of time steps in this beat cycle. The average power per cell for the sperm pair is
$\langle P\rangle =\sum _{i=1}^{2} \langle P_i\rangle /2$
.
Open access
