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

(3.1) \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

(4.1) \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

(4.2) \begin{equation} ^{k}\boldsymbol {u}_i(t)={^{k}\boldsymbol {v}_i(t)}+{^{k}{\boldsymbol{B}}_i(t)}\,\mathcal {U}(t), \end{equation}

where

(4.3) \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

(4.4) \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

(4.5) \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,

(4.6) \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

(4.7) \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

(4.8) \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 }}$ ,

(4.9) \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.

Algorithm 1. 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.

Figure 2. (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

Figure 3. 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

(4.10) \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

Figure 4. 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.

Figure 5. 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.

Figure 6. 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).

Figure 7. 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$ .