## Abstract

During bouts of evolutionary diversification, such as adaptive radiations, the emerging species cluster around different locations in phenotype space. How such multimodal patterns in phenotype space can emerge from a single ancestral species is a fundamental question in biology. Frequency-dependent competition is one potential mechanism for such pattern formation, as has previously been shown in models based on the theory of adaptive dynamics. Here, we demonstrate that also in models similar to those used in quantitative genetics, phenotype distributions can split into multiple modes under the force of frequency-dependent competition. In sexual populations, this requires assortative mating, and we show that the multimodal splitting of initially unimodal distributions occurs over a range of assortment parameters. In addition, assortative mating can be favoured evolutionarily even if it incurs costs, because it provides a means of alleviating the effects of frequency dependence. Our results reveal that models at both ends of the spectrum between essentially monomorphic (adaptive dynamics) and fully polymorphic (quantitative genetics) yield similar results. This underscores that frequency-dependent selection is a strong agent of pattern formation in phenotype distributions, potentially resulting in adaptive speciation.

## 1. Introduction

Explaining the origin of diversity is a core problem in evolutionary biology that continues to receive much attention from both empiricists and theoreticians (Coyne & Orr 2004; Dieckmann *et al*. 2004). The process of diversification can be described as evolutionary change taking place in phenotype space. If individual organisms are assessed for their phenotypes, populations can be represented by the corresponding phenotype distributions, giving information about the abundance of different phenotypes in the population. A single, ancestral population would typically yield a unimodal phenotype distribution, with the average phenotype being at or close to the distribution's peak. Processes of speciation can then often be described as the splitting of an ancestral and unimodal phenotype distribution into two (or more) peaks or modes, so that the descendent species emerging from speciation correspond to different peaks of the phenotype distribution. On the phenotypic level, speciation can thus cause pattern formation: during speciation, unimodal phenotype distributions may become multimodal.

Traditional explanations for such pattern formation through speciation are based on geographical isolation: different, but phenotypically similar subpopulations of an ancestral species come to occupy different and mutually isolated habitats in which they embark on different evolutionary trajectories. These trajectories may eventually take the populations evolving in different habitats to different locations in phenotype space so that the joint phenotype distribution of all descendent species becomes multimodal. It is important to appreciate that this allopatric mode of phenotypic pattern formation and speciation results from geographical isolation, rather than from ecological interactions within the ancestral population.

The situation is reversed for sympatric processes of speciation, which unfold due to ecological interactions within the ancestral population, rather than as a consequence of geographical isolation. For example, when phenotypes differ in their resource preference and when most individuals in an ancestral population prefer similar resources, selection may favour rare phenotypes with a different resource preference. In this case, diversification of the ancestral population may be an adaptive response to the detrimental effects of frequency-dependent competition. In general, adaptive speciation occurs when an ancestral lineage splits into phenotypically diverging descendent lineages due to disruptive selection caused by frequency-dependent interactions (Dieckmann *et al*. 2004). In this mode of speciation, pattern formation in phenotype space is caused by interactions that are intrinsic to the ancestral population. The theoretical framework of adaptive dynamics predicts that such adaptive diversification can occur under a wide variety of ecological scenarios (Metz *et al*. 1996; Geritz *et al*. 1998; Dieckmann & Doebeli 1999; Doebeli & Dieckmann 2000, 2003; Dieckmann *et al*. 2004; Kisdi & Gyllenberg 2005). In this framework, adaptive diversification is epitomized by the phenomenon of evolutionary branching, which occurs when frequency-dependent selection drives a population towards a point in phenotype space at which selection turns disruptive. Evolutionary branching can be characterized mathematically and is a generic outcome of adaptive dynamics models (Metz *et al*. 1992, 1996; Geritz *et al*. 1998; Kisdi & Gyllenberg 2005).

Most models of evolutionary branching are based on a number of seemingly significant simplifying assumptions. Chief among those are the assumptions that reproduction is asexual and that populations are essentially monomorphic at all times (except when branching occurs, after which each of the emerging lineages is assumed to be essentially monomorphic). Obviously, both of these assumptions are often violated in real populations. It is thus important that it has been shown that evolutionary branching is also a robust outcome in asexual models of polymorphic populations (Metz *et al*. 1996; Meszéna *et al*. 2005) and that a number of recent models have incorporated explicit genetics to study adaptive speciation in sexual populations (Doebeli 1996, 2005; Dieckmann & Doebeli 1999; Kisdi & Geritz 1999; Kondrashov & Kondrashov 1999; Drossel & McKane 2000; Geritz & Kisdi 2000; Doebeli & Dieckmann 2003; Dieckmann *et al*. 2004; Bürger & Schneider 2006; Schneider & Bürger 2006; Bürger *et al*. in press). In sexual populations under disruptive selection, random mating typically prevents speciation, so that diversification requires the presence of assortative mating mechanisms ensuring that individuals preferentially mate with similar phenotypes. Such mechanisms have been considered in models with genetic architectures based on small to intermediate numbers of loci with additive effects (see articles cited above). One general conclusion of such studies is that adaptive speciation, or pattern formation in phenotype space, is possible in sexual populations when mating is assortative.

