Moth Mating: Modeling Female Pheromone Calling and Male Navigational Strategies to Optimize Reproductive Success

Abstract

Male and female moths communicate in complex ways to search for and to select a mate. In a process termed calling, females emit small quantities of pheromones, generating plumes that spread in the environment. Males detect the plume through their antennae and navigate toward the female. The reproductive process is marked by female choice and male–male competition, since multiple males aim to reach the female but only the first can mate with her. This provides an opportunity for female selection on male traits such as chemosensitivity to pheromone molecules and mobility. We develop a mathematical framework to investigate the overall mating likelihood, the mean first arrival time, and the quality of the first male to reach the female for four experimentally observed female calling strategies unfolding over a typical one-week mating period. We present both analytical solutions of a simplified model as well as results from agent-based numerical simulations. Our findings suggest that, by adjusting call times and the amount of released pheromone, females can optimize the mating process. In particular, shorter calling times and lower pheromone titers at onset of the mating period that gradually increase over time allow females to aim for higher-quality males while still ensuring that mating occurs by the end of the mating period.

1. Introduction

Moths, of the order Lepidoptera, are one of the largest groups in the animal kingdom, with nearly 160,000 species tallied in almost all land habitats. Many more have yet to be described. These insects are highly adaptable to the local environment and develop close relationships with local food sources and other organisms; as a result, they are highly diverged in morphology and biology and display extraordinarily diverse colors, shapes, and sizes. The lifetime of a moth is divided into four periods: the egg, larva, pupa, and adult stages. The length of these periods depends on the specific species, habitat, climate, and available food sources, but in most cases, the adult stage varies between a few days (the yucca and silkworm moths) to a few months (the common brown house moth). A moth’s mating period mostly coincides with the adult stage (Lees and Zilli [1]).

The complex process of searching for a mate for copulation involves male and female moths communicating in complementary ways. Although some crepuscular and diurnal species have been identified, most moths are nocturnal. The predominant communication mechanism of nocturnal species is via the release of pheromones by females, termed calling, and the subsequent sensing of these chemical cues by males (Naka [2], Sarto i Monteys et al. [3]).

A continuous period of calling is known as a call. During a call, female moths emit pheromone molecules from dedicated glands at the end of their abdomens via short, discrete pulses called puffs that are dispersed in the atmosphere through convective and diffusive motion. Collectively, these puffs define a plume. Typical emission rates are on the order of a few puffs per second (Conner et al. [4]). Male moths detect the female-generated pheromone plumes through sensilla, sensory receptors located on their antennae and follow the scent toward the source (Stengl [5], Futrelle [6], Groot [7]). Due to the great olfactory sensitivity of male antennae, adaptive species-specific pheromone blends can be distinguished from other background, food, or plant odors even at low concentrations, and males are able to follow trails over long distances, up to several kilometers (Badeke et al. [8], Leary et al. [9]).

The reproductive process is marked by female choice and male–male competition, since typically multiple males aim to reach the female but only the first can mate with her (Greenfield [10]). This translates into a high opportunity for sexual selection on desirable male traits, such as those regarding body weight, antenna structure, wing length, genital clasper length, and symmetry (Koshio et al. [11]). However, maximizing these traits can come at a cost. For example, while larger antennae may be able to detect pheromone more efficiently, there will be more aerodynamic drag during flight and energy limitations due to more neural tissue (Niven et al. [12], Sane [13], Stockl et al. [14], Vogel et al. [15], Wang et al. [16]).

Females utilize several pheromone-producing strategies to maximize benefits and to minimize costs. If too little pheromone is released, the likelihood of being reached by a male is low and the call may result in no mating. Conversely, if the released pheromone is too abundant, then less fit males will be able to reach the female, resulting in low-quality offspring. Female moths thus adapt several strategies such as changing the onset, duration, frequency, and intensity of the mating call, though there is a lack of quantitative experimental data on pheromone production fitness and physiological costs (Umbers et al. [17], Harari et al. [18]).

Understanding how pheromones are secreted and followed may be useful in pest control, since ad hoc chemical blends may be synthesized to attract crop-damaging moths and may disrupt the mating process. From a neuroscience perspective, moths represent a useful model system for studying olfactory information processing, learning, and memory (Haupt et al. [19], Martinez et al. [20], Bau and Cardé [21], Grünbaum and Willis [22], Baker et al. [23], Loveless and Webb [24]). From an engineering perspective, moths may serve as a paradigm for odor source localization over long distances while subject to large perturbations from a highly dynamic environment, guiding the design of prototype plume tracing, hybrid robot–insect machines, and algorithms for the control of autonomous vehicles (Edwards et al. [25], Bau and Cardé [21], Ando et al. [26], Huang et al. [27], Liberzon et al. [28], Zhu et al. [29]).

While there are many variations of male navigational patterns, the key features are well understood: male moths search randomly until they encounter female-generated and wind-advected pheromone plumes, after which they surge, i.e., fly directly upwind, toward the odor source: the female. If the scent is lost, they undergo a meandering motion where the male adopts a zigzag trajectory perpendicular to the direction of the wind, termed casting, to reenter the plume (Kennedy [30], Mafra-Neto and Cardé [31], Cardé and Willis [32], Bau and Cardé [21], Liberzon et al. [28], Benelli et al. [33]). However, female pheromone release strategies are more diverse and vary in call length, amount and pattern of pheromone released, and onset time (Umbers et al. [17], Harari et al. [18]).

It is worth noting that, although the process described above is the most common, in some polyandrous species, females mate more than once, with the last male to arrive physically displacing any previously deposited spermatophore (Xu and Wang [34]). Female/male selectivity may also change dynamically, depending on seasonal effects or trade-offs between reproductive benefits and costs (Gao et al. [35]). Finally, in species where male reproductive efforts are costly, it may be the male that chooses which female to mate with. For example, Achroia grisella male moths select for more fecund, heavier females due to sperm production limitations (Goubault and Burlaud [36]).

To understand how female pheromone-release strategies affect mating, we perform large-scale agent-based simulations where a stationary female varies the characteristics of her calling strategies while male moths navigate toward her. More specifically, females will vary call length (the amount of time per night spent calling males) and the emitted pheromone titer (the amount of pheromone released per call) over the course of a typical one-week mating period. The dynamic pheromone plume will then attract male moths depending on their detection capability.

First, we present a simplified agent-based model where the female generated plume is assumed to be static and circular and where the male search is dominated by the time to first reach the stationary plume. We call this model the Male Random Flight Model since males perform a two-dimensional Brownian motion search until the plume is reached. Although many physical traits can be associated with male quality, in this work, we focus on two primary markers: flight speed/diffusivity and ability to sense small concentrations of pheromone, i.e., chemosensitivity. These two traits emerge as a combination of more specific body, antenna, sensilla, or wing characteristics from which females choose through sexual selection. We interpret these two traits as collective proxies and study how varying them affects mating likelihood. Second, we consider a more complete agent-based model where the female-generated plume is more realistically advected by wind currents and displays time-dependent behavior, that we have named the Plume Navigation Model. We use mathematical analysis to determine the mating probability and arrival time distribution in the Male Random Flight Model under various female calling strategies. We eschew some details, such as pheromone plume dynamics, in order to assess how male diffusivitiy and chemosensitivity affect the mating process during a single night and over the female’s lifespan. In the Plume Navigation Model (that does not permit mathematical analysis due to its complexity), we use extensive numerical simulations to explore how the complex interplay between female-produced pheromone plume dynamics, the male search for the plume, and subsequent male plume navigation affect the mating process. The source code for all numerical simulations in this paper is provided as Supplementary Materials.

Two broad findings emerge from our study. First, we observe that a female can select for male quality while ensuring that mating probability remains high by modifying her calling strategy: if unable to attract high-quality mates at onset of the mating period, she can then shift toward less restrictive calling at a later time. Second, we find that the average first mating time is significantly shorter than the mean arrival time, over large parameter ranges and in both models, confirming that pheromone-based signaling allows females to attract males of differing quality: the highest-quality males are able to reach her quickly, and less optimal ones reach her at later times. Our agent-based simulations also reveal that mating likelihood is less sensitive to pheromone titer amounts than to calling times. Females that release less pheromone preferentially select for males with lower detection thresholds, a trait that we can associate with longer or more elaborate antennae. This outcome is in agreement with field observations (Johnson et al. [37]).

An overview of the paper is as follows: in Section 2, we review female and male mating strategies; in Section 3, we introduce the Male Random Flight Model and present some analytical estimates. Related numerical results are discussed in Section 4. In Section 5, we present the Plume Navigation Model of moth mating with numerical results. A broader discussion and conclusions are presented in Section 6.

2. Moth Mating Mechanisms

In this section, we review experimental studies describing the calling methods employed by females (Section 2.1) and the navigational strategies used by males (Section 2.2) as they seek optimal mating outcomes. These strategies will be implemented in our simplified model (Section 3 and Section 4) and in our full agent-based simulation (Section 5).

2.1. Female Calling Strategies

The release of pheromone by female moths is known as calling and is characterized by the female raising her abdomen to better disperse pheromone from the scent gland. The subsequent travel of pheromone molecules through the air forms odor plumes that are dispersed by the wind (Murlis et al. [38]). For nocturnal species, calling occurs during scotophase, the dark phase of the light–dark cycle, although the timing of calling onset, duration, and pheromone titer varies with species. Adult moths attempt to mate over several days before perishing, depending on their life span. In this section, we identify patterns of female calling and determine ranges of values for calling duration and pheromone titer. We assume a one-week mating period, as typical (Groot [7], Lees and Zilli [1]).

The meta-analysis of Umbers et al. [17] reviewed 52 empirical studies of 44 Lepidoptera species and categorized the reported female calling patterns of each. The calling effort was quantified as the time females spent calling males and/or as the amount of female pheromone titer. While most of the studies measured titer as the daily amount of chemicals present in female pheromone glands, a handful of others considered the amount of released pheromone per night. We do not distinguish between the two classifications and pool all reported values into one range. We thus implicitly assume that the amount of pheromone stored in the female body per day is the same as the cumulative amount of pheromone released during her nocturnal call. Note that it is experimentally difficult to measure pheromone ambient concentration and to infer its exact relationship to pheromone gland content (Umbers et al. [17]). Furthermore, there is great variability in the amount of pheromone stored and released by females both within a single species and across species and over the spatial extent of a pheromone plume (Mankin and Mayer [39]).

Four possible categories of female calling patterns were considered in Umbers et al. [17]: constant (the calling effort is the same every night), increase (it is greater each night compared to the night before), decrease (it is less each night compared to the night before), or hat (it increases to a peak amount and then decreases). Among those studies that measured the time that females spent calling, approximately 70% reported that the effort followed the increase pattern, 15% reported the hat pattern, 9% reported the decrease pattern, and 6% reported the constant pattern. For studies that measured pheromone titer, approximately 40% reported that the effort followed the hat pattern, 34% reported the decrease pattern, 16% reported the increase pattern, and 10% reported the constant pattern. Due to variability, our models account for changes in both the time spent calling and the pheromone titer. A schematic is shown in Figure 1.

Since our model aims to give a general picture of the dynamics of Lepidoptera moth reproduction without focusing on a single species, we group all relevant estimates reviewed in Umbers et al. [17] for calling duration (Kou and Chow [40], Mazor and Dunkelblum [41], Ming et al. [42]) and pheromone titer (Mazor and Dunkelblum [41], Ming et al. [42], Giebultowicz et al. [43], Tang et al. [44], Foster et al. [45]), to estimate the following:

• calling duration: ranges from 0.5 to 8 h and

• emitted pheromone: ranges from 1 to 30 ng per night.

