A Two-State Random Walk Model of Sperm Search on Confined Domains

1. Introduction

Mammalian reproduction depends on the ability of sperm cells to locate an egg within the complex microenvironment of the female reproductive tract. Reproductive tracts impose numerous barriers to sperm success, ultimately reducing millions of potential fertilizers to only a few dozen cells that manage to find their way to the ampulla of the oviduct where fertilization takes place. Though rudimentary forms of navigation such as chemotaxis and rheotaxis have been investigated in mammals [1,2], the physiological relevance of these phenomena remains unclear, and may be more significant for externally fertilizing organisms that navigate in three dimensions [3,4]. For nearly all eukaryotes, the sperm search for an egg depends, at least in part, on intrinsically controlled motility patterns. Additionally, among species for which forms of taxis have been shown to be important, it is likely that sperm find themselves in regions or conditions where external guidance cues may be weak or absent - highlighting the importance of understanding how intrinsically controlled motility contributes to the search process.

Though intracellular biochemical mechanisms that control sperm post-ejaculatory behaviors are well-studied in rodents and humans [5], the stochasticity of sperm search owing to complex variation in intrinsically controlled motility patterns remains poorly understood - particularly at the scale of sperm cell populations which vary from ${10}^6$ to ${10}^{12}$ cells per ejaculate depending on the species.

Sperm cells are microswimmers that are subject to the physical limitations of life at low Reynolds number [6], i.e., inertia plays no role and the motion is overdamped. Mammalian sperm, furthermore, move on the epithelial surface of the reproductive tract and 2-dimensional modeling is a reasonable simplification [7]. Sperm motility is driven by active flagellar beating and can be modeled as a random walk [8]. Contemporary computer-aided sperm motility analysis generally focuses on averaged motility patterns on timescales of just a few seconds [9].

A human sperm cell has a width of about 2 $\mu$m. Human oviducts are called “fallopian tubes.” Fertilization generally occurs in the lining of a fallopian tube. After release in the vagina, only a fraction of less than one-thousandth of the about hundred million sperm cells in a human ejaculate reaches the fallopian tubes. The lining of the fallopian tubes has folds, but can still roughly be thought of as the inside of a cylinder. The cylinder is 10 to 14 cm long and about 1 cm in diameter. The egg cell is an ovoid with a diameter of about 0.15 mm. In other words, we are looking at something the size of a pea that is looking for a big beach ball in a search area that is the size of ten football fields.

An important parameter in this study is the persistence length. The persistence length characterizes the sperm’s trajectory and it can only be observed when the trajectory is followed on a timescale longer than just a few seconds. Persistence length is the distance of over which a curve remains directionally correlated and the concept is also frequently applied in the study of polymers [10,11,12]. We will present the formal definition of persistence length in Section 2. Suffice it to say at this point that a trajectory with a long persistence length appears as a relatively straight path. A smaller persistence length, on the other hand corresponds to more frequent changes in direction. Sperm within one ejaculate can exhibit persistence lengths that may differ by orders of magnitude (cf. Figure 1).

Recently, transitions between high-persistence-length movement and low-persistence length movement have been observed in bovine sperm [13]. This suggests that sperm may transition between mere dispersal and localized search behavior. In confined or semi-confined domains, such as an oviduct in-vivo or a culture dish surface in-vitro, effective sperm navigation likely requires a strategy that balances rapid traversal across searchable surface area with thorough exploration of local space [8]. The efficiency of sperm search may depend on transition rates among movement patterns. Transitions between different movement patterns have also been observed in E. coli bacteria [14]. In that case the transitions between the “movement states” were well-described with ordinary chemical kinetics (i.e., exponentially distributed dwelling times in states).

