Sperm chemotaxis in marine species is optimal at physiological flow rates according theory of filament surfing

Sperm of marine invertebrates have to find eggs cells in the ocean. Turbulent flows mix sperm and egg cells up to the millimeter scale; below this, active swimming and chemotaxis become important. Previous work addressed either turbulent mixing or chemotaxis in still water. Here, we present a general theory of sperm chemotaxis inside the smallest eddies of turbulent flow, where signaling molecules released by egg cells are spread into thin concentration filaments. Sperm cells ‘surf’ along these filaments towards the egg. External flows make filaments longer, but also thinner. These opposing effects set an optimal flow strength. The optimum predicted by our theory matches flow measurements in shallow coastal waters. Our theory quantitatively agrees with two previous fertilization experiments in Taylor-Couette chambers and provides a mechanistic understanding of these early experiments. ‘Surfing along concentration filaments’ could be a paradigm for navigation in complex environments in the presence of turbulent flow.

Shear α [s −1 ] with / without chemotaxis / − Experiment from Ref. [20] / − Simulation with chemotaxis, open green triangles: without chemotaxis; mean ± SD), using fertilizability p f = 60% in Eq. (5) as single fit parameter. Experiment and simulation again agree reasonably except for the data point without flow α = 0 s −1 , which corresponds to a different experimental protocol. While the simulations with co-rotation overestimate the reduction of P fert at high shear rate α > 6 s −1 , these high shear rates are less relevant for the spawning habitat of H. rufescens.
A Shear flow around freely-rotating egg and minimal case of ballistic swimmer For all simulations (except Fig 1A), we use a simple shear flow αy e x as idealized paradigm for small-scale turbulence. At the relevant shear rates α and typical egg radii r egg ∼ 100 µm, the Reynolds number . From the fit, we obtain p f P sperm:egg ≈ 9% for the product of fertilizability p f and encounter probability P sperm:egg . Assuming a ballistic swimmer that is captured at the egg surface (Eq. (5) with P sperm:egg (t max ) = 1 − exp(−qt max ) and rate q = πr 2 egg v h ρ egg = 0.02 s −1 ), we find p f ≈ 10% for exposure time t max = 120 s and ρ egg = 1.5 · 10 4 ml −1 . This value p f is used in Fig 4 and  Re = αr 2 egg /ν ≤ 0.1 is sufficiently small to justify the use of the analytical Stokes equation for viscous flow v ext (r). Throughout, we consider the co-moving frame of the egg allowing us to assume that the egg is at the origin r = 0. We introduce dimensionless coordinatesr = r/r egg and the dimensionless flow fieldv ext (r) = 2vext(r) αregg . The components of this flow field read ( [4], Eq. (12)) v ext,x = 2ŷ −ŷ 1 +Ω r −3 +r −5 − 5x 2ŷ r −5 −r −7 v ext,y =x 1 +Ω r −3 −r −5 − 5xŷ 2 r −5 −r −7 v ext,z = − 5xŷẑ r −5 −r −7 (S1) where no-slip boundary conditions on the surface |r| = 1 of the freely-rotating spherical egg are assumed. The egg rotates according to the undisturbed flow vorticity with the dimensionless rotation rateΩ = −1, corresponding to an rotation of the egg with angular velocity Ω = − α 2 e z . It is instructive to consider a ballistic swimmer in the above flow field v ext as a reference for the analysis of more complicated cases, such as swimmers performing chemotaxis. For instance, without flow or chemotaxis, sperm cells are considered to swim along a straight helix with helix radius r 0 much smaller than the egg radius. These sperm trajectories are well approximated by a ballistic swimmer moving along the helix axis h with net swimming speed v h . If the ballistic swimmers and the target eggs (with density ρ egg ) are uniformly distributed, the steady-state rate q at which a swimmer hits an egg is given by q = π(r egg + r 0 ) 2 v h ρ egg ≈ πr 2 egg v h ρ egg . If ballistic swimmers become trapped at the egg on encounter, this corresponds to the encounter probability P sperm:egg (t) = 1 − exp(−qt) (and fertilization probability P fert according to fertilization kinetics, see Eq. (5)). If ballistic swimmers are additionally Previous measurements of fertilization probability P fert ( ) for sea urchin S. purpuratus at strong turbulence, characterized by density-normalized dissipation rate (filled gray triangles) [1,25] and our corresponding simulations P fert (α) as function of shear rate α (open green triangles, mean ± SD) match well, using a single fit parameter a = 0.023 that relates dissipation rate and typical shear rate α (with the known relationship α( ) = a /ν [2,3]). Analogous to Fig 4, the case of low shear rates is well described by the limit case of a ballistic swimmer in the absence of flow α = 0 s −1 (dotted horizontal line, Eq. (5) with P sperm:egg (t) = 1 − exp(−qt) and rate q = πr 2 egg v h ρ egg ). The fertilizability p f = 10% is obtained from an independent experiment [25], see Fig B. From the experimental protocol, we estimate a high background concentration c bg = 500 − 4000 nM of chemoattractant, which renders sperm chemotaxis ineffective.
convected by an external fluid flow field v ext , we can characterize q (and thus P sperm:egg and P fert ) in terms of an universal curve: We introduce the dimensionless parameter f = αregg 2v h , which compares shear rate to net swimming speed. The combined velocity field of active swimming and fluid flow is now Note that without co-rotation, h does not change. Thus, for any h, all possible velocity fieldsû are given by a single one-parameter family parametrized by f . For each of these fields, the dimensionless rate of swimmersq (f, h) reaching the egg from |r| 1 specifies the actual rate q for any set of parameters α, r egg , v h , ρ egg with the same parameter f by We obtain a universal curve for q by computingq(f, h) numerically for all f and h and averagê q(f ) = q(f, h) h over all directions h, see Fig 4 for corresponding P fert . A prominent feature of the universal rate is that it vanishes at large shear ratesq (f → ∞) → 0. In the absence of flow α = 0, we haveq (f = 0) = π.
We compute the universal rateq efficiently by integrating a uniform grid of initial conditions on the surface of the egg, with |r| = 1 att = 0, backwards in time according to the velocity fieldû. Each initial condition is integrated until it either returns to the egg |r| (t) = 1 (fail) or leaves the outer boundarieŝ r(t) =r max (success) withr max 1. As the flow is volume conserving, the results are independent of the choice of the outer boundaryr max , as long asr max is sufficiently large to ensure the absence of closed orbits beyond it. We chooser max = 4 as numerics show that the in-and outflow on this sphere differs only by 4% between the Stokes flow around the freely-rotating sphere and the undisturbed simple shear flow, for which it is known no closed orbits exist. Based on the intersections with the outer boundary, the flow reaching the egg is interpolated. This is done for a grid of swim directions h. For efficiency, we exploit the symmetries of the Stokes flowv ext (x,ŷ,ẑ) ; thus, it is sufficient to consider h z ≥ 0 and h y ≥ 0, respectively.