Farrell et al. [46] estimated a pheromone release rate of 14 ng/h, corresponding to approximately 4 pg/s. As discussed above, we use these broad ranges as guiding estimates in our simulations but also consider larger values to capture the variation among species.

2.2. Male Navigational Strategies

Male moths from diverse species exhibit a set of well-known navigational behaviors in the presence of female pheromone plumes (Kennedy [30], Mafra-Neto and Cardé [31], Cardé and Willis [32], Bau and Cardé [21], Liberzon et al. [28], Benelli et al. [33]). While vertical motion is observed in the presence of trees and thick canopies and under conditions of high moth density (Girling et al. [47]), in most unobstructed cases, males stabilize their flight at a fixed altitude and motion along the vertical axis is limited (Cardé et al. [48], Elkinton et al. [49]). Male moths typically execute planar random walk flights during which their sensilla gather environmental information. Upon detection of a species-specific plume by antennal chemoreceptors, male moths orient upwind using antennal mechanoreceptors, and subsequently cease random flight and surge, flying directly upwind, towards the source.

Since turbulence of the ambient air makes pheromone plumes intermittent and filamentous, a male moth may lose contact with the plume as it surges toward the female (Martinez et al. [20]). Upon loss of pheromone signal, the moth adopts a zigzag trajectory perpendicular to the direction of wind flow at a fixed altitude. This behavior is called casting and is thought to increase the probability that the male moth would recapture the pheromone plume (Willis et al. [50]). If casting is unsuccessful for a given period of time, the male moth resumes random flight.

Estimates reported in the literature (Bau and Cardé [21], Belanger and Arabas [51], Kuenen and Cardé [52], Cardé et al. [48], Gaydecki [53]) lead us to the following male flight speed range and turning rate:

• male flight speed: ranges from 0.5 to 5 m/s

• male turning rate: ranges from 3 to 4 turns/s

Although flight speed and turning rates depend on whether male moths are engaged in a random flight, surging, or casting behavior, we assume that males broadly move and turn within the above ranges at all times. These values lead us to estimate typical random walk diffusivities $D\sim 1$ to ${10}^{-2}$ m²/s, as also confirmed by other studies on Lepidoptera diffusivity estimated at $D\sim 1$ to ${10}^{-1}$ m²/s (Watanabe [54], Kareiva [55].) For transitions between flight modes, the time necessary for males to measure pheromone concentrations is less than 100 ms, while the time to realize that the pheromone plume is lost is estimated as 300–500 ms (Stengl [5], de Bruyne and Baker [56], Baker and Vogt [57]). Since the typical time spent in any flight mode varies between seconds and hours, several orders of magnitude larger than the above awareness/decision times, we neglect the time lag between transitions. Finally, studies on pheromone-specific receptor neurons reveal that male moth sensilla are able to detect pheromones at very low concentrations, allowing for males to detect the plume at large distances (Angioy et al. [58], Tabuchi et al. [59]). If we assume a typical plume length of 2 km and estimate male flight speed at 1 m/s, we can also give the rough estimate of 30 min as the time a male moth spends in the female-generated plume. Most nocturnal male moths begin flying in the evening and until the predawn hours, indicating that male activity normally occurs away from the plume or when searching for it (Nowinszky et al. [60]).

3. Male Random Flight Model

We begin with modeling the male moth random search for the female-generated pheromone plume. This “foraging” stage of the mating process has often been considered independent of the surging/casting navigation processes, and few theoretical studies have focused on it (Cardé and Willis [32], Bau and Cardé [21], Sabelis and Schipper [61], Dusenbery [62,63]). The arrival time $t_{\mathrm{a}}$ to first reach the plume however is important, since males may spend a significant portion of the night in the foraging state and since shorter arrival times are indicative of moths that are better suited for mating. Although capturing the plume necessitates finding the female through surging/casting, for simplicity, we assume here that $t_{\mathrm{a}}$ is much larger than the time spent in the plume and study how the capture (and mating) probability depends on male traits relevant to random search, such as detection thresholds and diffusivity.

We utilize a two-dimensional setting with female and male moths located on the same plane at a fixed altitude based on the observed lack of male vertical motion in experimental studies, as discussed in Section 2.2. We also model the concentration of female-emitted plume $c(\mathbf{x},t)$ as a two-dimensional field. This quantity varies in space and time, and it is a function of wind currents and female emission activity. A male moth is assumed to be able to detect the plume if $c(\mathbf{x},t)$ is larger than a given threshold $C_{\mathrm{tol}}$, a value that indicates male chemosensitivity to pheromone and is a proxy for the length and/or elaborateness a male’s antenna, where longer/more elaborate antennae are able to sense smaller quantities of pheromone (Symonds et al. [64]).

The spatial region where $c(\mathbf{x},t)>C_{\mathrm{tol}}$ defines the trapping region $\Omega \left(t\right)$:

$\Omega \left(t\right)=\{\mathbf{x}:c(\mathbf{x},t)\ge C_{\mathrm{tol}}\}.$

Once the plume is detected, the male moth will cease random flight and will surge upwind toward the female. Since lower detection thresholds $C_{\mathrm{tol}}$ can be used as a proxy for male moths that are able to detect smaller amounts of pheromone and, thus, are higher quality, the trapping region $\Omega \left(t\right)$ is larger for higher quality males. Moths located outside the trapping region are not affected by the plume and are assumed to engage in a two-dimensional random flight with diffusion constant D, which serves as a proxy for flight speed. Similar to $C_{\mathrm{tol}}$, D can be considered an indicator of male moth quality: moths with larger D are more diffusive and are able to locate the plume earlier than moths with lower D.

We investigate how $C_{\mathrm{tol}}$ and D influence $t_{\mathrm{a}}$, especially when $C_{\mathrm{tol}}$ and D are both relatively large or when both are relatively small. These are interesting biological scenarios, since moths with longer and/or more elaborate antennae (smaller thresholds $C_{\mathrm{tol}}$) are associated to larger trapping regions but are also greater in body size and less motile (smaller diffusion constants D) (Symonds et al. [64], Wang et al. [16], Johnson et al. [37]).

The probability density function $g(\mathbf{x},t)$ that a moth is at $\mathbf{x}=(x,y)$ outside the trapping region $\Omega \left(t\right)$ at time t and undergoing random motion satisfies the exterior parabolic problem:

(1a)$\begin{array}{cc}\frac{\partial g}{\partial t}=D{\nabla }^{2}g,\hfill & \hfill \mathbf{x}\in {\mathbb{R}}^2\backslash \Omega \left(t\right),t>0,\end{array}$

(1b)$\begin{array}{cc}g=0,\hfill & \hfill \mathbf{x}\in \partial \Omega \left(t\right),t>0,\end{array}$

(1c)$\begin{array}{ccc}g=\delta (x-x_0)\delta (y-y_0),\hfill & & \hfill \mathbf{x}\in {\mathbb{R}}^2\backslash \Omega \left(t\right),t=0,\end{array}$

where ${\mathbf{x}}_0=(x_0,y_0)$ is the initial moth location, $\Omega \left(t\right)$ is the dynamic trapping region defined by the time evolution of $c(\mathbf{x},t)$, and $\partial \Omega \left(t\right)$ is its boundary. Equation (1c) implies that the moth is initially located outside the trapping region, in ${\mathbb{R}}^2\backslash \Omega \left(t\right)$. We use absorbing boundary conditions in Equation (1b) to indicate that random flight ceases once the moth has reached the plume’s edge. At this point, the male moth will navigate toward the female using well-known surging and casting behavior. The overall probability that the moth has not captured the plume at time t is thus given by

(2)$P\left(t\right)={\int }_{{\mathbb{R}}^2\backslash \Omega \left(t\right)}g(\mathbf{x},t)dxdy,$

while the overall probability that the moth has captured the plume at time t is defined as $C\left(t\right)=1-P\left(t\right)$. Since we equate reaching the plume with mating success, $C\left(t\right)$ can also be referred to as the moth mating probability.

One quantity of interest is the arrival time distribution $S\left(t\right)$ of a moth to the trapping region (Chou and D’Orsogna [65]), which is determined by calculating the total flux:

(3)$S\left(t\right)=-\frac{dP}{dt}=-\frac{d}{dt}\left({\int }_{{\mathbb{R}}^2\backslash \Omega \left(t\right)}g(\mathbf{x},t)dxdy\right).$

Using standard differentiation tools for the moving surface in Equation (1a) and the divergence theorem, we write the following:

(4)$S\left(t\right)=-\left(D{\int }_{\partial \Omega \left(t\right)}\nabla g·\mathbf{n}d\ell +{\int }_{\partial \Omega \left(t\right)}({\mathbf{v}}^{\mathrm{b}}·\mathbf{n})gd\ell \right),$

where $\mathbf{n}$ is the outward facing normal vector to the $\partial \Omega \left(t\right)$ boundary contour, ${\mathbf{v}}^{\mathrm{b}}$ is the velocity of the boundary contour element, and $d\ell$ is the arc length along $\partial \Omega \left(t\right)$. Equation (4) implies that the extent of trapping region $\Omega \left(t\right)$ is time-dependent and that its changing boundaries define non-isotropic shapes.

To follow the time evolution of $\Omega \left(t\right)$ we need to specify the dynamics of $c(\mathbf{x},t)$, which will be done in Section 5. Here, instead, we introduce some simplifications to Equation (4) which will allow for a more basic yet fundamental understanding of the mating process. We thus assume that $\Omega \left(t\right)$ is static and circular, centered at the origin, and characterized by a radius $r_{\mathrm{mate}}$ so that the first term in Equation (4) is an integral over a circumference and the second term vanishes. We associate higher quality moths with a larger $r_{\mathrm{mate}}$ (and larger D) and lower quality moths with a smaller $r_{\mathrm{mate}}$ (and smaller D).

In this scenario, closed form expressions for the quantities of $P\left(t\right)$ and $S\left(t\right)$ can be calculated (Wendel [66], Carslaw and Jaeger [67]). Detailed derivations are given in Appendix A. We find

(5a)$\begin{array}{cc}\hfill S\left(t\right)=& \frac{2D}{\pi r_{\mathrm{mate}}^2}{\int }_{\omega =0}^{\infty }\left[\frac{J_0\left(\omega \right)Y_0\left(\omega R^{-1}\right)-J_0\left(\omega R^{-1}\right)Y_0\left(\omega \right)}{Y_0^2\left(\omega \right)+J_0^2\left(\omega \right)}\right]\omega e^{-{\omega }^{2}Dt/r_{\mathrm{mate}}^2}d\omega ,\hfill \end{array}$

(5b)$\begin{array}{cc}\hfill P\left(t\right)=& \frac{2}{\pi }{\int }_{\omega =0}^{\infty }\left[\frac{J_0\left(\omega \right)Y_0\left(\omega R^{-1}\right)-J_0\left(\omega R^{-1}\right)Y_0\left(\omega \right)}{Y_0^2\left(\omega \right)+J_0^2\left(\omega \right)}\right]\frac{1}{\omega }{e}^{-{\omega }^{2}Dt/r_{\mathrm{mate}}^2}d\omega ,\hfill \end{array}$

where $J_0\left(z\right)$ and $Y_0\left(z\right)$ are Bessel functions of the first and second kinds, respectively, and the dimensionless R is defined as

(6)$R=\frac{r_{\mathrm{mate}}}{\sqrt{x_0^2+y_0^2}}=\frac{r_{\mathrm{mate}}}{r_0},$

and where $r_0$ is the male moth’s initial distance from the female. If we assume that the male moth is initially outside the circular plume, $r_0\ge r_{\mathrm{mate}}$, Equation (6) implies that $R\le 1$. We illustrate some characteristic forms for $P\left(t\right),C\left(t\right)$, and $S\left(t\right)$ in Figure 2A and typical moth trajectories obtained through numerical simulations in Figure 2B.