Here, we hypothesize that this two-state random walk — switching from a high-persistence length to a low-persistence length, can enhance search efficiency in domains where uniform coverage is essential. We develop a two-state active particle model to explore this hypothesis on a disk with an absorbing boundary where we let searchers emerge from the center of the disk. In Section 2, we show through numerical simulation how changes in persistence length can improve coverage of the search domain. As a baseline, we examine the case of a searcher with a constant persistence length. Next we study two enhancements. The first one involves a well-tuned transition rate (i.e. exponentially distributed waiting times) to a persistence length that is about an order of magnitude smaller. There is evidence that the change-of-state can be catalyzed by chemicals that are present in the environment [13,14]. Taking into account this positional guidance, we also study, as a second case, a modified version of our first model: we take the stochasticity out of the switch and let the searching sperm cell change its motility state upon entering the outer section of the disk. In Section 3, we use the diffusion equation to add mathematical rigor to results of Section 2.

2. A Two-State Random Searcher on a Disk

Our geometric setup is simple. Searching particles appear in the center of a disk and follow a random trajectory until they hit the absorbing boundary at the circumference. We set the disk’s radius at $r=1$ (cf. Figure 2). Each moving searcher starts out with a large persistence length that is on the order of magnitude of the radius of the disk. The switch is to a persistence length that is about an order of magnitude smaller. We gather statistics for three cases. (i) A persistence length that is fixed at half the disk’s radius and does not change. (ii) A persistence length that is initially set at one-and-a-half times the disk’s radius, but with a transition rate that is set such that the irreversible switch to the smaller persistence length occurs on average at about $r=1/2$. (iii) Persistence lengths as in the 2nd case, but with the irreversible switch occurring at the moment that $r=1/2$ is first crossed. We consider a search effective if there is no place for the target to “hide” on the disk. In practical terms this means that all parts of the disk are equally searched. We divide the disk into annuli and numerically investigate. In the next section we will also set up and analyze a diffusion equation to help us understand the results of the stochastic simulations.

We model the sperm cell as an active particle that moves at a constant speed $v_0$. Fluctuations in the direction of motion $\theta$ are accounted for with a diffusion coefficient $D_{\theta }$. This stochasticity is due to active flagellar fluctuations, rather than Brownian motion [15]. We let ${\mathbf{e}}_{\mathbf{u}}$ be the unit vector in the direction of motion and come to the following equations describing the sperm cell’s motion:

$$\dot{\mathbf{u}}=v_0{\mathbf{e}}_{\mathbf{u}},\dot{\theta }=\sqrt{2D_{\theta }}\xi \left(t\right).$$

(1)

At any time the unit vector in the direction of motion is given by ${\mathbf{e}}_{\mathbf{u}}=\left[\cos \theta ,\sin \theta \right]$. The term $\xi \left(t\right)$ represents Gaussian-distributed white noise. We use the Euler-Maruyama scheme for the actual numerical simulations: at every timestep the angle $\theta$ is updated through the addition of a small random number. The actual $\xi \left(t\right)$ is implemented by taking a set of Gaussian distributed numbers that are zero-average and unit variance. Subsequent values of $\xi \left(t_i\right)$ are obtained by dividing these numbers by the square root of the timestep $\Delta t$.

The trajectories generated by Equation (1) display a persistence length ${\mathsf{\ell }}_{\mathrm{p}}=v_0/D_{\theta }$. Persistence length is a notion that is also well known in polymer theory [10,11,12]. We parametrize the trajectory with the contour length ℓ. Obviously, the direction $\theta$ changes as a function of ℓ. We consider the orientations $\theta \left({\mathsf{\ell }}_0\right)$ and $\theta ({\mathsf{\ell }}_0+\mathsf{\ell })$, i.e., the orientations on two points that are a distance ℓ away from each other along the contour. A good measure for the correlation between $\theta \left({\mathsf{\ell }}_0\right)$ and $\theta ({\mathsf{\ell }}_0+\mathsf{\ell })$ is the average $\langle \cos \Delta \theta \left(\mathsf{\ell }\right)\rangle$, where $\Delta \theta \left(\mathsf{\ell }\right)=\theta ({\mathsf{\ell }}_0+\mathsf{\ell })-\theta \left({\mathsf{\ell }}_0\right)$ and the average is taken over all positions ${\mathsf{\ell }}_0$ along the entire trajectory such that the different segments of length ℓ do not overlap. For a very stiff almost straight trajectory, $\langle \cos \Delta \theta \left(\mathsf{\ell }\right)\rangle$ is close to one for all ℓ. The persistence length ${\mathsf{\ell }}_{\mathrm{p}}$ is defined through $\langle \cos \Delta \theta \left(\mathsf{\ell }\right)\rangle =\exp \left[-\mathsf{\ell }/{\mathsf{\ell }}_{\mathrm{p}}\right]$ [10,11,12]. A smaller persistence length makes a trajectory appear more tightly wound (cf. Figure 2).

