Mechanisms of reciprocity and diversity in social networks: a modeling and comparative approach

Mechanisms of reciprocity and diversity in social networks: a modeling and comparative approach Abstract Three mechanisms have been proposed to underlie reciprocation of social behaviors in gregarious animals: “calculated reciprocity,” “emotional bookkeeping,” and “symmetry-based reciprocity.” Among these explanations, emotional book-keeping has received the broadest support from experimental and observational studies. On the other hand, 3 individual-based models have shown that reciprocation may emerge via “symmetry-based reciprocity,” “emotional bookkeeping,” or a combination of both mechanisms. Here, we use these 3 models to assess their relative fit with empirical data on reciprocation and social network structure across different groups and species of macaques. We collected grooming data from 14 groups and 8 macaque species and simulated each group in each model. We analyzed and quantitatively compared social network metrics of the empirical and the models’ grooming networks. The 3 models captured fairly well the features of observed networks, and fitted data from wild groups better than captive ones. The emotional bookkeeping model seemed on average to fit slightly better the social networks metrics observed in empirical data, but failed to reproduce some grooming patterns. The symmetry-based models, on the other hand, fitted better other network parameters (e.g., modularity). No model generally fitted the data better than the others, and the fit with some metrics (e.g., modularity, centralization index) was low even after optimization. Thus, our analyses indicate that in the models social interactions may be simpler than in reality and models may miss social processes (e.g., third-party awareness). INTRODUCTION Group living is a common phenomenon across the animal kingdom, from invertebrates (e.g., ants, bees, spiders) to marine (e.g., dolphins, whales, fish) and terrestrial vertebrates (e.g., rodents, primates, felines, birds) (Krause and Ruxton 2002; Ward and Webster 2016). Group living may be widespread due to the many adaptive advantages it has for individuals, including the reduction of predation risks, better access to resources, and community breeding and rearing of offspring (Crook et al. 1976). In species of animals living in permanent, stable groups, individuals engage frequently in social interactions with other group members under many different circumstances and during periods of time that can span their whole lifetime, giving then rise to complex social systems (Freeberg et al. 2012). In such systems, individual recognition and repeated interactions between certain group members may promote the emergence of reciprocation, i.e., the exchange of beneficial behaviors over long-time periods. Patterns of reciprocation have been extensively studied in primate social systems—e.g., grooming and support in fights—and their ultimate functions and proximate mechanisms largely debated. At an ultimate level, reciprocation may be explained via kin selection (Hamilton 1964), reciprocal altruism (Trivers 1971), mutualism (Clutton-Brock 2009), and/or the biological market theory (Noë and Hammerstein 1994). At a proximate level, 3 mechanisms, differing in the degree of cognition required, have been proposed. At the highest cognitive level is “calculated reciprocity” (de Waal and Luttrell 1988). This mechanism implies keeping record of all services given and received from others, long-term memory of past events, and paying back accordingly (de Waal and Brosnan 2006). A cognitively simpler mechanism is “attitudinal reciprocity” (de Waal 2000), also named “emotional bookkeeping” (Schino and Aureli 2009). Emotional bookkeeping means that individuals develop, through social interactions, a specific emotionally mediated attitude toward each group member, and whether or not they exchange a service with a group member depends on the type of affiliative relationship associated with it (Schino and Aureli 2009). At the lowest cognitive level is “symmetry-based reciprocity.” This mechanism assumes that the balance between services given and received results from repeated interactions between individuals sharing a common or symmetrical variable—e.g., spatial proximity, kinship, rank, age, etc. (de Waal and Luttrell 1988). Currently, researchers find calculated reciprocity an implausible hypothesis because of the difficulty that represents keeping record of behaviors over long periods of time (Stevens and Hauser 2004; Stevens et al. 2011). Symmetry-based reciprocity has been disregarded because even when symmetrical variables are partialed out, the correlation between behaviors given and received remains significant, suggesting thus that the correlation is not a by-product of the symmetrical variable. The statistical test used to partial out a variable, however, may be insufficient to completely exclude the effects of the variable (Hemelrijk and Puga-Gonzalez 2012); symmetry-based reciprocity thus, has potentially been underestimated. Emotional bookkeeping is nowadays the leading hypothesis regarding the mechanisms giving rise to reciprocity of social behaviors. Indeed, several empirical studies have shown that individuals prefer to associate and reciprocate behaviors with some partners rather than with others (Tiddi et al. 2011; Sabbatini et al. 2012; Jaeggi et al. 2013) and that “social bonds” are important for the individuals’ survival and that of their offspring (Silk et al. 2003; Silk 2007; Schüelke et al. 2010; McFarland et al. 2015). An approach to the study of the proximate mechanisms underlying social behavior in nonhuman primates makes use of individual-based computer models and complex science—i.e., the study of phenomena emerging from the interactions among the elements of a system. Several different individual-based models have shown that simple behavioral rules and local interactions may suffice to explain many social patterns observed in primate societies (Evers et al. 2011; Campenni and Schino 2014; Evers et al. 2015; Hemelrijk et al. 2017). In this regard, 3 individual-based computer models—GrooFiWorld, FriendsWorld, and Reaper—have been able to reproduce not only social behaviors but also the behavioral contrasts observed between tolerant and intolerant societies of monkeys, especially macaques (Puga-Gonzalez et al. 2009; Hemelrijk and Puga-Gonzalez 2012; Puga-Gonzalez et al. 2015; Puga-Gonzalez and Sueur 2017a). The Macaca genus comprises 23 species in which females are the philopatric sex, which results in female-bonded societies. On the basis of the females’ social behavior, macaque species have been classified along a social style gradient going from extremely intolerant to extremely tolerant. Among females from species in grade 1 and 2 (extremely intolerant and intolerant), the dominance hierarchy is steep, conflicts are unidirectional, aggression is fierce, reconciliation is infrequent, and social grooming is directed up the dominance hierarchy and towards individuals of similar rank. On the contrary, among females from species in grades 3 and 4 (tolerant and extremely tolerant), the dominance hierarchy is more shallow, conflicts are more bidirectional, and aggression is mild (Thierry 2004). Further, in all 4 grades individuals cooperate by reciprocating (exchange of same behaviors) and interchanging (exchange of different behaviors) grooming and support in fights (Schino 2007; Schino and Aureli 2008a). All these behavioral patterns (reciprocation and interchange of grooming and support), as well as the behavioral differences between intolerant and tolerant societies (grooming up the dominance hierarchy and towards individuals of similar rank, uni- (bi-) directional aggression, high [low] steepness of the hierarchy), have been reproduced by the 3 individual-based models. The individual-based model GrooFiWorld proposes that spatial proximity (i.e., proximity-based reciprocity), is the main cause of the emergence of complex behavior in societies of macaques (Puga-Gonzalez et al. 2009). In GrooFiWorld, a spatial structure with dominant individuals in the centre and subordinates at the periphery emerges as a side-effect of aggression. Due to this spatial structure, individuals interact more with some than with others, and reciprocity of grooming and support in fights; and the interchange of grooming and support emerge. Interestingly, these patterns disappear in the model only when proximity is omitted by making individuals interact at random, but not when proximity is statistically partialed out with the partial Tau-kr test (Hemelrijk 1990a), which indicates that partial correlations do not sufficiently exclude the effect of proximity (Hemelrijk and Puga-Gonzalez 2012). Recently, GrooFiWorld was extended with a mechanism similar to emotional-bookkeeping (Puga-Gonzalez et al. 2015). In this model, called “FriendsWorld,” individuals have a tendency to follow their “friends”, i.e., those individuals with whom they groom the most. The model shows that this mechanism of “follow your friends,” reinforces patterns of reciprocity and interchange of grooming and support in fights. It remained unclear however, what factor was more important for the emergence of reciprocation and interchange, the spatial structure of the group, or the “follow your friends” rule. More recently, a new individual-based model has been developed. This model is identical to GrooFiWorld regarding the way individuals choose to interact with their partners. Contrary to GrooFiWorld, this model is not spatially-explicit; thus, interactions are not based on proximity. Instead, individuals interact first with others at random in this model, and as time goes by they develop a preference for interacting with those with whom they exchange more grooming, resembling thus emotional bookkeeping; we called this model Reaper (i.e., reinforcement of affiliative preferences) (Puga-Gonzalez and Sueur 2017a). These 3 models have been successful at reproducing empirical patterns of reciprocation and, at the same time, suggest different factors as potential causes. In these studies, however, no social network analysis (SNA) was done. Given that SNA provides us with a detailed picture of how individuals interact and distribute grooming behavior at the individual and group level; here, we used a SNA to distinguish which of these models and underlying mechanisms (proximity-based reciprocity, emotional bookkeeping, or a combination of both), if any, more accurately represents empirical data. We studied individual and global network metrics to assess whether in the model the individuals’ connections in the social network as well as the global features of the network were similar to empirical data. We focused on 4 different individual metrics: betweenness, closeness, eigenvector, and clustering coefficient; and 7 different global metrics: diameter, number of triangles, global clustering coefficient, modularity, average eigenvector, centrality index, and the correlation coefficient between individuals’ dominance rank and eigenvector (see Table 1 for an interpretation of the meaning of each network metric). We chose these metrics because they have been shown to be relevant in studies of primate social networks. For instance, Kasper and Voelkl (2009) suggest the use of betweenness and closeness (for sparse networks), and eigenvector (for near-to-all-connected networks) to understand the involvement, popularity, and sociability of an individual within the social network. Modularity is used to measure the degree of fragmentation of the network into subgroups, and clustering coefficient to characterize the cliquishness of the network by measuring the extent to which individuals linked to one individual are also linked to each other (Kasper and Voelkl 2009; Sueur, Petit et al. 2011; Pasquaretta et al. 2014). The correlation between eigenvector and dominance rank is used to measure the degree of popularity of dominant individuals in the group (Sueur, Petit, et al. 2011). To measure the centralization of the network around the most central individual, it is recommended to calculate an index obtained by dividing the centralization of the observed network by the maximum centralization possible—that of a star network (Kasper and Voelkl 2009; Sueur, Petit, et al. 2011; Pasquaretta et al. 2014). The approach taken in the present study was exploratory without clear predictions about which network measure should be best predicted by which model. Table 1 Interpretation of the network metrics used in the network analysis Metric Interpretation Individual: Betweenness The degree to which an individual lies in-between the paths of others, and thus the influence that exerts on the passing of information between others. Closeness How close an individual is to all others in the network, and thus how fast information originating from this individual may reach all others Eigenvector “Popularity” (centrality) of an individual; takes into account how well connected are the focal individual and its neighbors (individuals linked to the focal) Clustering coefficient The relative number of cliques (of 3 individuals) an individual forms with its neighbors. Global: Diameter How fast information travels between the 2 least well connected individuals in the network. Number of triangles The absolute number of cliques (of 3 individuals) Clustering coefficient The relative number of cliques (of 3 individuals) Eigenvector The average degree of “popularity” among individuals in the network Modularity The degree of subgrouping or clustering of the network Centrality index The extent to which a network is dominated by the most or most “popular” individuals in the network Correlation coefficient eigenvector-dominance rank The extent to which dominant individuals are the most “popular” individuals within the network Metric Interpretation Individual: Betweenness The degree to which an individual lies in-between the paths of others, and thus the influence that exerts on the passing of information between others. Closeness How close an individual is to all others in the network, and thus how fast information originating from this individual may reach all others Eigenvector “Popularity” (centrality) of an individual; takes into account how well connected are the focal individual and its neighbors (individuals linked to the focal) Clustering coefficient The relative number of cliques (of 3 individuals) an individual forms with its neighbors. Global: Diameter How fast information travels between the 2 least well connected individuals in the network. Number of triangles The absolute number of cliques (of 3 individuals) Clustering coefficient The relative number of cliques (of 3 individuals) Eigenvector The average degree of “popularity” among individuals in the network Modularity The degree of subgrouping or clustering of the network Centrality index The extent to which a network is dominated by the most or most “popular” individuals in the network Correlation coefficient eigenvector-dominance rank The extent to which dominant individuals are the most “popular” individuals within the network See Methods for a technical description of each network metric. View Large Table 1 Interpretation of the network metrics used in the network analysis Metric Interpretation Individual: Betweenness The degree to which an individual lies in-between the paths of others, and thus the influence that exerts on the passing of information between others. Closeness How close an individual is to all others in the network, and thus how fast information originating from this individual may reach all others Eigenvector “Popularity” (centrality) of an individual; takes into account how well connected are the focal individual and its neighbors (individuals linked to the focal) Clustering coefficient The relative number of cliques (of 3 individuals) an individual forms with its neighbors. Global: Diameter How fast information travels between the 2 least well connected individuals in the network. Number of triangles The absolute number of cliques (of 3 individuals) Clustering coefficient The relative number of cliques (of 3 individuals) Eigenvector The average degree of “popularity” among individuals in the network Modularity The degree of subgrouping or clustering of the network Centrality index The extent to which a network is dominated by the most or most “popular” individuals in the network Correlation coefficient eigenvector-dominance rank The extent to which dominant individuals are the most “popular” individuals within the network Metric Interpretation Individual: Betweenness The degree to which an individual lies in-between the paths of others, and thus the influence that exerts on the passing of information between others. Closeness How close an individual is to all others in the network, and thus how fast information originating from this individual may reach all others Eigenvector “Popularity” (centrality) of an individual; takes into account how well connected are the focal individual and its neighbors (individuals linked to the focal) Clustering coefficient The relative number of cliques (of 3 individuals) an individual forms with its neighbors. Global: Diameter How fast information travels between the 2 least well connected individuals in the network. Number of triangles The absolute number of cliques (of 3 individuals) Clustering coefficient The relative number of cliques (of 3 individuals) Eigenvector The average degree of “popularity” among individuals in the network Modularity The degree of subgrouping or clustering of the network Centrality index The extent to which a network is dominated by the most or most “popular” individuals in the network Correlation coefficient eigenvector-dominance rank The extent to which dominant individuals are the most “popular” individuals within the network See Methods for a technical description of each network metric. View Large We collected data from 14 groups belonging to 8 different macaque species, comprising the 4 grades of social style, and did a SNA on each group. Due to limitations in empirical data and because social styles are based on females’ social behavior, our analyses were restricted to grooming social networks of adult females. Then, we simulated the 14 macaque groups in each of the 3 models, analyzed the resulting social networks, and quantitatively compared the results of the models and those of the empirical groups. As a control, we built an extra individual-based model (the “null-model”) in which individuals randomly selected interaction partners, and also randomly chose whether to groom or fight. Before running simulations, we mimicked the following variables from empirical groups in each of the models: group size, sex ratio, intensity of aggression, and density of grooming network. After running the simulations in the models we realized that there was a poor fit between the models and empirical data regarding some global network metrics (e.g., modularity, centralization index). These network metrics may be influenced by the way individuals distribute grooming among group members: whether or not grooming is reciprocated, directed up the dominance hierarchy, and/or toward partners of similar rank. Hence, to optimize the fit between models and empirical data we decided to mimic in the models the distribution of grooming observed in the empirical groups (see Methods for details). Previous analyses of the grooming networks from GrooFiWorld and Reaper models showed that these networks resembled qualitatively those of nonhuman primates (Puga-Gonzalez and Sueur 2017a; Puga-Gonzalez and Sueur 2017b). Hence, we had no specific predictions as to which model may fit empirical data better, except for the metric of modularity. From modularity, we know that the Reaper model fails at reproducing the variation observed in primate groups (Puga-Gonzalez and Sueur 2017a). Therefore, we expected the Reaper model to perform poorly regarding this metric. Overall, our main aim was to investigate which model reproduced quantitatively better the grooming network structure of the macaque groups analyzed. In this way, we aimed to know which model and thus which underlying mechanism is more likely to explain the diversity of social styles observed in societies of macaques. METHODS Macaque species From the primate literature, we surveyed studies on macaque species which reported a matrix of directed grooming, dominance hierarchy, and sampling method (scan or focal). In case of focal sampling, social grooming matrices were corrected for observation time, i.e., individual’s rows were divided by the number of times the individual was observed. In these cases thus, grooming matrices contained relative grooming frequency. Seven studies comprising 10 groups from 5 species met these criteria (Supplementary Table SI1). Further, coauthors of this study provided data from 4 groups from 4 species of macaques. Thus, our data collection comprised a total of 14 groups from 8 macaque species representing all of the 4 categorical grades of social styles attributed to macaque societies (Thierry 2004). Since in some of these studies data were available only for adult female individuals, we restricted our analyses to this category of individuals in all the groups. Note that due to the small number of adult females (n = 4) in the group of M. tonkeana (Orangerie Zoo), we included 2 subadult females (>4 years) in the analysis of this group. A full description of the groups included in this study is presented in Supplementary Table SI1. Data available from the Dryad Digital Repository (Puga-Gonzalez et al. 2018). Distribution of social grooming Whether social grooming was reciprocated, directed up the dominance hierarchy, or toward individuals of similar dominance rank was evaluated by means of the Tau-Kr matrix correlation (Hemelrijk 1990b). The level of significance was calculated using 2000 permutations and P-values < 0.05. To check for reciprocation of grooming, we correlated the groom-given matrix with the groom-received matrix. A significant positive correlation implies reciprocity. To check for grooming directed up the dominance hierarchy and to females of similar dominance rank, we tested the correlation between the grooming given matrix with the partner-rank matrix and the similar-rank matrix respectively. The partner-rank-matrix had the dominance rank of individuals in the rows, and the similar-rank-matrix was filled with zeros apart from the partners closest and second closest in dominance rank which were indicated as 1’s. A significantly positive correlation with the partner-rank-matrix corresponds to grooming being directed up the hierarchy, and a significantly positive correlation with the similar-rank-matrix corresponds to grooming directed to females of similar rank. Social networks metrics SNA were performed using R statistical software, version 3.2.2 (Team-R-Core 2015) and package igraph (Csardi and Nepusz 2006). Matrices of grooming given were transformed to directed graphs with the igraph package; then, the analysis of the directed grooming network was performed. Analyses were performed at the individual and at the group level. Individual level metrics We calculated 4 different metrics: betweenness, eigenvector, clustering coefficient, and closeness. These are metrics commonly used in analyses of animal societies and they inform us about the patterning of the individuals’ social interactions (Table 1). The technical definition of each metric is described below. Note that in some cases we used binary instead of weighted matrices to calculate some of the network metrics. We used binary matrices when the metric calculated was affected by the weight in the matrix (absolute grooming frequency) or when the metric can only be calculated from binary matrices. To convert matrices with absolute or relative grooming frequency to binary ones, we set all cell values > 0 to 1 and all others were set to 0. Metrics calculated from binary networks were betweenness centrality, closeness centrality, and clustering coefficient. To calculate eigenvector, we used weighted networks (with absolute or relative grooming frequency) because this metric can cope with weights. Unweighted betweenness centrality The number of times a node (i.e., an individual) acts as a bridge along the shortest path between 2 other nodes. Thus, it is a measure of how well an individual connects different subgroups or clusters within the network. Unweighted closeness centrality A measure of the degree to which a node is near all other nodes in a network. It is the average length of the shortest path between the node and all other nodes in the graph. In other words, this metric calculates how fast a node can reach all others. Clustering coefficient It is the number of closed-triangles a node form with its neighbors divided by the total possible number of closed-triangles that could exist among them. Thus, it assesses the degree to which nodes tend to cluster together. This measure takes into account neither direction nor weight. Weighted eigenvector It takes into account not only the number and strength of connections of the individual but also those of the partners to which it is connected. Since we used directed networks, here, high eigenvector means that an individual receives and has strong connections to other individuals who themselves receive grooming frequently. Group level metrics We chose 7 metrics that have been shown to be relevant in studies of primate social networks (Kasper and Voelkl 2009; Sueur, Petit, et al. 2011; Pasquaretta et al. 2014). All of these metrics are used to characterize the overall network structure, and each of them captures a different feature of the network. We used these metrics in order to get a detail comparison of the features the models capture better from empirical data and thus get a better assessment of which model (and mechanism) is most accurate at fitting the overall network structure. Number of triangles The absolute number of triangles in the network. Neither direction nor weights are taken into account. Unweighted diameter It is the longest path among all the shortest path lengths between 2 nodes, i.e., the shortest distance between the 2 most distant nodes in the network. Modularity The difference between the proportion of the total association of individuals within clusters (i.e., subgroups) and the expected proportion, given the summed associations of the different individuals. This metric is calculated using the eigenvector of each individual (Newman 2006). A high value of modularity means a high number of contacts within a subgroup, but few contacts between subgroups and low modularity means a homogeneous distribution of contacts between all group members. Mean clustering coefficient The proportion between the number of triangles in the network divided by total possible number of triangles. Neither direction nor weights are taken into account. Mean eigenvector The average of individual’s eigenvector coefficients per group. Centralization of the network It was calculated according to Equation 1, where the numerator, Cmax, is the highest eigenvector in the group and the denominator, Max, is the value obtained assuming the highest centralization possible, i.e., if the network were a star (Pasquaretta et al. 2014). Hence, the higher the value the more centralized is the network around one or several individuals. CIB=100*∑in(cmax−ci)Max∑in(cmax−ci) (1) Correlation dominance-eigenvector Spearman rank correlation coefficient between the individual’s eigenvector and the individual’s hierarchical rank. A positive correlation indicates a centralization of dominant individuals in the network. Agent-based computational models Here, we only present a brief description of each model. For a full description of the processes of the model, we refer to previous publications. GrooFiWorld model GrooFiWorld (Grooming and Fighting) is an individual-based model that is spatially explicit (Puga-Gonzalez et al. 2009). The model comprises a continuous 2-dimensional “world” (without borders) in which individuals are able to move in all directions. Individuals have a fixed vision angle (VisionAngle, Supplementary Table SI2) and a maximum perception distance (MaxView, Supplementary Table SI2). At the beginning of each simulation, the individuals’ location is assigned randomly within a previously defined radius (InitRadius, Supplementary Table SI2) calculated by multiplying group size by an arbitrary constant. To regulate individuals’ activations, each individual is attributed a random waiting time drawn from a uniform distribution; the individual with the shortest waiting time gets activated first. These waiting times are combined with a kind of social facilitation (Galef 1988) implying that an individual’s waiting time is reduced when a dominance interaction occurs close by (radius of social facilitation, Supplementary Table SI2). Intensity of aggression is reflected by the StepDom value (Supplementary Table SI2). Fierce aggression (e.g., bites), as in intolerant macaque societies, is represented by high values, and mild aggression (e.g., threats, slaps), as in tolerant societies, is represented by low values. To represent sexual dimorphism, females have a StepDom value 80% lower than that of males (Hemelrijk et al. 2008a). Grouping rules Why individuals group (to avoid predators, search resources, etc.) is not specified in the model. Whenever an individual does not see another close by (i.e., within its personal space, PersSpace, Supplementary Table SI2); then, it looks for others at increasing distances (NearView and MaxView, Supplementary Table SI2). If, even then none is in the field of vision, then the individual turns over a SearchAngle (Supplementary Table SI2) in order to look for others. In this way individuals tend to remain in a group. If, however, an individual “sees” another within its personal space (PersSpace, Supplementary Table SI2), then a social interaction may take place. Interactions Upon meeting another in its personal space, an individual first chooses whether or not to perform a dominance interaction. This choice depends on the risks involved, whereby risk concerns the chance of losing a fight (i.e., “Risk-aversion strategy”, Hemelrijk 2000). If an individual chooses not to fight; then, it considers grooming the other. If the individual chooses not to groom, no interaction happens and the individual stays put. Dominance interactions Dominance interactions are modeled as in Hemelrijk (1999) and are an extension of the DoDom rules of Hogeweg (1988). Each individual has a dominance value, Dom (Supplementary Table SI2), which represents its capacity to win. A dominance interaction takes place only if an individual expects to be victorious, i.e. individuals avoid risks. These risks are estimated by means of a “mental battle.” During a “mental battle”, an individual i compares its dominance value (Domi) relatively to that of its opponent j (Domj) to a random value between zero and one (Equation 2). This process is repeated RiskSens times (Supplementary Table SI2). In the current simulation RiskSens is set to 2, thus individuals have to win a mental battle twice before engaging in a real dominance interaction. wi=[1 DOMiDOMi+DOMj>RND(0,1)0 else (2) If in both “mental battles” individual i expects to be victorious, then a real dominance interaction occurs. The outcome of the real dominance interaction is again decided according to Equation 2. To reflect the self-reinforcing effects of victory and defeat (Barchas and Mendoza 1984; Hsu and Wolf 1999; Setchell et al. 2008; Hemelrijk et al. 2008a; Franz et al. 2015), dominance values are updated by increasing the dominance value of the winner and decreasing that of the loser by the same amount (Equation 3). This positive feedback is “dampened” because a victory of a higher-ranking opponent increases its relative Dom-value only slightly, whereas an (unexpected) success of the lower-ranking individual increases its relative dominance value by a greater change. To keep Dom-values positive, their minimum value is, arbitrarily, set at 0.01. DOMi=DOMi+(wi−DOMiDOMi+DOMj)*STEPDOMDOMj=DOMj−(wi−DOMjDOMi+DOMj)*STEPDOM (3) The change in Dom-values is multiplied by a scaling factor, StepDom, which represents intensity of aggression (Hemelrijk 1999). StepDom values range from 0.08 to 0.8 in females and from 0.1 to 1 in males. High values of StepDom cause a bigger change in Dom-values than low values of StepDom. After a fight, the winner chases the loser over a distance of one unit and the loser responds by fleeing under a small random angle over a predefined FleeingDistance (Supplementary Table SI2). To reduce consecutive interactions between same opponents, after the chase the winner turns randomly 45° to right or left. Grooming interactions When an individual meets another in its PersSpace and “chooses” not to fight, then, it considers grooming its partner. Grooming is induced by the level of Anxiety (Supplementary Table SI2), which ranges from very relaxed to very tense on a scale from 0 to 1. Individuals groom their partners if their level of Anxiety is higher than a random number between 0 and 1; otherwise, they display “nonaggressive” proximity. As indicated by empirical studies, grooming (giving and receiving) reduces anxiety and thus the tendency to groom (Schino et al. 1988; Aureli et al. 1999; Shutt et al. 2007; but see Molesti and Majolo 2013; Semple et al. 2013). It does so more strongly in the groomee (AnxDcrGree) than in the groomer (AnxDcrGrmr) (Supplementary Table SI2). Further, during periods without grooming, individuals increase their Anxiety with AnxInc (Supplementary Table SI2) and thus, their motivation to groom as demonstrated in empirical studies (Keverne et al. 1989; Martel et al. 1995; Graves et al. 2002). Furthermore, Anxiety increases in both opponents after a fight in the model (AnxIncFight, Supplementary Table SI2), as reported in nonhuman primates (Aureli et al. 2002; Butovskaya 2008). To avoid consecutive interactions between same partners, after grooming both partners turn randomly 45° to the right or left. The default grooming parameters were tuned in order to match the percentage of time grooming observed in primate societies ~20% (Dunbar 1991). Among macaque species, there is some variation regarding the percentage of time allocated to grooming; however given that the most accurate percentage of grooming time in each macaque species is unknown, we decided to tune all models to ~20%. Nevertheless, grooming patterns (reciprocation and interchange of grooming and support; grooming up the dominance hierarchy and among individuals of similar rank) emerge in the models as long as the percentage of time allocated to grooming is ≥5% (Puga-Gonzalez et al., unpublished data). FriendsWorld In contrast to GrooFiWorld, in FriendsWorld, individuals categorize others as “friends” or not. We classified as “friends” the top 5 partners with whom an individual groomed the most (given and received). We choose 5 individuals because this is the average number of grooming partners found in empirical studies irrespective of group size (Kudo and Dunbar 2001; Lehmann et al. 2007; Sueur, Deneubourg, et al. 2011). Note that the classification of others as “friends” or not is done every time an individual is activated, and it is based on the accumulation of grooming interactions until that point. Thus, individuals considered “friends” at one point during the simulation may not be the same at a point later in the simulation. During the period of time in which we collect data from the model, however, friends remain stable and change little (Puga-Gonzalez et al. 2015). Grouping rules FriendsWorld has the same architecture and behavioral rules as GrooFiWorld, the only exception is that in FriendsWorld model individuals follow their “friends” (Puga-Gonzalez et al. 2015). Individuals have 3 different visual ranges, PersSpace, NearView, and MaxView. When an individual is activated, it searches for others within its PersSpace, if it does not perceive another in its close proximity, it acts according to the grouping rules. If an individual perceives one of its “friends” within its NearView, it will move one step towards it. If several “friends” are perceived, the individual moves towards the closest one. If none of its “friends” is perceived but others are, the individual just keeps on moving without modifying its original direction. When no others are perceived within NearView, the individual looks further away into MaxView. If other individuals are perceived within MaxView, the individual moves towards the closest “friend” if available; otherwise, it moves towards the closest group-member. If no individual is perceived within MaxView, the individual scans for others by turning over a SearchAngle. Social interactions In FriendsWorld, the social interaction rules are the same as in GrooFiWorld. If the individual perceives another one within its PersSpace, a dominance interaction may occur. Whether or not an individual chooses to fight depends on the outcome of a mental battle (Equation 2). If the individual wins the mental battle, a dominance interaction occurs (see “Dominance interactions” above). However, if the individual loses the mental battle, it may groom its partner depending on its level of Anxiety (see “Grooming interactions” above). Reaper Model This model is not spatially explicit and the selection of interaction partners is based solely on previous grooming interactions with other group-members; therefore, as already mentioned, we called it Reaper, i.e., reinforcement of affiliative preferences (Puga-Gonzalez and Sueur 2017a). Selection of interaction partners In contrast to GrooFiWorld, in this model individuals are not located in space and thus proximity has no influence on the distribution of their social interactions. At the start of the simulation, individuals are attributed a random waiting time drawn from a uniform distribution and the individual with the shortest waiting time gets activated first. Once activated, the probability of selecting a given group member as an interaction partner depends on the amount of grooming given and received exchanged with that specific partner. This probability is calculated according to Equation 4. Where Pij is the probability of individual i interacting with partner j; Gij is the number of grooming bouts exchanged between individual i and j; and the denominator is the sum of the number of grooming bouts individual i has exchanged with every other group member. α determines the degree of nonlinearity in the probability of selecting a given partner, the higher the value of α, the higher the tendency of individuals to interact with their most frequent grooming partners. This value of was set to 1.2. We choose this value of α because with it individuals allocated >75% of their grooming time to their top 5 interaction partners, as commonly observed in primate societies (Sueur, Deneubourg, et al. 2011)—for a full sensitivity analysis of the values of α, see Puga-Gonzalez and Sueur (2017a). At the beginning of the simulation, individuals have not yet interacted with any other group member and thus, all individuals have the same probability of being chosen as an interaction partner. As interactions go by, however, individuals will tend to select more frequently as interaction partners those individuals with whom they are more involved in social grooming. Note that Equation 4 does not imply reciprocal relationships, i.e., if individual i prefers to interact mostly with individual j, this does not necessarily mean that individual j also prefers to interact with individual of i. After selecting an interaction partner, individuals follow the same interactions rules as in GrooFiWorld. They first choose on the basis of a “mental fight” whether they perform or not a dominance interaction, if they choose not to fight, then they decide whether or not to groom (see dominance and grooming interactions above). Pij=Gijα∑a≠iN−1Giaα (4) Null-model In order to compare our results to a “null-model” model, we built a model in which we omit the effects of spatial proximity and preferential interactions with specific group members. In this model, individuals select interaction partners at random. Once an individual has randomly selected an interaction partner, a random number is drawn from a uniform distribution between [0, 1]. If this number is ≤ 0.5 the individual performs a dominance interaction, otherwise it grooms. Note that in this case individuals do not perform a “mental fight” to decide whether or not to attack, if the random number is ≤ to 0.5 individuals immediately attack. If the random number is > 0.5 individuals groom their partner independently of its value of anxiety. As in GrooFiWorld, individuals’ activations are regulated by attributing each individual a random waiting time drawn from a uniform distribution. The individual with the shortest waiting time gets activated first. Parametrization of the models In order to compare the empirical grooming networks with those of the models, in the simulations we mimicked the same group size and densities of each of the empirical networks. We did so because comparisons between networks differing in size and density (the number of observed edges divided by the number of possible edges) appear to be meaningless (Rankin et al. 2016). Indeed many network measures are dependent on these 2 parameters (Kasper and Voelkl 2009; Sueur, Petit, et al. 2011; Pasquaretta et al. 2014). Group size consisted of all adult individuals in the group; sex of individuals was also taken into consideration. Note that even though we only analyzed the behavior of females, we also included males in the simulations since male behavior may influence the behavior of real (Hemelrijk et al. 2008b) and computer simulated females (Hemelrijk et al. 2008a; Hemelrijk and Puga-Gonzalez 2012). In this way, we compared networks of females taken into consideration the presence of males in the group. To obtain female grooming networks with the same density as that observed in empirical data, first we collected data for long periods of time (up to 10,000 grooming interactions in some cases). Then, we constructed matrices in which the number of grooming interactions was increased sequentially following the order of data collection until the target density (± 2.5%) of the network was matched. Further, we also mimicked the grade of social style of each of the species of macaques. From previous analyses of the models, we knew that values of intensity of aggression of 1.0, 0.8, and 0.1, 0.08 for males and females, respectively, reproduce patterns that resemble those of intolerant and tolerant macaques, respectively (Hemelrijk and Puga-Gonzalez 2012; Puga-Gonzalez et al. 2015; Puga-Gonzalez and Sueur 2017a). Hence values for grades 1 and 4 of social style were set to 1.0, 0.8, and 0.1, 0.08 for males and females, respectively. Then, the values for grades 2 and 3 were chosen so they laid in between grades 1 and 4. In this way, we created 4 different sets of values of intensity of aggression “mimicking” the 4 different grades of social style described by Thierry (2004). These values were: grade 4, 0.1, 0.08; grade 3, 0.4, 0.32; grade 2, 0.7, 0.56; and grade one, 1.0, 0.8; for males and females respectively. The values of intensity of aggression assigned to each of the different groups are shown in Supplementary Table SI3. Preliminary analysis and optimization of the distribution of grooming in the models After running the models with the values of intensity of aggression described in the section “Parametrization of the models” (Supplementary Table SI3), we realized that the models did not match the presence/absence of some grooming patterns (grooming directed up the dominance hierarchy and/or towards individuals of similar rank) observed in the empirical groups. The presence/absence of these patterns may influence the value of some networks’ metrics. For instance, whether or not individuals direct grooming up the dominance hierarchy affects the correlation dominance rank—eigenvector and overall network centralization. Therefore, we decided to optimize the fit between the models’ grooming patterns and those observed in each group of empirical data. To optimize these patterns, we adjusted the parameter intensity of aggression (StepDom, Supplementary Table SI2). We specifically focused on this parameter because we knew it affects how grooming is distributed across group members whereas the other model’s parameters mainly affect cohesiveness of the group and/or the overall rate of grooming (Hemelrijk et al. 2017). The higher the value of intensity of aggression, the steeper the dominance hierarchy, the more individuals direct grooming up the hierarchy, at partners of similar rank and the more pronounced the spatial structure with dominant individuals at the center (Puga-Gonzalez et al. 2009). Hence, to match the presence/absence of these grooming patterns, we increased/decreased the values of intensity of aggression in the models. The new values of intensity of aggression assigned to each model and simulated group are shown in Supplementary Table SI8. Note that here we only show the results after the optimization of the grooming patterns. The results without optimization are found in the supporting information (see text therein, Supplementary Tables SI5–SI7; Supplementary SI12–SI14, and Supplementary Figure SI1–SI4). Data collection of the models For each empirical group (n = 14) and model (n = 4), we ran 20 different replicates, thus we ran a total of 1120 simulation. Each replicate consisted of up to 300 periods and each period consisted of (Group_Size * 20) individual’s activations. Data were collected from period 200 onwards to exclude any bias caused by transient values (Hemelrijk 1999). Data collection consisted of every social interaction: dyadic fights and grooming behavior. We recorded, for dominance interactions the identities of 1) the attacker and its opponent, 2) that of the winner/loser, 3) the updated Dom values of the individuals; and for grooming interactions, the identities of 4) groomer, 5) groomee, and 6) the updated anxiety values. From these interactions, we built grooming and aggression matrices. Note that the grooming matrices matched approximately the same density as that observed in empirical groups. Analyses of socio-behavioral patterns in the models For reasons of comparison with empirical data, the analyses of social behavior were performed only among females. As in empirical data we checked for grooming reciprocation, grooming up the hierarchy, and grooming towards partners of similar rank by means of the matrix Tau-Kr correlation (Hemelrijk 1990b). The level of significance was calculated using 2000 permutations and P-value < 0.05. The results shown are the average of 20 replicates with their combined probability using the improved Bonferroni procedure (Hochberg 1988). Dominance hierarchies were stable, and dominance ranks were significantly correlated (Kendall correlation every 30 periods, e.g., 200–230: tau > 0.88**). The asymmetry of dominance ranks was measured by the coefficient of variation of Dominance values (standard deviation divided by the mean). As in our previous studies, the higher the intensity of aggression the higher the asymmetry of dominance ranks (Hemelrijk and Puga-Gonzalez 2012; Puga-Gonzalez et al. 2015; Puga-Gonzalez and Sueur 2017a). Analysis of social networks of the models Like the analyses of social behavior, analyses were performed only among females. We performed the same social networks analysis as that performed in the social networks of the empirical groups (see above for metrics and definitions, sections “Individual level metrics” and “Group level metrics”). Note, however, that in the social networks derived from the models it was not necessary to correct for observation time since we kept track of all individuals’ interactions. Thus, network analyses were performed on directed binary/weighted social networks, depending on the network metric analyzed (see sections “Individual level metrics” and “Group level metrics”). In weighted networks, the weights represented the absolute grooming frequency given to another individual. Results of the global metrics shown are values of each of the 20 replicates per macaque group simulated. Likewise, results of the individuals’ metrics shown are values per individual in each of the 20 replicates per macaque group simulated. Thus, for each empirical data point (global or individual metric) there were 20 data points from each model. Comparison of social networks and model selection We performed linear regressions with the values of the observed metrics as the response variable, and the values of the metrics obtained from the simulated networks as the predictor variable. At the individual level, values were paired by sorting individuals by rank (within each group); at the global level, they were paired by groups. We compared effect sizes (R-squared), i.e., the percentage of variation explained by each model. However, to select which model fitted the observed data at best, we used the Akaike Information Criteria (AIC). We did so because we reasoned that if the social networks metrics from the models had a perfect fit with the observed metrics, then the observed and simulated values should be identical. In other words, if we plotted the values of the observed metrics (x axis) against the values from the model (y axis) and there was a perfect correspondence, then we would expect the values to lie along the diagonal with slope = 1 and intercept = 0 (note that the slope of the linear regressions was not always close to one). The AIC was calculated using the residual sum of squares (RSS) from this expectation according to Equation 5. AIC=n*ln(RSSn)+2*K (5) Where n is the number of data points, ln is the natural logarithm, RSS the residual sum of squares, and K the number of parameters in the model which was set to 2 (the model + error). The best model was determined by calculating the difference (Δ) between the model with the lowest AIC and all the other models, and by calculating the Akaike weight for each model (Equation 6). Where wi is the Akaike weight for model i, Δi is the difference between the AIC of the best fitting model and that of model i, and the denominator is the sum of the relative likelihood for all candidate models. Models with Δ values ≤2 were considered to have the same explanatory power. wi=exp(−0.5*Δi)∑r=1Rexp(−0.5*Δr) (6) RESULTS Grooming interactions and networks’ density The number of interactions necessary to match the target empirical density varied between the different models -GrooFiWorld, FriendsWorld, Reaper, and Null-model (Supplementary Table SI4). The null-model required the least number of grooming interactions to match the empirical network density. As regards the other 3 models, no apparent difference among them was found when the target density was below 60% and intensity of aggression was set to low (Supplementary Table SI4). However, for densities above 60% the Reaper model took many more grooming interactions than the two spatially explicit models to reach the target density. More interestingly, in all models except the null one, simulations belonging to some specific groups never reached a network density of 100% regardless of the number of grooming interactions. In these cases, the density reached a plateau around values of ~80%. In the spatially explicit models (GrooFiWorld and FriendsWorld), this only happened when intensity of aggression was set to too high values. When intensity of aggression was high the asymmetry of the dominance values among females was great and thus, when females at the top of the hierarchy interacted with those at the bottom, they chose to fight instead of grooming. In the Reaper model, however, this happened regardless of the values of intensity of aggression. In this case, females became so selective on their interaction partners that they never selected other group members as interaction partners. Patterns of grooming Reciprocation of social grooming In all groups of empirical data, grooming given was positively correlated with grooming received except in the 2 groups of M. tonkeana in which no significant correlation was found (Supplementary Table SI5). Concerning the computer models, regardless of the simulated group, we always found a significant positive matrix correlation between grooming given and received, except in the null-model. Correlations were never significant in this model (Supplementary Table SI5). Optimization of distribution of grooming We modified the values of intensity of aggression (StepDom) to optimize the match between the patterns emerging in the model (grooming directed up the dominance hierarchy and towards individuals of similar rank) and those observed in the groups of empirical data (see Methods and Supplementary Table SI8). Here we only present the results after optimization, for the results without optimization see supporting information (text therein, Supplementary Tables SI5–SI7; Supplementary Tables SI12–SI14, and Supplementary Figures SI1–SI4). Regarding grooming directed up the dominance hierarchy, after optimization of the parameters (Supplementary Table SI8), GrooFiWorld matched all observations in empirical data (14/14); Reaper model matched all except one group (13/14); and FriendsWorld did not match 3 groups (11/14); in the null-model, grooming up the dominance hierarchy never reached statistical significance (Supplementary Table S10). Regarding grooming directed towards individuals of similar rank, the optimization of parameters improved the match between models and empirical data, but there were still several mismatches. GrooFiWorld did not match 4 groups (10/14); FriendsWorld 5 groups (9/14); and the Reaper 7 groups (7/14). Note that in the Reaper and null models this pattern never reached significance; thus the cases that matched empirical data were those in which this pattern was absent (Supplementary Table S11). The reason for the mismatches is that in GrooFiWorld and FriendsWorld grooming up the hierarchy and directed towards partners of similar rank usually emerge (or not) together. In empirical data, however, this was not always the case. In the group of M. assamensis and M. radiata, grooming was directed up the hierarchy but not towards partners of similar rank. Conversely, in 2 groups of M. fuscata and one of M. fascicularis, grooming was directed towards partners of similar rank but not up the dominance hierarchy (Supplementary Tables S10–11). SNA Because the condition in which animals lived, captive or wild (free-ranging), can affect the way individuals interact with each other, we decided to compare the network metrics by dividing the empirical groups into wild/provisioned/free-ranging (n = 8) and captive (n = 6) (Supplementary Table SI1). Condition seemed to have an effect on whether the model matched or not the empirical data, particularly for the global metrics of modularity, eigenvector, correlation dominance—eigenvector, and centralization index (Supplementary Figures SI9–SI12). All models matched better these metrics when the simulated group was wild rather than captive (Supplementary Figures SI11–SI12). This is evident from the much higher R2 in wild than in the captive groups (Supplementary Table SI18) and from the negative slopes of the regression lines for the captive groups, i.e., opposite to expected (Supplementary Figure SI12). Thus, we only present here the comparison between the network metrics of the models with those of wild groups of macaques. For the network metric comparison without optimization of the grooming patterns, and for the comparison after optimizing the grooming patterns but without dividing groups into wild and captive, see supporting information. Individual metrics With regard to the individual network metrics of closeness, betweenness, and eigenvector all computer models seemed to predict the observed female’s values, particularly those of closeness (Figure 1). Effect sizes were similar among all models, explaining >30% of the variance of observed values of clustering coefficient; >60% of the variance of betweenness and eigenvector; and > 97% of the variance in values of closeness (Table 2). FriendsWorld had the highest effect size for betweenness; FriendsWorld and Reaper for closeness; and Reaper for eigenvector and clustering coefficient (Table 2). After calculating the Akaike information criteria (AIC), it appeared that patterns generated in FriendsWorld fitted best the empirical values of betweenness and closeness, and Reaper those of eigenvector and clustering coefficient (Table 2). In all cases, the null-model was the worst fitting model (Table 2). Figure 1 View largeDownload slide Plots of the values of the observed individual’s metric (x axis) against the values predicted by each theoretical model (y axis). Black diagonal represents values expected in a perfect match between model and empirical data. Values of betweenness centrality and closeness are shown in left and right columns, respectively. Colors represent macaque groups and each data point represents a single individual. Figure 1 View largeDownload slide Plots of the values of the observed individual’s metric (x axis) against the values predicted by each theoretical model (y axis). Black diagonal represents values expected in a perfect match between model and empirical data. Values of betweenness centrality and closeness are shown in left and right columns, respectively. Colors represent macaque groups and each data point represents a single individual. Table 2 Linear regressions between observed and predicted individuals’ network metric for each model, and its correspondent Akaike information criteria (AIC) (for details of the calculation, see Methods) Linear regression Akaike information Metric Model β R2 P-value AIC Delta (Δ) Weight Closeness FriendsWorld 0.97 0.98 <0.001 −21,120 0 1 Reaper 0.99 0.98 <0.001 −20,792 327 0 GrooFiWorld 0.94 0.97 <0.001 −20,131 989 0 Null 0.97 0.97 <0.001 −19,509 1611 0 Betweenness FriendsWorld 1.38 0.73 <0.001 10,090 0 1 Reaper 1.73 0.68 <0.001 10,730 640 0 GrooFiWorld 1.47 0.63 <0.001 10,745 655 0 Null 1.95 0.65 <0.001 10,989 899 0 Eigen Centrality Reaper 0.92 0.69 <0.001 −6611 0 1 FriendsWorld 0.80 0.65 <0.001 −5996 616 0 GrooFiWorld 0.77 0.61 <0.001 −5820 791 0 Null 1.11 0.60 <0.001 −4645 1966 0 Clustering Coeff. Reaper 0.58 0.33 <0.001 −8761 0 1 GrooFiWorld 0.42 0.29 <0.001 −8178 583 0 FriendsWorld 0.44 0.31 <0.001 −7875 886 0 Null 0.30 0.09 <0.001 −7013 1748 0 Linear regression Akaike information Metric Model β R2 P-value AIC Delta (Δ) Weight Closeness FriendsWorld 0.97 0.98 <0.001 −21,120 0 1 Reaper 0.99 0.98 <0.001 −20,792 327 0 GrooFiWorld 0.94 0.97 <0.001 −20,131 989 0 Null 0.97 0.97 <0.001 −19,509 1611 0 Betweenness FriendsWorld 1.38 0.73 <0.001 10,090 0 1 Reaper 1.73 0.68 <0.001 10,730 640 0 GrooFiWorld 1.47 0.63 <0.001 10,745 655 0 Null 1.95 0.65 <0.001 10,989 899 0 Eigen Centrality Reaper 0.92 0.69 <0.001 −6611 0 1 FriendsWorld 0.80 0.65 <0.001 −5996 616 0 GrooFiWorld 0.77 0.61 <0.001 −5820 791 0 Null 1.11 0.60 <0.001 −4645 1966 0 Clustering Coeff. Reaper 0.58 0.33 <0.001 −8761 0 1 GrooFiWorld 0.42 0.29 <0.001 −8178 583 0 FriendsWorld 0.44 0.31 <0.001 −7875 886 0 Null 0.30 0.09 <0.001 −7013 1748 0 The best fitting models according to the AIC weight are shown in bold. Note that for each metric the models are ordered according to increasing values of delta (Δ). View Large Table 2 Linear regressions between observed and predicted individuals’ network metric for each model, and its correspondent Akaike information criteria (AIC) (for details of the calculation, see Methods) Linear regression Akaike information Metric Model β R2 P-value AIC Delta (Δ) Weight Closeness FriendsWorld 0.97 0.98 <0.001 −21,120 0 1 Reaper 0.99 0.98 <0.001 −20,792 327 0 GrooFiWorld 0.94 0.97 <0.001 −20,131 989 0 Null 0.97 0.97 <0.001 −19,509 1611 0 Betweenness FriendsWorld 1.38 0.73 <0.001 10,090 0 1 Reaper 1.73 0.68 <0.001 10,730 640 0 GrooFiWorld 1.47 0.63 <0.001 10,745 655 0 Null 1.95 0.65 <0.001 10,989 899 0 Eigen Centrality Reaper 0.92 0.69 <0.001 −6611 0 1 FriendsWorld 0.80 0.65 <0.001 −5996 616 0 GrooFiWorld 0.77 0.61 <0.001 −5820 791 0 Null 1.11 0.60 <0.001 −4645 1966 0 Clustering Coeff. Reaper 0.58 0.33 <0.001 −8761 0 1 GrooFiWorld 0.42 0.29 <0.001 −8178 583 0 FriendsWorld 0.44 0.31 <0.001 −7875 886 0 Null 0.30 0.09 <0.001 −7013 1748 0 Linear regression Akaike information Metric Model β R2 P-value AIC Delta (Δ) Weight Closeness FriendsWorld 0.97 0.98 <0.001 −21,120 0 1 Reaper 0.99 0.98 <0.001 −20,792 327 0 GrooFiWorld 0.94 0.97 <0.001 −20,131 989 0 Null 0.97 0.97 <0.001 −19,509 1611 0 Betweenness FriendsWorld 1.38 0.73 <0.001 10,090 0 1 Reaper 1.73 0.68 <0.001 10,730 640 0 GrooFiWorld 1.47 0.63 <0.001 10,745 655 0 Null 1.95 0.65 <0.001 10,989 899 0 Eigen Centrality Reaper 0.92 0.69 <0.001 −6611 0 1 FriendsWorld 0.80 0.65 <0.001 −5996 616 0 GrooFiWorld 0.77 0.61 <0.001 −5820 791 0 Null 1.11 0.60 <0.001 −4645 1966 0 Clustering Coeff. Reaper 0.58 0.33 <0.001 −8761 0 1 GrooFiWorld 0.42 0.29 <0.001 −8178 583 0 FriendsWorld 0.44 0.31 <0.001 −7875 886 0 Null 0.30 0.09 <0.001 −7013 1748 0 The best fitting models according to the AIC weight are shown in bold. Note that for each metric the models are ordered according to increasing values of delta (Δ). View Large Global metrics All models seemed to accurately predict the observed values of number of triangles and clustering coefficient (Figure 3). Effect sizes were highest for GrooFiWorld (R2 = 0.82) for number of triangles, and Reaper and FriendsWorld (R2 = 0.89) for clustering coefficient (Table 3). According to the AIC, GroFiWorld and Reaper were the best fitting models for number of triangles and clustering coefficient, respectively (Table 3). Regarding diameter and the correlation dominance—eigenvector, all models, except the null one, moderately predicted the observed values (Figures 3 and 4). Effect sizes were highest for FriendsWorld (R2 = 0.49) for diameter and GrooFiWorld (R2 = 0.43) for the correlation dominance—eigenvector, and according to the AIC, these models were also the best fitting ones (Table 3). With respect to modularity of the network, the match between observed and predicted data was not as good as for the above metrics; and in the case of Reaper model, modularity seemed to not vary but to remain constant (Figure 4). Similarly, the values of centrality index were poorly matched by all models (Figure 4). Effect sizes were highest for GrooFiWorld (R2 = 0.34) and Reaper model (R2 = 0.36) for modularity and centrality index respectively (Table 3). And according to the AIC, the best fitting model were also GrooFiWorld and Reaper (Table 3). The average eigenvector value of the group was the network metric matched at worst by the models, all of them overestimated this value (Figure 4). Effect size was highest for FriendsWorld (R2 = 0.25), but according to the AIC, the Reaper model seemed to fit the data at best (Table 3). As for the individual metrics, in all cases, except for clustering coefficient, the null model was the worst fitting model (Table 3). Table 3 Linear regression between observed and predicted global network metric for each model, and its correspondent Akaike information (AIC) (for details of the calculation, see methods) Linear Regression Akaike information Metric Model Β R2 P-value AIC Delta (Δ) Weight Num Triangles GrooFiWorld 0.92 0.82 <0.001 1085 0 1 Reaper 0.82 0.81 <0.001 1124 39 0 FriendsWorld 0.81 0.81 <0.001 1140 55 0 Null 0.62 0.73 <0.001 1308 223 0 Diameter FriendsWorld 0.86 0.49 <0.001 52 0 1 Reaper 0.95 0.41 <0.001 77 25 0 GrooFiWorld 0.82 0.41 <0.001 82 30 0 Null 1.07 0.34 <0.001 105 53 0 Clustering Coeff. Reaper 1.00 0.89 <0.001 −780 0 1 FriendsWorld 1.28 0.89 <0.001 −657 123 0 Null 1.07 0.83 <0.001 −652 129 0 GrooFiWorld 1.17 0.81 <0.001 −636 144 0 Modularity GrooFiWorld 0.91 0.34 <0.001 −508 0 1 FriendsWorld 0.93 0.28 <0.001 −491 17 0 Reaper 0.13 0.00 0.627 −482 26 0 Null 1.74 0.28 <0.001 −389 118 0 Eigenvector Reaper 0.58 0.22 <0.001 −588 0 1 FriendsWorld 0.50 0.25 <0.001 −519 69 0 GrooFiWorld 0.36 0.14 <0.001 −496 93 0 Null 0.47 0.16 <0.001 −384 205 0 Corr Dom-Eigen GrooFiWorld 0.58 0.43 <0.001 −369 0 1 FriendsWorld 0.59 0.39 <0.001 −352 17 0 Reaper 0.44 0.17 <0.001 −349 20 0 Null 0.02 −0.01 0.779 −163 206 0 Centrality Index Reaper 0.86 0.36 <0.001 1057 0 1 FriendsWorld 0.74 0.26 <0.001 1133 76 0 GrooFiWorld 0.45 0.10 <0.001 1162 105 0 Null 0.60 0.15 <0.001 1263 206 0 Linear Regression Akaike information Metric Model Β R2 P-value AIC Delta (Δ) Weight Num Triangles GrooFiWorld 0.92 0.82 <0.001 1085 0 1 Reaper 0.82 0.81 <0.001 1124 39 0 FriendsWorld 0.81 0.81 <0.001 1140 55 0 Null 0.62 0.73 <0.001 1308 223 0 Diameter FriendsWorld 0.86 0.49 <0.001 52 0 1 Reaper 0.95 0.41 <0.001 77 25 0 GrooFiWorld 0.82 0.41 <0.001 82 30 0 Null 1.07 0.34 <0.001 105 53 0 Clustering Coeff. Reaper 1.00 0.89 <0.001 −780 0 1 FriendsWorld 1.28 0.89 <0.001 −657 123 0 Null 1.07 0.83 <0.001 −652 129 0 GrooFiWorld 1.17 0.81 <0.001 −636 144 0 Modularity GrooFiWorld 0.91 0.34 <0.001 −508 0 1 FriendsWorld 0.93 0.28 <0.001 −491 17 0 Reaper 0.13 0.00 0.627 −482 26 0 Null 1.74 0.28 <0.001 −389 118 0 Eigenvector Reaper 0.58 0.22 <0.001 −588 0 1 FriendsWorld 0.50 0.25 <0.001 −519 69 0 GrooFiWorld 0.36 0.14 <0.001 −496 93 0 Null 0.47 0.16 <0.001 −384 205 0 Corr Dom-Eigen GrooFiWorld 0.58 0.43 <0.001 −369 0 1 FriendsWorld 0.59 0.39 <0.001 −352 17 0 Reaper 0.44 0.17 <0.001 −349 20 0 Null 0.02 −0.01 0.779 −163 206 0 Centrality Index Reaper 0.86 0.36 <0.001 1057 0 1 FriendsWorld 0.74 0.26 <0.001 1133 76 0 GrooFiWorld 0.45 0.10 <0.001 1162 105 0 Null 0.60 0.15 <0.001 1263 206 0 The best fitting models according to the AIC weight are shown in bold. Note that for each metric the models are ordered according to increasing values of delta (Δ). View Large Table 3 Linear regression between observed and predicted global network metric for each model, and its correspondent Akaike information (AIC) (for details of the calculation, see methods) Linear Regression Akaike information Metric Model Β R2 P-value AIC Delta (Δ) Weight Num Triangles GrooFiWorld 0.92 0.82 <0.001 1085 0 1 Reaper 0.82 0.81 <0.001 1124 39 0 FriendsWorld 0.81 0.81 <0.001 1140 55 0 Null 0.62 0.73 <0.001 1308 223 0 Diameter FriendsWorld 0.86 0.49 <0.001 52 0 1 Reaper 0.95 0.41 <0.001 77 25 0 GrooFiWorld 0.82 0.41 <0.001 82 30 0 Null 1.07 0.34 <0.001 105 53 0 Clustering Coeff. Reaper 1.00 0.89 <0.001 −780 0 1 FriendsWorld 1.28 0.89 <0.001 −657 123 0 Null 1.07 0.83 <0.001 −652 129 0 GrooFiWorld 1.17 0.81 <0.001 −636 144 0 Modularity GrooFiWorld 0.91 0.34 <0.001 −508 0 1 FriendsWorld 0.93 0.28 <0.001 −491 17 0 Reaper 0.13 0.00 0.627 −482 26 0 Null 1.74 0.28 <0.001 −389 118 0 Eigenvector Reaper 0.58 0.22 <0.001 −588 0 1 FriendsWorld 0.50 0.25 <0.001 −519 69 0 GrooFiWorld 0.36 0.14 <0.001 −496 93 0 Null 0.47 0.16 <0.001 −384 205 0 Corr Dom-Eigen GrooFiWorld 0.58 0.43 <0.001 −369 0 1 FriendsWorld 0.59 0.39 <0.001 −352 17 0 Reaper 0.44 0.17 <0.001 −349 20 0 Null 0.02 −0.01 0.779 −163 206 0 Centrality Index Reaper 0.86 0.36 <0.001 1057 0 1 FriendsWorld 0.74 0.26 <0.001 1133 76 0 GrooFiWorld 0.45 0.10 <0.001 1162 105 0 Null 0.60 0.15 <0.001 1263 206 0 Linear Regression Akaike information Metric Model Β R2 P-value AIC Delta (Δ) Weight Num Triangles GrooFiWorld 0.92 0.82 <0.001 1085 0 1 Reaper 0.82 0.81 <0.001 1124 39 0 FriendsWorld 0.81 0.81 <0.001 1140 55 0 Null 0.62 0.73 <0.001 1308 223 0 Diameter FriendsWorld 0.86 0.49 <0.001 52 0 1 Reaper 0.95 0.41 <0.001 77 25 0 GrooFiWorld 0.82 0.41 <0.001 82 30 0 Null 1.07 0.34 <0.001 105 53 0 Clustering Coeff. Reaper 1.00 0.89 <0.001 −780 0 1 FriendsWorld 1.28 0.89 <0.001 −657 123 0 Null 1.07 0.83 <0.001 −652 129 0 GrooFiWorld 1.17 0.81 <0.001 −636 144 0 Modularity GrooFiWorld 0.91 0.34 <0.001 −508 0 1 FriendsWorld 0.93 0.28 <0.001 −491 17 0 Reaper 0.13 0.00 0.627 −482 26 0 Null 1.74 0.28 <0.001 −389 118 0 Eigenvector Reaper 0.58 0.22 <0.001 −588 0 1 FriendsWorld 0.50 0.25 <0.001 −519 69 0 GrooFiWorld 0.36 0.14 <0.001 −496 93 0 Null 0.47 0.16 <0.001 −384 205 0 Corr Dom-Eigen GrooFiWorld 0.58 0.43 <0.001 −369 0 1 FriendsWorld 0.59 0.39 <0.001 −352 17 0 Reaper 0.44 0.17 <0.001 −349 20 0 Null 0.02 −0.01 0.779 −163 206 0 Centrality Index Reaper 0.86 0.36 <0.001 1057 0 1 FriendsWorld 0.74 0.26 <0.001 1133 76 0 GrooFiWorld 0.45 0.10 <0.001 1162 105 0 Null 0.60 0.15 <0.001 1263 206 0 The best fitting models according to the AIC weight are shown in bold. Note that for each metric the models are ordered according to increasing values of delta (Δ). View Large DISCUSSION On basis of the effect sizes and the AIC, we could not distinguish one model (GrooFiWorld, FriendsWorld, or Reaper) fitting all network metrics better than the others. Each of these models fitted the empirical data in at least 3 network metrics better than the others (Table 4). The null model, on the other hand, performed worst for all metrics except one (Table 4). Patterns generated by Reaper and FriendsWorld seemed to fit the empirical network metrics slightly better than GrooFiWorld which had the worst fit among the non-null models in 7 out of 11 cases (Table 4). Regarding grooming reciprocation, directed up the dominance hierarchy, and towards individuals of similar rank, GrooFiWorld performed better, followed by FriendsWorld and Reaper models (Supplementary Tables SI9–SI11). Similarly, when the analysis were performed without optimization and/or without dividing the data into wild and captive, the overall ranks of the models changed slightly and no model fitted all patterns better than the others (Supplementary Tables SI14 and SI17). Our analysis thus does not allow to establish that one model was better than other. If anything, it indicates that both space, i.e., proximity-based reciprocity as modeled in GrooFiWorld and FriendsWorld, and social bonding, i.e., emotional bookkeeping as modeled in Reaper, are relevant factors influencing the distribution of social grooming. Table 4 Rank of the models according to the AIC weight and Delta (Δ) value for each network metric Computer Model Network metric FriendsWorld Groofiworld Reaper Null-model Individual: 1) Closeness 1 3 2 4 2) Betweenness 1 3 2 4 3) Eigenvector 2 3 1 4 4) Clustering coefficient 2.5 2.5 1 4 Global: 5) Number of triangles 3 1 2 4 6) Diameter 1 3 2 4 7) Average clustering coefficient 2 4 1 3 8) Modularity 2 1 3 4 9) Average eigenvector centrality 2 3 1 4 10) Corr. eigenvector - dominance rank 2 1 3 4 11) Centrality Index 2 3 1 4 Average Ranking 2.0 2.5 1.7 3.9 Computer Model Network metric FriendsWorld Groofiworld Reaper Null-model Individual: 1) Closeness 1 3 2 4 2) Betweenness 1 3 2 4 3) Eigenvector 2 3 1 4 4) Clustering coefficient 2.5 2.5 1 4 Global: 5) Number of triangles 3 1 2 4 6) Diameter 1 3 2 4 7) Average clustering coefficient 2 4 1 3 8) Modularity 2 1 3 4 9) Average eigenvector centrality 2 3 1 4 10) Corr. eigenvector - dominance rank 2 1 3 4 11) Centrality Index 2 3 1 4 Average Ranking 2.0 2.5 1.7 3.9 Models with a Delta (Δ) value difference ≤ 2 were given the same rank. View Large Table 4 Rank of the models according to the AIC weight and Delta (Δ) value for each network metric Computer Model Network metric FriendsWorld Groofiworld Reaper Null-model Individual: 1) Closeness 1 3 2 4 2) Betweenness 1 3 2 4 3) Eigenvector 2 3 1 4 4) Clustering coefficient 2.5 2.5 1 4 Global: 5) Number of triangles 3 1 2 4 6) Diameter 1 3 2 4 7) Average clustering coefficient 2 4 1 3 8) Modularity 2 1 3 4 9) Average eigenvector centrality 2 3 1 4 10) Corr. eigenvector - dominance rank 2 1 3 4 11) Centrality Index 2 3 1 4 Average Ranking 2.0 2.5 1.7 3.9 Computer Model Network metric FriendsWorld Groofiworld Reaper Null-model Individual: 1) Closeness 1 3 2 4 2) Betweenness 1 3 2 4 3) Eigenvector 2 3 1 4 4) Clustering coefficient 2.5 2.5 1 4 Global: 5) Number of triangles 3 1 2 4 6) Diameter 1 3 2 4 7) Average clustering coefficient 2 4 1 3 8) Modularity 2 1 3 4 9) Average eigenvector centrality 2 3 1 4 10) Corr. eigenvector - dominance rank 2 1 3 4 11) Centrality Index 2 3 1 4 Average Ranking 2.0 2.5 1.7 3.9 Models with a Delta (Δ) value difference ≤ 2 were given the same rank. View Large All models fitted relatively well the observed values of the individual metrics of betweenness, closeness, and eigenvector. In these cases, effect sizes were usually high and data points laid usually along the diagonal (Figures 1 and 2, Table 2). This was the case even for the null-model, which indicates that the value of these individual metrics may be partially constrained by group size and network density (i.e., the network properties matched in the simulations). The models did not explain more than 33% of the variance of clustering coefficient. This indicates that the grooming cliques formed by individuals were different in the model and empirical data. In the models, clustering coefficient was mostly overestimated, i.e., individuals and their neighbors were more connected than in reality; suggesting then, that in reality individuals are more selective in choosing their interaction partners. Figure 2 View largeDownload slide Plots of the values of the observed individual’s metric (x axis) against the values predicted by each model (y axis). Black diagonal represents values expected in a perfect match between model and empirical data. Values of clustering coefficient and eigenvector are shown in left and right columns, respectively. Colors represent macaque groups and each data point represents a single individual. Figure 2 View largeDownload slide Plots of the values of the observed individual’s metric (x axis) against the values predicted by each model (y axis). Black diagonal represents values expected in a perfect match between model and empirical data. Values of clustering coefficient and eigenvector are shown in left and right columns, respectively. Colors represent macaque groups and each data point represents a single individual. With respect to the global metrics, none of the models performed as good as for the individual metrics (Figures 3 and 4). All models seemed to overestimate the values of eigenvector centrality and of the correlation coefficient dominance-eigenvector; and underestimate those of the centrality index (Figure 4). Nonetheless, the null-model was consistently worse than the other 3 models (Table 3), indicating that the global network metrics depend more on the patterning of social interactions produced by the mechanisms implemented in the models. On average, “Reaper” predicted the global network metrics’ values slightly better than the other models. With regard to modularity, however, it produced values that did not vary but remained constant; whereas in the spatially explicit models modularity varied according to the values observed in the empirical data (Tables 3 and 4, Figure 4). This suggests that a spatial structure is necessary to produce variation in modularity. Figure 3 View largeDownload slide Plots of the values of the observed global metric (x axis) against the values predicted by each model (y axis). Black diagonal represents values expected in a perfect match between model and empirical data. Colors represent different models and their corresponding trend lines. Data points represent the metric’s value of a single simulation (n = 20) from each model. CC = clustering coefficient. Figure 3 View largeDownload slide Plots of the values of the observed global metric (x axis) against the values predicted by each model (y axis). Black diagonal represents values expected in a perfect match between model and empirical data. Colors represent different models and their corresponding trend lines. Data points represent the metric’s value of a single simulation (n = 20) from each model. CC = clustering coefficient. Figure 4 View largeDownload slide Plots of the values of the observed global metric (x axis) against the values predicted by each model (y axis). Black diagonal represents values expected in a perfect match between model and empirical data. Colors represent different models and their corresponding trend lines. Data points represent the metric’s value of a single simulation (n = 20) from each model. Figure 4 View largeDownload slide Plots of the values of the observed global metric (x axis) against the values predicted by each model (y axis). Black diagonal represents values expected in a perfect match between model and empirical data. Colors represent different models and their corresponding trend lines. Data points represent the metric’s value of a single simulation (n = 20) from each model. With respect to grooming reciprocation, up the dominance hierarchy, and towards partners of similar rank, none of the models was able to fit all of the patterns across all empirical groups (Supplementary Tables SI9–SI11). “GrooFiWorld” failed to match 6 out of 42 patterns, “Reaper” 9, and “FriendsWorld” 10. It is noteworthy that grooming directed towards individuals of similar rank never reached statistical significance in “Reaper” (Supplementary Table SI11): the absence of this pattern can be explained by a lack of spatial structure. In “GrooFiWorld” and “FriendsWorld,” individuals of similar rank are often close due to the spatial structure; this induces frequent interactions among them, and the emergence of this pattern (Puga-Gonzalez et al. 2009). GrooFiWorld, FriendsWorld, and Reaper models reproduced reciprocation. We were however unable to firmly confirm which model, and thus which underlying mechanism, “proximity-based” or “emotional bookkeeping,” may reproduce more accurately the diversity of networks features and distribution of grooming observed in macaque groups. It could be argued that exploring the parameter space (all possible combinations of values for all models’ parameters, Supplementary Table SI2) may produce a better fit with empirical data. We think this is not the case. Models made poor predictions in particular for values of modularity, centrality index, and average eigenvector (Figure 4). Modularity and centrality index assess the degree of subgrouping and the popularity of the most central individual within the network, respectively; and eigenvector is a measurement of how well the individuals and their neighbors are connected within the network. To reproduce in the models a high degree of modularity and centrality index as observed in empirical data, it would be necessary that individuals, somehow, interact more with their subgroup than with the rest of the group, and that all or most individuals in the group have the same preferred interaction partner respectively. These patterns, thus, seem a consequence of the way individuals distribute grooming among their partners. The grouping, grooming, and most of the dominance parameters (Supplementary Table SI2), modify the rate of interactions, grooming frequency, and cohesiveness of the group but not the way individuals distribute their grooming (Hemelrijk and Puga-Gonzalez 2012; Hemelrijk et al. 2017). Hence, the lack of a better fit between the models and these global metrics indicates that an additional process(es), affecting the way individuals distribute their grooming, is missing in the computer models here analyzed. None of the mechanisms in these models (spatial structure and proximity-based interactions: GrooFiWorld; social bonding: Reaper; or their combination: FriendsWorld) is sufficient to produce attraction to specific individuals and subgrouping effects as significant as those observed in nature. In nonhuman primates, individuals appear to be knowledgeable not only of their own social relationships but also of those between third-parties (Silk 1999; Perry et al. 2004; Schino et al. 2006; Paxton et al. 2010; Kubenova et al. 2017). Thus, a potential mechanism missing in the models may be that individuals can choose partners, taking into account not only how other individuals behave towards them but also towards other group members. Implementing such a mechanism in the models could have an impact on triadic closure, and thus on the centralization and modularity of the network. Another potential mechanism is that, in nonhuman primates, the degree with which individuals reduce their anxiety after being involved in grooming depends on the identity of the partner. Individuals that groom with “friends” seem to relax more than when they groom with “nonfriends” (Crockford et al. 2013; Wittig et al. 2016). In the model, however, the reduction in anxiety after grooming was the same regardless of the identity of the partner. Further, individuals may take into account not only grooming exchange but also aggression exchange when establishing social bonds, whereas in “Friendsworld” and “Reaper” grooming exchange was the only variable taken into account. Furthermore, the model omitted other individual attributes such as age and kinship, which are known to influence social bonding and subgrouping. Further research will be necessary to know whether the implementation of these or other mechanisms in the models will produce a better fit between models’ and empirical grooming networks. The models fitted the data from wild (provisioned) groups better than that of captive ones. Indeed, in some cases with the captive data the predictions of the models were opposite of expected (SI11-SI12). The models, thus, seem to represent the behavior of wild animals better than those maintained in captive conditions (see R2 in Supplementary Table SI18). The lack of a better fit of the models with the latter may be due to the spatial constraint imposed by captivity on these groups. In the wild individuals can avoid others by moving away (Heesen et al. 2014; Heesen et al. 2015), whereas in captivity individuals have to interact with most group members. This may cause an artificial distribution of grooming which would ultimately affect grooming networks. In the spatially explicit models, although individuals followed rules that made them stay in the group (by searching for others and moving towards them when no other group member was in view), they were free to move in any direction since there were no spatial boundaries. Whether imposing a spatial boundary in these models would result in a better fit with empirical data from captive groups should be studied in the future. To fit the grooming patterns observed in some groups of macaques, in the models the parameter of aggression intensity had to be optimized (Supplementary Tables SI8). To simulate the group of Barbary macaques, for instance, aggression intensity was set to high values because females directed grooming up the dominance hierarchy and towards partners of similar rank in this group (Supplementary Tables SI6–SI7). Similarly, to simulate some groups of Japanese macaques, aggression intensity was set to low values because in none of these groups grooming was directed up the hierarchy (Supplementary Table SI6). According to the biological market theory, grooming directed up the dominance hierarchy is a pattern expected and usually observed in intolerant primate societies and not in tolerant ones (Barrett et al. 1999; Schino and Aureli 2008b). Some of the macaque groups here analyzed, thus, presented a pattern opposite to the one expected according to their species’ classification of social style. We do not believe, however, that this observation is sufficient to consider a reclassification of macaque species’ social style because 1) although from the groups here studied some presented a pattern opposite to the expected, other groups of the same species did conform to the patterns expected according to its social style, and 2) the classification of social styles is based not only on the distribution of social grooming but on other patterns like rates of reconciliation, directionality of aggression, or permissiveness of mothers (Thierry 2004). We have no indication that in our data set living conditions—wild or captive—have influenced the observed grooming patterns; actually, the patterns of 5 wild and 3 captive groups did not conform to the patterns previously documented for social style (but see Balasubramaniam et al. 2017 for a significant result on the influence of living condition on grooming traits). Using the principle of parsimony, the simplest model should be favored. In the models, however, the cognitive demands on individuals do not seem exceptionally high. In all of them, individuals should be capable of identifying other group-members and recognizing their dominance status. And in “FriendsWorld” and “Reaper” models individuals should additionally differentiate others on the basis of the grooming given and received, a task that could be achieved through associative learning. Even though cognitive demands on individuals seem lower in GrooFiWorld, nonhuman primates have been shown to possess these cognitive capacities anyhow (Tomasello and Call 1997). Thus, it is difficult to favor one model over the other on the basis of cognitive limitations. The models presented here are not the only ones that reproduce patterns of reciprocation of grooming similar to those observed in primate societies. In the partner-choice model by Schino and Campenni (Campenni and Schino 2014) and the EMO-model by Evers and co-authors (2015), preferential interactions alone or in combination with spatial structure, respectively, suffice to reproduce grooming reciprocation. Nonetheless, no analyses of grooming social networks, grooming directed up the dominance hierarchy or towards individuals of similar rank, and of interchange of grooming for support have been done on these models. Thus, it is unknown if in these models the distribution of grooming and/or the grooming social networks are similar to those reported in macaques. In conclusion, our results indicate that 3 models capture to some degree the diversity of social styles observed in the grooming patterns and networks of macaques. Our comparative analysis, however, does not allow us to make a clear distinction between the models. Whereas GrooFiWorld fits the distribution of grooming among group members best, “Reaper” and “FriendsWorld” do so for the network metrics. Thus, it is not possible to specify which mechanism, proximity-based or emotional book-keeping, is more accurate at reproducing reciprocation, grooming distribution and social networks as observed in macaques. What clearly appears, however, is that a process is missing in these models. This process could be a complex one, such as individuals taking into account third-party relationships when deciding whom to interact with. Whether such a mechanism can more accurately reproduce the diversity of social styles observed in macaque societies we will study in future models. SUPPLEMENTARY MATERIAL Supplementary data are available at Behavioral Ecology online. FUNDING This work was supported by a post-doctoral research grant from the National Council of Science and Technology (CONACyT) of Mexico awarded to I.P.G. Long-term research on Assamese macaques is supported by the German Primate Center and Georg-August University of Gottingen. Data accessibility: Analyses reported in this article can be reproduced using the data provided by Puga-Gonzalez et al. (2018). J.O. and O.S. thank the National Research Council of Thailand (NRCT) and the Department of National Parks, Wildlife and Plant Conservation (DNP) for permission to conduct this study, A. Koening and C. Borries who developed the field site at PKWS, and K. Nitaya T. Wongsnak, M. Pongjantarasatien, K. Kreetiyutanont, M. Kumsuk, and W. Saenphala for their cooperation. B.T. thanks Arianna de Marco for data collection. S.S. thanks Ellen Merz for allowing him to lead his research in La Forêt des Singes and for providing the demographic data. REFERENCES Aureli F , Cords M , Van Schaik CP . 2002 . Conflict resolution following aggression in gregarious animals: a predictive framework . Anim Behav . 64 : 325 – 343 . Google Scholar CrossRef Search ADS Aureli F , Preston SD , de Waal FB . 1999 . Heart rate responses to social interactions in free-moving rhesus macaques (Macaca mulatta): a pilot study . J Comp Psychol . 113 : 59 – 65 . Google Scholar CrossRef Search ADS PubMed Balasubramaniam KN , Beisner BA , Berman CM , De Marco A , Duboscq J , Koirala S , Majolo B , MacIntosh AJ , McFarland R , Molesti S , et al. 2017 . The influence of phylogeny, social style, and sociodemographic factors on macaque social network structure . Am J Primatol. 80 : e22727 . Barchas PR , Mendoza SD . 1984 . Emergent hierarchical relationships in rhesus macaques: An application of chase’s model . In: Barchas PR , editor. Social hierarchies: essays towards a sociophysiological perspective . Westport (CT) : Greenwood Press . p. 81 . Barrett L , Henzi SP , Weingrill T , Lycett JE , Hill RA . 1999 . Market forces predict grooming reciprocity in female baboons . Proc R Soc Lond B . 266 : 665 – 670 . Google Scholar CrossRef Search ADS Butovskaya ML . 2008 . Reconciliation, dominance and cortisol levels in children and adolescents (7–15-year-old boys) . Behaviour . 145 : 1557 – 1576 . Google Scholar CrossRef Search ADS Campennì M , Schino G . 2014 . Partner choice promotes cooperation: the two faces of testing with agent-based models . J Theor Biol . 344 : 49 – 55 . Google Scholar CrossRef Search ADS PubMed Clutton-Brock T . 2009 . Cooperation between non-kin in animal societies . Nature . 462 : 51 – 57 . Google Scholar CrossRef Search ADS PubMed Crockford C , Wittig RM , Langergraber K , Ziegler TE , Zuberbuehler K , Deschner T . 2013 . Urinary oxytocin and social bonding in related and unrelated wild chimpanzees . Proc Biol Sci . 280 : 20122765 . Google Scholar CrossRef Search ADS PubMed Crook JH , Ellis JE , Gosscustard JD . 1976 . Mammalian social-systems - structure and function . Anim Behav . 24 : 261 – 274 . Google Scholar CrossRef Search ADS Csardi G , Nepusz T . 2006 . The igraph software package for complex network research . InterJournal Complex Systems . 1695 . Dunbar RIM . 1991 . Functional significance of social grooming in primates . Folia Primatol . 57 : 121 – 131 . Google Scholar CrossRef Search ADS Evers E , de Vries H , Spruijt BM , Sterck EH . 2011 . Better safe than sorry–socio-spatial group structure emerges from individual variation in fleeing, avoidance or velocity in an agent-based model . PLoS One . 6 : e26189 . Google Scholar CrossRef Search ADS PubMed Evers E , de Vries H , Spruijt BM , Sterck EH . 2015 . Emotional bookkeeping and high partner selectivity are necessary for the emergence of partner-specific reciprocal affiliation in an agent-based model of primate groups . PLoS One . 10 : e0118921 . Google Scholar CrossRef Search ADS PubMed Franz M , McLean E , Tung J , Altmann J , Alberts SC . 2015 . Self-organizing dominance hierarchies in a wild primate population . Proc Biol Sci . 282 . Freeberg TM , Dunbar RI , Ord TJ . 2012 . Social complexity as a proximate and ultimate factor in communicative complexity . Philos Trans R Soc Lond B Biol Sci . 367 : 1785 – 1801 . Google Scholar CrossRef Search ADS PubMed Galef BG , Jr . 1988 . Imitation in animals: History, definitions, and interpretation of data from the psychological laboratory . In: Zentall T , Galef B , editors. Social learning: psychobiological and biological perspectives . Hillsdale (NJ) : Erlbaum . p. 3 . Graves FC , Wallen K , Maestripieri D . 2002 . Opioids and attachment in rhesus macaque (Macaca mulatta) abusive mothers . Behav Neurosci . 116 : 489 – 493 . Google Scholar CrossRef Search ADS PubMed Hamilton WD . 1964 . The genetical evolution of social behaviour. I . J Theor Biol . 7 : 1 – 16 . Google Scholar CrossRef Search ADS PubMed Heesen M , Macdonald S , Ostner J , Schülke O . 2015 . Ecological and social determinants of group cohesiveness and within-group spatial position in wild assamese macaques . Ethology . 121 : 270 – 283 . Google Scholar CrossRef Search ADS Heesen M , Rogahn S , Macdonald S , Ostner J , Schulke O . 2014 . Predictors of food-related aggression in wild assamese macaques and the role of conflict avoidance . Behav Ecol Sociobiol . 68 : 1829 – 1841 . Google Scholar CrossRef Search ADS Hemelrijk CK . 1990a . A matrix partial correlation test used in investigations of reciprocity and other social interaction patterns at group level . J Theor Biol . 143 : 405 – 420 . Google Scholar CrossRef Search ADS Hemelrijk CK . 1990b . Models of, and tests for, reciprocity, unidirectional and other social interaction patterns at a group level . Anim Behav . 39 : 1013 – 1029 . Google Scholar CrossRef Search ADS Hemelrijk CK . 1999 . An individual-oriented model on the emergence of despotic and egalitarian societies . P Roy Soc Lond B Bio . 266 : 361 – 369 . Google Scholar CrossRef Search ADS Hemelrijk CK . 2000 . Towards the integration of social dominance and spatial structure . Anim Behav . 59 : 1035 – 1048 . Google Scholar CrossRef Search ADS PubMed Hemelrijk CK , Kappeler PM , Puga-Gonzalez I . 2017 . The self-organization of social complexity in group-living animals: Lessons from the DomWorld model . In: Naguib M , Podos J , Simmons LW , et al. , editors. Advances in the study of behavior . p. 361 . Google Scholar CrossRef Search ADS Hemelrijk CK , Puga-Gonzalez I . 2012 . An individual-oriented model on the emergence of support in fights, its reciprocation and exchange . PLoS One . 7 : e37271 . Google Scholar CrossRef Search ADS PubMed Hemelrijk CK , Wantia J , Isler K . 2008a . Female dominance over males in primates: self-organisation and sexual dimorphism . PLoS One . 3 : e2678 . Google Scholar CrossRef Search ADS Hemelrijk CK , Wantia J , Isler K . 2008b . The more males, the more dominant are female primates . Folia Primatol . 79:337-337 . Hochberg Y . 1988 . A sharper bonferroni procedure for multiple tests of significance . Biometrika . 75:800-802 . Hogeweg P . 1988 . MIRROR beyond MIRROR, puddles of LIFE . In: Langton C , editor. Artificial life, SFI studies in the sciences of complexity . Redwood City (CA) : Adisson-Wesley Publishing Company . p. 297 Hsu Y , Wolf LL . 1999 . The winner and loser effect: integrating multiple experiences . Anim Behav . 57 : 903 – 910 . Google Scholar CrossRef Search ADS PubMed Jaeggi AV , De Groot E , Stevens JMG , Van Schaik CP . 2013 . Mechanisms of reciprocity in primates: testing for short-term contingency of grooming and food sharing in bonobos and chimpanzees . Evol Hum Behav . 34 : 69 – 77 . Google Scholar CrossRef Search ADS Kasper C , Voelkl B . 2009 . A social network analysis of primate groups . Primates . 50 : 343 – 356 . Google Scholar CrossRef Search ADS PubMed Keverne EB , Martensz ND , Tuite B . 1989 . Beta-endorphin concentrations in cerebrospinal fluid of monkeys are influenced by grooming relationships . Psychoneuroendocrinology . 14 : 155 – 161 . Google Scholar CrossRef Search ADS PubMed Krause J , Ruxton GD . 2002 . Living in groups . Oxford : Oxford University Press . Kubenova B , Konecna M , Majolo B , Smilauer P , Ostner J , Schülke O . 2017 . Triadic awareness predicts partner choice in male-infant-male interactions in Barbary macaques . Anim Cogn . 20 : 221 – 232 . Google Scholar CrossRef Search ADS PubMed Kudo H , Dunbar RIM . 2001 . Neocortex size and social network size in primates . Anim Behav . 62 : 711 – 722 . Google Scholar CrossRef Search ADS Lehmann J , Korstjens AH , Dunbar RIM . 2007 . Group size, grooming and social cohesion in primates . Anim Behav . 74 : 1617 – 1629 . Google Scholar CrossRef Search ADS Martel FL , Nevison CM , Simpson MJ , Keverne EB . 1995 . Effects of opioid receptor blockade on the social behavior of rhesus monkeys living in large family groups . Dev Psychobiol . 28 : 71 – 84 . Google Scholar CrossRef Search ADS PubMed McFarland R , Fuller A , Hetem RS , Mitchell D , Maloney SK , Henzi SP , Barrett L . 2015 . Social integration confers thermal benefits in a gregarious primate . J Anim Ecol . 84 : 871 – 878 . Google Scholar CrossRef Search ADS PubMed Molesti S , Majolo B . 2013 . Grooming increases self-directed behaviour in wild barbary macaques, macaca sylvanus . Anim Behav . 86 : 169 – 175 . Google Scholar CrossRef Search ADS Newman MEJ . 2006 . Modularity and community structure in networks . Proc Natl Acad Sci . 103 : 8577 – 8582 . Google Scholar CrossRef Search ADS PubMed Noë R , Hammerstein P . 1994 . Biological markets: supply and demand determine the effect of partner choice in cooperation, mutualism and mating . Behav Ecol Sociobiol . 35 : 1 – 11 . Google Scholar CrossRef Search ADS Pasquaretta C , Levé M , Claidière N , van de Waal E , Whiten A , MacIntosh AJ , Pelé M , Bergstrom ML , Borgeaud C , Brosnan SF , et al. 2014 . Social networks in primates: smart and tolerant species have more efficient networks . Sci Rep . 4 : 7600 . Google Scholar CrossRef Search ADS PubMed Paxton R , Basile BM , Adachi I , Suzuki WA , Wilson ME , Hampton RR . 2010 . Rhesus monkeys (Macaca mulatta) rapidly learn to select dominant individuals in videos of artificial social interactions between unfamiliar conspecifics . J Comp Psychol . 124 : 395 – 401 . Google Scholar CrossRef Search ADS PubMed Perry S , Barret CH , Manson JH . 2004 . White-faced capuchin monkeys show triadic awareness in their choices of allies . Animal Behav . 67 : 165 – 170 . Google Scholar CrossRef Search ADS Puga-Gonzalez I , Hildenbrandt H , Hemelrijk CK . 2009 . Emergent patterns of social affiliation in primates, a model . PLoS Comput Biol . 5 : e1000630 . Google Scholar CrossRef Search ADS PubMed Puga-Gonzalez I , Hoscheid A , Hemelrijk CK . 2015 . Friendships, reciprocation and interchange in an individual-based model . Behav Ecol Sociobiol . 69 : 383 – 394 . Google Scholar CrossRef Search ADS Puga-Gonzalez I , Ostner J , Schuelke O , Sosa S , Thierry B , Sueur C . 2018 . Data from: mechanisms of reciprocity and diversity in social networks: a modeling and comparative approach . Dryad Digital Repository . http://dx.doi.org/10.5061/dryad.6m220kd Puga-Gonzalez I , Sueur C . 2017a . Friendships and social networks in an individual-based model of primate social behaviour . J Artif Soc Soc Simul . 20 : 10 . Google Scholar CrossRef Search ADS Puga-Gonzalez I , Sueur C . 2017b . Emergence of complex social networks from spatial structure and rules of thumb: a modelling approach . Ecol Complex . 31 : 189 – 200 . Google Scholar CrossRef Search ADS Rankin RW , Mann J , Singh L , Patterson EM , Krzyszczyk E , Bejder L . 2016 . The role of weighted and topological network information to understand animal social networks: a null model approach . Anim Behav . 113 : 215 – 228 . Google Scholar CrossRef Search ADS Sabbatini G , De Bortoli Vizioli A , Visalberghi E , Schino G . 2012 . Food transfers in capuchin monkeys: an experiment on partner choice . Biol Lett . 8 : 757 – 759 . Google Scholar CrossRef Search ADS PubMed Schino G . 2007 . Grooming and agonistic support: a meta-analysis of primate reciprocal altruism . Behav Ecol . 18 : 115 – 120 . Google Scholar CrossRef Search ADS Schino G , Aureli F . 2008a . Grooming reciprocation among female primates: a meta-analysis . Biol Lett . 4 : 9 – 11 . Google Scholar CrossRef Search ADS Schino G , Aureli F . 2008b . Trade-offs in primate grooming reciprocation: testing behavioural flexibility and correlated evolution . Biol J Linn Soc . 95 : 439 – 446 . Google Scholar CrossRef Search ADS Schino G , Aureli F . 2009 . Reciprocal altruism in primates: partner choice, cognition, and emotions . Adv Stud Behav . 39 : 45 – 69 . Google Scholar CrossRef Search ADS Schino G , Scucchi S , Maestripieri D , Turillazzi PG . 1988 . Allogrooming as a tension-reduction mechanism: a behavioral approach . Am J Primatol . 16 : 43 – 50 . Google Scholar CrossRef Search ADS Schino G , Tiddi B , Di Sorrentino EP . 2006 . Simultaneous classification by rank and kinship in japanese macaques . Anim Behav . 71 : 1069 – 1074 . Google Scholar CrossRef Search ADS Schülke O , Bhagavatula J , Vigilant L , Ostner J . 2010 . Social bonds enhance reproductive success in male macaques . Curr Biol . 20 : 2207 – 2210 . Google Scholar CrossRef Search ADS PubMed Semple S , Harrison C , Lehmann J . 2013 . Grooming and anxiety in barbary macaques . Ethology . 119 : 779 – 785 . Google Scholar CrossRef Search ADS Setchell JM , Smith T , Wickings EJ , Knapp LA . 2008 . Social correlates of testosterone and ornamentation in male mandrills . Horm Behav . 54 : 365 – 372 . Google Scholar CrossRef Search ADS PubMed Shutt K , MacLarnon A , Heistermann M , Semple S . 2007 . Grooming in Barbary macaques: better to give than to receive ? Biol Lett . 3 : 231 – 233 . Google Scholar CrossRef Search ADS PubMed Silk JB . 1999 . Male bonnet macaques use information about third-party rank relationships to recruit allies . Anim Behav . 58 : 45 – 51 . Google Scholar CrossRef Search ADS PubMed Silk JB . 2007 . Social components of fitness in primate groups . Science . 317 : 1347 – 1351 . Google Scholar CrossRef Search ADS PubMed Silk JB , Alberts SC , Altmann J . 2003 . Social bonds of female baboons enhance infant survival . Science . 302 : 1231 – 1234 . Google Scholar CrossRef Search ADS PubMed Stevens JR , Hauser MD . 2004 . Why be nice? Psychological constraints on the evolution of cooperation . Trends Cogn Sci . 8 : 60 – 65 . Google Scholar CrossRef Search ADS PubMed Stevens JR , Volstorf J , Schooler LJ , Rieskamp J . 2011 . Forgetting constrains the emergence of cooperative decision strategies . Front Psychol . 1 : 235 . Google Scholar CrossRef Search ADS PubMed Sueur C , Deneubourg JL , Petit O , Couzin ID . 2011 . Group size, grooming and fission in primates: a modeling approach based on group structure . J Theor Biol . 273 : 156 – 166 . Google Scholar CrossRef Search ADS PubMed Sueur C , Petit O , De Marco A , Jacobs AT , Watanabe K , Thierry B . 2011 . A comparative network analysis of social style in macaques . Anim Behav . 82 : 845 – 852 . Google Scholar CrossRef Search ADS Team-R-Core . 2015 . R: a language and environment for statistical computing . Vienna (Austria) : R Foundation for Statistical Computing . Thierry B . 2004 . Social epigenesis . In: Thierry B , Singh M , Kaumanns W , editors. Macaque societies: A model for the study of social organisation . Cambridge : Cambridge University Press . p. 267 . Tiddi B , Aureli F , di Sorrentino EP , Janson CH , Schino G . 2011 . Grooming for tolerance? two mechanisms of exchange in wild tufted capuchin monkeys . Behav Ecol . 22 : 663 – 669 . Google Scholar CrossRef Search ADS Tomasello M , Call J . 1997 . Primate cognition . New York : Oxford University Press . Trivers RL . 1971 . The evolution of reciprocal altruism . Q Rev Biol . 46 : 35 – 57 . Google Scholar CrossRef Search ADS de Waal FBM . 2000 . Attitudinal reciprocity in food sharing among brown capuchin monkeys . Animal Behav . 60 : 253 – 261 . Google Scholar CrossRef Search ADS de Waal FBM , Brosnan SF . 2006 . Simple and complex reciprocity in primates . In: Kappeler PM , van Schaick CP , editors. Cooperation in primates and humans: Mechanisms and evolution . Berlin (Germany) : Springer . p. 85 . Google Scholar CrossRef Search ADS de Waal FBM , Luttrell LM . 1988 . Mechanisms of social reciprocity in three primate species: symmetrical relationship characteristics or cognition ? Ethol Sociobiol . 9 : 101 – 118 . Google Scholar CrossRef Search ADS Ward A , Webster M . 2016 . Sociality: the behaviour of group-living animals . 1st ed . Switzerland: Springer International Publishing . Google Scholar CrossRef Search ADS Wittig RM , Crockford C , Weltring A , Langergraber KE , Deschner T , Zuberbühler K . 2016 . Social support reduces stress hormone levels in wild chimpanzees across stressful events and everyday affiliations . Nat Commun . 7 : 13361 . Google Scholar CrossRef Search ADS PubMed © The Author(s) 2018. Published by Oxford University Press on behalf of the International Society for Behavioral Ecology. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Behavioral Ecology Oxford University Press

