Fission-fusion dynamics and group-size dependent composition in heterogeneous populations

Fission-fusion dynamics and group-size dependent composition in heterogeneous populations

Gokul G. Nair Department of Physics, Indian Institute of Science, Bangalore, 560 012, India    Athmanathan Senthilnathan Department of Mathematics, Indian Institute of Science, Bangalore, 560 012, India    Srikanth K. Iyer Department of Mathematics, Indian Institute of Science, Bangalore, 560 012, India    Vishwesha Guttal Centre for Ecological Sciences, Indian Institute of Science, Bangalore, 560 012, India
July 15, 2019

Many animal groups are heterogeneous and may even consist of individuals of different species, called mixed-species flocks. Mathematical and computational models of collective animal movement behaviour, however, typically assume that groups and populations consist of identical individuals. In this paper, using the mathematical framework of the coagulation-fragmentation process, we develop and analyse a model of merge and split group dynamics, also called fission-fusion dynamics, for heterogeneous populations that contain two types (or species) of individuals. We assume that more heterogeneous groups experience higher split rates than homogeneous groups, forming two daughter groups whose compositions are drawn uniformly from all possible partitions. We analytically derive a master equation for group size and compositions and find mean-field steady-state solutions. We predict that there is a critical group size below which groups are more likely to be homogeneous and contain the abundant type/species. Despite the propensity of heterogeneous groups to split at higher rates, we find that groups are more likely to be heterogeneous but only above the critical group size. Monte-Carlo simulation of the model show excellent agreement with these analytical model results. Thus, our model makes a testable prediction that composition of flocks are group-size dependent and do not merely reflect the population level heterogeneity. We discuss the implications of our results to empirical studies on flocking systems.

I Introduction

Collective phenomena and self-organisation are widespread in the animal kingdom Parrish and Edelstein-Keshet (1999); Camazine (2003); Sumpter (2010); Couzin and Krause (2003). Theory as well as empirical works suggest that these macroscopic behaviours often emerge from simple microscopic interactions among individuals Vicsek et al. (1995). Much of collective behaviour theory and models assume that individuals in populations are identical Vicsek et al. (1995); del Mar Delgado et al. (2018). Animal populations in nature, however, are rarely homogeneous. Within conspecific social groups, heterogeneity may arise from differences in age, size, or sex. Social groups may also have dominance hierarchies including differences in behavioural tendencies such as boldness and shyness Dyer et al. (2008); Ioannou et al. (2008); Michelena et al. (2010); del Mar Delgado et al. (2018). Heterogeneity also arises when individuals of different species interact to form groups, also called mixed-species flocks Morse (1970); Diamond (1981); Sridhar et al. (2009); Stensland et al. (2003); Greenberg (2000); Lukoschek and McCormick (2000). Given the wide prevalence of individual variations among grouping species, it is pertinent to investigate how heterogeneity among individuals influences macroscopic features of collective animal behaviour Gueron et al. (1996); Nagy et al. (2010); Biro et al. (2006); Romey and Vidal (2013); Herbert-Read et al. (2013); Aplin et al. (2014); Farine et al. (2017); del Mar Delgado et al. (2018).

Most of the previous studies that incorporate heterogeneity focus on emergent properties of single groups Gueron et al. (1996); Couzin et al. (2002, 2005); Jolles et al. (2013); Romey and Vidal (2013); Farine et al. (2014); Herbert-Read et al. (2013); Farine et al. (2017); del Mar Delgado et al. (2018). Computational studies show that differences among individuals in phenotypes such as mobility, local cohesion or environmental sensing ability can lead to spontaneous assortment of phenotypes within groups Couzin et al. (2002). For example, individuals with higher speed, or ‘leaders’ who sense environmental gradients, are often at the leading edge of groups despite the absence of any communication or signalling among group members Couzin et al. (2005). Furthermore, even a relatively small proportion of such leaders can facilitate consensus decision making and transfer of information within groups Couzin et al. (2005); Conradt and Roper (2005). Recently, spin-based models have been used to show analytically the existence of phase-transition like behaviour in the consensus decision making in heterogeneous groups Pinkoviezky et al. (2018).