Related to the persistence length is the Kuhn length. Imagine a point P along the trajectory. Points that are ${\mathsf{\ell }}_{\mathrm{p}}$ to the left and right of P come with orientations that are not correlated to the orientation at P. With this realization we can next imagine the trajectory as being built up from straight segments with a length of ${\mathsf{\ell }}_{\mathrm{k}}=2{\mathsf{\ell }}_{\mathrm{p}}$. At each junction, the angle between the linked segments is a random number from a flat distribution between 0 and $2\pi$. The polymer or trajectory conceived as such is called a Freely Jointed Chain (FJC) and ${\mathsf{\ell }}_{\mathrm{k}}$ is known as the Kuhn length.

The FJC comes with concise formulae. For a chain of a total length L and with N Kuhn segments of length ${\mathsf{\ell }}_{\mathrm{k}}$, we have

$$L=N{\mathsf{\ell }}_{\mathrm{k}}\mathrm{and}\langle S^2\rangle =N{\mathsf{\ell }}_{\mathrm{k}}^2$$

(2)

for the total contour length and the average square end-to-end distance, respectively. These formulae are well-known, much applied, and explained in many authoritative textbooks, though mostly in the context of polymer theory [10,11,12].

For a searching particle that emerges at the center of the disk, it makes sense to search in such a way that any random point on the disk has an equal likelihood of being “found.” An annulus at a distance r from the center of width $\Delta r$ covers an area $2\pi r(\Delta r)$ (cf. Figure 2). A search by a moving particle is thus more likely to be successful if the dwelling time at radius r is proportional to r itself. For an effective search, the trajectory should get tighter as the rim is approached, i.e., the Kuhn length ${\mathsf{\ell }}_{\mathrm{k}}$ should decrease during the search. A simple physiological course of action to achieve such decrease of ${\mathsf{\ell }}_{\mathrm{k}}$ is to emerge from the center with an initial Kuhn length ${\mathsf{\ell }}_{\mathrm{k},1}$ and a well-tuned transition rate to a smaller value ${\mathsf{\ell }}_{\mathrm{k},2}$.

Consider the disk with radius $r=1$ as depicted in Figure 2. Figure 2a shows a search with $v_0=1$ and $D_{\theta }=2.0$, which implies a persistence length of ${\mathsf{\ell }}_{\mathrm{p}}=0.5$. Figure 2b depicts the realization of a search with an initial $D_{\theta ,1}=0.67$ and a transition rate $k=2.0$ to a final $D_{\theta ,2}=6.0$. Note that the $D_{\theta }=2.0$ of Figure 2a is the geometric average of $D_{\theta ,1}$ and $D_{\theta ,2}$.

The unit disk has three times as much area outside $r={\displaystyle \frac{1}{2}}$ as that it has inside $r={\displaystyle \frac{1}{2}}$. It thus makes sense to explore the disk’s outer part ($r>{\displaystyle \frac{1}{2}}$) with three times the contour length of the inner part. After some algebra, it follows from Equation (2) that with a Kuhn length that is one-ninth of the original one, approximately the same end-to-end distance is covered with three times the contour length. It is therefore that we take $D_{\theta ,2}=9D_{\theta ,1}$. Finally, it needs to be noted that for all generated trajectories we took $\Delta t=0.001$ for the timestep.