Mechanisms of reciprocity and diversity in social networks: a modeling and comparative approach

Loading next page...
 
/lp/ou_press/mechanisms-of-reciprocity-and-diversity-in-social-networks-a-modeling-yYeLroSmwG
Publisher
Oxford University Press
Copyright
© The Author(s) 2018. Published by Oxford University Press on behalf of the International Society for Behavioral Ecology. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com
ISSN
1045-2249
eISSN
1465-7279
D.O.I.
10.1093/beheco/ary034
Publisher site
See Article on Publisher Site

Abstract

Abstract Three mechanisms have been proposed to underlie reciprocation of social behaviors in gregarious animals: “calculated reciprocity,” “emotional bookkeeping,” and “symmetry-based reciprocity.” Among these explanations, emotional book-keeping has received the broadest support from experimental and observational studies. On the other hand, 3 individual-based models have shown that reciprocation may emerge via “symmetry-based reciprocity,” “emotional bookkeeping,” or a combination of both mechanisms. Here, we use these 3 models to assess their relative fit with empirical data on reciprocation and social network structure across different groups and species of macaques. We collected grooming data from 14 groups and 8 macaque species and simulated each group in each model. We analyzed and quantitatively compared social network metrics of the empirical and the models’ grooming networks. The 3 models captured fairly well the features of observed networks, and fitted data from wild groups better than captive ones. The emotional bookkeeping model seemed on average to fit slightly better the social networks metrics observed in empirical data, but failed to reproduce some grooming patterns. The symmetry-based models, on the other hand, fitted better other network parameters (e.g., modularity). No model generally fitted the data better than the others, and the fit with some metrics (e.g., modularity, centralization index) was low even after optimization. Thus, our analyses indicate that in the models social interactions may be simpler than in reality and models may miss social processes (e.g., third-party awareness). INTRODUCTION Group living is a common phenomenon across the animal kingdom, from invertebrates (e.g., ants, bees, spiders) to marine (e.g., dolphins, whales, fish) and terrestrial vertebrates (e.g., rodents, primates, felines, birds) (Krause and Ruxton 2002; Ward and Webster 2016). Group living may be widespread due to the many adaptive advantages it has for individuals, including the reduction of predation risks, better access to resources, and community breeding and rearing of offspring (Crook et al. 1976). In species of animals living in permanent, stable groups, individuals engage frequently in social interactions with other group members under many different circumstances and during periods of time that can span their whole lifetime, giving then rise to complex social systems (Freeberg et al. 2012). In such systems, individual recognition and repeated interactions between certain group members may promote the emergence of reciprocation, i.e., the exchange of beneficial behaviors over long-time periods. Patterns of reciprocation have been extensively studied in primate social systems—e.g., grooming and support in fights—and their ultimate functions and proximate mechanisms largely debated. At an ultimate level, reciprocation may be explained via kin selection (Hamilton 1964), reciprocal altruism (Trivers 1971), mutualism (Clutton-Brock 2009), and/or the biological market theory (Noë and Hammerstein 1994). At a proximate level, 3 mechanisms, differing in the degree of cognition required, have been proposed. At the highest cognitive level is “calculated reciprocity” (de Waal and Luttrell 1988). This mechanism implies keeping record of all services given and received from others, long-term memory of past events, and paying back accordingly (de Waal and Brosnan 2006). A cognitively simpler mechanism is “attitudinal reciprocity” (de Waal 2000), also named “emotional bookkeeping” (Schino and Aureli 2009). Emotional bookkeeping means that individuals develop, through social interactions, a specific emotionally mediated attitude toward each group member, and whether or not they exchange a service with a group member depends on the type of affiliative relationship associated with it (Schino and Aureli 2009). At the lowest cognitive level is “symmetry-based reciprocity.” This mechanism assumes that the balance between services given and received results from repeated interactions between individuals sharing a common or symmetrical variable—e.g., spatial proximity, kinship, rank, age, etc. (de Waal and Luttrell 1988). Currently, researchers find calculated reciprocity an implausible hypothesis because of the difficulty that represents keeping record of behaviors over long periods of time (Stevens and Hauser 2004; Stevens et al. 2011). Symmetry-based reciprocity has been disregarded because even when symmetrical variables are partialed out, the correlation between behaviors given and received remains significant, suggesting thus that the correlation is not a by-product of the symmetrical variable. The statistical test used to partial out a variable, however, may be insufficient to completely exclude the effects of the variable (Hemelrijk and Puga-Gonzalez 2012); symmetry-based reciprocity thus, has potentially been underestimated. Emotional bookkeeping is nowadays the leading hypothesis regarding the mechanisms giving rise to reciprocity of social behaviors. Indeed, several empirical studies have shown that individuals prefer to associate and reciprocate behaviors with some partners rather than with others (Tiddi et al. 2011; Sabbatini et al. 2012; Jaeggi et al. 2013) and that “social bonds” are important for the individuals’ survival and that of their offspring (Silk et al. 2003; Silk 2007; Schüelke et al. 2010; McFarland et al. 2015). An approach to the study of the proximate mechanisms underlying social behavior in nonhuman primates makes use of individual-based computer models and complex science—i.e., the study of phenomena emerging from the interactions among the elements of a system. Several different individual-based models have shown that simple behavioral rules and local interactions may suffice to explain many social patterns observed in primate societies (Evers et al. 2011; Campenni and Schino 2014; Evers et al. 2015; Hemelrijk et al. 2017). In this regard, 3 individual-based computer models—GrooFiWorld, FriendsWorld, and Reaper—have been able to reproduce not only social behaviors but also the behavioral contrasts observed between tolerant and intolerant societies of monkeys, especially macaques (Puga-Gonzalez et al. 2009; Hemelrijk and Puga-Gonzalez 2012; Puga-Gonzalez et al. 2015; Puga-Gonzalez and Sueur 2017a). The Macaca genus comprises 23 species in which females are the philopatric sex, which results in female-bonded societies. On the basis of the females’ social behavior, macaque species have been classified along a social style gradient going from extremely intolerant to extremely tolerant. Among females from species in grade 1 and 2 (extremely intolerant and intolerant), the dominance hierarchy is steep, conflicts are unidirectional, aggression is fierce, reconciliation is infrequent, and social grooming is directed up the dominance hierarchy and towards individuals of similar rank. On the contrary, among females from species in grades 3 and 4 (tolerant and extremely tolerant), the dominance hierarchy is more shallow, conflicts are more bidirectional, and aggression is mild (Thierry 2004). Further, in all 4 grades individuals cooperate by reciprocating (exchange of same behaviors) and interchanging (exchange of different behaviors) grooming and support in fights (Schino 2007; Schino and Aureli 2008a). All these behavioral patterns (reciprocation and interchange of grooming and support), as well as the behavioral differences between intolerant and tolerant societies (grooming up the dominance hierarchy and towards individuals of similar rank, uni- (bi-) directional aggression, high [low] steepness of the hierarchy), have been reproduced by the 3 individual-based models. The individual-based model GrooFiWorld proposes that spatial proximity (i.e., proximity-based reciprocity), is the main cause of the emergence of complex behavior in societies of macaques (Puga-Gonzalez et al. 2009). In GrooFiWorld, a spatial structure with dominant individuals in the centre and subordinates at the periphery emerges as a side-effect of aggression. Due to this spatial structure, individuals interact more with some than with others, and reciprocity of grooming and support in fights; and the interchange of grooming and support emerge. Interestingly, these patterns disappear in the model only when proximity is omitted by making individuals interact at random, but not when proximity is statistically partialed out with the partial Tau-kr test (Hemelrijk 1990a), which indicates that partial correlations do not sufficiently exclude the effect of proximity (Hemelrijk and Puga-Gonzalez 2012). Recently, GrooFiWorld was extended with a mechanism similar to emotional-bookkeeping (Puga-Gonzalez et al. 2015). In this model, called “FriendsWorld,” individuals have a tendency to follow their “friends”, i.e., those individuals with whom they groom the most. The model shows that this mechanism of “follow your friends,” reinforces patterns of reciprocity and interchange of grooming and support in fights. It remained unclear however, what factor was more important for the emergence of reciprocation and interchange, the spatial structure of the group, or the “follow your friends” rule. More recently, a new individual-based model has been developed. This model is identical to GrooFiWorld regarding the way individuals choose to interact with their partners. Contrary to GrooFiWorld, this model is not spatially-explicit; thus, interactions are not based on proximity. Instead, individuals interact first with others at random in this model, and as time goes by they develop a preference for interacting with those with whom they exchange more grooming, resembling thus emotional bookkeeping; we called this model Reaper (i.e., reinforcement of affiliative preferences) (Puga-Gonzalez and Sueur 2017a). These 3 models have been successful at reproducing empirical patterns of reciprocation and, at the same time, suggest different factors as potential causes. In these studies, however, no social network analysis (SNA) was done. Given that SNA provides us with a detailed picture of how individuals interact and distribute grooming behavior at the individual and group level; here, we used a SNA to distinguish which of these models and underlying mechanisms (proximity-based reciprocity, emotional bookkeeping, or a combination of both), if any, more accurately represents empirical data. We studied individual and global network metrics to assess whether in the model the individuals’ connections in the social network as well as the global features of the network were similar to empirical data. We focused on 4 different individual metrics: betweenness, closeness, eigenvector, and clustering coefficient; and 7 different global metrics: diameter, number of triangles, global clustering coefficient, modularity, average eigenvector, centrality index, and the correlation coefficient between individuals’ dominance rank and eigenvector (see Table 1 for an interpretation of the meaning of each network metric). We chose these metrics because they have been shown to be relevant in studies of primate social networks. For instance, Kasper and Voelkl (2009) suggest the use of betweenness and closeness (for sparse networks), and eigenvector (for near-to-all-connected networks) to understand the involvement, popularity, and sociability of an individual within the social network. Modularity is used to measure the degree of fragmentation of the network into subgroups, and clustering coefficient to characterize the cliquishness of the network by measuring the extent to which individuals linked to one individual are also linked to each other (Kasper and Voelkl 2009; Sueur, Petit et al. 2011; Pasquaretta et al. 2014). The correlation between eigenvector and dominance rank is used to measure the degree of popularity of dominant individuals in the group (Sueur, Petit, et al. 2011). To measure the centralization of the network around the most central individual, it is recommended to calculate an index obtained by dividing the centralization of the observed network by the maximum centralization possible—that of a star network (Kasper and Voelkl 2009; Sueur, Petit, et al. 2011; Pasquaretta et al. 2014). The approach taken in the present study was exploratory without clear predictions about which network measure should be best predicted by which model. Table 1 Interpretation of the network metrics used in the network analysis Metric Interpretation Individual: Betweenness The degree to which an individual lies in-between the paths of others, and thus the influence that exerts on the passing of information between others. Closeness How close an individual is to all others in the network, and thus how fast information originating from this individual may reach all others Eigenvector “Popularity” (centrality) of an individual; takes into account how well connected are the focal individual and its neighbors (individuals linked to the focal) Clustering coefficient The relative number of cliques (of 3 individuals) an individual forms with its neighbors. Global: Diameter How fast information travels between the 2 least well connected individuals in the network. Number of triangles The absolute number of cliques (of 3 individuals) Clustering coefficient The relative number of cliques (of 3 individuals) Eigenvector The average degree of “popularity” among individuals in the network Modularity The degree of subgrouping or clustering of the network Centrality index The extent to which a network is dominated by the most or most “popular” individuals in the network Correlation coefficient eigenvector-dominance rank The extent to which dominant individuals are the most “popular” individuals within the network Metric Interpretation Individual: Betweenness The degree to which an individual lies in-between the paths of others, and thus the influence that exerts on the passing of information between others. Closeness How close an individual is to all others in the network, and thus how fast information originating from this individual may reach all others Eigenvector “Popularity” (centrality) of an individual; takes into account how well connected are the focal individual and its neighbors (individuals linked to the focal) Clustering coefficient The relative number of cliques (of 3 individuals) an individual forms with its neighbors. Global: Diameter How fast information travels between the 2 least well connected individuals in the network. Number of triangles The absolute number of cliques (of 3 individuals) Clustering coefficient The relative number of cliques (of 3 individuals) Eigenvector The average degree of “popularity” among individuals in the network Modularity The degree of subgrouping or clustering of the network Centrality index The extent to which a network is dominated by the most or most “popular” individuals in the network Correlation coefficient eigenvector-dominance rank The extent to which dominant individuals are the most “popular” individuals within the network See Methods for a technical description of each network metric. View Large Table 1 Interpretation of the network metrics used in the network analysis Metric Interpretation Individual: Betweenness The degree to which an individual lies in-between the paths of others, and thus the influence that exerts on the passing of information between others. Closeness How close an individual is to all others in the network, and thus how fast information originating from this individual may reach all others Eigenvector “Popularity” (centrality) of an individual; takes into account how well connected are the focal individual and its neighbors (individuals linked to the focal) Clustering coefficient The relative number of cliques (of 3 individuals) an individual forms with its neighbors. Global: Diameter How fast information travels between the 2 least well connected individuals in the network. Number of triangles The absolute number of cliques (of 3 individuals) Clustering coefficient The relative number of cliques (of 3 individuals) Eigenvector The average degree of “popularity” among individuals in the network Modularity The degree of subgrouping or clustering of the network Centrality index The extent to which a network is dominated by the most or most “popular” individuals in the network Correlation coefficient eigenvector-dominance rank The extent to which dominant individuals are the most “popular” individuals within the network Metric Interpretation Individual: Betweenness The degree to which an individual lies in-between the paths of others, and thus the influence that exerts on the passing of information between others. Closeness How close an individual is to all others in the network, and thus how fast information originating from this individual may reach all others Eigenvector “Popularity” (centrality) of an individual; takes into account how well connected are the focal individual and its neighbors (individuals linked to the focal) Clustering coefficient The relative number of cliques (of 3 individuals) an individual forms with its neighbors. Global: Diameter How fast information travels between the 2 least well connected individuals in the network. Number of triangles The absolute number of cliques (of 3 individuals) Clustering coefficient The relative number of cliques (of 3 individuals) Eigenvector The average degree of “popularity” among individuals in the network Modularity The degree of subgrouping or clustering of the network Centrality index The extent to which a network is dominated by the most or most “popular” individuals in the network Correlation coefficient eigenvector-dominance rank The extent to which dominant individuals are the most “popular” individuals within the network See Methods for a technical description of each network metric. View Large We collected data from 14 groups belonging to 8 different macaque species, comprising the 4 grades of social style, and did a SNA on each group. Due to limitations in empirical data and because social styles are based on females’ social behavior, our analyses were restricted to grooming social networks of adult females. Then, we simulated the 14 macaque groups in each of the 3 models, analyzed the resulting social networks, and quantitatively compared the results of the models and those of the empirical groups. As a control, we built an extra individual-based model (the “null-model”) in which individuals randomly selected interaction partners, and also randomly chose whether to groom or fight. Before running simulations, we mimicked the following variables from empirical groups in each of the models: group size, sex ratio, intensity of aggression, and density of grooming network. After running the simulations in the models we realized that there was a poor fit between the models and empirical data regarding some global network metrics (e.g., modularity, centralization index). These network metrics may be influenced by the way individuals distribute grooming among group members: whether or not grooming is reciprocated, directed up the dominance hierarchy, and/or toward partners of similar rank. Hence, to optimize the fit between models and empirical data we decided to mimic in the models the distribution of grooming observed in the empirical groups (see Methods for details). Previous analyses of the grooming networks from GrooFiWorld and Reaper models showed that these networks resembled qualitatively those of nonhuman primates (Puga-Gonzalez and Sueur 2017a; Puga-Gonzalez and Sueur 2017b). Hence, we had no specific predictions as to which model may fit empirical data better, except for the metric of modularity. From modularity, we know that the Reaper model fails at reproducing the variation observed in primate groups (Puga-Gonzalez and Sueur 2017a). Therefore, we expected the Reaper model to perform poorly regarding this metric. Overall, our main aim was to investigate which model reproduced quantitatively better the grooming network structure of the macaque groups analyzed. In this way, we aimed to know which model and thus which underlying mechanism is more likely to explain the diversity of social styles observed in societies of macaques. METHODS Macaque species From the primate literature, we surveyed studies on macaque species which reported a matrix of directed grooming, dominance hierarchy, and sampling method (scan or focal). In case of focal sampling, social grooming matrices were corrected for observation time, i.e., individual’s rows were divided by the number of times the individual was observed. In these cases thus, grooming matrices contained relative grooming frequency. Seven studies comprising 10 groups from 5 species met these criteria (Supplementary Table SI1). Further, coauthors of this study provided data from 4 groups from 4 species of macaques. Thus, our data collection comprised a total of 14 groups from 8 macaque species representing all of the 4 categorical grades of social styles attributed to macaque societies (Thierry 2004). Since in some of these studies data were available only for adult female individuals, we restricted our analyses to this category of individuals in all the groups. Note that due to the small number of adult females (n = 4) in the group of M. tonkeana (Orangerie Zoo), we included 2 subadult females (>4 years) in the analysis of this group. A full description of the groups included in this study is presented in Supplementary Table SI1. Data available from the Dryad Digital Repository (Puga-Gonzalez et al. 2018). Distribution of social grooming Whether social grooming was reciprocated, directed up the dominance hierarchy, or toward individuals of similar dominance rank was evaluated by means of the Tau-Kr matrix correlation (Hemelrijk 1990b). The level of significance was calculated using 2000 permutations and P-values < 0.05. To check for reciprocation of grooming, we correlated the groom-given matrix with the groom-received matrix. A significant positive correlation implies reciprocity. To check for grooming directed up the dominance hierarchy and to females of similar dominance rank, we tested the correlation between the grooming given matrix with the partner-rank matrix and the similar-rank matrix respectively. The partner-rank-matrix had the dominance rank of individuals in the rows, and the similar-rank-matrix was filled with zeros apart from the partners closest and second closest in dominance rank which were indicated as 1’s. A significantly positive correlation with the partner-rank-matrix corresponds to grooming being directed up the hierarchy, and a significantly positive correlation with the similar-rank-matrix corresponds to grooming directed to females of similar rank. Social networks metrics SNA were performed using R statistical software, version 3.2.2 (Team-R-Core 2015) and package igraph (Csardi and Nepusz 2006). Matrices of grooming given were transformed to directed graphs with the igraph package; then, the analysis of the directed grooming network was performed. Analyses were performed at the individual and at the group level. Individual level metrics We calculated 4 different metrics: betweenness, eigenvector, clustering coefficient, and closeness. These are metrics commonly used in analyses of animal societies and they inform us about the patterning of the individuals’ social interactions (Table 1). The technical definition of each metric is described below. Note that in some cases we used binary instead of weighted matrices to calculate some of the network metrics. We used binary matrices when the metric calculated was affected by the weight in the matrix (absolute grooming frequency) or when the metric can only be calculated from binary matrices. To convert matrices with absolute or relative grooming frequency to binary ones, we set all cell values > 0 to 1 and all others were set to 0. Metrics calculated from binary networks were betweenness centrality, closeness centrality, and clustering coefficient. To calculate eigenvector, we used weighted networks (with absolute or relative grooming frequency) because this metric can cope with weights. Unweighted betweenness centrality The number of times a node (i.e., an individual) acts as a bridge along the shortest path between 2 other nodes. Thus, it is a measure of how well an individual connects different subgroups or clusters within the network. Unweighted closeness centrality A measure of the degree to which a node is near all other nodes in a network. It is the average length of the shortest path between the node and all other nodes in the graph. In other words, this metric calculates how fast a node can reach all others. Clustering coefficient It is the number of closed-triangles a node form with its neighbors divided by the total possible number of closed-triangles that could exist among them. Thus, it assesses the degree to which nodes tend to cluster together. This measure takes into account neither direction nor weight. Weighted eigenvector It takes into account not only the number and strength of connections of the individual but also those of the partners to which it is connected. Since we used directed networks, here, high eigenvector means that an individual receives and has strong connections to other individuals who themselves receive grooming frequently. Group level metrics We chose 7 metrics that have been shown to be relevant in studies of primate social networks (Kasper and Voelkl 2009; Sueur, Petit, et al. 2011; Pasquaretta et al. 2014). All of these metrics are used to characterize the overall network structure, and each of them captures a different feature of the network. We used these metrics in order to get a detail comparison of the features the models capture better from empirical data and thus get a better assessment of which model (and mechanism) is most accurate at fitting the overall network structure. Number of triangles The absolute number of triangles in the network. Neither direction nor weights are taken into account. Unweighted diameter It is the longest path among all the shortest path lengths between 2 nodes, i.e., the shortest distance between the 2 most distant nodes in the network. Modularity The difference between the proportion of the total association of individuals within clusters (i.e., subgroups) and the expected proportion, given the summed associations of the different individuals. This metric is calculated using the eigenvector of each individual (Newman 2006). A high value of modularity means a high number of contacts within a subgroup, but few contacts between subgroups and low modularity means a homogeneous distribution of contacts between all group members. Mean clustering coefficient The proportion between the number of triangles in the network divided by total possible number of triangles. Neither direction nor weights are taken into account. Mean eigenvector The average of individual’s eigenvector coefficients per group. Centralization of the network It was calculated according to Equation 1, where the numerator, Cmax, is the highest eigenvector in the group and the denominator, Max, is the value obtained assuming the highest centralization possible, i.e., if the network were a star (Pasquaretta et al. 2014). Hence, the higher the value the more centralized is the network around one or several individuals. CIB=100*∑in(cmax−ci)Max∑in(cmax−ci) (1) Correlation dominance-eigenvector Spearman rank correlation coefficient between the individual’s eigenvector and the individual’s hierarchical rank. A positive correlation indicates a centralization of dominant individuals in the network. Agent-based computational models Here, we only present a brief description of each model. For a full description of the processes of the model, we refer to previous publications. GrooFiWorld model GrooFiWorld (Grooming and Fighting) is an individual-based model that is spatially explicit (Puga-Gonzalez et al. 2009). The model comprises a continuous 2-dimensional “world” (without borders) in which individuals are able to move in all directions. Individuals have a fixed vision angle (VisionAngle, Supplementary Table SI2) and a maximum perception distance (MaxView, Supplementary Table SI2). At the beginning of each simulation, the individuals’ location is assigned randomly within a previously defined radius (InitRadius, Supplementary Table SI2) calculated by multiplying group size by an arbitrary constant. To regulate individuals’ activations, each individual is attributed a random waiting time drawn from a uniform distribution; the individual with the shortest waiting time gets activated first. These waiting times are combined with a kind of social facilitation (Galef 1988) implying that an individual’s waiting time is reduced when a dominance interaction occurs close by (radius of social facilitation, Supplementary Table SI2). Intensity of aggression is reflected by the StepDom value (Supplementary Table SI2). Fierce aggression (e.g., bites), as in intolerant macaque societies, is represented by high values, and mild aggression (e.g., threats, slaps), as in tolerant societies, is represented by low values. To represent sexual dimorphism, females have a StepDom value 80% lower than that of males (Hemelrijk et al. 2008a). Grouping rules Why individuals group (to avoid predators, search resources, etc.) is not specified in the model. Whenever an individual does not see another close by (i.e., within its personal space, PersSpace, Supplementary Table SI2); then, it looks for others at increasing distances (NearView and MaxView, Supplementary Table SI2). If, even then none is in the field of vision, then the individual turns over a SearchAngle (Supplementary Table SI2) in order to look for others. In this way individuals tend to remain in a group. If, however, an individual “sees” another within its personal space (PersSpace, Supplementary Table SI2), then a social interaction may take place. Interactions Upon meeting another in its personal space, an individual first chooses whether or not to perform a dominance interaction. This choice depends on the risks involved, whereby risk concerns the chance of losing a fight (i.e., “Risk-aversion strategy”, Hemelrijk 2000). If an individual chooses not to fight; then, it considers grooming the other. If the individual chooses not to groom, no interaction happens and the individual stays put. Dominance interactions Dominance interactions are modeled as in Hemelrijk (1999) and are an extension of the DoDom rules of Hogeweg (1988). Each individual has a dominance value, Dom (Supplementary Table SI2), which represents its capacity to win. A dominance interaction takes place only if an individual expects to be victorious, i.e. individuals avoid risks. These risks are estimated by means of a “mental battle.” During a “mental battle”, an individual i compares its dominance value (Domi) relatively to that of its opponent j (Domj) to a random value between zero and one (Equation 2). This process is repeated RiskSens times (Supplementary Table SI2). In the current simulation RiskSens is set to 2, thus individuals have to win a mental battle twice before engaging in a real dominance interaction. wi=[1 DOMiDOMi+DOMj>RND(0,1)0 else (2) If in both “mental battles” individual i expects to be victorious, then a real dominance interaction occurs. The outcome of the real dominance interaction is again decided according to Equation 2. To reflect the self-reinforcing effects of victory and defeat (Barchas and Mendoza 1984; Hsu and Wolf 1999; Setchell et al. 2008; Hemelrijk et al. 2008a; Franz et al. 2015), dominance values are updated by increasing the dominance value of the winner and decreasing that of the loser by the same amount (Equation 3). This positive feedback is “dampened” because a victory of a higher-ranking opponent increases its relative Dom-value only slightly, whereas an (unexpected) success of the lower-ranking individual increases its relative dominance value by a greater change. To keep Dom-values positive, their minimum value is, arbitrarily, set at 0.01. DOMi=DOMi+(wi−DOMiDOMi+DOMj)*STEPDOMDOMj=DOMj−(wi−DOMjDOMi+DOMj)*STEPDOM (3) The change in Dom-values is multiplied by a scaling factor, StepDom, which represents intensity of aggression (Hemelrijk 1999). StepDom values range from 0.08 to 0.8 in females and from 0.1 to 1 in males. High values of StepDom cause a bigger change in Dom-values than low values of StepDom. After a fight, the winner chases the loser over a distance of one unit and the loser responds by fleeing under a small random angle over a predefined FleeingDistance (Supplementary Table SI2). To reduce consecutive interactions between same opponents, after the chase the winner turns randomly 45° to right or left. Grooming interactions When an individual meets another in its PersSpace and “chooses” not to fight, then, it considers grooming its partner. Grooming is induced by the level of Anxiety (Supplementary Table SI2), which ranges from very relaxed to very tense on a scale from 0 to 1. Individuals groom their partners if their level of Anxiety is higher than a random number between 0 and 1; otherwise, they display “nonaggressive” proximity. As indicated by empirical studies, grooming (giving and receiving) reduces anxiety and thus the tendency to groom (Schino et al. 1988; Aureli et al. 1999; Shutt et al. 2007; but see Molesti and Majolo 2013; Semple et al. 2013). It does so more strongly in the groomee (AnxDcrGree) than in the groomer (AnxDcrGrmr) (Supplementary Table SI2). Further, during periods without grooming, individuals increase their Anxiety with AnxInc (Supplementary Table SI2) and thus, their motivation to groom as demonstrated in empirical studies (Keverne et al. 1989; Martel et al. 1995; Graves et al. 2002). Furthermore, Anxiety increases in both opponents after a fight in the model (AnxIncFight, Supplementary Table SI2), as reported in nonhuman primates (Aureli et al. 2002; Butovskaya 2008). To avoid consecutive interactions between same partners, after grooming both partners turn randomly 45° to the right or left. The default grooming parameters were tuned in order to match the percentage of time grooming observed in primate societies ~20% (Dunbar 1991). Among macaque species, there is some variation regarding the percentage of time allocated to grooming; however given that the most accurate percentage of grooming time in each macaque species is unknown, we decided to tune all models to ~20%. Nevertheless, grooming patterns (reciprocation and interchange of grooming and support; grooming up the dominance hierarchy and among individuals of similar rank) emerge in the models as long as the percentage of time allocated to grooming is ≥5% (Puga-Gonzalez et al., unpublished data). FriendsWorld In contrast to GrooFiWorld, in FriendsWorld, individuals categorize others as “friends” or not. We classified as “friends” the top 5 partners with whom an individual groomed the most (given and received). We choose 5 individuals because this is the average number of grooming partners found in empirical studies irrespective of group size (Kudo and Dunbar 2001; Lehmann et al. 2007; Sueur, Deneubourg, et al. 2011). Note that the classification of others as “friends” or not is done every time an individual is activated, and it is based on the accumulation of grooming interactions until that point. Thus, individuals considered “friends” at one point during the simulation may not be the same at a point later in the simulation. During the period of time in which we collect data from the model, however, friends remain stable and change little (Puga-Gonzalez et al. 2015). Grouping rules FriendsWorld has the same architecture and behavioral rules as GrooFiWorld, the only exception is that in FriendsWorld model individuals follow their “friends” (Puga-Gonzalez et al. 2015). Individuals have 3 different visual ranges, PersSpace, NearView, and MaxView. When an individual is activated, it searches for others within its PersSpace, if it does not perceive another in its close proximity, it acts according to the grouping rules. If an individual perceives one of its “friends” within its NearView, it will move one step towards it. If several “friends” are perceived, the individual moves towards the closest one. If none of its “friends” is perceived but others are, the individual just keeps on moving without modifying its original direction. When no others are perceived within NearView, the individual looks further away into MaxView. If other individuals are perceived within MaxView, the individual moves towards the closest “friend” if available; otherwise, it moves towards the closest group-member. If no individual is perceived within MaxView, the individual scans for others by turning over a SearchAngle. Social interactions In FriendsWorld, the social interaction rules are the same as in GrooFiWorld. If the individual perceives another one within its PersSpace, a dominance interaction may occur. Whether or not an individual chooses to fight depends on the outcome of a mental battle (Equation 2). If the individual wins the mental battle, a dominance interaction occurs (see “Dominance interactions” above). However, if the individual loses the mental battle, it may groom its partner depending on its level of Anxiety (see “Grooming interactions” above). Reaper Model This model is not spatially explicit and the selection of interaction partners is based solely on previous grooming interactions with other group-members; therefore, as already mentioned, we called it Reaper, i.e., reinforcement of affiliative preferences (Puga-Gonzalez and Sueur 2017a). Selection of interaction partners In contrast to GrooFiWorld, in this model individuals are not located in space and thus proximity has no influence on the distribution of their social interactions. At the start of the simulation, individuals are attributed a random waiting time drawn from a uniform distribution and the individual with the shortest waiting time gets activated first. Once activated, the probability of selecting a given group member as an interaction partner depends on the amount of grooming given and received exchanged with that specific partner. This probability is calculated according to Equation 4. Where Pij is the probability of individual i interacting with partner j; Gij is the number of grooming bouts exchanged between individual i and j; and the denominator is the sum of the number of grooming bouts individual i has exchanged with every other group member. α determines the degree of nonlinearity in the probability of selecting a given partner, the higher the value of α, the higher the tendency of individuals to interact with their most frequent grooming partners. This value of was set to 1.2. We choose this value of α because with it individuals allocated >75% of their grooming time to their top 5 interaction partners, as commonly observed in primate societies (Sueur, Deneubourg, et al. 2011)—for a full sensitivity analysis of the values of α, see Puga-Gonzalez and Sueur (2017a). At the beginning of the simulation, individuals have not yet interacted with any other group member and thus, all individuals have the same probability of being chosen as an interaction partner. As interactions go by, however, individuals will tend to select more frequently as interaction partners those individuals with whom they are more involved in social grooming. Note that Equation 4 does not imply reciprocal relationships, i.e., if individual i prefers to interact mostly with individual j, this does not necessarily mean that individual j also prefers to interact with individual of i. After selecting an interaction partner, individuals follow the same interactions rules as in GrooFiWorld. They first choose on the basis of a “mental fight” whether they perform or not a dominance interaction, if they choose not to fight, then they decide whether or not to groom (see dominance and grooming interactions above). Pij=Gijα∑a≠iN−1Giaα (4) Null-model In order to compare our results to a “null-model” model, we built a model in which we omit the effects of spatial proximity and preferential interactions with specific group members. In this model, individuals select interaction partners at random. Once an individual has randomly selected an interaction partner, a random number is drawn from a uniform distribution between [0, 1]. If this number is ≤ 0.5 the individual performs a dominance interaction, otherwise it grooms. Note that in this case individuals do not perform a “mental fight” to decide whether or not to attack, if the random number is ≤ to 0.5 individuals immediately attack. If the random number is > 0.5 individuals groom their partner independently of its value of anxiety. As in GrooFiWorld, individuals’ activations are regulated by attributing each individual a random waiting time drawn from a uniform distribution. The individual with the shortest waiting time gets activated first. Parametrization of the models In order to compare the empirical grooming networks with those of the models, in the simulations we mimicked the same group size and densities of each of the empirical networks. We did so because comparisons between networks differing in size and density (the number of observed edges divided by the number of possible edges) appear to be meaningless (Rankin et al. 2016). Indeed many network measures are dependent on these 2 parameters (Kasper and Voelkl 2009; Sueur, Petit, et al. 2011; Pasquaretta et al. 2014). Group size consisted of all adult individuals in the group; sex of individuals was also taken into consideration. Note that even though we only analyzed the behavior of females, we also included males in the simulations since male behavior may influence the behavior of real (Hemelrijk et al. 2008b) and computer simulated females (Hemelrijk et al. 2008a; Hemelrijk and Puga-Gonzalez 2012). In this way, we compared networks of females taken into consideration the presence of males in the group. To obtain female grooming networks with the same density as that observed in empirical data, first we collected data for long periods of time (up to 10,000 grooming interactions in some cases). Then, we constructed matrices in which the number of grooming interactions was increased sequentially following the order of data collection until the target density (± 2.5%) of the network was matched. Further, we also mimicked the grade of social style of each of the species of macaques. From previous analyses of the models, we knew that values of intensity of aggression of 1.0, 0.8, and 0.1, 0.08 for males and females, respectively, reproduce patterns that resemble those of intolerant and tolerant macaques, respectively (Hemelrijk and Puga-Gonzalez 2012; Puga-Gonzalez et al. 2015; Puga-Gonzalez and Sueur 2017a). Hence values for grades 1 and 4 of social style were set to 1.0, 0.8, and 0.1, 0.08 for males and females, respectively. Then, the values for grades 2 and 3 were chosen so they laid in between grades 1 and 4. In this way, we created 4 different sets of values of intensity of aggression “mimicking” the 4 different grades of social style described by Thierry (2004). These values were: grade 4, 0.1, 0.08; grade 3, 0.4, 0.32; grade 2, 0.7, 0.56; and grade one, 1.0, 0.8; for males and females respectively. The values of intensity of aggression assigned to each of the different groups are shown in Supplementary Table SI3. Preliminary analysis and optimization of the distribution of grooming in the models After running the models with the values of intensity of aggression described in the section “Parametrization of the models” (Supplementary Table SI3), we realized that the models did not match the presence/absence of some grooming patterns (grooming directed up the dominance hierarchy and/or towards individuals of similar rank) observed in the empirical groups. The presence/absence of these patterns may influence the value of some networks’ metrics. For instance, whether or not individuals direct grooming up the dominance hierarchy affects the correlation dominance rank—eigenvector and overall network centralization. Therefore, we decided to optimize the fit between the models’ grooming patterns and those observed in each group of empirical data. To optimize these patterns, we adjusted the parameter intensity of aggression (StepDom, Supplementary Table SI2). We specifically focused on this parameter because we knew it affects how grooming is distributed across group members whereas the other model’s parameters mainly affect cohesiveness of the group and/or the overall rate of grooming (Hemelrijk et al. 2017). The higher the value of intensity of aggression, the steeper the dominance hierarchy, the more individuals direct grooming up the hierarchy, at partners of similar rank and the more pronounced the spatial structure with dominant individuals at the center (Puga-Gonzalez et al. 2009). Hence, to match the presence/absence of these grooming patterns, we increased/decreased the values of intensity of aggression in the models. The new values of intensity of aggression assigned to each model and simulated group are shown in Supplementary Table SI8. Note that here we only show the results after the optimization of the grooming patterns. The results without optimization are found in the supporting information (see text therein, Supplementary Tables SI5–SI7; Supplementary SI12–SI14, and Supplementary Figure SI1–SI4). Data collection of the models For each empirical group (n = 14) and model (n = 4), we ran 20 different replicates, thus we ran a total of 1120 simulation. Each replicate consisted of up to 300 periods and each period consisted of (Group_Size * 20) individual’s activations. Data were collected from period 200 onwards to exclude any bias caused by transient values (Hemelrijk 1999). Data collection consisted of every social interaction: dyadic fights and grooming behavior. We recorded, for dominance interactions the identities of 1) the attacker and its opponent, 2) that of the winner/loser, 3) the updated Dom values of the individuals; and for grooming interactions, the identities of 4) groomer, 5) groomee, and 6) the updated anxiety values. From these interactions, we built grooming and aggression matrices. Note that the grooming matrices matched approximately the same density as that observed in empirical groups. Analyses of socio-behavioral patterns in the models For reasons of comparison with empirical data, the analyses of social behavior were performed only among females. As in empirical data we checked for grooming reciprocation, grooming up the hierarchy, and grooming towards partners of similar rank by means of the matrix Tau-Kr correlation (Hemelrijk 1990b). The level of significance was calculated using 2000 permutations and P-value < 0.05. The results shown are the average of 20 replicates with their combined probability using the improved Bonferroni procedure (Hochberg 1988). Dominance hierarchies were stable, and dominance ranks were significantly correlated (Kendall correlation every 30 periods, e.g., 200–230: tau > 0.88**). The asymmetry of dominance ranks was measured by the coefficient of variation of Dominance values (standard deviation divided by the mean). As in our previous studies, the higher the intensity of aggression the higher the asymmetry of dominance ranks (Hemelrijk and Puga-Gonzalez 2012; Puga-Gonzalez et al. 2015; Puga-Gonzalez and Sueur 2017a). Analysis of social networks of the models Like the analyses of social behavior, analyses were performed only among females. We performed the same social networks analysis as that performed in the social networks of the empirical groups (see above for metrics and definitions, sections “Individual level metrics” and “Group level metrics”). Note, however, that in the social networks derived from the models it was not necessary to correct for observation time since we kept track of all individuals’ interactions. Thus, network analyses were performed on directed binary/weighted social networks, depending on the network metric analyzed (see sections “Individual level metrics” and “Group level metrics”). In weighted networks, the weights represented the absolute grooming frequency given to another individual. Results of the global metrics shown are values of each of the 20 replicates per macaque group simulated. Likewise, results of the individuals’ metrics shown are values per individual in each of the 20 replicates per macaque group simulated. Thus, for each empirical data point (global or individual metric) there were 20 data points from each model. Comparison of social networks and model selection We performed linear regressions with the values of the observed metrics as the response variable, and the values of the metrics obtained from the simulated networks as the predictor variable. At the individual level, values were paired by sorting individuals by rank (within each group); at the global level, they were paired by groups. We compared effect sizes (R-squared), i.e., the percentage of variation explained by each model. However, to select which model fitted the observed data at best, we used the Akaike Information Criteria (AIC). We did so because we reasoned that if the social networks metrics from the models had a perfect fit with the observed metrics, then the observed and simulated values should be identical. In other words, if we plotted the values of the observed metrics (x axis) against the values from the model (y axis) and there was a perfect correspondence, then we would expect the values to lie along the diagonal with slope = 1 and intercept = 0 (note that the slope of the linear regressions was not always close to one). The AIC was calculated using the residual sum of squares (RSS) from this expectation according to Equation 5. AIC=n*ln(RSSn)+2*K (5) Where n is the number of data points, ln is the natural logarithm, RSS the residual sum of squares, and K the number of parameters in the model which was set to 2 (the model + error). The best model was determined by calculating the difference (Δ) between the model with the lowest AIC and all the other models, and by calculating the Akaike weight for each model (Equation 6). Where wi is the Akaike weight for model i, Δi is the difference between the AIC of the best fitting model and that of model i, and the denominator is the sum of the relative likelihood for all candidate models. Models with Δ values ≤2 were considered to have the same explanatory power. wi=exp(−0.5*Δi)∑r=1Rexp(−0.5*Δr) (6) RESULTS Grooming interactions and networks’ density The number of interactions necessary to match the target empirical density varied between the different models -GrooFiWorld, FriendsWorld, Reaper, and Null-model (Supplementary Table SI4). The null-model required the least number of grooming interactions to match the empirical network density. As regards the other 3 models, no apparent difference among them was found when the target density was below 60% and intensity of aggression was set to low (Supplementary Table SI4). However, for densities above 60% the Reaper model took many more grooming interactions than the two spatially explicit models to reach the target density. More interestingly, in all models except the null one, simulations belonging to some specific groups never reached a network density of 100% regardless of the number of grooming interactions. In these cases, the density reached a plateau around values of ~80%. In the spatially explicit models (GrooFiWorld and FriendsWorld), this only happened when intensity of aggression was set to too high values. When intensity of aggression was high the asymmetry of the dominance values among females was great and thus, when females at the top of the hierarchy interacted with those at the bottom, they chose to fight instead of grooming. In the Reaper model, however, this happened regardless of the values of intensity of aggression. In this case, females became so selective on their interaction partners that they never selected other group members as interaction partners. Patterns of grooming Reciprocation of social grooming In all groups of empirical data, grooming given was positively correlated with grooming received except in the 2 groups of M. tonkeana in which no significant correlation was found (Supplementary Table SI5). Concerning the computer models, regardless of the simulated group, we always found a significant positive matrix correlation between grooming given and received, except in the null-model. Correlations were never significant in this model (Supplementary Table SI5). Optimization of distribution of grooming We modified the values of intensity of aggression (StepDom) to optimize the match between the patterns emerging in the model (grooming directed up the dominance hierarchy and towards individuals of similar rank) and those observed in the groups of empirical data (see Methods and Supplementary Table SI8). Here we only present the results after optimization, for the results without optimization see supporting information (text therein, Supplementary Tables SI5–SI7; Supplementary Tables SI12–SI14, and Supplementary Figures SI1–SI4). Regarding grooming directed up the dominance hierarchy, after optimization of the parameters (Supplementary Table SI8), GrooFiWorld matched all observations in empirical data (14/14); Reaper model matched all except one group (13/14); and FriendsWorld did not match 3 groups (11/14); in the null-model, grooming up the dominance hierarchy never reached statistical significance (Supplementary Table S10). Regarding grooming directed towards individuals of similar rank, the optimization of parameters improved the match between models and empirical data, but there were still several mismatches. GrooFiWorld did not match 4 groups (10/14); FriendsWorld 5 groups (9/14); and the Reaper 7 groups (7/14). Note that in the Reaper and null models this pattern never reached significance; thus the cases that matched empirical data were those in which this pattern was absent (Supplementary Table S11). The reason for the mismatches is that in GrooFiWorld and FriendsWorld grooming up the hierarchy and directed towards partners of similar rank usually emerge (or not) together. In empirical data, however, this was not always the case. In the group of M. assamensis and M. radiata, grooming was directed up the hierarchy but not towards partners of similar rank. Conversely, in 2 groups of M. fuscata and one of M. fascicularis, grooming was directed towards partners of similar rank but not up the dominance hierarchy (Supplementary Tables S10–11). SNA Because the condition in which animals lived, captive or wild (free-ranging), can affect the way individuals interact with each other, we decided to compare the network metrics by dividing the empirical groups into wild/provisioned/free-ranging (n = 8) and captive (n = 6) (Supplementary Table SI1). Condition seemed to have an effect on whether the model matched or not the empirical data, particularly for the global metrics of modularity, eigenvector, correlation dominance—eigenvector, and centralization index (Supplementary Figures SI9–SI12). All models matched better these metrics when the simulated group was wild rather than captive (Supplementary Figures SI11–SI12). This is evident from the much higher R2 in wild than in the captive groups (Supplementary Table SI18) and from the negative slopes of the regression lines for the captive groups, i.e., opposite to expected (Supplementary Figure SI12). Thus, we only present here the comparison between the network metrics of the models with those of wild groups of macaques. For the network metric comparison without optimization of the grooming patterns, and for the comparison after optimizing the grooming patterns but without dividing groups into wild and captive, see supporting information. Individual metrics With regard to the individual network metrics of closeness, betweenness, and eigenvector all computer models seemed to predict the observed female’s values, particularly those of closeness (Figure 1). Effect sizes were similar among all models, explaining >30% of the variance of observed values of clustering coefficient; >60% of the variance of betweenness and eigenvector; and > 97% of the variance in values of closeness (Table 2). FriendsWorld had the highest effect size for betweenness; FriendsWorld and Reaper for closeness; and Reaper for eigenvector and clustering coefficient (Table 2). After calculating the Akaike information criteria (AIC), it appeared that patterns generated in FriendsWorld fitted best the empirical values of betweenness and closeness, and Reaper those of eigenvector and clustering coefficient (Table 2). In all cases, the null-model was the worst fitting model (Table 2). Figure 1 View largeDownload slide Plots of the values of the observed individual’s metric (x axis) against the values predicted by each theoretical model (y axis). Black diagonal represents values expected in a perfect match between model and empirical data. Values of betweenness centrality and closeness are shown in left and right columns, respectively. Colors represent macaque groups and each data point represents a single individual. Figure 1 View largeDownload slide Plots of the values of the observed individual’s metric (x axis) against the values predicted by each theoretical model (y axis). Black diagonal represents values expected in a perfect match between model and empirical data. Values of betweenness centrality and closeness are shown in left and right columns, respectively. Colors represent macaque groups and each data point represents a single individual. Table 2 Linear regressions between observed and predicted individuals’ network metric for each model, and its correspondent Akaike information criteria (AIC) (for details of the calculation, see Methods) Linear regression Akaike information Metric Model β R2 P-value AIC Delta (Δ) Weight Closeness FriendsWorld 0.97 0.98 <0.001 −21,120 0 1 Reaper 0.99 0.98 <0.001 −20,792 327 0 GrooFiWorld 0.94 0.97 <0.001 −20,131 989 0 Null 0.97 0.97 <0.001 −19,509 1611 0 Betweenness FriendsWorld 1.38 0.73 <0.001 10,090 0 1 Reaper 1.73 0.68 <0.001 10,730 640 0 GrooFiWorld 1.47 0.63 <0.001 10,745 655 0 Null 1.95 0.65 <0.001 10,989 899 0 Eigen Centrality Reaper 0.92 0.69 <0.001 −6611 0 1 FriendsWorld 0.80 0.65 <0.001 −5996 616 0 GrooFiWorld 0.77 0.61 <0.001 −5820 791 0 Null 1.11 0.60 <0.001 −4645 1966 0 Clustering Coeff. Reaper 0.58 0.33 <0.001 −8761 0 1 GrooFiWorld 0.42 0.29 <0.001 −8178 583 0 FriendsWorld 0.44 0.31 <0.001 −7875 886 0 Null 0.30 0.09 <0.001 −7013 1748 0 Linear regression Akaike information Metric Model β R2 P-value AIC Delta (Δ) Weight Closeness FriendsWorld 0.97 0.98 <0.001 −21,120 0 1 Reaper 0.99 0.98 <0.001 −20,792 327 0 GrooFiWorld 0.94 0.97 <0.001 −20,131 989 0 Null 0.97 0.97 <0.001 −19,509 1611 0 Betweenness FriendsWorld 1.38 0.73 <0.001 10,090 0 1 Reaper 1.73 0.68 <0.001 10,730 640 0 GrooFiWorld 1.47 0.63 <0.001 10,745 655 0 Null 1.95 0.65 <0.001 10,989 899 0 Eigen Centrality Reaper 0.92 0.69 <0.001 −6611 0 1 FriendsWorld 0.80 0.65 <0.001 −5996 616 0 GrooFiWorld 0.77 0.61 <0.001 −5820 791 0 Null 1.11 0.60 <0.001 −4645 1966 0 Clustering Coeff. Reaper 0.58 0.33 <0.001 −8761 0 1 GrooFiWorld 0.42 0.29 <0.001 −8178 583 0 FriendsWorld 0.44 0.31 <0.001 −7875 886 0 Null 0.30 0.09 <0.001 −7013 1748 0 The best fitting models according to the AIC weight are shown in bold. Note that for each metric the models are ordered according to increasing values of delta (Δ). View Large Table 2 Linear regressions between observed and predicted individuals’ network metric for each model, and its correspondent Akaike information criteria (AIC) (for details of the calculation, see Methods) Linear regression Akaike information Metric Model β R2 P-value AIC Delta (Δ) Weight Closeness FriendsWorld 0.97 0.98 <0.001 −21,120 0 1 Reaper 0.99 0.98 <0.001 −20,792 327 0 GrooFiWorld 0.94 0.97 <0.001 −20,131 989 0 Null 0.97 0.97 <0.001 −19,509 1611 0 Betweenness FriendsWorld 1.38 0.73 <0.001 10,090 0 1 Reaper 1.73 0.68 <0.001 10,730 640 0 GrooFiWorld 1.47 0.63 <0.001 10,745 655 0 Null 1.95 0.65 <0.001 10,989 899 0 Eigen Centrality Reaper 0.92 0.69 <0.001 −6611 0 1 FriendsWorld 0.80 0.65 <0.001 −5996 616 0 GrooFiWorld 0.77 0.61 <0.001 −5820 791 0 Null 1.11 0.60 <0.001 −4645 1966 0 Clustering Coeff. Reaper 0.58 0.33 <0.001 −8761 0 1 GrooFiWorld 0.42 0.29 <0.001 −8178 583 0 FriendsWorld 0.44 0.31 <0.001 −7875 886 0 Null 0.30 0.09 <0.001 −7013 1748 0 Linear regression Akaike information Metric Model β R2 P-value AIC Delta (Δ) Weight Closeness FriendsWorld 0.97 0.98 <0.001 −21,120 0 1 Reaper 0.99 0.98 <0.001 −20,792 327 0 GrooFiWorld 0.94 0.97 <0.001 −20,131 989 0 Null 0.97 0.97 <0.001 −19,509 1611 0 Betweenness FriendsWorld 1.38 0.73 <0.001 10,090 0 1 Reaper 1.73 0.68 <0.001 10,730 640 0 GrooFiWorld 1.47 0.63 <0.001 10,745 655 0 Null 1.95 0.65 <0.001 10,989 899 0 Eigen Centrality Reaper 0.92 0.69 <0.001 −6611 0 1 FriendsWorld 0.80 0.65 <0.001 −5996 616 0 GrooFiWorld 0.77 0.61 <0.001 −5820 791 0 Null 1.11 0.60 <0.001 −4645 1966 0 Clustering Coeff. Reaper 0.58 0.33 <0.001 −8761 0 1 GrooFiWorld 0.42 0.29 <0.001 −8178 583 0 FriendsWorld 0.44 0.31 <0.001 −7875 886 0 Null 0.30 0.09 <0.001 −7013 1748 0 The best fitting models according to the AIC weight are shown in bold. Note that for each metric the models are ordered according to increasing values of delta (Δ). View Large Global metrics All models seemed to accurately predict the observed values of number of triangles and clustering coefficient (Figure 3). Effect sizes were highest for GrooFiWorld (R2 = 0.82) for number of triangles, and Reaper and FriendsWorld (R2 = 0.89) for clustering coefficient (Table 3). According to the AIC, GroFiWorld and Reaper were the best fitting models for number of triangles and clustering coefficient, respectively (Table 3). Regarding diameter and the correlation dominance—eigenvector, all models, except the null one, moderately predicted the observed values (Figures 3 and 4). Effect sizes were highest for FriendsWorld (R2 = 0.49) for diameter and GrooFiWorld (R2 = 0.43) for the correlation dominance—eigenvector, and according to the AIC, these models were also the best fitting ones (Table 3). With respect to modularity of the network, the match between observed and predicted data was not as good as for the above metrics; and in the case of Reaper model, modularity seemed to not vary but to remain constant (Figure 4). Similarly, the values of centrality index were poorly matched by all models (Figure 4). Effect sizes were highest for GrooFiWorld (R2 = 0.34) and Reaper model (R2 = 0.36) for modularity and centrality index respectively (Table 3). And according to the AIC, the best fitting model were also GrooFiWorld and Reaper (Table 3). The average eigenvector value of the group was the network metric matched at worst by the models, all of them overestimated this value (Figure 4). Effect size was highest for FriendsWorld (R2 = 0.25), but according to the AIC, the Reaper model seemed to fit the data at best (Table 3). As for the individual metrics, in all cases, except for clustering coefficient, the null model was the worst fitting model (Table 3). Table 3 Linear regression between observed and predicted global network metric for each model, and its correspondent Akaike information (AIC) (for details of the calculation, see methods) Linear Regression Akaike information Metric Model Β R2 P-value AIC Delta (Δ) Weight Num Triangles GrooFiWorld 0.92 0.82 <0.001 1085 0 1 Reaper 0.82 0.81 <0.001 1124 39 0 FriendsWorld 0.81 0.81 <0.001 1140 55 0 Null 0.62 0.73 <0.001 1308 223 0 Diameter FriendsWorld 0.86 0.49 <0.001 52 0 1 Reaper 0.95 0.41 <0.001 77 25 0 GrooFiWorld 0.82 0.41 <0.001 82 30 0 Null 1.07 0.34 <0.001 105 53 0 Clustering Coeff. Reaper 1.00 0.89 <0.001 −780 0 1 FriendsWorld 1.28 0.89 <0.001 −657 123 0 Null 1.07 0.83 <0.001 −652 129 0 GrooFiWorld 1.17 0.81 <0.001 −636 144 0 Modularity GrooFiWorld 0.91 0.34 <0.001 −508 0 1 FriendsWorld 0.93 0.28 <0.001 −491 17 0 Reaper 0.13 0.00 0.627 −482 26 0 Null 1.74 0.28 <0.001 −389 118 0 Eigenvector Reaper 0.58 0.22 <0.001 −588 0 1 FriendsWorld 0.50 0.25 <0.001 −519 69 0 GrooFiWorld 0.36 0.14 <0.001 −496 93 0 Null 0.47 0.16 <0.001 −384 205 0 Corr Dom-Eigen GrooFiWorld 0.58 0.43 <0.001 −369 0 1 FriendsWorld 0.59 0.39 <0.001 −352 17 0 Reaper 0.44 0.17 <0.001 −349 20 0 Null 0.02 −0.01 0.779 −163 206 0 Centrality Index Reaper 0.86 0.36 <0.001 1057 0 1 FriendsWorld 0.74 0.26 <0.001 1133 76 0 GrooFiWorld 0.45 0.10 <0.001 1162 105 0 Null 0.60 0.15 <0.001 1263 206 0 Linear Regression Akaike information Metric Model Β R2 P-value AIC Delta (Δ) Weight Num Triangles GrooFiWorld 0.92 0.82 <0.001 1085 0 1 Reaper 0.82 0.81 <0.001 1124 39 0 FriendsWorld 0.81 0.81 <0.001 1140 55 0 Null 0.62 0.73 <0.001 1308 223 0 Diameter FriendsWorld 0.86 0.49 <0.001 52 0 1 Reaper 0.95 0.41 <0.001 77 25 0 GrooFiWorld 0.82 0.41 <0.001 82 30 0 Null 1.07 0.34 <0.001 105 53 0 Clustering Coeff. Reaper 1.00 0.89 <0.001 −780 0 1 FriendsWorld 1.28 0.89 <0.001 −657 123 0 Null 1.07 0.83 <0.001 −652 129 0 GrooFiWorld 1.17 0.81 <0.001 −636 144 0 Modularity GrooFiWorld 0.91 0.34 <0.001 −508 0 1 FriendsWorld 0.93 0.28 <0.001 −491 17 0 Reaper 0.13 0.00 0.627 −482 26 0 Null 1.74 0.28 <0.001 −389 118 0 Eigenvector Reaper 0.58 0.22 <0.001 −588 0 1 FriendsWorld 0.50 0.25 <0.001 −519 69 0 GrooFiWorld 0.36 0.14 <0.001 −496 93 0 Null 0.47 0.16 <0.001 −384 205 0 Corr Dom-Eigen GrooFiWorld 0.58 0.43 <0.001 −369 0 1 FriendsWorld 0.59 0.39 <0.001 −352 17 0 Reaper 0.44 0.17 <0.001 −349 20 0 Null 0.02 −0.01 0.779 −163 206 0 Centrality Index Reaper 0.86 0.36 <0.001 1057 0 1 FriendsWorld 0.74 0.26 <0.001 1133 76 0 GrooFiWorld 0.45 0.10 <0.001 1162 105 0 Null 0.60 0.15 <0.001 1263 206 0 The best fitting models according to the AIC weight are shown in bold. Note that for each metric the models are ordered according to increasing values of delta (Δ). View Large Table 3 Linear regression between observed and predicted global network metric for each model, and its correspondent Akaike information (AIC) (for details of the calculation, see methods) Linear Regression Akaike information Metric Model Β R2 P-value AIC Delta (Δ) Weight Num Triangles GrooFiWorld 0.92 0.82 <0.001 1085 0 1 Reaper 0.82 0.81 <0.001 1124 39 0 FriendsWorld 0.81 0.81 <0.001 1140 55 0 Null 0.62 0.73 <0.001 1308 223 0 Diameter FriendsWorld 0.86 0.49 <0.001 52 0 1 Reaper 0.95 0.41 <0.001 77 25 0 GrooFiWorld 0.82 0.41 <0.001 82 30 0 Null 1.07 0.34 <0.001 105 53 0 Clustering Coeff. Reaper 1.00 0.89 <0.001 −780 0 1 FriendsWorld 1.28 0.89 <0.001 −657 123 0 Null 1.07 0.83 <0.001 −652 129 0 GrooFiWorld 1.17 0.81 <0.001 −636 144 0 Modularity GrooFiWorld 0.91 0.34 <0.001 −508 0 1 FriendsWorld 0.93 0.28 <0.001 −491 17 0 Reaper 0.13 0.00 0.627 −482 26 0 Null 1.74 0.28 <0.001 −389 118 0 Eigenvector Reaper 0.58 0.22 <0.001 −588 0 1 FriendsWorld 0.50 0.25 <0.001 −519 69 0 GrooFiWorld 0.36 0.14 <0.001 −496 93 0 Null 0.47 0.16 <0.001 −384 205 0 Corr Dom-Eigen GrooFiWorld 0.58 0.43 <0.001 −369 0 1 FriendsWorld 0.59 0.39 <0.001 −352 17 0 Reaper 0.44 0.17 <0.001 −349 20 0 Null 0.02 −0.01 0.779 −163 206 0 Centrality Index Reaper 0.86 0.36 <0.001 1057 0 1 FriendsWorld 0.74 0.26 <0.001 1133 76 0 GrooFiWorld 0.45 0.10 <0.001 1162 105 0 Null 0.60 0.15 <0.001 1263 206 0 Linear Regression Akaike information Metric Model Β R2 P-value AIC Delta (Δ) Weight Num Triangles GrooFiWorld 0.92 0.82 <0.001 1085 0 1 Reaper 0.82 0.81 <0.001 1124 39 0 FriendsWorld 0.81 0.81 <0.001 1140 55 0 Null 0.62 0.73 <0.001 1308 223 0 Diameter FriendsWorld 0.86 0.49 <0.001 52 0 1 Reaper 0.95 0.41 <0.001 77 25 0 GrooFiWorld 0.82 0.41 <0.001 82 30 0 Null 1.07 0.34 <0.001 105 53 0 Clustering Coeff. Reaper 1.00 0.89 <0.001 −780 0 1 FriendsWorld 1.28 0.89 <0.001 −657 123 0 Null 1.07 0.83 <0.001 −652 129 0 GrooFiWorld 1.17 0.81 <0.001 −636 144 0 Modularity GrooFiWorld 0.91 0.34 <0.001 −508 0 1 FriendsWorld 0.93 0.28 <0.001 −491 17 0 Reaper 0.13 0.00 0.627 −482 26 0 Null 1.74 0.28 <0.001 −389 118 0 Eigenvector Reaper 0.58 0.22 <0.001 −588 0 1 FriendsWorld 0.50 0.25 <0.001 −519 69 0 GrooFiWorld 0.36 0.14 <0.001 −496 93 0 Null 0.47 0.16 <0.001 −384 205 0 Corr Dom-Eigen GrooFiWorld 0.58 0.43 <0.001 −369 0 1 FriendsWorld 0.59 0.39 <0.001 −352 17 0 Reaper 0.44 0.17 <0.001 −349 20 0 Null 0.02 −0.01 0.779 −163 206 0 Centrality Index Reaper 0.86 0.36 <0.001 1057 0 1 FriendsWorld 0.74 0.26 <0.001 1133 76 0 GrooFiWorld 0.45 0.10 <0.001 1162 105 0 Null 0.60 0.15 <0.001 1263 206 0 The best fitting models according to the AIC weight are shown in bold. Note that for each metric the models are ordered according to increasing values of delta (Δ). View Large DISCUSSION On basis of the effect sizes and the AIC, we could not distinguish one model (GrooFiWorld, FriendsWorld, or Reaper) fitting all network metrics better than the others. Each of these models fitted the empirical data in at least 3 network metrics better than the others (Table 4). The null model, on the other hand, performed worst for all metrics except one (Table 4). Patterns generated by Reaper and FriendsWorld seemed to fit the empirical network metrics slightly better than GrooFiWorld which had the worst fit among the non-null models in 7 out of 11 cases (Table 4). Regarding grooming reciprocation, directed up the dominance hierarchy, and towards individuals of similar rank, GrooFiWorld performed better, followed by FriendsWorld and Reaper models (Supplementary Tables SI9–SI11). Similarly, when the analysis were performed without optimization and/or without dividing the data into wild and captive, the overall ranks of the models changed slightly and no model fitted all patterns better than the others (Supplementary Tables SI14 and SI17). Our analysis thus does not allow to establish that one model was better than other. If anything, it indicates that both space, i.e., proximity-based reciprocity as modeled in GrooFiWorld and FriendsWorld, and social bonding, i.e., emotional bookkeeping as modeled in Reaper, are relevant factors influencing the distribution of social grooming. Table 4 Rank of the models according to the AIC weight and Delta (Δ) value for each network metric Computer Model Network metric FriendsWorld Groofiworld Reaper Null-model Individual: 1) Closeness 1 3 2 4 2) Betweenness 1 3 2 4 3) Eigenvector 2 3 1 4 4) Clustering coefficient 2.5 2.5 1 4 Global: 5) Number of triangles 3 1 2 4 6) Diameter 1 3 2 4 7) Average clustering coefficient 2 4 1 3 8) Modularity 2 1 3 4 9) Average eigenvector centrality 2 3 1 4 10) Corr. eigenvector - dominance rank 2 1 3 4 11) Centrality Index 2 3 1 4 Average Ranking 2.0 2.5 1.7 3.9 Computer Model Network metric FriendsWorld Groofiworld Reaper Null-model Individual: 1) Closeness 1 3 2 4 2) Betweenness 1 3 2 4 3) Eigenvector 2 3 1 4 4) Clustering coefficient 2.5 2.5 1 4 Global: 5) Number of triangles 3 1 2 4 6) Diameter 1 3 2 4 7) Average clustering coefficient 2 4 1 3 8) Modularity 2 1 3 4 9) Average eigenvector centrality 2 3 1 4 10) Corr. eigenvector - dominance rank 2 1 3 4 11) Centrality Index 2 3 1 4 Average Ranking 2.0 2.5 1.7 3.9 Models with a Delta (Δ) value difference ≤ 2 were given the same rank. View Large Table 4 Rank of the models according to the AIC weight and Delta (Δ) value for each network metric Computer Model Network metric FriendsWorld Groofiworld Reaper Null-model Individual: 1) Closeness 1 3 2 4 2) Betweenness 1 3 2 4 3) Eigenvector 2 3 1 4 4) Clustering coefficient 2.5 2.5 1 4 Global: 5) Number of triangles 3 1 2 4 6) Diameter 1 3 2 4 7) Average clustering coefficient 2 4 1 3 8) Modularity 2 1 3 4 9) Average eigenvector centrality 2 3 1 4 10) Corr. eigenvector - dominance rank 2 1 3 4 11) Centrality Index 2 3 1 4 Average Ranking 2.0 2.5 1.7 3.9 Computer Model Network metric FriendsWorld Groofiworld Reaper Null-model Individual: 1) Closeness 1 3 2 4 2) Betweenness 1 3 2 4 3) Eigenvector 2 3 1 4 4) Clustering coefficient 2.5 2.5 1 4 Global: 5) Number of triangles 3 1 2 4 6) Diameter 1 3 2 4 7) Average clustering coefficient 2 4 1 3 8) Modularity 2 1 3 4 9) Average eigenvector centrality 2 3 1 4 10) Corr. eigenvector - dominance rank 2 1 3 4 11) Centrality Index 2 3 1 4 Average Ranking 2.0 2.5 1.7 3.9 Models with a Delta (Δ) value difference ≤ 2 were given the same rank. View Large All models fitted relatively well the observed values of the individual metrics of betweenness, closeness, and eigenvector. In these cases, effect sizes were usually high and data points laid usually along the diagonal (Figures 1 and 2, Table 2). This was the case even for the null-model, which indicates that the value of these individual metrics may be partially constrained by group size and network density (i.e., the network properties matched in the simulations). The models did not explain more than 33% of the variance of clustering coefficient. This indicates that the grooming cliques formed by individuals were different in the model and empirical data. In the models, clustering coefficient was mostly overestimated, i.e., individuals and their neighbors were more connected than in reality; suggesting then, that in reality individuals are more selective in choosing their interaction partners. Figure 2 View largeDownload slide Plots of the values of the observed individual’s metric (x axis) against the values predicted by each model (y axis). Black diagonal represents values expected in a perfect match between model and empirical data. Values of clustering coefficient and eigenvector are shown in left and right columns, respectively. Colors represent macaque groups and each data point represents a single individual. Figure 2 View largeDownload slide Plots of the values of the observed individual’s metric (x axis) against the values predicted by each model (y axis). Black diagonal represents values expected in a perfect match between model and empirical data. Values of clustering coefficient and eigenvector are shown in left and right columns, respectively. Colors represent macaque groups and each data point represents a single individual. With respect to the global metrics, none of the models performed as good as for the individual metrics (Figures 3 and 4). All models seemed to overestimate the values of eigenvector centrality and of the correlation coefficient dominance-eigenvector; and underestimate those of the centrality index (Figure 4). Nonetheless, the null-model was consistently worse than the other 3 models (Table 3), indicating that the global network metrics depend more on the patterning of social interactions produced by the mechanisms implemented in the models. On average, “Reaper” predicted the global network metrics’ values slightly better than the other models. With regard to modularity, however, it produced values that did not vary but remained constant; whereas in the spatially explicit models modularity varied according to the values observed in the empirical data (Tables 3 and 4, Figure 4). This suggests that a spatial structure is necessary to produce variation in modularity. Figure 3 View largeDownload slide Plots of the values of the observed global metric (x axis) against the values predicted by each model (y axis). Black diagonal represents values expected in a perfect match between model and empirical data. Colors represent different models and their corresponding trend lines. Data points represent the metric’s value of a single simulation (n = 20) from each model. CC = clustering coefficient. Figure 3 View largeDownload slide Plots of the values of the observed global metric (x axis) against the values predicted by each model (y axis). Black diagonal represents values expected in a perfect match between model and empirical data. Colors represent different models and their corresponding trend lines. Data points represent the metric’s value of a single simulation (n = 20) from each model. CC = clustering coefficient. Figure 4 View largeDownload slide Plots of the values of the observed global metric (x axis) against the values predicted by each model (y axis). Black diagonal represents values expected in a perfect match between model and empirical data. Colors represent different models and their corresponding trend lines. Data points represent the metric’s value of a single simulation (n = 20) from each model. Figure 4 View largeDownload slide Plots of the values of the observed global metric (x axis) against the values predicted by each model (y axis). Black diagonal represents values expected in a perfect match between model and empirical data. Colors represent different models and their corresponding trend lines. Data points represent the metric’s value of a single simulation (n = 20) from each model. With respect to grooming reciprocation, up the dominance hierarchy, and towards partners of similar rank, none of the models was able to fit all of the patterns across all empirical groups (Supplementary Tables SI9–SI11). “GrooFiWorld” failed to match 6 out of 42 patterns, “Reaper” 9, and “FriendsWorld” 10. It is noteworthy that grooming directed towards individuals of similar rank never reached statistical significance in “Reaper” (Supplementary Table SI11): the absence of this pattern can be explained by a lack of spatial structure. In “GrooFiWorld” and “FriendsWorld,” individuals of similar rank are often close due to the spatial structure; this induces frequent interactions among them, and the emergence of this pattern (Puga-Gonzalez et al. 2009). GrooFiWorld, FriendsWorld, and Reaper models reproduced reciprocation. We were however unable to firmly confirm which model, and thus which underlying mechanism, “proximity-based” or “emotional bookkeeping,” may reproduce more accurately the diversity of networks features and distribution of grooming observed in macaque groups. It could be argued that exploring the parameter space (all possible combinations of values for all models’ parameters, Supplementary Table SI2) may produce a better fit with empirical data. We think this is not the case. Models made poor predictions in particular for values of modularity, centrality index, and average eigenvector (Figure 4). Modularity and centrality index assess the degree of subgrouping and the popularity of the most central individual within the network, respectively; and eigenvector is a measurement of how well the individuals and their neighbors are connected within the network. To reproduce in the models a high degree of modularity and centrality index as observed in empirical data, it would be necessary that individuals, somehow, interact more with their subgroup than with the rest of the group, and that all or most individuals in the group have the same preferred interaction partner respectively. These patterns, thus, seem a consequence of the way individuals distribute grooming among their partners. The grouping, grooming, and most of the dominance parameters (Supplementary Table SI2), modify the rate of interactions, grooming frequency, and cohesiveness of the group but not the way individuals distribute their grooming (Hemelrijk and Puga-Gonzalez 2012; Hemelrijk et al. 2017). Hence, the lack of a better fit between the models and these global metrics indicates that an additional process(es), affecting the way individuals distribute their grooming, is missing in the computer models here analyzed. None of the mechanisms in these models (spatial structure and proximity-based interactions: GrooFiWorld; social bonding: Reaper; or their combination: FriendsWorld) is sufficient to produce attraction to specific individuals and subgrouping effects as significant as those observed in nature. In nonhuman primates, individuals appear to be knowledgeable not only of their own social relationships but also of those between third-parties (Silk 1999; Perry et al. 2004; Schino et al. 2006; Paxton et al. 2010; Kubenova et al. 2017). Thus, a potential mechanism missing in the models may be that individuals can choose partners, taking into account not only how other individuals behave towards them but also towards other group members. Implementing such a mechanism in the models could have an impact on triadic closure, and thus on the centralization and modularity of the network. Another potential mechanism is that, in nonhuman primates, the degree with which individuals reduce their anxiety after being involved in grooming depends on the identity of the partner. Individuals that groom with “friends” seem to relax more than when they groom with “nonfriends” (Crockford et al. 2013; Wittig et al. 2016). In the model, however, the reduction in anxiety after grooming was the same regardless of the identity of the partner. Further, individuals may take into account not only grooming exchange but also aggression exchange when establishing social bonds, whereas in “Friendsworld” and “Reaper” grooming exchange was the only variable taken into account. Furthermore, the model omitted other individual attributes such as age and kinship, which are known to influence social bonding and subgrouping. Further research will be necessary to know whether the implementation of these or other mechanisms in the models will produce a better fit between models’ and empirical grooming networks. The models fitted the data from wild (provisioned) groups better than that of captive ones. Indeed, in some cases with the captive data the predictions of the models were opposite of expected (SI11-SI12). The models, thus, seem to represent the behavior of wild animals better than those maintained in captive conditions (see R2 in Supplementary Table SI18). The lack of a better fit of the models with the latter may be due to the spatial constraint imposed by captivity on these groups. In the wild individuals can avoid others by moving away (Heesen et al. 2014; Heesen et al. 2015), whereas in captivity individuals have to interact with most group members. This may cause an artificial distribution of grooming which would ultimately affect grooming networks. In the spatially explicit models, although individuals followed rules that made them stay in the group (by searching for others and moving towards them when no other group member was in view), they were free to move in any direction since there were no spatial boundaries. Whether imposing a spatial boundary in these models would result in a better fit with empirical data from captive groups should be studied in the future. To fit the grooming patterns observed in some groups of macaques, in the models the parameter of aggression intensity had to be optimized (Supplementary Tables SI8). To simulate the group of Barbary macaques, for instance, aggression intensity was set to high values because females directed grooming up the dominance hierarchy and towards partners of similar rank in this group (Supplementary Tables SI6–SI7). Similarly, to simulate some groups of Japanese macaques, aggression intensity was set to low values because in none of these groups grooming was directed up the hierarchy (Supplementary Table SI6). According to the biological market theory, grooming directed up the dominance hierarchy is a pattern expected and usually observed in intolerant primate societies and not in tolerant ones (Barrett et al. 1999; Schino and Aureli 2008b). Some of the macaque groups here analyzed, thus, presented a pattern opposite to the one expected according to their species’ classification of social style. We do not believe, however, that this observation is sufficient to consider a reclassification of macaque species’ social style because 1) although from the groups here studied some presented a pattern opposite to the expected, other groups of the same species did conform to the patterns expected according to its social style, and 2) the classification of social styles is based not only on the distribution of social grooming but on other patterns like rates of reconciliation, directionality of aggression, or permissiveness of mothers (Thierry 2004). We have no indication that in our data set living conditions—wild or captive—have influenced the observed grooming patterns; actually, the patterns of 5 wild and 3 captive groups did not conform to the patterns previously documented for social style (but see Balasubramaniam et al. 2017 for a significant result on the influence of living condition on grooming traits). Using the principle of parsimony, the simplest model should be favored. In the models, however, the cognitive demands on individuals do not seem exceptionally high. In all of them, individuals should be capable of identifying other group-members and recognizing their dominance status. And in “FriendsWorld” and “Reaper” models individuals should additionally differentiate others on the basis of the grooming given and received, a task that could be achieved through associative learning. Even though cognitive demands on individuals seem lower in GrooFiWorld, nonhuman primates have been shown to possess these cognitive capacities anyhow (Tomasello and Call 1997). Thus, it is difficult to favor one model over the other on the basis of cognitive limitations. The models presented here are not the only ones that reproduce patterns of reciprocation of grooming similar to those observed in primate societies. In the partner-choice model by Schino and Campenni (Campenni and Schino 2014) and the EMO-model by Evers and co-authors (2015), preferential interactions alone or in combination with spatial structure, respectively, suffice to reproduce grooming reciprocation. Nonetheless, no analyses of grooming social networks, grooming directed up the dominance hierarchy or towards individuals of similar rank, and of interchange of grooming for support have been done on these models. Thus, it is unknown if in these models the distribution of grooming and/or the grooming social networks are similar to those reported in macaques. In conclusion, our results indicate that 3 models capture to some degree the diversity of social styles observed in the grooming patterns and networks of macaques. Our comparative analysis, however, does not allow us to make a clear distinction between the models. Whereas GrooFiWorld fits the distribution of grooming among group members best, “Reaper” and “FriendsWorld” do so for the network metrics. Thus, it is not possible to specify which mechanism, proximity-based or emotional book-keeping, is more accurate at reproducing reciprocation, grooming distribution and social networks as observed in macaques. What clearly appears, however, is that a process is missing in these models. This process could be a complex one, such as individuals taking into account third-party relationships when deciding whom to interact with. Whether such a mechanism can more accurately reproduce the diversity of social styles observed in macaque societies we will study in future models. SUPPLEMENTARY MATERIAL Supplementary data are available at Behavioral Ecology online. FUNDING This work was supported by a post-doctoral research grant from the National Council of Science and Technology (CONACyT) of Mexico awarded to I.P.G. Long-term research on Assamese macaques is supported by the German Primate Center and Georg-August University of Gottingen. Data accessibility: Analyses reported in this article can be reproduced using the data provided by Puga-Gonzalez et al. (2018). J.O. and O.S. thank the National Research Council of Thailand (NRCT) and the Department of National Parks, Wildlife and Plant Conservation (DNP) for permission to conduct this study, A. Koening and C. Borries who developed the field site at PKWS, and K. Nitaya T. Wongsnak, M. Pongjantarasatien, K. Kreetiyutanont, M. Kumsuk, and W. Saenphala for their cooperation. B.T. thanks Arianna de Marco for data collection. S.S. thanks Ellen Merz for allowing him to lead his research in La Forêt des Singes and for providing the demographic data. REFERENCES Aureli F , Cords M , Van Schaik CP . 2002 . Conflict resolution following aggression in gregarious animals: a predictive framework . Anim Behav . 64 : 325 – 343 . Google Scholar CrossRef Search ADS Aureli F , Preston SD , de Waal FB . 1999 . Heart rate responses to social interactions in free-moving rhesus macaques (Macaca mulatta): a pilot study . J Comp Psychol . 113 : 59 – 65 . Google Scholar CrossRef Search ADS PubMed Balasubramaniam KN , Beisner BA , Berman CM , De Marco A , Duboscq J , Koirala S , Majolo B , MacIntosh AJ , McFarland R , Molesti S , et al. 2017 . The influence of phylogeny, social style, and sociodemographic factors on macaque social network structure . Am J Primatol. 80 : e22727 . Barchas PR , Mendoza SD . 1984 . Emergent hierarchical relationships in rhesus macaques: An application of chase’s model . In: Barchas PR , editor. Social hierarchies: essays towards a sociophysiological perspective . Westport (CT) : Greenwood Press . p. 81 . Barrett L , Henzi SP , Weingrill T , Lycett JE , Hill RA . 1999 . Market forces predict grooming reciprocity in female baboons . Proc R Soc Lond B . 266 : 665 – 670 . Google Scholar CrossRef Search ADS Butovskaya ML . 2008 . Reconciliation, dominance and cortisol levels in children and adolescents (7–15-year-old boys) . Behaviour . 145 : 1557 – 1576 . Google Scholar CrossRef Search ADS Campennì M , Schino G . 2014 . Partner choice promotes cooperation: the two faces of testing with agent-based models . J Theor Biol . 344 : 49 – 55 . Google Scholar CrossRef Search ADS PubMed Clutton-Brock T . 2009 . Cooperation between non-kin in animal societies . Nature . 462 : 51 – 57 . Google Scholar CrossRef Search ADS PubMed Crockford C , Wittig RM , Langergraber K , Ziegler TE , Zuberbuehler K , Deschner T . 2013 . Urinary oxytocin and social bonding in related and unrelated wild chimpanzees . Proc Biol Sci . 280 : 20122765 . Google Scholar CrossRef Search ADS PubMed Crook JH , Ellis JE , Gosscustard JD . 1976 . Mammalian social-systems - structure and function . Anim Behav . 24 : 261 – 274 . Google Scholar CrossRef Search ADS Csardi G , Nepusz T . 2006 . The igraph software package for complex network research . InterJournal Complex Systems . 1695 . Dunbar RIM . 1991 . Functional significance of social grooming in primates . Folia Primatol . 57 : 121 – 131 . Google Scholar CrossRef Search ADS Evers E , de Vries H , Spruijt BM , Sterck EH . 2011 . Better safe than sorry–socio-spatial group structure emerges from individual variation in fleeing, avoidance or velocity in an agent-based model . PLoS One . 6 : e26189 . Google Scholar CrossRef Search ADS PubMed Evers E , de Vries H , Spruijt BM , Sterck EH . 2015 . Emotional bookkeeping and high partner selectivity are necessary for the emergence of partner-specific reciprocal affiliation in an agent-based model of primate groups . PLoS One . 10 : e0118921 . Google Scholar CrossRef Search ADS PubMed Franz M , McLean E , Tung J , Altmann J , Alberts SC . 2015 . Self-organizing dominance hierarchies in a wild primate population . Proc Biol Sci . 282 . Freeberg TM , Dunbar RI , Ord TJ . 2012 . Social complexity as a proximate and ultimate factor in communicative complexity . Philos Trans R Soc Lond B Biol Sci . 367 : 1785 – 1801 . Google Scholar CrossRef Search ADS PubMed Galef BG , Jr . 1988 . Imitation in animals: History, definitions, and interpretation of data from the psychological laboratory . In: Zentall T , Galef B , editors. Social learning: psychobiological and biological perspectives . Hillsdale (NJ) : Erlbaum . p. 3 . Graves FC , Wallen K , Maestripieri D . 2002 . Opioids and attachment in rhesus macaque (Macaca mulatta) abusive mothers . Behav Neurosci . 116 : 489 – 493 . Google Scholar CrossRef Search ADS PubMed Hamilton WD . 1964 . The genetical evolution of social behaviour. I . J Theor Biol . 7 : 1 – 16 . Google Scholar CrossRef Search ADS PubMed Heesen M , Macdonald S , Ostner J , Schülke O . 2015 . Ecological and social determinants of group cohesiveness and within-group spatial position in wild assamese macaques . Ethology . 121 : 270 – 283 . Google Scholar CrossRef Search ADS Heesen M , Rogahn S , Macdonald S , Ostner J , Schulke O . 2014 . Predictors of food-related aggression in wild assamese macaques and the role of conflict avoidance . Behav Ecol Sociobiol . 68 : 1829 – 1841 . Google Scholar CrossRef Search ADS Hemelrijk CK . 1990a . A matrix partial correlation test used in investigations of reciprocity and other social interaction patterns at group level . J Theor Biol . 143 : 405 – 420 . Google Scholar CrossRef Search ADS Hemelrijk CK . 1990b . Models of, and tests for, reciprocity, unidirectional and other social interaction patterns at a group level . Anim Behav . 39 : 1013 – 1029 . Google Scholar CrossRef Search ADS Hemelrijk CK . 1999 . An individual-oriented model on the emergence of despotic and egalitarian societies . P Roy Soc Lond B Bio . 266 : 361 – 369 . Google Scholar CrossRef Search ADS Hemelrijk CK . 2000 . Towards the integration of social dominance and spatial structure . Anim Behav . 59 : 1035 – 1048 . Google Scholar CrossRef Search ADS PubMed Hemelrijk CK , Kappeler PM , Puga-Gonzalez I . 2017 . The self-organization of social complexity in group-living animals: Lessons from the DomWorld model . In: Naguib M , Podos J , Simmons LW , et al. , editors. Advances in the study of behavior . p. 361 . Google Scholar CrossRef Search ADS Hemelrijk CK , Puga-Gonzalez I . 2012 . An individual-oriented model on the emergence of support in fights, its reciprocation and exchange . PLoS One . 7 : e37271 . Google Scholar CrossRef Search ADS PubMed Hemelrijk CK , Wantia J , Isler K . 2008a . Female dominance over males in primates: self-organisation and sexual dimorphism . PLoS One . 3 : e2678 . Google Scholar CrossRef Search ADS Hemelrijk CK , Wantia J , Isler K . 2008b . The more males, the more dominant are female primates . Folia Primatol . 79:337-337 . Hochberg Y . 1988 . A sharper bonferroni procedure for multiple tests of significance . Biometrika . 75:800-802 . Hogeweg P . 1988 . MIRROR beyond MIRROR, puddles of LIFE . In: Langton C , editor. Artificial life, SFI studies in the sciences of complexity . Redwood City (CA) : Adisson-Wesley Publishing Company . p. 297 Hsu Y , Wolf LL . 1999 . The winner and loser effect: integrating multiple experiences . Anim Behav . 57 : 903 – 910 . Google Scholar CrossRef Search ADS PubMed Jaeggi AV , De Groot E , Stevens JMG , Van Schaik CP . 2013 . Mechanisms of reciprocity in primates: testing for short-term contingency of grooming and food sharing in bonobos and chimpanzees . Evol Hum Behav . 34 : 69 – 77 . Google Scholar CrossRef Search ADS Kasper C , Voelkl B . 2009 . A social network analysis of primate groups . Primates . 50 : 343 – 356 . Google Scholar CrossRef Search ADS PubMed Keverne EB , Martensz ND , Tuite B . 1989 . Beta-endorphin concentrations in cerebrospinal fluid of monkeys are influenced by grooming relationships . Psychoneuroendocrinology . 14 : 155 – 161 . Google Scholar CrossRef Search ADS PubMed Krause J , Ruxton GD . 2002 . Living in groups . Oxford : Oxford University Press . Kubenova B , Konecna M , Majolo B , Smilauer P , Ostner J , Schülke O . 2017 . Triadic awareness predicts partner choice in male-infant-male interactions in Barbary macaques . Anim Cogn . 20 : 221 – 232 . Google Scholar CrossRef Search ADS PubMed Kudo H , Dunbar RIM . 2001 . Neocortex size and social network size in primates . Anim Behav . 62 : 711 – 722 . Google Scholar CrossRef Search ADS Lehmann J , Korstjens AH , Dunbar RIM . 2007 . Group size, grooming and social cohesion in primates . Anim Behav . 74 : 1617 – 1629 . Google Scholar CrossRef Search ADS Martel FL , Nevison CM , Simpson MJ , Keverne EB . 1995 . Effects of opioid receptor blockade on the social behavior of rhesus monkeys living in large family groups . Dev Psychobiol . 28 : 71 – 84 . Google Scholar CrossRef Search ADS PubMed McFarland R , Fuller A , Hetem RS , Mitchell D , Maloney SK , Henzi SP , Barrett L . 2015 . Social integration confers thermal benefits in a gregarious primate . J Anim Ecol . 84 : 871 – 878 . Google Scholar CrossRef Search ADS PubMed Molesti S , Majolo B . 2013 . Grooming increases self-directed behaviour in wild barbary macaques, macaca sylvanus . Anim Behav . 86 : 169 – 175 . Google Scholar CrossRef Search ADS Newman MEJ . 2006 . Modularity and community structure in networks . Proc Natl Acad Sci . 103 : 8577 – 8582 . Google Scholar CrossRef Search ADS PubMed Noë R , Hammerstein P . 1994 . Biological markets: supply and demand determine the effect of partner choice in cooperation, mutualism and mating . Behav Ecol Sociobiol . 35 : 1 – 11 . Google Scholar CrossRef Search ADS Pasquaretta C , Levé M , Claidière N , van de Waal E , Whiten A , MacIntosh AJ , Pelé M , Bergstrom ML , Borgeaud C , Brosnan SF , et al. 2014 . Social networks in primates: smart and tolerant species have more efficient networks . Sci Rep . 4 : 7600 . Google Scholar CrossRef Search ADS PubMed Paxton R , Basile BM , Adachi I , Suzuki WA , Wilson ME , Hampton RR . 2010 . Rhesus monkeys (Macaca mulatta) rapidly learn to select dominant individuals in videos of artificial social interactions between unfamiliar conspecifics . J Comp Psychol . 124 : 395 – 401 . Google Scholar CrossRef Search ADS PubMed Perry S , Barret CH , Manson JH . 2004 . White-faced capuchin monkeys show triadic awareness in their choices of allies . Animal Behav . 67 : 165 – 170 . Google Scholar CrossRef Search ADS Puga-Gonzalez I , Hildenbrandt H , Hemelrijk CK . 2009 . Emergent patterns of social affiliation in primates, a model . PLoS Comput Biol . 5 : e1000630 . Google Scholar CrossRef Search ADS PubMed Puga-Gonzalez I , Hoscheid A , Hemelrijk CK . 2015 . Friendships, reciprocation and interchange in an individual-based model . Behav Ecol Sociobiol . 69 : 383 – 394 . Google Scholar CrossRef Search ADS Puga-Gonzalez I , Ostner J , Schuelke O , Sosa S , Thierry B , Sueur C . 2018 . Data from: mechanisms of reciprocity and diversity in social networks: a modeling and comparative approach . Dryad Digital Repository . http://dx.doi.org/10.5061/dryad.6m220kd Puga-Gonzalez I , Sueur C . 2017a . Friendships and social networks in an individual-based model of primate social behaviour . J Artif Soc Soc Simul . 20 : 10 . Google Scholar CrossRef Search ADS Puga-Gonzalez I , Sueur C . 2017b . Emergence of complex social networks from spatial structure and rules of thumb: a modelling approach . Ecol Complex . 31 : 189 – 200 . Google Scholar CrossRef Search ADS Rankin RW , Mann J , Singh L , Patterson EM , Krzyszczyk E , Bejder L . 2016 . The role of weighted and topological network information to understand animal social networks: a null model approach . Anim Behav . 113 : 215 – 228 . Google Scholar CrossRef Search ADS Sabbatini G , De Bortoli Vizioli A , Visalberghi E , Schino G . 2012 . Food transfers in capuchin monkeys: an experiment on partner choice . Biol Lett . 8 : 757 – 759 . Google Scholar CrossRef Search ADS PubMed Schino G . 2007 . Grooming and agonistic support: a meta-analysis of primate reciprocal altruism . Behav Ecol . 18 : 115 – 120 . Google Scholar CrossRef Search ADS Schino G , Aureli F . 2008a . Grooming reciprocation among female primates: a meta-analysis . Biol Lett . 4 : 9 – 11 . Google Scholar CrossRef Search ADS Schino G , Aureli F . 2008b . Trade-offs in primate grooming reciprocation: testing behavioural flexibility and correlated evolution . Biol J Linn Soc . 95 : 439 – 446 . Google Scholar CrossRef Search ADS Schino G , Aureli F . 2009 . Reciprocal altruism in primates: partner choice, cognition, and emotions . Adv Stud Behav . 39 : 45 – 69 . Google Scholar CrossRef Search ADS Schino G , Scucchi S , Maestripieri D , Turillazzi PG . 1988 . Allogrooming as a tension-reduction mechanism: a behavioral approach . Am J Primatol . 16 : 43 – 50 . Google Scholar CrossRef Search ADS Schino G , Tiddi B , Di Sorrentino EP . 2006 . Simultaneous classification by rank and kinship in japanese macaques . Anim Behav . 71 : 1069 – 1074 . Google Scholar CrossRef Search ADS Schülke O , Bhagavatula J , Vigilant L , Ostner J . 2010 . Social bonds enhance reproductive success in male macaques . Curr Biol . 20 : 2207 – 2210 . Google Scholar CrossRef Search ADS PubMed Semple S , Harrison C , Lehmann J . 2013 . Grooming and anxiety in barbary macaques . Ethology . 119 : 779 – 785 . Google Scholar CrossRef Search ADS Setchell JM , Smith T , Wickings EJ , Knapp LA . 2008 . Social correlates of testosterone and ornamentation in male mandrills . Horm Behav . 54 : 365 – 372 . Google Scholar CrossRef Search ADS PubMed Shutt K , MacLarnon A , Heistermann M , Semple S . 2007 . Grooming in Barbary macaques: better to give than to receive ? Biol Lett . 3 : 231 – 233 . Google Scholar CrossRef Search ADS PubMed Silk JB . 1999 . Male bonnet macaques use information about third-party rank relationships to recruit allies . Anim Behav . 58 : 45 – 51 . Google Scholar CrossRef Search ADS PubMed Silk JB . 2007 . Social components of fitness in primate groups . Science . 317 : 1347 – 1351 . Google Scholar CrossRef Search ADS PubMed Silk JB , Alberts SC , Altmann J . 2003 . Social bonds of female baboons enhance infant survival . Science . 302 : 1231 – 1234 . Google Scholar CrossRef Search ADS PubMed Stevens JR , Hauser MD . 2004 . Why be nice? Psychological constraints on the evolution of cooperation . Trends Cogn Sci . 8 : 60 – 65 . Google Scholar CrossRef Search ADS PubMed Stevens JR , Volstorf J , Schooler LJ , Rieskamp J . 2011 . Forgetting constrains the emergence of cooperative decision strategies . Front Psychol . 1 : 235 . Google Scholar CrossRef Search ADS PubMed Sueur C , Deneubourg JL , Petit O , Couzin ID . 2011 . Group size, grooming and fission in primates: a modeling approach based on group structure . J Theor Biol . 273 : 156 – 166 . Google Scholar CrossRef Search ADS PubMed Sueur C , Petit O , De Marco A , Jacobs AT , Watanabe K , Thierry B . 2011 . A comparative network analysis of social style in macaques . Anim Behav . 82 : 845 – 852 . Google Scholar CrossRef Search ADS Team-R-Core . 2015 . R: a language and environment for statistical computing . Vienna (Austria) : R Foundation for Statistical Computing . Thierry B . 2004 . Social epigenesis . In: Thierry B , Singh M , Kaumanns W , editors. Macaque societies: A model for the study of social organisation . Cambridge : Cambridge University Press . p. 267 . Tiddi B , Aureli F , di Sorrentino EP , Janson CH , Schino G . 2011 . Grooming for tolerance? two mechanisms of exchange in wild tufted capuchin monkeys . Behav Ecol . 22 : 663 – 669 . Google Scholar CrossRef Search ADS Tomasello M , Call J . 1997 . Primate cognition . New York : Oxford University Press . Trivers RL . 1971 . The evolution of reciprocal altruism . Q Rev Biol . 46 : 35 – 57 . Google Scholar CrossRef Search ADS de Waal FBM . 2000 . Attitudinal reciprocity in food sharing among brown capuchin monkeys . Animal Behav . 60 : 253 – 261 . Google Scholar CrossRef Search ADS de Waal FBM , Brosnan SF . 2006 . Simple and complex reciprocity in primates . In: Kappeler PM , van Schaick CP , editors. Cooperation in primates and humans: Mechanisms and evolution . Berlin (Germany) : Springer . p. 85 . Google Scholar CrossRef Search ADS de Waal FBM , Luttrell LM . 1988 . Mechanisms of social reciprocity in three primate species: symmetrical relationship characteristics or cognition ? Ethol Sociobiol . 9 : 101 – 118 . Google Scholar CrossRef Search ADS Ward A , Webster M . 2016 . Sociality: the behaviour of group-living animals . 1st ed . Switzerland: Springer International Publishing . Google Scholar CrossRef Search ADS Wittig RM , Crockford C , Weltring A , Langergraber KE , Deschner T , Zuberbühler K . 2016 . Social support reduces stress hormone levels in wild chimpanzees across stressful events and everyday affiliations . Nat Commun . 7 : 13361 . Google Scholar CrossRef Search ADS PubMed © The Author(s) 2018. Published by Oxford University Press on behalf of the International Society for Behavioral Ecology. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

Journal

Behavioral EcologyOxford University Press

Published: Apr 9, 2018

There are no references for this article.

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