On the evolution of sexual receptivity in female primates

There has been much interest in the evolutionary forces responsible for, and underlying the diversity in, female primate reproductive cycles. While there has been limited research on sexual receptivity in primates, this has been one recurring topic of interest. Some primate species are like humans, sexually receptive to mating throughout their entire estrus cycle, while other species are the opposite, receptive for mere hours out of their several-week cycles. Why is there such prominent variation in sexual receptivity length among primate species? Here we examine the evolutionary trade-offs associated with sexual receptivity length using mathematical modeling. We investigate how various factors, including having ovulation signs present versus concealed ovulation, female physiological costs, and group size, each influence the length of females’ receptive periods. We find that both continuous receptivity and very short lengths of receptivity are able to evolve. Our model predicts that increasing the impacts of infanticide will increase the length of the female receptive period, emphasizing the possible importance of paternity confusion. Similar effects can also be achieved by increasing the non-genetic benefits provided by males. Overall, our work offers a theoretical framework for understanding the evolution and diversity of mating traits in female primates.

The purpose of the model in this article is to propose a general framework for modeling the evolution of sexual receptivity in female primates. Using this model, we will investigate the 712 following questions: How does length of sexual receptivity contribute to female reproductive success? What factors facilitate the evolution of continuous receptivity and/or very short lengths 714 of receptivity? What role could infanticide have in the evolution of receptivity? How is receptivity linked to the evolution of sexual attractiveness, namely a female having obvious visual 716 ovulation signs or concealed ovulation?

718
We consider a population of individuals interacting in G groups, each comprised of N males and N females. Females can differ genetically in both their visible ovulation signs present (with 720 magnitude m and length ) and their length of receptivity (r), while males differ in their quality to females, meaning any benefit from the male to a female and/or her offspring. We do not 722 consider evolution in males, assuming instead that male traits are at a [stochastic] evolutionary equilibrium.

724
Male-provided genetic (y g ) and non-genetic (y ng ) benefits are randomly drawn from the bivariate normal distribution, each with meanŷ and standard deviation b, and correlation param-726 eter ρ; parameterŷ characterizes mean male quality and b > 0 the extent of variation. Males are ranked according to the value of their genetic quality: y g,j such that y g,1 > y g,2 > ... for each male 728 j. Small b implies small additional benefits to females from investing in ovulation signaling.
Our model explicitly accounts for the female cycle by using D discrete units of time. Without 730 loss of generality, we refer to these units of time as 'days' (e.g., D = 29 days). For each day of the cycle, every female will have an associated probability of fertilization. We assume each 732 fertile period to last C ≤ D days (e.g., C = 7). For all days lying outside these C fertile days, females are assumed to have zero probability of fertilization but can still mate if receptive on 734 those days. We assume a triangular shape for the fertility function to be identical in all females, but with midpoints randomly distributed so as cycle synchrony occurs only probabilistically.

