# Modeling the Origin of Eusociality

A 2010 paper***** that challenged a decades-old explanation for the evolution of “eusocial organization” in many insects (ants, wasps, bees, termites, thrips, aphids) and certain other animals (snapping shrimp and naked mole rats) caused a stir among evolutionary biologists. The distinguishing feature of a eusocial population is that some individuals reduce their own lifetime reproductive potential to help raise the offspring of others. Eusociality is rare (occurring in only 2% of the 900,000 known species), but successful (social insects account for more than 50% of the known biomass of insects). The controversy engendered by the paper was no doubt enhanced by the identity of one of the authors: Edward O. Wilson, one of the early and strongest proponents of the challenged explanation, which is called “inclusive fitness.” Wilson’s co-authors, Martin Nowak and Corina Tarnita, are mathematicians with appointments in the Program for Evolutionary Dynamics at Harvard. Nowak, in the keynote address at the Mathematical Biosciences Institute’s 10th-anniversary celebration, discussed some of the ideas and the mathematical model presented in the 2010 paper.

Altruistic behavior among bees was recognized by Darwin as a conundrum for his theory of evolution; he suggested that the whole colony was, in some sense, the unit of selection. Slightly more than a century later, the biologist William D. Hamilton, building on work of J.B.S. Haldane, Ronald Fisher, and Sewall Wright, provided an explanation and formalized a foundation for a theory of “inclusive fitness” using kin selection. Hamilton defined inclusive fitness as the personal fitness expressed by an individual in its production of offspring, stripped of all components that could be attributed to its environment. This fitness is augmented by certain quantities of harm and benefit that the individual causes to the fitness of neighbors by his actions. The quantities, expressed in fractions, are coefficients of relatedness—one for clonal individuals, one-half for siblings, one-quarter for half-siblings, one-eighth for cousins, and so forth. Hamilton codified the theory in an inequality: \(R > c/b\), where \(R\) is relatedness, \(c\) the cost, and \(b\) the benefit. An especially trenchant case is that of haploid–diploid societies, in which females are produced by fertilization of eggs, males by unfertilized eggs. In these societies, sisters are related more closely to each other (\(R\) = three-quarters) than to their mother (\(R\) = one-half). The classic example of such societies is found in the *Hymenoptera *(an order of ants, bees, and wasps).

In their paper, Nowak et al. show that inclusive fitness holds only in special cases and describe the mathematical model they developed. Based on natural selection, the model is more general than inclusive fitness and agrees with the results for inclusive fitness in the special cases in which it works. They identify weaknesses of inclusive fitness: It holds only in the weak selection limit (all individuals have approximately the same fitness, and the two social strategies (eusocial and solitary) are roughly equally abundant); the interactions between individuals are additive and pairwise; and the theory can deal only with special population structures. (This paper is not the first to point out the limitations of inclusive fitness based on mathematical arguments. In a 1984 paper in the *Proceedings of the National Academy of Sciences*, for example, biologist Carlos Matessi and mathematician Samuel Karlin concluded that Hamilton’s rule has quantitative validity only in the special case of linear fitness functions.)

Nowak et al. begin with a simple example, which they then extend to the haploid–diploid case. The simpler model, for asexual reproduction, is a deterministic model for evolutionary dynamics described by an ordinary differential equation. In their solitary model

\[\begin{equation}\tag{1} \frac{dx_0}{dt} = (b_0 - d_0)x_0,\end{equation}\]

\(x_0\) is the number of solitary queens, and \(b_0\) is the reproduction rate and \(d_0\) the death rate of solitary queens.

They compare that model with a eusocial model,

\[\begin{equation}\tag{2}

\frac{dx_1}{dt}= \sum^{\infty}_{i=1}{b_i}(1-q)x_i-b_1qx_1-d_1x_1\\

\frac{dx_i}{dt}= b_{(i-1)}qx_{(i-1)}-b_iqx_i-d_ix_i,

\end{equation}\]

in which \(x_i\) is the number of colonies of size \(i\); \(b_i\) is the reproductive rate and \(d_i\) the death rate of a queen in a colony of size \(i\); and \(q\) is the probability that a daughter stays in the colony, which can be considered an effect of a mutation.

