A Computational Study of Hydrodynamic Interactions Between Pairs of Sperm With Planar and Quasi-Planar Beat Forms

Although hydrodynamic interactions and cooperative swimming of mammalian sperm are observed, the key factors that lead to attraction or repulsion in different confined geometries are not well understood. In this study, we simulate the 3-dimensional fluid-structure interaction of pairs of swimmers utilizing the Method of Regularized Stokeslets, accounting for a nearby wall via a regularized image system. To investigate emergent trajectories of swimmers, we look at different preferred beat forms, planar or quasi-planar (helical with unequal radii). We also explored different initializations of swimmers in either the same plane (co-planar) or with centerlines in parallel planes. In free space, swimmers with quasi-planar beat forms and those with planar beat forms that are co-planar exhibit stable attraction. The swimmers reach a maintained minimum distance apart that is smaller than their initial distance apart. In contrast, for swimmers initialized in parallel beat planes with a planar beat form, we observe alternating periods of attraction and repulsion. When the pairs of swimmers are perpendicular to a nearby wall, for all cases considered, they approach the wall and reach a constant distance between swimmers. Interestingly, we observe sperm rolling in the case of swimmers with preferred planar beat forms that are initialized in parallel beat planes and near a wall.


INTRODUCTION
The tumultuous journey of the mammalian sperm involves navigating the female reproductive tract. Even though millions of sperm are deposited at the beginning of the tract, only a select few are able to traverse the long distances and overcome all of the hurdles to make it to the egg [1,2]. Using a single flagellum, sperm progress through a wide range of environments and their motility patterns must change in response to various chemical and physical cues; this could act to group sperm together or separate them [3][4][5]. A sperm senses other nearby sperm and surfaces via hydrodynamic interactions, which results in a wide range of collective motion, from alignment in trains and vortices to synchronization and attraction [6][7][8].
Sperm are able to propel themselves through different fluid environments by bending their elastic and flexible flagellum or tail. Within the flagellum, dynein is the molecular motor that actively generates force along the length of the flagellum, which results in a bending wave [2,[9][10][11]. The emergent beat form and trajectory will depend on the local fluid environment and the particular species of sperm. In experiments, emergent planar, helical, and quasi-planar flagellar beat forms have been observed, with planar beat forms more likely in higher viscosity fluids [10,11]. Similarly, tracking of sperm trajectories has shown linear trajectories as well as helical trajectories [12,13].
Surface interactions play an important role of many microorganisms, including sperm swimming in the reproductive tract [6,14]. Once sperm reach the oviduct, they generally bind to the epithelial cell wall and remain in this sperm reservoir until the time of ovulation [15]. At ovulation, signaling molecules such as heparin or progesterone will increase in concentration and will aid in the release of the sperm from the oviductal wall [15][16][17]. Likely, these signaling molecules will either act to break bonds of the cell body with the wall or bind to the flagellum and initiate a different beat form that will aid in generating increased force to break away from the wall [18,19]. There are many interesting and unanswered questions as to how sperm get to the walls, how long they stay at the walls, as well as how the sperm reservoir could act as a filter to sort out healthy sperm [20,21]. Investigations with microfluidic channels have revealed that sperm guidance can be achieved with surface topography and microchannels [22,23].
There are many different modeling approaches that can be taken to model sperm motility, but one of the key ingredients is how the sperm cell itself is being represented. Simplified approaches neglect the head or cell body [24][25][26]; headless sperm can still swim and previous analysis has shown that neglecting the cell body results in similar dynamics [27,28]. Other studies have focused on capturing accurate cell body geometries to investigate the role on swimming and interactions [29][30][31]. Similarly, when accounting for the bending of the flagellum, there are a few options. One can actuate or drive the dynamics of the flagellum with forces or torques corresponding to a curvature wave [24,32], prescribe a preferred curvature that is utilized in an energy functional that determines forces [25,26,33], or exactly prescribe the beat form [29][30][31]. When the beat form is prescribed exactly, the swimmer may rotate, but the flagellum will always have its entire centerline in the same plane and it will always have the same beat form parameters (e.g., amplitude and beat frequency). Preferred beat forms will have emergent beat forms based on interactions with other swimmers and/or surfaces.
To date, there have been many studies that have investigated the 3-dimensional (3D) dynamics of a single sperm near a wall [10,19,34]. Wall attraction was observed when approaching at specific angles, including perpendicular to the wall. Similarly, there have been studies of pairs of swimmers in free space to study the dynamics of attraction since biologically, most sperm are not swimming in isolation [24][25][26]31]. Since there has not been a detailed study of the interactions of pairs of swimmers near a wall, our goal is to further investigate and characterize how pairwise interactions vary from free space to the case where the sperm are swimming in close proximity and perpendicular to a planar wall. Using the method of regularized Stokeslets with an image system to account for the wall, we study swimmers propagating both preferred planar and quasi-planar beat forms. The dynamics of attraction and repulsion of swimmers is explored through different configurations; we vary both the initial distance between the sperm as well as the planes in which the swimmers are initialized. Our results show that sperm will attract to and stay near the wall while phenomenon such as sperm rolling will occur for a subset of sperm configurations. Our results further contextualize divergent results for pairs of swimmers in previous studies and provides insight into relevant interactions that can be utilized in the development of artificial micro-swimmers [35][36][37].