It has recently been suggested by Polechová & Barton (2005) that occurrences of adaptive speciation in sexual populations could often be a consequence of the particular genetic models used and that other genetic models would not generate diversification in sexual populations even with assortative mating. One reason for this caveat might be that in genetically explicit models with a finite number of loci and with finite allelic effects, a population's variance is automatically constrained, leading to more intense intraspecific competition and thus strengthening disruptive selection. In models with more flexible genetic architectures, intraspecific competition might simply result in increased population variance. In particular, for populations described by a continuous phenotype distribution (rather than by a single monomorphic type or by the frequencies of a finite number of types), one might have the intuitive expectation that frequency-dependent competition merely flattens unimodal phenotype distributions, thus compensating for the effects of competition. According to this intuition, frequency-dependent selection would not result in pattern formation in phenotype space, i.e. in a bimodal or multimodal split of the phenotype distribution and hence would not result in adaptive speciation.

To describe the dynamics of continuous phenotype distributions under frequency-dependent competition, Polechová & Barton (2005) used the ‘infinitesimal’ model of quantitative genetics, which assumes a large (infinite) number of unlinked loci with additive effects (Roughgarden 1979; Bulmer 1980). Polechová & Barton (2005) claim that in such models, frequency dependence never leads to adaptive speciation even if mating is assortative. This would support the intuitive notion that frequency dependence can generate increased population variance, but not pattern formation in the form of bimodal or multimodal phenotype distributions. These questions are very interesting and deserve further study. In this paper, we use a more general class of models to show that frequency-dependent competition in sexual populations indeed leads to pattern formation in phenotype space under many circumstances.

The intuitive notion that in models for continuous phenotype distributions, frequency dependence only leads to increased variance, but not to phenotypic clusters, thus turns out to be wrong in general. Instead, if mating is assortative, frequency-dependent competition often generates multiple phenotypic modes also in infinitesimal models. Since the splitting of a population reduces the strength of disruptive selection, assortative mating facilitates the evolutionary response to frequency dependence. Consequently, there is selection for assortative mating in initially randomly mating populations, in which segregation and recombination would otherwise prevent the emergence of multiple modes. This is why pattern formation in phenotype space is a possible outcome of frequency-dependent competition in infinitesimal models of sexual populations.

Our results show that, with regard to adaptive diversification, the outcomes of asexual adaptive dynamics models at one end of the spectrum, and of infinitesimal sexual models at the other end, are surprisingly congruent. In the sexual models, assortative mating is required for adaptive speciation to occur, but in both types of model, the emergence of distinct phenotypic clusters out of unimodal or even monomorphic ancestral populations can readily be caused by frequency-dependent ecological interactions. This pattern formation in sexual models could be an important mechanism underlying the instability and disruption of the sexual continuum of phenotypes (Maynard Smith & Szathmáry 1995; Noest 1997) and hence could help address the fundamental question why life forms appear to cluster phenotypically (Coyne & Orr 2004).

## 2. Model description

Below, we introduce the dynamics of the density distribution *ϕ*(*x*) of a quantitative character *x* in a sexual population.

### (a) Ecological dynamics