Animal populations across taxa, from insects to mammals, often form a large number of groups that frequently merge (fusion) and split (fission) among themselves Couzin and Laidre (2009). Microbial populations too exhibit such dynamics either because of their self-propulsion or by being driven by their environment Joshi et al. (2017); Durham and Stocker (2012). Previous studies have focussed on deriving the emergence of group size distributions in such fission-fusion populations Gueron and Levin (1995); Gueron (1998); Durrett et al. (1999); Niwa (2003); Ma et al. (2011). The role of heterogeneity, which as discussed above is widely prevalent in natural systems, has not attracted much attention in the literature on fission-fusion systems Sueur et al. (2011); Greenberg (2000). Evolutionary models of collective behaviour predict the emergence of heterogeneity in social, navigational or cooperative traits in fission-fusion populations Guttal and Couzin (2010); Joshi et al. (2017). In such heterogeneous populations, each group needs to be characterised by an additional property that describes the degree of heterogeneity (referred to as group composition). In the literature on mixed-species flocks, group composition patterns are in fact used to infer species-level interactions in ecological communities Sridhar et al. (2012); Berry and Widder (2014); Graves and Gotelli (1993); Ulrich and Gotelli (2010); Sridhar and Shanker (2014). However, group compositions are highly dynamic due to the underlying fission-fusion process among groups.

In this paper, we develop and analyse a model of fission-fusion dynamics of heterogeneous populations. Coagulation-Fragmentation processes provide an excellent mathematical framework to model such flocking dynamics Gueron and Levin (1995); Gueron (1998); Durrett et al. (1999); Niwa (2003); Ma et al. (2011). One such important model, proposed by Niwa Niwa (2003), assumes homogeneous groups on a fixed number of discrete sites. The two most important parameters governing the group movement between sites are the split and move rates. The former determines the rate at which a group splits into two smaller groups (fission), while the latter determines the rate at which a group moves to a new site, merging with any group (fusion) present at the new site. This fission-fusion dynamic model predicts that, in populations of identical individuals, group size distribution is approximately logarithmic. These models have been successful in predicting qualitative features of empirically observed group size distributions from the field Bonabeau et al. (1999); Niwa (2003); Ma et al. (2011). In our study, we employ this framework and generalise it to account for heterogeneity among individuals.

For simplicity, we assume that the population consists of two types of individuals (or species). Unlike homogeneous populations, here we need to keep track of group compositions in addition to the group size distribution. We incorporate the effect of heterogeneity via increased split-rate for groups of heterogeneous composition. The resulting two daughter groups are drawn randomly from all possible partitions of the parent group. We discuss alterations to these assumptions later, but these help keep the model analytically tractable while offering interesting insights on real-world heterogeneous flocks. We first derive master equations for the group sizes and composition and obtain approximate steady-state solutions in the large population limit. We also carry out Monte-Carlo simulations of the model which show considerable agreement with the analytical solution.

Our main finding is that the composition of the flocks depends on the group size. This is despite the merge and split rates being independent of the group size. In particular, we show that there exists a critical group size below which they are more likely to be homogeneous and contain the abundant type/species. However, groups larger than the critical size are representative of the population heterogeneity. The prevalence of heterogeneous groups is surprising, given our assumption that heterogeneous groups exhibit a higher propensity to split. In the Discussion section, we provide a reasoning for this phenomenon. We also discuss some interesting implications of our results for current methods used to infer interspecies interactions from mixed-species flock compositions.

Ii Merge-Split model for Heterogeneous populations

Our formulation of the problem in heterogeneous populations is based on the merge-split model for homogeneous populations, originally conceived by Niwa Niwa (2003) and later analysed by Ma et al Ma et al. (2011) (we present a review of this model in Appendix A). We generalise Niwa’s model to accommodate two species instead of one and derive master equations from the underlying stochastic process. The derivation is non-trivial and includes some assumptions and approximations. We present key steps here, leaving the detailed algebraic steps to the Appendix B.

ii.1 Key Assumptions

We assume sites without geometry and a population consisting of type-I individuals and type-II individuals which can occupy these sites with total population size, . Individuals of the same type are indistinguishable. A group is defined to be the set of individuals occupying the same site at any point in time.

As in the previous model Niwa (2003); Ma et al. (2011), groups move from their current site to a randomly chosen site at a rate that is independent of the size of the groups. If the group lands at a site that is already occupied by another group, they merge to form a larger group with size equal to the sum of the smaller groups.

Unlike the previous model, the groups can be heterogeneous. A group with size of which are of type-I, referred to as the ‘composition’ of a group, will be denoted by the ordered pair . We incorporate the role of heterogeneity via the following assumption: heterogeneous groups have a higher split rate than homogeneous ones. This assumption is justified based on previous individual-based models that indeed predict that more heterogeneous groups are less stable Gueron et al. (1996); del Mar Delgado et al. (2018). More specifically, we assume a group-size-independent split rate which is a function only of the proportion of each type in the group ( and ). The split rate of an -group is given by:


In Eq (1), is the base split rate that is experienced by homogeneous groups (i.e. when or ). The parameter determines excess split rate experienced by heterogeneous groups. The function, is concave down with respect to the proportion , i.e. heterogeneous groups have a higher split rate than homogeneous ones. Groups with proportion experience the maximum split rate, .