First Arrival and Mating Time

The capture/mating probability $C\left(t\right)=1-P\left(t\right)$ that an individual male moth reaches a female by time t and the corresponding arrival time distribution $S\left(t\right)$, evaluated via Equation (5a,b), are important overall metrics in determining the mating success of a single moth. However, competition arises when multiple males attempt to reach the same female and only the first one can mate. In this section, we consider N competing moths and adapt some known results on first arrival times (Lawley [68]) to give analytical estimates for the average time it takes the first male to reach the female. These results are valid in the $N\gg 1$ limit.

In accordance with our previous discussion, we identify the male’s first arrival time to the female-generated plume as the first mating time. Lawley [68] analyzes the probability $P_{\mathrm{a}}(t_{\mathrm{a}}<t)$ that the first arrival time $t_{\mathrm{a}}$ out of a cohort of N independent random walkers is less than t and shows that if $P_{\mathrm{a}}(t_{\mathrm{a}}<t)$ can be approximated in the short time limit as

(7)$P_{\mathrm{a}}(t_{\mathrm{a}}<t)\sim A{t}^q{e}^{-\frac{B}{t}},t\to 0^{+},$

for a given q, $A\ne 0$, and $B>0$, then the first arrival time $t_{\mathrm{a}}$ is distributed according to Gumbel distribution:

(8)$P_{\mathrm{a}}\left(t_{\mathrm{a}}\right)=\frac{1}{a_N}exp\left[\frac{t_{\mathrm{a}}-b_N}{a_N}-exp\left(\frac{t_{\mathrm{a}}-b_N}{a_N}\right)\right],$

where

(9)$a_N=\frac{b_N}{q(1+W_N)},b_N=\frac{B}{q{W}_N},W_N=W_{*}\left(\frac{B}{q}{\left(AN\right)}^{\frac{1}{q}}\right),$

and where $W_{*}\left(z\right)$ is the principal branch of the LambertW function, defined as the inverse function of $f\left(z\right)=z{e}^z$. The first two moments of the Gumbel distribution are

(10)$\mathbb{E}\left[t_{\mathrm{a}}\right]=b_N-{\gamma }_e{a}_N,\mathrm{Var}\left[t_{\mathrm{a}}\right]=\frac{{\pi }^2}{6}{a}_N^2,$

where ${\gamma }_e\approx 0.5772$ is the Euler–Mascheroni constant.

We can use these results by noting that, within our model, $P_{\mathrm{a}}(t_a<t)$ can be identified with the capture/mating probability $C\left(t\right)$: if the first arrival time is less than t, then at time t, the moth will have reached the plume. In Equation (A20) of Appendix B, we derive the short time behavior for $C\left(t\right)$ and verify that it follows Equation (7) as $t\to 0^{+}$. We find

(11)$C\left(t\right)\sim A{t}^q{e}^{-\frac{B}{t}},t\to 0^{+},$

with

(12)$q=\frac{1}{2},A=\frac{2}{(1-R)r_{\mathrm{mate}}}\sqrt{\frac{DR}{\pi }},B=\frac{r_0^2{(1-R)}^2}{4D},$

from which we conclude that the average first arrival time $t_{\mathrm{a}}$ for $N\to \infty$ competing male moths is given by

(13a)$\mathbb{E}\left[t_{\mathrm{a}}\right]\sim \frac{r_0^2{(1-R)}^2}{2D{W}_N}\left[1-\frac{2{\gamma }_e}{1+W_N}\right],$

where

(13b)$W_N=W_{*}\left(\frac{2N^{2}R}{\pi }\right)\sim lnN.$

The mean first arrival time $\mathbb{E}\left[t_{\mathrm{a}}\right]$ is thus a decreasing function of D and a quadratically increasing function of $r_0(1-R)={(x_0^2+y_0^2)}^{1/2}-r_{\mathrm{mate}}$, as can be expected from random searchers in two dimensions. Also note that, since $W_N\sim lnN$, in the $N\to \infty$ limit, $\mathbb{E}\left[t_{\mathrm{a}}\right]$ scales as ${(lnN)}^{-1}$. Finally, note that these estimates assume that the time frame for male moths to reach the female is unbounded, so that all values of $t_{\mathrm{a}}$, including large ones, contribute to $\mathbb{E}\left[t_{\mathrm{a}}\right]$ as evaluated in Equation (13a). In practice, the female call may last for a limited period $t_{\mathrm{call}}$, so that arrival times $t_{\mathrm{a}}>t_{\mathrm{call}}$ would not be tallied.

4. Numerical Results for the Male Random Flight Model

We now use a two-dimensional Brownian motion-based particle model with N independent individual male moths to numerically simulate the simplified mating process described so far. We begin with a single female call in Section 4.1 and numerically evaluate the capture and mating probabilities as well as the mean first mating time in Section 4.2. Where possible, we compare our numerical results with analytical estimates. Later, in Section 4.3 and Section 4.4, we will consider multiple calling events unfolding over a typical seven-night period with the female being able to select her mating strategy. Table 1 lists all parameters used.

4.1. Single Female Calling Period

In this section, we consider male moths responding to a single female calling period. Moths are initially placed in an equidistant manner on a circle of radius $r_0={(x_0^2+y_0^2)}^{1/2}>r_{\mathrm{mate}}$, so that $R<1$, and undergo a random walk until the stationary, circular plume of radius $r_{\mathrm{mate}}$ is reached or the total simulation time $t_{\mathrm{sim}}$ elapses. We use an adaptive time stepping scheme to ensure that sufficiently small time steps (hence sufficiently small spatial steps) are taken near the plume circumference to reduce boundary errors in our simulation. If the time step is too large, the corresponding spatial step may take a free moth that starts at the exterior of the plume at time t across a short secant through the plume to an updated position that remains on the exterior of the plume at time $t+\Delta t$. The plume flyover would incorrectly tally this moth as unsuccessful in reaching the female and would result in an underestimate of the mating probability.

The position $\mathbf{x}\left(t\right)$ of a given moth at time t is updated as

(14)$\mathbf{x}(t+\Delta t)=\mathbf{x}\left(t\right)+\sqrt{4D\Delta t}{r}_{\mathrm{s}}\widehat{\mathbf{r}},$

where D is the diffusion constant, $r_{\mathrm{s}}$ is randomly drawn from the standard normal distribution, and $\widehat{\mathbf{r}}$ is a unit vector isotropically oriented. This corresponds to a position-jump process where the mean square displacement after time T is $4DT$, as expected for two-dimensional random movement with diffusion coefficient D. As discussed above, the time step $\Delta t$ is chosen in an adaptive manner depending on the moth’s distance to the circular plume $d\left(t\right)={(x{\left(t\right)}^2+y{\left(t\right)}^2)}^{1/2}-r_{\mathrm{mate}}$. We pose

(15)$\Delta t\left(d\right)=(\Delta t_{\max }-\Delta t_{\min })tanh\left(kd\right)+\Delta t_{\min },$

so that $\Delta t\left(d\right)$ interpolates between $\Delta t_{\min }$ (the shortest time step) and $\Delta t_{\max }$ (the largest time step). In the vicinity of the plume, $\Delta t(d\to 0)$ = $\Delta t_{\min }$, whereas at sufficiently large distances, $\Delta t(d\to \infty )$ = $\Delta t_{\max }$. The parameter k controls the steepness of the transition between $\Delta t_{\min }$ and $\Delta t_{\max }$. We set k = 1 m⁻¹. A male moth is assumed to have captured the plume and to have successfully mated once $d\left(t\right)<0$.

We follow the trajectories of $N=5000$ moths engaged in a random flight for a range of diffusivities D and mating radii $r_{\mathrm{mate}}$. All moths are initiated at an initial radius $r_0=1$ m. Simulations are halted once the simulation time $t_{\mathrm{sim}}=100$ s has elapsed. Since the N male moths are independent, the fraction of moths to have reached the female by $t_{\mathrm{sim}}$ can reasonably approximate the capture and mating probability $C\left(t_{\mathrm{sim}}\right)=1-P\left(t_{\mathrm{sim}}\right)$ as analytically evaluated in Equation (5b).

We find that the analytical expression in Equation (5b) and our simulation results are in good agreement for all surveyed values of D and $r_{\mathrm{mate}}$, as shown in Figure 3A, for the specific choice $t_{\mathrm{sim}}=100$ s. The dimensionless $0\le R\le 1$ is plotted on the horizontal axis of Figure 3A instead of $r_{\mathrm{mate}}$ for compactness. Larger values of $R\to 1$ indicate moths that are initiated closer to the stationary plume, whereas smaller values of $R\to 0$ indicate moths that are initiated further away from the female. A contour plot for the numerical capture and mating probability at time $t_{\mathrm{sim}}$ as a function of D and R is shown in Figure 3B. Increasing the simulation time $t_{\mathrm{sim}}$ yields capture and mating probabilities that approach one.

As can be seen from Figure 3A, and as can be expected, male moths with very low diffusivity are not able to reach the female within $t_{\mathrm{sim}}$ unless they are initially very close to her and $R\simeq 1$. At higher diffusivity, the mating probability is an increasing function of R. Figure 3B reveals the existence of a trade-off between the diffusivity and a moth’s initial distance from the female. For the capture and mating probability $C\left(t\right)$ to remain fixed, decreasing R (or equivalently $r_{\mathrm{mate}}$) requires a concurrent increase in the moth’s diffusivity D. Thus, to be competitive, moths that are less sensitive and that cannot yet sense the plume must be more agile in exploring the physical space.

4.2. First Arrival Time

We now evaluate the expected time for a male moth to first arrive at the plume. The mean first arrival time $\mathbb{E}\left[t_{\mathrm{a}}\right]$ is analytically estimated in Equation (13a) and is a function of male diffusivity D, plume radius $r_{\mathrm{mate}}$, and N. Numerically, it will also depend on the duration of the calling time, which poses an upper limit for the time available to males to reach the female. Here, we we set $t_{\mathrm{call}}$ to be the total simulation time, $t_{\mathrm{call}}=t_{\mathrm{sim}}$. We repeatedly simulate independent male moths performing random flights with diffusivity D and record the time $t_{\mathrm{a}}$ for each to reach the female-generated plume at $r_{\mathrm{mate}}$ starting from an initial radius $r_0$until the number of arrivals is N or the maximum number M of simulations is reached. For those males who are unable to reach the female within $t_{\mathrm{call}}$, no entry is tallied and the attempt is discarded and not used when calculating the mean first arrival time. If $R\simeq 1$ and $t_{\mathrm{call}}$ is large enough, we expect the number of discarded events to be few, and a male reaching the female at times larger than $t_{\mathrm{call}}$ is expected to be a rare event. However, if $R\ll 1$ and $t_{\mathrm{call}}$ is not large enough, many males will fail to reach the female within the allotted time. The mean first arrival time we evaluate here is thus conditioned on the male reaching the female within $t_{\mathrm{call}}$.

For most D and R parameter values, simulating $N=1000$ arrivals is straightforward since the capture/mating probability $C\left(t\right)$ is relatively large; however, for those parameter values with small $C\left(t\right)$, many simulations are required. When $t_{\mathrm{call}}=1800$ s, ${min}_{R,D}C\left(t\right)\approx 0.0175$. Thus, we choose the maximum number of simulations $M=60000$ to obtain approximately $MC\left(t\right)=1050>N=1000$ arrivals. From these $N=1000$ arrival times, we then simply obtain the minimum and average values.