B Equations of motion for navigating sperm cells
We simulate the swimming path r(t) of a sperm cell in a concentration field c(r) of chemoattractant in the presence of an external fluid flow field v ext (r). For this, we extend a previous theory of chemotaxis of marine sperm cells along helical paths [5][6][7][8] by incorporating convection and co-rotation by flow: The sperm cell is described in terms of the time-dependent center position r(t), averaged over one flagellar beat cycle, and the set of ortho-normal vectors e 1 (t), e 2 (t), e 3 (t) of the co-moving coordinate frame, where the vector e 1 (t) points in the direction of active swimming with speed v 0 . The equations of motion reaḋ The two angular velocities, Ω h and Ω f , describe the rotation of the coordinate frame due to helical chemotaxis and external flow, respectively. For Eq. (S4) a constant swim speed is assumed and motility noise is neglected; the persistence length of sperm swimming paths in the absence of chemoattractant cues was estimated as 3 − 25 mm [9] which validates this assumption. Note that Eq. (S4) is also valid for time-dependent concentration and flow fields. Note further that the quantitative comparison of the two angular velocities in Eq. (S4) already suggests that the rotation due to external flow is negligible, as the rate of change due to external flow Ω f ∼ α = 0 − 1 s −1 (see Eq. (S8)) is always smaller than due to the helical motion Ω h ∼ τ 0 v 0 = 3 − 13 s −1 (see Eq. (S5) and parameters in Table A). Without external flow or chemotaxis, cells swim along a helical path with constant path curvature κ(t) = κ 0 and torsion τ (t) = τ 0 . The angular velocity Ω h is defined by the Frenet-Serret equations where the coordinate frame e 1 , e 2 , e 3 corresponds to the Frenet-Serret frame of r(t), i.e., tangent, normal and bi-normal vector. During chemotactic steering, sperm cells dynamically regulate curvature κ(t) and torsion τ (t) of active swimming according to the output a(t) of a chemotactic signaling system Here, the sensori-motor gain factor ρ characterizes the amplitude of chemotactic steering responses. The chemotactic signaling system takes as input the local concentration c(r(t)) at the position of the cell This minimal signaling system comprises sensory adaption with sensitivity threshold c b and relaxation with time scale µ to a rest state a = 1 for any constant stimulus c(r(t)) = c 0 . The variable p describes an dynamic sensitivity which is regulated down when the stimulus is high, or regulated up when the stimulus is low (a loose analogy would be that p corresponds to the opening of our eye's pupils as adaption to brightness). In principle, p and a could have different time-scales [5]. However, equal time-scales automatically ensure that the phase-lag between small-amplitude oscillations of the input signal c(r(t)) and resulting oscillation of the output signal a(t) attains the value π/2 optimal for helical chemotaxis [6]. This special case is sufficient for the purpose of a minimal model. The gain factor ρ sets the rate of chemotactic steering. While ρ could depend on the chemotactic signal by a feedback mechanism [10], we assume here a constant gain factor ρ = 5 for simplicity. The values of all parameters are listed and discussed in Sec G. We approximate the angular velocity Ω f for co-rotation by external flow using the Jeffery equation for a small prolate spheroid with major axis along e 1 [11,12] , with the flow vorticity ω, the strain rate tensor E, and a geometric factor G = g 2 −1 g 2 +1 , which depends on the aspect ratio g ≥ 1 of major to minor axis of the spheroid. Together with Eq. (S4), Eq. (S8) describes the cell rotation due the flow, i.e., The first term in the first line of Eq. (S8) describes rotation of a spherical body due to flow vorticity and the second term the correction for non-sperical bodies that can be approximated as spheroids. For a swimming sperm cell, we take the swim direction e 1 as effective major axis, and employ an effective aspect ratio, g = 5, reflecting the ratio of the length of the flagellum and a typical beat amplitude [13]. Note that in general instead of e 1 , the major axis could be any co-moving vector.
We numerically integrate the equations of motion, i.e., Eqs. (S4,S7), using an Euler scheme with fixed small time step dt. For efficient computation, Rodrigues rotation formula [14] with respect to the co-moving coordinate frame is used to integrate e 1 , e 2 , e 3 , resulting in faster computation compared to the algorithm used in [10].

C Analysis of concentration filaments
Turbulent flows cause turbulent mixing of diffusing chemicals and generate filamentous concentration fields. As a minimal model, we simplify the turbulent flow and the filamentous concentration field by the case of a simple shear flow. We consider a spherical egg located at the origin r = 0 releasing chemoattractant with diffusion coefficient D at a constant rateQ in the presence of shear flow v ext (r) given by Eq. (S1). We compute the time-dependent concentration field c(r, t) of chemoattractant numerically using Lagrangian particle tracking, see Sec G. We empirically find that the far-field at distances r r egg is well approximated by a generic profile, see Fig 1B for which describes a concentration filament with time-dependent parameters c 0 (t), k(t), a y , as well as timeand position-dependent variance and midline profiles σ(x, t) and y 0 (x, t), respectively. This formula for the concentration filament is consistent with results obtained using the analytic solution for an instantaneous point source in a shear flow αy e x , see below. We present and discuss scaling laws for the parameters in the following. While these dependencies are not explicitly required for our theory, they demonstrate the universality of our theory. Finally, we use these scaling laws to quantify how the filaments become longer and thinner with increasing shear rate α.
From numerical simulations, we empirically find the following scaling laws of the parameters from where all parameters except p 2,t0 are positive. Note that in a turbulent flow the time t in which the filament is formed may scale with the Kolmogorov time t ∼ τ Kol ; in this case σ would scale as the Batchelor length σ ∼ √ Dτ Kol . We also found power-law dependencies for the coefficients p 0,t0 (t), p 1,t0 (t), and p 2,t0 (t). The factor a y appears to be constant for sufficiently large t. These numerical observations become plausible by analysis of a point source in shear flow. The Fokker-Planck equation for this case can be written in dimensionless form by using the Batchelor scale √ Dτ Kol ∼ D/α to re-scale to dimensionless coordinateš with shear rate α, and release rateQ of the source (i.e. Consequently, the solutionč ř,ť of this equation can be re-scaled to the solution c(r, t) for any set of parameters α, D,Q. For the above form of the far-field of the filament, this implies that the parameters δ k , δ c0 , p 1,y0 and p 2,t0 are universal as they are invariant under the re-scaling Eq. (S12). The analytical solution for the dimensionless concentrationč reads ( [15], Eq. (18)) with Greens functionǦ, i.e., the solution for an instantaneous source at the origin ( [16], Eq. (26)) (S14) While the integral Eq. (S13) cannot be solved analytically, it explains the empirical scaling for the parameters in Eq. (S9) heuristically: It is reasonable to assume that for anyx, the parametery 0 (ť) is close to the pointy max of the maximal concentration ofǦ(ř,ť). From ∂yǦ(x,y, in accordance with the fitted power-law. The power law c 0 (t) ∼ t − 3 2 , as suggested by numerics, is plausible sinceǦ (0,š 1) ∼š − 5 2 , which impliesč 0,ť 1 ∼ ť 0 dšš − 5 2 ∼ť − 3 2 . We introduce the concentrationč max at the centerline of the filamentč max x,ť =č x,y 0 x,ť ,ž = 0 . We make the ansatzč max x,ť =č 0 (ť) exp −ǩ ť |x| and derive a power-law forǩ(ť) in the following.
We expect thatč max scales proportional to the summed contributions of the Greens functions at the timedependent centerline, hence we estimate (assumingť 1, we approximateť 2 + 3 →ť 2 , 1 +ť 2 /12 . (S16) We are interested in the shape of the concentration filament up to a maximal distancex max at which the concentration at the centerline decayed to a fraction ι ofč 0 ,č max (x max ,ť) = ιc 0 (ť). Any asymptotic tails beyond this distance will likely not be relevant for chemotaxis. Since the decay ofč max as function ofx is dominated by the numerator in Eq. (S16), the distancex max has a time-dependencyx max ť ∼ť 3 2 according to the argument of the complementary error-function erfc. Using Eq. (S16), we estimate the time-dependency ofǩ(ť) from where we introduced a(x,ť) = 3/4ť − 3 2x . The crucial point is that for 0 ≤x ≤x max (ť), the variable a varies only in a finite interval 0 ≤ a ≤ a max with upper bound a max = a(x max (ť),ť) ∼ť − 3 2x max (ť) independent of timeť. This allows us to approximate ln erfc (a) by its Taylor expansion for small a 1 in the last step of Eq. (S17). We concludeǩ ť ∼ a max /x max (ť) ∼ť − 3 2 , as suggested by numerics. From the above considerations follows that for a constant exposure time t max the filaments become longer and thinner with increasing shear rate α: From Eq. (S12) follows for the dimensionless exposure timě t max = t max α and thus the exponentǩ of the dimensionless version of Eq. (S9) scales withǩ ∼ (t max α) − 3 2 , which implies according toǩx = kx a scaling of the effective decay length This means that for constant exposure time t max the effective length of the filament increases with shear rate α. Analogously, fromσ ∼ √ t max α andč 0 ∼ (t max α) − 3 2 follows with the re-scaling Eq. (S12) for the effective decay length σ away from the center of the filament and the base concentration c 0 (S19) The combination of the effective decay length σ being independent of α and the base concentration c 0 decreasing with increasing α means that the effective width of the filament decreases with increasing α.

D Chemotactic navigation within filament
We derive an effective equation of motion for chemotactic navigation within a typical concentration filament. For simplicity, we initially ignore interaction with the flow and assume that the motion is effectively two-dimensional, i.e., in the xy-plane. Additionally, we employ a two-dimensional version of Eq. (S9) for the concentration filament, setting a y = 1, We introduce the centerline r h (t) = (x(t), y(t), 0) of the helical swimming path r(t), withṙ h = v h h. From a previously established equation for r h [5,6], we havė describing the alignment of the helix axis h with the local gradient ∇c(r h (t)) of a concentration field c(r). The first equation corresponds to ballistic motion along the helix axis h(ϕ) = cos ϕ e x + sin ϕ e y with net swimming speed v h = v 0 τ 0 / κ 2 0 + τ 2 0 . The second equation describes chemotactic turning of the orientation angle ϕ, where Ψ denotes the angle enclosed by h and the local gradient ∇c(r h (t)). Here, c b denotes the adaption threshold and v ϕ the chemotactic turning speed, v ϕ = ρv h κ 2 0 / κ 2 0 + τ 2 0 , v ϕ > 0, with the gain factor ρ and helix parameters κ 0 , τ 0 . We apply this general theory, Eq. (S21), to the filamentous profile Eq. (S20) and obtain a single dimensionless ODË withẊ 2 +Ẏ 2 = 1,Ẋ = 0 and a single dimensionless parameter Here, we introduce a characteristic time-scale τ , as well as re-scaled coordinates Dots denote differentiation with respect to re-scaled time T = t/τ , e.g.,Ẏ = dY /dT . The time scale τ is the geometric mean of a characteristic time-scale σ/v ϕ of chemotactic steering and a typical time σ/v h for traversing the cross-sectional width σ of the filament if steering was absent. We have an equation for X analogous to Eq. (S22) (which requiresẎ = 0 and covers the caseẊ = 0), The factor c/(c + c b ) in the effective equations of motion, Eqs. (S22, S25), represents a 'dimmer switch' that attenuates chemotactic navigation at low concentration c. Thus, it is reasonable to define the filament as the region where c(r) ≥ c b . In the following, we focus on the dynamics within the filament and approximate c/(c + c b ) ≈ 1.
The effective equation of motion, Eq. (S22), describes a damped, non-linear oscillator: The first terṁ X 2 Y originates from the perpendicular component ∇ ⊥ c = (e y · ∇c) e y of the concentration filament and governs the observed oscillations of sperm cells around the centerline Y = 0 of the filament. Heuristically, these oscillations result from sperm cells slowly aligning their helix axis h parallel to ∇ ⊥ c while approaching Y = 0. At Y = 0, ∇ ⊥ c changes its direction, yet sperm cells overshoot due to their finite chemotactic turning speed v ϕ < ∞, before they eventually make a 'U-turn'. The second term sgn (X) γẊẎ in Eq. (S22) originates from the exponential decay of concentration along the centerline of the filament and changes the amplitude of the oscillation. In particular, for sgn XẊ < 0, i.e., sperm cells surfing towards the egg, the oscillation is damped, whereas for sgn XẊ > 0, i.e., sperm cells surfing away from the egg, it is amplified. This increase in amplitude can cause sperm cells that are surfing away from the egg to eventually turn around, redirecting them towards the egg. A linear stability analysis of Eq. (S22) around the case of a non-oscillating trajectory Y,Ẏ = (0, 0) yields the eigenvalues ω 1,2 of the Jacobian of the linearization, which define a harmonic oscillator with dimensionless damping ratio ζ and dimensionless oscillation frequency 1 − ζ 2 . This analytic result agrees with full simulations of helical chemotaxis in threedimensional space, see Fig D. Note that the predicted exponential decay of oscillation amplitude, exp (ζT ) = exp (γ/2 · t/τ ), is independent of x since γ/τ is independent of σ 2 (x). Interestingly, both for Eq. (S22) and full simulations, the angle at which trajectories intersect the centerline Y = 0 of the concentration filament is essentially independent of the angle at which they first entered the filament at Y (c = c b ), provided Y (c = c b ) is sufficiently large: For smaller Y (c = c b ), i.e., outer and thus thinner parts of the filament, trajectories will simply pass through the filament, unable to execute a successful turn before they have left the filament again. As the width of the filament decreases away from the egg, this implies that filament surfing will be operative, at most, up to a maximal distance from the egg (which depends on the entry angle), characterized by p in . If we account for convection by shear flow v ext = αy e x , Eq. (S25) changes toẊ →Ẋ + ατ (Y + y 0 (X)/L). Note that due to sgn (y 0 (x)) = sgn (x), sperm cells that surf within the filament towards the egg swim on average against the external flow. Surfing along filaments can be described as damped oscillation. Distance d = (y − y 0 ) 2 + z 2 from the centerline of concentration filament Eq. (S9) superimposed for n = 9 sperm trajectories simulated according to Sec B (black). Trajectories are shown after they entered the surface of the concentration filament, defined by c(r(t)) = c b , and shifted in time to align the first oscillation peak at t 0 = 0. Remarkably, all trajectories display stereotypic oscillations that overlap perfectly, despite the fact that trajectories entered the filament at different x-positions and initial direction angles. The observed damped oscillation are well reproduced by a minimal analytical theory for the centerline of the helical swimming path, which predicts damping ratio and oscillation period (dashed red line, see Eq. (S26)). Parameters as in Fig 1B, corresponding to A. punctuala.

E Minimal theory for sperm-egg-encounter probability
We provide an estimate for the encounter probability P sperm:egg , building on the effective equation of motion of the helix axis derived in Sec D. The fertilization probability P fert is obtained then from P sperm:egg using fertilization kinetics, Eq. (5). For P sperm:egg , we decompose the search problem for the egg into an outer search problem of finding the concentration filament and an inner search problem of surfing along the filament. We obtain (exploiting the symmetry between the two branches of the filament for x < 0 and x > 0) Here, we introduce the following quantities: • the cross-sectional area A(x) at the position x of the filament, which is defined by c(r) ≥ c b , i.e., A(x) = ∞ −∞ dy dz Θ(c(x, y, z) − c b ), with the Heavyside-function Θ (Θ(c > 0) = 1 and Θ(c ≤ 0) = 0), • the circumference S(x) corresponding to the cross-section, • the average probability p in (x, t max ) that a trajectory entering the filament at x > 0 will surf along it and reach the egg within exposure time t max , • the mean steady-state flux j out of trajectories arriving at the surface of the filament, and • the time limit t out for the outer search problem.
These quantities are explained in detail below. The first term in Eq. (S27) accounts for sperm cells found inside the concentration filament already at t = 0, assuming a random uniform distribution of initial positions. The second term in Eq. (S27) accounts for trajectories, which first search for the filament and, after encountering the filament, surf along it towards the egg. We compute the probability p in (x, t max ) of successful inner search numerically using the effective equation of motion for the helix axis Eq. (S22) as function of entry position x and exposure time t max . Specifically, we average over simulations of Eq. (S22) with uniformly distributed initial entry points and isotropic initial directions, i.e., entry angles. In order to account for the ellipsoidal cross-section of the concentration filament with σ y = σa y , σ z = σ, we average results for σ y and σ z . From the successful trajectories, we also obtain the mean travel time t in within the filament, which represents a conditional mean first passage time. Accordingly, we set the maximal time t out allowed for the outer search t out (x, t max ) = t max − t in (x, t max ) if p in > 0 and t out = 0 else.
Note that the first term in Eq. (S27) can be written as V eff ρ egg with an effective volume V eff = 2 ∞ 0 dx A(x)p in (x, t max ) of the concentration filament, weighted by the probability p in of successful chemotaxis to the egg. This contribution is negligible compared to the second term for long exposure times t max and low egg densities ρ egg .
The flux j out of trajectories arriving at the surface of the concentration filament can be determined by a fit to P sperm:egg (α) from simulations at different shear rates α. Alternatively, we can estimate j out by treating sperm cells outside of the filament as ballistic swimmers with net swimming speed v h and uniformly distributed random positions r and orientations h with probability distribution Assuming that the filament is convex, each point on its surface is reached at time t from initial conditions on a surface of a half-sphere with radius v h t. The flux of trajectories with direction h into the filament at r 0 is j out (r 0 , h) = −n · v h hp sperm (r 0 , h) for n · h < 0 and j out (r 0 , h) = 0 else, where n denotes the outer surface normal vector at r 0 . For the constant density p sperm (r 0 , h) = p sperm the total flux of sperm cells into the filament is j out = 2π 0 dϕ π/2 0 dθ sin θj out (r 0 , h(ϕ, θ)) = p sperm πv h , where we use spherical coordinates ϕ, θ with e z = n to express h. Note that an isotropic distribution of orientations h is a simplification, since co-rotation by flow alters this distribution, see Sec F.
Despite the simplifications made, Eq. (S27) can quantitatively account for the encounter probability in full simulations, see Fig 2. In particular, we find that the numerical fit for j out = 0.063 m −2 s −1 is close to our simple estimate for a ballistic swimmer j out = ρ egg v h /4 = 0.04 m −2 s −1 . Of course, our simple theory has limitations: First, trajectories are three-dimensional, not two-dimensional, and are characterized by oscillations both in y-and z-direction. As a result, sperm trajectories are super-helical, which reduces the effective speed along the filament. Second, our theory does not account for the fact that some sperm cells may miss the egg on the first attempt, and find it only after reversing their motion in x-direction, which increases the mean time t in to find the egg. Preliminary simulations suggest that the difference between simulations and theory in Fig 2 indeed originate from this effect. Finally, co-rotation is neglected in the simple theory. However, this is justified for α τ −1 , see Eq. (S24), i.e., when rotation due to navigation is much faster than co-rotation due to flow. Note that simulations with neither convection nor co-rotation exhibit also an optimal shear rate α * , but at higher shear rate and different encounter probability. The reason is that convection implies a flow opposing surfing towards the egg, which increases t in compared to the case without convection. Thus, P sperm:egg increases for large α when convection is not included, resulting in a shift of α * .
For the experiment of Zimmer and Riffell (data reproduced in Fig 3 and Fig A), we estimate a high background concentration of chemotattractant c bg ∼ 4 nM, see Sec G. Adding a background concentration c → c + c bg in Eq. (S20) leads to an effective, higher threshold c b,eff = c b + c bg in Eq. (S22). Consequently, the volume of the filament with sufficiently high concentration c(r) ≥ c b,eff is situated only in the vicinity of the egg. While our far-field theory of filament surfing does not apply directly to this special near-field case, we can make a simple estimate: We assume that sperm cells always swim directly towards the egg within the concentration plume defined by c(r) ≥ c b,eff due to the close-to-spherical shape of the plume. Thus, sperm cells entering the plume at x 0 = 0 approach it with net radial speed v h , as the external flow only convects the sperm cells parallel to the egg surface, see Eq. (S1). A second, alternative calculation applies if sperm cells enter the plume at x 0 r egg : In this case, we can estimate the net speed towards the egg byẋ = αy 0 (x) − v h . This yields for the distance x(t) from the egg, αb exp (αbt) (using y 0 (x) ≈ bx, see Sec C). We use these two limit cases to compute p in and t out for Eq. (S27) and obtain similar fertilization probabilities P fert (α) in both cases. For these limit cases, P fert (α) displays a similar decay as function of α as the simulation results without co-rotation, see Fig 3. In particular, the fitted flux j out = 4.8 · 10 3 m −2 s −1 is consistent with the theoretical value j out = ρ egg v h /4 = 7.5 · 10 3 m −2 s −1 .

F Analytic solution of Jeffery equation in shear flow
As shear flow is a fundamental paradigm for small-scale turbulence, we present here the analytic solution to the Jeffery equation, Eq. (S8), for particles suspended in simple shear flow. The application to helical swimmers is discussed. The results provide the distribution of helix orientations h on the periodic boundary used in the simulations, i.e., p sperm in Eq. (6). In particular, the results quantify the common notion that non-spherical swimmers align their major axis parallel to the flow direction. In fact, these swimmers rotate all the time, but with non-constant rotation rate, causing these swimmers to spend more time aligned with the flow axis. Consequently, the time-average of the orientation vector is not zero, but aligned with the flow axis. Note that analytic results for Poiseuille flow can be found in [17,18].
While ρ e is periodic in time with period T by Eq. (S31), we can compute a time-average over one period, starting with a uniform distribution of directions e at t = 0. The time-averaged density displays a maximum at the axis of flow e = ±e x and a minimum at the shear axis e = ±e y . These extrema vanish for a sphere (G = 0) and become more pronounced with increasing G. While the above results are derived for the case of a suspended particle, numerical simulations show that they also approximately apply to the centerline r h (t) of a helical swimmer with helix axis h (without chemotaxis) if we use an effective aspect ratio g eff . Specifically, the dynamics of the helix axis h resembles the above solutions with a smaller aspect ratio 1 ≤ g eff ≤ g. This approximation is valid for small times t and at small α, i.e., as long as the helix period is much smaller than the period T (g). For instance, we fit g eff = 1.3 ± 0.1 (G eff = 0.26) for the sea urchin helix parameters and g = 5 (G = 0.92). This effective parameter is a result of averaging the instantaneous co-rotation for the swimming direction e 1 with parameter G over one period of helical swimming. Generally, g eff depends on the angle between h and e 1 . For larger α, complicated behavior of h is observed with limit cycles and stable fixed points, which is consistent with recent results for Jeffery equation in perturbed shear flow [19]. We use the value g eff in all simulations to determine the periodic boundary conditions at the boundary of the simulation domain.

G Choice of parameters
Parameters used throughout the three simulation scenarios (Arabacia punctuala for Fig 1B, Fig 2, Fig D, Strongylocentrotus purpuratus from [1,21,22,25] for Fig 4, Fig C and Fig B, Haliotis rufescens from [13,20] for Fig 3 and Fig A) are listed in Table A and discussed in the following.
Mean path curvature κ 0 = 0.065 µm −1 and mean path torsion τ 0 = 0.067 µm −1 of the helical paths are set according to three-dimensional tracking of A. punctuala sperm cells [8]. Three-dimensional tracking for S. purpuratus give similar values [23], though with larger error intervals. Moreover, the sperm morphology for A. punctuala [8,10], S. purpuratus [24], and H. rufescens [26] is similar, which justifies the use of the same helix parameters for all three species. Likewise, the effective aspect ratio g = 5 between major and minor axis of a sperm cell, i.e., length of flagellum divided by typical beat amplitude, suggested for H. rufescens [13] is employed for all three species in the Jeffery equation Eq. (S8). We observe that simulation results are largely independent of the precise value of g. The signaling time-scale µ = 1/ v 0 κ 2 0 + τ 2 0 is chosen to ensure the optimal phase-lag between concentration input c(r(t)) and motor response a(t) [5,27], see Eq. (S7), consistent with experimental observations [8]. For all three species, the gain factor is set as ρ = 5, corresponding to the mean of the values used in [10]. This value reproduces typical bending rates of helical swimming paths as observed in experiments [8]. The threshold of sensory adaption c b = 10 pM is chosen as suggested in [28]. At the concentration c b , about 20 chemoattractant molecules would diffuse to a sperm cell during one helical turn. Note that sea urchin sperm cells respond to single chemoattractant molecules [29]; the change in intra-cellular calcium concentration caused by the binding of chemoattractant molecules as function of stimulus strength becomes sub-linear already for chemoattractant concentrations on the order of c b [28]. For A. punctuala, other parameters were also tested, i.e., ρ = 2 and c b = 1 pM, which yielded qualitatively similar simulation results and again agreement of theory and simulations. Note that the experimental protocol used in [20] for H. rufescens results in a substantial background concentration of chemoattractant, which we estimate as c bg ∼ 4 nM (experiments are conducted 10 − 30 min after spawning at a high density of eggs ρ egg = 10 3 ml −1 with the known release rateQ = 0.18 fmol min −1 of chemoattractant [30]). According to our theory, such a background concentration causes effectively a higher sensitivity threshold c b,eff = c b + c bg (see Sec E), which may the be reason for the higher behavioral threshold 300 pM observed in [20]. In the case of S. purpuratus, we estimate an even higher background concentration, c bg ∼ 500 − 4000 nM, which renders chemotaxis ineffective. For this estimate, we use that experiments were conducted 1 − 8 h after spawning at a high egg density ρ egg = 1.5 · 10 4 ml −1 [21,25] and assume a release rateQ = 0.46 fmol min −1 of chemoattractant as for A. punctuala [28].
For the swimming speed v 0 of sperm cells along helical paths for both sea urchin species, we use the measured value v 0 = 200 µm s −1 from [8]. Note that some experiments effectively measure the net swimming speed along the helix axis v h = v 0 τ 0 / κ 2 0 + τ 2 0 , which is smaller than v 0 . For H. rufescens, we use the speed v h measured during the same experiment [20]. Note that this experiment also indicated chemokinesis, i.e., higher swimming speeds at elevated chemoattractant concentration, an effect which we neglect here for simplicity. For A. punctuala, we use the diffusion coefficient D = 239 µm 2 s −1 and release rateQ = 0.46 fmol min −1 of chemoattractant [28]. For this simulation, we assume a low egg density ρ egg = 10 −3 ml −1 , which yields the radius r max = 6 · 10 4 µm of the outer boundary centered around the egg according to ρ egg = 4πr 3 max /3 −1 . For this reference case, the filament is completely included inside the simulation domain for all considered shear rates α. The exposure time t max = 360 s is chosen comparable to the experiment in [25], where t max = 120 s. While for this work the exposure time t max is set by the protocol of the considered experiment, in a generic turbulent flow t max corresponds to the time-scale of flow changes, i.e. scale with the Kolmogorov time t max ∼ τ Kol . For comparison with the experiments with S. purpuratus and H. rufescens, the radius r max is computed directly from the stated egg density ρ egg = 4πr 3 max /3 −1 . From the 5 vol% solution with r egg = 40 − 55 µm ( [21], pg. 161), we infer a range ρ egg = 0.9 − 3.4 · 10 4 ml −1 for the experiments with S. purpuratus. This estimate already takes into account that, according to the experimental protocol, the above egg solution is mixed 9 : 1 with sperm solution [21,25]. Likewise, from the range of sperm densities ρ sperm = 1.9 − 3.1 · 10 6 ml −1 in Fig. 4 of [25] and the estimate ρ sperm = 4 · 10 6 ml −1 from pg. 59 of [21], both before 9 : 1-dilution, we infer a final concentration ρ sperm = 3.9 · 10 5 ml −1 . We use the kinematic viscosity ν = 10 −6 m 2 s −1 of sea water at room temperature.

H Numerical Simulation
The equations of motion are integrated using an Euler scheme with fixed time step dt . For all time integrations, a time step dt = 10 −3 s is used. Integration with smaller dt = 10 −4 s for some test cases gave consistent results. The number N sperm of sperm cells simulated in each case is 10 5 , except for S. purpuratus, where N sperm = 10 4 is used. The concentration field is computed from Lagrangian particle tracking with Euler-Maruyama method for the Fokker-Planck equation with v ext from Eq. (S1). Test particles were released at random points of the surface of sphere of radius r egg located at the origin. In total, we used 4 · 10 6 test particles, which corresponds effectively to 1.6 · 10 7 particles by exploiting symmetries of the flow field. Concentrations are evaluated on a cubic 50 × 50 × 50 grid, spanning in each dimension from −r max to r max , and then interpolated by a spline interpolation of order 3. This grid is sufficiently fine to resolve the details of the concentration filaments.
The rapid convergence to a near-steady state allows to use a static concentration field corresponding to exposure time t max for each simulation. We checked for test cases that full simulations with time-varying concentration field do not yield different results. The implementation of an unsteady shear flow for a shear rate α used as illustration in Fig 1A is inspired by [31]: We use the flow field v ext (r, t) = α (r, t) r · e y (t) e x (t), where the shear axis e y (t) and the flow axis e x (t) are subject to a three-dimensional random walk on the unit sphere with rotational diffusion coefficient D rot = πα. The shear rate profile is given by α (r, t) = √ 2 α sin (2πt/T α ) h(r). The shear rate α (r, t) decays as h(r) with distance r away from the center. This decay h(r) mimics the decay of velocity from the center of a vortex. We use the decay of an Lamb-Oseen vortex h(r) =  Fig 4, Fig B,

I Parameter study
In order to demonstrate the sensitivity of the quantitative results, shown in Fig 2 of the main text, on the parameters, we computed the encounter probability P sperm:egg (α) for a range of exposure times t max , egg densities ρ egg expressed in terms of boundary radii r max , threshold of sensory adaption c b , and gain factors ρ as shown in Fig E, Fig F, Fig G, and Fig H, respectively. In all cases, there is a pronounced optimum present at some intermediate shear rate α * , where the position of the optimum α * is only slightly affected by the parameter variations. The parameters mostly affect the height of the optimum, in particular t max and r max , and its ratio to the flow-less α = 0 case, see Table B. This suggests that the existence of an optimum is quite insensitive to parameter variations, i.e., regardless of how the parameters are adapted, fertilization is optimal at an intermediate shear flow for a broad physiological range of parameter values.
For each parameter study all parameters but one are kept constant on the values reported for A. punctuala in Table A. The only exception are t max and r max whose base values are lowered to t max = 90 s and r max = 15 mm for numerical efficiency, i.e., the blue dots in Fig E, Fig F, Fig H, and Fig G correspond always to the same parameters. The increase of the exposure time t max in Fig E from 45 to 90 s causes only a slight decrease of the optimal shear rate α * from 0.3 to 0.1 s −1 , but increases the absolute encounter probability P sperm:egg (α * ) by an order of magnitude. Such an increase of P sperm:egg (α * ) is also observed for the increase of the egg density, i.e., the decrease of the boundary radius from 30 to 10 mm in Fig F. These increases are in accordance with the simple argument that longer search time or smaller search volume increases the chances of finding the egg. The advantage of the optimum to the flow-less case P sperm:egg (α * )/P sperm:egg (α = 0) varies for both parameters between a factor 2 and 14, see Table B. In contrast, the variation of c b and ρ hardly affects the optimum in terms of α * and P sperm:egg (α * ) but rather alters the probability in the absence of flow: Increasing c b or decreasing ρ increases P sperm:egg (α = 0 s −1 ). This is probably an effect of signal-noise, originating from the computed concentration field which, due to the very nature of Lagrangian particle tracking, can exhibit low signal-to-noise ratio at low concentrations, i.e. at the surface of the concentration plume. (Note that our model does not explicitely account for sensing noise [34].) This noise results in an effective reflection of incoming sperm trajectories at the surface of the plume for increasing sensitivity of the concentration measurement, expressed by c b , or increasing reaction to signal stimulus, expressed by ρ, see also discussion in [10,35]. The effect is expected to be much smaller for concentration filaments at α > 0 s −1 as the concentration gradient towards the center of the filament is higher and thus the signal-to-noise ratio generally higher as for a concentration plume solely established by diffusion in the flow-less case. Flow-dependent sperm-egg encounter probability for different exposure times t max . Encounter probabilities P sperm:egg (α) as function of external shear rate α for three values of sperm-egg exposure time t max obtained from simulations with co-rotation. (symbols according to legend, mean ± SD; flow-less results P sperm:egg (α = 0 s −1 ) displayed by dashed horizontal lines in respective color). Parameters taken for A. punctuala, see Table A, except boundary radius r max = 15 mm for numerical efficiency. Chemotactic gain factor: ρ = 2, 5, 10 3.8 13.6 18.6 Table B. Relative amplitude of optimum in sperm-egg encounter probability. Ratios of encounter probability P sperm:egg (α = α * ) at optimal shear rate α * normalized by encounter probability P sperm:egg (α = 0) in the absence of flow for parameter study displayed in Fig E, Fig F, Fig G, and Fig H.