When groups do split, they do so uniformly at random, i.e. every possibility that results in two daughter groups is equally probable. Hence, heterogeneous groups are more likely to split but the mechanism of the split does not favour any type of group. A group splits into two groups and such that and , where is the uniform distribution on the integers in the interval . The random variables are sampled conditional on , which ensures that there is a split. After splitting, the two daughter groups occupy random sites.

ii.2 Transition events

The number of -groups at time , denoted by , is the primary random variable of interest. We derive an equation for the rate of change of expected value of this random variable, defined as . This is done by considering all events that will lead to a change in in a small time interval. All such events, along with the resulting change to the number of focal groups , are listed below in Table 1. There are six such events that can lead to a change in the number of focal groups— three merge and three split events.

We denote the rates at which these events occur as follows, with exact expressions for these rates derived in Appendix B.1,

  1. : A merge event changes by .

  2. : A split event changes by .

ii.3 Dynamical equations

Using the above notations for the rates of various events, we obtain the following equation that determines how the expected number of groups of composition changes with time,


This is also known as the master equation. In the large limit it is reasonable to assume that the random variables are pairwise independent. This allows us to rewrite the master equation as (for details see Appendix B.1)


where we remind the reader that is the split rate, is the move rate, and is the number of sites. The notations and represent the maximum and the minimum of and , respectively. Finally is an indicator function defined for a statement as


We also write a mean-field equation by generalising the one in Ma et al. (2011) (Refer to Eq (9b) in Appendix A). The expected total number of groups, , obeys the following equation


ii.4 Steady-state equations

In steady state, i.e.  and , we derive equations relating , the expected total number of groups to , the expected proportion of -groups. When the system size is large (), it is natural to assume that also grows such that the ratio of the two, the fraction of occupied sites, also converges to a constant (). This finally results in the following two equations:


Using an iterative scheme, we solve Eq (6) and Eq (7) to obtain . A detailed description of the derivation, including all the approximations and the iterative technique is provided in Appendix B.

ii.5 Monte-Carlo Simulations

Using a Monte-Carlo algorithm, we simulate the system described above. We maintain a two dimensional counter, that stores the number of groups of size with type-I individuals. At discrete time points, a Bernoulli random variable with appropriately calculated parameter was used to decide between the occurrence of a split and merge. In the case of a split event, the group that undergoes splitting is decided using a bivariate random variable whose probability mass function, is proportional to , where is the split rate. When a group splits, the number of type-I and type-II individuals in the daughter groups is uniformly distributed between 0 and the value for the parent group. Merge events are simulated in an analogous way. The initial condition for the simulation is obtained by placing type-I and type-II individuals uniformly at random on the s sites. After the system reaches steady-state we sample the counter at regular intervals to produce the distribution.

Data and Codes

The source codes for the Monte-Carlo simulation and iterative solutions to (6) and (7) can be found at the following link: Detailed documentation for the Monte-Carlo simulation is also provided in this repository. Data used to plot the figures can also be found here.

Iii Results

To begin, we consider populations with equal proportion of type-I individuals and type-II individuals. From Fig 1, we observe that the results of Monte-Carlo simulations (top row) and iterative solutions to Eq (6)-(7) (bottom row) are in qualitative agreement. We find that groups smaller than a critical size, denoted by , are more likely to be homogeneous, that is, they are dominated by either one of type-I or type-II individuals. Groups larger than the critical size () are usually mixed in equal proportions. We infer this by studying the probabilities


of relative composition of type-I for various group sizes (). For small groups () the distribution is bimodal with modes close to or , suggesting a largely homogeneous composition of groups. As the group size increases to the critical size , the two modes converge to a single mode at (Fig 1), representing the greater tendency of groups to contain both types in equal proportions. The distribution remains unimodal, and thus heterogeneous, for all group sizes (Fig 1). This is surprising given that heterogeneous groups, for all group sizes, have a higher split rate.

It must be noted that the value of is dependent on the excess split rate () and increases as we increase (Fig 1 and 4).

To demonstrate the above transition from homogeneous to heterogeneous groups at a critical group size, in Fig 2, we plot the location of the modes of as a function of group size (). The transition from bimodality to unimodality appears qualitatively similar to a pitchfork bifurcation Strogatz (1994). In this bifurcation, two stable and one unstable fixed points converge to give a single stable fixed point. In our system, the modes (maxima) of the distribution can be viewed as stable fixed points and minima as unstable fixed points.

We show the plots for two cases of unequal abundances of the two types/species in the population in Fig 3. First, when the proportion of type-I () is closer to 0.5, we find that the above results broadly hold true (top row of Fig 3): As the group size () increases, the distribution changes from a bimodal to a unimodal distribution. Unlike the equal proportion scenario where extreme modes merge to form a unimodal distribution (Fig 1), in this case, the two extreme modes vanish with increasing and a mode at the population proportion emerges. Second, when the proportion of type-I () is much smaller than 0.5, the distribution remains unimodal for all group sizes. However, the mode of the distribution gradually moves from an extreme end representing homogeneous groups composed of the abundant species to one representing the population proportion ().