The ecological model underlying our analysis is an extension of Lotka–Volterra competition equations to polymorphic populations, in which the competitive impact of phenotype *y* on a phenotype *x* is measured by the competition kernel *α*(*x*−*y*). For a focal phenotype *x*, the total competitive impact experienced in a population described by the distribution *ϕ* is given by the convolution(2.1)In the asexual case, the dynamics of the distribution *ϕ* is then given by the following partial differential equation:(2.2)Here, *r* is the intrinsic growth rate, which we assume to be independent of the phenotype *x*, and *K*(*x*) determines the carrying capacity as a function of *x*. Thus, −*rϕ*.*α***ϕ*/*K* corresponds to the usual competition term in Lotka–Volterra models, whereas *rϕ* describes exponential population growth.

### (b) Mating and reproduction

We incorporate sexual reproduction following standard procedures (Roughgarden 1979; Bulmer 1980; see also Polechová & Barton 2005). We assume that matings are initiated bilaterally. The probability of mating between two phenotypes *u* and *v* is therefore proportional to the product of two preference functions, which we assume to be Gaussian,(2.3)where *σ*_{A} is a measure for the degree of assortment: large *σ*_{A} correspond to random mating, while small *σ*_{A} correspond to assortative mating (note that we will always assume here that assortative mating occurs with respect to the quantitative character that determines the ecological interactions).

In accordance with the assumptions underlying the infinitesimal model of quantitative genetics (Bulmer 1980), we assume that a mating between phenotypes *u* and *v* produces a Gaussian offspring distribution , with a mean equalling the midparent value (*u*+*v*)/2 and a variance of .

To establish a baseline case, we assume that all phenotypes have the same per capita birth rate. This means that the relative contribution a mating with phenotype *v* makes to the offspring pool of a given phenotype *u* must be normalized by the total amount of mating that phenotype *u* participates in,(2.4)Then, the distribution of offspring with phenotypes *x* produced by phenotype *u* is given by(2.5)Finally, the total density of offspring at phenotype *x* resulting from all possible matings is given by(2.6)The mating scheme just described for the infinitesimal model is a direct extension of the one used in Dieckmann & Doebeli (1999) for genetically explicit multilocus models.

Putting everything together, we obtain the following equation for the dynamics of phenotype distributions in sexual populations:(2.7)The essential parameters in this dynamical system are *σ*_{A} (degree of assortment) and *σ*_{f} (width of the so-called segregation kernel (Roughgarden 1979) that describes the offspring distribution of a given mating pair), as well as the functional forms of the ecological functions *α* and *K*. For numerical simulations of the partial differential equation (2.7), we always used carrying-capacity functions *K* with finite variance, which implies that phenotypes far from the optimal phenotype are not viable. This allows the numerical simulations to be restricted to a finite interval without creating artefacts.

### (c) Competition kernel and carrying-capacity function

It is already very interesting to study the dynamics of the asexual model, equation (2.2), which is determined by the functions *α* and *K*. In particular, one can ask whether, for given functions *α* and *K*, equilibrium distributions of the asexual model exhibit phenotypic clustering in the form of multiple modes. For example, if the competition kernel *α* and the carrying capacity *K* are both of Gaussian type with variances and , respectively, then the model has an equilibrium density distribution that is also Gaussian, with variance (if is negative, the equilibrium distribution has all its density concentrated at the maximum of *K*). In particular, with Gaussian *α* and *K*, equilibrium distributions of the asexual model never exhibit more than one phenotypic cluster.

However, it is known that the asexual model with Gaussian ecological functions is structurally unstable (Sasaki & Ellner 1995; Sasaki 1997) and that generic choices for the ecological functions often lead to pattern formation with distinct phenotypic clusters (Meszéna *et al*. 2005). We therefore use competition kernels of the form(2.8)and carrying-capacity functions of the form(2.9)Here, the shape parameters *ϵ*_{α} and *ϵ*_{K} measure deviations from the Gaussian case.

### (d) Equilibrium distributions

For *ϵ*_{α}=*ϵ*_{K}=2 (the ‘quartic’ case in which the competition kernel and the carrying capacity are both platykurtic), it can easily be shown numerically that equilibrium distributions of the asexual model (2.2) have multiple peaks whenever *σ*_{α} is small enough.

By contrast, for the sexual model (2.7) with Gaussian ecological functions *α* and *K* with variances and , one can show, by carrying out the various integrals introduced above, that a Gaussian equilibrium distribution exists whose variance satisfies the following equation:(2.10)For example, in the case of random mating, *σ*_{A}=∞, the variance of the Gaussian equilibrium distribution satisfies(2.11)Similarly, in the case of extreme assortative mating, *σ*_{A}=0, there is a Gaussian equilibrium distribution with a variance satisfying(2.12)The existence of Gaussian equilibrium distributions in infinitesimal models in which the ecological functions have Gaussian form may be perceived as supporting the claim that frequency-dependent competition in polymorphic populations does not usually generate multimodal phenotype distributions. However, two important caveats need to be kept in mind. First, even though a Gaussian equilibrium distribution exists, it may not be stable under the dynamics given by equation (2.7). Second, the existence of the Gaussian equilibrium given by equation (2.10) depends on the assumption that the ecological functions *α* and *K* have Gaussian form. As mentioned above, the asexual model with Gaussian ecological functions is structurally unstable, and hence there is no reason to believe that sexual models with non-Gaussian ecological functions and assortative mating would generally admit unimodal equilibrium distributions.

The use of Gaussian functions for *α* and *K* has a long tradition in the literature (Roughgarden 1979). Unfortunately, other than for the fact that a Gaussian decrease in competitive effects and in carrying capacities appears to be heuristically appealing, there is no reason for using these particular functional forms. In fact, Ackermann & Doebeli (2004) have shown that the case in which both the competition kernel and the carrying capacity are Gaussian with finite variance cannot be derived from the underlying mechanistic consumer-resource model introduced by MacArthur (1972), which lies at the basis of most competition models for continuous characters (Roughgarden 1979). This in itself does not mean that the Gaussian case is biologically implausible, but it means that there is no biological reason why this case should receive preferential treatment over other, more general functions, such as those given by equations (2.8) and (2.9). In fact, the mathematical simplicity of the Gaussian case, which sometimes allows analytical equilibrium solutions, may lead to an undesirable bias towards drawing conclusions from a structurally unstable scenario (Meszéna *et al*. 2005). More general models, such as those based on equations (2.8) and (2.9), will generally yield more robust results, even though one typically has to resort to numerical simulations for solving the corresponding dynamical equations for the phenotype distribution.

## 3. Results

Before we turn our attention to the effects of assortative mating on the dynamics of phenotype distributions in sexual populations, we mention two general conditions that are necessary for pattern formation in the form of multimodal distributions. First, the width of the offspring distribution of a given mating pair, *σ*_{f}, must be small enough compared with the width of the carrying-capacity function, *σ*_{K}. Wide offspring distributions tend to homogenize populations and hence to prevent pattern formation. Second, the force of frequency-dependent selection needs to be strong enough for the emergence of multiple phenotypic clusters. For our purposes, this means that in the ecological functions given by equations (2.8) and (2.9), the width of the competition kernel, *σ*_{α}, must be small enough compared with the width of the carrying-capacity function, *σ*_{K}. Wide competition kernels weaken frequency-dependent disruptive selection and hence prevent pattern formation.

### (a) Implications of assortative mating

Even with these necessary conditions being satisfied, we never observed phenotypic pattern formation when mating was random, in which case the equilibrium distributions were invariably unimodal. However, strikingly different outcomes resulted when mating was assortative, i.e. for small enough *σ*_{A}. This is illustrated in figure 1, which shows stable equilibrium distributions of the infinitesimal model for different values of *σ*_{A} for the case in which the competition kernel and the carrying capacity are both Gaussian. As we pointed out in the previous section, this model admits Gaussian equilibrium distributions with variances given by equation (2.10). These equilibrium distributions are stable for high *σ*_{A} (random mating; figure 1*a*) as well as for very low *σ*_{A} (very strong assortment; figure 1*d*). In these cases, the numerical simulations are in exact agreement with the analytical predictions for the variances of the equilibrium distribution given by equation (2.11) for *σ*_{A}=∞ and by equation (2.12) for *σ*_{A}=0.

However, there is a range of intermediate values of *σ*_{A} for which the Gaussian equilibrium distributions are unstable, and instead the dynamics converges to an equilibrium distribution exhibiting distinct phenotypic modes as shown in figure 1*b*,*c*. Since mating is assortative, the phenotypic clusters emerging through such pattern formation represent incipient species: the resultant clusters are reproductively isolated to a large degree, with little gene flow occurring between them. To illustrate the niche partitioning between the incipient species, the grey lines in figure 1 show the carrying-capacity function *K*, indicating the total available niche space. For figure 1*b*,*c*, the initial phenotype distributions were chosen to be very close to the Gaussian equilibrium distribution, but, rather than approaching this Gaussian equilibrium, the system diverges from these unimodal distributions and exhibits pattern formation. Our numerical simulations indicate that when the multimodal equilibrium distributions are stable, they are attractors for a large range of initial conditions. This is illustrated in figure 2 for the case shown in figure 1*b*.

We note that the fact that the Gaussian equilibrium is stable for very small *σ*_{A} (figure 1*d*) is a consequence of the special and non-robust characteristics of the Gaussian case for the asexual model, in which Gaussian ecological functions always generate unimodal solutions (see §2): for very strong assortative mating, the sexual model becomes similar to the asexual model (albeit even in the limit of *σ*_{A}=0, the sexual model is not exactly equivalent to the asexual model unless *σ*_{f}=0).

Figure 3*a*–*d* shows examples of equilibrium distributions for quartic ecological functions, i.e. for *ϵ*_{α}=*ϵ*_{K}=2 in equations (2.8) and (2.9). Again, random mating results in unimodality (figure 3*a*), but assortative mating readily results in multimodal phenotype distributions (figure 3*b*–*d*). In this case, diversification occurs even for very strong assortative mating (figure 3*d*), corresponding to the fact that models with quartic ecological functions admit multimodal solutions even in the asexual case. In contrast to the case of Gaussian ecological functions, the existence of unimodal equilibrium distributions (stable or unstable) cannot be asserted when ecological functions are non-Gaussian. Even if such equilibrium distributions exist in the quartic case, our simulations indicate that they are never stable when assortment is strong enough. In particular, for the values of *σ*_{A} used for figure 3*a*–*d*, the dynamics converges to the shown multimodal equilibrium distributions independent of the various initial conditions that we tested.

In the quartic case, our extensive numerical simulations indicate that the dependence of pattern formation on the various parameters can be roughly summarized as follows. First, for multimodal pattern formation, we have the basic requirement that *σ*_{α} must be small enough to produce frequency-dependent disruptive selection, i.e. *σ*_{α}<*σ*_{K}. Second, both *σ*_{f} and *σ*_{A} need to be small compared with *σ*_{α} and *σ*_{K}. We have found that this can be approximately summarized by the two conditions *σ*_{f}+*σ*_{A}<*σ*_{α} and *σ*_{f}+*σ*_{A}<*σ*_{K}/3. Our simulations indicate that these conditions generally imply pattern formation in the quartic case. These conditions also apply in the case of Gaussian ecological functions, except that with Gaussian functions, we have the additional condition *σ*_{f}<*σ*_{A}. If this condition is not satisfied, the sexual system behaves like the Gaussian asexual model and possesses a stable unimodal distribution (figure 1*d*). On theoretical grounds, it is difficult to assess the biological relevance of the above conditions. There is at least some empirical support for the ecological condition *σ*_{α}<*σ*_{K} (Bolnick *et al*. 2003), and situations in which the genetic kernels (described by *σ*_{f} and *σ*_{A}) are narrower than the ecological kernels (described by *σ*_{α} and *σ*_{K}) do not appear to be unrealistic.

Figure 4 further illustrates the generality of the phenomenon of diversification through pattern formation in phenotype space in the presence of assortative mating. In figure 4*a*, we considered different forms of the carrying-capacity function by varying the shape parameter *ϵ*_{K}, while assuming a Gaussian form for the competition kernel (*ϵ*_{α}=0). For a given carrying-capacity function *K*, we varied the assortative mating parameter *σ*_{A} from values corresponding to random mating (right) to values representing strong assortment (left). For each parameter combination (*σ*_{A}, *ϵ*_{K}), the figure indicates whether the resulting equilibrium phenotype distribution had a single or multiple modes. Analogously, in figure 4*b*, we considered different forms of the competition kernel *α* by varying the shape parameter *ϵ*_{α}, while assuming a Gaussian form for the carrying-capacity function (*ϵ*_{K}=0).

To produce figure 4, we used uniform initial phenotype distributions to start the dynamics for each tested parameter combination. However, the results were virtually identical when Gaussian initial distributions with unit variance were used. The fact that these very different initial conditions yielded the same results underscores that the long-term dynamics of the models considered is largely independent of the initial conditions. Thus, figure 4 shows that diversification resulting in multimodal phenotype distribution occurs for a wide range of assortative mating parameters, and for general classes of competition kernels and carrying-capacity functions.

### (b) Evolution of assortative mating

Given that assortative mating can facilitate phenotypic diversification due to frequency-dependent interactions, as evidenced in figures 1, 3 and 4, it is natural to ask whether there is selection pressure on assortment itself to evolve in initially randomly mating populations. We analyse the selection acting on assortment in two steps. We first assume that the degree of assortment is asexually inherited (one could think of it as being maternally inherited), which permits an adaptive dynamic analysis. We then implement the sexual inheritance of the assortment trait based on standard quantitative genetics in an individual-based model.

For the adaptive dynamic analysis, we extended equation (2.7) to two types differing in their degree of assortment. This allows us to follow the dynamics of the phenotype distributions of the two different types, and in particular to determine when one type can invade the other. With *ϕ*_{1}(*x*) and *ϕ*_{2}(*x*) denoting the phenotype distributions of the two types with assortative mating parameters and , respectively, the resulting dynamics is given by(3.1)(3.2)Since the two types are ecologically equivalent, their per capita death rates *r*.*α**(*ϕ*_{1}+*ϕ*_{2})/*K* are equal, while their birth rates *β*_{1}(*x*) and *β*_{2}(*x*) may differ as a result of differential assortment. These birth rates are derived in Appendix A.

To understand the evolutionary dynamics of assortative mating, we used equations (3.1) and (3.2) to generate pairwise invasibility plots (Metz *et al*. 1996; Geritz *et al*. 1998). These are two-dimensional plots in which possible resident phenotypes are shown on the horizontal axis and possible mutant phenotypes on the vertical axis. For each resident–mutant pair (*σ*_{A,res}, *σ*_{A,mut}), we first let a population consisting only of the resident type reach equilibrium, and then introduced a mutant type at small total density, in order to evaluate whether the mutant's growth rate was positive or negative. The mutant's initial phenotype distribution was assumed to have the same shape as the resident's equilibrium distribution, but with a much reduced total density. Using equations (3.1) and (3.2), the mutant's growth rate was measured as the change in total density over a number of subsequent generations. This procedure generates a partitioning of the pairwise invasibility plot into plus regions, indicating that for such resident–mutant pairs, the mutant can increase when rare and hence will potentially invade the resident, and minus regions, indicating that the mutant cannot invade the corresponding resident, but instead will go extinct.

Figure 5 shows examples of such pairwise invasibility plots that were obtained using the same ecological functions as used in figures 1 and 3. In figure 5, regions in which the mutant can invade the corresponding resident are black, while regions in which the mutant cannot invade the corresponding resident are light grey. In both figure 5*a* (Gaussian ecological functions) and figure 5*b* (quartic ecological functions), the area below the diagonal is black, whereas the area above the diagonal is light grey (note that the diagonal itself belongs neither to the plus nor to the minus region, because a rare mutant with the same assortment phenotype as the resident will neither grow nor decline in total density, since the resident is at equilibrium). For very small resident values of *σ*_{A}, mutant growth rates are very close to zero. This is because in such resident populations, any rare mutant has a strong effective assortment very similar to the resident, which is a consequence of the assumption that the probability of mating between two types is determined by the product of their respective preferences; see equation (A 1) in Appendix A. Thus, for very small values of *σ*_{A}, selection as measured by initial mutant growth rates is nearly neutral, which is indicated by medium grey shading in figure 5*a*,*b*. Nevertheless, the figure shows that there is directional selection for decreased *σ*_{A} and hence for increased assortment. This is not surprising: selection favours increased assortment because assortative mating is a mechanism that facilitates the evolutionary response to frequency-dependent competition (Dieckmann & Doebeli 1999). This mitigation of frequency dependence manifests itself as pattern formation in phenotype space.

There are various ways in which assortative mating could incur fertility costs. One straightforward way to incorporate such costs in the models studied here is to assume that the intrinsic growth rate *r* is negatively affected by increased assortment, i.e. by decreased *σ*_{A}. For example, we can replace the birth terms *β*_{i}(*x*) in equations (3.1) and (3.2) by(3.3)so that the new cost parameter *c* determines the maximal fertility cost, incurred for very strong assortment (i.e. for *σ*_{A}→0). With costs of assortment, the pairwise invasibility plots change qualitatively as shown in figure 5*c*,*d*. For low resident values of *σ*_{A}, the plus and minus regions are now reversed across the diagonal, so that the plus region is above the diagonal and the minus region is below the diagonal. This means that for low resident values of *σ*_{A}, mutants with higher values of *σ*_{A} than the resident, i.e. less-assortative mutants, can invade, while more assortative mutants cannot. Thus, at low values of *σ*_{A}, there is directional selection for less assortative mating. However, at high values of *σ*_{A}, there is still directional selection for increased assortment (i.e. for lower *σ*_{A}). The point at which the two regimes of directional selection meet on the horizontal axis is an evolutionary attractor for the degree of assortment. Once the population has reached the corresponding degree of assortment, either from above or from below, no further invasion of nearby mutants occurs. As expected, costs of assortative mating thus move the evolutionary attractor for the trait *σ*_{A} away from 0. Figure 5*c*,*d* shows that even for moderately high costs of assortative mating, the degree of assortment is still expected to evolve to substantial levels.

Finally, we used an individual-based model to investigate the full evolutionary dynamics of assortment. In such a model, individuals are described by their ecological trait *x* and their assortment trait *σ*_{A}. At each point in time, every individual experiences a per capita death rate and a per capita birth rate. The per capita death rate is determined by the ecological trait and is calculated according to the death term in equation (2.7) (integrals are replaced by sums over all individuals in the population). The per capita birth rate incorporates potential costs of assortment and is given by equation (3.3). At each point in time, individual rates are summed up to give the total birth and death rates *B* and *D*, respectively. The waiting time until the next birth or death event is drawn from an exponential probability distribution with mean 1/(*B*+*D*), and a birth or death event is then chosen with probabilities *B*/(*B*+*D*) and *D*/(*B*+*D*), respectively. If a death event occurs, one individual is chosen probabilistically according to its relative contribution to the total death rate. The chosen individual is removed, and the birth and death rates of all other individuals are adjusted accordingly. If a birth event occurs, one individual is chosen probabilistically according to its relative contribution to the total birth rate. The chosen individual then selects a mating partner probabilistically according to the mate choice function given by equation (A 1) in Appendix A, evaluated for all other individuals in the population (as before, mate choice is based on the ecological trait). The resulting mating pair produces an offspring whose phenotypes are drawn from two Gaussian distributions with means given by the midparent values of the two traits and with standard deviations *σ*_{f} for the ecological trait and *σ*_{f,ass} for the assortment trait. The offspring individual is inserted, and the birth and death rates of all other individuals are adjusted accordingly. This stochastic model naturally extends to finite populations the deterministic models introduced and analysed above.

