From female choice to social structure: Modeling harem formation in camelids

Herbivorous wild species constantly strive to optimize the trade-off between energy and nutrient intake and predation risk during foraging. This has led to the selection of several evolutionary traits – such as diet, habitat selection, and behavior – which are simultaneously shaped by the spatio-temporal variability of the habitat. Among camelid species, polygyny is a prevalent behavioral strategy that encompasses both mating and foraging activities. This group-level behavior has multiple interacting dimensions, contributing to an interesting ecological and evolutionary complexity. We developed an individual-based stochastic model in which camelid females transition between different familial groups in response to their environmental conditions, aiming to maximize individual fitness. Our results indicate that the behavioral strategy of individual females can shape, by itself, emergent population-level properties, including group size and fitness distribution. Furthermore, these properties are modulated, in a non-additive manner, by other factors such as population density, sex ratio and system heterogeneity.

💡 Research Summary

The paper presents an individual‑based stochastic model that captures how female camelids (guanacos and vicuñas) choose the harem they belong to in order to maximize their fitness under varying ecological and social conditions. The authors first review the ecological background of camelid societies, emphasizing the trade‑offs that females face between resource acquisition, predation risk, and harassment by bachelor males. They argue that while male reproductive success generally increases with harem size, females experience diminishing returns as group size grows, creating a complex decision‑making landscape that has been difficult to predict from field observations alone.

To formalize this landscape, the authors define a female payoff function Q(F_i, B) that incorporates (i) the quality of a male’s territory (σ_i R_i) divided by the cost of sharing it among F_i females (1 + c F_i), and (ii) the costs associated with environmental hazards (H) and harassment by B bachelor males (β B), both diluted by group size (1 + F_i). The payoff equation is normalized by the male’s subsistence cost, allowing all parameters to be expressed in relative units. A switching penalty η is added to model the energetic or social cost of leaving a current harem.

The model assumes a closed population with fixed numbers of females (F_T) and males (M_T). Each Monte‑Carlo step (MCS) represents a decision round: every female evaluates the potential payoff of joining each male’s harem (including the scenario of forming a new harem with a bachelor male) by inserting herself as the sole new member (F_i + 1). Because payoffs can be negative, the authors shift the entire payoff distribution by adding the absolute value of the minimum payoff, ensuring all values are non‑negative. These shifted payoffs are then normalized to produce a probability distribution P_i for each male, which drives a multinomial selection of the new group for each female. The process is iterated, allowing the system to evolve toward a dynamic equilibrium.

Simulation experiments explore how three key factors shape emergent population‑level patterns: (1) the sex ratio M_T/F_T, (2) overall population density, and (3) heterogeneity in group intrinsic quality (σ_i R_i). In a homogeneous scenario where all harems have identical quality, the proportion of males that retain a harem quickly stabilizes, and fluctuations around this steady state depend mainly on total numbers of males and females. Changing the sex ratio reveals non‑linear effects: a higher proportion of bachelor males (larger β B term) increases the harassment cost for all females, pushing them toward larger harems where the per‑individual risk is diluted. Conversely, when females outnumber males, harem sizes shrink and the distribution of female fitness becomes broader, reflecting increased competition for limited male protection.

Introducing heterogeneity in σ_i R_i produces a pronounced clustering of females around high‑quality territories. In these simulations, a few “super‑harems” dominate, while many low‑quality males remain bachelors, generating a bimodal distribution of harem sizes. The model therefore reproduces the empirically observed pattern that resource‑rich areas attract disproportionate numbers of females, even though the underlying decision rule is purely individual‑based.

The authors discuss the implications of these findings for understanding camelid social organization. By isolating female decision‑making and treating males as passive providers of resources (or sources of harassment), the model demonstrates that many emergent properties of the system—stable numbers of harems, the shape of the fitness distribution, and the sensitivity to sex ratio—can arise without invoking complex male strategies. However, the simplifications also constitute limitations. Real male camelids actively defend territories, engage in coalition formation, and experience mortality and recruitment, all of which are omitted. Moreover, the model parameters (β, H, c, σ_i R_i) are not calibrated against field measurements, so the results are qualitative rather than predictive.

Future work is suggested to incorporate dynamic male behavior, explicit birth‑death processes, and spatially explicit resource maps, as well as to validate the model against longitudinal field data on harem composition, resource use, and predation pressure. Such extensions would allow the framework to move from a conceptual demonstration of emergent complexity toward a quantitative tool for wildlife management and conservation planning.

In conclusion, the study provides a clear example of how individual‑level optimization can generate population‑level structure in a socially complex species. It highlights the importance of sex‑specific trade‑offs, the role of harassment by non‑breeding males, and the impact of environmental heterogeneity on social organization. The modeling approach offers a valuable bridge between ecological theory and the observed dynamics of camelid harems, and it sets the stage for more comprehensive, data‑driven investigations of polygynous mammalian societies.