Mathematical Model
We utilize a fluid-structure interaction framework where the sperm is assumed to be neutrally buoyant and immersed in the fluid. Since sperm swim in viscosity dominated environments and are often in close proximity to a wall or boundary [2,38], we will consider the 3-dimensional (3D) incompressible Stokes equations for the fluid velocity u above a planar wall at x W: subject to the boundary condition u(x)| x W 0. Here, μ is the viscosity, p is the pressure, and f b is the sperm body force density exerted on the fluid as the flagellum actively bends. That is, the hydrodynamic stresses are coupled to the bending of the flagella; the fluid "feels" the swimmers through f b and the swimmers "feel" each other through the fluid.
To determine the body forces, we will study a simplified representation of a sperm cell where we neglect the cell body, similar to previous studies of [32,39]. Additionally, we assume the sperm flagellum is isotropic and homogeneous with constant radius much smaller than length L so that we can represent it using the Kirchhoff Rod (KR) model. This allows for each cylindrical elastic flagellum ι to be represented via a space curve X ι (t, q) for arc length parameter q (0 ≤ q ≤ L). Local twisting is accounted for via a right-handed orthonormal triad, We briefly summarize the model here and refer the reader to [39][40][41][42] for additional details.
In the KR model, through a variational energy argument, we can define the internal force F ι and torque M ι on a cross section of flagellum ι as F ι respectively. The scalar components of the force and torque in the Darboux frame are defined as: where δ 3i is the Kronecker delta and F i ι , M i ι are defined for any cyclic permutation of (i, j, k). The physical properties of the flagellum determine the moduli (a 1 , a 2 -bending, a 3 -twisting, b 1 , b 2 -shearing, b 3 -stretching); the stiffness of sperm flagella can be measured experimentally and related to these moduli [43,44]. The preferred kinematics can be described by the preferred twist κ 3 , normal curvatureκ 1 , and geodesic curvatureκ 2 [45,46]. In this formulation, there will be no torque when the flagellum is in its preferred shape, which we will now define.
Sperm propagate bending along the length of the flagellum due to the coordinated action of molecular motors inside the flagellum [2,47]. We can capture a range of beat forms representative of experimental observations [9,10] by assuming that the flagella have the preferred configuration of with beat frequency f ω/2π, wavelength 2π/η, and beat amplitudes α and β. Planar bending of the flagellum occurs when either α or β are zero whereas quasi-planar bending occurs when α, β ≠ 0 and either α < β or β < α. We note that we are not prescribing a rotation of the flagellum; we simply actuate or drive it to beat (and bend) in this preferred planar or quasiplanar configuration that is a function of arc length parameter q and time t. Similarly, these are the initialized and preferred configurations that can deviate later in the simulation due to interactions with other swimmers and/or the wall. Utilizing this preferred configuration, we can calculate the preferred orthonormal triadD(t, q) by settingD 3 as the tangent toX, and choosingD 1 as the normal andD 2 as the binormal vector.
Then, the preferred curvature and twist can be calculated aŝ which is a spatiotemporal function based onX(q, t) in Eq. 3.
For this configuration, we impose force and torque balance along the length of the flagellum. That is, the fluid feels the sperm through a force per unit length f ι and torque per unit length m ι [only defined on X ι (q, t)], The body force f b ι in Eq. 1 will be equal and opposite to the internal forces of the flagellum and will also need to be a smooth force field on the fluid domain with finite velocities on X. Hence, we will set where Γ ι X ι (q, t), x is any point in the 3D domain with x > W and r x − X . The regularization functions ψ ε and ϕ ε are radially symmetric, satisfying R 3 ψ ε (r)dx R 3 ϕ ε (r)dx 1, and smooth the singular force field with parameter ε governing the region where most of the force is spread [48,49]. In this model, each sperm flagellum will have a preferred configurationX and D(q, t), but its ability to achieve that configuration will depend on the local fluid flow through a no-slip condition, where w 1 2 ∇ × u is the angular velocity.

Numerical Method
We discretize each of the m flagella into p points and use standard second order central finite difference approximations to determine forces and torques in Eqs 2, 5. The body force in Eq. 6 is approximated using the trapezoidal rule. In free-space, utilizing the linearity of Eq. 1, the flow can be written as a superposition of regularized fundamental solutions (regularized due to the smoothing of the body force in Eq. 6). Since we are interested in the flow above a wall at x W, we utilize a regularized image system that first accounts for the resulting flow in R 3 (no wall) and then cancels the flow at the wall x W through combinations of fundamental solutions at additional image points on the other side of the wall (outside of the half space of interest). Let X ℓ ι (x ℓ ι , y ℓ ι , z ℓ ι ) be the ℓ-th discretized point on the ι-th flagellum, h ℓ ι x ℓ ι − W the height above the wall, and X im ℓι (−h ℓι , y ℓι , z ℓι ) the image point of X ℓι . The linear velocity u is where r ℓι x − X ℓι , r im ℓ ι x − X im ℓ ι , ξ n,ℓ ι −ξm ℓι , ξ f,ℓ ι −ξf ℓ ι , and ξ is the quadrature weight. S and R denote the fundamental solutions for a point force and point torque, corresponding to the Stokeslet and Rotlet, respectively, and the subindex refers to the particular regularization function that is being utilized to smoothly spread the force or torque in Eq. 6. S im ψ ε ,ϕ ε and R im ψ ε ,ϕ ε are the regularized image systems for the Stokeslet and Rotlet. Via a direct calculation, w 1 2 ∇ × u for u in Eq. 8. As detailed previously in [50,51]; [34], we set the regularization functions as which ensures that the boundary condition at the wall is satisfied, u| x W 0. Additionally, this ensures that the regularized image system in Eq. 8 satisfies the property that in the limit as ε → 0, the singular image system is recovered.
To solve this coupled fluid-structure interaction, at each time point t, we complete the above steps to determine the resulting linear and angular velocity due to a given flagellar configuration (assuming the preferred configuration in Eqs 3,4). Once the velocity of the flagellum, u (X ι ), and angular velocity of the triads, w (D ι ), are known, we determine the new flagellar location and triads at time t + τ by solving the no-slip equations in Eq. 7 with the forward Euler method. This process is repeated, solving for the new instantaneous flow due to each time-dependent force balance equation, which depends on the emergent flagellar configurations. We emphasize that this is a Lagrangian method; only the flagella are discretized. Through the regularized image system in Eq. 8, we automatically satisfy that the wall at x W has zero velocity. Additionally, we can evaluate Eq. 8 at any point of interest in the fluid domain.

