A two-sex matrix population model to represent harem structure.

Population dynamic models often include males in the calculation of population change, but even in those cases males have rarely been introduced to represent polygyny (harem social structure), where it is particularly important to include males in the reproductive performance of the population. In this article we develop an adaptable matrix population modeling framework for species that have a harem-like social structure under an assumption that the transitions from newborn to juvenile and juvenile to adult both take one time step. We are able to calculate not only the growth rates and stable stage distributions, but also the mathematical expressions for harem size for this model. We then provide applications of this model to two mammal species with slightly different harem behavior.

1. Introduction. The use of two-sex models is becoming more common in demographic studies where there is either an unbalanced sex ratio or a unique social structure. Examples of two-sex models include sperm whales [14], bengal tigers [9], perennial seagrass [23] and red deer [13,20], to name a few. One of the primary difficulties with two-sex models is determining how to incorporate the interaction between the males and females, which is often modeled through a birth function [3,4,25]. In this article, we present a novel stage-structured two-sex model that includes parameters that model the transition between stages that are determined by both age and social structure.
Although the inclusion of males in demographic modeling is not always necessary, the representation of males in matrix models has been increasing. Even though there are various reasons that necessitate including males in the modeling effort, they are generally related to gender specific considerations like predation or survival [18] or population density considerations like sperm limitation and infanticide [22].
Another common reason to include males is when the sex ratio is biased in one direction [18,23].
On the other hand, there are models useful for conservation planning and for the comparison of management strategies that require explicit social categories and/or the representation of both sexes [11,24]. In the case of population management of mammals or birds that have a social structure involving breeding groups, a common strategy is based on the harvest of individuals of one sex or stage/age-group; in many cases when males are regarded as a "surplus" the harvest of only males is usually considered [12]. In these cases, two-sex models with a harem structure allow an optimization of alternative management strategies.
Harem formation is a social structure employed by numerous species (polygyny). Harem size can also be used as a quantitative estimate of the opportunity for sexual selection [27] and has been shown to correlate with the degree of sexual dimorphism in body size [26]. Although there may be an effect on female reproductive strategy [2,17], the focus of this article is on social species with a harem structure where there is a dominant male that has exclusive access to a group of females.
Matrix population models are widely used in ecological modeling [3]. Stagestructured models classify individuals into stages and aggregate all of the individuals within each stage for analyses. Although it is possible to collapse a matrix population model, as suggested in [15], the value of performing and understanding the mathematical analysis can be extremely beneficial.
One of the difficulties of matrix population models is the assumption of a constant sex ratio. A more recent model for populations with a variable sex ratio is frequencydependent matrix models [8]. These models elicit long-term growth rates like matrix population models with a long-term stable stage vector. However, these models are difficult to implement and their sensitivity with respect to the model parameters is extremely difficult to calculate.
Since one-sex population matrices are ubiquitous in population modeling, the focus of this paper is building a two-sex stage-structured matrix population model that incorporates the social structure of harems, while at the same time reproducing the same growth rate as the one-sex matrix population model. In so doing, we will establish properties and behavior of harems while discussing the most common data needed to perform this type of analysis. In addition, we will be able to estimate both harem size and transition probabilities between stages from data.
2. Model development. For the following matrix population models, we assume that the recursion equation is given by n t+1 = P n t , where n t denotes the number of individuals in each stage at time t and P is the particular population matrix. For the two applications in this paper and in order to highlight the effects of the harem formation mechanisms (by simplifying the analysis), we will assume that t is measured in years and the transition rates for the newborn-to-juvenile and juvenileto-adult transitions are of the same one year length. Each of the following matrix population models contains an initial condition n 0 that consists of an initial number of individuals in each of the stages. Given the recursion equation and the initial condition, we can iterate the recursion equation to approximate the population for any time t. For linear population matrices, exponential growth, exponential decay, or constant population are the three possible behaviors. Analysis of these matrices involves finding the maximum eigenvalue, which corresponds to the growth rate of the population, and the corresponding eigenvector, which corresponds to the stable age distribution.
2.1. The classic 1-sex population projection matrix. The population projection matrix for females, P f , is given by where A s is the adult female survival, J s is the juvenile survival, N s is the newborn survival, A f is the female fecundity expressed as female offspring per female adult, and J f is the juvenile fecundity.