Theoretical results are shown in Figure 4A, where we present a contour plot of the mean first mating time $\mathbb{E}\left[t_{\mathrm{a}}\right]$ as a function of D and R based on Equation (13a). Simulation results for the mean first mating time conditioned on arrival prior to $t_{\mathrm{call}}=t_{\mathrm{sim}}$ are presented in Figure 4B. There is good agreement between the two, indicating that moths reaching the female plume at times larger than $t_{\mathrm{call}}$ are few and do not affect $\mathbb{E}\left[t_{\mathrm{a}}\right]$. Ranges are $\mathbb{E}\left[t_{\mathrm{a}}\right]\simeq {10}^{-5}$ to ${10}^3$, depending on parameter choices. In addition to the first arrival time, we also calculate the mean arrival time which is displayed in Figure 4C. As can be expected for all parameters, this quantity is larger than the expected mean first arrival time, with typical values ranging between $\mathbb{E}\left[t\right]\approx 1$ to ${10}^3$.

Figure 4 reveals that the mean first arrival and mating time $\mathbb{E}\left[t_{\mathrm{a}}\right]$ is lowest for moths with large diffusivity D and large values of $R\simeq 1$, with mating radius $r_{\mathrm{mate}}\simeq r_0$. Similar trends are observed for the mean arrival time $\mathbb{E}\left[t\right]$. One key observation is that the difference between the mean first mating time and the average arrival time can be several orders of magnitude in $\{R,D\}$ phase space.

4.3. Multiple Female Calling Periods

So far, females have been modeled to call males via a single calling period of time $t_{\mathrm{sim}}$. In this section, we model calling over a typical one-week period where calling times can be adjusted nightly according to the chosen constant, increase, decrease, or hat strategies described in Section 2.1 and Figure 1. The female still generates a stationary, circular plume that attracts male moths, and male quality will still be associated with large diffusivities D and large sensitivity radii $r_{\mathrm{mate}}$.

For each of the k nights, with $1\le k\le 7$, we denote the time spent calling by $t_{\mathrm{call}}\left(k\right)$, while the radius of the circular plume per call is fixed at $r_{\mathrm{mate}}$. Regardless of female strategy, the values of $t_{\mathrm{call}}\left(k\right)$ are chosen between a minimum of 0.5 and maximum of 8 h as reported in Section 2.1 from previous experimental observations, such that the total allotted total calling time $T_{\mathrm{call}}={\sum }_{k=1}^7{t}_{\mathrm{call}}\left(k\right)$ does not change among strategies. For $1\le k\le 7$, we thus consider

$\begin{array}{cc}\hfill \mathrm{Increase}:& t_{\mathrm{call}}\left(k\right)=\{0.5000,1.7500,3.0000,4.2500,5.5000,6.7500,8.0000\},\hfill \\ \hfill \mathrm{Decrease}:& t_{\mathrm{call}}\left(k\right)=\{8.0000,6.7500,5.5000,4.2500,3.0000,1.7500,0.5000\},\hfill \\ \hfill \mathrm{Constant}:& t_{\mathrm{call}}\left(k\right)=4.2500,\hfill \\ \hfill \mathrm{Hat}:& t_{\mathrm{call}}\left(k\right)=\{2.5870,3.8804,5.1739,6.4674,5.1739,3.8804,2.5870\}.\hfill \end{array}$

We set $r_0=1$ m and utilize two possible values for $r_{\mathrm{mate}}$: $r_{\mathrm{mate}}=3/4r_0$ for a higher quality male and $r_{\mathrm{mate}}=1/4r_0$ for a lower quality male. Finally, male diffusivities are set at $D={10}^{-3/2}$ or $D={10}^{-7/2}$ m²/s, resulting in four possibilities:

$\{R,D\}=\{0.75,{10}^{-1.5}{\mathrm{m}}^2/\mathrm{s}\},\{0.75,{10}^{-3.5}{\mathrm{m}}^2/\mathrm{s}\},\{0.25,{10}^{-1.5}{\mathrm{m}}^2/\mathrm{s}\},\mathrm{or}\{0.25,{10}^{-3.5}{\mathrm{m}}^2/\mathrm{s}\}.$

These parameter values are chosen as typical values from Figure 3.

To calculate the capture and mating probabilities $C_k$ on each night, we evaluated $C_k\equiv C\left(t_{\mathrm{call}}\left(k\right)\right)$ from Equation (5) with $t=t_{\mathrm{call}}\left(k\right)$ and the four $\{R,D\}$ parameters specified above. We find that the capture and mating probability $C_k$ for male moths with the highest diffusivities and the largest mating radius $r_{\mathrm{mate}}$ remains elevated for all k values and relatively constant across all strategies, as can be seen in Figure 5. Moths with lower diffusivities and where $r_{\mathrm{mate}}$ is smaller display capture and mating probabilities that vary with k and that follow the same trends as calling time strategy used, so that increasing calling time generally increases the mating probability from one night to the next. Intermediate capture and mating probabilities emerge for small diffusivities D but high $r_{\mathrm{mate}}$, and vice versa for large diffusivities D and low $r_{\mathrm{mate}}$, with the latter scenario yielding slightly more favorable outcomes. Figure 5 suggests that the increase or hat strategies may be best for a female to successfully attract a high-quality male through the least calling effort while still finding a mate within the one-week mating period. At onset of the mating process, both the increase or hat strategies yield a large capture and mating probability for high-quality males, with high D and large $r_{\mathrm{mate}}$, and a relatively lower probability for low-quality males, with low D and small $r_{\mathrm{mate}}$. At this stage, the female employs short calling times $t_{\mathrm{call}}\left(k\right)$, ensuring that she is being selective while minimizing her calling effort. As the mating period continues and $t_{\mathrm{call}}\left(k\right)$ increases, the likelihood of attracting lower quality males increases, so that overall mating success with males of any quality increases. This observation is consistent with field observations where it is shown that female moths preferentially employ the increase strategy to optimize their mating likelihood (Umbers et al. [17]).

We confirmed that the increase and hat strategies are the most effective for females. To do so, we defined the effectiveness of the strategy on the kth night, $E_k$, as the ratio of the cumulative total mating probability and total time spent calling (Figure 6). Using the calling times $t_k$ specified above for each strategy, we calculated the mating probability $C_k=C\left(t_k\right)$ in the $\{R,D\}$ parameter plane as in Figure 3B. This defines the total mating probability $M_k$, which is calculated as the sum of $C_k$ over all R and D values. The effectiveness is thus calculated as the ratio of the cumulative total mating probability to the total time spent calling

(16)$\begin{array}{c}\hfill \begin{array}{c}\hfill E_k=\frac{{\sum }_{i=1}^k{M}_i}{{\sum }_{i=1}^k{t}_i}.\end{array}\end{array}$

In Figure 6, we plot $({\sum }_{i=1}^k{t}_i,E_k)$ for each night $k=1,\dots ,7$ and for each strategy. The results show that the increase strategy is the most effective, since a female employing this strategy has high total mating probability obtained with a relatively short calling time. On the other hand, the decrease strategy is the least effective, since a female has a higher total mating probability than in the increase strategy at the cost of a long calling time.

4.4. Night-by-Night Mating Probabilities

For a more nuanced understanding of the mating process, we study the night-by-night mating probability for each of the four strategies shown in Figure 5. We use the same calling time sequences $t_{\mathrm{call}}\left(k\right)$ as in Section 4.3 but consider a wider range of values for the male characteristics R and D. In particular, we evaluate the mating probability ${\mathcal{M}}_{\ell }$ that a female mates on night ℓ given that no male reached her on prior nights, defined as follows:

(17)$\begin{array}{cc}\hfill {\mathcal{M}}_1& =C_1,\hfill \end{array}$

(18)$\begin{array}{cc}\hfill {\mathcal{M}}_{\ell }& =C_{\ell }\prod _{k=1}^{\ell -1}(1-C_k),\ell =2,\dots 7.\hfill \end{array}$

We also calculate the total mating probability $\mathcal{M}$ over the one-week mating period given by

(19)$\begin{array}{c}\hfill \mathcal{M}=C_1+C_2(1-C_1)+\cdots +C_7\prod _{k=1}^6(1-C_k)=\sum _{j=1}^7{C}_j\prod _{k=1}^{j-1}(1-C_k).\end{array}$

The results are shown in Figure 7, Figure 8, Figure 9 and Figure 10. We find that all four calling strategies result in qualitatively similar trends for ${\mathcal{M}}_{\ell }$. Specifically, the mating probability for the first night, ${\mathcal{M}}_1$, is an increasing function of D and R and is largest for high-quality moths, as can be expected. However, the maximum probability on successive nights, ${\mathcal{M}}_{\ell }$ with $\ell >1$, declines in magnitude and shifts to lower diffusivity D and smaller relative radius R. Compared to high-quality moths, those with medium $\{R,D\}$ do not display a large mating likelihood on the first night; however, their mating likelihood remains relatively large on successive nights. This implies that, if only intermediate quality males are present, the female can increase her mating likelihood and partner with one of them if the call extends further than one night. The last panels of Figure 7, Figure 8, Figure 9 and Figure 10 show the total mating probability, $\mathcal{M}$, at the end of the one-week mating process, confirming that mating occurs with large probability for almost all moth qualities considered. This suggests that a female moth can employ a conservative strategy such as the increase strategy to maximize the quality of a mate early in the calling process without risking not copulating at all after the one-week mating period.

5. Plume Navigation Model

We now present an agent-based model of moth mating, called the Plume Navigation Model, where a stationary female emits pheromones and males navigate to find her. As discussed in Section 2.2, plume dynamics and moth flight are modeled in two dimensions since males typically fly at constant altitudes, whether engaged in random flight, surging, or casting behavior. Rather than considering a stationary, circular plume, as done so far, the static female is now assumed to emit a series of pheromone puffs, defining a plume that evolves according to known atmospheric dispersion models, while males search for the plume and fly through it (Cardé [69], Arbas et al. [70], Baker et al. [71], Willis et al. [72], Witzgall [73], Kaissling [74]).

Our Plume Navigation Model is based on previous work by Liberzon et al. [28] and Benelli et al. [33], who developed MothPy [75], an open-source package to simulate moth mating. While the focus of their work was to study male navigational strategies, we developed our Plume Navigation Model to investigate female calling strategies, in particular, by varying the amount of emitted pheromone and its distribution over a typical mating period of seven nights. A schematic of the Plume Navigation Model with typical male trajectories and female pheromone plume dispersion is given in Figure 11, and Table 1 lists all parameters used. Details of the Plume Navigation Model setup are given in Section 5.1 and Section 5.2 and the results are discussed in Section 5.3.

5.1. Male Navigation Algorithm

As originally formulated, MothPy [75] focuses on the latter stage of the male search for a female. Male moths are initially stationary, positioned downwind of the female, and do not actively search for the plume. Rather, the plume will diffuse and be advected toward them. Once a sufficient amount of pheromone reaches the male moths, they start moving upwind, toward the female. Thus, the random flight mode described in Section 2.2 is not employed in MothPy. Conversely, our Plume Navigation Model considers a full night of calling where the male random flight mode is included and where several female strategies are implemented.

We define our coordinate system $\mathbf{x}=(x,y)$ with the female moth at the origin and several, noninteracting, male moths initially positioned downwind of her as shown in Figure 11A. In addition to male moth movement, we track the concentration of pheromone $c(\mathbf{x},t)$ released by the female moth as it diffuses and is convected by wind. We assume that the wind velocity field $\mathbf{U}$ is uniform along the positive x direction so that $\mathbf{U}=\left(U\right(\mathbf{x},t),V(\mathbf{x},t\left)\right)=(U,0)$, where U is a constant. Male moths are assumed to move at a constant flight speed v so that the motion of a moth at location $\mathbf{x}=\left(x\right(t),y(t\left)\right)$ is updated at every time step $\Delta t$ via

(20a)$\begin{array}{cc}\hfill x(t+\Delta t)& =x\left(t\right)+v\Delta tcos\theta ,\hfill \end{array}$