Quantifying Interactions of Swimmers
The actual beat form, swimming speed, and trajectory that the sperm achieves is an emergent property of the coupled system that depends Frontiers in Physics | www.frontiersin.org November 2021 | Volume 9 | Article 735438 on both the geometry of the swimmers and the wall, fluid parameters such as viscosity, beat form parameters, as well as stiffness or moduli of the flagella. We will explore these dynamic and nonlinear relations through simulations. Here we highlight some of the metrics and different parameters that are studied. The sperm number is a non-dimensional number that characterizes the ratio of viscous fluid effects to elastic effects of the bending flagellum, computed as Here, L is the flagellum length, f is the beat frequency (inverse time), ε is the regularization parameter that approximates the flagellum radius, μ is fluid viscosity, and a 1 is the bending modulus. To compare to previous studies, we utilize the resistive force theory coefficient χ to capture viscous effects, similar to [32]. For mammalian sperm, previous estimates have had S p in the range of 2-17 [2]. Using the values reported in Table 1, the baseline sperm number for our simulations is S p 5.11. We explore increasing S p to a value of 7.64 by increasing μ to 5μ and increasing S p to 9.08 by increasing μ to 10μ. Similarly, we also explore increasing S p by setting a i to a i /5 or a i /10 for i 1, 2, 3. Even though this leads to the same sperm number, we explore the differences since a change in μ changes the magnitude of all terms in Eq. 8 whereas a change in a i only changes the magnitude of the terms involving the torque, which can be seen from Eq. 2.
Through interactions, the beating planes of swimmers may continue to deviate from the plane they were initialized in. A diagram of the directions of pitching and/or rolling out of the beating plane is given in Figure 1. The beating plane is calculated as the plane that passes through the center of mass X 1 p p ℓ 1 X ℓ of the swimmer and minimizes the orthogonal distances between the points on the flagellum and the plane [26,52]. The pitching angle θ of the swimmer's beating plane is the angle between the plane and the unitary vector in the x-direction e x , where n is the normal vector to the beating plane. The pitching angle θ is in the range [ −90°, 90°] and θ > 0 (θ < 0) corresponds to the plane pitching upward (downward). The rolling angle c of the swimmer's beating plane is the angle between the plane and the unitary vector in the y-direction e y , c 90°− arccos n · e y n .
The rolling angle c is also in the range [ −90°, 90°] and c > 0 (c < 0) corresponds to the plane rolling right (left) with respect to the direction of motion.

RESULTS AND DISCUSSION
The aim of this study is to further quantify and understand the hydrodynamic interactions of pairs of swimmers close to or far away from a wall. To do this, we will consider a few different scenarios that include planar or quasi-planar beat forms initialized a distance d apart. For the case of planar beat forms, we also consider the case of flagellar beating with centerlines initialized in parallel planes or the same plane (coplanar). The wall is either initialized at x −5 (near a wall) or at x −10,000, which we denote as the free space solution since the wall has negligible effects on swimming at this location. The baseline parameters used for the numerical methods and preferred beat form are given in Table 1; we assume that in a given simulation, all sperm flagella have the same preferred configuration given by Eqs 3, 4, with values in the range of mammalian sperm [10]. Swimmers are separated by a distance d in either the y or z-directions. We explore distances d in the range of 3-30 μm, i.e., distances where there is non-negligible flow effects from the nearby swimmer. Here, 30 μm is half the length of the swimmer (L 60 μm) and 3 μm is equal to the beat amplitude α.

Co-Planar Swimmers
This case involves two planar swimmers (α 3 and β 0 in Eq. 3) initialized such that their initial beating planes are in the plane z 0, shifted a distance d apart in the y-direction. That is, the average y-value, y, is set at y 0 for the bottom swimmer and y d for the top swimmer. Figure 2A shows the configuration of two co-planar sperm cells at time t 0 (gray dashed lines) and at time t 0.6 s (gray solid lines) for swimmers initialized d 6 μm apart. The first point is represented by a solid gray circle and denotes the cell body and swimming direction. The trajectories of the first point during this time frame (blue and black solid lines) show the hydrodynamic attraction of the two sperm cells; the top sperm starts to go down and attract to the bottom sperm, which is swimming with an upward trajectory. We classify this movement as yaw since these motions remain within the plane z 0, which has been observed in other modeling work for a single swimmer [29,30,34]. A solo swimmer has an upward yaw [34], but attraction dominates and TABLE 1 | Computational parameters for preferred planar (P) and quasi-planar (Q) beat forms.
Figures 2B,C track the minimum distance between the two sperm (solid lines) and the distance between the first points of the two sperm (circles) as a function of time, for the initial distances d considered. In both cases, for d 6 and d 30 μm, we observe a monotonic decrease in the distance between the sperm for this initial time period of attraction, similar to results of [26]. Even though the sperm are initialized at the same distance d apart along the entire flagellum, the minimum distance between the sperm cells occurs at the head (first point). For longer time simulations, we continue to observe attraction, reaching a steady configuration of flagella that are a very small distance apart. We do not show simulations or report this distance since a repulsion term is required to keep the filaments from occupying the same space and hence, the steady state distance between the flagella will depend on the strength and form of the repulsion term.
We also varied the sperm number S p (given in Eq. 10); similar dynamics of attraction were observed but on a slower time scale (results not shown). As previously noted in [32], for a solo swimmer, increasing S p results in a decreased swimming speed. We observe a similar decrease for a pair of swimmers with increased S p , which leads to the increased time scale for attraction to occur.
Previous studies with planar beat forms required to stay in the plane have reported an increase in swimming speeds for pairs of swimmers relative to the case of a solo swimmer when filaments propagate planar beat forms, are co-planar, and the flagella are sufficiently stiff [25]. Other studies that allow beat forms to deviate from the plane have reported a slowdown of swimming speed while swimmers attract [26]. For the parameters utilized in Figure 2 (reported in Table 1), we observe a decrease in swimming speed for a pair of swimmers initialized d 6 μm apart and a small increase in swimming speed for a pair of swimmers initialized d 30 μm. The speeds reported in Table 2 are looking at the magnitude of velocity of the first  Frontiers in Physics | www.frontiersin.org November 2021 | Volume 9 | Article 735438 5 point on the flagellum by utilizing locations at time t 0 and t 0.6 s; for two swimmers, the speed is the average over both swimmers. We do not observe a significant increase in swimming speed during the attraction phase for a range of initial separation distances. However, on longer time intervals after attraction has occurred, we do observe an increase in swimming speed (results not shown).