Figure 6 shows examples of the joint evolutionary dynamics of the ecological phenotype and the assortment phenotype in the individual-based model. The initial conditions for these dynamics were chosen such that populations were mating approximately randomly. As a consequence, the phenotype distribution for the ecological trait was initially unimodal (figure 6*a*). However, despite costs of assortment, assortative mating readily evolved to a degree that allowed the formation of phenotypic clusters and hence diversification (figure 6*b*,*c*).

## 4. Discussion

Our results show that even in infinitesimal genetic models for the dynamics of continuous phenotype distributions in sexual populations, frequency-dependent selection can split the population into separate phenotypic clusters when mating is assortative. If such pattern formation in phenotype space occurs, the emerging phenotypic clusters represent incipient species, as they are at least partially reproductively isolated due to assortative mating. Thus, frequency-dependent selection can cause adaptive speciation in these models.

Apart from assortative mating, two requirements need to be met for diversification to occur. The width of the offspring distribution produced by a given mating pair must be small enough, and frequency dependence must be strong enough. Our extensive numerical explorations of thousands of different cases revealed that when these conditions are satisfied, pattern formation occurs for a wide range of assortative mating parameters and for a wide range of forms of the competition kernel and the carrying-capacity function. Moreover, assortative mating can readily evolve in initially randomly mating populations even if it comes at considerable cost. This is in accordance with the recent results from explicit multilocus models showing that costs to assortment do not prevent adaptive sympatric speciation unless such costs are high (Doebeli 2005; Doebeli & Dieckmann 2005; Bürger & Schneider 2006; Schneider & Bürger 2006; Bürger *et al*. in press). In this study, we have focused on assortative mating based on the ecological trait under frequency-dependent selection. Such assortment models are usually called one-allele models (Kirkpatrick & Ravigné 2002), in contrast with two-allele models, in which assortment is based on a selectively neutral display trait. The evolution of assortment in two-allele models for adaptive speciation has been studied in models with explicit multilocus genetics (Dieckmann & Doebeli 1999; Doebeli 2005), with the conclusion that while recombination between the display and the ecological traits hinders adaptive diversification, speciation is nevertheless possible in such scenarios, even if there are costs to assortment (Doebeli 2005). It would clearly be interesting to study two-allele scenarios in infinitesimal models incorporating frequency-dependent competition.