In Figure 2 the unit disk is divided up into nineteen annuli each with a width of 1/20. In the center there is one small circle with radius 1/20. A search is most effective if each location on the disk has an equal chance of being visited by the searcher. In terms of the annuli this means that in each annulus, the searching particle spends the same amount of time per unit of area, i.e., $\Delta T/\Delta A$ (cf. Figure 3) should be constant. Each of the three bar charts in Figure 3 is the result of 5000 trajectories that were generated with Equation (1). For each trajectory the amount of time in each annulus was recorded and then divided by the area of that annulus. Figure 3 represents the average over 5000 trajectories for three cases. For Figure 3a the radial diffusion coefficient was held constant at $D_{\theta }=2$. For Figure 3b each of the 5000 searchers starts at the origin with $D_{\theta ,1}=0.67$ and a transition rate $k=2.0$ to $D_{\theta ,2}=6.0$. Figure 3c depicts the result for searchers that can detect their own position and irreversibly switch from $D_{\theta ,1}=0.67$ to $D_{\theta ,2}=6.0$ when $r=0.5$ is first traversed.

Going from Figure 3a to Figure 3b we see an improved search. The average and standard deviation of $\Delta T/\Delta A$ over the twenty regions go from 1.1 and 1.5 to 1.7 and 1.7, which implies that the coefficient of variation (the standard deviation divided by the average) goes from 1.3 to 0.98. Furthermore, the escape time $T_{\mathrm{esc}}$, i.e. the average time it takes to go from the center of the disk to $r=1$, increases from $T_{\mathrm{esc}}=1.8$ to $T_{\mathrm{esc}}=3.1$ in the setup with the switch of $D_{\theta }$.

There is another small improvement as we go from Figure 3b to Figure 3c, i.e. go to a setup where the searcher “reads” its position and makes the switch contingent on the position. For the search in Figure 2c, the average and standard deviation are 1.6 and 1.5, resulting in a coefficient of variation of 0.94. The value of $T_{\mathrm{esc}}$ is the same as for the non-position guided switch, i.e., $T_{\mathrm{esc}}=3.1$. Particularly striking in Figure 3c is how $\Delta T/\Delta A$ stays constant in the region between about $r=0.2$ and $r=0.7$.

The improvements going from Figure 3a to Figure 3b and next to Figure 3c may appear small. However, in the 1930s the “New Evolutionary Synthesis” combined genetics and population dynamics to put evolution on a more quantitative and scientifically rigorous foundation. It was found that even small selective advantages drive evolution. The authoritative textbook by D. Futuyma [16] summarizes it as follows: “. . . a character state with even a miniscule advantage will be fixed by natural selection. Hence even very slight differences among species, in seemingly trivial characters such as the distribution of hairs on a fly or veins on a leaf, could conceivably have evolved as adaptations.”

3. The Evolution of the Probability Distribution on the Disk

The previous section involved an ordinary differential equation with a stochastic input at every timestep. The ensuing stochastic differential equation describes a trajectory of the diffusing particle. It is possible to formulate an equivalent partial differential equation that describes the time evolution of the particle’s probability distribution $P(x,y,t)$. The underlying theory is covered in many standard textbooks [17,18,19].

The equation that rules the 2D diffusion on the disk is

$${\partial }_{t}P=D\Delta P,\mathrm{where}P=P(x,y,t).$$

(3)

Here D and $\Delta$ stand for the diffusion coefficient and the Laplace operator, respectively. The diffusion coefficient D in this equation has a simple relation to the $D_{\theta }$ in Equation (1). For a particle diffusing in 2D, the mean square displacement and the elapsed time are related through $\langle \Delta r^2\rangle =4D\Delta \tau$. The $\Delta r$ can be identified with S in Equation (2). The contour length L can be identified with $v_0\Delta \tau$. This leads first to $D={\displaystyle \frac{1}{4}}{\mathsf{\ell }}_{\mathrm{k}}{v}_0$ and, after invoking ${\mathsf{\ell }}_{\mathrm{p}}={\displaystyle \frac{1}{2}}{\mathsf{\ell }}_{\mathrm{k}}$ and ${\mathsf{\ell }}_{\mathrm{p}}=v_0/D_{\theta }$, we find

