Mating, births, and transitions: a flexible two-sex matrix model for evolutionary demography

Mating, births, and transitions: a flexible two-sex matrix model for evolutionary demography Models of sexually-reproducing populations that consider only a single sex cannot capture the effects of sex-specific demo- graphic differences and mate availability. We present a new framework for two-sex demographic models that implements and extends the birth-matrix mating-rule approach of Pollak. The model is a continuous-time matrix model that explicitly includes the processes of mating (which is nonlinear but homogeneous), offspring production, and demographic transitions and survival. The resulting nonlinear model converges to exponential growth with an equilibrium population composition. The model can incorporate age- or stage-structured life histories and flexible mating functions. As an example, we apply the model to analyze the effects of mating strategies (polygamy or monogamy, and mated unions composed of males and females, of variable duration) on the response to sex-biased harvesting. The combination of demographic complexity with the interaction of the sexes can have major population dynamic effects and can change the outcome of evolution on sex- related characters. Keywords Birth matrix-mating rule · BMMR · Demography · Matrix population models · Sex-biased harvest · Two-sex models Introduction mortalities (Kuczynski 1932; Pollak 1990; Jenouvrier et al. 2010), developmental schedules (Caswell 2001), behavioral Models of sexually-reproducing populations that consider interactions (Rankin and Kokko 2007), dispersal patterns a only single sex (typically females) implicitly assume that (Miller et al. 2011), and selective harvest pressures (Ginsberg both sexes have identical vital rates and that the availability and Milner-Gulland 1994). Additionally, ecological, environ- of the neglected sex (typically males) does not affect fertility mental, and evolutionary changes may alter the most limiting (Pollard 1974; Caswell 2001; Iannelli et al. 2005). In reality, sex over time (Hardy 2002; Miller and Inouye 2011). One- both these assumptions are frequently violated. Males and sex models miss, and cannot explore, important components females often differ significantly in terms of fertilities and of population dynamics in all these cases. Here, we introduce a flexible framework that can explore many of these factors. In response to growing concerns over discrepancies in male and female reproductive rates (Karmel 1947), early Electronic supplementary material The online version of this article (https ://doi.org/10.1007/s1014 4-018-0615-8) contains dynamical models with sex structure were introduced in the supplementary material, which is available to authorized users. late 1940s. Pollard (1948) used coupled Lotka integral equa- tions, considering female births to males and male births to * Hal Caswell females in order to reconcile the growth rates of both sexes. h.caswell@uva.nl Kendall (1949) introduced a system of ordinary differential Biology Department MS-34, Woods Hole Oceanographic equations for males and females (and later, married cou- Institution, Woods Hole, MA 02543, USA ples), and was the first to incorporate nonlinear interactions Institute for Biodiversity and Ecosystem Dynamics, between the sexes via a mating term. Subsequent models University of Amsterdam, PO Box 94248, considered other nonlinear mating functions (Pollard 1974; 1090 GE Amsterdam, The Netherlands Vol.:(0123456789) 1 3 22 Population Ecology (2018) 60:21–36 Yellin and Samuelson 1974) and couple dissolution through Model development death and divorce (Hadeler et al. 1988). Extensions to age-structured populations were made by For reasons that will become apparent, we have incorporated Fredrickson (1971) and Hoppensteadt (1975), who allowed sex structure, stage structure, and life cycle processes into birth and death rates, as well as couple formation and divorce a continuous-time, rather than a discrete-time matrix popu- rates, to depend on age and sex. Hadeler later included the lation model. The mating, birth, and transition processes age (i.e., duration) of married pairs (Hadeler et al. 1988) from the BMMR framework are described by separate rate and maturation delays (Hadeler 1993). Age-structured mat- matrices. The mating process introduces nonlinearity into ing functions have similarly been proposed (Martcheva and the model because reproduction depends on the relative pro- Milner 2001). Many of these models use continuous-time portions of males and females, of appropriate stages, in the equations that incorporate age structure through coupled population. This dependence can be flexibly modeled by McKendrick–von Foerster partial differential equations (Fre- generalized weighted mean mating functions, which satisfy drickson 1971; Keyfitz 1972; Hoppensteadt 1975; Hadeler the biological criteria of sexual reproduction (Iannelli et al. 1989, 1993), though discrete-time two sex matrix models 2005). The resulting BMMR matrix models are nonlinear have also been developed (Caswell and Weeks 1986). but homogeneous (i.e., frequency-dependent). Models of A powerful conceptual approach to two-sex models was this form generally converge to an exponential growth rate proposed by Pollak (1986, 1987, 1990) and called the birth and a stable stage frequency distribution (Martcheva 1999). matrix-mating rule (BMMR) model. It proposes that mating, births, and life cycle transition processes repeat periodically, Incorporating sex and stage structure one after the other. It contains three main components: The model classifies individuals into stages based on age, 1. A mating rule function that gives the number of matings developmental state, sex, reproductive status, or other vari- u between males of age i and females of age j. ables of interest. Stage densities are projected forward in ij 2. A birth matrix whose entries b are the expected number time by a projection matrix that contains the demographic ij of male and female offspring produced by a mating of a rates characterizing survival, reproduction, and transitions male of age i and a female of age j. between stages. The properties of this projection matrix 3. Sex-specific mortality schedules, which project surviv - provide information about the population as a whole, ing individuals to the next age class, or, in our generali- thereby linking individual-level life cycle information (i.e., zation, include other stage-specific life cycle transitions. the stage-specific vital rates in the matrix entries) to popu- lation-level properties important for ecology and evolution BMMR is a useful approach for describing two-sex popu- (e.g., growth rates or stage distributions). lations because it can specify age (and, more generally, A population with s stages is described by a s × 1 popula- stage) structure over all parts of the life cycle. This struc- tion vector (t) , the entries of which are the densities of each ture, in turn, can have significant effects on two-sex popula- stage at time t. In a two-sex population, (t) would contain tion dynamics (e.g., Sundelöf and Åberg 2006, where the male stages, female stages, and mated stages (unions) that addition of size-specific birth functions affects growth rate could include married couples or breeding harems. and reproductive output) and recommended management For example, the population vector for a two-sex popula- strategies (e.g., Ginsberg and Milner-Gulland 1994, where tion with mating adults and nonmating juveniles could have incorporating age-specific fecundity changes the outcomes the form: of sex-biased harvest). Conceptually appealing as it is, the BMMR model has yet ⎛ ⎞ juvenile males ⎜ ⎟ to be fully incorporated into the framework of stage-classie fi d adult males ⎜ ⎟ matrix population models. Our goal here is to do so, providing (t)= juvenile females (1) ⎜ ⎟ a general model that can be adapted to a wide range of life ⎜ ⎟ adult females ⎜ ⎟ cycles and mating strategies. We do so by a novel extension adult unions ⎝ ⎠ of periodic matrix models to continuous time, based on transi- tion rate matrices. In this paper, we focus on ecological impli- cations of the model, but the approach can also be applied to The dynamics of the population vector are given by a system of ordinary differential equations study the evolution of sex-related traits using methods from adaptive dynamics (Shyu and Caswell 2016a, b). = [] (t) (2) dt 1 3 Population Ecology (2018) 60:21–36 23 where the entries of  are either transition rates or rates of offspring production, and we have indicated that they may depend on the population vector. Incorporating the BMMR processes The BMMR framework incorporates mating, birth, and transition processes. We describe each of these processes by a separate matrix: 1. The mating (union formation) process, where adult Fig. 1 Mating functions from the generalized weighted mean family males and females organize into reproductive unions, is (Eq.  6) with a check to indicate which of the biologically desirable described by the matrix . mating function criteria they satisfy 2. The birth process, where unions produce new offspring, is described by the matrix . 3. The transition process, which includes other life cycle M() events like mortality, maturation, or divorce, is described U ()= (4) by the matrix . M() Other life cycle processes can be included with additional U ()= . (5) matrices. A discrete-time model would incorporate these processes The total population mating rate is M = U m = U f . m f as a periodic matrix product (Caswell and Shyu 2012). For Many commonly used mating functions are generalized example, the product  would describe mating, followed weighted means (Hölder means) of the form by the production of offspring from the matings, followed 1∕ by survival and transitions of the resulting individuals. In M()= f +(1 − )m (6) continuous-time, we conceptualize the processes as occur- where  and  are constants; 0 ≤  ≤ 1 , 𝛼< 0 (Hadeler 1989; ring simultaneously. It can be shown (Appendix 1) that the Bessa-Gomes et al. 2010; Iannelli et al. 2005). Figure  1 projection matrix  in the continuous-time matrix model shows several generalized weighted mean mating functions (Eq. 2) is the average of the transition rate matrices, e.g.,: and biologically desirable criteria that they satisfy (McFar- d 1 land 1972; Pollard 1974; Yellin and Samuelson 1974). = ( +  + )(t) dt 3 (3) If multiple male and female stages interbreed to form = (t) different types of unions, stage-specific mating preferences can also be integrated into this mating function (Shyu and Caswell 2016b; see  Appendix 4). Modeling the mating process It is difficult to distinguish between mating functions in populations where the sex ratio does not vary signifi- The mating process, as described by the union formation matrix , depends on the relative numbers of males and females in the cantly (Keyfitz 1972). However, recent empirical studies (Miller and Inouye 2011) support the harmonic mean as a population, not all of which may be mature enough or available for breeding (Pollard 1974; Iannelli et al. 2005). As a result,  mating function. Because the harmonic mean satisfies rea- sonable biological criteria for mating functions (Caswell depends on the population’s sex and stage composition, making a function of the population vector (t). and Weeks 1986; Iannelli et al. 2005), and captures the qualitative properties of other generalized means, we will The total mating function M() gives the rate of union forma- tion (total number of unions formed per unit time in a population use a harmonic mean mating function for our analyses. of composition  ). This embodies the mating rule in the BMMR framework. Here, “unions” refers to any mated, reproducing Frequency‑dependent dynamics units in the population, including both one-to-one male-female pairs and harems with multiple individuals of the one sex. The mating process often depends on the relative frequen- cies, rather than absolute abundances, of males and females. We convert the total mating function into the per capita mating rates U () and U () (the average mating rates per As a result, although the mating function is nonlinear, it is m f homogeneous of degree 1 in  . That is: available males m or females f respectively), 1 3 24 Population Ecology (2018) 60:21–36 M(c)= cM() (7) for any positive constant c. As a result, the per capita mating functions (Eqs. 4 and 5) are homogeneous of degree 0 in  , so that: U (c)= U () m m (8) U (c)= U () f f If all entries in the projection matrix  in Eq.  3 are also Fig. 2 Life cycle diagram for a 5-stage population with juvenile males m and juvenile females f , adult males m and adult females homogeneous of degree 0, the system is said to be  fre- 1 1 2 f , and reproducing unions u. The functions and parameters shown quency-dependent. This means that  can be written as a here, as described in Table 1, appear in the union formation matrix function of the population frequency vector: (Eq. 14) (red), birth matrix  (Eq. 15) (green), or transition matrix (Eq. 16) (blue). From Shyu and Caswell (2016a) (9) ‖‖ where ‖‖ is the 1-norm of . 2m f 2 2 Frequency-dependent models of this type usually con- M()= m + f 2 2 verge asymptotically to an equilibrium population structure 2f ̂ (the stable stage distribution) in which all stage frequen- U ()= (12) m + f 2 2 cies are constant (e.g., see Yellin and Samuelson 1974; 2m Hadeler 1989; Martcheva 1999). The population then grows 2 U ()= m + f or decays exponentially at a long-term growth rate given by 2 2 the dominant eigenvalue  of [ ̂ ]. To find the equilibrium stage distribution  and popula- Again, we consider the life cycle in terms of mating, birth, tion growth rate  , it is sufficient to consider the dynamics and transition processes, which are described by matrices  , of (t) . It can be shown (Appendix 2) that: , and  respectively. =  −  []  (10) dt 1. The union formation matrix  contains the per capita mating functions from Eq. 12. where  is a s × s identity matrix and  is a 1 × s vector of ones. One can integrate Eq. 10 with a numerical differential ⎛000 0 0⎞ equation solver until population frequencies converge to  ̂ . ⎜ ⎟ 0 − U () 0 00 ̂ m Then  is the dominant eigenvalue of [] . The dominant ⎜ ⎟ 000 0 0 ()= (13) right eigenvector  of [ ̂ ] is proportional to the stable stage ⎜ ⎟ ⎜000 − U () 0⎟ distribution . ⎜ 1 1 ⎟ 0 U () 0 U () 0 ⎝ m f ⎠ 2 2 A 5‑stage BMMR matrix model ⎛0 00 00⎞ M() ⎜ ⎟ 0 − 0 00 We now present an example of a BMMR matrix model with ⎜ ⎟ five stages: juvenile males m and juvenile females f , adult ⎜ ⎟ 0 00 00 1 1 (14) ⎜ M() ⎟ males m and adult females f , and reproducing unions u that 0 00 − 0 2 2 ⎜ ⎟ consist of one adult male and one adult female. Single adult M() M() ⎜ ⎟ 0 0 0 ⎝ ⎠ 2m 2f males and females interact to form unions, which then produce 2 2 new juvenile offspring (Fig.  2). A summary of the variables, parameters, and matrices in this model is provided in Table 1. Note that U and U must halved in the last row of  to m f As in Eq. 1, we write the population vector as avoid double counting the unions formed from both male and female partners. (t)= m m f f u (11) 1 2 1 2 2. The birth matrix  contains k, the average reproductive rate of a union, and the primary sex ratio s , the propor- Using the harmonic mean mating function, the total and per tion of offspring that are male. capita mating functions are 1 3 Population Ecology (2018) 60:21–36 25 a b 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 01 23 45 67 8 01 23 45 67 8 time time Fig. 3 Dynamics of the 5-stage BMMR model with monogamous (11), where dynamics are given by (3). b) Convergence of the pop- unions and no harvest. The population consists of juvenile males m ulation frequency vector  (9), where dynamics are given by (17). and juvenile females f , adult males m and adult females f , and Parameters are fixed at  =  = 0.5,  =  = 0.1,  =  = 1 2 2 m1 f 1 m2 f 2 m f reproducing unions u. a) Growth of the population density vector  0.5, s = 0.5, k = 20, d = 0.1, h = 1 , E = 0. Table 1 A summary of the Matrices and vectors variables, parameters, matrices, Projection matrix and population properties in the 5-stage BMMR matrix model Birth matrix (15) Transition matrix (16) Union matrix (14) Population density vector (11) Population frequency vector (9) ̂ or  Equilibrium stage structure Population properties Long-term population growth rate, dominant eigenvalue of [ ̂ ] s Primary sex ratio (proportion of offspring that are born male) s Secondary sex ratio (proportion of adults that are male) f , f Juvenile, adult female density 1 2 m ,m Juvenile, adult male density 1 2 u Union (mated pair) density Life cycle parameters ,  Male, female maturation rates m f d Divorce rate (rate at which a male-female pair bond breaks) h Average harem size k Union reproductive rate ,  Juvenile, adult female mortality rates f 1 f 2 ,  Juvenile, adult male mortality rates m1 m2 M Total mating rate (total unions formed per time) U , U Per capita mating rates (4) and (5) m f Harvest parameters E Total harvest rate in (18) s Harvested sex ratio (proportion of harvest that targets males) in Eq. 18 1 3 stage densities stage frequencies 26 Population Ecology (2018) 60:21–36 Mating systems and the effects of sex‑biased ⎛ ⎞ 00 00 ks harvest ⎜ ⎟ 00 00 0 ⎜ ⎟ = 00 00 k(1 − s ) (15) As an example of the use of the two-sex model, we consider ⎜ 1 ⎟ ⎜ ⎟ 00 00 0 the effectes of sex-biased harvesting. In sport or trophy hunt- ⎜ ⎟ 00 00 0 ing, harvest is often male-biased and significantly exceeds ⎝ ⎠ natural mortality (Festa-Bianchet 2003; Milner et al. 2007). Age or size bias is also common, as larger or older males 3. The transition matrix  contains the male mortality rates with well-developed adult characteristics (e.g., large antlers (  for juveniles,  for adults) and female mortal- m1 m2 or horns) are usually the most desirable targets. This unnat- ity rates (  for juveniles,  for adults), the rates of f 1 f 2 ural selection may alter population structure, reproductive maturation from juveniles to adults (  for males,  for m f strategies, body morphology, and developmental timing (e.g., females), and the divorce rate d (rate at which the male- Ashley et al. 2003; Festa-Bianchet 2003; Allendorf and Hard female pair bond breaks). Note that unions may also 2009). The population response depends on multiple demo- dissolve due to partner death, and that union dissolution graphic factors, including the mating system. through both death and divorce returns surviving males The mating system determines how males and females and females to the single adult stages. organize for breeding and is thus a key component of two- sex population structure (Emlen and Oring 1977). Some species form monogamous, one-to-one pair bonds between ⎛−( +  ) 0 00 0 ⎞ m1 m ⎜ ⎟ males and females. Other species have polygynous mating −  00  + d ⎜ m m2 f 2 ⎟ systems in which one male mates with multiple females ⎜ ⎟ = 00 −( +  ) 00 f 1 f ⎜ ⎟ (e.g., lions, deer, seals), or, more rarely, polyandrous systems ⎜ ⎟ 00  −   + d f f 2 m2 ⎜ ⎟ where one female mates with multiple males (e.g., jacanas, ⎝ 0 0 00 −( +  + d)⎠ m2 f 2 pipefish). Unions formed by mating may be transient and (16) limited to a single breeding episode (e.g., lek systems) or The advantages of the continuous-time formulation are may persist over multiple breeding seasons (e.g., lion har- apparent at this point. The rates of maturation, mortality, ems) or even until partner death (e.g., albatrosses and other and divorce in Eq. 16 combine in a simple additive fashion. species with high mate fidelity) (Cézilly and Danchin 2008). Transition probabilities would not combine additively; they These factors motivate the use of a demographic two-sex would involve sums of products of conditional probabilities, model in analyzing harvest effects. To this end, we will use over all the possible events. The number of these products our BMMR matrix framework to explore the effect of mating increases dramatically when more stages, and hence more systems on the response to sex-biased harvest. A range of mating combinations, are included. See Shyu and Caswell mating systems will be approximated by varying two model (2016b) and Appendix 4 for incorporation of multiple mat- parameters, d and h (Table 2). ing stages in the continuous-time model. As per Eq. 3, the average of these three matrices is the – The divorce rate d. This is a measure of union transience. continuous-time projection matrix [] . The corresponding Unions with higher values of d are more likely to dissolve equation for the proportional structure (Eq. 10) is thus: after a given mating, while unions with lower values of d are more likely to persist over multiple breeding interactions. =  −  ( +  + []) (17) – The harem size h. This is a measure of polygamy. Unions dt 3 with h = 1 are monogamous and consist only of one-to- where  is given by Eq. 14,  is given by Eq. 15, and  is one male-female pair bonds, while unions with h > 1 are given by Eq. 16. polygamous groups of size h + 1 . As polyandrous mating As shown in Fig. 3, the population vector  ultimately systems are relatively rare (Cézilly and Danchin 2008), grows exponentially, while the frequency vector  ultimately we will consider only the polygynous form of polygamy, converges to the constant distribution of stages  ̂ . To deter- where one male mates with h females. mine the equilibrium stage distribution  ̂ , we integrated Eq. 17 with the Matlab ODE45 differential equation solver Harvest strategies are characterized by the overall har- until population frequencies converged (e.g., until vector vest rate E and the harvested sex ratio s (proportion of entries do not change significantly over consecutive inte- the total harvest rate that targets males). Assuming that gration intervals). only adults are harvested, the adult mortality rates are The Matlab code used for all the following examples is modified by harvest as follows: provided in the Electronic Supplementary Material. 1 3 Population Ecology (2018) 60:21–36 27 Table 2 Mating systems Low d (persistent unions) High d (transient unions) corresponding to different values of the divorce rate d and h = 1 (monogamy) Persistent pair bonds, high mate fidelity (e.g., Serial pair bonds (e.g., harem size h albatross) humans, Emperor penguins) h > 1 (polygyny) Persistent harems (e.g., lion prides) Leks, scramble competi- tion (e.g., grouse, cod, horseshoe crabs) c change in λ from harvest d change in s from harvest -0.125 0.4 0.4 a 0.9 0.9 -0.13 0.3 0.8 0.8 0.3 -0.135 0.2 0.7 0.7 -0.14 0.1 0.2 0.6 0.6 -0.145 0.5 0.5 -0.15 0.4 0.4 0.3 b -0.1 -0.155 0.3 0.3 0.25 -0.2 -0.16 0.2 0.2 0.2 -0.3 0.1 0.1 -0.165 0.15 0 0 -0.4 01 23 45 01 23 45 01 23 45 divorce rate d divorce rate d divorce rate d Fig. 4 Population growth rates, structure, and responses to adult har- the divorce rate d. c) The change in  and d) the change in second- vest in the monogamous ( h = 1 ) model. a) Proportion of mated adults ary sex ratio s , when a proportion s of harvest targets males. With- 2 h (adults in stage u, rather than m or f ) for an unharvested population, out harvest, s = 0.5 for all values of d. Other parameters are fixed at 2 2 2 and the b) corresponding population growth rates  , as functions of  =  = 0.5,  =  = 0.1,  =  = 0.5, s = 0.5, k = 20, E = 1 m1 f 1 m2 f 2 m f 1 As shown in Fig.  4a, the proportion of adults in the →  + Es m2 m2 h reproductive union stage (mated adults) decreases as d (18) →  + E(1 − s ) increases. The unharvested population growth rate  simi- f 2 f 2 h larly decreases (Fig. 4b), because populations with more To determine how various mating systems, as characterized transient couples (fewer mated adults) cannot produce off- by different values of h and d, respond to sex-biased harvest, spring as rapidly as populations with more persistent cou- we will examine harvest effects on the long-term population ples (more mated adults). Because males and females have growth rate  and the secondary sex ratio s (proportion of the same baseline vital rates, the secondary sex ratio s 2 2 all adults that are male). We will assume that males and remains unbiased ( s = 0.5 , not shown) for all values of d. females have the same baseline vital rates, and that the pri- Figure 4c shows how increasingly sex-biased harvest mary sex ratio (proportion of males at birth) is 0.5. Thus, strategies affect population growth. Both biased and unbi- the main sex-specic fi die ff rences we consider are sex-biased ased harvest strategies most strongly reduce growth in harvest pressures and male versus females roles within the populations with lower divorce rates, as adult mortality polygynous mating systems. will also disrupt pair bonds. The greatest decreases in occur when harvest is strongly sex-biased, i.e., s is close Monogamy ( h = 1) to 0 (only females are harvested) or 1 (only males are har- vested). This suggests that monogamous populations with Consider a monogamous two-sex population with juve- high fidelity pair bonds will be the most impacted by sex- niles and adults. The mating process forms unions that are biased harvest, and that concentrating harvest on a single one-to-one pair bonds of adult males and females. This sex will more greatly reduce population growth. scenario uses the rate matrices  ,  , and  as given by Figure 4d shows how the same harvest strategies decrease Eqs. 14, 15, and 16 respectively, and the mating functions the secondary sex ratio s from equality. Unsurprisingly, male- in Eq. 12. We vary the divorce rate d to explore the effects biased strategies reduce s , female-biased strategies increase of transient vs. persistent pair bonds. s , and unbiased strategies ( s = 0.5 ) have relatively minimal 2 h 1 3 % adults mated λ, no harvest harvested sex ratio s harvested sex ratio s h 28 Population Ecology (2018) 60:21–36 effect. Populations with greater divorce rates experience larger reductions in s , possibly due to their lower growth rates (Fig. 4a) being unable to replenish harvested individuals as rapidly. Polygyny ( h > 1) Now consider a polygynous population that forms unions con- sisting of one male with a harem of females. Because the death or departure of a single female changes the harem’s size and reproductive rate, we must now account for multiple union (harem) stages. The stage u represents polygynous unions consisting of one male and a harem of i females. When h is the maximum harem size, the population vector contains h union stages, which Fig. 5 Reproduction and three possible transitions for u , a union with range from a union with a harem of size 1 ( u , equivalent to 1 harem size i a monogamous couple) to a union with a harem of size h (u , the largest possible union): (t)= m m f f u u … u (19) 1 2 1 2 1 2 h Assume that when males and females mate, they form the largest possible union u . Their union formation rate is still given by the harmonic mean mating function Eq. 12, but with the number of single females now replaced by the num- ber of prospective harems: Fig. 6 Stages in a population with maximum harem size h, which include juveni le males m and juvenile females f , adult males m and 1 1 2 2m adult females f , and reproducing unions u (one male with a harem 2 i M()= of i females). Adults form harems of size h when mating, and these m + harems can shrink in size, but not grow, over time 2f (20) U ()= hm + f 2 2 2m 2. A female harem member dies (with mortality rate  ). f 2 U ()= hm + f This shrinks the union from u to u . For the union 2 2 i i−1 u , which has only one female, u = u corresponds to 1 i−1 0 the single adult male stage m (i.e., the death of a wife returns her husband to the pool of unmated singles). Note that the total union formation rate M() is maximized 3. A female harem member departs from the union, with when the sex ratio of single adults is divorce rate d. We assume that only females leave, which 2 seems biologically plausible. Her departure returns one (21) m + f female to f and shrinks the union from u to u . For 2 2 1 + h 2 i i−1 the union u , divorce dissolves the union and returns one male to m and one female to f . Thus, as the harem size h increases, a higher proportion of 2 2 single females is needed to maximize the mating rate. As a result, a union may shrink (but not grow) in size due to If an individual female has a reproductive rate of k, the departure or death of its members (Fig. 6). After a union a harem with i females has a total reproductive rate of shrinks to the smallest possible size h = 1 , or if the male harem ik; larger harems are thus more productive. Each union leader dies, the union dissolves and its members return to the u , regardless of harem size, can change in three possible stages for unmated adults. ways (Fig. 5): Appendix 3 shows how to write the rate matrices  ,  , and for a polygynous system with maximum harem size h. The 1. The male harem leader dies (with mortality rate  ). m2 population vector and matrices for the case where h = 3 are This returns i adult females to the stage f . as follows: 1 3 Population Ecology (2018) 60:21–36 29 a λ before harvest b s before harvest 0.497 14 14 0.4 0.496 0.38 12 12 0.495 0.36 10 10 0.34 0.494 0.32 8 8 0.493 0.3 6 6 0.492 0.28 0.26 0.491 4 4 0.24 0.49 2 2 01 23 45 01 23 45 divorce rate d divorce rate d Fig. 7 Unharvested population dynamics in the polygynous ( h > 1 ) union model, as functions of the divorce rate d and harem size h. a) Popula- tion growth rate  . b) Secondary sex ratio s (proportion of males in the adult population). Other parameters are the same as in Fig. 4 Figure 7 shows how the unharvested population rate of (t)= m m f f u u u (22) 1 2 1 2 1 2 3 growth and secondary sex ratio vary across different values of d and h. As in the monogamous case, lower divorce rates ⎛ ⎞ 000 0 0 0 0 result in more mated reproducing adults and, thus, higher ⎜ ⎟ 0 − U 0 0 0 0 0 population growth. Larger harems lead to more rapid popu- ⎜ ⎟ 000 0 0 0 0 ⎜ ⎟ lation growth, possibly because of their higher total repro- ⎜ ⎟ 000 − 3U 00 0 f (23) ductive rates. Unions with a high maximum harem size are ⎜ ⎟ 000 0 0 0 0 more resilient to divorce and female mortality, because they ⎜ ⎟ ⎜000 0 0 0 0⎟ can lose more females before shrinking to a u and then 1 1 ⎜ ⎟ 0 U 0 U 00 0 dissolving. m f ⎝ ⎠ 2 2 Even without sex-biased harvest, the secondary sex ratio is slightly female-biased ( s ≈ 0.494 ), but varies only a few tenths of a percentage point across a wide range of h and d. 00 00 ks 2ks 3ks ⎛ 1 1 1 ⎞ Populations with high h and low d (large, persistent harems) ⎜00 00 0 0 0 ⎟ are the most biased. ⎜ ⎟ 00 00 k(1 − s ) 2k(1 − s ) 3k(1 − s ) 1 1 1 Figure 8 demonstrate how female-biased ( s = 0 ), unbi- ⎜ ⎟ h = 00 00 0 0 0 (24) ⎜ ⎟ ased ( s = 0.5 ), and male-biased ( s = 1 ) harvest strategies h h ⎜ ⎟ 00 00 0 0 0 affect population growth rates and secondary sex ratios. ⎜ ⎟ 00 00 0 0 0 ⎜ ⎟ ⎝ ⎠ 00 00 0 0 0 1. Female-biased harvest (Fig. 8a) most strongly reduces growth in populations with large d and small h, the same ⎛ ⎞ −( +  ) 0 00 0 0 0 m1 m ⎜ ⎟ −  00  + d 00 m m2 f 2 ⎜ ⎟ 00 −( +  ) 00 0 0 ⎜ ⎟ f 1 f ⎜ ⎟ (25) = 00  −   + d 2 + d 3 + d f f 2 m2 m2 m2 ⎜ ⎟ 0 0 00 −( +  + d)  + d 0 m2 f 2 f 2 ⎜ ⎟ ⎜ 0 0 00 0 −( +  + d)  + d ⎟ m2 f 2 f 2 ⎜ ⎟ 0 0 00 0 0 −( +  + d) m2 f 2 ⎝ ⎠ 1 3 harem size h harem size h 30 Population Ecology (2018) 60:21–36 change in λ from harvest change in s from harvest a female-biased harvest (s = 0) 0.28 14 14 -0.08 0.26 -0.09 12 12 0.24 0.22 -0.1 10 10 0.2 -0.11 0.18 -0.12 0.16 6 6 0.14 -0.13 0.12 4 4 -0.14 0.1 -0.15 2 2 01 23 45 01 23 45 divorce rate d divorce rate d b unbiased harvest (s = 0.5) -0.08 14 14 -0.085 -0.015 -0.09 -0.02 -0.095 10 10 -0.1 -0.025 -0.105 -0.03 -0.11 -0.115 -0.035 4 4 -0.12 -0.125 -0.04 01 23 45 01 23 45 divorce rate d divorce rate d c male-biased harvest (s = 1) -0.15 14 -0.06 14 -0.16 -0.17 12 12 -0.07 -0.18 10 10 -0.19 -0.08 -0.2 8 8 -0.09 -0.21 -0.22 6 6 -0.1 -0.23 4 4 -0.24 -0.11 -0.25 2 2 01 23 45 01 23 45 divorce rate d divorce rate d Fig. 8 Population responses to harvest that is a female-biased rate  . (Right) The change in secondary sex ratio s . Other parameters ( s = 0 ), b unbiased ( s = 0.5 ), and c male-biased ( s = 1 ) in the are the same as in Fig. 4 h h h polygynous ( h > 1 ) model. (Left) The change in population growth 1 3 harem size h harem size h harem size h harem size h harem size h harem size h Population Ecology (2018) 60:21–36 31 0.4 Fig. 9 Growth rates  as a male-biased harvest (s = 1) function of the total harvest unbiased harvest (s = 0.5) 0.3 rate E for populations with female-biased harvest (s = 0) various mating systems. The 0.2 four types of mating systems shown correspond to different 0.1 harem sizes h and divorce rates low h, low d low h, high d d (Table 2); in this example, low h = 2 , high h = 10 , low d = 0 , 0.4 and high d = 2 . Harvest may be female-biased ( s = 0 ), unbi- 0.3 ased ( s = 0.5 ), or male-biased ( s = 1 ). Other parameters are h 0.2 the same as in Fig. 4 0.1 high h, low d high h, high d 00.2 0.40.6 0.81 1.21.4 1.61.8 2 00.2 0.40.6 0.81 1.21.4 1.61.8 2 E E populations with the lowest unharvested growth rates normally would (similar to the low d, low h case for (Fig. 7a). Smaller harems are less resilient to female- unbiased harvest). As in the previous scenarios, the biased harvest for the same reason they are less resilient growth of large h populations is less affected by har - to divorce and female mortality - because they cannot vest. Even though male-biased mortality causes unions lose as many females before dissolving. Female-biased to break up more frequently, it also returns (potentially harvest may reduce the average harem size, and also many) females to the f pool. This may be beneficial makes it difficult for large harems to form or reform when h is high, as a higher proportion of females is after breaking up. When h is high, a higher proportion needed to maximize the mating rate in accordance with of females is needed to maximize the mating rate, in Eq. 21. accordance with Eq. 21, but these females are depleted by harvest. Increasing the divorce rate increases the rate As expected, the secondary sex ratio s increases during at which harems dissolve. This more drastically reduces female-biased harvest, decreases during male-biased har- growth for larger harems (contours are steeper at large vest, and undergoes only minimal changes when harvest is h), because they need more females to reform. unbiased (Fig. 8, right). Populations with high d and low h 2. Unbiased harvest (Fig.  8b) yields similar qualitative experience the largest sex ratio shifts under biased harvest. trends. Again, populations with higher divorce rates This may be because the smaller growth rates of high d, low experience greater reductions in growth, and larger har- h populations are less effective in offsetting harvest-induced ems are more affected by divorce. The effect of increas- sex ratio biases. ing divorce is not as pronounced as with female-biased Figure 9 shows how the growth rates of the mating sys- harvest (contours are flatter overall), as less female tems in Table 2 vary with harvest bias and intensity. Again, harvest makes it easier for harems to reform. At low we see that high h, low d populations (large, persistent har- h, however, populations with lower d are actually more ems) have the largest growth rates of all the mating systems, impacted by harvest. Low h unions have only a few even under harvest. Low h, high d populations (small, tran- females and are more likely to dissolve from increased sient harems) have the smallest growth rates. mortality. Unions with high divorce rates are already Increasing the total harvest rate E in Eq.  18 amplifies dissolving quickly, regardless of harvest mortality. the differences between female-biased, unbiased, and Unions with low divorce rates, in contrast, break up male-biased harvest strategies. Female-biased harvest (red) much more frequently once harvest mortality occurs. decreases population growth more severely than male-biased As a result, low d, low h populations experience the harvest (blue) does, even across populations with different h largest decreases in growth. and d. This may be because there is an excess of single males 3. Male-biased harvest (Fig.  8c) reverses the effects of waiting to become harem leaders, whereas single females increased divorce rate. Focusing harvest on males is are usually in shorter supply (especially when h is large). more likely to dissolve harems by killing their male Additionally, the death of a male harem leader immediately leaders. Populations with low d experience the larg- dissolves his union; this can return many females to the sin- est reductions in growth, because male-biased harvest gles stage, allowing new, full-sized harems to reform with makes these unions dissolve more frequently than they a new male leader. The death of female harem members, 1 3 32 Population Ecology (2018) 60:21–36 in contrast, does not necessarily cause the union to dis- and Danchin 2008). Sex-biased harvest may have different solve, and may instead generate small, stunted harems with consequences for offspring survival in these mating systems. reduced productivity. Other species have additional nuances in how they respond Depending on the mating system, unbiased harvest to sex-biased harvest; African lions, for instance, commit (black) can decrease growth rates more or less than female- infanticide when male harem leaders are killed (Whitman biased harvest does. In populations with low h (small har- et al. 2004), which exacerbates the effects of male harvest ems), female-biased harvest has the most drastic impacts of on population growth. all the harvest strategies — again, perhaps, because small How populations respond to selective harvest also has harems cannot afford to lose as many females before dis- important consequences for evolution. Growing evidence solving. In populations with high h (large harems), however, suggests that evolutionary considerations are relevant to sus- unbiased harvest may be just as, if not more, detrimental to tainable long-term management (Ashley et al. 2003), and that population growth. human-induced selection is especially important for harvested species. As harvest mortalities are often more severe and selective than natural mortalities, they may drive evolution Discussion in directions that would not occur under natural conditions. Because it integrates life cycle structure, sex ratio, and The framework we present here is a tool for studying the a mating function, the approach introduced here is ideally effects of sexual reproduction, mating systems, and life his- suited to studying the evolution of sex ratio as a compo- tories on population dynamics. Our implementation of the nent of life history evolution. Sex ratio evolution has a long BMMR approach places no limitation on the complexity of and distinguished history in evolutionary demography (e.g., the life cycle, and is flexible enough to accommodate age or Darwin 1871; Fisher 1930; Trivers 1972; Charnov 1982, stage structure, diverse mating systems, and mating preferences. among many others). Many of these discussions invoke Because it is formulated as a matrix model, powerful sensitivity demographic factors, such as relative mortality of males analysis tools are available (see Shyu and Caswell 2016a, b). and females, but do so without a demographic model that There are interesting open problems to which the incorporates those factors. approach can be applied. As formulated here, the model Our approach provides such a model, and can be analyzed describes a constant environment. Time-varying (seasonal using the framework of adaptive dynamics (e.g., Geritz et al. or stochastic) models would be interesting and challenging 1998), which accounts for the frequency-dependence inher- extensions. Spatial models connecting multiple subpopula- ent in the mating process. We have applied this to the evolu- tions might have important effects on mate limitation. And, tion of sex ratio as influenced by sex-biased offspring costs because the model is formulated in terms of pair formation, (Shyu and Caswell 2016a) and multiple maternal conditions it will be interesting to see how it might apply to species (Shyu and Caswell 2016b). In these studies, the dynamics that do not form pairs, but in which sex ratio may influence of the population frequency vector  are used to evaluate population dynamics (e.g., pollen limitation in plants). the possibility for a mutant phenotype to invade a resident In our application to sex-biased harvest we found that phenotype; phenotypes that are not invasible are singular mating factors, including harem size and union duration, strategies; their stability properties lead to the identification affect not only unharvested population growth, but also the of evolutionarily stable (ESS) or convergence stable strate- responses of growth rate and sex ratio to sex-biased harvest. gies. See Shyu and Caswell (2016a, b) for details. In unharvested populations, high rates of divorce, which lead Acknowledgements This work was supported by a National Sci- to more transient unions, tend to reduce population growth, ence Foundation Graduate Research Fellowship to ES, under Grant especially when harem size is small. Sex-biased harvest 1122374. HC acknowledges support from NSF Grants DEB1145017 affects not only population sex ratios, but also long-term and DEB1257545 and support from the European Research Coun- growth rates, with effects depending on sex bias, harem size, cil under the European Union’s Seventh Framework Programme (FP7/2007-2013), ERC Advanced Grant 322989. ES acknowledges and divorce rate. These complex, and sometimes counterin- support from the Academic Programs office of the Woods Hole Ocean- tuitive, nuances would be impossible to capture without a ographic Institution. demographic two-sex model like this, motivating the use of such models in ecological studies. Open Access This article is distributed under the terms of the Crea- tive Commons Attribution 4.0 International License (http://creat iveco Our two-variable depiction of the mating system, defined mmons.or g/licenses/b y/4.0/), which permits unrestricted use, distribu- by harem size and union duration, can be extended to other tion, and reproduction in any medium, provided you give appropriate factors. Mating systems may differ in parental investment by credit to the original author(s) and the source, provide a link to the males and females. While polygynous males tend to provide Creative Commons license, and indicate if changes were made. minimal parental care, monogamous males invest on par with their female partners (Emlen and Oring 1977; Cézilly 1 3 Population Ecology (2018) 60:21–36 33 Appendix 1: Sequential processes the derivative is dt in continuous time d d‖‖ ‖‖ − dt dt The continuous-time two sex model is obtained by imagining dt ‖‖ (32) that the processes of mating, birth, and transitions take place d‖‖ 1 d = − as a periodic sequence, each over a discrete time interval, and ‖‖ dt dt ‖‖ then shrinking the time interval so that, in the limit, all three processes are operating at every moment. The resulting contin- Let  be a 1 × s vector of ones. Because  is a uous-time projection matrix  is the average of the individual non-negative vector, transition rate matrices in Eq. 3. Recall that the solution of a linear continuous-time system is ‖‖ = (33) d‖‖ d d = (t) → (t)= e (0) (26) dt dt dt Δt Thus, the matrix exponential e projects the population Then Eq. 32 can be rewritten as over a discrete period of length Δt: 1 d  d Δt = − (t +Δt)= e (t). (27) dt   dt ( ) dt Suppose that, e.g., three constant matrices  ,  , and  act 1 d =  −  (34) in sequence, over periods of length Δt . Then dt Δt Δt Δt =  −  (t) (t + 3Δt)= e e e (t) (28) s To first order, as Δt → 0, (t + 3Δt)≈( +Δt)( +Δt)( +Δt)(t) where  is a s × s identity matrix. (29) ≈ [ +Δt( +  + )](t) If the population is initialized with a population vector =  so that ‖‖ =   = 1 , then Eq. 34 can be rewritten as: Dividing by 3Δt and taking the limit as Δt → 0 gives the derivative, d =  −  (t) (35) dt (t + 3Δt)− (t) = lim dt Δt→0 3Δt as in Eq. 10. (30) = ( +  + )(t) Appendix 3: BMMR matrices = (t) for a polygynous mating system Consider a polygynous system with a maximum harem size of as in Eq. 3. Heuristically, this can be thought of as a periodic h. When h is large, it is cumbersome to write the rate matrices product operating so fast that all three processes are operat- ,  ,  in full, especially since many of their entries will be ing simultaneously. In our case, the matrix  is a function zeros. Instead, we will consider these matrices in terms of their of (t) , which changes over the interval [t, t +Δt] . However, contributions to these nine regions of the projection matrix: as Δt gets sufficiently small, the result is equivalent to evaluated at (t). ⎛ ⎞ m→m f→m u→m Technically, this construction is known as the product ⎜ ⎟ (36) m→f f→f u→f integral (Gantmacher 1959). It originated in the mathemati- ⎜ ⎟ ⎝ ⎠ m→u f→u u→u cal work of the same Vito Volterra who is recognized by ecologists for his work on population dynamics (Volterra and For the union formation matrix  , the only regions with Hostinksy 1938). nonzero contributions are: Appendix 2: Dynamics of the relative frequency vector   = 2 × 2 (37) m→m 0 − U Given the population vector of stage frequencies 0 0 = 2 × 2 (38) f→f 0 − hU = f (31) ‖‖ 1 3 34 Population Ecology (2018) 60:21–36 a particular advantage to work with a continuous time model ⎛ ⎞ because it specifies the rates of mating rather than prob- ⎜ ⎟ ⋮⋮ = h × 2 abilities; to specify probabilities would require accounting (39) m→u ⎜ 1 ⎟ 0 U for probabilities of events that did not happen. ⎝ ⎠ The mating process, where adult males and females pair into reproducing unions, is described by the union forma- tion matrix  . Mating functions in  give the rates of union ⎛ ⎞ ⎜ ⎟ formation as functions of the relative frequencies of males ⋮⋮ = h × 2 (40) f→u ⎜ ⎟ and females available to mate, and are thus functions of the 0 U ⎝ ⎠ stage frequency vector . Mating preferences in the mating functions describe the For the birth matrix  , the only regions with nonzero con- probabilities of favoring partners of certain conditions. The tributions are: female preference distribution g (i) gives the proportion ks 2ks … khs of Condition j females that mate with Condition i males. 1 1 1 = 2 × h (41) u→m Similarly, the male preference distribution h (j) gives the 00 … 0 proportion of Condition i males that mate with Condition j females. Summing these distributions over all male and k(1 − s ) 2k(1 − s )… kh(1 − s ) 1 1 1 = 2 × h u→f female conditions respectively yields a total probability of 1: 00 … 0 (42) g (i)= 1 ∀ j (47) For the transition matrix  , the only regions with nonzero contributions are: h (j)= 1 ∀ i −( +  ) 0 (48) m1 m = 2 × 2 (43) m→m m m2 Examples of mating preference distributions include: −( +  ) 0 f 1 f = 2 × 2 (44) f→f 1. Fully assortative mating, where individuals only mate f f 2 with partners in the same condition: 0 0 … 0 g (i)= 1 if i = j, 0 else = 2 × h (45) u→m j + d 0 … 0 f 2 (49) h (j)= 1 if i = j, 0 else 00 … 0 2. Random mating, where individuals pick partners based = 2 × h (46) u→f + d 2 + d … h + d m2 m2 m2 on their relative abundances in the population: The h × h submatrix  is also nonzero. It contains entries u→u g (i)= of −( +  + d) all along its diagonal, and entries of m2 f 2 + d all along its first superdiagonal. f 2 (50) As an example, the  ,  ,  matrices in the case where h (j)= h = 3 are provided in Eqs. 23, 24, and 25 respectively. 3. Biased mating, where individuals prefer partners of cer- Appendix 4: Multiple mating stages tain conditions. An attractiveness or competitiveness and mating preferences factor c weighs the abundance of each partner condi- tion, e.g.,: Suppose that there are q adult male stages and r adult c m female stages. There are qr types of unions possible. The i i g (i)= rate of production of (i, j) matings depends on males of stage c m i i i i and females of stage j, and on preferences. In this case, it is (51) c f j j h (j)= ∑ . c f j j j The material in Appendix  4 is extracted and modified from Shyu and Caswell 2016b under a Creative Commons Attribution license. 1 3 Population Ecology (2018) 60:21–36 35 Cézilly F, Danchin É (2008) Mating systems and parental care. In: Partners with larger c are more preferable mates. If all Danchin É, Girladeau LA, Cézilly F (eds) Behavioural ecology. c are equal, Eq. 51 reduces to the random mating case Oxford University Press, Oxford, pp 429–465 Eq. 50. If c = 0 , individuals of stage i do not mate. Charnov EL (1982) The theory of sex allocation. Princeton Univer- sity Press, Princeton Darwin C (1871) The descent of man and selection in relation to sex. The total mating function M () gives total unions u ij ij John Murray, London (Condition i males mated with Condition j females) formed Emlen ST, Oring LW (1977) Ecology, sexual selection, and the evo- per time. The most general and flexible mating functions lution of mating systems. Science 197:215–223 are based on generalized weighted means (Hölder means). Festa-Bianchet M (2003) Exploitative wildlife management as a selective pressure for life-history evolution of large mammals. These have the general form: In: Festa-Blanchet M, Apollonio M (eds) Animal behavior and wildlife conservation. Island Press, Washington, DC, pp 1∕ a a M ()= b[f g (i)] +(1 − b)[m h (j)] (52) ij j j i i 191–210 Fisher RA (1930) The genetical theory of natural selection. Oxford where 0 ≤ b ≤ 1 and a < 0 (Hadeler 1989; Martcheva and University Press, Oxford Fredrickson A (1971) A mathematical theory of age structure in Milner 2001; Caswell 2001). Note that M () is calculated ij sexual populations: random mating and monogamous marriage only over individuals that are available to mate (i.e., adult models. Math Biosci 20:117–143 single male stages m and adult single female stages f ). As a i j Gantmacher FR (1959) The theory of matrices. Chelsea, New York result, the mating function does not depend on the males and Geritz SAH, Mesze G, Metz JA (1998) Evolutionarily singular strate- gies and the adaptive growth and branching of the evolutionary females in non-mating stages, such as immature juveniles or tree. Evol Ecol 12:35–57 adults already in unions. Ginsberg JS, Milner-Gulland EJ (1994) Sex-based harvesting and The harmonic mean mating function is one of the most population dynamics in ungulates: implications for conservation widely used, because it satisfies the biological criteria for and sustainable use. Conserv Biol 18:157–166 Hadeler KP (1989) Pair formation in age-structured populations. two sex models and captures the qualitative properties of Acta Appl Math 14:91–102 a wide range of Holder means (Caswell and Weeks 1986; Hadeler KP (1993) Pair formation with maturation period. J Math Iannelli et al. 2005). Hence, we use a harmonic mean mat- Biol 32:1–15 ing function where a =−1, b = 1∕2 , so that: Hadeler KP, Waldstätter R, Wörz-Busekros A (1988) Models for pair formation in bisexual populations. J Math Biol 26:635–649 Hardy ICW (2002) Sex ratios: concepts and methods. Cambridge 2m h (j)f g (i) i i j j M ()= . University Press, Cambridge (53) ij m h (j)+ f g (i) i i j j Hoppensteadt F (1975) Mathematical theories of populations: demo- graphics, genetics, and epidemics. Society for Industrial and Applied Mathematics, Philadelphia The corresponding male and female per capita mating func- Iannelli M, Martcheva M, Milner F (2005) Gender-structured tions are: population modeling: mathematical methods, numerics, and simulations. Society for Industrial and Applied Mathematics, M () ij Philadelphia U ()= m,ij m Jenouvrier S, Caswell H, Barbraud C, Weimerskirch H (2010) Mat- (54) ing behavior, population growth, and the operational sex ratio: a M () ij periodic two-sex model approach. Am Nat 175:739–752 U ()= . f ,ij f Karmel PH (1974) The relations between male and female reproduction rates. Popul Stud 1:249–274 Kendall DG (1949) Stochastic processes and population growth. J R Stat Soc B 11:230–282 References Keyfitz N (1972) The mathematics of sex and marriage, vol 4. In: Proceedings of the sixth Berkeley symposium on mathematical Allendorf FW, Hard JJ (2009) Human-induced evolution caused by statistics and probability, pp 89–108 unnatural selection through harvest of wild animals. Proc Natl Kuczynski RR (1932) Fertility and reproduction. Falcon Press, New Acad Sci USA 106:9987–9994 York Ashley MV, Willson MF, Pergams ORW, O’Dowd DJ, Gende SM, Martcheva M (1999) Exponential growth in age-structured two-sex Brown JS (2003) Evolutionarily-enlightened management. Biol populations. Math Biosci 157:1–22 Conserv 111:115–123 Martcheva M, Milner F (2001) The mathematics of sex and marriage, Bessa-Gomes C, Legendre S, Clobert J (2010) Discrete two-sex mod- revisited. Math Popul Stud 9:123–141 els of population dynamics: on modeling the mating function. McFarland DD (1972) Comparison of alternative marriage models. Acta Oecol 36:439–445 In: Grenville TNE (ed) Population dynamics. Academic Press, Caswell H (2001) Matrix population models. Sinauer Associates, New York, pp 89–106 Sunderland Miller TEX, Inouye BD (2011) Confronting two-sex demographic Caswell H, Shyu E (2012) Sensitivity analysis of periodic matrix models with data. Ecol 92:2141–2151 population models. Theor Popul Biol 82:329–339 Miller TEX, Shaw AK, Inouye BD, Neubert MG (2011) Sex-biased dis- Caswell H, Weeks DE (1986) Two-sex models: chaos, extinction, persal and the speed of two-sex invasions. Am Nat 177:549–561 and other dynamic consequences of sex. Am Nat 128:707–735 1 3 36 Population Ecology (2018) 60:21–36 Milner JM, Nilsen EB, Andreaseen HP (2007) Demographic side Shyu E, Caswell H (2016a) A demographic model for sex ratio evo- effects of hunting ungulates and carnivores. Conserv Biol lution and the effects of sex-biased offspring costs. Ecol Evol 21:36–47 6:1470–1492 Pollak RA (1986) A reformulation of the two-sex problem. Demog- Shyu E, Caswell H (2016b) Frequency-dependent two-sex models: a raphy 23:247–259 new approach to sex ratio evolution with multiple maternal condi- Pollak RA (1987) The two-sex problem with persistent unions: a gen- tions. Ecol Evol 6:6855–6879 eralization of the birth matrix-mating rule model. Theor Popul Trivers RL (1972) Parental investment and sexual selection. In: Camp- Biol 32:176–187 bell B (ed) Sexual selection and the descent of man, 1871–1971. Pollak RA (1990) Two-sex demographic models. J Polit Econ Aldine, Chicago, pp 136–179 98:399–420 Volterra V, Hostinksy B (1938) Opérations infinitésimales linéaires. Pollard AH (1948) The measurement of reproductivity. J Insti Actuar Gauthier-Villars, Paris (in French) 74:288–337 Whitman K, Starfield AM, Quadling HS, Packer C (2004) Sustainable Pollard JH (1974) Mathematical models for the growth of human popu- trophy hunting of African lions. Nature 428:175–177 lations. Cambridge University Press, Cambridge Yellin J, Samuelson PA (1974) A dynamical model for human popula- Rankin DJ, Kokko H (2007) Do males matter? The role of males in tion. Proc Natl Acad Sci USA 71:2813–2817 population dynamics. Oikos 116:335–348 Sundelöf A, Åberg P (2006) Birth functions in stage structured two-sex models. Ecol Model 193:787–795 1 3 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Population Ecology Springer Journals

Mating, births, and transitions: a flexible two-sex matrix model for evolutionary demography

Free
16 pages

Loading next page...
 
/lp/springer_journal/mating-births-and-transitions-a-flexible-two-sex-matrix-model-for-kV0Jt50eo9
Publisher
Springer Journals
Copyright
Copyright © 2018 by The Author(s)
Subject
Life Sciences; Ecology; Zoology; Plant Sciences; Evolutionary Biology; Behavioral Sciences; Forestry
ISSN
1438-3896
eISSN
1438-390X
D.O.I.
10.1007/s10144-018-0615-8
Publisher site
See Article on Publisher Site

Abstract

Models of sexually-reproducing populations that consider only a single sex cannot capture the effects of sex-specific demo- graphic differences and mate availability. We present a new framework for two-sex demographic models that implements and extends the birth-matrix mating-rule approach of Pollak. The model is a continuous-time matrix model that explicitly includes the processes of mating (which is nonlinear but homogeneous), offspring production, and demographic transitions and survival. The resulting nonlinear model converges to exponential growth with an equilibrium population composition. The model can incorporate age- or stage-structured life histories and flexible mating functions. As an example, we apply the model to analyze the effects of mating strategies (polygamy or monogamy, and mated unions composed of males and females, of variable duration) on the response to sex-biased harvesting. The combination of demographic complexity with the interaction of the sexes can have major population dynamic effects and can change the outcome of evolution on sex- related characters. Keywords Birth matrix-mating rule · BMMR · Demography · Matrix population models · Sex-biased harvest · Two-sex models Introduction mortalities (Kuczynski 1932; Pollak 1990; Jenouvrier et al. 2010), developmental schedules (Caswell 2001), behavioral Models of sexually-reproducing populations that consider interactions (Rankin and Kokko 2007), dispersal patterns a only single sex (typically females) implicitly assume that (Miller et al. 2011), and selective harvest pressures (Ginsberg both sexes have identical vital rates and that the availability and Milner-Gulland 1994). Additionally, ecological, environ- of the neglected sex (typically males) does not affect fertility mental, and evolutionary changes may alter the most limiting (Pollard 1974; Caswell 2001; Iannelli et al. 2005). In reality, sex over time (Hardy 2002; Miller and Inouye 2011). One- both these assumptions are frequently violated. Males and sex models miss, and cannot explore, important components females often differ significantly in terms of fertilities and of population dynamics in all these cases. Here, we introduce a flexible framework that can explore many of these factors. In response to growing concerns over discrepancies in male and female reproductive rates (Karmel 1947), early Electronic supplementary material The online version of this article (https ://doi.org/10.1007/s1014 4-018-0615-8) contains dynamical models with sex structure were introduced in the supplementary material, which is available to authorized users. late 1940s. Pollard (1948) used coupled Lotka integral equa- tions, considering female births to males and male births to * Hal Caswell females in order to reconcile the growth rates of both sexes. h.caswell@uva.nl Kendall (1949) introduced a system of ordinary differential Biology Department MS-34, Woods Hole Oceanographic equations for males and females (and later, married cou- Institution, Woods Hole, MA 02543, USA ples), and was the first to incorporate nonlinear interactions Institute for Biodiversity and Ecosystem Dynamics, between the sexes via a mating term. Subsequent models University of Amsterdam, PO Box 94248, considered other nonlinear mating functions (Pollard 1974; 1090 GE Amsterdam, The Netherlands Vol.:(0123456789) 1 3 22 Population Ecology (2018) 60:21–36 Yellin and Samuelson 1974) and couple dissolution through Model development death and divorce (Hadeler et al. 1988). Extensions to age-structured populations were made by For reasons that will become apparent, we have incorporated Fredrickson (1971) and Hoppensteadt (1975), who allowed sex structure, stage structure, and life cycle processes into birth and death rates, as well as couple formation and divorce a continuous-time, rather than a discrete-time matrix popu- rates, to depend on age and sex. Hadeler later included the lation model. The mating, birth, and transition processes age (i.e., duration) of married pairs (Hadeler et al. 1988) from the BMMR framework are described by separate rate and maturation delays (Hadeler 1993). Age-structured mat- matrices. The mating process introduces nonlinearity into ing functions have similarly been proposed (Martcheva and the model because reproduction depends on the relative pro- Milner 2001). Many of these models use continuous-time portions of males and females, of appropriate stages, in the equations that incorporate age structure through coupled population. This dependence can be flexibly modeled by McKendrick–von Foerster partial differential equations (Fre- generalized weighted mean mating functions, which satisfy drickson 1971; Keyfitz 1972; Hoppensteadt 1975; Hadeler the biological criteria of sexual reproduction (Iannelli et al. 1989, 1993), though discrete-time two sex matrix models 2005). The resulting BMMR matrix models are nonlinear have also been developed (Caswell and Weeks 1986). but homogeneous (i.e., frequency-dependent). Models of A powerful conceptual approach to two-sex models was this form generally converge to an exponential growth rate proposed by Pollak (1986, 1987, 1990) and called the birth and a stable stage frequency distribution (Martcheva 1999). matrix-mating rule (BMMR) model. It proposes that mating, births, and life cycle transition processes repeat periodically, Incorporating sex and stage structure one after the other. It contains three main components: The model classifies individuals into stages based on age, 1. A mating rule function that gives the number of matings developmental state, sex, reproductive status, or other vari- u between males of age i and females of age j. ables of interest. Stage densities are projected forward in ij 2. A birth matrix whose entries b are the expected number time by a projection matrix that contains the demographic ij of male and female offspring produced by a mating of a rates characterizing survival, reproduction, and transitions male of age i and a female of age j. between stages. The properties of this projection matrix 3. Sex-specific mortality schedules, which project surviv - provide information about the population as a whole, ing individuals to the next age class, or, in our generali- thereby linking individual-level life cycle information (i.e., zation, include other stage-specific life cycle transitions. the stage-specific vital rates in the matrix entries) to popu- lation-level properties important for ecology and evolution BMMR is a useful approach for describing two-sex popu- (e.g., growth rates or stage distributions). lations because it can specify age (and, more generally, A population with s stages is described by a s × 1 popula- stage) structure over all parts of the life cycle. This struc- tion vector (t) , the entries of which are the densities of each ture, in turn, can have significant effects on two-sex popula- stage at time t. In a two-sex population, (t) would contain tion dynamics (e.g., Sundelöf and Åberg 2006, where the male stages, female stages, and mated stages (unions) that addition of size-specific birth functions affects growth rate could include married couples or breeding harems. and reproductive output) and recommended management For example, the population vector for a two-sex popula- strategies (e.g., Ginsberg and Milner-Gulland 1994, where tion with mating adults and nonmating juveniles could have incorporating age-specific fecundity changes the outcomes the form: of sex-biased harvest). Conceptually appealing as it is, the BMMR model has yet ⎛ ⎞ juvenile males ⎜ ⎟ to be fully incorporated into the framework of stage-classie fi d adult males ⎜ ⎟ matrix population models. Our goal here is to do so, providing (t)= juvenile females (1) ⎜ ⎟ a general model that can be adapted to a wide range of life ⎜ ⎟ adult females ⎜ ⎟ cycles and mating strategies. We do so by a novel extension adult unions ⎝ ⎠ of periodic matrix models to continuous time, based on transi- tion rate matrices. In this paper, we focus on ecological impli- cations of the model, but the approach can also be applied to The dynamics of the population vector are given by a system of ordinary differential equations study the evolution of sex-related traits using methods from adaptive dynamics (Shyu and Caswell 2016a, b). = [] (t) (2) dt 1 3 Population Ecology (2018) 60:21–36 23 where the entries of  are either transition rates or rates of offspring production, and we have indicated that they may depend on the population vector. Incorporating the BMMR processes The BMMR framework incorporates mating, birth, and transition processes. We describe each of these processes by a separate matrix: 1. The mating (union formation) process, where adult Fig. 1 Mating functions from the generalized weighted mean family males and females organize into reproductive unions, is (Eq.  6) with a check to indicate which of the biologically desirable described by the matrix . mating function criteria they satisfy 2. The birth process, where unions produce new offspring, is described by the matrix . 3. The transition process, which includes other life cycle M() events like mortality, maturation, or divorce, is described U ()= (4) by the matrix . M() Other life cycle processes can be included with additional U ()= . (5) matrices. A discrete-time model would incorporate these processes The total population mating rate is M = U m = U f . m f as a periodic matrix product (Caswell and Shyu 2012). For Many commonly used mating functions are generalized example, the product  would describe mating, followed weighted means (Hölder means) of the form by the production of offspring from the matings, followed 1∕ by survival and transitions of the resulting individuals. In M()= f +(1 − )m (6) continuous-time, we conceptualize the processes as occur- where  and  are constants; 0 ≤  ≤ 1 , 𝛼< 0 (Hadeler 1989; ring simultaneously. It can be shown (Appendix 1) that the Bessa-Gomes et al. 2010; Iannelli et al. 2005). Figure  1 projection matrix  in the continuous-time matrix model shows several generalized weighted mean mating functions (Eq. 2) is the average of the transition rate matrices, e.g.,: and biologically desirable criteria that they satisfy (McFar- d 1 land 1972; Pollard 1974; Yellin and Samuelson 1974). = ( +  + )(t) dt 3 (3) If multiple male and female stages interbreed to form = (t) different types of unions, stage-specific mating preferences can also be integrated into this mating function (Shyu and Caswell 2016b; see  Appendix 4). Modeling the mating process It is difficult to distinguish between mating functions in populations where the sex ratio does not vary signifi- The mating process, as described by the union formation matrix , depends on the relative numbers of males and females in the cantly (Keyfitz 1972). However, recent empirical studies (Miller and Inouye 2011) support the harmonic mean as a population, not all of which may be mature enough or available for breeding (Pollard 1974; Iannelli et al. 2005). As a result,  mating function. Because the harmonic mean satisfies rea- sonable biological criteria for mating functions (Caswell depends on the population’s sex and stage composition, making a function of the population vector (t). and Weeks 1986; Iannelli et al. 2005), and captures the qualitative properties of other generalized means, we will The total mating function M() gives the rate of union forma- tion (total number of unions formed per unit time in a population use a harmonic mean mating function for our analyses. of composition  ). This embodies the mating rule in the BMMR framework. Here, “unions” refers to any mated, reproducing Frequency‑dependent dynamics units in the population, including both one-to-one male-female pairs and harems with multiple individuals of the one sex. The mating process often depends on the relative frequen- cies, rather than absolute abundances, of males and females. We convert the total mating function into the per capita mating rates U () and U () (the average mating rates per As a result, although the mating function is nonlinear, it is m f homogeneous of degree 1 in  . That is: available males m or females f respectively), 1 3 24 Population Ecology (2018) 60:21–36 M(c)= cM() (7) for any positive constant c. As a result, the per capita mating functions (Eqs. 4 and 5) are homogeneous of degree 0 in  , so that: U (c)= U () m m (8) U (c)= U () f f If all entries in the projection matrix  in Eq.  3 are also Fig. 2 Life cycle diagram for a 5-stage population with juvenile males m and juvenile females f , adult males m and adult females homogeneous of degree 0, the system is said to be  fre- 1 1 2 f , and reproducing unions u. The functions and parameters shown quency-dependent. This means that  can be written as a here, as described in Table 1, appear in the union formation matrix function of the population frequency vector: (Eq. 14) (red), birth matrix  (Eq. 15) (green), or transition matrix (Eq. 16) (blue). From Shyu and Caswell (2016a) (9) ‖‖ where ‖‖ is the 1-norm of . 2m f 2 2 Frequency-dependent models of this type usually con- M()= m + f 2 2 verge asymptotically to an equilibrium population structure 2f ̂ (the stable stage distribution) in which all stage frequen- U ()= (12) m + f 2 2 cies are constant (e.g., see Yellin and Samuelson 1974; 2m Hadeler 1989; Martcheva 1999). The population then grows 2 U ()= m + f or decays exponentially at a long-term growth rate given by 2 2 the dominant eigenvalue  of [ ̂ ]. To find the equilibrium stage distribution  and popula- Again, we consider the life cycle in terms of mating, birth, tion growth rate  , it is sufficient to consider the dynamics and transition processes, which are described by matrices  , of (t) . It can be shown (Appendix 2) that: , and  respectively. =  −  []  (10) dt 1. The union formation matrix  contains the per capita mating functions from Eq. 12. where  is a s × s identity matrix and  is a 1 × s vector of ones. One can integrate Eq. 10 with a numerical differential ⎛000 0 0⎞ equation solver until population frequencies converge to  ̂ . ⎜ ⎟ 0 − U () 0 00 ̂ m Then  is the dominant eigenvalue of [] . The dominant ⎜ ⎟ 000 0 0 ()= (13) right eigenvector  of [ ̂ ] is proportional to the stable stage ⎜ ⎟ ⎜000 − U () 0⎟ distribution . ⎜ 1 1 ⎟ 0 U () 0 U () 0 ⎝ m f ⎠ 2 2 A 5‑stage BMMR matrix model ⎛0 00 00⎞ M() ⎜ ⎟ 0 − 0 00 We now present an example of a BMMR matrix model with ⎜ ⎟ five stages: juvenile males m and juvenile females f , adult ⎜ ⎟ 0 00 00 1 1 (14) ⎜ M() ⎟ males m and adult females f , and reproducing unions u that 0 00 − 0 2 2 ⎜ ⎟ consist of one adult male and one adult female. Single adult M() M() ⎜ ⎟ 0 0 0 ⎝ ⎠ 2m 2f males and females interact to form unions, which then produce 2 2 new juvenile offspring (Fig.  2). A summary of the variables, parameters, and matrices in this model is provided in Table 1. Note that U and U must halved in the last row of  to m f As in Eq. 1, we write the population vector as avoid double counting the unions formed from both male and female partners. (t)= m m f f u (11) 1 2 1 2 2. The birth matrix  contains k, the average reproductive rate of a union, and the primary sex ratio s , the propor- Using the harmonic mean mating function, the total and per tion of offspring that are male. capita mating functions are 1 3 Population Ecology (2018) 60:21–36 25 a b 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 01 23 45 67 8 01 23 45 67 8 time time Fig. 3 Dynamics of the 5-stage BMMR model with monogamous (11), where dynamics are given by (3). b) Convergence of the pop- unions and no harvest. The population consists of juvenile males m ulation frequency vector  (9), where dynamics are given by (17). and juvenile females f , adult males m and adult females f , and Parameters are fixed at  =  = 0.5,  =  = 0.1,  =  = 1 2 2 m1 f 1 m2 f 2 m f reproducing unions u. a) Growth of the population density vector  0.5, s = 0.5, k = 20, d = 0.1, h = 1 , E = 0. Table 1 A summary of the Matrices and vectors variables, parameters, matrices, Projection matrix and population properties in the 5-stage BMMR matrix model Birth matrix (15) Transition matrix (16) Union matrix (14) Population density vector (11) Population frequency vector (9) ̂ or  Equilibrium stage structure Population properties Long-term population growth rate, dominant eigenvalue of [ ̂ ] s Primary sex ratio (proportion of offspring that are born male) s Secondary sex ratio (proportion of adults that are male) f , f Juvenile, adult female density 1 2 m ,m Juvenile, adult male density 1 2 u Union (mated pair) density Life cycle parameters ,  Male, female maturation rates m f d Divorce rate (rate at which a male-female pair bond breaks) h Average harem size k Union reproductive rate ,  Juvenile, adult female mortality rates f 1 f 2 ,  Juvenile, adult male mortality rates m1 m2 M Total mating rate (total unions formed per time) U , U Per capita mating rates (4) and (5) m f Harvest parameters E Total harvest rate in (18) s Harvested sex ratio (proportion of harvest that targets males) in Eq. 18 1 3 stage densities stage frequencies 26 Population Ecology (2018) 60:21–36 Mating systems and the effects of sex‑biased ⎛ ⎞ 00 00 ks harvest ⎜ ⎟ 00 00 0 ⎜ ⎟ = 00 00 k(1 − s ) (15) As an example of the use of the two-sex model, we consider ⎜ 1 ⎟ ⎜ ⎟ 00 00 0 the effectes of sex-biased harvesting. In sport or trophy hunt- ⎜ ⎟ 00 00 0 ing, harvest is often male-biased and significantly exceeds ⎝ ⎠ natural mortality (Festa-Bianchet 2003; Milner et al. 2007). Age or size bias is also common, as larger or older males 3. The transition matrix  contains the male mortality rates with well-developed adult characteristics (e.g., large antlers (  for juveniles,  for adults) and female mortal- m1 m2 or horns) are usually the most desirable targets. This unnat- ity rates (  for juveniles,  for adults), the rates of f 1 f 2 ural selection may alter population structure, reproductive maturation from juveniles to adults (  for males,  for m f strategies, body morphology, and developmental timing (e.g., females), and the divorce rate d (rate at which the male- Ashley et al. 2003; Festa-Bianchet 2003; Allendorf and Hard female pair bond breaks). Note that unions may also 2009). The population response depends on multiple demo- dissolve due to partner death, and that union dissolution graphic factors, including the mating system. through both death and divorce returns surviving males The mating system determines how males and females and females to the single adult stages. organize for breeding and is thus a key component of two- sex population structure (Emlen and Oring 1977). Some species form monogamous, one-to-one pair bonds between ⎛−( +  ) 0 00 0 ⎞ m1 m ⎜ ⎟ males and females. Other species have polygynous mating −  00  + d ⎜ m m2 f 2 ⎟ systems in which one male mates with multiple females ⎜ ⎟ = 00 −( +  ) 00 f 1 f ⎜ ⎟ (e.g., lions, deer, seals), or, more rarely, polyandrous systems ⎜ ⎟ 00  −   + d f f 2 m2 ⎜ ⎟ where one female mates with multiple males (e.g., jacanas, ⎝ 0 0 00 −( +  + d)⎠ m2 f 2 pipefish). Unions formed by mating may be transient and (16) limited to a single breeding episode (e.g., lek systems) or The advantages of the continuous-time formulation are may persist over multiple breeding seasons (e.g., lion har- apparent at this point. The rates of maturation, mortality, ems) or even until partner death (e.g., albatrosses and other and divorce in Eq. 16 combine in a simple additive fashion. species with high mate fidelity) (Cézilly and Danchin 2008). Transition probabilities would not combine additively; they These factors motivate the use of a demographic two-sex would involve sums of products of conditional probabilities, model in analyzing harvest effects. To this end, we will use over all the possible events. The number of these products our BMMR matrix framework to explore the effect of mating increases dramatically when more stages, and hence more systems on the response to sex-biased harvest. A range of mating combinations, are included. See Shyu and Caswell mating systems will be approximated by varying two model (2016b) and Appendix 4 for incorporation of multiple mat- parameters, d and h (Table 2). ing stages in the continuous-time model. As per Eq. 3, the average of these three matrices is the – The divorce rate d. This is a measure of union transience. continuous-time projection matrix [] . The corresponding Unions with higher values of d are more likely to dissolve equation for the proportional structure (Eq. 10) is thus: after a given mating, while unions with lower values of d are more likely to persist over multiple breeding interactions. =  −  ( +  + []) (17) – The harem size h. This is a measure of polygamy. Unions dt 3 with h = 1 are monogamous and consist only of one-to- where  is given by Eq. 14,  is given by Eq. 15, and  is one male-female pair bonds, while unions with h > 1 are given by Eq. 16. polygamous groups of size h + 1 . As polyandrous mating As shown in Fig. 3, the population vector  ultimately systems are relatively rare (Cézilly and Danchin 2008), grows exponentially, while the frequency vector  ultimately we will consider only the polygynous form of polygamy, converges to the constant distribution of stages  ̂ . To deter- where one male mates with h females. mine the equilibrium stage distribution  ̂ , we integrated Eq. 17 with the Matlab ODE45 differential equation solver Harvest strategies are characterized by the overall har- until population frequencies converged (e.g., until vector vest rate E and the harvested sex ratio s (proportion of entries do not change significantly over consecutive inte- the total harvest rate that targets males). Assuming that gration intervals). only adults are harvested, the adult mortality rates are The Matlab code used for all the following examples is modified by harvest as follows: provided in the Electronic Supplementary Material. 1 3 Population Ecology (2018) 60:21–36 27 Table 2 Mating systems Low d (persistent unions) High d (transient unions) corresponding to different values of the divorce rate d and h = 1 (monogamy) Persistent pair bonds, high mate fidelity (e.g., Serial pair bonds (e.g., harem size h albatross) humans, Emperor penguins) h > 1 (polygyny) Persistent harems (e.g., lion prides) Leks, scramble competi- tion (e.g., grouse, cod, horseshoe crabs) c change in λ from harvest d change in s from harvest -0.125 0.4 0.4 a 0.9 0.9 -0.13 0.3 0.8 0.8 0.3 -0.135 0.2 0.7 0.7 -0.14 0.1 0.2 0.6 0.6 -0.145 0.5 0.5 -0.15 0.4 0.4 0.3 b -0.1 -0.155 0.3 0.3 0.25 -0.2 -0.16 0.2 0.2 0.2 -0.3 0.1 0.1 -0.165 0.15 0 0 -0.4 01 23 45 01 23 45 01 23 45 divorce rate d divorce rate d divorce rate d Fig. 4 Population growth rates, structure, and responses to adult har- the divorce rate d. c) The change in  and d) the change in second- vest in the monogamous ( h = 1 ) model. a) Proportion of mated adults ary sex ratio s , when a proportion s of harvest targets males. With- 2 h (adults in stage u, rather than m or f ) for an unharvested population, out harvest, s = 0.5 for all values of d. Other parameters are fixed at 2 2 2 and the b) corresponding population growth rates  , as functions of  =  = 0.5,  =  = 0.1,  =  = 0.5, s = 0.5, k = 20, E = 1 m1 f 1 m2 f 2 m f 1 As shown in Fig.  4a, the proportion of adults in the →  + Es m2 m2 h reproductive union stage (mated adults) decreases as d (18) →  + E(1 − s ) increases. The unharvested population growth rate  simi- f 2 f 2 h larly decreases (Fig. 4b), because populations with more To determine how various mating systems, as characterized transient couples (fewer mated adults) cannot produce off- by different values of h and d, respond to sex-biased harvest, spring as rapidly as populations with more persistent cou- we will examine harvest effects on the long-term population ples (more mated adults). Because males and females have growth rate  and the secondary sex ratio s (proportion of the same baseline vital rates, the secondary sex ratio s 2 2 all adults that are male). We will assume that males and remains unbiased ( s = 0.5 , not shown) for all values of d. females have the same baseline vital rates, and that the pri- Figure 4c shows how increasingly sex-biased harvest mary sex ratio (proportion of males at birth) is 0.5. Thus, strategies affect population growth. Both biased and unbi- the main sex-specic fi die ff rences we consider are sex-biased ased harvest strategies most strongly reduce growth in harvest pressures and male versus females roles within the populations with lower divorce rates, as adult mortality polygynous mating systems. will also disrupt pair bonds. The greatest decreases in occur when harvest is strongly sex-biased, i.e., s is close Monogamy ( h = 1) to 0 (only females are harvested) or 1 (only males are har- vested). This suggests that monogamous populations with Consider a monogamous two-sex population with juve- high fidelity pair bonds will be the most impacted by sex- niles and adults. The mating process forms unions that are biased harvest, and that concentrating harvest on a single one-to-one pair bonds of adult males and females. This sex will more greatly reduce population growth. scenario uses the rate matrices  ,  , and  as given by Figure 4d shows how the same harvest strategies decrease Eqs. 14, 15, and 16 respectively, and the mating functions the secondary sex ratio s from equality. Unsurprisingly, male- in Eq. 12. We vary the divorce rate d to explore the effects biased strategies reduce s , female-biased strategies increase of transient vs. persistent pair bonds. s , and unbiased strategies ( s = 0.5 ) have relatively minimal 2 h 1 3 % adults mated λ, no harvest harvested sex ratio s harvested sex ratio s h 28 Population Ecology (2018) 60:21–36 effect. Populations with greater divorce rates experience larger reductions in s , possibly due to their lower growth rates (Fig. 4a) being unable to replenish harvested individuals as rapidly. Polygyny ( h > 1) Now consider a polygynous population that forms unions con- sisting of one male with a harem of females. Because the death or departure of a single female changes the harem’s size and reproductive rate, we must now account for multiple union (harem) stages. The stage u represents polygynous unions consisting of one male and a harem of i females. When h is the maximum harem size, the population vector contains h union stages, which Fig. 5 Reproduction and three possible transitions for u , a union with range from a union with a harem of size 1 ( u , equivalent to 1 harem size i a monogamous couple) to a union with a harem of size h (u , the largest possible union): (t)= m m f f u u … u (19) 1 2 1 2 1 2 h Assume that when males and females mate, they form the largest possible union u . Their union formation rate is still given by the harmonic mean mating function Eq. 12, but with the number of single females now replaced by the num- ber of prospective harems: Fig. 6 Stages in a population with maximum harem size h, which include juveni le males m and juvenile females f , adult males m and 1 1 2 2m adult females f , and reproducing unions u (one male with a harem 2 i M()= of i females). Adults form harems of size h when mating, and these m + harems can shrink in size, but not grow, over time 2f (20) U ()= hm + f 2 2 2m 2. A female harem member dies (with mortality rate  ). f 2 U ()= hm + f This shrinks the union from u to u . For the union 2 2 i i−1 u , which has only one female, u = u corresponds to 1 i−1 0 the single adult male stage m (i.e., the death of a wife returns her husband to the pool of unmated singles). Note that the total union formation rate M() is maximized 3. A female harem member departs from the union, with when the sex ratio of single adults is divorce rate d. We assume that only females leave, which 2 seems biologically plausible. Her departure returns one (21) m + f female to f and shrinks the union from u to u . For 2 2 1 + h 2 i i−1 the union u , divorce dissolves the union and returns one male to m and one female to f . Thus, as the harem size h increases, a higher proportion of 2 2 single females is needed to maximize the mating rate. As a result, a union may shrink (but not grow) in size due to If an individual female has a reproductive rate of k, the departure or death of its members (Fig. 6). After a union a harem with i females has a total reproductive rate of shrinks to the smallest possible size h = 1 , or if the male harem ik; larger harems are thus more productive. Each union leader dies, the union dissolves and its members return to the u , regardless of harem size, can change in three possible stages for unmated adults. ways (Fig. 5): Appendix 3 shows how to write the rate matrices  ,  , and for a polygynous system with maximum harem size h. The 1. The male harem leader dies (with mortality rate  ). m2 population vector and matrices for the case where h = 3 are This returns i adult females to the stage f . as follows: 1 3 Population Ecology (2018) 60:21–36 29 a λ before harvest b s before harvest 0.497 14 14 0.4 0.496 0.38 12 12 0.495 0.36 10 10 0.34 0.494 0.32 8 8 0.493 0.3 6 6 0.492 0.28 0.26 0.491 4 4 0.24 0.49 2 2 01 23 45 01 23 45 divorce rate d divorce rate d Fig. 7 Unharvested population dynamics in the polygynous ( h > 1 ) union model, as functions of the divorce rate d and harem size h. a) Popula- tion growth rate  . b) Secondary sex ratio s (proportion of males in the adult population). Other parameters are the same as in Fig. 4 Figure 7 shows how the unharvested population rate of (t)= m m f f u u u (22) 1 2 1 2 1 2 3 growth and secondary sex ratio vary across different values of d and h. As in the monogamous case, lower divorce rates ⎛ ⎞ 000 0 0 0 0 result in more mated reproducing adults and, thus, higher ⎜ ⎟ 0 − U 0 0 0 0 0 population growth. Larger harems lead to more rapid popu- ⎜ ⎟ 000 0 0 0 0 ⎜ ⎟ lation growth, possibly because of their higher total repro- ⎜ ⎟ 000 − 3U 00 0 f (23) ductive rates. Unions with a high maximum harem size are ⎜ ⎟ 000 0 0 0 0 more resilient to divorce and female mortality, because they ⎜ ⎟ ⎜000 0 0 0 0⎟ can lose more females before shrinking to a u and then 1 1 ⎜ ⎟ 0 U 0 U 00 0 dissolving. m f ⎝ ⎠ 2 2 Even without sex-biased harvest, the secondary sex ratio is slightly female-biased ( s ≈ 0.494 ), but varies only a few tenths of a percentage point across a wide range of h and d. 00 00 ks 2ks 3ks ⎛ 1 1 1 ⎞ Populations with high h and low d (large, persistent harems) ⎜00 00 0 0 0 ⎟ are the most biased. ⎜ ⎟ 00 00 k(1 − s ) 2k(1 − s ) 3k(1 − s ) 1 1 1 Figure 8 demonstrate how female-biased ( s = 0 ), unbi- ⎜ ⎟ h = 00 00 0 0 0 (24) ⎜ ⎟ ased ( s = 0.5 ), and male-biased ( s = 1 ) harvest strategies h h ⎜ ⎟ 00 00 0 0 0 affect population growth rates and secondary sex ratios. ⎜ ⎟ 00 00 0 0 0 ⎜ ⎟ ⎝ ⎠ 00 00 0 0 0 1. Female-biased harvest (Fig. 8a) most strongly reduces growth in populations with large d and small h, the same ⎛ ⎞ −( +  ) 0 00 0 0 0 m1 m ⎜ ⎟ −  00  + d 00 m m2 f 2 ⎜ ⎟ 00 −( +  ) 00 0 0 ⎜ ⎟ f 1 f ⎜ ⎟ (25) = 00  −   + d 2 + d 3 + d f f 2 m2 m2 m2 ⎜ ⎟ 0 0 00 −( +  + d)  + d 0 m2 f 2 f 2 ⎜ ⎟ ⎜ 0 0 00 0 −( +  + d)  + d ⎟ m2 f 2 f 2 ⎜ ⎟ 0 0 00 0 0 −( +  + d) m2 f 2 ⎝ ⎠ 1 3 harem size h harem size h 30 Population Ecology (2018) 60:21–36 change in λ from harvest change in s from harvest a female-biased harvest (s = 0) 0.28 14 14 -0.08 0.26 -0.09 12 12 0.24 0.22 -0.1 10 10 0.2 -0.11 0.18 -0.12 0.16 6 6 0.14 -0.13 0.12 4 4 -0.14 0.1 -0.15 2 2 01 23 45 01 23 45 divorce rate d divorce rate d b unbiased harvest (s = 0.5) -0.08 14 14 -0.085 -0.015 -0.09 -0.02 -0.095 10 10 -0.1 -0.025 -0.105 -0.03 -0.11 -0.115 -0.035 4 4 -0.12 -0.125 -0.04 01 23 45 01 23 45 divorce rate d divorce rate d c male-biased harvest (s = 1) -0.15 14 -0.06 14 -0.16 -0.17 12 12 -0.07 -0.18 10 10 -0.19 -0.08 -0.2 8 8 -0.09 -0.21 -0.22 6 6 -0.1 -0.23 4 4 -0.24 -0.11 -0.25 2 2 01 23 45 01 23 45 divorce rate d divorce rate d Fig. 8 Population responses to harvest that is a female-biased rate  . (Right) The change in secondary sex ratio s . Other parameters ( s = 0 ), b unbiased ( s = 0.5 ), and c male-biased ( s = 1 ) in the are the same as in Fig. 4 h h h polygynous ( h > 1 ) model. (Left) The change in population growth 1 3 harem size h harem size h harem size h harem size h harem size h harem size h Population Ecology (2018) 60:21–36 31 0.4 Fig. 9 Growth rates  as a male-biased harvest (s = 1) function of the total harvest unbiased harvest (s = 0.5) 0.3 rate E for populations with female-biased harvest (s = 0) various mating systems. The 0.2 four types of mating systems shown correspond to different 0.1 harem sizes h and divorce rates low h, low d low h, high d d (Table 2); in this example, low h = 2 , high h = 10 , low d = 0 , 0.4 and high d = 2 . Harvest may be female-biased ( s = 0 ), unbi- 0.3 ased ( s = 0.5 ), or male-biased ( s = 1 ). Other parameters are h 0.2 the same as in Fig. 4 0.1 high h, low d high h, high d 00.2 0.40.6 0.81 1.21.4 1.61.8 2 00.2 0.40.6 0.81 1.21.4 1.61.8 2 E E populations with the lowest unharvested growth rates normally would (similar to the low d, low h case for (Fig. 7a). Smaller harems are less resilient to female- unbiased harvest). As in the previous scenarios, the biased harvest for the same reason they are less resilient growth of large h populations is less affected by har - to divorce and female mortality - because they cannot vest. Even though male-biased mortality causes unions lose as many females before dissolving. Female-biased to break up more frequently, it also returns (potentially harvest may reduce the average harem size, and also many) females to the f pool. This may be beneficial makes it difficult for large harems to form or reform when h is high, as a higher proportion of females is after breaking up. When h is high, a higher proportion needed to maximize the mating rate in accordance with of females is needed to maximize the mating rate, in Eq. 21. accordance with Eq. 21, but these females are depleted by harvest. Increasing the divorce rate increases the rate As expected, the secondary sex ratio s increases during at which harems dissolve. This more drastically reduces female-biased harvest, decreases during male-biased har- growth for larger harems (contours are steeper at large vest, and undergoes only minimal changes when harvest is h), because they need more females to reform. unbiased (Fig. 8, right). Populations with high d and low h 2. Unbiased harvest (Fig.  8b) yields similar qualitative experience the largest sex ratio shifts under biased harvest. trends. Again, populations with higher divorce rates This may be because the smaller growth rates of high d, low experience greater reductions in growth, and larger har- h populations are less effective in offsetting harvest-induced ems are more affected by divorce. The effect of increas- sex ratio biases. ing divorce is not as pronounced as with female-biased Figure 9 shows how the growth rates of the mating sys- harvest (contours are flatter overall), as less female tems in Table 2 vary with harvest bias and intensity. Again, harvest makes it easier for harems to reform. At low we see that high h, low d populations (large, persistent har- h, however, populations with lower d are actually more ems) have the largest growth rates of all the mating systems, impacted by harvest. Low h unions have only a few even under harvest. Low h, high d populations (small, tran- females and are more likely to dissolve from increased sient harems) have the smallest growth rates. mortality. Unions with high divorce rates are already Increasing the total harvest rate E in Eq.  18 amplifies dissolving quickly, regardless of harvest mortality. the differences between female-biased, unbiased, and Unions with low divorce rates, in contrast, break up male-biased harvest strategies. Female-biased harvest (red) much more frequently once harvest mortality occurs. decreases population growth more severely than male-biased As a result, low d, low h populations experience the harvest (blue) does, even across populations with different h largest decreases in growth. and d. This may be because there is an excess of single males 3. Male-biased harvest (Fig.  8c) reverses the effects of waiting to become harem leaders, whereas single females increased divorce rate. Focusing harvest on males is are usually in shorter supply (especially when h is large). more likely to dissolve harems by killing their male Additionally, the death of a male harem leader immediately leaders. Populations with low d experience the larg- dissolves his union; this can return many females to the sin- est reductions in growth, because male-biased harvest gles stage, allowing new, full-sized harems to reform with makes these unions dissolve more frequently than they a new male leader. The death of female harem members, 1 3 32 Population Ecology (2018) 60:21–36 in contrast, does not necessarily cause the union to dis- and Danchin 2008). Sex-biased harvest may have different solve, and may instead generate small, stunted harems with consequences for offspring survival in these mating systems. reduced productivity. Other species have additional nuances in how they respond Depending on the mating system, unbiased harvest to sex-biased harvest; African lions, for instance, commit (black) can decrease growth rates more or less than female- infanticide when male harem leaders are killed (Whitman biased harvest does. In populations with low h (small har- et al. 2004), which exacerbates the effects of male harvest ems), female-biased harvest has the most drastic impacts of on population growth. all the harvest strategies — again, perhaps, because small How populations respond to selective harvest also has harems cannot afford to lose as many females before dis- important consequences for evolution. Growing evidence solving. In populations with high h (large harems), however, suggests that evolutionary considerations are relevant to sus- unbiased harvest may be just as, if not more, detrimental to tainable long-term management (Ashley et al. 2003), and that population growth. human-induced selection is especially important for harvested species. As harvest mortalities are often more severe and selective than natural mortalities, they may drive evolution Discussion in directions that would not occur under natural conditions. Because it integrates life cycle structure, sex ratio, and The framework we present here is a tool for studying the a mating function, the approach introduced here is ideally effects of sexual reproduction, mating systems, and life his- suited to studying the evolution of sex ratio as a compo- tories on population dynamics. Our implementation of the nent of life history evolution. Sex ratio evolution has a long BMMR approach places no limitation on the complexity of and distinguished history in evolutionary demography (e.g., the life cycle, and is flexible enough to accommodate age or Darwin 1871; Fisher 1930; Trivers 1972; Charnov 1982, stage structure, diverse mating systems, and mating preferences. among many others). Many of these discussions invoke Because it is formulated as a matrix model, powerful sensitivity demographic factors, such as relative mortality of males analysis tools are available (see Shyu and Caswell 2016a, b). and females, but do so without a demographic model that There are interesting open problems to which the incorporates those factors. approach can be applied. As formulated here, the model Our approach provides such a model, and can be analyzed describes a constant environment. Time-varying (seasonal using the framework of adaptive dynamics (e.g., Geritz et al. or stochastic) models would be interesting and challenging 1998), which accounts for the frequency-dependence inher- extensions. Spatial models connecting multiple subpopula- ent in the mating process. We have applied this to the evolu- tions might have important effects on mate limitation. And, tion of sex ratio as influenced by sex-biased offspring costs because the model is formulated in terms of pair formation, (Shyu and Caswell 2016a) and multiple maternal conditions it will be interesting to see how it might apply to species (Shyu and Caswell 2016b). In these studies, the dynamics that do not form pairs, but in which sex ratio may influence of the population frequency vector  are used to evaluate population dynamics (e.g., pollen limitation in plants). the possibility for a mutant phenotype to invade a resident In our application to sex-biased harvest we found that phenotype; phenotypes that are not invasible are singular mating factors, including harem size and union duration, strategies; their stability properties lead to the identification affect not only unharvested population growth, but also the of evolutionarily stable (ESS) or convergence stable strate- responses of growth rate and sex ratio to sex-biased harvest. gies. See Shyu and Caswell (2016a, b) for details. In unharvested populations, high rates of divorce, which lead Acknowledgements This work was supported by a National Sci- to more transient unions, tend to reduce population growth, ence Foundation Graduate Research Fellowship to ES, under Grant especially when harem size is small. Sex-biased harvest 1122374. HC acknowledges support from NSF Grants DEB1145017 affects not only population sex ratios, but also long-term and DEB1257545 and support from the European Research Coun- growth rates, with effects depending on sex bias, harem size, cil under the European Union’s Seventh Framework Programme (FP7/2007-2013), ERC Advanced Grant 322989. ES acknowledges and divorce rate. These complex, and sometimes counterin- support from the Academic Programs office of the Woods Hole Ocean- tuitive, nuances would be impossible to capture without a ographic Institution. demographic two-sex model like this, motivating the use of such models in ecological studies. Open Access This article is distributed under the terms of the Crea- tive Commons Attribution 4.0 International License (http://creat iveco Our two-variable depiction of the mating system, defined mmons.or g/licenses/b y/4.0/), which permits unrestricted use, distribu- by harem size and union duration, can be extended to other tion, and reproduction in any medium, provided you give appropriate factors. Mating systems may differ in parental investment by credit to the original author(s) and the source, provide a link to the males and females. While polygynous males tend to provide Creative Commons license, and indicate if changes were made. minimal parental care, monogamous males invest on par with their female partners (Emlen and Oring 1977; Cézilly 1 3 Population Ecology (2018) 60:21–36 33 Appendix 1: Sequential processes the derivative is dt in continuous time d d‖‖ ‖‖ − dt dt The continuous-time two sex model is obtained by imagining dt ‖‖ (32) that the processes of mating, birth, and transitions take place d‖‖ 1 d = − as a periodic sequence, each over a discrete time interval, and ‖‖ dt dt ‖‖ then shrinking the time interval so that, in the limit, all three processes are operating at every moment. The resulting contin- Let  be a 1 × s vector of ones. Because  is a uous-time projection matrix  is the average of the individual non-negative vector, transition rate matrices in Eq. 3. Recall that the solution of a linear continuous-time system is ‖‖ = (33) d‖‖ d d = (t) → (t)= e (0) (26) dt dt dt Δt Thus, the matrix exponential e projects the population Then Eq. 32 can be rewritten as over a discrete period of length Δt: 1 d  d Δt = − (t +Δt)= e (t). (27) dt   dt ( ) dt Suppose that, e.g., three constant matrices  ,  , and  act 1 d =  −  (34) in sequence, over periods of length Δt . Then dt Δt Δt Δt =  −  (t) (t + 3Δt)= e e e (t) (28) s To first order, as Δt → 0, (t + 3Δt)≈( +Δt)( +Δt)( +Δt)(t) where  is a s × s identity matrix. (29) ≈ [ +Δt( +  + )](t) If the population is initialized with a population vector =  so that ‖‖ =   = 1 , then Eq. 34 can be rewritten as: Dividing by 3Δt and taking the limit as Δt → 0 gives the derivative, d =  −  (t) (35) dt (t + 3Δt)− (t) = lim dt Δt→0 3Δt as in Eq. 10. (30) = ( +  + )(t) Appendix 3: BMMR matrices = (t) for a polygynous mating system Consider a polygynous system with a maximum harem size of as in Eq. 3. Heuristically, this can be thought of as a periodic h. When h is large, it is cumbersome to write the rate matrices product operating so fast that all three processes are operat- ,  ,  in full, especially since many of their entries will be ing simultaneously. In our case, the matrix  is a function zeros. Instead, we will consider these matrices in terms of their of (t) , which changes over the interval [t, t +Δt] . However, contributions to these nine regions of the projection matrix: as Δt gets sufficiently small, the result is equivalent to evaluated at (t). ⎛ ⎞ m→m f→m u→m Technically, this construction is known as the product ⎜ ⎟ (36) m→f f→f u→f integral (Gantmacher 1959). It originated in the mathemati- ⎜ ⎟ ⎝ ⎠ m→u f→u u→u cal work of the same Vito Volterra who is recognized by ecologists for his work on population dynamics (Volterra and For the union formation matrix  , the only regions with Hostinksy 1938). nonzero contributions are: Appendix 2: Dynamics of the relative frequency vector   = 2 × 2 (37) m→m 0 − U Given the population vector of stage frequencies 0 0 = 2 × 2 (38) f→f 0 − hU = f (31) ‖‖ 1 3 34 Population Ecology (2018) 60:21–36 a particular advantage to work with a continuous time model ⎛ ⎞ because it specifies the rates of mating rather than prob- ⎜ ⎟ ⋮⋮ = h × 2 abilities; to specify probabilities would require accounting (39) m→u ⎜ 1 ⎟ 0 U for probabilities of events that did not happen. ⎝ ⎠ The mating process, where adult males and females pair into reproducing unions, is described by the union forma- tion matrix  . Mating functions in  give the rates of union ⎛ ⎞ ⎜ ⎟ formation as functions of the relative frequencies of males ⋮⋮ = h × 2 (40) f→u ⎜ ⎟ and females available to mate, and are thus functions of the 0 U ⎝ ⎠ stage frequency vector . Mating preferences in the mating functions describe the For the birth matrix  , the only regions with nonzero con- probabilities of favoring partners of certain conditions. The tributions are: female preference distribution g (i) gives the proportion ks 2ks … khs of Condition j females that mate with Condition i males. 1 1 1 = 2 × h (41) u→m Similarly, the male preference distribution h (j) gives the 00 … 0 proportion of Condition i males that mate with Condition j females. Summing these distributions over all male and k(1 − s ) 2k(1 − s )… kh(1 − s ) 1 1 1 = 2 × h u→f female conditions respectively yields a total probability of 1: 00 … 0 (42) g (i)= 1 ∀ j (47) For the transition matrix  , the only regions with nonzero contributions are: h (j)= 1 ∀ i −( +  ) 0 (48) m1 m = 2 × 2 (43) m→m m m2 Examples of mating preference distributions include: −( +  ) 0 f 1 f = 2 × 2 (44) f→f 1. Fully assortative mating, where individuals only mate f f 2 with partners in the same condition: 0 0 … 0 g (i)= 1 if i = j, 0 else = 2 × h (45) u→m j + d 0 … 0 f 2 (49) h (j)= 1 if i = j, 0 else 00 … 0 2. Random mating, where individuals pick partners based = 2 × h (46) u→f + d 2 + d … h + d m2 m2 m2 on their relative abundances in the population: The h × h submatrix  is also nonzero. It contains entries u→u g (i)= of −( +  + d) all along its diagonal, and entries of m2 f 2 + d all along its first superdiagonal. f 2 (50) As an example, the  ,  ,  matrices in the case where h (j)= h = 3 are provided in Eqs. 23, 24, and 25 respectively. 3. Biased mating, where individuals prefer partners of cer- Appendix 4: Multiple mating stages tain conditions. An attractiveness or competitiveness and mating preferences factor c weighs the abundance of each partner condi- tion, e.g.,: Suppose that there are q adult male stages and r adult c m female stages. There are qr types of unions possible. The i i g (i)= rate of production of (i, j) matings depends on males of stage c m i i i i and females of stage j, and on preferences. In this case, it is (51) c f j j h (j)= ∑ . c f j j j The material in Appendix  4 is extracted and modified from Shyu and Caswell 2016b under a Creative Commons Attribution license. 1 3 Population Ecology (2018) 60:21–36 35 Cézilly F, Danchin É (2008) Mating systems and parental care. In: Partners with larger c are more preferable mates. If all Danchin É, Girladeau LA, Cézilly F (eds) Behavioural ecology. c are equal, Eq. 51 reduces to the random mating case Oxford University Press, Oxford, pp 429–465 Eq. 50. If c = 0 , individuals of stage i do not mate. Charnov EL (1982) The theory of sex allocation. Princeton Univer- sity Press, Princeton Darwin C (1871) The descent of man and selection in relation to sex. The total mating function M () gives total unions u ij ij John Murray, London (Condition i males mated with Condition j females) formed Emlen ST, Oring LW (1977) Ecology, sexual selection, and the evo- per time. The most general and flexible mating functions lution of mating systems. Science 197:215–223 are based on generalized weighted means (Hölder means). Festa-Bianchet M (2003) Exploitative wildlife management as a selective pressure for life-history evolution of large mammals. These have the general form: In: Festa-Blanchet M, Apollonio M (eds) Animal behavior and wildlife conservation. Island Press, Washington, DC, pp 1∕ a a M ()= b[f g (i)] +(1 − b)[m h (j)] (52) ij j j i i 191–210 Fisher RA (1930) The genetical theory of natural selection. Oxford where 0 ≤ b ≤ 1 and a < 0 (Hadeler 1989; Martcheva and University Press, Oxford Fredrickson A (1971) A mathematical theory of age structure in Milner 2001; Caswell 2001). Note that M () is calculated ij sexual populations: random mating and monogamous marriage only over individuals that are available to mate (i.e., adult models. Math Biosci 20:117–143 single male stages m and adult single female stages f ). As a i j Gantmacher FR (1959) The theory of matrices. Chelsea, New York result, the mating function does not depend on the males and Geritz SAH, Mesze G, Metz JA (1998) Evolutionarily singular strate- gies and the adaptive growth and branching of the evolutionary females in non-mating stages, such as immature juveniles or tree. Evol Ecol 12:35–57 adults already in unions. Ginsberg JS, Milner-Gulland EJ (1994) Sex-based harvesting and The harmonic mean mating function is one of the most population dynamics in ungulates: implications for conservation widely used, because it satisfies the biological criteria for and sustainable use. Conserv Biol 18:157–166 Hadeler KP (1989) Pair formation in age-structured populations. two sex models and captures the qualitative properties of Acta Appl Math 14:91–102 a wide range of Holder means (Caswell and Weeks 1986; Hadeler KP (1993) Pair formation with maturation period. J Math Iannelli et al. 2005). Hence, we use a harmonic mean mat- Biol 32:1–15 ing function where a =−1, b = 1∕2 , so that: Hadeler KP, Waldstätter R, Wörz-Busekros A (1988) Models for pair formation in bisexual populations. J Math Biol 26:635–649 Hardy ICW (2002) Sex ratios: concepts and methods. Cambridge 2m h (j)f g (i) i i j j M ()= . University Press, Cambridge (53) ij m h (j)+ f g (i) i i j j Hoppensteadt F (1975) Mathematical theories of populations: demo- graphics, genetics, and epidemics. Society for Industrial and Applied Mathematics, Philadelphia The corresponding male and female per capita mating func- Iannelli M, Martcheva M, Milner F (2005) Gender-structured tions are: population modeling: mathematical methods, numerics, and simulations. Society for Industrial and Applied Mathematics, M () ij Philadelphia U ()= m,ij m Jenouvrier S, Caswell H, Barbraud C, Weimerskirch H (2010) Mat- (54) ing behavior, population growth, and the operational sex ratio: a M () ij periodic two-sex model approach. Am Nat 175:739–752 U ()= . f ,ij f Karmel PH (1974) The relations between male and female reproduction rates. Popul Stud 1:249–274 Kendall DG (1949) Stochastic processes and population growth. J R Stat Soc B 11:230–282 References Keyfitz N (1972) The mathematics of sex and marriage, vol 4. In: Proceedings of the sixth Berkeley symposium on mathematical Allendorf FW, Hard JJ (2009) Human-induced evolution caused by statistics and probability, pp 89–108 unnatural selection through harvest of wild animals. Proc Natl Kuczynski RR (1932) Fertility and reproduction. Falcon Press, New Acad Sci USA 106:9987–9994 York Ashley MV, Willson MF, Pergams ORW, O’Dowd DJ, Gende SM, Martcheva M (1999) Exponential growth in age-structured two-sex Brown JS (2003) Evolutionarily-enlightened management. Biol populations. Math Biosci 157:1–22 Conserv 111:115–123 Martcheva M, Milner F (2001) The mathematics of sex and marriage, Bessa-Gomes C, Legendre S, Clobert J (2010) Discrete two-sex mod- revisited. Math Popul Stud 9:123–141 els of population dynamics: on modeling the mating function. McFarland DD (1972) Comparison of alternative marriage models. Acta Oecol 36:439–445 In: Grenville TNE (ed) Population dynamics. Academic Press, Caswell H (2001) Matrix population models. Sinauer Associates, New York, pp 89–106 Sunderland Miller TEX, Inouye BD (2011) Confronting two-sex demographic Caswell H, Shyu E (2012) Sensitivity analysis of periodic matrix models with data. Ecol 92:2141–2151 population models. Theor Popul Biol 82:329–339 Miller TEX, Shaw AK, Inouye BD, Neubert MG (2011) Sex-biased dis- Caswell H, Weeks DE (1986) Two-sex models: chaos, extinction, persal and the speed of two-sex invasions. Am Nat 177:549–561 and other dynamic consequences of sex. Am Nat 128:707–735 1 3 36 Population Ecology (2018) 60:21–36 Milner JM, Nilsen EB, Andreaseen HP (2007) Demographic side Shyu E, Caswell H (2016a) A demographic model for sex ratio evo- effects of hunting ungulates and carnivores. Conserv Biol lution and the effects of sex-biased offspring costs. Ecol Evol 21:36–47 6:1470–1492 Pollak RA (1986) A reformulation of the two-sex problem. Demog- Shyu E, Caswell H (2016b) Frequency-dependent two-sex models: a raphy 23:247–259 new approach to sex ratio evolution with multiple maternal condi- Pollak RA (1987) The two-sex problem with persistent unions: a gen- tions. Ecol Evol 6:6855–6879 eralization of the birth matrix-mating rule model. Theor Popul Trivers RL (1972) Parental investment and sexual selection. In: Camp- Biol 32:176–187 bell B (ed) Sexual selection and the descent of man, 1871–1971. Pollak RA (1990) Two-sex demographic models. J Polit Econ Aldine, Chicago, pp 136–179 98:399–420 Volterra V, Hostinksy B (1938) Opérations infinitésimales linéaires. Pollard AH (1948) The measurement of reproductivity. J Insti Actuar Gauthier-Villars, Paris (in French) 74:288–337 Whitman K, Starfield AM, Quadling HS, Packer C (2004) Sustainable Pollard JH (1974) Mathematical models for the growth of human popu- trophy hunting of African lions. Nature 428:175–177 lations. Cambridge University Press, Cambridge Yellin J, Samuelson PA (1974) A dynamical model for human popula- Rankin DJ, Kokko H (2007) Do males matter? The role of males in tion. Proc Natl Acad Sci USA 71:2813–2817 population dynamics. Oikos 116:335–348 Sundelöf A, Åberg P (2006) Birth functions in stage structured two-sex models. Ecol Model 193:787–795 1 3

Journal

Population EcologySpringer Journals

Published: Jun 1, 2018

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off