Additionally, in this paper, we have considered only the costs to assortment that differentiate between different levels of assortative mating. For a population with a fixed degree of assortment, these costs are the same for all individuals. However, it is important to consider also scenarios in which there is a cost to assortment due to Allee effects. In this case, individuals may differ in fertility even in populations with a fixed degree of assortment, because rare phenotypes will encounter fewer preferred mates than common phenotypes, and so may have lower fertility. Allee effects can be incorporated into infinitesimal models (Noest 1997), and results for the dynamics of pattern formation in such models will be reported elsewhere. Implementing the individual-based model introduced at the end of the previous section is a straightforward exercise, and we invite readers to explore the dynamics of phenotypic pattern formation based on their own models and/or implementations. Whether the regions in parameter space for which diversification can be observed are biologically realistic is a question that needs to be addressed in empirical studies, but mathematically speaking, it is clear that phenotypic pattern formation is a robust outcome of infinitesimal models.

In fact, diversification involving assortative mating may be easier in infinitesimal models than in other more explicit genetic models based on a finite number of loci with finite effects, such as those investigated by Dieckmann & Doebeli (1999), Kirkpatrick & Nuismer (2004), Schneider & Bürger (2006), Bürger & Schneider (2006) and Bürger *et al*. (in press). One obstacle to speciation in such models is that genetic variation in the ecological trait can be exhausted if mating is strongly assortative (Kirkpatrick & Nuismer 2004; Bürger *et al*. in press). This cannot happen in deterministic infinitesimal models, in which offspring distributions always range across the whole spectrum of phenotypes (albeit with very low frequencies at most phenotypes). In infinitesimal models, all phenotypes are thus present at all times, and hence any loss of phenotypes on which selection can act is not a problem. At any rate, the results reported here for infinitesimal models are in surprisingly good overall agreement with models for adaptive diversification based on adaptive dynamics and multilocus genetics (Dieckmann & Doebeli 1999; Doebeli & Dieckmann 2000; Bürger *et al*. in press), supporting the understanding that adaptive speciation due to frequency-dependent interactions can safely be considered a theoretically plausible scenario.