$$D={\displaystyle \frac{1}{2}}{\displaystyle \frac{v_0^2}{D_{\theta }}}.$$

(4)

In polar coordinates and with circular symmetry, Equation (3) takes the form

$${\partial }_{t}P=D\left({\partial }_r^2+{\displaystyle \frac{1}{r}}{\partial }_r\right)P,\mathrm{where}P=P(r,t).$$

(5)

Before further analysis, it is worth reflecting on Equation (5) in the light of what was established in the previous section. The first term on the right-hand-side of Equation (5) is a diffusion term that does not discriminate between the directions r and $-r$. The second term on the right-hand-side is a drift term that contains a characteristic speed $D/r$. This may be puzzling as the diffusion should have no apparent direction. In his classic text, Feller [20] explains the apparent anisotropy as follows (Volume 2, Chapter 10, Section 6): “The existence of a drift away from the origin can be understood if one considers a plane Brownian motion starting at the point r of the x-axis. For reasons of symmetry its abscissa at epoch $h>0$ is equally likely to be $>r$ or $<r$. In the first case certainly $R\left(h\right)>r$, but this relation can occur also in the second case. Thus the relation $R\left(h\right)>r$ has probability $>{\displaystyle \frac{1}{2}}$, and on the average R is bound to increase.”

With an initial condition $P(r,0)=\delta \left(r\right)/\left(2\pi r\right)$ and on a disk with infinite radius, Equation (5) is readily solved by the spreading Gaussian:

$$P(r,t)={\displaystyle \frac{1}{4\pi Dt}}\exp \left[-{\displaystyle \frac{r^2}{4Dt}}\right].$$

(6)

This solution comes with a root-mean-square distance to the origin, $\tilde{r}$, that increases following ${\tilde{r}}^2=4Dt$. Taking the derivative w.r.t. time on both sides, we find for the speed of the root-mean-square distance $\dot{\tilde{r}}=2D/\tilde{r}$.

For a particle that drifts away from the origin with a speed $v=2D/r$, the time that is spent in an annulus of area $\Delta A=2\pi r\Delta r$, is found to be $\Delta T=(r/2D)\Delta r$. So the search time per area, $\Delta T/\Delta A=1/\left(4\pi D\right)$, appears to be independent of r. This uniform coverage of the search area occurs with any drift speed that is proportional to $1/r$.

However, it is not hard to see that such r-independence of $\Delta T/\Delta A$ only works for a disk with a radius $R\to \infty$. If a diffusing particle is absorbed upon hitting the boundary at a finite $r=R$, any trajectory leading to a position $r<R$ after passing through an $r>R$ region is eliminated. This leads to a decreasing $\Delta T/\Delta A$ as r increases. Such decrease is indeed what we see in Figure 3. Below we will add mathematical rigor to these insights and observations.

An absorbing boundary at $r=R$ means $P(R,t)=0$. Solutions for setups with an absorbing boundary condition can commonly be found through placing an “image initial condition,” $\overline{P}(r,t=0)$, outside the disk. The image is negative and the positioning of the image initial conditions has to be such that the ensuing flow towards $r=R$ (governed also by Equation (5)) leads to $\overline{P}(R,t)+P(R,t)=0$. The full solution $(P+\overline{P})(r,t)$ then satisfies Equation (5) in the area $r<R$.