We remark that despite differences in the way the modes of behave for different population proportions of two types, our model predicts a consistent pattern of group-size dependent composition, i.e. small group sizes are likely to be homogeneous with the abundant species whereas larger groups contain two species reflecting population proportion. These surprising qualitative features arise despite simple assumptions of the model such as group-size independent merge and split rates and an excess split rate associated with heterogeneous groups. We provide an intuitive explanation for this in the Discussion section below.

On a similar note, we study as a function of the excess split rate () due to group heterogeneity. We find that when is less than a critical value, the distribution has a single mode at (Fig 4)), representing heterogeneous groups. For above that critical value, however, the distribution becomes bimodal with modes occurring close to and , indicating higher likelihood of homogeneous groups. The location of the modes plotted as a function of also shares qualitative features of a pitchfork bifurcation (Fig 5).

Iv Discussion

In summary, we develop and analyse a heterogeneous flocking model with two types (or species) of individuals. To the best of our knowledge, this is the first model of merge-split dynamics for heterogeneous populations. We use a first principles approach to derive an analytical description of group sizes and composition. We assumed that heterogeneous groups split at higher rates than homogeneous ones but the rates are independent of the size of the groups. Merge rates are independent of both group size and composition. Our key prediction is that composition of small groups is likely to be skewed towards the abundant type. Above a critical group size, , groups reflect the relative composition of species in the population, i.e. they are more likely to be heterogeneous. This is despite the assumption that heterogeneous groups split more often.

We offer an intuitive explanation of the result via two opposing ‘forces’ at play in this model. The first being chance, driven by the number of combinatorial ways a group can be realised by randomly choosing individuals from the population. Given a heterogeneous population, the combinations for the formation of heterogeneous groups far outweigh that of homogeneous ones; this effect is pronounced when the group size is large. The second force which opposes this formation (or maintenance) of heterogeneous group arises from the model assumption that heterogeneous groups are more likely to split into two daughter groups. A single split, however, is not biased towards formation of homogeneous groups. Nevertheless, successive splits have a cumulative effect of homogenising, and reducing the size of daughter groups. Therefore, this homogenising force manifests strongly for smaller group sizes. These forces put together, we find the occurrence of homogeneous groups are dominant up to a critical group size , beyond which the combinatorial forces result in heterogeneous groups.

Our assumption that heterogeneous groups are more likely to split or fragment is supported by previous agent-based simulation models of group movement Gueron et al. (1996); Couzin et al. (2005); del Mar Delgado et al. (2018). Furthermore, these models suggest that, unlike our model assumption, groups do not split uniformly randomly into any of possible partitions; rather, fission events are more likely to cause homogeneous daughter groups. We suspect that incorporating this additional feature, for example in analyses leading to Fig 2, will increase critical group size (), at which the group compositions transition from homogeneous to heterogeneous groups. In other words, we are likely to find groups dominated by the abundant type for much larger groups than predicted by our analyses. Nevertheless, the qualitative features of our results are unlikely to change. Incorporation of these additional features makes our model analytically intractable. Therefore, to confirm these speculations, we suggest studies based on individual-based simulations of fission-fusion group dynamics.

A natural generalisation of our model is one that incorporates species, with . The split rate function for groups could be extrapolated from Eq (1) in a way that preserves its qualitative aspects, i.e. heterogeneous groups having higher split rates than homogeneous ones. For such a system, we expect to find qualitatively similar behaviour to that exhibited by the two species model, i.e. smaller groups are likely to be dominated by one of the species but groups beyond a critical size to be mixed in ratios that are representative of the population composition. To investigate the type of bifurcations and the behaviour of the system near critical points, we require a formal analysis of the generalised model.

We now discuss implications of our results to ecological studies on mixed-species flocks, one of the most widely studied type of heterogeneous flocks. Our model predicts that a study of mixed-species flocks focussing on groups smaller than critical group size of the system will yield observation of flocks that are largely homogeneous; this is despite the fact that the population is heterogeneous. On the other hand, a study on large groups will find flock compositions that represent the population heterogeneity. Therefore, empirical study designs must account for group-size dependent composition of flocks.

