Mass action in two-sex population models: encounters, mating encounters and the associated numerical correction

Ideal gas models are a paradigm used in Biology for the phenomenological modelling of encounters between individuals of diﬀerent types. These models have been used to approximate encounter rates given densities, velocities and distance within which an encounter certainly occurs. When using mass action in two-sex populations, however, it is necessary to recognize the diﬀerence between encounters and mating encounters . While the former refers in general to the (possibly simultaneous) collisions between particles, the latter represents pair formation that will produce oﬀspring. The classical formulation of the law of mass action does not account this diﬀerence. In this short paper, we present an alternative derivation of the law of mass action that uses dimensional reduction together with simulated data. This straightforward approach allows to correct the expression for the rate of mating encounters between individuals in a two-sexpopulationwithrelativeease.Inaddition,variabilityinmatingencounterrates(duetoenvironmentalstochasticity)isnumericallyexploredthroughrandomﬂuctuationsonthenewmassactionproportionalityconstant.ThesimulationsshowhowtheconditionedtimetoextinctioninapopulationsubjecttoareproductiveAlleeeﬀectisaﬀected.


Introduction
In chemistry, the law of mass action states that the rate of a reaction is proportional to the product of the concentrations of the reactants. Since Alfred Lotka pioneered its use in Lotka (1925) to justify the encounter term in his predator-prey system of differential equations, the law has become ubiquitous in mathematical ecology for modelling the interactions between individuals of different groups. Lotka's arguments for the use of mass action in biological encounters were motivated by the analogy to the kinetic theory of gases (ideal gas models) and have been applied to describe a variety of phenomena, including fertilization kinetics, search theory and mate finding, see Hutchinson and Waser (2007), Voit, Martens, and Omholt (2015) for thorough reviews. Before Lotka, the chemical law of mass action was used for the first time by A.G. McKendrick for describing the interactions among susceptible and infected individuals in epidemiological contexts (Heesterbeek, 2005;McKendrick, 1912).
In the theory of molecular collisions in gases, the collision frequency among particles of different types is expressed in terms of the velocities and radiuses of particles. The quantitative law is deduced directly from geometrical abstractions using the mean of a Poisson process that models the number of collisions a particle receives from others and where simultaneous collisions are allowed (Hutchinson & Waser, 2007;Kauzmann, 2012). At relatively low population densities, the same idea is used as a phenomenological approach to approximate the encounter rates between males and females, with the birth rate taken proportional to the product of their densities (Bazykin, 1998;Fauvergue, 2013). We remark that two-sex population models become relevant when sexual dimorphism in vital rates is present, which has been observed in several species (Caswell, 2001).
While ideal gas models, which assume a heterogeneous population of linearly moving individuals, are used for modelling encounters of individuals, dispersal and movement of homogeneous populations in ecology are modelled following diffusion models, i.e. where individuals move randomly in space (Codling, Plank, & Benhamou, 2008;Turchin, 1998). However, due to the complexity involved, most of the studies that include encounters in models of randomly moving organisms are based on computer simulations (see for instance Bartumeous, Catalan, Viswanathan, Raposo, & da Luz, 2008;Gurarie & Ovaskainen, 2011;Brown, 2008, 2010, where encounters are defined in general animal search problems). This paper addresses the modelling of encounters that lead to reproduction, which requires ruling out simultaneous encounters that occur in the ideal gas model, i.e. encounters where a single female mates two or more males simultaneously, otherwise the rate at which new offspring appear would be overestimated. Thus, we differentiate between counting encounters, possibly simultaneous and counting mating encounters, which are understood here as the formation of female-male pairs from which offspring are successfully produced.
Our aim here is to present a correction for the constant used in the mass action term corresponding to the ideal gas model that accounts for the pair formation. To achieve this goal, we first build a functional relation among the variables using dimensional reduction and simulated data of individuals' movement. This allows us to approximate the value of the proportionality constant for the mass action with relative precision in comparison to theoretical results. Then, we generate new data through computer simulations that only count female-male pairs, and with this we approximate the new value for the constant. Finally, we use the new constant in the mass action law to explore the effects of environmental stochasticity on the conditioned time to extinction for a population model via the variability on encounter rates. The stochastic model is derived from a deterministic model that uses mass action at low population densities and thus shows a reproductive Allee effect (Courchamp, Berec, & Gascoigne, 2008). The simulations reveal the effects that simultaneous random fluctuations around the new constant (and therefore the mating encounter rate) and demographic stochasticity have on the extinction time. This elementary example justifies having reasonable approximations for the non-linear term that models encounters: it demonstrates how variability in the environment could play an essential role in regulating the time to extinction distribution.