For 1D diffusion along an axis, the image method is easily implemented and readily leads to a solution. But complications arise in the 2D case. As was mentioned before in the context of Equation (5), the directions towards the disk’s center and away from the disk’s center are not equivalent. For an image that is initially positioned at $r=2R$, i.e. $\overline{P}(r,t=0)=-\delta (r-2R)/\left(4\pi R\right)$, a solution according to Equation (5) can be derived after some nontrivial mathematical operations [21]: $\overline{P}(r,t)=-1/\left(4\pi Dt\right)I_0\left(rR/\left(Dt\right)\right)\exp \left[-(r^2+4R^2)/\left(4Dt\right)\right]$, where $I_0$ represents the Modified Bessel function of the first kind. The problem with this solution is that we do not have $\overline{P}(R,t)+P(R,t)=0$ at all times.

Below we derive some relevant features of the diffusion on the 2D disk without an actual closed-form solution. Through Fick’s Law, the outflow from the disk can be identified with $J(R,t)=-D\nabla P(R,t)=-D{\partial }_{r}P(R,t)$. If we set ${\partial }_{r}P(R,t)=-\gamma$ for the slope of the $P(r,t)$ curve at $r=R$ (cf. Figure 4), then we have $J(r,t)=D\gamma$ at $r=R$. $J(r,t)=D\gamma$ is also a good approximation on $r<R$ near the outer edge of the disk.

In our 2D case, the flux $J(r_0,t)$ is the amount of probability that crosses at $r=r_0$ per unit of time and per unit of distance along the circumference of the circle with radius $r_0$. Next it is the ratio $J(r_0,t)/P(r_0,t)$ that represents the average speed that is necessary to bring about this flux. With $P(r,t)=\gamma (R-r)$, we then obtain $v\left(r\right)=D/(R-r)$ for the drift speed of the diffusing particle. We observe that, unlike with the $R\to \infty$ case discussed earlier, the drift speed now goes up as the absorbing boundary of the disk gets closer. This can be readily intuited upon the realization that with $r=R$ more nearby, ever more possible trajectories hit $r\ge R$ and disappear from the domain.

As before, we divide up our disk in annuli of a width $\Delta r$. With the drift speed that we derived in the last paragraph, we find for the dwelling time per annulus $\Delta T/\Delta r=(R-r)/D$. From there we derive a search time per area that depends on r:

$${\displaystyle \frac{\Delta T}{\Delta A}}\left(r\right)={\displaystyle \frac{R-r}{2\pi Dr}}.$$

(7)

It is obvious that $\Delta T/\Delta A$ vanishes as $r\to R$. Figure 5 shows the results of a simulation. We let 1000 particles diffuse from the center of the disk to the absorbing boundary at $r=1$. The motion of the particles is again ruled by Equation (1) with $v_0=1$ and $D_{\theta }=200$. Through Equation (4) this leads to $D=2.5\times {10}^{-3}$. The timestep used for the Euler scheme was again $\Delta t=0.001$. As before, the unit disk was segmented into a small $r=0.05$ disk and nineteen annuli of a width of 0.05. Like before, the twenty bars in the histogram (Figure 5a) give the average time per unit area spent in each of the twenty segments. Figure 5b gives $\Delta T/\Delta A$ for the four outermost annuli together with the best fitting straight line. Equation (7) gives 63.7 for this slope. The slope of the fitted line is 67.7. The small discrepancy is not just due to the neglected curvature, but also to the finite persistence length that is implicit in simulating Equation (1). In our case ${\mathsf{\ell }}_{\mathrm{p}}=5\times {10}^{-3}$. The partial differential equations that are central in this section, i.e., Equations (3) and (5), actually correspond to a fractal trajectory with a zero persistence length.

A diffusion coefficient that decreases as the particle diffuses away from the center leads to a D that effectively decreases with r. As was pointed out in the previous section and as now also should be clear from Equation (7), this flattens out $\Delta T/\Delta A$ as a function of r and leads to a more efficient search of the disk’s area.