Near a Wall
We have shown that our model can capture a range of known phenomena for a pair of swimmers in free space, but sperm are not swimming in isolation and they are swimming in close proximity to walls as they navigate the female reproductive tract [2,5]. Hence, we wish to go beyond previous studies of a pair of swimmers in free space or a solo swimmer near a wall to understand whether the wall will help or hinder attraction of swimmers. For a solo swimmer, simulations have shown attraction to a wall is immediate but no yawing motion or vertical translation (up or down with respect to the centerline of the flagellum) is observed [34]. We note that for a solo swimmer initialized perpendicular to the wall, it will always attract to the wall regardless of whether it starts 2 μm or 50 μm away (the wall does not cause the swimmer to tilt and escape the wall). The questions are now 1) will additional sperm change the dynamics of attraction to a wall; 2) will swimmers near a wall still be able to attract to each other; 3) if attraction occurs,  Frontiers in Physics | www.frontiersin.org November 2021 | Volume 9 | Article 735438 6 will it happen with heads decreasing their distance apart at a faster rate (similar to Figures 2B,C)?
To study these questions, we consider the case of two coplanar swimmers initialized at a distance d apart and perpendicular to a planar wall at x −5 μm. Figure 3 shows the configurations of the swimmers at time t 0 (gray dashed lines) and at three snapshots in time (gray solid lines), for the initial distance of d 3 μm in A, d 6 μm in B, and d 11 μm in C. The first point of the swimmer is represented by a solid gray circle and the trace of this point (from t 0 until the time point given in that panel) is shown with the solid lines. For all values of d considered, the swimmers start with the head 5 μm away from the wall and rapidly attract to the wall, remaining close to the wall for the full time of the simulation (over 2.5 s).
In order to explore the rich dynamics of swimmers near a wall, we also report the minimum (solid line) and maximum (dotted line) distance between the two swimmers, along with the average distance between the first points of the two swimmers (circles) in Figure 3, for initial distance d 3 μm in D, d 6 μm in E and d 11 μm in F. The snapshots in time in A-C are identified with vertical gray lines in the corresponding distance graphs in D-F, and were chosen to highlight specific points of interest. The top panel in Figures 3A-C is the last time point for which the distance between the first points equals the minimum distance between the swimmers, i.e., the end of the first attraction period between the two swimmers where the heads are moving closer. The middle panel in Figures 3A-C is where the maximum distance between the swimmers is at a maximum. The stable configuration achieved by the swimmers is shown in the bottom panel in Figures 3A-C. Here, we are defining stable as attaining an average distance between points on the swimmer that persists in time. The same preferred beat form is given for all swimmers in all of these simulations; the presence of the wall and the initial separation distance is what causes the different beat forms to emerge. We note that two swimmers in free space attract, with the top swimmer yawing down and the bottom swimmer yawing up ( Figure 2B). The dynamics of the nearby wall prevent the swimmers to attract with equal yawing due to the emergent flagellar beat forms.
For all the values of d considered, the stable configuration achieved toward the end of the simulation (third snapshot) consists of an average first point distance of approximately 5.5 μm, in between the maximum and minimum average distance between the swimmers. As shown in Figures 3A,D, for an initial distance d smaller than 6 μm, the swimmer's first points initially attract, reaching a minimum distance on the order of 1 μm at t 0.34 s, then repulse reaching the maximum average distance of approximately 8 μm at t 0.62 s (with the top swimmer moving up). Due to hydrodynamic interactions of the swimmer and the close wall, even though the swimmers have centerlines that are parallel, we observe a dramatic yaw in the top swimmer as the head of the top swimmer is pushed up. Later, the swimmers reach a stable configuration after 2.50 s where flagellar centerlines are again parallel. For an initial distance d equal to 6 μm, shown in Figures 3B,E, the head of the swimmers (first point) attract up until t 0.33 s, reaching a distance of approximately 4.5 μm. After this initial attraction, the heads of the swimmer then repulse, reaching a maximum distance 9 μm apart at t 0.71 s, and reach the stable configuration after 2.50 s. For an initial distance d greater than 6 μm, Figures 3C,F, the maximum distance of 14 μm between the swimmers is reached at t 0.86 s and the distance between the first points show a continuous decrease in time until reaching a plateau and the stable configuration after 5.69 s. In this last scenario, the top swimmer, after reaching the wall, is moving downward to get closer to the bottom swimmer, clearly visible in Figure 3C at t 5.69 s. The asymmetry between the swimmer's behavior is due to the direction of motion chosen. When the direction of motion is reversed (wall at x 5, swimmers reflected about the y-axis, and preferred beat form propagating a wave in the opposite direction with − ω), we observe the bottom swimmer moving upward to get closer to the top swimmer (results not shown).
To explore competing effects of wall attraction and swimmer dynamics, we vary the sperm number S p in Eq. 10. Figure 4 shows the swimmer configurations obtained for the different values of S p considered in the case of two co-planar swimmers near a wall, varying μ in A,B,E and varying a i in C,D,F. Increasing S p means that the viscous effects are increasing relative to the elastic effect. As a result, in both Figures 4A-D, increasing S p exhibits a beat form with decreased achieved amplitude. In A,B, this is due to the increased viscosity of the fluid creating additional drag on the bending flagellum whereas in C,D, this is due to the decrease in the bending moduli making the flagellum less stiff and less able to propagate bending at the preferred amplitude. Even though the achieved amplitude of the flagellar beat form decreases, this does not prevent or slow down attraction to the wall (shown in zoomed views in B and D). Through a close examination of the distance between the swimmers at different time points, we can see the subtle differences between increasing S p via increasing μ in Figure 4E and decreasing a i in Figure 4F. The baseline value of S p 5.11 exhibits attraction and then repulsion with the heads being further apart than the rest of the flagellum at later time points (Figures 4E,F) and reaching a steady state configuration by t 2 s. In contrast, if an increase in S p to 7.64 or 9.08 is obtained by increasing μ by a factor 5 or 10, the swimmers do not maintain a stable configuration near the wall, but they continue to attract with the heads being the closest points at later time points ( Figure 4E). Increased S p in Figure 4E results in continued attraction whereas increased S p in Figure 4F results in a quasi steady state configuration at t 2 s.
We emphasize that for all simulations presented here for coplanar swimmers near a wall, the swimmers had a preferred planar configuration. With swimmer interactions and the presence of the wall, there was some out of plane motion, but on average, the swimmers tended to maintain a beat form in the same plane (due to torques in Eq. 2 that penalize deviations from the preferred motion). Additionally, for all simulations, the head or first point of the swimmers attract to the wall and maintains a small distance away from the wall. This occurs for a range of parameters in the case of two swimmers and also occurs on a similar time scale to that of a solo swimmer [34]. For a sperm number of S p 5.11, in the range for mammalian sperm, we observe swimmers reaching a somewhat steady state distance of Frontiers in Physics | www.frontiersin.org November 2021 | Volume 9 | Article 735438 ∼5.5 μm within a few seconds when starting in the same plane at a distance of 3-11 μm apart. Since the beat frequency of the swimmers is set to f 20 Hz; attraction on the order of seconds near the wall requires hundreds of beats of the flagellum. The dynamics of attraction will depend on the sperm number and increased S p (increased viscous forces) results in sperm being able to attract closer at the first point or cell body. This is important to consider since the epithelial cells on the oviductal walls may be secreting proteins and/or fluids that will change the viscosity near the walls [54] and potentially control or dominate emergent interactions and motility of sperm [55]. represents the wall. The minimum distance between the two swimmers (solid lines) and distance between the first points of the two swimmers (circles) are shown as averages over a beat period, varying μ in (E) and a i in (F).