It should, of course, be analysed whether infinitesimal models can indeed provide a sufficiently accurate approximation of multilocus dynamics under frequency-dependent disruptive selection. This may well be the case over shorter time spans, involving mostly standing genetic variation (Bulmer 1980), but the robustness of this approximation becomes more uncertain when one takes into account mutation and substitution of allelic effects at the loci. It could happen that variation typically becomes concentrated on just one or a few loci (van Doorn & Dieckmann in press), which is similar to the outcome predicted in adaptive dynamics models. Alternatively, variation might increase at all loci, potentially increasing *σ*_{f}, which could in principle result in a stable unimodal equilibrium distribution. In addition, there are other possible evolutionary responses to frequency-dependent competition, such as sexual dimorphism and a widening of individual niche widths (Bolnick & Doebeli 2003; Ackermann & Doebeli 2004; Rueffler *et al*. 2006) that could be considered. Although these are relevant issues, here we have focused our treatment on the infinitesimal model used by Polechová & Barton (2005) because, at the very least, it represents a conceptually interesting and traditionally well-received case.

Our results are in contrast to those reported by Polechová & Barton (2005) for infinitesimal models in which both the competition kernel and the carrying capacity are Gaussian functions. Rather than focusing on the actual dynamics of phenotype distributions, results presented by Polechová & Barton (2005) only concern the variance of Gaussian equilibrium distributions: these authors seem to have implicitly assumed that the dynamics of the infinitesimal model always converges to such Gaussian solutions. In particular, Polechová & Barton (2005) did neither consider the possibility that the Gaussian equilibrium could be unstable nor investigate models with non-Gaussian ecological functions, which might not admit unimodal equilibrium solutions in the first place. Our results show that, in general, the dynamics of infinitesimal models does not converge towards unimodal equilibrium distributions when mating is assortative. While they do not appear to have numerically solved the infinitesimal model, Polechová & Barton (2005) mention that their simulations of a related model, the ‘symmetric’ model with explicit multilocus genetics, suggest that dynamics in that symmetric model always converges to Gaussian equilibrium solutions, thus apparently lending support to their implicit assumption of stable Gaussian equilibrium distributions for the infinitesimal model. However, evaluating dynamical stability of one model in terms of another model is obviously not possible. Moreover, it has already been shown by Doebeli (1996) that the symmetric model also often exhibits pattern formation in the form of bimodal equilibrium distributions when mating is assortative.