Finally, it is worth pointing out that Equation (7) and the drift speed $v\left(r\right)=D/(R-r)$ are independent of the actual slope $-\gamma$. This is a nice feature, because as the particle diffuses away from $r=0$, that slope actually changes over time. Immediately after $t=0$, the slope increases. Later in the relaxation, the slope starts decreasing again. It is remarkable that Equation (7) and the drift speed $v\left(r\right)$ remain the same throughout that process.

4. Discussion

We have studied a simple, bare-bone model to keep the analysis intuitive and to have manageable mathematics and numerics. Figure 6 presents a model with more realism. In the oviduct the searching sperm cells enter on one side. Here we have them enter at an arbitrary point on the left edge. We account for being on the inside of a cylinder through a periodicity of the vertical coordinate. The two-state strategy that we studied on the disk would also likely be effective in this case. A good dispersal of searchers could be achieved by taking off with a large persistence length in an arbitrary direction from an arbitrary point on the left edge. Each individual sperm cell would next do its more detailed search after it switches to the smaller persistence length.

Some of the trajectories in Figure 1 exhibit chirality. In our minimal model we did not not consider chirality. The chirality could readily be added to the model in Equation (1) through the addition of an angular velocity $\omega$ on the right-hand-side of the second equation. Such study was undertaken in Ref. [13]. It is obvious that chirality, or any curvature for that matter, is like a smaller Kuhn length in that it leads to an area being more thoroughly searched.

Routine semen analysis in clinical practice relies on microscopic inspection of sperm characteristics and statistical comparison with fertile samples [22,23]. Contrary to what is commonly thought, fast motility alone is not necessarily an indicator of high fertilizing potential due to the complexity of the sperm search task and statistical variation among the large number of sperm present in an ejaculate [24]. Efficient search likely requires a combination of motility patterns including high-persistence-length movement for rapid dispersal, and low-persistence-length movement for thorough local exploration. Our modeling suggests that a well-tuned combination of both modes may enhance search efficiency on restricted domains such as those that occur in the complex environment of the female reproductive tract, where sperm must balance wide-area coverage with detailed exploration to locate the egg.

The two-state random walkers (Figure 2b and Figure 3b,c) that we studied represent improvements over the regular random walk (Figure 2a and Figure 3a). In closing, it is instructive to look at our subject matter and results from the wider perspective of computation and problem solving in biological systems. Searches are also part of macroscopic everyday reality. When a search-and-rescue team goes looking for a lost or missing person in a large wilderness area, a common strategy is to cut up the search area. The smaller areas get labeled and next different teams or individuals are assigned to different areas in the grid for detailed searches. Obviously, sperm cells cannot achieve such a level of organization and coordination.

The random walk with identical searchers with constant persistence length (cf. Figure 2a) is the most brute-force approach to the search problem. It is the Monte-Carlo modus operandi and it requires the mobilization of many searchers to achieve a nonegligible probability of fertilization. The random walk of Figure 2b and Figure 3b contains a transition to a smaller persistence length. We have shown that this feature leads to improved efficacy. Through the transition rate k, this search involves a characteristic time $1/k$. So the improvement is due to the usage of a very crude clock and the search being split up into a “before” and “after.” Figure 3c is the barchart that corresponds to a random-walk searcher that does not have a clock, but instead uses a very crude GPS. The searcher only “knows” whether it is in Region 1 or 2, and switches persistence length the first time it “sets foot” in Region 2.

Biologically, a well-tuned k or an effective partitioning into two regions would be an outcome of natural selection. This improvement over the pure Monte-Carlo algorithm can be viewed as enhanced agency and as something “learned” - a genetic legacy built up after much past “trial and error.” Neither of the two improvements involve any sperm-to-sperm communication.

The analogy with a team on a search-and-rescue mission underscores an important point. Popular culture often views ejaculate as a swimming contest. However, it is more appropriate to see ejaculate as tissue. A heart pumps, an ear to detect sound, and ejaculate searches. Successful ejaculate is ejaculate that yields a high probability that one sperm cell finds and fertilizes an egg cell. It is obvious that there is selectional pressure towards more successful ejaculate.