Fluid Mixing
It is well known that as mammalian sperm navigate through the female reproductive tract, their motility patterns change in response to different signaling molecules, which may be released from the wall or be present in the entire fluid [2,3,5]. When sperm are trapped on the oviductal walls, signaling molecules such as heparin or progesterone initiate changes in motility by binding to the cell body or the flagellum [15,16]. In turn, changes to asymmetrical beating occur, which aids in the release of these trapped sperm [15][16][17]19]. The ability for these molecules to reach the sperm and bind to receptors is highly dependent on the local fluid flows. In Figure 5, we look at the mixing of the fluid by the swimmers in the case of an initial separation of d 3 μm in A,B and a separation of d 11 μm in C,D. The left hand side (A,C) is the case of the wall at x −5 μm where as the right hand side is the free space case (B,D). For each plot, the initialization is a plane of passive fluid markers (in the plane of the swimmers at t 0), similar to the inset in A. We initialize these locations and track their movement by solving Eq. 8 for their velocity due to the moving swimmers and update the fluid marker locations using Eq. 7. Similar to the swimmers in Figures 2, 3, the fluid particles initialized in this plane remain mostly in the plane. In all of the cases in Figure 5, we observe a region at the end of the flagella that is empty. These fluid particles have traveled with the swimmer and can be seen along the flagella. Due to the beat form of the flagella, we observe fluid markers that were originally a vertical stripe that have moved and now taken on a range of x-values. For example, the yellow fluid particles were initialized at x ∼ 9-10 μm and are now in the approximate range of x −5 to 20 μm in D. We observe significant movement in both directions for this x − location; this means that signaling molecules can get close to the flagella if being passively advected by the flow that the flagella are generating.

Parallel Swimmers
We continue to consider the case of two planar swimmers (α 3 and β 0 in Eq. 3), but now these swimmers are initialized with parallel beating planes. The distance d between the swimmers is initialized by placing the bottom swimmer in the plane z 0 and the top swimmer in the plane z d. However, the emergent beat plane may pitch upward or downward in the z-direction and/or roll left or right around the lateral axis (refer to Figure 1 for a schematic). Similar to the co-planar case, we wish to first benchmark our model and further explore the case of swimmers far away (free space) or close to a wall to understand how these dynamics change.