The above prediction of our model has further implications for empirical studies that try to infer interspecies interactions from the frequency of their co-occurrence in mixed-species groups. In such studies, typically, a high frequency of co-occurrence beyond what is expected of a null association is typically interpreted as evidence for positive interspecies interactions Sridhar et al. (2012); Berry and Widder (2014); Graves and Gotelli (1993); Ulrich and Gotelli (2010); Sridhar and Shanker (2014). A study that samples flocks that are of size smaller than may rarely find mixed-species associations, thus leading to the conclusion that two species have no or weak positive interspecies interactions. In contrast, a study that samples flocks that are larger than will find many groups with mixed-species associations and thus may arrive at the opposite conclusion of positive interspecies interactions. Therefore, our study highlights that the merge-split dynamics of flocks must be accounted for when making inferences on interspecies interactions.


Our model analysis yields interesting predictions about the composition of heterogeneous groups, suggesting that groups below a certain threshold do not reflect the population level composition. An interesting direction for further study would be to generalise the model of multiple species to allow for differential interactions between species. The differential interactions may arise because the degree of affinities for different pairs of species are not the same. For example, some species may like to be associated with each other while some may avoid each other. Our model provides a starting point to investigate such complex interactions via suitably modified merge and split rate functions. In conclusion, our study highlights the importance of investigating mechanistic models of how individual level interactions between species results in heterogeneous flock dynamics and compositions.

V Acknowledgements

We thank Hari Sridhar for comments on the manuscript. VG acknowledges support from DBT-IISc partnership program, DST Centre for Mathematical Biology at IISc Phase II (SR/S4/MS:799/12) and infrastructure support from DST-FIST.

Appendix A Fission-fusion dynamics of homogeneous populations

Our formulation of the problem in heterogeneous populations is based on the model originally conceived by Niwa Niwa (2003) and later analysed by Ma et al Ma et al. (2011). To keep the paper self contained as well as keep it easier to understand the more involved derivations for the heterogeneous case, we review the homogeneous model and its assumptions in this appendix. The model assumes sites with no geometry and a population of indistinguishable individuals which can occupy these sites. A group is defined to be the set of individuals occupying the same site at any point of time. All groups move at rate (which will often be referred to as merge rate) and split at rate . These rates are independent of group size. When a group moves to an occupied site, they merge to form a larger group with size equal to the sum of smaller groups. A split results in the formation of smaller groups that move to random empty sites. The model can be thought of as the coarse-grained version of a microscopic model with local interactions. Ma et al. derived deterministic evolution equations for the merge-split model from first principles. They did so by considering the various changes that could happen to , the expected number of groups of size at time , in a small time interval . Their analysis of the above described merge and split processes resulted in the following coupled differential equations, given below, where denotes the total number of groups at a given time.


where the symbol denotes the indicator function, which is 1 when ‘expr’ is satisfied and 0 otherwise. We describe the dynamics that each term of Eq (9a) and Eq (9b) represent and a few minor modifications that we propose for better accuracy. The first term captures the event where groups of size and merge, but in the case of , we need to account for over-counting. Furthermore, since a group cannot merge with itself, the term corresponding to has to be .

The second term in Eq (9a) corresponds to larger groups splitting to form groups of size . The factor of 2 accounts for the fact that a group of size can split in two equally probable ways to yield a group of size . Since the system is finite, the upper limit in the sum cannot be infinity.

The third term represents the probability of groups with size merging with other groups. This can happen in two ways- either the group of size moves to an occupied site or a group moves to a site occupied by an -sized group, hence the factor of two. This term, however, includes merger with groups of size greater than , which is impossible. When we generalise this model to heterogeneous groups, we resolve this issue.

The fourth term is the decrease in , due to a group of size splitting, and does not require modification.

Eq (9b), also called the mean-field equation can be obtained by considering the processes that lead to a change in the total number of groups, . Each split event can increase by 1. Since all groups of split at rate , the term is is the rate at which increases. Each merge event, on the other hand, decreases by 1. Since merge events happen when groups move (at rate ) to already occupied sites, the total rate associated with merge events is . The solutions to these equations shows that fission-fusion dynamics approximately yields a logarithmic group size distribution. We generalise this model and its analytical formulation to heterogeneous populations in the main text.

Appendix B Deriving the steady state distribution function

Here, we present the complete derivation of the steady state equations, Eq (6) and Eq (7) from the main text.

b.1 The Master Equation

Groups move from site to site (and effectively merge) at rate . They split at rate . The primary random variable of concern is , the number of groups with size , of which are of type-I. Let } denote the -algebra generated by events up to time (which contains all the information about the dynamics of the system up to time ). Now, we find the expected number of groups of size with individuals of type-1 at time given information up to time ,


where .

