Sociality
From Evolution and Games
(→Editorial Information) |
(→Run the simulations) |
||
| (23 intermediate revisions not shown) | |||
| Line 1: | Line 1: | ||
[[Category:It takes grouping and cooperation to get sociality]] | [[Category:It takes grouping and cooperation to get sociality]] | ||
| - | + | The software we provide and describe here reimplements simulations first described in [http://www.pnas.org/content/99/22/14268.abstract Aviles 2002] using Java . Our simulation program also has alternative fitness functions and some additional features. If you have questions or require further information about this work or this tutorial please do email [mailto:garciaATevolbio.mpg.de Julian Garcia], [mailto:c.m.veelenATuva.nl Matthijs van Veelen] or [mailto:lavilesATzoology.ubc.ca Leticia Avilés]. | |
| - | + | = Editorial Information = | |
| + | [http://www1.fee.uva.nl/creed/people/veelen/index.shtml van Veelen, Matthijs], [http://staff.feweb.vu.nl/sociality/ García Julián] and [http://www.zoology.ubc.ca/person/laviles Avilés Leticia]. '''It takes grouping and cooperation to get sociality'''. ''Journal of Theoretical Biology'' 264(4) Pages 1240-1253. [http://dx.doi.org/10.1016/j.jtbi.2010.02.043 Download the article] | ||
| - | = | + | = Abstract = |
| - | + | Cooperation and grouping are regularly studied as separate traits. The evolution of sociality however requires both that individuals get together in groups and that they cooperate within them. Because the level of cooperation can influence selection for group size, and vice versa, it is worth studying how these traits coevolve. Using a generally applicable two-trait optimization approach, we provide analytical solutions for three specific models. These solutions describe how cooperative associations of nonrelatives evolve, and predict how large and how cooperative they will be. The analytical solutions help understand how changes in parameter values, such as the group carrying capacity and the costs of cooperation, affect group size and the level of cooperation in equilibrium. Although the analytical model makes a few simplifying assumptions - populations are assumed to be monomorphic for grouping as well as for cooperative tendencies, and group size is assumed to be deterministic - simulations show that its predictions are matched quite closely by results for settings where these assumptions do not hold. | |
| - | + | = Run the simulations = | |
| - | + | Click on the link corresponding to the model you want to run, and ''choose download and double click, or just open with Java Runtime''. You should set the parameters and press the ''play button'' to start a run. If you want to change the parameters press stop, and then refresh so that new values can be introduced. Feel free to re-arrange the windows or zoom into the graphs. | |
| - | + | ||
| - | + | * [http://www.evolbio.mpg.de/~garcia/jars/LimitedOnlineVersion.jar Start the Limited Nesting Sites model] | |
| - | + | * [http://www.evolbio.mpg.de/~garcia/jars/GlobalOnlineVersion.jar Start the Global Density Dependence model] | |
| - | Individuals are characterized by a diploid genome coding for two separate traits; grouping and cooperative tendencies. The number of loci for each of these traits can be chosen by the user of the program. The number of loci for grouping will be denoted by n, and the number of loci for cooperation will be denoted by m. We followed [http://www.pnas.org/content/99/22/14268.abstract Avilés et al. (2002)] by choosing 15 as a default value for the number of loci for both traits. The loci of an individual will be numbered from < | + | * [http://www.evolbio.mpg.de/~garcia/jars/ConstantOnlineVersion.jar Start the Constant Population Size model] |
| + | |||
| + | = The genetic system = | ||
| + | |||
| + | == What individuals are made of == | ||
| + | |||
| + | Individuals are characterized by a diploid genome coding for two separate traits; grouping and cooperative tendencies. The number of loci for each of these traits can be chosen by the user of the program. The number of loci for grouping will be denoted by n, and the number of loci for cooperation will be denoted by m. We followed [http://www.pnas.org/content/99/22/14268.abstract Avilés et al. (2002)] by choosing 15 as a default value for the number of loci for both traits. The loci of an individual will be numbered from <math>1</math> to <math>n+m</math>, where the first n loci code for grouping tendencies and the loci from n+1 to n+m code for cooperative tendencies. The alleles will be indicated as left and right. At each locus, two alleles (1 or 0) can occur. The tendency of individual i to form groups is denoted by gi, and equals the proportion of 1’s in the grouping portion of its genome, or the number of 1’s divided by 2n. The tendency of individual i to cooperate, denoted by Íi, is the proportion of 1’s in the cooperation part of its genome, or the number of 1’s in the second part of its genome, divided by 2m. This implies that gi and Íi vary between 0 and 1, and are multiples of 1 over twice the number of loci. | ||
[[File:GenomeSociality.png|500px|center|thumb|Example genome as represented in the computer]] | [[File:GenomeSociality.png|500px|center|thumb|Example genome as represented in the computer]] | ||
| Line 23: | Line 29: | ||
<center> | <center> | ||
{| | {| | ||
| - | |[[File: | + | | [[File:SlowCooperation.png|400px|center|thumb|Cooperation is slow relative to grouping tendencies]] ||[[File:FastCooperation.png|390px|center|thumb|Cooperation is fast relative to grouping tendencies]]|| |
|} | |} | ||
</center> | </center> | ||
| + | |||
| + | == Mutations; inversions and translocations first == | ||
| + | |||
| + | |||
| + | At the mutation stage (see below for when that is) one Bernouilli trial determines whether or not a mutation occurs in the grouping part of the genome of an individual, and another Bernouilli trial determines whether or not a mutation occurs at the cooperation part. These trials are independently run for every population member again. The probability <math>p</math> with which a mutation occurs is chosen by the user of the program, and <math>1-p</math> obviously is the probability with which no mutation occurs. If the outcome of the Bernouilli trial is, for instance, that a mutation will occur in the grouping part of the genome, then first a random allele is selected from that part of the genome, with all alleles having equal probability to be chosen. Then the program tries to find an allele that is dissimilar to it (that is: a <math>0</math> if it is a <math>1</math> itself, and a <math>1</math> if it is a <math>0</math> itself) in order to swap places. It starts looking for a dissimilar allele at the same side of the genome. Suppose the randomly determined first allele is the <math>left</math> allele at locus <math>j \in \{1,...,n\}</math> The program then starts checking the <math>left</math> allele at locus <math>j+ 1</math>. If that is the same as the <math>left</math> allele at locus <math>j</math>, then it moves on to the <math>left</math> allele at locus <math>j+2</math>, and so on, and after locus <math>n</math> it continues looking from locus <math>1</math> and goes on until locus <math>j-1</math>. If it has found a first different allele, then the two are interchanged, which is equivalent to inverting the string in between and including them. If it has not found a different allele, then it continues looking for a different allele at the other side; in this case, it starts at the right allele at locus <math>j</math>, and then repeats the search for an allele that is different. If the program finds one, the <math>0</math> and the <math>1</math> are interchanged, imitating translocation. If the program does not find one here either, then there are either only <math>1</math>’s or only <math>0</math>’s to be found, and then a random allele is chosen to be changed into the other value. | ||
| + | |||
| + | |||
| + | This procedure is repeated for all individuals in the pool. Only at the edges, where a substantial share of the individuals will be either all be all <math>1</math>’s or all be all <math>0</math>’s, there will be mutation pressure with a bias against the most common allele, but the priority of inversion and translocation reduces that as much as possible when not almost at the extremes. | ||
| + | |||
| + | |||
| + | == Mating == | ||
| + | |||
| + | If an individual reproduces, then a random other member of the population is chosen with which it will reproduce. This implies that there is sex, but there are no different sexes. Both individuals then produce a gamete, and these together form the new individual. | ||
| + | |||
| + | == Producing gametes == | ||
| + | |||
| + | If an individual produces a gamete, then locus 1 of the gamete will with probability one half be a copy of the <math>right</math> allele at locus one of the individual, and with probability one half it will be a copy of the <math>left</math> allele at locus one of the individual. The next locus of the gamete will be a copy from the same side of the genome with probability <math>0.5 + 0.5 P_l</math>, and with probability <math>0.5 - 0.5 P_l</math> it will be a copy from the other side of the genome. Here, <math>P_l \in [0,1]</math> is the linkage level. In the current version of the program, the linkage is set to <math>0</math>. | ||
| + | |||
| + | == Asexual reproduction == | ||
| + | |||
| + | In the program, the user can also uncheck sexual reproduction. In that case individuals just produce exact copies of themselves. Automatically, the mutation procedure will then also change; it will not consist of inversion and translocation, but only change <math>0</math>’s in <math>1</math>’s and <math>1</math>’s in <math>0</math>’s. If asexual reproduction were to be combined with translocation and inversion first, then once all individuals share the same number of loci for grouping and/or for cooperation, this would remain stable forever, if that number is positive, because translocation and inversion keep the number of <math>1</math>’s and <math>0</math>’s constant. | ||
| + | |||
| + | == Initial population == | ||
| + | |||
| + | The default initial population consists of individuals with grouping and cooperative tendency <math>0</math>. The user can choose to initialize the population randomly with positive grouping and cooperative tendencies by choosing positive values for parameters <math>ProbabilityActiveLociCooperation</math> and <math>ProbabilityActiveLociGroupishness</math>. Then the initial population is determined by drawing every allele in the genome of every individual randomly, with these probabilities as parameters. | ||
| + | |||
| + | = Lifecycle = | ||
| + | |||
| + | == Global pool and group formation == | ||
| + | |||
| + | A central and shared feature of all simulation models encompassed in the program is that there are two distinct stages in a life cycle. At one stage, all individuals are gathered together without structure in a global population pool. At the other stage, individuals form groups, in which reproduction takes place. The transition of a population from the global pool to the groups, and the reproduction within the groups before the offspring returns to the global pool is where the simulation program offers different possibilities. Mutation takes place when individuals are in the global pool. | ||
| + | |||
| + | <center> | ||
| + | [[File:LifecycleSociality.png|600px|thumb|center|Lifecycle]] | ||
| + | </center> | ||
| + | |||
| + | == From the global pool to the groups == | ||
| + | Individuals leave the global pool in random order. When they leave, they go to visit sites, where they either initiate new groups, or join already formed groups. If a site is empty, they ''join'' (initiate a new group) with probability <math>1</math>. If a site is not empty, then the program offers two possibilities for the probability with which they join, from which the user can choose. | ||
| + | |||
| + | For the default option, the probability that individual <math>i</math>, that itself has grouping tendency <math>g_i</math>, will join a group whose average grouping tendency of <math>\bar{g}</math> is | ||
| + | <center> | ||
| + | <math>g_i \times \bar{g}</math> | ||
| + | </center> | ||
| + | The second option for the joining probability is called ''truncated joining'' in the computer program, and it truncates the possible group sizes at <math>n^{*}</math>, and reduces the probability of joining more for larger groups than for smaller groups. Here the probability of joining a group that consists of n individuals is | ||
| + | |||
| + | <center> | ||
| + | <math>g_i \times \bar{g} \left(1- \frac{n}{n^{*}}\right)</math> | ||
| + | </center> | ||
| + | |||
| + | for <math>n < n^{*}</math>, where <math>n^{*}</math> is the expected stable group size. The stable group size is the fixed point of the fitness function (see below), assuming that <math>\gamma_i = \bar{\gamma}</math>. Whenever <math>n\geq n^{*}</math>, the group acceptance probability is 0. | ||
| + | |||
| + | |||
| + | The second option matches the model in [http://www.journals.uchicago.edu/doi/abs/10.1086/324792 Avilés et al. (2002)]. Given a distribution of <math>g</math>’s in the population, the distribution of sizes will be somewhat different between the two. Also, if selection would produce the same average group size, the shape of the distribution of group sizes will differ a bit. Running the two options however shows that this does not really make a difference for the group size in the fixed point. Because the computation of <math>n^{*}</math> requires information that individuals in realistic situations may not have, one can therefore, without changing the model results, resort to the first version that requires less information. The only advantage of the second, truncated version is that it gets to the fixed point faster. | ||
| + | |||
| + | The program allows for three options for the order in which the sites are visited. One option (called ''Random'' in the program) is that every individual visits the sites in random order. Another option (''Forward'') is that everyone starts at site number ''1'', then goes to site number ''2'', then to number ''3'' and so on. The third option (''Backwards'') is that everyone starts trying to be admitted at the non-empty site with the highest number. Say that this is site number ''n''. Then it goes backwards, first to site ''n-1'', then to ''n-2'' and so on, until it is at site number ''1''. If the individual is then still not admitted, it starts a new group at site ''n+1''. | ||
| + | |||
| + | These different orders of visiting will produce different size distributions; with ''Forward'', the variance in sizes will be the largest, with ''Backwards'' it will be the smallest, and with ''Random'' it will be somewhere in between. | ||
| + | |||
| + | How many sites there will be and when the flow from the global population pool to the groups stops, depends on the way in which the size of the global population is balanced; there are three options (limited nesting sites, global density dependence, limited total population) that are discussed below. | ||
| + | |||
| + | == From the groups to the global pool == | ||
| + | Within the group, a hump-shaped function will determine the expected number of offspring that an individual contributes to the global pool. In the program, we coded three such functions, called <math>f_i</math>, <math>g_i</math> and <math>h_i</math> respectively. The expected number of offspring of individual <math>i</math> with cooperativeness <math>\gamma_i</math>, in a group with <math>n</math> individuals in it and average cooperativeness <math>\overline{\gamma}</math> is then | ||
| + | |||
| + | |||
| + | <center> | ||
| + | <math>f_{i}\left(n,\overline{\gamma },\gamma _{i}\right) =e^{r}e^{-cn}n^{\overline{\gamma}}\left(1+\beta\left(\overline{\gamma}-\gamma_{i}\right)\right)</math> | ||
| + | |||
| + | |||
| + | <math>g_{i}\left(n,\overline{\gamma },\gamma _{i}\right)=e^{r}e^{-cn}\left( 1+n\overline{\gamma}\right) \left( 1+\beta \left(\overline{\gamma }-\gamma_{i}\right) \right)</math> | ||
| + | |||
| + | |||
| + | <math>h_{i}\left( n,\overline{\gamma },\gamma _{i}\right)=e^{r}e^{-cn}n\overline{\gamma }\left( 1+\beta \left( \overline{\gamma }-\gamma _{i}\right)\right)</math> | ||
| + | </center> | ||
| + | |||
| + | where <math>r</math> is an intrinsic rate of growth parameter, <math>c</math> is a local negative density-dependent parameter inversely proportional to a group carrying capacity, and <math>\beta</math> determines the cost of cooperation. The functions <math>g_i</math> and <math>h_i</math> can be seen as special cases of a function | ||
| + | |||
| + | <center> | ||
| + | <math>e^{r}e^{-cn}\left( \alpha +n\overline{\gamma }\right) \left(1+\beta \left(\overline{\gamma }-\gamma _{i}\right) \right) </math>, | ||
| + | </center> | ||
| + | with <math>\alpha</math> determining the relative importance of cooperation. | ||
| + | |||
| + | The default initialization of the population is one where all individuals have grouping tendency ''0'' and cooperative tendency ''0''. For <math>f_i</math> and <math>h_i</math>, not too low choices for <math>r</math> can guarantee that the population will not go extinct from that initial population. For <math>g_i</math>, that will be different; when <math>\gamma = 0</math>, <math>g_i(n,\gamma,\gamma_i)</math> will be <math>0</math> too, and the population will go extinct. One can however initialize the population with higher values of grouping and cooperative tendencies, by choosing positive probabilities with which alleles in the initial population will be a <math>1</math> instead of a <math>0</math>. | ||
| + | |||
| + | == Global population control: Limited nesting sites, global density dependence and limited total population == | ||
| + | |||
| + | The ecology can put a limit on the growth of the overall population in different ways. One can assume that there is only a limited number of nesting sites. In Avilés et al. ([http://www.journals.uchicago.edu/doi/abs/10.1086/324792 2002], [http://www.journals.uchicago.edu/doi/abs/10.1086/422263 2004]) this is implemented by choosing a fixed number of sites, and assuming that the flow from the global pool towards the groups stops as soon as the last nesting site has been occupied by one individual. The Limited Nesting Sites version in this program also does exactly that. If the user however unchecks the Limited Season option, then the flow from the population pool does not stop when the last nesting site is occupied, but goes on until the global pool is empty. Because individuals that do not get admitted to any of the groups are discarded, this implies that selective pressure for higher grouping tendencies is increased. The possibility of not being admitted at all therefore heavily distorts the link between grouping tendency and how that brings the expected group size closer to the optimal group size; not getting into a group at all will always be (much) worse than making it into a group that has a suboptimal size. Grouping tendencies and group size will therefore explode. | ||
| + | |||
| + | One can also assume that the global density negatively affects individual fitnesses. In the Global Density Dependence version of the simulation model, as in Avilés et al. (2002, 2004), the local per capita growth function (in our case: <math>f_i</math>,<math>g_i</math> or <math>h_i</math>) is multiplied by a factor | ||
| + | |||
| + | <center> | ||
| + | <math> | ||
| + | e^{-w\left( N-n\right) } | ||
| + | </math> | ||
| + | </center> | ||
| + | |||
| + | where <math>w</math> is a global negative density-dependence parameter, inversely proportional to a global carrying capacity, and <math>N</math> is the total number of individuals in the global population. | ||
| + | |||
| + | A third option is to assume that every generation, the global pool is reduced to a fixed number of individuals before they leave in order to form groups. This option (Fixed Population Size) is similar to the Limited Nesting Sites; the difference is that with the Limited Nesting Sites, the number of groups is fixed, and the total population size is a random variable, with a very moderate variance, while with a Fixed Population Size, the total population size is fixed, and the number of groups is a random variable, with a very moderate variance. The average group size will be a random variable for both cases, and will be pretty comparable. | ||
| + | |||
| + | The most important distinctive feature of the global density dependence version is that, in combination with particular choices of the parameters, this model can allow for cyclical or chaotic behaviour of the size of the total population. Group sizes will cycle along, and this implies that the optimal grouping and cooperative tendencies become moving targets; if the overall size of the population is large, then low grouping and cooperative tendencies can do relatively well, but the next period the overall population size might already be very small, implying that high grouping and cooperative tendencies do better. | ||
| + | |||
| + | === Equlibrium values for different models === | ||
| + | The figures below compare equilibrium values for a Limited Nesting Sites simulation (''first figure'') and a Global Density Dependence simulation (''second figure''). Equilibrium values are shifted to the right in the Global density dependence as a result of moving optimal targets. For further details see the paper. | ||
| + | |||
| + | [[File:LimitedNestingEquilibriaSociality.png|center|thumb|450px|Limited nesting sites. Function G.]] | ||
| + | [[File:GlobalDensityEquilibriaSociality.png|center|thumb|450px|Global density dependence. Function G. R equals 1.]] | ||
Current revision as of 09:25, 21 October 2011
The software we provide and describe here reimplements simulations first described in Aviles 2002 using Java . Our simulation program also has alternative fitness functions and some additional features. If you have questions or require further information about this work or this tutorial please do email Julian Garcia, Matthijs van Veelen or Leticia Avilés.
Contents |
Editorial Information
van Veelen, Matthijs, García Julián and Avilés Leticia. It takes grouping and cooperation to get sociality. Journal of Theoretical Biology 264(4) Pages 1240-1253. Download the article
Abstract
Cooperation and grouping are regularly studied as separate traits. The evolution of sociality however requires both that individuals get together in groups and that they cooperate within them. Because the level of cooperation can influence selection for group size, and vice versa, it is worth studying how these traits coevolve. Using a generally applicable two-trait optimization approach, we provide analytical solutions for three specific models. These solutions describe how cooperative associations of nonrelatives evolve, and predict how large and how cooperative they will be. The analytical solutions help understand how changes in parameter values, such as the group carrying capacity and the costs of cooperation, affect group size and the level of cooperation in equilibrium. Although the analytical model makes a few simplifying assumptions - populations are assumed to be monomorphic for grouping as well as for cooperative tendencies, and group size is assumed to be deterministic - simulations show that its predictions are matched quite closely by results for settings where these assumptions do not hold.
Run the simulations
Click on the link corresponding to the model you want to run, and choose download and double click, or just open with Java Runtime. You should set the parameters and press the play button to start a run. If you want to change the parameters press stop, and then refresh so that new values can be introduced. Feel free to re-arrange the windows or zoom into the graphs.
The genetic system
What individuals are made of
Individuals are characterized by a diploid genome coding for two separate traits; grouping and cooperative tendencies. The number of loci for each of these traits can be chosen by the user of the program. The number of loci for grouping will be denoted by n, and the number of loci for cooperation will be denoted by m. We followed Avilés et al. (2002) by choosing 15 as a default value for the number of loci for both traits. The loci of an individual will be numbered from 1 to n + m, where the first n loci code for grouping tendencies and the loci from n+1 to n+m code for cooperative tendencies. The alleles will be indicated as left and right. At each locus, two alleles (1 or 0) can occur. The tendency of individual i to form groups is denoted by gi, and equals the proportion of 1’s in the grouping portion of its genome, or the number of 1’s divided by 2n. The tendency of individual i to cooperate, denoted by Íi, is the proportion of 1’s in the cooperation part of its genome, or the number of 1’s in the second part of its genome, divided by 2m. This implies that gi and Íi vary between 0 and 1, and are multiples of 1 over twice the number of loci.
The choice for the genome size has an effect on the speed of selection. Reducing the genome size for any one of the tendencies increases the speed of selection for that tendency. This is because the marginal effect of one allele changing from a 0 into a 1 on fitness is larger with a smaller number of loci. The speed of selection of one of the tendencies relative to the other can obviously be tuned by adjusting the one genome size and keeping the other genome size constant.
 |  |
Mutations; inversions and translocations first
At the mutation stage (see below for when that is) one Bernouilli trial determines whether or not a mutation occurs in the grouping part of the genome of an individual, and another Bernouilli trial determines whether or not a mutation occurs at the cooperation part. These trials are independently run for every population member again. The probability p with which a mutation occurs is chosen by the user of the program, and 1 − p obviously is the probability with which no mutation occurs. If the outcome of the Bernouilli trial is, for instance, that a mutation will occur in the grouping part of the genome, then first a random allele is selected from that part of the genome, with all alleles having equal probability to be chosen. Then the program tries to find an allele that is dissimilar to it (that is: a 0 if it is a 1 itself, and a 1 if it is a 0 itself) in order to swap places. It starts looking for a dissimilar allele at the same side of the genome. Suppose the randomly determined first allele is the left allele at locus 
The program then starts checking the left allele at locus j + 1. If that is the same as the left allele at locus j, then it moves on to the left allele at locus j + 2, and so on, and after locus n it continues looking from locus 1 and goes on until locus j − 1. If it has found a first different allele, then the two are interchanged, which is equivalent to inverting the string in between and including them. If it has not found a different allele, then it continues looking for a different allele at the other side; in this case, it starts at the right allele at locus j, and then repeats the search for an allele that is different. If the program finds one, the 0 and the 1 are interchanged, imitating translocation. If the program does not find one here either, then there are either only 1’s or only 0’s to be found, and then a random allele is chosen to be changed into the other value.
This procedure is repeated for all individuals in the pool. Only at the edges, where a substantial share of the individuals will be either all be all 1’s or all be all 0’s, there will be mutation pressure with a bias against the most common allele, but the priority of inversion and translocation reduces that as much as possible when not almost at the extremes.
Mating
If an individual reproduces, then a random other member of the population is chosen with which it will reproduce. This implies that there is sex, but there are no different sexes. Both individuals then produce a gamete, and these together form the new individual.
Producing gametes
If an individual produces a gamete, then locus 1 of the gamete will with probability one half be a copy of the right allele at locus one of the individual, and with probability one half it will be a copy of the left allele at locus one of the individual. The next locus of the gamete will be a copy from the same side of the genome with probability 0.5 + 0.5Pl, and with probability 0.5 − 0.5Pl it will be a copy from the other side of the genome. Here, ![P_l \in [0,1]](/wiki/images/math/2/7/0/2704206d7fbaff15e16f559add667c69.png)
is the linkage level. In the current version of the program, the linkage is set to 0.
Asexual reproduction
In the program, the user can also uncheck sexual reproduction. In that case individuals just produce exact copies of themselves. Automatically, the mutation procedure will then also change; it will not consist of inversion and translocation, but only change 0’s in 1’s and 1’s in 0’s. If asexual reproduction were to be combined with translocation and inversion first, then once all individuals share the same number of loci for grouping and/or for cooperation, this would remain stable forever, if that number is positive, because translocation and inversion keep the number of 1’s and 0’s constant.
Initial population
The default initial population consists of individuals with grouping and cooperative tendency 0. The user can choose to initialize the population randomly with positive grouping and cooperative tendencies by choosing positive values for parameters ProbabilityActiveLociCooperation and ProbabilityActiveLociGroupishness. Then the initial population is determined by drawing every allele in the genome of every individual randomly, with these probabilities as parameters.
Lifecycle
Global pool and group formation
A central and shared feature of all simulation models encompassed in the program is that there are two distinct stages in a life cycle. At one stage, all individuals are gathered together without structure in a global population pool. At the other stage, individuals form groups, in which reproduction takes place. The transition of a population from the global pool to the groups, and the reproduction within the groups before the offspring returns to the global pool is where the simulation program offers different possibilities. Mutation takes place when individuals are in the global pool.
From the global pool to the groups
Individuals leave the global pool in random order. When they leave, they go to visit sites, where they either initiate new groups, or join already formed groups. If a site is empty, they join (initiate a new group) with probability 1. If a site is not empty, then the program offers two possibilities for the probability with which they join, from which the user can choose.
For the default option, the probability that individual i, that itself has grouping tendency gi, will join a group whose average grouping tendency of 
is
The second option for the joining probability is called truncated joining in the computer program, and it truncates the possible group sizes at n * , and reduces the probability of joining more for larger groups than for smaller groups. Here the probability of joining a group that consists of n individuals is
for n < n * , where n * is the expected stable group size. The stable group size is the fixed point of the fitness function (see below), assuming that 
. Whenever 
, the group acceptance probability is 0.
The second option matches the model in Avilés et al. (2002). Given a distribution of g’s in the population, the distribution of sizes will be somewhat different between the two. Also, if selection would produce the same average group size, the shape of the distribution of group sizes will differ a bit. Running the two options however shows that this does not really make a difference for the group size in the fixed point. Because the computation of n * requires information that individuals in realistic situations may not have, one can therefore, without changing the model results, resort to the first version that requires less information. The only advantage of the second, truncated version is that it gets to the fixed point faster.
The program allows for three options for the order in which the sites are visited. One option (called Random in the program) is that every individual visits the sites in random order. Another option (Forward) is that everyone starts at site number 1, then goes to site number 2, then to number 3 and so on. The third option (Backwards) is that everyone starts trying to be admitted at the non-empty site with the highest number. Say that this is site number n. Then it goes backwards, first to site n-1, then to n-2 and so on, until it is at site number 1. If the individual is then still not admitted, it starts a new group at site n+1.
These different orders of visiting will produce different size distributions; with Forward, the variance in sizes will be the largest, with Backwards it will be the smallest, and with Random it will be somewhere in between.
How many sites there will be and when the flow from the global population pool to the groups stops, depends on the way in which the size of the global population is balanced; there are three options (limited nesting sites, global density dependence, limited total population) that are discussed below.
From the groups to the global pool
Within the group, a hump-shaped function will determine the expected number of offspring that an individual contributes to the global pool. In the program, we coded three such functions, called fi, gi and hi respectively. The expected number of offspring of individual i with cooperativeness γi, in a group with n individuals in it and average cooperativeness 
is then
where r is an intrinsic rate of growth parameter, c is a local negative density-dependent parameter inversely proportional to a group carrying capacity, and β determines the cost of cooperation. The functions gi and hi can be seen as special cases of a function
,
with α determining the relative importance of cooperation.
The default initialization of the population is one where all individuals have grouping tendency 0 and cooperative tendency 0. For fi and hi, not too low choices for r can guarantee that the population will not go extinct from that initial population. For gi, that will be different; when γ = 0, gi(n,γ,γi) will be 0 too, and the population will go extinct. One can however initialize the population with higher values of grouping and cooperative tendencies, by choosing positive probabilities with which alleles in the initial population will be a 1 instead of a 0.
Global population control: Limited nesting sites, global density dependence and limited total population
The ecology can put a limit on the growth of the overall population in different ways. One can assume that there is only a limited number of nesting sites. In Avilés et al. (2002, 2004) this is implemented by choosing a fixed number of sites, and assuming that the flow from the global pool towards the groups stops as soon as the last nesting site has been occupied by one individual. The Limited Nesting Sites version in this program also does exactly that. If the user however unchecks the Limited Season option, then the flow from the population pool does not stop when the last nesting site is occupied, but goes on until the global pool is empty. Because individuals that do not get admitted to any of the groups are discarded, this implies that selective pressure for higher grouping tendencies is increased. The possibility of not being admitted at all therefore heavily distorts the link between grouping tendency and how that brings the expected group size closer to the optimal group size; not getting into a group at all will always be (much) worse than making it into a group that has a suboptimal size. Grouping tendencies and group size will therefore explode.
One can also assume that the global density negatively affects individual fitnesses. In the Global Density Dependence version of the simulation model, as in Avilés et al. (2002, 2004), the local per capita growth function (in our case: fi,gi or hi) is multiplied by a factor
where w is a global negative density-dependence parameter, inversely proportional to a global carrying capacity, and N is the total number of individuals in the global population.
A third option is to assume that every generation, the global pool is reduced to a fixed number of individuals before they leave in order to form groups. This option (Fixed Population Size) is similar to the Limited Nesting Sites; the difference is that with the Limited Nesting Sites, the number of groups is fixed, and the total population size is a random variable, with a very moderate variance, while with a Fixed Population Size, the total population size is fixed, and the number of groups is a random variable, with a very moderate variance. The average group size will be a random variable for both cases, and will be pretty comparable.
The most important distinctive feature of the global density dependence version is that, in combination with particular choices of the parameters, this model can allow for cyclical or chaotic behaviour of the size of the total population. Group sizes will cycle along, and this implies that the optimal grouping and cooperative tendencies become moving targets; if the overall size of the population is large, then low grouping and cooperative tendencies can do relatively well, but the next period the overall population size might already be very small, implying that high grouping and cooperative tendencies do better.
Equlibrium values for different models
The figures below compare equilibrium values for a Limited Nesting Sites simulation (first figure) and a Global Density Dependence simulation (second figure). Equilibrium values are shifted to the right in the Global density dependence as a result of moving optimal targets. For further details see the paper.