Free Space
The results in the case of free space are reported in Figure 6. The short-term simulations, up to t 3.25 s, are shown in Figures  6A-D. A non-monotonic behavior depending on the initial distance d between the swimmer's beating planes is observed. For d 5 μm, the distance between the first points of the two swimmers increases in time ( Figure 6A) while the rest of the swimmers slightly attract and then slowly increase their distance.
The beating plane of the bottom swimmer pitches downward while the top swimmer pitches upward by angles of ±2.3°at t 3 s ( Figures 6B,C), corresponding to the swimmers pushing away from each other. In contrast, for d 30 μm, the distance between the first points of the two swimmers decrease in time ( Figure 6A), attracting with the beating plane of the bottom swimmer pitching upward and the top swimmer pitching downward with angles of ±1.6° (Figures 6B,D). In the shorter time simulations ( Figures  6A-D), the beating plane of both swimmers have a minimal rolling motion, alternating between left and right rolling with −0.26°≤ c ≤ 0.26°for d 5 μm and −0.14°≤ c ≤ 0.14°for d 30 μm.
In terms of the swimming speeds for these shorter term dynamics over 3 s, the swimming speed of two parallel swimmers was similar to the corresponding solo swimmer ( Table 2). This is similar to previous results for pairs of swimmers separated by a distance of at least half their length, swimming speeds are similar to that of a solo swimmer [24,26,31]. However, we observe marked differences between swimmers that are co-planar and those that are in parallel planes. At an initial separation distance of d 6 μm, the swimmers in parallel planes are significantly faster (∼28 μm/s) than the case of the coplanar swimmers (∼20 μm/s).
To investigate the long-term dynamics, we look at an initial separation of d 3 μm for t 0-11 s in Figures 6E-G. The zoomed in portions of Figure 6E show the top and bottom swimmer's beating planes and corresponding pitching angles θ for the three snapshots in time, delineating the switches among near-field, mid-field and far-field dynamics. The times of the snapshot are identified with vertical gray lines in Figures 6F,G. The swimmers show near-field repulsion until t 3.5 s, i.e., with heads or first points increasing in separation ( Figure 6F) and beating planes pitching away from each other ( Figure 6G and right-most zoomed portion of Figure 6E). The top and bottom swimmer's beating planes obtain their maximum pitching angle at t 3.5 s and after t 3.5 s, the swimmers enter the mid-field regime where they will continue to repel each other ( Figure 6F) but the pitching angles will decrease in magnitude ( Figure 6G) and reach θ ≃ 0 at t 9.3 s (central zoomed portion of Figure 6E). After t 9.3 s, the swimmers show far-field attraction, i.e., decreasing distance between the swimmers ( Figure 6F) with beating planes pitching toward each other ( Figure 6G and leftmost zoomed portion of Figure 6E). In the long-term simulations, the rolling of the beating planes is also minimal, with −0.38°≤ c ≤ 0.38°over 11 s.
In summary, swimmers that are close to each other will initially show near-field repulsion and then eventually, after reaching a significant distance between each other, will transition to far-field attraction ( Figure 6E). Conversely, if the swimmers are initialized relatively far away from each other, the swimmers will initially show far-field attraction and then eventually, when getting too close to each other, will transition to near-field repulsion (results not shown for d 14 μm). Hence, dynamics of swimmers with planar beat forms in initially parallel beating planes will not reach a stable configuration of attraction and will continue to oscillate between attraction and repulsion. Our far-field attraction results differ from previous results of [26] where only repulsion was observed and [31] where only attraction is observed; this is likely due to different modeling assumptions with regards to the preferred planar beat form, how out of plane beating is penalized, and geometry of the cell body. We note that rotations of swimmers with respect to θ and c are also on par with previous studies [24].

Near a Wall
Similar to the previous cases, we wish to understand whether pairs of swimmers initialized in parallel planes will attract or repulse when near a wall. Results for the case of a wall at x −5 μm are highlighted in Figure 7. For all the values of d considered, the swimmers also attract to the wall (similar to the case of initially co-planar swimmers in Figures 3, 4). When the swimmers are initialized d 3 μm apart, they push apart and then quickly reach a constant distance apart ( Figure 7A) whereas in free space, they continued to push apart initially ( Figure 6A) and then oscillate between attraction and repulsion in the long term ( Figure 6E). With the wall, the beating plane of the bottom swimmer pitches downward while the top swimmer pitches upward, but at a very small angle (Figures 7B,C). For the full ∼ 4 s simulation, the beating planes show minimal rolling behavior with −0.91°≤ c ≤ 0.91°. For the case of an initial distance of d 30 μm, we observe constant attraction in Figure 7A, similar to the free space case in Figure 6A. In Figure 7A, the average distance between the two swimmers decreases monotonically in time for t ≥ 2.5 s. To understand why the first point or head distance is in between the minimum and maximum distance, we can see in Figure 7D that the entire length of the flagellum is not remaining in the same plane. The beating plane of the bottom swimmer pitches upward while the top swimmer pitches downward (i.e., toward each other), with a more preeminent rolling motion of the swimmer's beating planes ( Figures 7B,D), with −3.29°≤ c ≤ 3.29°for d 30 μm.
To investigate the long-term behavior, we report in Figures  7E-I the results for d 11 μm for ∼ 10 s. Figure 7E shows the Frontiers in Physics | www.frontiersin.org November 2021 | Volume 9 | Article 735438  Figure 7F. When initialized at d 11 μm, the swimmers show farfield attraction until t ≃ 4 s, reaching an average distance apart of ∼5 μm ( Figure 7F top panel). The top and bottom swimmer's beating planes obtain their maximum average pitching angle at t 1.2 s ( Figure 7G) and after t ≃ 4 s, the swimmers enter in the nearfield stability regime with the swimmer's beating planes pitching away from each other ( Figure 7F mid panel). However, after t ≃ 8 s, the swimmer's dynamics change drastically since the average rolling angle c, for both swimmers, increases and reaches the maximum value of ∼84°( Figure 7F bottom panel), i.e., both swimmer's beating planes roll to the right and the beating planes become almost perpendicular to the z 0 plane ( Figure 7H). Then, both swimmer's beating planes roll to the left with an almost 180°m otion ( Figure 7I), the average rolling angle c decreases and reaches the minimum value of ∼ −86°( Figure 7F bottom panel). After this second rotation, the swimmers reach a configuration similar to the one obtained in the two co-planar swimmers case in Figures 3C,G for the same initial distance d 11 μm.
In summary, if the swimmers are initialized close to a wall and relatively far away from each other, the swimmers will initially attract and then eventually, when getting too close to each other, will transition to a short-term near-field stability. After a certain period of time, the stability is broken by variations in the rolling angle c that cause the swimmer's planes to rotate (twice) and reach a final configuration in which the swimmers are almost coplanar, with beating planes almost perpendicular to the z 0 plane.
The rolling of sperm has been observed in experiments where the frequency of rolling is correlated with the beat frequency of the flagellum [10,56]. In this longer term simulation, we observe two rotations in ∼10 s with a beat frequency of 20 Hz ( Table 2), so this is at a higher rate than the beat frequency. Simulations observe rolling with a very low frequency but we hypothesize that additional perturbations to the flow from additional swimmers would increase the rolling rate; this is backed up by a recent study that showed a nonplanar component of the beat form is necessary to see rolling [57]. Indeed, rolling was previously observed in free space with a pair of swimmers when they were initialized as a perturbation to the coplanar configuration [26]. This will be important to further investigate as it has been proposed that sperm rolling plays an important role in selection of sperm as well as in the organization of sperm in the female reproductive tract [56]. In our simulations, the rolling episode is what enables the swimmers to fully align, allowing for cooperative movement of sperm swimming in close proximity and near a wall.