Based on their analysis of the infinitesimal model, Polechová & Barton (2005) concluded that the process of assortment itself, irrespective of any frequency-dependent competition, was the most important driver of divergence in sexual models of sympatric speciation. This conclusion was based on the observation that in the infinitesimal model with *σ*_{K}=∞ and sufficiently strong assortment, the variance of a solution can increase without bound, even in the absence of frequency-dependent competition (i.e. if *σ*_{α}=∞). The possible role of assortment in permitting genetic divergence is of course a relevant issue. For the infinitesimal model, equilibrium solutions with infinite genetic variance exist. For instance, for an infinite width of the carrying-capacity function, and with very strong assortative mating (*σ*_{A}=0) but without frequency dependence, equation (2.12) for determining the equilibrium variance *σ*_{eq} reduces to , which only admits *σ*_{eq}=∞ as a solution. This solution might be regarded as an artefact of the assumption of the infinitesimal model that there is an unlimited supply of genetic variation in the population. Nevertheless, the qualitative conclusion that assortment sometimes can relax genetic constraints and thus enable an increase in genetic variation seems valid.

An interesting question is then whether assortment itself, without frequency-dependent competition, can lead to pattern formation. For our model, making the assumptions of very strong assortment (*σ*_{A}=0), no frequency dependence (*σ*_{α}=*∞*) and a Gaussian carrying capacity with finite *σ*_{K}, we infer from equation (2.12) that there is a Gaussian equilibrium distribution with finite variance , which is approximately equal to *σ*_{f}*σ*_{K} for small *σ*_{f}. Numerical simulations indicate that this unimodal solution is always stable. Similarly, when the carrying capacity is not of Gaussian form, but still unimodal, even very strong assortative mating never generated pattern formation in the absence of frequency dependence. Thus, in our numerical analysis of the infinitesimal model, speciation was never observed when frequency dependence was absent, independent of the strength of assortment.