Dimensional reduction for mass action
We consider a two-sex population and assume that individuals' movement is done in two spatial dimensions, i.e. we assume that the total variation in the displacement of each individual on the plane is much larger than changes made in their altitude. We also assume that individuals of both sexes are initially homogeneously mixed and uniformly distributed in space, and move over a terrain with area A under the following assumptions: • The velocity, v, is constant and the same for both sexes, • Individuals move independently from each other, • Movement is in straight lines, • The initial direction of movement for each individual is independently chosen at random from the interval [0, 2π), • Individuals' sizes are negligible.
Let n m and n f be the number of males and females, respectively, and c the average number of encounters that one female has with males during the observation time t. We assume that c is related to (i) the velocity v, (ii) the density of males n m /A and (iii) the size of a small area surrounding the female where males are attracted to mate. This is thought as a circular area with radius R. We remark that these assumptions are the same used in the theory of molecular collisions to deduce the law of mass action for gases of two different types, with just the words 'particle' replacing 'individual' and 'type' instead of 'sex' (Hutchinson & Waser, 2007;Kauzmann, 2012). We write the relation between the system parameters in terms of some (unknown) function F, By the -Theorem, see Barenblatt (2006) or Logan (2006) for instance, Equation (1) is equivalent to a relation that involves only the dimensionless quantities with f yet to be determined. Next, we use data generated from agent-based simulations that count the number of contacts (with males) per female. For the simulations, the individuals were programmed to follow the rules stated above and, for a single female, we counted a contact when its distance to a male is less than R. In the computations, the units chosen for length and time were meters and hours (1 h = 1 time step). We used fixed values for the time of observation, t = t * = 24 (h), the radius R = R * = .05 (m) and the number of males, n m * = 100, while varying area size A and individuals' velocity v. On the plane v − A, we arbitrarily choose the strip [25, 4 × 10 4 ] × [25, 2 × 10 3 ], and within this domain, we fixed 8 values of areas: (25,35,50,400,625,900,1200,2000), and 28 values of velocities: (25,50,75,100,200,300,400,500,600,700,800,900,1000,1200,1400,1600,1800,3000,3200,5400,6400,9000,11000,16000,18000,22400,25200,40000). These values were chosen with the aim of capturing the characteristics of the surface at low and high parameter values. Using periodic boundary conditions, we then produced 10 simulations for each corresponding pair of parameters and finally computed the contact averages. The simulations were run in NetLogo 1 Railsback & Grimm (2012) and the code is available upon request to K. Snyder. 2 It is well known that the average number of contacts observed during a fixed period of time should increase with larger velocities and decrease with larger areas, therefore suggesting a relation of the form where p, q and a are constants. This expression can then be rewritten in terms of the dimensionless quantities 1 and 2 , Using least squares to fit relation (3) to the averages of the data points collected gave approximate values of p ≈ 1, q ≈ 1 and K = 1.2887. Denoting with n f the number of females and with C the total averaged number of encounters that females have with males, i.e. C = cn f , gives where x and y are the densities of males and females, respectively. The right-hand term in (4) is traditionally obtained from the theory of molecular collisions, also known as the ideal gas model. That theory provides a constant value of K = 4/π = 1.2732..., which is in good agreement with the value obtained through our simulations (the relative error is less than 2%).