Particle Trajectories
We also investigated trajectories of passive fluid markers to further understand how signaling molecules near or around the swimmer would be advected by the flow. In Figures 8A,B, we look at particles initialized at different x-locations with swimmers separated by distances of d 3, 5 μm. For reference, even though the swimmer is not shown, it is similar to that of Figure 6 where the swimmers are 60 μm in length and the top swimmer is in the plane z d and the bottom swimmer is in the plane z 0. In Figure 8A, the fluid particles start below the bottom swimmer (z −2.9) and mid-way along the swimmer in the x-direction. We observe that for both cases, trajectories below the bottom swimmer are the same at this location and that particles are being pushed down and further back, similar to the pitching angles of the swimmers. Signaling molecules initialized in this region will not be able to reach and bind to the flagellum. When looking at a particle trajectory initialized slightly above the plane z 3 ( Figure 8B), we observe that the particles are moving in an upward trajectory and are being advected in the negative x-direction, corresponding to the direction of swimming. At these time points, the fluid particle is also pitching at a similar angle to that of the top swimmer, while attracting to the swimmer. The movement of these particles is also interesting in that the amplitude of their movement is growing in time in the y-direction. In comparison, we look at fluid particles with a wall at x −5 and compare it to the free space case for an initial separation distance of d 3 μm in Figures 8C,D. Mid-way along the swimmer in Figure 8C, at these early time points, we observe that trajectories of fluid markers behave in the same way, being pushed down in the z − direction regardless of whether the wall is present close by. This fluid particle is initialized close to but below the top swimmer and it is being pushed downward in the direction of the bottom swimmer. In contrast, we observe the effects of the wall in Figure 8D when initializing a fluid particle close to the wall and directly below the swimmer. In the case of the wall, the particle trajectory is following and getting close to the swimmer, increasing in the z-direction with hardly any progression in the x − direction due to the wall. The free space case shows the extensive movement in the x-direction, but little movement in the z-direction.

Quasi-Planar Swimmers
Everything presented thus far has been for swimmers with a planar preferred beat form. Due to interactions with a swimmer or the wall, nonplanar beat forms have emerged (e.g., Figure 7D). Since different species of sperm exhibit a variety of nonplanar beat forms [2,9], we now consider here the case of two quasiplanar swimmers where α 3 μm and β 1 μm in the preferred beat form in Eq. 3. The bottom swimmer is initialized with its centerline lying on the plane z 0 and the top swimmer's centerline is on the plane z d.