736
Visual ovulation signs and receptivity are both correlated with these days of fertility (in that a female's day of peak ovulation signs will align with both her median day of receptivity and her 738 peak day of fertility), but can each still evolve freely.
Although the form of such fertile periods will be identical in all females, the timing of each 740 fertile period will be randomly distributed among females. For instance, when C = 7 (calling each of these fertile days C1, C2, ... , C7) and D = 29 (calling each of these days D1, D2, ... , D29), 742 C1 in every female has an equal probability of landing on any of D1, D2, ... , D29. Since C is a cycle, we similarly assume C7 can land on any such day. For example, if C1 landed on D24, then 744 C7 would 'wrap around' the cycle to land on D1. Moreover, for all females we assume ovulation happens directly in the middle of this C cycle. Hence, for C = 7 we assume ovulation happens 746 on C4 and the probabilities of fertilization increase from C1 to C4 and then decrease from C4 to C7.

748
Each female i will have r i ≤ D days of receptivity. These days of receptivity line up with the female's days of fertility, such that a female's median day of fertility will equal her median 750 day of receptivity. Similarly, since having visible ovulation signs at least loosely correlates with a female's fertility, we assume the day(s) of having maximum visual ovulation signs also align with 752 the day(s) of that female's maximum fertilization probability, and align these cycles accordingly (note in most cases, C = r i = i ).

754
Let r be the number of days of the cycle a female is receptive to mating (a non-negative, integer value). We treat visual ovulation signs, x(d) for each day d of the cycle, as overlapping 756 graded curves. Each female's curve is characterized by two traits: magnitude (m) and length ( ). m is the maximum amount of ovulation signs a female has visible during her cycle (a non-758 negative, continuous value), while is the number of days a female has some amount of ovulation signs visible (a non-negative, integer value).

760
Given female traits m i ≥ 0 (ovulation signs magnitude) and i ∈ Z, i ∈ [0, D] (ovulation signs length), the amount x i (d) of ovulation signs visible on any day d of the cycle is defined as for parameter γ > 0 and d measured from the start of ovulation signs being visible. Mutation effects in m i are randomly chosen from the normal distribution N(0, σ 2 ), whereas mutations in 762 i and r i are discrete, representing adding or subtracting exactly one day to or from a female's time of having ovulation signs visible or being receptive, respectively.

764
Note that the exact function x i (d) is chosen such that (i) m i will denote the peak of the curve, for m i = 1, i = 7, and C = 7. Ovulation signs on each day are depicted by the red shaded bars, while fertility probabilities are given by the black lines.

A.1.3 Process Overview and Scheduling
The model proceeds in discrete, non-overlapping generations. Males and receptive females com-772 pete for mating opportunities. Mating is followed by offspring production and dispersal. Males can also engage in infanticide. Within each generation, reproduction occurs based on female fit-774 ness, which in turn is based on mating success during the female reproductive cycle. The relative contribution of a female to the offspring generation is proportional to her relative fitness payoff.

776
Hence, "fitness" can be thought of as the expected number of offspring surviving to the age of reproduction, or each female's relative reproductive success. For each trait, mutations occur with 778 a small probability. Females are also able to migrate between groups.
We consider three different sets of results regarding selection:

786
Individual differences in female-expressed traits (ovulation signs magnitude m, ovulation signs length , and receptivity r) emerge, as individuals that perform poorly are eliminated over time 788 and replaced by individuals adopting the most efficient strategy. The initial gene pool is initiated at random in the first generation.
790 Figure S1: Sample receptivity lengths are displayed for r i = 3, 7, 11 (from top to bottom) via the light blue shading. Each graph also includes ovulation signs present (with m i = 1, i = 7) via the red bars and fertilization probabilities (with C = 7) via the black lines. In the top graph, we see that a female can be fertile and/or have ovulation signs visible on days where she is not receptive. In the bottom graph, we instead see a female can be receptive on days where she does not have any ovulation signs visible and/or is not fertile. The middle graph depicts the situation where visible ovulation signs, fertility, and receptivity all line up with lengths of 7 days.

A.2.2 Adaptation
Each female is characterized by up to three evolving traits. r is the length of time a female is 792 receptive to mating (a non-negative, integer value). In addition, females may also be characterized by their visual ovulation signs, denoted x(d) for each day d of the cycle, thought of as 794 overlapping curves. Visual ovulation signs are characterized by two genetically-controlled traits.
Ovulation signs magnitude m is the maximum amount of ovulation signs a female has visible 796 during her cycle (a non-negative, continuous value), while ovulation signs length is the number of days a female has some amount of ovulation signs visible (a non-negative, integer value).

798
The evolution of each of evolving trait is governed by female fitness, or each female's relative reproductive success.

A.2.3 Fitness
The fitness function for female i in the model is as follows: Here, w 0 is baseline fitness, y g,i the genetic benefit provided by the male who fertilizes female i, 802ỹ ng = 1 D ∑ (all mates j) y ng,j the average non-genetic benefit provided by all female i's mates, p i,j the perceived paternity probability of male j for female i's offspring, c ·x i the costs to a female of 804 supporting her visual ovulation signs withx i = 1 D ∑ D d=1 x i (d), a female's average visual ovulation signs, and c r r i D the costs to a female of being receptive to mating for r i D proportion of her cycle.

806
Notice for the case when α = 0 and, thus, g = 0 (i.e., without the effects of infanticide), the fitness function above collapses into Note in this α = 0 (no infanticide) case, parameters β, κ, τ, and ω all become unnecessary.

A.2.4 Learning and Prediction
Individuals in the model do not change their behavior in response to experience. In addition, individuals in the model do not predict future consequences of their decisions.

A.2.5 Sensing
All sensing is local. Males and females can identify the opposite sex and can discriminate during 812 mate choice. Males will preferentially mate with females with more ovulation signs visible, and females will only mate with males on days when they are receptive. Males will also know which 814 females they mated with during a cycle, in order to "estimate" their probability of paternity.

816
Interaction between individuals occurs through both mating and infanticide. See "Submodels" for more details on infanticide. 818 We assume each female to mate once on every day on which she is receptive, and each male to be able to mate with no more than one female on any given day. This means if there are R 820 receptive females on any given day of the cycle, there will only be R males able to mate with a female on that day. Biologically, a male after mating with a female may guard her to ensure 822 no other male is able to mate with her (i.e., mate guard) and/or a male must do other activities besides mating/searching for mates (e.g., hunting/eating). Mating is followed by offspring 824 production and dispersal.
We assume that, for every day of the cycle, males of higher GC will preferentially mate with 826 females with more ovulation signs visible. Parameters m and f control the amount of stochasticity in the process of mating pair formation. On each day, mating pairs are formed by first 828 randomly perturbing both the male trait y g and female trait x by adding to each an independent, normally distributed random variable with standard deviation m ≥ 0 and f ≥ 0, respectively 830 (i.e., x = x + e x , y g = y g + e y , where e x ∼ N(0, 2 f ), e y ∼ N(0, 2 m )). Males and receptive females are then sorted according to these perturbed values and mating occurs between individuals of the same order. With m , f → ∞, mating pairs are formed completely independently of the values of x i and y g,j .

A.2.7 Stochasticity
Stochasticity is present in mating, reproduction, and mutation. Stochasticity in mating pairs 836 comes from m , f , as described above. Stochasticity in reproduction comes from each male mate having some probability of paternity, as described below. Both the occurrence and the effect

842
The population and gene pool are the two collective levels. The population is important, since female fitness is impacted by both her male mate(s) and infanticide. The gene pool acts through 844 the model's population of sexually-reproducing individuals.

846
In each generation, the evolving female-expressed traits -that is, r, m, and -of all individuals in the population are recorded.

850
For this model, recall we consider 3 different sets of simulations regarding selection: (1) Evolution of receptivity (r) given fixed concealed ovulation (m, = 0),

A.3.2 Input Data
This simulation model does not use input data.

A.3.3 Submodels
There are two submodels considered in the overall model. The first is calculating probabilities of 866 paternity for all males, and the second is infanticide by males.
To calculate the probability of paternity for each male, first the fertility of each female with 868 which that male has mated is summed up for every day a particular male mates with her. These values are normalized across males for each female in the group, meaning a male who mates 870 with the female on a day where she has a higher fertility probability will have a higher paternity probability. The actual father for that female's offspring is then determined randomly,

872
proportional to each male's probability of paternity for her offspring.
Note these actual probabilities of paternity are different from the perceived probabilities of 874 paternity, as detailed in the main text. We make this distinction because males will not know any female's probability of fertility at the time he mates with her; a male will only know how many 876 visible ovulation signs she has present. Separating these quantities allows actual probabilities of paternity (which males do not know) to be calculated using females' probabilities of fertility, and 878 perceived probabilities of paternity (which males do know) to be calculated using females' visible ovulation signs.

880
In addition, each male in a female's group can affect her fitness directly, due to infanticide, although males will differ in their actions' effectiveness (e.g., due to strength, rank, size, alliances, etc.). We assume each male's effectiveness to be proportional to f (j) for each male j (sorted by males' y g ) via an exponential function: f (j) ∼ e −ωj with parameter ω > 0 controlling the amount of disparity among males in their corresponding effectiveness within the group. Note a larger value of ω indicates more disparity, and ω = 0 equality. We define this exponential function of male effectiveness to be: A male's contribution g(p j ) to female fitness (positive via helping protect the offspring, or negative via not helping and/or harming the offspring) depends on his perceived probability of 882 paternity p j . We define g(