A eusocial strategy will be selected over a solitary strategy provided that the maximum eigenvalue, \(\lambda\), of the matrix associated with the right-hand side of the eusocial model satisfies \(\lambda > b_0 – d_0\). This depends, in turn, on how the demographic parameters of the queen change as the size of the colony increases. If for a colony of size \(m\), \(i < m\), the reproductive rate and death rate of the eusocial queen remain the same as those of the solitary queen, but for \(i > m\) the eusocial queen has a reproductive rate of at least \(b > b_0\) and a death rate of at most \(d < d_0\), then a necessary condition for the selection of eusociality is

\[\begin{equation}\tag{3}

b > k_m (b_0 + (d – d_0)).

\end{equation}\]

The value \(k_m\) grows exponentially with \(m\). For small colonies to evolve eusociality, the reproductive rate of the queen must be increased.

This model can be extended to account for density limitations on reproductive rate and worker mortality—in the former, by multiplying the values of \(b\) by a factor \(\gamma = 1/(x_0 + \eta X)\), where \(X\) is the total population size, and in the latter, by adding the following terms to the equations (2):

\[\begin{equation}

\alpha x_2~~ \mathrm{for}~~ i = 1; \\

–\alpha(i – 1)x_i + \alpha ix_{(i+1)} ~~\mathrm{for}~~ i > 1,

\end{equation}\]

where \(\alpha\) is the worker death rate.

For the parameter values \(b_0 = .5\), \(d_0 = .1\), \(m = 3\), \(b = 4\), \(d = .01\), \(\eta = .01\), and \(\alpha = .1\), eusociality is selected for an intermediate range of values \((.36 < q < .9)\). For low values of \(q\), there is a small chance that a colony will reach critical size; for large values of \(q\), the colonies produce too few new queens.

The authors build a model for the key case of sexual reproduction with haploid–diploid genetics. They consider a dynamic evolution model to study competition between two alleles, a wildtype allele denoted “A” and a mutant allele denoted “a” that causes daughters to stay in the colony with some probability. There are then three different types of females, AA, Aa, and aa, with probabilities of staying in the colony 0, \(q_2\), and \(q_3\), respectively. If the mutant allele a is dominant, \(q = q_2 = q_3 > 0\); if it is recessive, only \(q = q_3\) is positive. There are two types of males, A and a, and as a result six types of fertilized queens, AA-A, AA-a, Aa-A, Aa-a, aa-A, and aa-a, where the last letter indicates a wildtype or mutant mate.

The resulting system of equations is similar to the one for the simple case. There are equations for the number of colonies of size \(i\) associated with each type of queen, supplemented by three equations, one for the number of each type of virgin female and two for the number of each type of male. These latter equations require three parameters \((g, h, \beta)\) that specify the death rate of females, the death rate of males, and the rate of successful mating. It is assumed that the males do not die after mating, but have a short lifespan and are thus unlikely to mate more than once.

Assuming that the initial state is dominated by the solitary allele A, Nowak et al. simulate the system for the same parameter values considered in the above example and for values \(g = h = .1\) and \(\beta = .1\). At equilibrium, they ask, which allele wins the competition as a function of \(q\)? If a is dominant, then for \(.26 < q <.88\), the eusocial allele wins and for \(.88 < q\) there is coexistence of alleles. If a is recessive, the solitary allele A wins the competition as long as \(q < .7\); for \(.7 < q < .88\), the eusocial allele a wins, and for \(.88 < q\), the two alleles coexist.

Based on their model, the authors make the following observations: First, it is difficult to evolve eusociality because of the dependence on favorable parameters—e.g., an increased birth rate. Second, the haploid–diploid model is bistable, and is more easily maintained than evolved. Third, inclusive fitness is not needed to explain eusociality. Finally, once eusociality is evolved, relatedness becomes a consequence of eusociality, rather than a cause of it.

Evolutionary biologists continue to pose objections to the paper, none of which address the model. The future would seem to hold opportunities for refining evolutionary models and gathering the experimental data to support model parameters.

*****M.A. Nowak, C.E. Tarnita, and E.O. Wilson, *The evolution of eusociality*, Nature, 46 (2010), 1057–1062.