Free Space
We now wish to further characterize interactions of two swimmers with quasi-planar beat forms and understand how they are similar or different to swimmers with planar beat forms. In Figure 9, for an initial distance d 6 μm apart, the two swimmer's trajectories rotate around each other creating a bundle formed by the two flagella. At the same time, the two swimmers are attracting to each other, as shown in Figure 9B, where all distance metrics considered are oscillating and decreasing in time.
Here, there are no signs of repulsion between the swimmers as they reach a minimum distance between the swimmers on the order of 3-4 μm for t ∈ [4,5]s, similar to the swimmers with planar beat forms that were initialized as co-planar in Figure 2B and in contrast to those initialized in parallel planes in Figures  6A,E where repulsion was observed when starting d 3, 5 μm apart. For the quasi-planar case, as expected, the trace of the first point shown in the zoomed portions of Figure 9A exhibit a more complicated trajectory, known as the flagelloid curve (or f-curve), as previously recorded in experiments and simulations for a single sperm [32,39,58,59]. The flagelloid curve is shown in the yz-plane over the time interval from 2 to 3.2 s, where the curvature is higher at the bottom of the bundle and lower at the top of the bundle; this trend in curvature is consistent for the full time of the simulation. We have also considered the case of two quasi-planar swimmers initialized at a distance of d 30 μm in Figure 9C. In this case, the swimmer's trajectories show clear attraction between the swimmers. That is, the minimum distance between the swimmers in Figure 9D is monotonically decreasing. Here, the average minimum distance and the average distance between the first points coincide for the full simulation. The flagelloid curves for d 30 μm are also reported in the zoomed in portions of Figure 9C and exhibit a similar pattern to those in the zoomed in portions of Figure 9A.
The results reported in Figure 9 suggest that the fundamental dynamics in free space of two quasi-planar swimmers, in terms of attraction and repulsion, is similar to the dynamics of two coplanar swimmers reported in Section 3.1.1 and Figure 2. We also quantified the swimming speeds of a solo quasi-planar swimmer as well as a pair of quasi-planar swimmers ( Table 2). Again, similar to the dynamics of attraction, the swimming speed trends were similar to that of the co-planar swimmers. Relative to the swimming speed of a solo swimmer, a pair of swimmers 5 μm apart had a decrease in swimming speed whereas swimmers initially 30 μm apart had a very small increase in swimming speed (at earlier time points). For the preferred configurations studied, the quasi-planar swimmers were slower than the planar swimmers (by a few μm/s). We also emphasize that no difference in the results are obtained if the second swimmer was initialized with a centerline lying on the plane y d, instead of z d, i.e., translating on the y-axis instead of the z-axis. Figure 10A shows the dynamics near a planar wall for a pair of swimmers initialized a distance d 6 μm apart; the two swimmers attract to the wall and start rotating around each other. Similar to the free space case in Figure 9A, the swimmers continue to circle each other. However, with the wall in Figure 10A, they do not progress forward but stay a constant distance away from the wall, remaining perpendicular to the wall. The zoomed portion of Figure 10A shows the flagelloid curves traced by the first point on the swimmers. The curvature is approximately the same whether the swimmer is at the top or at the bottom of the bundle. This trend in curvature is consistent for the full time of the simulation and in contrast to quasi-planar swimmers in free space ( Figures  9A,C). In terms of the dynamics of attraction, after an initial transient period of approximately 1 s where the first points of the swimmer attract and then repulse, the swimmers reach an almost-constant average distance between the heads at ∼6 μm apart ( Figure 10C). Similarly, swimmers initialized 3 μm apart reach a constant distance apart around 1 s, but the heads repulse initially and level off at a distance ∼5 μm apart ( Figure 10B).

Near a Planar Wall
The case of two quasi-planar swimmers initialized at a distance of 30 μm apart and also near the wall at x −5 is shown in Figure 10D. In this case, the swimmer's trajectories show clear attraction, i.e., monotonic decrease of the average minimum distance between the swimmers ( Figure 10E). The results reported in Figure 10 suggest that the fundamental dynamics near a wall of two quasi-planar swimmers, in terms of attraction and repulsion, is similar to the dynamics of two co-planar swimmers near a wall reported in Section 3.1.1 and Figure 3.
In particular, we point out the strong similarity between Figures 3D-F and Figures 10B,C,E. We also emphasize that no difference in the results were obtained if the second swimmer was initialized with a centerline lying on the plane y d, instead of z d.

CONCLUSION
The ability of mammalian sperm to reach and fertilize the egg is aided by a multitude of dynamic interactions between swimmers, signaling molecules in the fluid, and walls of the female reproductive tract. In this work, we provide a detailed look at pairs of swimmers to characterize conditions that lead to emergent phenomena such as attraction or repulsion of swimmers. In free space, we observe long-term attraction of two swimmers in the case of initially co-planar sperm with preferred planar beat forms and sperm initially with centerlines in parallel planes with preferred quasi-planar beat forms. In contrast, sperm initially with centerlines in parallel planes and preferred planar beat forms exhibit oscillatory dynamics, alternating between periods of attraction and repulsion. For both of these, we emphasize that these classifications are for separation distances on a length scale smaller than the length of the flagella and greater than or equal to the beat amplitude. When sperm are swimming in close proximity to a wall, we observe attraction to the wall for planar and quasi-planar beat forms, i.e., even if swimmers are in close proximity when near the wall, they are still trapped at the wall. For sperm initialized in parallel planes with a planar preferred beat form and near a wall, due to the instability in the angle of the attracted swimmers, we observe significant rolling episodes that allow the swimmers to attain a beat form and distance apart that can then be maintained. The observation of this rolling behavior is important as it is proposed to be an important mechanism in sperm selection [56].
The results presented in this work further clarify and contextualize divergent results in the literature. For example, in the case of parallel sperm in free space, our farfield attraction results differ from previous results of [26] where only repulsion was observed and [31] where only attraction is observed. We are able to show that the swimmers in this configuration will not reach a stable configuration of attraction and will continue to oscillate between attraction and repulsion. Zooming in on a particular time frame and/or different parameter choice leads to these divergent behaviors. Understanding the complex interactions of beat form and elasticity of the flagella can also be utilized to design artificial microswimmers that navigate in complex environments [35][36][37].
The modeling framework used did not account for background fluid flows and limited the study to the case of two swimmers initialized in the same plane or with centerlines in parallel planes. In all of the cases where attraction is observed, we emphasize that perturbations to the flow would likely cause additional rolling and pitching. Similarly, we focused on the case of a purely homogeneous fluid with a constant viscosity. It is well known that the viscosity of the fluid in the female reproductive tract varies and will often exhibit nonlinear properties with respect to stress and strain [2,10]. We observed changes in the dynamics of attraction with increases in the viscosity and hence, we expect that nonhomogeneous or nonlinear fluid contributions will also act to change the frequency of rolling and the degree of pitching in the swimmers. This will be the focus of future studies.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
LC and SO conceptualized this research and completed the writing. LC and DD completed simulations and created figures. All authors have contributed to this article and have given approval for this submission.

FUNDING
The work of SDO and DD was supported, in part, by NSF DMS 1455270. SDO was also supported, in part, by the Fulbright Research Scholar Program. Simulations were run on a cluster at WPI acquired through NSF grant 1337943 and on a cluster provided by Research Computing at the Rochester Institute of Technology (60).