894
Parameters used in the model are outlined in Table S1.

C Effects of all Parameters
See Table S2 for a summary of effects for all main parameters used in the model. These effects are 900 all illustrated for simulations with concealed ovulation fixed ( Figures S2,S3,S4), visible ovulation signs fixed ( Figures S5,S6,S7), and each of receptivity length, ovulation signs magnitude, and 902 ovulation signs length evolving (Figures S8-S16).  Figure: Value(s) of ρ: Other Parameters Investigated:

C.1 ANOVA Results
We ran an analysis of variance (ANOVA) to determine which of receptivity length (r), ovulation 906 signs magnitude m, and ovulation signs length are affected most by which parameters, as introduced in the main text.

908
The following tables (Tables S3, S4, S5) give the detailed results of these tests. Each table below reflects the results from different simulation sets. The numbers in each of these tables 910 correspond to percentage of variance, with the sign (±) corresponding to the direction of the effect. Any table entry with a zero means that effect is not significant (i.e., p > 0.05). For 912 example, in Table S3, 31% of the variation in receptivity length (r) in this simulation set can be explained by parameter c r , and another 38% by parameter η. Since the effect of c r is negative and 914 η positive, we know that as c r increases, r decreases, and as η increases, r instead increases. Each of the following tables can be interpreted in this fashion.          Figure S11: The effects of parameters α (maximum benefit of infanticide), β (proportional to the maximum cost of infanticide), κ (weight males put on females having ovulation signs visible), and η (relative weighting of NGC) on the average equilibria values of receptivity length (r), ovulation signs magnitude (m), and ovulation signs length ( ) with ρ = 0.5 (correlation between males' GC and NGC). Equilibria are obtained by averaging over 16 initial condition runs (with standard deviation indicated by error bars). All other parameters were held constant: N = 8, b = 0.2, c = 0.2, c r = 0.6, m = 0.25, f = 0.01. Figure S12: The effects of parameters α (maximum benefit of infanticide), β (proportional to the maximum cost of infanticide), κ (weight males put on females having ovulation signs visible), and η (relative weighting of NGC) on the average equilibria values of receptivity length (r), ovulation signs magnitude (m), and ovulation signs length ( ) with ρ = 0 (correlation between males' GC and NGC). Equilibria are obtained by averaging over 16 initial condition runs (with standard deviation indicated by error bars). All other parameters were held constant: N = 8, b = 0.2, c = 0.2, c r = 0.6, m = 0.25, f = 0.01. Figure S13: The effects of parameters α (maximum benefit of infanticide), β (proportional to the maximum cost of infanticide), κ (weight males put on females having ovulation signs visible), and η (relative weighting of NGC) on the average equilibria values of receptivity length (r), ovulation signs magnitude (m), and ovulation signs length ( ) with ρ = −0.5 (correlation between males' GC and NGC). Equilibria are obtained by averaging over 16 initial condition runs (with standard deviation indicated by error bars). All other parameters were held constant: N = 8, b = 0.2, c = 0.2, c r = 0.6, m = 0.25, f = 0.01. Figure S14: The effects of parameters m (reproductive stochasticity among males), f (reproductive stochasticity among females), c r (cost of receptivity length), and η (relative weighting of NGC) on the average equilibria values of receptivity length (r), ovulation signs magnitude (m), and ovulation signs length ( ) with ρ = 0.5 (correlation between males' GC and NGC). Equilibria are obtained by averaging over 16 initial condition runs (with standard deviation indicated by error bars). All other parameters were held constant: N = 8, b = 0.2, c = 0.2, α = 0. Figure S15: The effects of parameters m (reproductive stochasticity among males), f (reproductive stochasticity among females), c r (cost of receptivity length), and η (relative weighting of NGC) on the average equilibria values of receptivity length (r), ovulation signs magnitude (m), and ovulation signs length ( ) with ρ = 0 (correlation between males' GC and NGC). Equilibria are obtained by averaging over 16 initial condition runs (with standard deviation indicated by error bars). All other parameters were held constant: N = 8, b = 0.2, c = 0.2, α = 0. Figure S16: The effects of parameters m (reproductive stochasticity among males), f (reproductive stochasticity among females), c r (cost of receptivity length), and η (relative weighting of NGC) on the average equilibria values of receptivity length (r), ovulation signs magnitude (m), and ovulation signs length ( ) with ρ = −0.5 (correlation between males' GC and NGC). Equilibria are obtained by averaging over 16 initial condition runs (with standard deviation indicated by error bars). All other parameters were held constant: N = 8, b = 0.2, c = 0.2, α = 0.