2.2.
A 2-sex population projection matrix. We want to design a 2-sex matrix to model harem formation. For the remainder of this paper we will assume that J f = 0, although we could develop a similar model that includes juvenile fecundity.
We are going to begin by taking the only category with fecundity, in this case adults, and make one category for males and the other category for females. The projection matrix for the population of males and females, P m , is given by where A s is the adult survival for males and females, J s is the juvenile survival, N s is the newborn survival, A f is the number of female newborns per female, and γ is the proportion of newborn females. Also recall that we assume that the categories of newborns and juveniles contain both males and females. Since γ is defined as the proportion of newborn females, we have γ = , where A m is number of male newborns. Notice that by cross-multiplying we arrive at where the right-hand side is the female fecundity in terms of γ.

2.3.
The harem 2-sex population projection matrix. In this section we assume that there are two groups of males. One group, called the breeding males, are the dominant (reproductive) males and the other group, called bachelor males, do not mate. In a similar fashion, we can derive a two-sex model where the projection matrix is given by where α is the transition proportion from bachelor male to breeding male, β is the transition proportion from breeding male to bachelor male, and γ is the proportion of females in the population. We again assume that the categories of newborns and juveniles contain both males and females. We also assume that the survival for the entire bachelor male population is A s , which is broken up into α that transition to become breeding male and A s − α that remain as bachelor males. The parameter β behaves in the same way for breeding males. A lower bound for β is 0; β = 0 means that the loss in breeding males is only by mortality. In general, 0 ≤ α, β ≤ A s . Note that we again assume that the female fecundity is the same as (2) to account for births of both males and females.
2.4. The harem 2-sex population projection matrix with male and female transition. In this section, we assume that there are two groups of males and females. In a similar fashion, we can derive a two-sex model where the projection matrix is given by Ns where α m is the transition proportion from bachelor male to breeding male, β m is the transition proportion from breeding male to bachelor male, α f is the transition proportion from bachelorette female to breeding female, β f is the transition proportion from breeding female to bachelorette female, and γ is the proportion females in the population. We also assume that the categories of newborns and juveniles contain both males and females.
2.5. The density-dependent 1-sex population projection matrix. We also examine the impact of nonlinearities on population growth. We define the projection matrix to be where F ∆ is a nonlinear function. Although density dependence may affect survival or fecundity, in this article we only apply density dependence to fecundity, which will be justified when we apply the density-dependent model to guanacos.
2.6. The harem 2-sex population projection matrix with density-dependence. We consider the following two-sex population model with the following projection matrix where F ∆ is a nonlinear function.

Eigenvalue analysis.
In this section, we use the general definition of P f given by (1) with J f = 0, P m given by (2), P 2m given by (3), and P 2mf given by (4) to examine the behavior of the eigenvalues and eigenvectors of the projection matrices. The three eigenvalues of the female projection matrix P f are given by the solution of for λ, referred to as the characteristic equation. The three eigenvalues of P f are and . Notice that V is always greater than or equal to zero. For projection matrices, we are interested in the dominant or largest eigenvalue, which is λ 1 for P f . This dominant eigenvalue is often referred to as the asymptotic growth rate for the population. It is straightforward to show that λ 1 is always bigger than λ 2,3 by comparing the magnitude of the complex eigenvalue λ 2,3 with the square of λ 1 .
The four eigenvalues of the harem matrix P m are given by the solution of Notice that the cubic term is exactly the same equation as (7) in P f , so the eigenvalues λ 1 , λ 2 , and λ 3 of P m are the same as P f . The new eigenvalue λ 4 = A s is from the linear piece. Notice that the eigenvalues do not depend on γ. The appendix shows that λ 1 is always greater than the additional eigenvalue A s , which means that λ 1 is the dominant eigenvalue for P m . The five eigenvalues of the harem matrix P 2m are given by the solution of Notice that the cubic term is again exactly the same equation as in P m , so the eigenvalues λ 1 , λ 2 , and λ 3 of P 2m are the same as P f with the addition of two eigenvalues, namely λ 4 = A s and λ 5 = A s − α − β. Similarly, the eigenvalue λ 1 is the dominant eigenvalue for P 2m , as was the case with P f and P m . This result confirms that indeed P f , P m , and P 2m have the same dominant eigenvalue. Also, notice that the dominant eigenvalue does not depend on either α, β, or γ. Therefore, we have developed a two-sex matrix whose growth rate is exactly the same as the one-sex matrix, but has added complexities that will become important in the stable stage distribution (next section).
The eigenvalues of P 2mf start to get more complicated and we will not write down explicit eigenvalue results. However, we can still obtain information about the matrix by examining the characteristic equation, which is given by This result is very similar to the previous case, except now the maximum eigenvalue will come from the quartic term, but the important similarity is that the dominant eigenvalue does not depend on α m , β m or γ, but it does depend on α f and β f . Therefore, the dominant eigenvalue depends only on the female transitions, α f and β f .