We can gain an intuitive understanding of the reason for this conclusion by noting that equation (2.7) in the limit of *σ*_{A}→0 approaches the asexual case in equation (2.2), except that the offspring distribution acts like a mutation kernel, smoothing the distribution *ϕ*. For such an asexual model, with very high rate of mutation and no frequency-dependent competition, there is no reason to expect clustering of phenotypes, since there are no forces that could counteract the homogenization of a multimodally clustered distribution. Thus, under very strong assortment, the sexual production of offspring, involving segregation and recombination, can increase genetic variation in the infinitesimal model, in a manner analogous to the process of mutation. This source of variation, however, cannot by itself drive pattern formation, just as mutation cannot by itself drive pattern formation. In the infinitesimal models considered here and in Polechová & Barton (2005), speciation is thus impossible without frequency-dependent competition.

The results reported here are in complete agreement with those of Noest (1997), who presented an analytical study of a special class of infinitesimal models, in which the carrying capacity was assumed to be uniform (i.e. independent of *x*), and the competition kernel was assumed to be Gaussian. These models admit a uniform equilibrium phenotype distribution, and Noest (1997) investigated the conditions under which this uniform solution becomes unstable in the presence of assortative mating. His analytical results match our numerical results in essential aspects: if the offspring distribution is sufficiently narrow and frequency dependence is sufficiently strong, then the uniform solution can become unstable when mating is assortative. Noest (1997) did not study more general forms of the ecological functions, for which analytical results are not feasible, and he did not consider the evolution of assortative mating. However, his results already clearly showed that frequency dependence and assortative mating can break up the sexual continuum (Maynard Smith & Szathmáry 1995) through the formation of multimodal distributions in phenotype space. Our results lead to the same conclusion for a more general class of models and evolutionary scenarios: adaptive speciation can occur as a result of pattern formation in phenotype space due to frequency-dependent selection and the evolution of assortative mating.

## Acknowledgments

M.D. was supported by NSERC (Canada) and the James S. McDonnell Foundation (USA). O.L. was supported by the Swedish Research Council. U.D. acknowledges financial support by the Vienna Science and Technology Fund, WWTF.

## Footnotes

- Received August 12, 2006.
- Accepted September 3, 2006.

- © 2006 The Royal Society

## References

## Notice of correction

Figures 2 and 6 are now presented in their correct form.

16 November 2006