(20b)$\begin{array}{cc}\hfill y(t+\Delta t)& =y\left(t\right)+v\Delta tsin\theta ,\hfill \end{array}$

where $\theta \in [0,2\pi )$ is the angular direction of moth movement, as measured from the x axis. Upwind movement thus occurs along $\theta =\pi$, while casting occurs along $\theta =\pm \pi /2$. We also assume that the pheromone plume can be sensed by a male moth only if it is locally present at a large enough threshold concentration $C_{\mathrm{tol}}$. The algorithm for updating $\theta$ for a moth at position $\mathbf{x}$ is as follows:

• Random flight: If the male moth does not detect pheromone ($c(\mathbf{x},t)<C_{\mathrm{tol}}$) and has never detected pheromone at earlier times $t^{\prime }<t$, the moth chooses a random direction $\theta$ drawn uniformly from $[0,2\pi )$. Note that this motion implies that, overall, the moth executes a random flight with effective diffusivity $D\simeq v^2\Delta t/4$. For the experimentally measured male flight speeds reported in Section 2.2, $v=0.5$ to 5 m/s, and for $\Delta t=1$ s, this corresponds to $D\simeq 1$ to ${10}^{-2}$ m²/s.

• Surging: Upon detection of pheromone signals for $c(\mathbf{x},t)\ge C_{\mathrm{tol}}$, the moth will align with the upwind direction of airflow ($\theta =\pi$) with a margin of error $ϵ$. The updated direction of moment $\theta$ will be selected uniformly from $(\pi -ϵ,\pi +ϵ)$ at each timestep.

• Casting: Upon loss of contact with the pheromone plume for $c(\mathbf{x},t)<C_{\mathrm{tol}}$, the moth will search for the plume perpendicular to the direction of airflow ($\theta =\pm \pi /2$) with a margin of error $ϵ$. The updated direction of movement $\theta$ will be selected uniformly in $(-\pi /2-ϵ,-\pi /2+ϵ)$ or in $(\pi /2-ϵ,\pi /2+ϵ)$ at each timestep. If the male moth is unable to find the lost pheromone plume after casting for a given time $t_{\mathrm{cast}}$, it returns to random flight mode. Otherwise, it returns to surging. Typical total casting times are on the order of 10 s (Martinez et al. [20]), so we choose $t_{\mathrm{cast}}=10$ s.

• Mating: Once the male moth is within a radius $r_{\mathrm{cap}}$ of a female, the two successfully mate.

Example male moth trajectories exhibiting each movement mode (random flight, surging, casting, and mating) are shown as colored trajectories in Figure 11B.

5.2. Female-Generated Pheromone Plume Dynamics

The dynamics of the female-generated pheromone plume dispersing in turbulent wind is evaluated in MothPy [75] via the open source package PomPy (Graham [76]), based on the model of Farrell et al. [46]. We implement a simplified version of its dynamics in our Plume Navigation Model. Here, a stationary female moth, located at the origin, emits a series of i pheromone puffs at time $t_i$, where $i=1,\dots ,n$. Each puff contains $m_i$ pheromone, contributes to the overall plume as a spreading Gaussian centered at ${\mathbf{x}}_p\left(t\right)=(x_p\left(t\right),y_p\left(t\right))$, and is characterized by a time-dependent diffusivity ${\sigma }_i\left(t\right)={\sigma }_0(t-t_i)$ relative to the turbulent diffusivity. The ${\mathbf{x}}_p\left(t\right)$ center of each Gaussian may change in time, as advected by the wind velocity field $\mathbf{U}=(U,0)$. The overall plume concentration $c(\mathbf{x},t)$ at position $\mathbf{x}=(x,y)$ and at time t may thus be written as a superposition of diffusing terms originating from each of the i puffs as follows:

(21)$c(\mathbf{x},t)=\sum _{i=1}^{n}H(t-t_i)\frac{m_i}{2\pi {\sigma }_{\begin{array}{c}\hfill i\end{array}}{\left(t\right)}^2}exp\left(-\frac{|\mathbf{x}-{\mathbf{x}}_p(t-t_i){|}^2}{2{\sigma }_{\begin{array}{c}\hfill i\end{array}}{\left(t\right)}^2}\right),$

where H is the heaviside function and ${\mathbf{x}}_p\left(t\right)$ is advected by the wind velocity field so that, at each time step,

(22)${\mathbf{x}}_p(t+\Delta t)={\mathbf{x}}_p\left(t\right)+(U,0)\Delta t.$

The various female calling strategies can be incorporated by selecting the proper $t_i$ sequences as well as the amount of pheromone released per puff $m_i$.

5.3. The Interplay of Calling Time and Pheromone Amount on Male Fitness

In this section, we consider the full mating problem with a stationary female generating the dynamically evolving plume $c(\mathbf{x},t)$ presented in Equation (21). The N independent males seeking the female follow the plume according to Equation (20) and are characterized by a detection threshold $C_{\mathrm{tol}}$ and a speed v; low $C_{\mathrm{tol}}$ and large v represent high-quality males. We study the mating process for a wide range of calling times and released pheromone amounts.

In our simulations, the female is located at the origin and the $N=10000$ non-interacting male moths are initialized 500 m downwind of her within a 100-m line perpendicular to the direction of the wind flow as shown in Figure 11B. We assume mating occurs once the male is within $r_{\mathrm{cap}}=5$ m from the female. Once this distance is reached, the moth arrival time is recorded and the moth is removed from the simulation. Given the relatively large value of $r_{\mathrm{cap}}$, we do not adaptively change the time step $\Delta t$ at the boundary of the female-generated plume (as in the simplified model) but keep it fixed. To represent male phenotypic variability, each male is assigned a $C_{\mathrm{tol}}$ value uniformly sampled in logarithmic space between ${10}^{-10}$ and ${10}^{-6}$pg/m² and a flight speed v uniformly sampled between 0.5 and 2 m/s.

To simulate different female calling strategies, calling times $t_{\mathrm{call}}$ were varied from 1 to 8 h. Within each strategy, the emitted pheromone per puff was kept fixed at $m_i=m_{\mathrm{p}}$ and the puffs were emitted throughout the simulation at the rate of 1 puff per second. The $m_{\mathrm{p}}$ range we used across strategies was varied between 1 and 10 pg, corresponding to overall emitted pheromone of roughly 2 ng (1 pg emitted over a 30 min call) and 300 ng (10 pg emitted over an 8 h call), exceeding the range of typical values presented in Section 2.1. Each simulation represents one night of calling. Note that, since male moths are independent, concurrently simulating the arrival of N moths is equivalent to N simulations to study the arrival of a single moth. Finally, whereas moth quality in Section 3 was represented by $\{R,D\}$, here, moth-quality indicators are $\{C_{\mathrm{tol}},v\}$.

In Figure 12A, we show how the mating probability (the fraction of the N males that reach the female for all $\{C_{\mathrm{tol}},v\}$ values), depends on the calling time $t_{\mathrm{call}}$ and emitted pheromone $m_{\mathrm{p}}$. Simulations are run over the maximum calling time $t_{\mathrm{sim}}=8$ h when the female call ends. We then tally the total number of males that arrived from the onset of the simulation up to time $t_{\mathrm{call}}$ where $t_{\mathrm{call}}=1,\cdots ,8$ h and compute all relevant statistics associated to $t_{\mathrm{call}}$, such as mating probability, mean speed, and mean chemosensitivity threshold. The results shown in Figure 12A thus represent cumulative subsamples of the maximum $t_{\mathrm{call}}=8$ h call. As can be seen, the mating likelihood is much more sensitive to $t_{\mathrm{call}}$ than to $m_{\mathrm{p}}$, as it increases only by approximately 10% over a tenfold increase in $m_{\mathrm{p}}$. This is in agreement with experimental observations where doubling the effective titer on a moth trap does not substantially affect the number of males that are captured (Johnson et al. [37]).

Figure 12B shows the mean speed of successful males, regardless of $C_{\mathrm{tol}}$. This panel reveals that, by tuning $t_{\mathrm{call}}$, females may preferentially select for males with different speeds: calling for a relatively short time each night (0.5 h) selects for faster males on average, while calling for full night (8 h) gives slower males the opportunity to mate (Figure 12B). Similarly, Figure 12C displays the mean detection threshold of successful males, regardless of v. This panel shows that modulating pheromone titer $m_{\mathrm{p}}$ is an avenue for females to select for males with lower $C_{\mathrm{tol}}$ that are able to sense smaller pheromone concentrations. If we assume that lower $C_{\mathrm{tol}}$ is a proxy for longer antennae in moths, our findings are consistent with experimental observations where doubling the effective titer on a moth trap leads to capture of males with smaller antennae (Johnson et al. [37]).

In Figure 13A we show the mean first arrival times in $\{C_{\mathrm{tol}},v\}$ phase space for a fixed pheromone titer $m_{\mathrm{p}}=1$ pg and calling time $t_{\mathrm{call}}=8$ h. As can be seen, the fastest and most chemosensitive males display the shortest mean first arrival times; vice versa, the slowest and least sensitive males display the longest mean first arrival times. Similar trends are seen in Figure 13B, where instead, we plot the mean arrival and mating time. These findings are in qualitative agreement with the simplified male random flight model where moths with the largest diffusivity D and largest ratio $R=r_{\mathrm{mate}}/r_0$ are the first to reach the female and display the lowest mean arrival times.

Finally, to compare our Plume Navigation Model to existing models, we used MothPy [75] to quantify the fraction of successful males and minimum arrival times as a function of female calling strategies (pheromone puff release rate and molecular amount) and male quality parameters (speed and pheromone detection threshold). The main trends regarding the relationship between female calling and male quality with reproductive success are similar to those obtained via our Plume Navigation Model as discussed in the previous three paragraphs; details of the MothPy simulations can be found in Appendix C.

6. Discussion

In this work, we studied different sexual selection methods used by female moths to attract males with high-quality traits: faster flight speed/diffusivity and better ability to sense small concentrations of pheromone, i.e., chemosensitivity. Specifically, we considered four experimentally observed calling strategies where the female varies calling times and amounts of emitted pheromone over a one-week period in an attempt to optimize both mating success and quality of the first male to mate. Since male moth navigation and pheromone dynamics typically occur at fixed altitude, we only considered two-dimensional representations.

We developed an initial simplified agent-based model (the Male Random Flight Model) where the female-generated pheromone plume is assumed to be static and circular and male moths are assumed to be random walkers before reaching the plume, which allowed us to find analytical estimates for the mean first arrival time. We considered a cohort of N independent males, with diffusion constant D, who can sense the female-generated plume at a distance $r_{\mathrm{mate}}$ that started at a distance $r_0>r_{\mathrm{mate}}$ from the female. For a single calling period of fixed duration, the mean first arrival is found to scale as $\mathbb{E}\left[t_{\mathrm{a}}\right]\simeq {(r_0-r_{\mathrm{mate}})}^2/\begin{array}{c}\hfill (\end{array}DlnN\begin{array}{c}\hfill )\end{array}$. We identified $r_{\mathrm{mate}}$ and D with male quality: large values of $r_{\mathrm{mate}}$ indicate moths that sense the female at larger distances, and large values of D indicate moths that can more readily explore the environment. An interesting avenue for future research would be to extend these results to a variety of random walk models such as those analyzed in Petrovskii et al. [77] and in Bläßle and Tyson [78]. We also considered multiple calling periods distributed over seven nights according to an increase, decrease, hat, or constant pattern such that the total possible calling time, which is a measure of effort on the female’s part, remains constant across strategies. We find that the optimal calling strategy that minimizes expenditure while maximizing male quality and overall mating likelihood is the increase strategy, as has also been shown by experimental studies.