Eigenvector analysis.
Recall that the dominant eigenvalue represents the growth rate of a finite population and the dominant eigenvector represents the stable age distribution, which are shown through linearization techniques. Since both P f and P 2m have the same dominant eigenvalue, they have same growth rate. The dominant eigenvector v 1 f corresponding to the dominant eigenvector λ 1 for P f is given by where a = (1 + Ns Recall that the dominant eigenvector gives the proportion of the population that are newborns, juveniles, and adults, respectively.
The dominant eigenvector of P 2m is given by where a is the same we defined for (9). Recall that the dominant eigenvector gives the proportion of the population that are newborns, juveniles, bachelor males, breeding males, and females, respectively. Notice that the only effect of α and β is to change the stable stage distribution between bachelor males and breeding males and the only effect of γ is on the proportion of males and females.
3.3. Harem size. We can also use the eigenvector (10) to obtain the harem size, h, which is the proportion of adult females divided by the proportion of breeding males. We find that Another method to calculate harem size is given in [27]. Define p 0 to be the fraction of non-mating males out of the total population of adult males, which we refer to as bachelor males. We first calculate the total number of adult males to be (1 − γ) aλ 1 A f by adding the number of bachelor and breeding males in (10). We find that p 0 = λ1−As+β λ1−As+α+β . Then define p s to be the fraction of successful males, which we call breeding males and is given by p s = 1 − p 0 = α λ1−As+α+β . Define R as the sex ratio expressed as the ratio of the total number of females to the number of males. In terms of γ, we define R = γ 1−γ . Finally, define the harem size to be R ps , which gives the same result as (11).
One result of (11) is the relationship between the fraction of juveniles that are female, γ, and the harem size for different values of α. Using values from Table 1, Figure 1 is very similar to a graph in [27] that plots harem size against the fraction of sexual selection for a survey of 27 different harem-type species over 17 taxa. Note that Figure 1 has harem size on the y-axis, whereas in [27] harem size is on the xaxis; however, (11) displays the power of this modeling methodology. Furthermore, notice that when γ = 0.5, decreasing α from 0.35 to 0.05 increases the harem size from approximately 1.5 to approximately 4.5.   Table 1 using (11). 3.4. Sensitivity and elasticity. Another advantage of having a closed form solution for the dominant eigenvalue is the ability to calculate sensitivities exactly.
Since we have exact values for the eigenvalues and eigenvectors, it is much easier to explore the impact of changing the parameter values. For instance, in order to calculate the impact of changing adult survival A s on the dominant eigenvalue λ 1 , simply find ∂λ1 ∂As or differentiate (8) with respect to A s . As a more tangible example, consider P f in (1) with juvenile fecundity and examine the impact of including juvenile fecundity on the dominant growth rate. We find that the characteristic equation is given by Using the parameter values in Table 1, we find that by increasing J f from 0 to 0.1, increases the dominant eigenvalue from λ 1 = 1.044786 to λ 1 = 1.052971. Therefore, a 0.1 increase in juvenile fecundity only increases the asymptotic growth rate by 0.8 %. Note that since we have an exact representation of the eigenvalue, we can get an exact value of the impact of changing J f on the eigenvalue without the need for the typical sensitivity analysis used with projection matrices.

Applications.
3.5.1. Application to Bengal tigers. A two-sex matrix population model was used to model the Bengal tiger population for fixed values of α m and α f [9]. We can apply the model in (4), with β m = β f = 0, for the Bengal tiger population. Notice in (7) that the only effect of N s is in the third term. So, in [9] they combine N s into the term for A f . For comparison sake, we will follow this convention. Both females and males transition to the breeding group, which we will refer to as α f and α m , respectively. The stages for the following model are juveniles, transient (bachelor) males, breeding males, transient females, and breeding females. The population projection matrix is given by where the values of α m = 0.256, α f = 0.288, A f = 1.333, N s = 0.6, J s = 0.9, and γ = 0.5 were estimated from [9]. Notice that the probability of survival is different for transient and breeding males and females. All of the previous parameter values come directly from [9] and are based on field data. The dominant eigenvalue, which is given by 1.0534, and the stable age distribution, which is given by (0.2278, 0.1554, 0.1570, 0.1598, 0.3000), both match the exact values from [9]. Notice that the estimated ratio between breeding males and females is approximately 1:2, but it was found to be 1:3 in the field [9]. The strength of our model is the ability to examine how changing both α m and α f impacts the ratio between breeding males and females, and hence the harem size. Recall that the dominant eigenvalue only depends on α f and does not depend on α m . Therefore, by changing the parameter value α m , we have the same growth rate, but different harem sizes. If we fix α f = 0.288, we find that when α m = 0.132 the observed harem size of 1:3 is found. This result suggests that if there is a high confidence in the female transition proportion, then the actual male transition proportion may be smaller than the calculated value.
However, a more plausible problem is that both α f and α m are not correct. Figure 2 shows how the harem size changes as a function of both α m and α f for (12). Any combination of values for α m and α f on the red line has a harem size of 3, which is the experimental ratio. So, any of these combinations of α m and α f are more likely candidates for the actual transition proportions.
3.5.2. Application to guanacos. In this model, we have chosen to apply the nonlinear function only to the adult fecundity. Very little is known about the population regulation mechanisms of guanaco populations, and a comparative analysis of the inter-annual variation in fecundity rates in three species of ungulates (Red deer, Saiga antelope, and Soay sheep) demonstrated [5] that in all three species fecundity  rates are associated with both density and weather in similar manners. We decided to apply density-dependence to fecundity alone while no climatic effects were used in this model.
We will show that the results of the nonlinear models are much easier to understand due to the results of the prior mathematical analysis. We will use the data for guanacos given in Table 1 [21]. We will first analyze the results of (5) where F ∆ is a nonlinear function with a Ricker-type nonlinearity such that F ∆ = (1 + exp((ρ − a) · b)) −1 , where a > 0 and b > 0. For this example, ρ is the density of the total population as a proportion of the carrying capacity and both a and b are prescribed from data. Although the following results are independent of a and b in the nonlinear model, we will use a = 0.7939, b = 3.74469, which are fit from data [21]. We will also assume that the population is near the carrying capacity. We find that the stable age distribution of 10.5% newborns, 6.9% juveniles, and 82.7% adults agrees with (9). For the female nonlinear model given in (5), the stable age distribution is the same as the linear model. Now we will analyze the two-sex nonlinear model given by (6). There is a nice relationship between the harem size and α as shown in (11), which prompts the  Table 1. Parameter values for a guanaco population as required by the two-sex nonlinear model given by (6). Note that both β and γ were not discussed in [21]. question of whether it is possible to let α vary while keeping the harem size constant. Figure 3 shows the results when α is allowed to vary, β = 0 and we fix the harem size to 10. In this case, we can use (11) in order to predict the steady state value of α. We can use the mathematical results in order to obtain the final stage distributions for the nonlinear model. Since this particular nonlinear problem goes to a steady state that means that the long-term eigenvalue given in (8) is equal to 1. The nonlinear component of P ∆ f is in the upper right entry, so we want to solve (8)=1 for A f , because it is in the same place as the nonlinear component. Performing this substitution gives A f = 0.245. In order to calculate the stable age distribution we need to use (9) with λ 1 = 1 and A f = 0.245, which agrees with the earlier numerical result. Although this mathematical analysis gives us some insight into this particular problem, this type of analysis is not often possible.
Similarly to the one-sex population projection matrix, we are able to calculate the stable age distribution at steady state. However, we can also use (11) to find the harem size for a given value of α, where we once again assume that λ 1 = 1 at steady state. For instance, when α = 0.035 the long term harem size is 4.714.
We can also analyze the impact of sex ratio on harem size by using the nonlinear two-sex population projection matrix given in (6), where γ is the proportion of females in the population. We find that the proportion of newborns and juveniles remains the same as the one-sex model and the adult proportion is divided among bachelor males, breeding males, and females. The proportion of females is γ times the proportion in the one-sex model and α and β control the ratios between bachelor and breeding males. New values of α, β, and h emerge due to the nonlinearity in the fecundity. Figure 4 demonstrates the impact of changing α and γ on the harem size. We are again able to see that as the sex ratio, γ, changes, so does the value of α in order to keep a constant harem size.  4. Discussion. We begin with an all-female matrix population model (1) and construct two-sex populations matrices, P m (2) and P 2m (3), that have the same growth rate, and do not depend on the primary parameters α (the transition proportion from bachelor male to breeding male), β (the transition proportion from breeding male to bachelor male), and γ (the proportion of females). However, these parameters do affect the social category distribution. We also construct a two-sex, male and female transition matrix P 2mf (4) with the same growth rate as P f , the parameter values not only change the stable age distributions, but α f and β f also affect the population growth rate. Our mathematical analysis leads to an expression that describes the harem size under the assumption that newborn-to-juvenile and juvenile-to-adult transitions both take one time step. Furthermore, we can write the harem size in terms of the sex ratio, γ 1−γ , α, β and the other demographic parameters of the species. Figure 1 shows the impact that changes in α m , (the transition proportion to become a breeding male), has on the harem size. We notice that α m must increase for a fixed harem size if the proportion of females, γ , increases. This result shows that more bachelor males are needed to lead harems in order to keep the same fixed harem size. We then compared our results with results from [27], which show that in either formulation, we arrive at the same relationship for harem size.
Even though the mathematical expressions are complex, the ability to take derivatives and obtain exact values for sensitivities is a definite advantage. We examined one example where we were able to calculate the impact of including juvenile fecundity on the growth rate, but this type of analysis can be done with any of the parameter values. However, even a slight modification of the matrices, like in P 2mf , will require using traditional means to calculate sensitivity and elasticity [3].
Nonlinear models are typically difficult to understand, even though they demonstrate more realistic results. We proceed to demonstrate that we can better understand our nonlinear models, because of our thorough understanding of the mathematical analysis of the matrix population models. A key finding with these nonlinear models is the ability to observe the behavior of α m for a fixed harem size.
Our new population matrix model provided new insights in its application to the two species used as examples; regarding the Bengal tigers, it is interesting to see that the observed harem size (equivalent to territory size of breeding males) of three mature females per mature male can be obtained by different combinations of α f and α m ; of course, there will probably be a limit to the possible combinations of these two parameters in terms of feasible biological maturation rates, population densities, proportion of breeders and transients, or territory size of breeding males.
We applied P 2mf (4) to a Bengal tiger model in order to reconcile the difference in the calculated harem size and the field value [9]. Recall that we included a transition term for both males and females, α m and α f , respectively. Since the stable age distribution does not depend on α f , we are able to change the value of α m until the calculated ratio agrees with the field value. This result demonstrates the usefulness of the two-sex matrix population model, given in (4), because we are able to reconcile the difference in computed and field male/female ratios, while focusing on the transition proportions of both the males and females.
Also, interesting conclusions and insights can be obtained from the application of the new model to the guanaco populations; with our average α m value of 0.035 the harem size is 5.6 females/FG (FG= family group), which conforms well with the value of 6 females/FG (N= 925) estimated by Jefferson [10] (cited in [6]) for guanacos of Tierra del Fuego (Chile); additionally Puig and Videla [19], based on six different bibliographic sources, presented mean female/FG value of 5.6 (stand. dev. 1.4, coef. var 24.1%). On the other hand, under steady state conditions (λ = 1) and after 10-15 generations the parameter α stabilizes around a value of 0.01, and when we apply equation (11) with this α value, the predicted harem size is around 14 females/FG; despite this relatively high number of females /FG, such extremes have been found in the field [16] although field values of this size are rare (only 5% of 637 FG observed had values greater than 10 females/FG); however, this result conforms well with our model's prediction at steady state conditions, which probably are seldom reached under field populations.
Remember that all of the results in this article are based on the assumption that the transition from newborn to juvenile and juvenile to adult are the same length. The mathematical analysis was greatly simplified due to the presence of numerous zeros in the population matrices. This assumption can be relaxed, but most of the typical analysis methods would need to be used. For instance, instead of finding the sensitivity of the eigenvalue to changes in a particular parameter by taking derivatives with respect to that parameter, one would need to use traditional methods to calculate the sensitivity [3].
These models can be extended to study the development and optimization of conservation and management plans/program for a species; e.g., our model can show, in cases where population size is greatly reduced, that there may be a risk of inbreeding (a process of concern in biological conservation), as Ballou and Ralls [1] showed for small populations of ungulates. In terms of management, Lancia et al. [12] and Jensen [11] studied a sustained and intense male-only hunting in white tailed deer populations, and the former claimed that the male-only hunting may affect future recruitment. Probably, the results of this male-only hunting strategy will be influenced by the social features or structure of the species; currently, for many species with social structures similar to that of the guanaco, management programs rarely take into account this structure; this is due, in part, to the lack of models to represent it, and so we are confident that our harem model could be a step in correcting this situation. For instance, we could extend (4) to include harvesting using both nonlinear modeling and incorporating harem size, as well as other neglected demographic features. For example, the Allee effect must be considered when studying the impact of proportional harvesting on population growth. With our models we could evaluate how this harvest affects population dynamics: for a given α m parameter, a large decrease in the male group (bachelor males) may lead to a decrease in reproductive males (breeding males), and thus may cause either an increase in the size of the reproductive groups (harem or family groups) or a decrease in the number of these groups, decreasing the population growth rate. On the other hand, we could use our model to estimate the number of males that can be hunted without affecting the population growth rate, given a more robust demographic model.
Other possible applications of our model can be extended to the subject of evolutionary studies; e.g., Wittenberger [29] analyzed the evolution of polygamy and social behavior in mammals based on the costs and benefits to female survival and reproductive rates. That study also evaluates how these costs and benefits may determine the optimal size of female only groups. Nevertheless, Wittenberger [29] did not carry out a quantitative analysis, something we can carry out with our model, and determine quantitatively how the size of harem and other groups affect the demographic parameters, and estimate the optimal harem and bachelor size that maximize fitness. In particular our model can also be used to analyze possible positive frequency-dependent selection as an explanation for differences in different social-group behaviors in different species with similar population structures, as Franz et al [7] have analyzed for cooperative behavior. Furthermore, Wade et al. [28] had to resort to an explicit population genetic model for the evolution of a maternal-effect that biases offspring sex ratio, while such a model based on the demographic parameters is still lacking.
Appendix. We need to show that λ 1 > A s , so that the dominant eigenvalue for (2) and (3) is λ 1 . Substituting (8) gives Rewriting the right hand side while assuming that A s = 0 gives After some algebra, this equation becomes which is true as long as 2A s = V . Recall that or V = (8A 3 s + stuff) 1/3 , where stuff is a positive number. Now it is easy to see that unless all the other variables are zero, V > 2A s and therefore λ 1 is the dominant eigenvalue for (2) and (3).