Mating encounters
We generate different data by repeating the simulations of individual movement with the same assumptions as above but counting at most one mating encounter per female at each time step ruling out simultaneous mating. Fitting the model to the new data produces p ≈ 1, q ≈ 1 and K ≈ .1231, which is less than 10% of the previous K value.
Our interest now is to include changes in movement direction at every step in time, depending on the previous direction rather than restricting individual movement to straight lines as in the assumptions. A limitation in the theoretical ideal gas model, as traditionally conceived, is that it does not capture the effects of varying correlation in individuals' movement (i.e. the autocorrelation of the directions of subsequent individual moves) on contact rates. From the dimensional analysis, however, we conclude that the degree of correlation should appear as a functional dependence between the value K and the range of possible directions for individuals' movement, i.e. K = K(θ), where at each time step each individual changes to a direction chosen uniformly at random from [−θ, θ], 0 ≤ θ ≤ π, independently of the other individuals. Repeating the numerical experiments but now allowing individuals to change movement direction with different degrees of correlation, we compute the average of K for different values of θ. Unsurprisingly, these values appear almost constant, as first pointed out by Skellam (1958) for the classical  (3) to the original average data points ± standard deviation, respectively. Bottom: Mating encounter rates (encounters/h) for moving individuals as function of the population density (individuals/m 2 ) and velocity (m/h). For illustration purposes, the detection radius R was fixed and chosen equal to .02 m and the velocities are in the range from 1 to 40 Km/h, which includes estimates for several insect species, see Table 1. The area A is one square meter. mass action, see Figure 1 (top). This information can thus be used to approximate average mating encounter rates for individuals, given estimates of velocities, see Figure 1 (bottom). We emphasize, however, that the assumptions made on the movement constitute a rough simplification of reality: males and females do not necessarily move at the same speed and the details of mating mechanisms have been deliberately left out.
To assess the variability due to movement correlation on the parameter K, we re-fit Equation (3) for different values of θ to two different data-sets defined by taking the simulated data at each corresponding point of the v − A plane and then computing average+standard deviation and average-standard deviation. With these sets of points, we obtained two new constants (that depend on θ) and denoted by K + and K − . The data suggested that for values of θ close to zero, the differences are consistently larger, giving larger values for K + and K − (See Figure 1). But with increasing θ, although the movement is more irregular, large variability in the values of K + and K − is absent. Our interpretation of this outcome is that, although individual movement is apparently more convoluted for larger θ, at each trial the (stochastic) process is the same. On the contrary, for values of θ close to zero, the initial random directions for each individual define paths that look like the deterministic trajectories (θ = 0) in the domain (in this case, a torus). Those might differ completely each time the experiment is repeated because of the randomization of initial directions. This fact was found to be independent whether the boundary conditions are periodic or reflective after running additional simulations for the latter.

Effects of environmental stochasticity
The rate at which offspring are generated by a two-sex population in real scenarios is likely to be subject to random fluctuations due to environmental factors, like rainfall and temperature, but with sensitivity that is species specific. If mass action is chosen to model mating encounters, variability on the parameter K will appear as consequence of those fluctuations. This translates directly into variability of birth and death rates. As an illustration, let us initially consider the simplest deterministic population model ephemeral mating interactions, where x and y are the densities of males and females, respectively, P(x, y) is a (symmetric) birth rate with even sex ratio for all births and μ is the death rate that is assumed equal for both sexes. Suppose also that the initial sex ratio is 1:1. It is natural that we would like to use the law of mass action for P(x, y), i.e. proportional to xy, but as D.G. Kendall first pointed out, this leads to solutions that blow up in finite time (Kendall, 1949). One way, this trouble can be fixed by taking into consideration the average refractory time τ of Figure 2. Cumulative probability of extinction (conditioned to extinction) obtained from the stochastic model defined in Section 3. The effect of variability on the number of contacts, K ∼ Normal(K , σ 2 ), is shown for different values of σ 2 = 0 (continuous line), and .1 (dash-dots). The initial population is taken equal to half the Allee threshold. The simulations suggest that an increasing variance in the random fluctuations corresponding to mating encounters pushes the probability mass to the left, i.e. extinction is likely to happen sooner than otherwise expected. The parameters for the simulations were chosen arbitrarily 1/μ = 10 (days), v = 1.33 (Km/h), R = .02 (m) and b = 3, although they resemble the characteristics in some relatively small insects. The value p = .01, which in practice depends on the complexity of the mating mechanisms, was also set arbitrarly.
females (Bazykin, 1998), during which a female avoids further sexual encounters just after successfully mating with a male. Let r denote the rate of mating encounters per female, r = KvRx = αx, then 1/r is the average time between encounters for one female. If the population densities are low enough so that 1/r is very much more than τ , then the population growth will depend on the number of successful mating encounters made. In this case, the average progeny produced by a female has to saturate when 1/r ↓ τ as the density of males increases. Therefore, the average birth rate per female is more conveniently approximated with bp × αx × 1 N+x , where b is the average number of offspring per female per encounter that survive to adulthood, p is the probability that an encounter produces offspring and N is the average male population density at which half of the females are able to reproduce. This can be seen by making x = N, so the average birth rate per female is bpα/2 or equivalently that half of the females produce offspring at the per capita rate bpα. For simplicity, we have set N = 1, suggested by the reasonable approximation of the mean number of pairs formed at low densities (Gordillo, 2015). Given that initially we assumed a 1:1 sex ratio, the whole population dynamics can be described by the equation where z = x + y and γ = (bpα − 2μ)/2 > 0. Thus, if z 2, the population growth will approximately correspond to a mass action regime while if z 2, it grows exponentially. More precisely, z < 0 when z < 2μ/γ , giving rise to positive density dependence with a critical population size.
We construct a stochastic version of the model by seeing the birth and death rates as stochastic rates (Gillespie, 1977); that is, we integrate demographic random fluctuations by considering only the two types of events, birth and death of individuals, happening at exponentially distributed times with the overall rate given by At a given time, a death happens with probability q = μz/ , or a birth with probability 1− q. In addition, environmental stochasticity is introduced by letting K fluctuate randomly, with K ∼ Normal(K, σ 2 ) (K = .1231). The effects of varying σ 2 are displayed in Figure 2, which shows the cumulative probabilities of extinction up to a fixed time, conditioned on extinction. The initial population density is taken at half of the Allee threshold value. The results suggest that an increasing variance in the random fluctuations associated with the mating encounters would increase the probability of earlier extinction events.