The terms in Eq (10) arise from the following considerations:

  1. : A group, moves to a site occupied by a group and they merge to produce a larger group , thus increasing by 1. The rate for the event wherein a group moves is . The probability that a site is occupied by a group is . Hence the rate for the entire event is given by the product summed over all the possible values of and . When and are even, there is a corner case where two identical groups merge. Since a group cannot merge with itself the rate will be . This yields


    The second term is accompanied by (where means that is divisible by ) imposes the condition that only groups with even and can be formed by the merger of identical groups (otherwise would not make sense, since and would not be integers.

    The limits of the sums have been chosen very carefully to account for the finite size of the population. The outer sum goes from to , since merge events always result in an increase in group size, groups with size cannot merge and form a group of size .

    The inner sum starts at (where ). If , then , which would make meaningless (since a group can’t have more type-I individuals than its size). So when we start at instead of .

    Similarly, the upper limit is (where ) because groups with more than type-I individuals can never produce through a merge event. As long as we consider mergers involving groups with up to type-I individuals, but for we restrict the sum to .

  2. : A group merges with a group of another size and composition, and decreases the count by one. This term is calculated in a similar way to except that the group under consideration is merging. It is important to note that the merge can happen in two ways- the group moves to a site occupied by , and vice-versa. The multiplicative factor of 2 in the first term is to account for this. As reasoned earlier, the second term in the parenthesis accounts for the fact that groups cannot merge with themselves. When a group does merge with an identical group, the change in is -2 and not -1, so we need to exclude these kind of events from . The last term that is subtracted accounts for this fact.


    Groups of size larger than cannot merge with other groups of the same size, since there cannot be more than one group of size greater than . Even in groups of size lesser than , two groups cannot have more than or less than type-I individuals. The two indicator functions multiplied with the second term in parenthesis ensures that these constraints are not ignored.

    We are considering events where groups merge with other groups, therefore groups cannot have sizes more than , (the population is finite). The limits of the inner sum, analogously to the one in , restricts the compositions of the groups that can merge. Say , then , which is impossible, since it would imply that the right hand side of the last inequality, which is the number of type-II individuals in the resulting group is more than total type-II population, .

  3. : Two identical groups merge to give a group , which results in decreasing by 2. A group moves at rate and lands on a site occupied by an identical group with probability , to yield the expression

  4. : A larger group splits into and increasing by one. The mechanism of splitting is uniformly at random; consider a group ; the type-I individuals arranged linearly can be split at different points and similarly the type-II individuals can be split in ways. Combining one part from each of these groups yields the two groups resulting from splitting. However, splits that produce empty groups are not allowed, and this can happen in two extreme cases. Thus the subtraction of 2 from the total is necessary, resulting in combinations. Since splits into and in two symmetric ways, a multiplicative factor of 2 will arise in the expression.


    The second term accounts for the corner case where a group splits equally to give -groups. In this case, increases by 2 and not 1, and is accounted for in the term below.

    In the first term, the index of the outer sum goes from to , since only groups of size larger than can split to give groups of size . The inner sum starts at because when , due to the finite number of type-II species, the parent group has to have at least type-I individuals. The upper limit is , because , which would imply that the daughter group has more type-II individuals than the parent , which is impossible.

  5. : The event where a group splits to give two identical groups , increasing by 2. The second term in Eq (14) prevents double counting of this event.

  6. : Group splits. This is the most straightforward of all the events, with a rate given by


We can rewrite Eq (10) as follows:


Now take the expectation on both sides of Eq (17) and in the continuum limit (i.e. letting ), we get


b.2 Steady-State Equation

To estimate the expected number of groups at steady state, we set , to obtain


In Eq (19), and , are defined in the same way as before, except that is replaced by , the stationary distribution of the Continuous time Markov Chain .

To proceed further, we assume that is large enough so that the random variables are pairwise independent. So if we set , then we obtain, from Eq (11)-(16),Eq (19),


b.3 Mean-Field Equation

We can also derive an equation for the total number of groups, by generalising Eq (9b). Each split event increases by 1 and each merge event decreases by 1.

The rate of increase of due to all the current groups is , hence the total contribution of split events to is this term summed over all valid ’s and ’s.
At a given time, groups merge with at a rate , so this term summed over all possible ’s, ’s, ’s, and ’s gives the total rate of decrease of .


An additional term is subtracted in the second sum to account for the fact that groups cannot merge with themselves.

b.4 Scaling Limit

Dividing Eq (20) by we can rewrite the equation in terms of (the proportion of -groups), we get


We also divide the steady-state () version of Eq (21) by ,


Since is large and , we must speed up the split and merge rates to obtain a non-trivial limit and consequently set the split and merge rates to and in Eq (22) and Eq (23). In this limit, we can use the a priori knowledge that the mass of the group distribution is concentrated at small , to remove the constraints on the last two terms of Eq (22). Consequently, the double sum yields (since is normalised). We drop the constraint in the last term as well since as we will show below this term will be negligible in the limit and can be ignored.

Using these approximations in Eq (22) gives


We assume that when the system scales, also scales such that , (the fraction of occupied sites in steady-state) and consequently we can write the steady state equation for as


We also incorporate the scaling assumption for and into the steady-state mean-field equation, Eq (23) which gives


In the large and limit, the second term in parentheses goes to 0. Relaxing the finite size constraint on the first term in parentheses will simplify the equation considerably, since , leaving only . Again assuming that , we can write


Remarks: It is important to note that, although for , the split rate reduces to the homogeneous case, will not obey the homogeneous master equation Eq (9a). This is because in Eq (2), from the main text, the terms and are derived by assuming that individuals of the same type are indistinguishable, but of different types are distinguishable. However, in the homogeneous case, all individuals are indistinguishable. Hence, if both the equations were derived assuming that all individuals are distinguishable, will obey the homogeneous equation also.

b.5 Iterative Solution

It is a challenging task to solve Eq (26) and Eq (27) directly, so we use an iterative scheme to obtain the solutions. Eq (26) can be used to express the proportion of groups of a particular size and composition in terms of the proportions of groups of other sizes and compositions.


is the index of iteration. As the initial condition for the fixed point iteration we use the solution to the homogeneous equation (). We distribute uniformly to , for . This set of iterative equations goes to steady state, and produces stable solutions. However, as a result of the drastic approximations we adopted, the total population size is not conserved during the iteration.

Appendix C Group size distributions

Earlier studies that adopted merge-split dynamics Niwa (2003); Ma et al. (2011) primarily investigated the group size distributions in homogeneous populations. We found it instructive to look at the group size distribution for heterogeneous populations too. The probability of a group having size is obtained by summing the composition dependent proportion, over all possible compositions, resulting in a group size probability defined by . Fig 6 shows as a function of on log-log scale from simulations and iterative solution to the analytical equations. The plots shows qualitative match with the earlier predicted distributions and is approximately logarithmic. Although the likelihood of occurrence of small groups is nearly the same for different values of delta, the decays much faster for larger values of delta. This means that large groups are rarer for higher values of .


  • Parrish and Edelstein-Keshet (1999) J. K. Parrish and L. Edelstein-Keshet, Science 284, 99 (1999).
  • Camazine (2003) S. Camazine, Self-organization in biological systems (Princeton University Press, 2003).
  • Sumpter (2010) D. J. Sumpter, Collective animal behavior (Princeton University Press, 2010).
  • Couzin and Krause (2003) I. D. Couzin and J. Krause, Advances in the Study of Behavior 32, 1 (2003).
  • Vicsek et al. (1995) T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen,  and O. Shochet, Physical review letters 75, 1226 (1995).
  • del Mar Delgado et al. (2018) M. del Mar Delgado, M. Miranda, S. J. Alvarez, E. Gurarie, W. F. Fagan, V. Penteriani, A. di Virgilio,  and J. M. Morales, Phil. Trans. R. Soc. B 373, 20170008 (2018).
  • Dyer et al. (2008) J. R. Dyer, D. P. Croft, L. J. Morrell,  and J. Krause, Behavioral Ecology 20, 165 (2008).
  • Ioannou et al. (2008) C. Ioannou, M. Payne,  and J. Krause, Oecologia 157, 177 (2008).
  • Michelena et al. (2010) P. Michelena, R. Jeanson, J.-L. Deneubourg,  and A. M. Sibbald, Proceedings of the Royal Society of London B: Biological Sciences 277, 1093 (2010).
  • Morse (1970) D. H. Morse, Ecological monographs 40, 119 (1970).
  • Diamond (1981) J. M. Diamond, Nature 292, 408 (1981).
  • Sridhar et al. (2009) H. Sridhar, G. Beauchamp,  and K. Shanker, Animal Behaviour 78, 337 (2009).
  • Stensland et al. (2003) E. Stensland, A. Angerbjörn,  and P. Berggren, Mammal Review 33, 205 (2003).
  • Greenberg (2000) R. Greenberg, On the move: How and why animals travel in groups , 521 (2000).
  • Lukoschek and McCormick (2000) V. Lukoschek and M. McCormick, in Proceeding 9th International Coral Reef Symposium, Vol. 1 (Ministry of Environment, Indonesian Institute of Sciences and International Society for Reef Studies, 2000) pp. 467–474.
  • Gueron et al. (1996) S. Gueron, S. A. Levin,  and D. I. Rubenstein, Journal of Theoretical Biology 182, 85 (1996).
  • Nagy et al. (2010) M. Nagy, Z. Akos, D. Biro,  and T. Vicsek, Nature 464, 890 (2010).
  • Biro et al. (2006) D. Biro, D. J. Sumpter, J. Meade,  and T. Guilford, Current Biology 16, 2123 (2006).
  • Romey and Vidal (2013) W. L. Romey and J. M. Vidal, Ecological modelling 258, 9 (2013).
  • Herbert-Read et al. (2013) J. E. Herbert-Read, S. Krause, L. Morrell, T. Schaerf, J. Krause,  and A. Ward, Proc. R. Soc. B 280, 20122564 (2013).
  • Aplin et al. (2014) L. M. Aplin, D. R. Farine, R. P. Mann,  and B. C. Sheldon, Proc. R. Soc. B 281, 20141016 (2014).
  • Farine et al. (2017) D. R. Farine, A. Strandburg-Peshkin, I. D. Couzin, T. Y. Berger-Wolf,  and M. C. Crofoot, Proc. R. Soc. B 284, 20162243 (2017).
  • Couzin et al. (2002) I. D. Couzin, J. Krause, R. James, G. D. Ruxton,  and N. R. Franks, Journal of theoretical biology 218, 1 (2002).
  • Couzin et al. (2005) I. D. Couzin, J. Krause, N. R. Franks,  and S. A. Levin, Nature 433, 513 (2005).
  • Jolles et al. (2013) J. W. Jolles, A. J. King, A. Manica,  and A. Thornton, Animal Behaviour 85, 743 (2013).
  • Farine et al. (2014) D. R. Farine, L. M. Aplin, C. J. Garroway, R. P. Mann,  and B. C. Sheldon, Animal behaviour 95, 173 (2014).
  • Conradt and Roper (2005) L. Conradt and T. J. Roper, Trends in ecology & evolution 20, 449 (2005).
  • Pinkoviezky et al. (2018) I. Pinkoviezky, I. D. Couzin,  and N. S. Gov, Physical Review E 97, 032304 (2018).
  • Couzin and Laidre (2009) I. D. Couzin and M. E. Laidre, Current biology 19, R633 (2009).
  • Joshi et al. (2017) J. Joshi, I. D. Couzin, S. A. Levin,  and V. Guttal, PLoS computational biology 13, e1005732 (2017).
  • Durham and Stocker (2012) W. M. Durham and R. Stocker, Annual review of marine science 4, 177 (2012).
  • Gueron and Levin (1995) S. Gueron and S. A. Levin, Mathematical biosciences 128, 243 (1995).
  • Gueron (1998) S. Gueron, Journal of Mathematical Biology 37, 1 (1998).
  • Durrett et al. (1999) R. Durrett, B. L. Granovsky,  and S. Gueron, Journal of Theoretical Probability 12, 447 (1999).
  • Niwa (2003) H.-S. Niwa, Journal of Theoretical Biology 224, 451 (2003).
  • Ma et al. (2011) Q. Ma, A. Johansson,  and D. J. Sumpter, Journal of theoretical biology 283, 35 (2011).
  • Sueur et al. (2011) C. Sueur, A. J. King, L. Conradt, G. Kerth, D. Lusseau, C. Mettke-Hofmann, C. M. Schaffner, L. Williams, D. Zinner,  and F. Aureli, Oikos 120, 1608 (2011).
  • Guttal and Couzin (2010) V. Guttal and I. D. Couzin, Proceedings of the national academy of sciences 107, 16172 (2010).
  • Sridhar et al. (2012) H. Sridhar, U. Srinivasan, R. A. Askins, J. C. Canales-Delgadillo, C.-C. Chen, D. N. Ewert, G. A. Gale, E. Goodale, W. K. Gram, P. J. Hart, et al., The American naturalist 180, 777 (2012).
  • Berry and Widder (2014) D. Berry and S. Widder, Frontiers in microbiology 5 (2014).
  • Graves and Gotelli (1993) G. R. Graves and N. J. Gotelli, Proceedings of the National Academy of Sciences 90, 1388 (1993).
  • Ulrich and Gotelli (2010) W. Ulrich and N. J. Gotelli, Ecology 91, 3384 (2010).
  • Sridhar and Shanker (2014) H. Sridhar and K. Shanker, Behavioral ecology and sociobiology 68, 185 (2014).
  • Bonabeau et al. (1999) E. Bonabeau, L. Dagorn,  and P. Fréon, Proceedings of the National Academy of Sciences 96, 4472 (1999).
  • Strogatz (1994) S. H. Strogatz, Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering (Westview press, 1994).


Figure 1: Flock composition depends on flock size in heterogeneous populations. We find qualitative agreement between simulation results ((a)-(e), top row) and the iterative solutions ((f)-(j); bottom row) of the analytic model equations (6) and (7). We represent flock composition for each group size by the probabilities (; see Eq (8)) as a function of the relative composition of type-I individuals (). Here, we assumed that two types are equally abundant in the population (). For small , shown in (a)-(c) and (f)-(h), is a bimodal function with modes occurring away from the population ratio of type-I to type-II (which is 0.5); this suggests that small group sizes are dominated by one or the other type/species. For large above a critical value, shown in (e) and (j), is unimodal at