Our more complex agent-based model (the Plume Navigation Model) allowed us to include realistic female-generated plumes that evolve dynamically due to advection by wind and atmospheric dispersion over timescales of several hours. Male navigation is modeled in a more detailed manner: prior to sensing the plume, males are assumed to perform Brownian motion; once the pheromone is detected, they surge upwind toward the female. If the plume is lost, the moth will engage in casting behavior to relocate the plume. At the end of the casting phase, if the pheromone remains undetected, the moth will return to performing Brownian motion, else it will surge toward the female. We used the pheromone detection threshold of male moths, $C_{\mathrm{tol}}$, and their diffusive velocity v as proxies for male quality. Finally, we allowed male moths to search for the female plume over realistic periods of up to eight hours per night.

Where possible, we used experimental studies to inform parameter values. Ideally, a female reproduces at the beginning of adulthood with a high-quality male; however this may not always be possible and the female may need to adjust her calling strategy to attract males of lesser quality before the mating period ends. Our models confirm that, if a female initially releases limited quantities of pheromone $m_{\mathrm{p}}$ per puff and the calling time $t_{\mathrm{call}}$ is brief, only the more motile and chemosensitive males may reach her. Mating likelihoods for moths of lesser quality improve once the amount of released pheromone or calling time increases. We found that the mean first capture and mating times and the mean arrival times can differ by several orders of magnitude, indicating that pheromone-based signaling may help balance the goal of selecting a high-quality male versus the need to mate within a finite time.

We find that overall mating likelihood is more sensitive to $t_{\mathrm{call}}$ than to $m_{\mathrm{p}}$ but that modulating $m_{\mathrm{p}}$ allows the female to select for higher quality males with lower pheromone detection thresholds $C_{\mathrm{tol}}$. If we associate more chemosensitive males with lower $C_{\mathrm{tol}}$ levels to males with longer antennae, these findings are consistent with experimental observations where doubling the effective titer on a moth trap does not change the overall abundance of captured moths, but the antennae of males captured by single-female traps were observed to be longer than the antennae of males captured by double-female traps (Johnson et al. [37]). Simulations using MothPy performed over a 20-s timeframe when males are in close proximity to the female also confirm these findings, as shown in Appendix C.

Antenna structure is an evolutionary aspect of moth anatomy. While some species display large, branching antennae to maximize surface area, sensilla count, and pheromone detection (Symonds et al. [64], Wang et al. [16]), phylogeny shows that a majority of moths have simple non-branching antennae to minimize aerodynamic drag during flight (Symonds et al. [64]). In this case, sensilla and scales are carefully arranged so that chemical detection, pheromone retention, and flight capabilities are concurrently optimized (Wang et al. [16]). Greater power and strength allow moths to sustain large antennae and elevated speeds, so that moth body and antenna sizes are often correlated (Johnson et al. [37]). Our results indicate that maintaining the same mating probability when pheromone detection thresholds decrease (so that the male is more chemosensitive) requires the male flight speed to decrease.

Our models include several simplifying assumptions about moth behavior and interactions. For example, we consider the behavior of many independent males and a single calling female. These models could be expanded to study multiple females simultaneously emitting pheromones and their effects on male selection. Would female–female competition alter the optimal female strategies? There is also evidence that, in some species of moths, female calling attracts conspecific females as well as males (Su et al. [79]). This may lead to spatial clusters of females, which in turn produce larger aggregate plumes capable of attracting more males and males with higher $C_{\mathrm{tol}}$. This is another aspect of female–female behavioral interaction that may alter optimal female strategies. While our models assume that females mate only once and that the first male to reach the female mates with her, additional male–male competitive behaviors (e.g., multiple mating and spermatophore displacement) could be incorporated (Xu and Wang [34]). For species that engage in multiple mating, these behaviors may substantially change the selection pressures and female strategies as the last male to mate, rather than the first, would be considered the successful male. Additionally, while our models were built with a general moth in mind, species-specific parameter values could be used in cases where sufficient experimental data is available.

When studying the Male Random Flight Model, we made the simplifying assumption for mathematical convenience that the plume is a static circle of radius $r_{\mathrm{mate}}$. We thus focus on the “foraging” part of the mating process where the male searches for this static circular plume. This simplifies the mathematical analysis carried out in Section 3 but does ignores the physical fact that the plume is dynamic and is affected by the environment. An extension to our work could be to consider a dynamic trap radius that increases over time $r_{\mathrm{mate}}\left(t\right)$ that coarsely captures the plume’s spread. Although our results are focused on the random movement of individual moths and female calling strategies, it may also be insightful to compare the mean results of many stochastic simulations with the solution to a continuum model describing the chemotaxis-diffusion of male moths and reaction-convection-diffusion of the female pheromone plume. The construction of such a continuum model would require detailed consideration of the distinct modes of male movement (random, surging, casting, and mating) and whether the males undergo position or velocity jumps (Othmer et al. [80]).

There are many aspects of moth life dynamics that may also play a role in the mating process. We studied the mating process over seven nights of calling—a typical moth lifespan. However, other life dynamics such as death rates may be important considerations in future work. Ecologically, moths play a fundamental role near the bottom of the food chain, acting as prey for bats, birds, and other small mammals. Tcheslavskaia et al. [81] estimated that the mortality of female gypsy moths due to predation ranges between 0 to 35% per night, with an average of 14.2 ± 2.5% fatalities per night.

Another possibility is to study hybrid mating. Here, multiple moth strains could be introduced, each responding to a specific pheromone blend, for which the chemical components are expressed in tightly controlled ratios that define blend-specific diffusion and decay rates. These blends allow only moths of the same strain to recognize one another. Although males will initially follow the blend from their own strain, when plumes overlap they may eventually relax their specificity and begin following plumes from other strains, leading to hybrid mating (Kárpáti et al. [82]).

While moths are efficient pollinators (Walton et al. [83]), some species are considered major agricultural pests, and environmentally safe pest control procedures are wanting. For example, mating disruption is a method of releasing formulated synthetic copies of female pheromone that drive males to traps and curb reproduction (Cardé and Minks [84], Barclay and Judd [85], Gordon et al. [86], Lucchi et al. [87]). There is an opportunity for computational simulations, such as those presented in this work, to provide additional information that experiments cannot, as comprehensive field studies are difficult and existing data is sparse. To improve mating disruption methods, minimizing the amount of synthetic pheromone released from automatic devices can reduce costs, both economic- and health-related. Our computational Plume Navigation Model provides information on the effects of female calling strategies which can help inform the timing and length of synthetic pheromone release in mating disruption.

Supplementary Materials

The computer code used to generate the results for this article is available through GitHub at https://github.com/tstepien/moth-mating [v1.0.0]. The code is platform independent and written in MATLAB, which requires a license, and Python.

Author Contributions

Conceptualization, M.R.D. and A.E.L.; methodology, T.L.S., C.Z., J.B.H., R.M.C.R., M.R.D., and A.E.L.; software, T.L.S., C.Z., J.B.H., and A.E.L.; validation, T.L.S., C.Z., J.B.H., M.R.D., and A.E.L.; formal analysis, T.L.S., C.Z., J.B.H., M.R.D., and A.E.L.; writing—original draft preparation, T.L.S., C.Z., J.B.H., R.M.C.R., M.R.D., and A.E.L.; writing—review and editing, T.L.S., C.Z., J.B.H., R.M.C.R., M.R.D., and A.E.L.; visualization, T.L.S., C.Z., J.B.H., and A.E.L.; supervision, M.R.D. and A.E.L.; project administration, T.L.S., M.R.D., and A.E.L.; funding acquisition, M.R.D. and A.E.L. All authors have read and agreed to the published version of the manuscript.

Funding

This material is based upon work supported by the National Science Foundation under grant numbers DMS 1641020 (AMS MRC program), DMS 1562078 (C.Z.), DMS 1815216 (A.E.L.), and DMS 1814090 (M.R.D.); by the Army Research Office under grant numbers W911NF-18-1-0345 and W911NF-16-1-0165 (M.R.D.); and by a Natural Sciences and Engineering Research Council of Canada (NSERC) Postdoctoral Fellowship Award (C.Z.). This work was conducted in part using the resources of the Advanced Computing Center for Research and Education (ACCRE) at Vanderbilt University, Nashville, TN.

Acknowledgments

We wish to acknowledge the American Mathematical Society’s Mathematical Research Community (MRC) in Agent-Based Modeling in Biological and Social Systems, which brought our collaboration together and supported our work.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Derivation of Arrival Time to Circular Plumes

In this appendix, we evaluate the arrival time distribution for a two-dimensional random walker to reach a trapping disk of radius $r_{\mathrm{mate}}$, i.e., $\Omega =\{(x,y)\mid x^2+y^2\le r_{\mathrm{mate}}^2\}$. In Cartesian coordinates, the probability density function $g(x,y,t)$ for a walker to occupy position $\mathbf{x}=(x,y)$ solves the parabolic exterior problem:

(A1a)

$\begin{array}{ll}\frac{\partial g}{\partial t}=D{\nabla }^{2}g,& x^2+y^2>r_{\mathrm{mate}}^2,t>0,\end{array}$

(A1b)

$\begin{array}{ll}g=0,& \hfill x^2+y^2=r_{\mathrm{mate}}^2,t>0,\end{array}$

(A1c)

$\begin{array}{ll}g=\delta (x-x_0)\delta \left(y\right),& \hfill x^2+y^2>r_{\mathrm{mate}}^2,t=0.\end{array}$

Without loss of generality, in Equation (A1b) we set the initial location of the moth to be at $(x_0,0)$, so that its distance from the female is $r_0=x_0$. To solve Equation (A1), it is convenient to first transform to dimensionless variables and a polar coordinate system. Accordingly, we write

$x=r_{\mathrm{mate}}rcos\theta ,y=r_{\mathrm{mate}}rsin\theta ,t=\frac{r_{\mathrm{mate}}^2}{D}\tau ,R=\frac{r_{\mathrm{mate}}}{r_0},g(x,y,t)=p(r,\theta ,\tau ),$

so that the problem for $p(r,\theta ,\tau )$ satisfies

(A2a)

$\begin{array}{cc}{p}_{\tau }=p_{rr}+\frac{1}{r}{p}_r+\frac{1}{r^2}{p}_{\theta \theta },\hfill & \hfill r>1,\theta \in [0,2\pi ),\tau >0,\end{array}$

(A2b)

$\begin{array}{cc}p=0,\hfill & \hfill r=1,\theta \in [0,2\pi ),\tau >0,\end{array}$

(A2c)

$\begin{array}{cc}\hfill p=R\delta (r-R^{-1})\delta \left(\theta \right),& \hfill r>1,\theta \in [0,2\pi ),\tau =0.\end{array}$

The solution to Equation (A2) is now determined with the goal of developing closed form expressions for the following three quantities:

$\begin{array}{cc}\hfill \mathrm{Probability} \mathrm{moth} \mathrm{still} \mathrm{free} \mathrm{at} \mathrm{time}\tau :& {\displaystyle P\left(\tau \right)={\int }_{\theta =0}^{2\pi }{\int }_{r=1}^{\infty }p(r,\theta ,\tau )rdrd\theta ,}\hfill \\ \hfill \mathrm{Probability} \mathrm{moth} \mathrm{trapped} \mathrm{by} \mathrm{time}\tau :& C\left(\tau \right)=1-P\left(\tau \right),\hfill \\ \hfill \mathrm{Arrival} \mathrm{time} \mathrm{distribution} \mathrm{of} \mathrm{moths} \mathrm{to} \mathrm{trap}:& S\left(\tau \right)=-P^{\prime }\left(\tau \right).\hfill \end{array}$

The solution begins by defining the Laplace transform:

(A3)

$\widehat{p}(r,\theta ;s)={\int }_{\tau =0}^{\infty }p(r,\theta ,\tau )e^{-s\tau }d\tau ,s\in \mathbb{C}.$

The quantity $\widehat{p}(r,\theta ;s)$ satisfies the associated elliptic equation:

(A4a)

$\begin{array}{c}{\widehat{p}}_{rr}+\frac{1}{r}{\widehat{p}}_r+\frac{1}{r^2}{\widehat{p}}_{\theta \theta }-c^2\widehat{p}=-R\delta (r-R^{-1})\delta \left(\theta \right),r>1,\theta \in [0,2\pi ),\end{array}$

(A4b)

$\begin{array}{c}\widehat{p}(1,\theta )=0,\theta \in [0,2\pi ),\end{array}$

where $c^2=s$. In practice we are only interested in the radial component of the solution which describes the arrival time distribution. For completeness, we describe the full solution obtained by seeking a separable solution of Equation (A4) that is continuous, periodic, even in $\theta$, bounded as $r\to \infty$, and satisfies the boundary condition $\widehat{p}=0$ on $r=1$. This series is

(A5)

$\widehat{p}(r,\theta )=\left\{\begin{array}{c}{\displaystyle \sum _{n=0}^{\infty }{A}_n\left[I_n\left(cr\right)-\frac{I_n\left(c\right)}{K_n\left(c\right)}{K}_n\left(cr\right)\right]cos\left(n\theta \right),r<R^{-1},}\hfill \\ {\displaystyle \sum _{n=0}^{\infty }{A}_n\left[\frac{I_n\left(c{R}^{-1}\right)}{K_n\left(c{R}^{-1}\right)}-\frac{I_n\left(c\right)}{K_n\left(c\right)}\right]K_n\left(cr\right)cos\left(n\theta \right),r\ge R^{-1},}\hfill \end{array}\right.$

where $A_n$ are constants to be found from incorporation of the Dirac source term. The functions $I_n\left(z\right)$ and $K_n\left(z\right)$ are modified Bessel functions of the first and second kinds, respectively. Multiplying Equation (A4) by $cos\left(m\theta \right)$ and integrating over $(r,\theta )\in (R^{-1}-{\displaystyle \epsilon },R^{-1}+{\displaystyle \epsilon })\times [0,2\pi )$, we find in the limit as ${\displaystyle \epsilon }\to 0$ that

(A6)

$\underset{{\displaystyle \epsilon }\to 0}{lim}{\int }_{r=R^{-1}-{\displaystyle \epsilon }}^{R^{-1}+{\displaystyle \epsilon }}{\int }_{\theta =0}^{2\pi }cos\left(m\theta \right)\frac{\partial }{\partial r}\left(r\frac{\partial \widehat{p}}{\partial r}\right)drd\theta =-1.$

The orthogonality of ${\{cosm\theta \}}_{m=0}^{\infty }$ over $[0,2\pi )$ yields the jump condition

(A7)

$\begin{array}{c}\frac{A_{n}c}{R}\left[{\int }_{\theta =0}^{2\pi }{cos}^2\left(n\theta \right)d\theta \right]\left[\left(\frac{I_n\left(c{R}^{-1}\right)}{K_n\left(c{R}^{-1}\right)}-\frac{I_n\left(c\right)}{K_n\left(c\right)}\right)K_n^{\prime }\left(c{R}^{-1}\right)\right.\hfill \\ \hfill \left.-I_n^{\prime }\left(c{R}^{-1}\right)+\frac{I_n\left(c\right)}{K_n\left(c\right)}{K}_n^{\prime }\left(c{R}^{-1}\right)\right]=-1.\end{array}$

The jump condition (A7) can be simplified with the Wronskian identity

(A8)

$K_n\left(z\right)I_n^{\prime }\left(z\right)-I_n\left(z\right)K_n^{\prime }\left(z\right)=\frac{1}{z},$

to find that

(A9)

$A_n=K_n\left(c{R}^{-1}\right){\left[{\int }_{\theta =0}^{2\pi }{cos}^2\left(n\theta \right)d\theta \right]}^{-1}=\left\{\begin{array}{cc}{\displaystyle \frac{K_0\left(c{R}^{-1}\right)}{2\pi },}\hfill & n=0,\hfill \\ {\displaystyle \frac{K_n\left(c{R}^{-1}\right)}{\pi },}\hfill & n\ge 1.\hfill \end{array}\right.$

The Laplace transform of the flux over the disk is then given by

(A10)

$\frac{\partial \widehat{p}}{\partial r}{|}_{r=1}=\sum _{n=0}^{\infty }\frac{A_n}{K_n\left(c\right)}cos\left(n\theta \right)=\frac{1}{2\pi }\frac{K_0\left(c{R}^{-1}\right)}{K_0\left(c\right)}+\frac{1}{\pi }\sum _{n=1}^{\infty }\frac{K_n\left(c{R}^{-1}\right)}{K_n\left(c\right)}cosn\theta .$

Returning to the variable to $\sqrt{s}=c$, we now apply the inverse Laplace transform of Equation (A10) and evaluate the contour integrals:

(A11)

$\frac{1}{2\pi i}{\int }_{\Gamma }\frac{K_n\left(\sqrt{s}{R}^{-1}\right)}{K_n\left(\sqrt{s}\right)}{e}^{s\tau }ds,n=0,1,2,\dots ,$

where the integral is evaluated over the complex curve $\Gamma =\{\gamma +iy|y\in (-\infty ,\infty )\}$. The parameter $\gamma =Re(\Gamma )$ is chosen such that singularities of the integrand lie to the left of $\Gamma$. Since the integrand has no poles and its only singularity is a branch cut on the negative real axis, we deform the contour to a hairpin along $Re\left(s\right)<0$ and introduce the substitution $s=-w^2$. The integral becomes

(A12)

$\begin{array}{c}\frac{1}{\pi i}{\int }_{\omega =0}^{\infty }\left[\frac{K_n\left(i\omega R^{-1}\right)}{K_n\left(i\omega \right)}-\frac{K_n(-i\omega R^{-1})}{K_n(-i\omega )}\right]\omega e^{-{\omega }^2\tau }d\omega \hfill \\ \hfill =\frac{2}{\pi }{\int }_{\omega =0}^{\infty }\left[\frac{J_n\left(\omega \right)Y_n\left(\omega R^{-1}\right)-J_n\left(\omega R^{-1}\right)Y_n\left(\omega \right)}{Y_n^2\left(\omega \right)+J_n^2\left(\omega \right)}\right]\omega e^{-{\omega }^2\tau }d\omega ,\end{array}$

where we have used the identity (NIST DLMF ([88] § 10.27))

(A13)

$K_n(-iz)=\frac{\pi }{2}\left[-Y_n\left(z\right)+i{J}_n\left(z\right)\right],$

to write the integrand in terms of the Bessel functions of the first and second kinds, $J_n\left(z\right)$ and $Y_n\left(z\right)$, respectively.

The expression for the flux $\mathcal{J}(\tau ,\theta )=p_r{|}_{r=1}$ is now

(A14)

$\begin{array}{c}\mathcal{J}(\tau ,\theta )=\frac{1}{{\pi }^2}{\int }_{\omega =0}^{\infty }\left[\frac{J_0\left(\omega \right)Y_0\left(\omega R^{-1}\right)-J_0\left(\omega R^{-1}\right)Y_0\left(\omega \right)}{Y_0^2\left(\omega \right)+J_0^2\left(\omega \right)}\right]\omega e^{-{\omega }^2\tau }d\omega \hfill \\ \hfill +\frac{2}{{\pi }^2}\sum _{n=1}^{\infty }cos\left(n\theta \right){\int }_{\omega =0}^{\infty }\left[\frac{J_n\left(\omega \right)Y_n\left(\omega R^{-1}\right)-J_n\left(\omega R^{-1}\right)Y_n\left(\omega \right)}{Y_n^2\left(\omega \right)+J_n^2\left(\omega \right)}\right]\omega e^{-{\omega }^2\tau }d\omega .\end{array}$

The total flux to the trapping region is given by

(A15)

$S\left(\tau \right)={\int }_{\theta =0}^{2\pi }\mathcal{J}(\tau ,\theta )d\theta =\frac{2}{\pi }{\int }_{\omega =0}^{\infty }\left[\frac{J_0\left(\omega \right)Y_0\left(\omega R^{-1}\right)-J_0\left(\omega R^{-1}\right)Y_0\left(\omega \right)}{Y_0^2\left(\omega \right)+J_0^2\left(\omega \right)}\right]\omega e^{-{\omega }^2\tau }d\omega .$

This gives the distribution of arrival times to the inner disk. For a fixed time $\tau ={\tau }^{*}$, the distribution of arrival angles $\theta$ is then given by $\mathcal{J}({\tau }^{*},\theta )$.

The free probability $P\left(\tau \right)$ (2) that a moth has not yet reached the circular trapping region by time t satisfies $S\left(\tau \right)=-P^{\prime }\left(\tau \right)$ and therefore

(A16)

$P\left(\tau \right)={\int }_{\tau }^{\infty }S\left(\eta \right)d\eta =1-{\int }_0^{\tau }S\left(\eta \right)d\eta .$

The captured fraction is then given by the complement $C\left(\tau \right)=1-P\left(\tau \right)$. The survival probability is given by

(A17)

$\begin{array}{cc}\hfill P\left(\tau \right)& ={\int }_t^{\infty }S\left(\eta \right)d\eta \hfill \\ \hfill & =\frac{2}{\pi }{\int }_{\omega =0}^{\infty }\left[\frac{J_0\left(\omega \right)Y_0\left(\omega R^{-1}\right)-J_0\left(\omega R^{-1}\right)Y_0\left(\omega \right)}{Y_0^2\left(\omega \right)+J_0^2\left(\omega \right)}\right]\frac{1}{\omega }{e}^{-{\omega }^2\tau }d\omega .\hfill \end{array}$

where $S\left(\eta \right)$ is given in (A15). The integrals (A15) and (A17) are evaluated numerically in MATLAB. The results in dimensionless time $\tau$ are rescaled to physical time t through the relationship

$\tau =\frac{D}{r_{\mathrm{mate}}^2}t.$

Appendix B. Asymptotic of the Survival Probability at Short Times

To obtain the limiting distribution of first arrival times, we estimate the short time asymptotic for survival probability $P\left(t\right)$ derived in (A17). The first step is to apply the Laplace transform to the relationship $S\left(t\right)=-P^{\prime }\left(t\right)$, which yields

(A18)

$\widehat{S}\left(s\right)=1-s\widehat{P}\left(s\right).$

By combining (A10) and (A18), we obtain

(A19)

$\widehat{S}\left(s\right)={\int }_0^{2\pi }\frac{\partial \widehat{p}}{\partial r}{|}_{r=1}d\theta =\frac{K_0\left(\sqrt{s}{R}^{-1}\right)}{K_0\left(\sqrt{s}\right)},\widehat{P}\left(s\right)=\frac{1}{s}-\frac{K_0\left(\sqrt{s}{R}^{-1}\right)}{s{K}_0\left(\sqrt{s}\right)}.$

Our goal is to obtain the limiting behavior of $P\left(t\right)$ as $t\to 0$. This corresponds to the limiting behavior of $\widehat{P}\left(s\right)$ as $s\to \infty$. By applying the large argument expression $K_0\left(z\right)\sim \sqrt{\frac{\pi }{2z}}{e}^{-z}$ as $z\to \infty$, we have

(A20a)

$\begin{array}{c}\widehat{P}\left(s\right)\sim \frac{1}{s}-\frac{\sqrt{R}}{s}{e}^{-\frac{(1-R)\sqrt{s}}{R}}\left[1+\mathcal{O}\left(s^{-\frac{3}{2}}\right)\right],s\to \infty ,\end{array}$

(A20b)

$\begin{array}{c}P\left(t\right)=1-\frac{2\sqrt{Dt}}{r_0(1-R)\sqrt{\pi R}}{e}^{-\frac{r_0^2{(1-R)}^2}{4Dt}}\left[1+\mathcal{O}\left(t^{-\frac{3}{2}}\right)\right],\tau \to 0.\end{array}$

Appendix C. MothPy Simulations

We used MothPy [75] to run simulations to compare to the results of our Plume Navigation Model (Section 5). We initialized a single stationary female moth at position $(x,y)=(25,500)$ cm with 200 male moths uniformly distributed within the rectangle $\left\{\right(x,y)\mid x\in [393,450]\mathrm{cm},y\in [490,499\left]\mathrm{cm}\right\}$. The male moth navigational strategy we chose to simulate is called “Strategy A” in Benelli et al. [33] and is the same as that implemented in Liberzon et al. [28]: upon detecting the pheromone plume, male moths surge upwind, and upon losing the plume, they cast. The wind and plume models used are from the software package PomPy (Graham [76]), which implements the models in Farrell et al. [46] and as described in Equations (21) and (22). Contrary to our Plume Navigation Model implementation, in MothPy [75], the pheromone puff release is modeled as a stochastic Poisson process and the wind velocity includes random white noise. Finally, MothPy [75] focuses on very short time scales, of the order of tens of seconds, whereas our Plume Navigation Model unfolds over several hours. We use the default parameters as given in version 0.0.1 of MothPy [75] but vary MothPy’s male moth quality and female calling indicators through a parameter sweep of four variables as follows:

• female pheromone puff release rate:

  $f_{\mathrm{r}}=\{50,66.66,83.33,100,116.66,133.33,150,166.66,183.33,200\}{\mathrm{s}}^{-1}$

• female pheromone puff molecular amount:

  $m_{\mathrm{p}}=\{0.1,0.68,1.26,1.84,2.42,3.00,3.58,4.16,4.74,5.32,6.48,8.80,13.44,22.72,41.28\}\begin{array}{c}\hfill \mathrm{pg}\end{array}$

• male ground flight speed:

  $v=\{50,140,230,320,410,500,590\}\mathrm{cm}/\mathrm{s}$

• male pheromone detection threshold:

  $C_{\mathrm{tol}}=\{120,500,880,1260,1640,2020,2400\}\begin{array}{c}\hfill {\mathrm{pg}/\mathrm{cm}}^2\end{array}$

Each of the 7350 simulations is run for 20 s, and a male moth is considered to have reached the female when located within a radius of $r_{\mathrm{cap}}=15$ cm from her. We show the percentage of successful mating in Figure A1 and the first arrival time of male moths in Figure A2.

Although run over a different timescale than our Plume Navigation Model simulations, these MothPy simulations confirm the trends discussed in the main text. Specifically, increasing female effort, here represented by pheromone release rate $f_{\mathrm{r}}$ and titer amount $m_{\mathrm{p}}$, increases the mating likelihood of males traveling at speed v and characterized by pheromone detection threshold $C_{\mathrm{tol}}$ while at the same time decreasing the mean first arrival time. This is also observed for our long-term Plume Navigation Model simulations. Also note that the probability of successful mating for males with relatively low speed v is is more sensitive to changes in the female’s calling effort ($m_{\mathrm{p}}$ and $f_{\mathrm{r}}$) as can be seen from the broader color range in the bottom rows of panels in Figure A1 (e.g., for $v=140,230,320$ cm/s) as opposed to the top rows (e.g., for $v=410,500,590$ cm/s), another feature observed in our Plume Navigation Model simulations and discussed in the main text.

Figure A1

Successful mating percentage using MothPy [75] as described in Appendix C. Velocities v are measured in cm/s; female pheromone amounts $m_{\mathrm{p}}$ are measured in pg; pheromone puff release rates $f_{\mathrm{r}}$ are measured is s${}^{-1}$; and pheromone detection thresholds $C_{\mathrm{tol}}$ are measured in pg/cm${}^2$.

[Figure omitted. See PDF]

Figure A2

First male arrival time to the female, normalized by the total simulation time $t_{\mathrm{sim}}=20$ s using MothPy [75] as described in Appendix C. White regions correspond to no male moths reaching the female within $t_{\mathrm{sim}}$, velocities v are measured in cm/s; female pheromone amounts $m_{\mathrm{p}}$ are measured in pg; pheromone puff release rates $f_{\mathrm{r}}$ are measured is s${}^{-1}$; and pheromone detection thresholds $C_{\mathrm{tol}}$ are measured in pg/cm${}^2$.

[Figure omitted. See PDF]

Figures and Table

Figure 1. Schematic of four female moth calling strategies, expressed as female calling effort over a typical seven night call. Calling effort is quantified as the time females spent calling males and/or the pheromone titer.

Figure 2. (A) Typical profiles for P(t),C(t)=1−P(t), and S(t) as derived from Equation (5), for D = 1 m/s2, rmate=1 m, and r0=5 m. The arrival time distribution S(t) is plotted on an axis with a smaller scale for clarity (red axis, right). (B) Example numerical simulation for 6 independent male moths (various colors; D = 1 m/s2), with rmate=1 m (solid circle) and r0=5 m (dotted circle).

Figure 3. The probability for a male moth to reach a stationary, circular female-generated plume after tsim=100 s in the Male Random Flight Model. (A) Capture and mating probability C(t) as a function of R=rmate/r0, calculated from Equation (5b) for male diffusivities D=10−1, 10−2, 10−3, and 10−4. Numerical results obtained by simulating N=5000 independent moths as described in Section 4.1 are overlaid. (B) Capture and mating probability C(t) plotted via contours in {R,D} space.

Figure 4. Male arrival time distributions evaluated in {R,D} phase space, where R=rmate/r0. (A) Mean first arrival times E[ta] calculated from Equation (13a) for N=1000 moths. (B) Mean first arrival times E[ta] calculated from numerical simulations with N=1000 arrivals and tcall=1800 s (Section 4.1). (C) Mean arrival times E[t] calculated from numerical simulations. Mean first arrival times E[ta] are the lowest for large diffusivities D and for R→1, as can be expected. Note the nine orders of magnitude variation for E[ta] across the {R,D} ranges considered. The mean arrival time E[t] is generally much larger than E[ta] and only spans three orders of magnitude in the same {R,D} space.

Figure 5. Effectiveness of each of the four female calling strategies described in Section 2.1 and Figure 1. Capture and mating probability C(t) for four representative moths with {R,D}={0.75,10−1.5 m2/s} (blue curve), {0.75,10−3.5 m2/s} (red curve), {0.25,10−1.5 m2/s} (yellow curve), and {0.25,10−3.5 m2/s} (purple curve). Calling times tcall(k) for 1≤k≤7 nights are specified in Section 4.3. Male moths of the highest quality display similar capture and mating probability across female strategies, with C(t)≃1 in all cases. For intermediate quality moths C(t) depends on the chosen female strategy.

Figure 6. Effectiveness (mating probability/time spend calling) of the four female calling strategies. Each dot corresponds to a night. The increase strategy is the most effective for small calling times, while all strategies have similar effectiveness for large calling times.

Figure 7. Night-by-night and total mating probabilities for the increase calling strategy. Note the different likelihood scales for each panel. On the first night, the mating likelihood M1 is highest for high-quality moths (large {R,D}), and lower for intermediate and low-quality moths (medium and low {R,D}). The latter reaches the female with relatively large likelihoods on the second and successive nights, as can be seen for Mℓ with ℓ>1. The last panel displays the total capture and mating likelihood at the end of the one-week mating period.

Figure 8. Night-by-night and total mating probabilities under decrease calling. Note the different likelihood scales for each panel. Similar to Figure 7, intermediate-quality moths (medium {R,D}) display relatively large mating likelihoods Mℓ on the second and successive nights.

Figure 9. Night-by-night and total mating probabilities under constant calling. Note the different likelihood scales for each panel. Similar to Figure 7, intermediate-quality moths (medium {R,D}) display relatively large mating likelihoods Mℓ on the second and successive nights.

Figure 10. Night-by-night and total mating probabilities for the hat calling strategy. Note the different likelihood scales for each panel. Similar to Figure 7, intermediate-quality moths (medium {R,D}) display relatively large mating likelihoods Mℓ on the second and successive nights.

Figure 11. (A) Schematic of the agent-based Plume Navigation Model as described in Section 5. A stationary female at the origin releases pheromone puffs while male moths fly upwind toward her. Note the capture radius rcap and the presence of an advecting wind. (B) Representative female pheromone plume and male flight trajectories in our Plume Navigation Model simulations as described in Section 5.

Figure 12. Simulation results from the agent-based Plume Navigation Model described in Section 5 for N=10,000 moths. The range of the male quality trait Ctol is 10−10 to 10−6pg/m2; for v, it is 0.5 to 2 m/s. (A) Contour plot of mating probability as a function of female calling effort in {tcall,mp} phase space, regardless of male quality: The panel tallies moths with all levels of Ctol and v that reach the female. (B) Contour plot of the mean speed E[v] of males upon reaching the female in {tcall,mp} phase space and for any level of male Ctol. (C) Contour plot of the mean detection threshold E[Ctol] of males upon reaching the female in {tcall,mp} phase space for any level of male v. Male arrival likelihood is highly dependent on tcall: modulating tcall allows the female to select for male speed v, and modulating mp allows the female to select for Ctol.

Figure 13. Simulation results from the agent-based Plume Navigation Model described in Section 5 for N=10,000 moths attempting to reach a female calling for tcall=8 h via the continued release of mp=1 pg per second. (A) Contour plot of mean first arrival time E[ta] of males in {Ctol,v} phase space. (B) Contour plot of the mean arrival time of males in {Ctol,v} phase space. Both mean first arrival times and mean arrival times are measured in hours.

Table 1

Summary of variables used.

Symbol	Representation	Value/Units
$\mathbf{x}$	Position of moth	m
$t_{\mathrm{a}}$	First arrival and mating time	s
	Male Random Flight Model
$\Delta t_{\max }$	Largest simulation time step	${10}^{-2}$ s
$\Delta t_{\min }$	Smallest simulation time step	${10}^{-5}$ s
$t_{\mathrm{sim}}$	Total simulation time	100 s
N	Number of male moths	1000 to 5000
$r_{\mathrm{mate}}$	Mating radius	0 to 1 m
$r_0={(x_0^2+y_0^2)}^{1/2}$	Initial position of male moth	1 m
$r_{\mathrm{mate}}$	Mating radius	0 to 1 m
$R=r_0/r_{\mathrm{mate}}$	Male initial position to mating radius ratio	dimensionless
D	Male moth diffusion constant	${10}^{-1}$ to ${10}^{-4}$ m2/s
	Plume Navigation Model
$\Delta t$	Time step	1 s
N	Number of male moths	10,000
${\mathbf{x}}_p$	Center of wind-advected plume	m
v	Speed of male moths	0.5 to 2 m/s
$C_{\mathrm{tol}}$	Pheromone detection threshold	${10}^{-10}$ to ${10}^{-6}$ pg/m2
$ϵ$	Casting and surging margin of error	$\pi /6$
$t_{\mathrm{cast}}$	Maximum time spent casting	10 s
$m_{\mathrm{p}}$	Pheromone amount per puff	1 to 10 pg
$f_{\mathrm{r}}$	Puff release rate	1 puff/s
${\sigma }_0$	Puff growth parameter	0.1 m/s
U	Wind velocity	5 m/s
$r_{\mathrm{cap}}$	Capture distance	5 m