Discussion
Mating processes in two-sex insect populations are species specific and generally involve vast complexity (Bonduriansky, 2001;Choe & Crespi, 1997). This necessitates an approximated description of the non-linear process for mating encounters. Mass action has been continuously used for this end, even in cases with convoluted mating mechanisms (Courchamp et al., 2008). The multiplicative constant for the mating encounter rate might be difficult to estimate in many real applications. The novelty in this note is that, under general assumptions, we have approximated the value of that constant with the use of dimensional reduction and simulated data. Before now, this constant has not been computed by such analytical arguments. The dimensional reduction approach also offers a way to examine the effects of correlation in the movement. The simulations suggest that in general, the variances of the values of K are small and similar when compared among results obtained with different degrees of correlation. Thus, the valueK can be used as the multiplicative constant in the expression for the encounter rate,KvRAxy, independently of movement correlation. Relative larger variance values observed for highly correlated movement (θ close to zero) are attributed to the boundary conditions in the simulations. We remark that the programme used to produce the data only counts ephemeral encounters and does not consider further association between paired individuals, as well as other complexities involved in mating mechanisms.
While our analysis assumes at most one mating encounter per female at each time step, the approach may be generalized to systems where individuals must spend more or less time between mating encounters. For the purposes of modelling sex role evolution, H. Kokko and M. Jennions presented in Kokko and Jennions (2008) a scheme that accounts for time delays following encounters when individuals cannot mate. This framework could be adopted to modify and extend our work.
Small random fluctuations in the number of mating encounters can appear nevertheless, caused by non-permanent random environmental changes. These effects might alter substantially the encounter rate (now random). As an illustration, we simulated a basic stochastic birth and death process to obtain the probability distribution of the time to extinction, conditioned on extinction. The model has a deterministic analogue for which extinction due to Allee effects is possible. The simulations suggest that small fluctuations on the parameterK might induce relatively large increments on the (conditioned) cumulative probability of extinction, see Figure 2. For conservation or pest control efforts in which it is critical to assess how rapidly a population might become extinct, the approximation for K presented here might help shape more accurate quantitative predictions based on the law of mass action.
We would like to emphasize that, after all, mass action is a highly idealized model where physical and biological details are overlooked, including those corresponding to mating mechanisms. For instance, velocities associated with aerial mating in insects may be sensitive to the relative mass of flight muscle (Dudley, 2000), while (the radius of) attraction may depend on chemical, visual or auditory signals (Resh & Cardé, 2009;Roelofs, 1995). Trying to introduce these species-specific details into the mass action scheme could be a challenging task that we have not considered. In the case that velocity and radius of attraction were both scaled to individual total mass, for example, it would be possible to insert these relations into the main formulae as long as the assumptions for the mass action still hold. If the case is that individual mass implies non-negligible (relative) size, then the simulations cannot be used as we set them, and new computational experiments together with a re-assessment of the whole model have to be done. In contrast to the crude assumptions made for mass action, we refer the reader to the interesting papers (Gurarie & Ovaskainen, 2011, 2013, where the authors present a refined continuous timecontinuous space theoretical framework that shows how encounter rates depend on the interplay of spatial distributions, scales of movement, individuals densities and inherent dynamics. However, we believe the simple model presented here can inform the theoretical modelling of ephemeral (short-lived) mating encounters for insects, important in pest control management (Boukal & Berec, 2008;Fauvergue, 2013;Gordillo, 2015).