The role of demographic stochasticity in a speciation model with sexual reproduction
aa r X i v : . [ q - b i o . P E ] M a r The role of demographic stochasticity in a speciation model with sexual reproduction
Luis F. Lafuerza and Alan J. McKane
Theoretical Physics Division, School of Physics and Astronomy,The University of Manchester, Manchester M13 9PL, UK (Dated: October 2, 2018)Recent theoretical studies have shown that demographic stochasticity can greatly increase thetendency of asexually reproducing phenotypically diverse organisms to spontaneously evolve intolocalised clusters, suggesting a simple mechanism for sympatric speciation. Here we study the roleof demographic stochasticity in a model of competing organisms subject to assortative mating. Wefind that in models with sexual reproduction, noise can also lead to the formation of phenotypicclusters in parameter ranges where deterministic models would lead to a homogeneous distribution.In some cases, noise can have a sizeable effect, rendering the deterministic modelling insufficient tounderstand the phenotypic distribution.
PACS numbers: 05.40.-a, 87.23.Cc, 87.10.Ca
I. INTRODUCTION
Establishing the determinants of biological diversity isa fundamental question in biology. A particular aspect ofthis question that has attracted a great deal of attentionis the distribution (in genotype or phenotype space) of apopulation of interacting organisms, and to what extentand under what conditions clusters of similar individualstend to arise [1–3]. A better understanding of this ques-tion could shed some light onto the process of sympatricspeciation, whereby a ‘mother’ species splits into two ormore other species without geographic isolation [4, 5].There is a history of mathematical models proposedto elucidate the mechanisms whereby species (as clustersof phenotypically similar organisms) tend to form spon-taneously due to competition [6–10]. Models initiallyconsidered asexual populations, to facilitate analyticaltractability. It was shown that spontaneous clusteringcan occur [7, 9, 11], but the phenomena is somewhatsensitive to particular assumptions about the functionalform of the interaction kernels used [12, 13], question-ing the biological relevance of the findings. Models inwhich reproduction is sexual tend to lead to spontaneousclustering with less restrictive assumptions [8, 14], butthe biological realism of the conditions required has alsobeen questioned [5]. Recently, it was shown that stochas-tic effects can greatly increase the range of parametersfor which species are formed in asexual models [15, 16],presenting stochastic pattern formation as a novel mech-anism for speciation, which has particular biological im-plications [17]. In this paper we will study the effects ofstochasticity in a speciation model with sexual reproduc-tion. We will show that demographic stochasticity canincrease the parameter range in which species cluster-ing is observed and that, in some cases, noise can havea sizeable effect, rendering the deterministic modellinginsufficient to understand the phenotype distribution.The rest of the paper is organised as follows. We willpresent the model as well as its mathematical formula-tion in Sec. II. In Sec. III we will present the analysis inthe deterministic (large population) limit, summarising some previously known results and presenting some newones. We will then perform the analysis of the stochasticeffects in Sec. IV, highlighting the cases in which noiseeffects are largest. We will conclude with a summaryand conclusions. Some technical details are left for twoappendixes.
II. DESCRIPTION OF THE MODEL
The model consists of a population of individuals whichreproduce sexually and die through competition witheach other. Each individual is described by its pheno-type, which determines the extent to which the organ-ism competes with other organisms, and the likelihoodwith which it can mate with a given other organism. Forsimplicity, we will assume that the phenotype is well de-scribed by a single scalar variable which can take on allpossible real values. A simple example in which a singlescalar variable is the main determinant of the strengthof competition between organisms could be the beak sizein birds or the jaw size in lizards [18] (which determinesthe extent to which they feed on the same resources), orin general, body size. The model is stochastic in thatdeaths and births are random events (which, however,take place with probabilities determined by the state ofthe system). There are two basic processes:(i)
Death . The death rate (probability per unittime) of individual i , d i , is given by d i =[ Kf ( x i )] − P j g ( x i − x j ), where g ( x ) is the com-petition kernel that quantifies the strength withwhich two individuals with phenotype distance x compete, K is a constant that controls the overallcarrying capacity of the ecosystem and f ( x ) is afunction that determines the relative intrinsic ad-vantage of phenotype x .(ii) Reproduction . Each individual reproduces at a rateone. The probability that individual i mates withindividual j is proportional to m ( x i − x j ), with m a function determining the strength of assortment(preference for individuals which are alike). If indi-viduals i and j mate, an offspring is generated withphenotype given by x = ( x i + x j ) / ζ , where ζ is a random variable with probability density func-tion r ( x ). This models mutation about the averageof the parents phenotype ( x i + x j ) / ζ ,is a Gaussian random variable with zero mean and vari-ance given by σ m , that, for simplicity, will be assumed tobe phenotype-independent (the influence of this quantitywill be one of the main aspects of our investigation). Thefunction f ( x ) modulates the relative advantage of thephenotypes; it typically decays to zero for large | x | , con-straining the phenotypes to a given region of interest. Weare also assuming that the organisms are hermaphroditic,but this assumption will be relaxed later. The model isa generalisation of the asexual model studied in [15].Following [15], we will describe the state of the systemby the density in phenotype space, Kφ ( x ), where φ ( x ) = 1 K N ( t ) X i =1 δ ( x − x i ) . (1)Here N ( t ) is the number of organisms at time t and δ ( x )is the Dirac delta function. We have introduced the car-rying capacity, K , to obtain a function φ ( x ) that has awell-defined K → ∞ limit.When an organism with phenotype y dies, φ ( x ) is mod-ified by the subtraction of delta function centred at y .Similarly, if a new organism with phenotype y is born, φ ( x ) is modified by the addition of a delta function at y .With this in mind, we define the operators ∆ ± y by theiraction on a generic functional F [ φ ( x )] as:∆ ± y F [ φ ( x )] = F [ φ ( x ) ± K δ ( x − y )] . (2)Now suppose that γ ( x, φ ) is the density rate at which anindividual with phenotype x dies, so that γ ( x, φ ) dxdt isthe probability that an individual with phenotype in theinterval ( x, x + dx ) dies in the time interval ( t, t + dt ),given that the state of the system is given by φ at time t .The definition of the process implies that γ ( φ, t ) is givenby: γ ( x, φ ) = Kφ ( x ) 1 Kf ( x ) N ( t ) X i =1 g ( x − x i )= K φ ( x ) f ( x ) Z φ ( y ) g ( x − y ) dy ≡ K φ ( x ) f ( x ) φ ∗ g ( x ) . (3)The probability that individual i , if it mates, does sowith individual j is m ( x i − x j ) / P k m ( x i − x k ) = m ( x i − x j ) /K R φ ( y ) m ( x i − y ) dy = m ( x i − x j ) /Kφ ∗ m ( x i ). Theprobability density that their offspring has phenotype x is given by r ( x − ( x i + x j ) / r ( x ) is the Gaussianprobability density with variance σ m . With this in mind,we see that the rate density at which a new individual iscreated with phenotype x , β ( x, φ ), is: β ( x, φ ) = N ( t ) X i,j =1 r ( x − ( x i + x j ) / m ( x i − x j ) Kφ ∗ m ( x i ) (4)= K Z φ ( y ) Z φ ( z ) m ( y − z ) φ ∗ m ( y ) r ( x − ( y + z ) / dydz. Combining the two contributions to the change in φ , theprobability density of finding the system at state φ attime t , P ( φ, t ), changes in time according to the followingfunctional master equation [15]: ∂∂t P ( φ, t ) = Z [(∆ − x − β ( φ, x ) P ( φ, t )+ (∆ + x − γ ( φ, x ) P ( φ, t )] dx. (5)Due to the non-linearity of the system (that arises dueto the interactions) we are unable to obtain an exact so-lution and some approximation scheme is needed to pro-ceed. Expanding the ∆ ± x operators in Eq. (5) to secondorder in K − , we can derive a functional Fokker-Planckequation, given in Appendix A. For our purposes, it isclearer to work with the equivalent stochastic differentialequation which takes the form (see Eqs. (A4) and (A5)) ∂φ ( x, t ) ∂t = − Z φ ( x, t ) φ ( y, t ) g ( x − y ) f ( x ) dy + Z φ ( y, t ) φ ( z, t ) φ ∗ m ( y ) m ( y − z ) r ( x − ( y + z ) / dy dz + η ( x, t ) √ K , (6)where η ( x, t ) is a Gaussian white noise with zero meanand with a correlator which is given by Eq. (A5) of theAppendix. It is interesting to compare this equation withthe analogous equation in the asexual case. There m ( x − y ) = δ ( x − y ) and so φ ∗ m ( x ) = φ ( x, t ). Then Eq. (6)becomes ∂φ ( x, t ) ∂t = − Z φ ( x, t ) φ ( y, t ) g ( x − y ) f ( x ) dy + Z φ ( y, t ) r ( x − y ) dy + η ( x, t ) √ K , (7)which is the equation found in Ref. [15], apart from thefunction f ( x ), which was not included in the form ofthe model previously analysed (note that we are allowingself-fertilization since, for simplicity, in Eq. (4) we donot exclude i = j ; forbidding the i = j case would add O (1 /K ) terms). We will now analyse Eq. (6), first ofall in the deterministic limit, and then in the generalstochastic setting. III. DETERMINISTIC ANALYSIS
The deterministic limit corresponds to taking K → ∞ ,and so the governing equation is simply Eq. (6), butwith the last (noise) term absent. Some progress maybe made analytically if we assume that the ecologicalfunctions, namely the competition kernel [ g ( x )], the mat-ing function [ m ( x )], the function modulating the carry-ing capacity [ f ( x )] and the offspring distribution [ r ( x )],are all Gaussian functions. We denote their variancesby σ c , σ a , σ f and σ m , respectively. In this case, the de-terministic equation has a Gaussian stationary solution, φ ( x ) st = Ce − x / (2 σ ) , with σ and C both satisfyingcomplicated algebraic equations; the equation for σ be-ing first derived by Doebeli et. al. [14]. These authorsalso found that numerical integration of the determin-istic equation showed that when σ c and σ m are smallcompared with σ f , there is an intermediate range of σ a , σ m . σ a . σ f / − σ m , for which the Gaussian solutionbecomes unstable and a multi-modal stationary solutionis obtained.The random mating case ( σ a → ∞ ) is particularlyinteresting to analyse. In this limit, the equation for σ reduces to a cubic equation σ + α σ + α σ + α = 0,where α = 2 σ m + σ c , α = σ f σ c − σ m (2 σ f − σ c ) ,α = − σ m σ c σ f . (8)Since α > α <
0, the cubic equation alwayshas a single positive solution. This then is the requiredsolution. In this limit, the constant C is found to be C = s σ + σ c σ + 2 σ m , (9)(we have not taken f ( x ) to be normalised; f ( x ) =exp( − x / (2 σ f )), so that σ f controls the size of pheno-type space available). If we now, in addition, investigatethe σ f → ∞ limit (the most relevant when the phenotypedistribution is not too constrained by the external fitnesslandscape), we find two very different regimes, dependingon the relative width of the competition kernel and thereproduction noise distribution. When the reproductionnoise, σ m , is smaller than a critical value given by σ c / σ f , σ = 2 σ m σ c / ( σ c − σ m ). If σ dependson σ f in a way which means that it diverges as σ f → ∞ ,then the term σ dominates over the term quartic in σ st and the term quadratic in σ st dominates over α . There-fore, σ ≃ q σ f (4 σ m − σ c ), and this implies that if thereproduction noise exceeds the critical value, σ m > σ c / σ f . In both cases, numerical integration of thedeterministic equation shows that the Gaussian distribu-tion is always stable and the phenotype distribution isalways uni-modal, a consequence of the random mating. The factors determining the transition into a multi-modal distribution can be understood in a simpler way ifwe assume that the range of phenotype space is finite, forinstance − π < x ≤ π , and assume that the deterministicequation (Eq. (6) without the noise term) satisfies peri-odic boundary conditions. The competition function, g ,the assortment function, m , and the offspring distribu-tion will be assumed to be periodic functions of period2 π . Since there is now no need for the function f ( x )to regulate the competition process at large | x | , we set f ( x ) = 1. While the periodic boundary condition as-sumption is biologically unrealistic, we expect it to havea small impact when the scales of the competition, mat-ing and offspring distributions are all much smaller thanthe available region of phenotype space (as determinedby f ).Under these conditions, the deterministic equation hasa uniform solution φ ( x, t ) = constant. Using the normal-isation R π − π g ( x ) dx = 1, one finds that φ st = 1. A linearstability analysis of the deterministic equation aroundthe solution φ = 1 is performed in Appendix A. Thisshows that the uniform solution is unstable if:2 r k m − k/ − r k m k m − k/ − − g k > , (10)for some k , with r k , m k and g k the Fourier modes ofthe reproduction noise, the mating and the competitionkernel, respectively. An equivalent expression to (10) ina slightly different model was derived in [8].From the inequality (10) one can obtain the bifurcationdiagram of the system that shows the regions of parame-ter space in which the deterministic equation [Eq. (A6)]leads to patterns. Figure 1 portrays several bifurcationdiagrams for the case in which the ecological functionsare uniform and Gaussian functions (projected onto the( σ a , σ m ) plane). The figure shows that phenotypic clus-ters are observed for low values of the reproduction noise σ m and for intermediate values of the assortativity scale σ a , in accordance with what was found in [14]. Intu-itively, clusters tend to form with a distance, d , betweenthem so that they barely compete (so d ' σ c ), but theseclusters will only be durable if individuals from differentclusters do not mate (so d > σ a ), suggesting that clus-ters will not form if σ a ≫ σ c . Sexual reproduction tendsto concentrate the phenotype distribution, so a moder-ate value of σ a can promote the formation of clusters.Moreover, the clusters can be well-defined only if the re-production noise is not too large ( σ m < d ).In the asexual case [15, 16] it was found that the rangeof parameters for which patterns occur was much greaterthan that predicted by a stability analysis of the kindcarried out above. So motivated by the expectation thatthe same may be true in the case of sexual reproduction,we go on to analyse our model for finite values of K ,when stochastic fluctuations will be present. . . . . . . s a2 s m s c2 = s c2 = s c2 = s c2 = . . . . s a2 s m s c2 = s c2 = s c2 = s c2 = FIG. 1: Region of stability of the homogeneous solution ofthe deterministic equation in σ a − σ m space (mating range -mutation range), for several values of σ c (competition range),for uniform (upper panel) and Gaussian (lower panel) formsof the ecological functions. The homogeneous solution is un-stable below the lines, leading to the appearance of patterns.The kinks are produced when the maximum value of the left-hand side of the inequality (10) changes from one value of k to another. The black dot marks parameter values used inFig. 2. IV. STOCHASTIC EFFECTSA. Weak noise effects
Numerical simulations of the individual-based modelshow that, for moderate values of the carrying capacity, K , the observed phenotype distribution does not alwaysagree with the results obtained from the deterministicequation (A6). For example, at points in parameter spacecorresponding to a stable uniform solution, but not toofar from the instability boundary, one can observe clus-ters in phenotype space forming for small values of K (see Fig. 2).The origin of these stochastic patterns can be under-stood by analysing the stochastic differential equation t x 0.51.01.52.0 t x FIG. 2: Time evolution of the phenotype distribution for σ c =0 . , σ a = 0 . , σ m = 0 .
043 and uniform ecological functions,in the region of deterministic stability of the homogeneoussolution (marked as a black dot in Fig. 1). The carryingcapacity is K = 20 (upper panel) and K = 200 (lower panel).Note that for the smaller carrying capacity clear phenotypicclusters form. (6). Specifically, we apply the van Kampen system-size expansion [20] by expanding φ about the homoge-neous solution found in Sec. III and writing φ ( x, t ) =1 + K − / ξ ( x, t ). The factor K − / reflects the nature ofthe fluctuations at large K , and ξ ( x, t ) is a new stochas-tic field. The bulk of the calculation is exactly the sameas that carried out when performing the linear stabil-ity analysis in the deterministic case, except in this casethe existence of the K − / factor ensures that the noiseterm in Eq. (6) is retained. Thus, going over to Fouriervariables and using Eq. (A8), one finds that ddt ξ k ( t ) = (cid:2) r k m − k/ − − g k − r k m k m − k/ (cid:3) ξ k ( t )+ η k ( t ) , (11)where η k ( t ) is the Fourier transform of η ( x, t ). The corre-lation function of this noise is only calculated to leadingorder within the linear noise approximation, that is, set-ting φ ( x, t ) = 1. This gives, using Eqs. (A2), (A3) and(A5), R B ( φ, x, y ) dy = 2, and therefore h η k ( t ) η − k ( t ′ ) i = 2 (2 π ) δ ( t − t ′ ) , (12)showing that, as usual, the noise is additive within thelinear noise approximation.The size of the stochastic patterns can be quantifiedby looking at the spatial covariance of the phenotypedistribution in the stationary state:Cov( φ ( x ) , φ ( x + ∆)) ≡h ( φ ( x ) − h φ ( x ) i ) ( φ ( x + ∆) − h φ ( x + ∆) i ) i . (13)If the distribution in phenotype space shows high-densityregions separated by low-density ones with a well-definedaverage distance, the spatial covariance will display a siz-able spatial modulation in ∆. In the linear noise regime,the covariance becomes Z Cov( φ ( x ) , φ ( x + ∆)) dx = 1 K Z h ξ ( x ) ξ ( x + ∆) i dx = 12 πK X n h ξ n ξ − n i e in ∆ . (14)Using Eq. (A13), and assuming that g n , m n and r n areall even functions of n (which follows if the ecologicalfunctions are symmetric), we may find h ξ n ξ − n i , and thenfinally obtain: Z Cov( φ ( x ) , φ ( x + ∆)) dx = 1 K ∞ X n = −∞ e in ∆ (cid:2) g n + r n m n m n/ − r n m n/ (cid:3) . (15)Expression (15) is compared against numerical simula-tions in Fig. 3 for the case in which the ecological func-tions are uniform distributions. For low values of thecarrying capacity, K , the covariance shows a clear spatialmodulation, revealing the presence of stochastic patterns.The theoretical expression agrees qualitatively, and thequantitative agreement becomes better as K increases,where the linear noise assumption is expected to be abetter approximation to the stochastic dynamics. The re-sults are very similar when the ecological functions haveother forms, for instance if they are Gaussians or belongto the family exp( −| x | l / (2 σ l )) with varying l .In summary, the linear noise approximation shows thatstochastic pattern formation can lead to the formationof clusters in phenotype space in situations when the de-terministic description predicts a uniform distribution.Interestingly, the stochastic patterns (as opposed to thedeterministic ones) decrease with the carrying capacity(that controls the abundance of individuals), which leadsto particular biological implications [17] when this mech-anism is dominant.When comparing with the asexual case [15], thestochastic patterns now appear on a relatively narrow − − D C o v ( f ( x + D ) , f ( x )) K=20K=50K=200
FIG. 3: Spatial covariance of the phenotype density distribu-tion. The ecological functions are uniform distributions withvariances given by σ c = 0 . σ a = 0 . σ m = 0 .
043 (markedas a black dot in Fig. 1). Numerical results (symbols) wereaveraged over at least 100 measurements, after a transient of t = 2000. The lines correspond to the theoretical expression,Eq. (15). zone close to the deterministic transition line. We willnow show that stochastic effects can also play a moreprominent role in the sexual case. B. Strong noise effects
A stronger noise-induced effect occurs when the mat-ing is random. As discussed in Sec. III, in the randommating case there are two regimes: one leading to a nar-row phenotype distribution, for low reproduction noise(4 σ m < σ c ), and one leading to a broad phenotype distri-bution in the large reproduction noise case. In the broadphenotype regime (4 σ m > σ c ), simulations show that thephenotype distribution observed is much narrower thatthe one predicted by the deterministic analysis, see Fig. 4.This effect happens rather generically in this regime.In the random mating case, since the mating range isthe largest scale of the system, assuming periodic bound-ary conditions gives results that are quite different fromthose of the system with open boundaries. In this lastcase there is a strong bias towards the centre of the phe-notypic space. For this reason we focus on the openboundary conditions case, which is the more biologicallyrelevant. This, however, greatly complicates the math-ematical analysis of the noise effects. Simulations show(Fig. 4) that stochastic effects lead to the formation ofa phenotypic cluster with a width that approaches thedeterministic prediction only for very large values of thecarrying capacity, K . We, therefore, find that in this −10 −5 0 5 10 . . . x f ( x , t ) K=100K=400K=1000Theory t x FIG. 4: Upper panel shows the observed instantaneous pheno-type distribution in the steady state for different values of thecarrying capacity, together with the deterministic prediction.Lower panel shows the evolution of the phenotype distribu-tion for K = 400, starting with the stationary distributionpredicted by the deterministic analysis. Parameter values are σ c = 0 . , σ m = 0 . , σ f = 6400. regime demographic stochasticity has a quite large ef-fect, again leading to the formation of tight phenotypicclusters when the deterministic analysis predicts a broaddistribution. Similar results are obtained when the eco-logical functions do not have a Gaussian form, showingthe robustness of this phenomenon. V. THE CASE WHEN INDIVIDUALS BELONGTO ONE OF TWO SEXES
So far we have assumed, for simplicity, that any two or-ganisms can mate, i.e. the organisms are hermaphroditic.We now go on to model the situation with two explicitlydifferent sexes.We will denote the phenotype of female organism by x i and that of the male organism by y α . We use differentindices, since the number of male and female organismsat a given time will typically be different, and so the range of these indices will be different. The model is asin Sec. II, but with the following modifications:(a) Only females reproduce and at rate one.(b) The probability that female organism i mates withmale organism α is proportional to m ( x i − y α ).(c) The probability that the offspring is male or femaleis 1 / i is d i = P j g ( x i − x j ) + P α g ( x i − y α ).With these assumptions, one can follow through theanalysis given in Sec. III and find the analogues of Eq. (6): ∂φ f ∂t = − Z φ f ( x, t ) [ φ f ( y, t ) + φ m ( y, t )] g ( x − y ) f ( x ) dy + 12 Z φ f ( y, t ) φ m ( z, t ) φ m ∗ m ( y ) m ( y − z ) r ( x − ( y + z ) / dy dz + η f ( x, t ) √ K , (16)and ∂φ m ∂t = − Z φ m ( x, t ) [ φ f ( y, t ) + φ m ( y, t )] g ( x − y ) f ( x ) dy + 12 Z φ f ( y, t ) φ m ( z, t ) φ m ∗ m ( y ) m ( y − z ) r ( x − ( y + z ) / dy dz + η m ( x, t ) √ K , (17)where the subscripts f and m denote female and malerespectively. It is convenient to work with the sum anddifferences of φ f and φ m : S ( x, t ) = φ f ( x, t ) + φ m ( x, t )and D ( x, t ) = φ f ( x, t ) − φ m ( x, t ). If we additionally take K → ∞ to obtain the deterministic equations, we findthat ∂S ( x, t ) ∂t = − Z S ( x, t ) S ( y, t ) g ( x − y ) f ( x ) dy + 12 Z [ S ( y, t ) + D ( y, t )][ S ( z, t ) − D ( z, t )][ S − D ] ∗ m ( y ) × [ m ( y − z ) r ( x − ( y + z ) / dy dz ] , (18)and ∂D ( x, t ) ∂t = − Z D ( x, t ) S ( y, t ) g ( x − y ) f ( x ) dy. (19)From Eq. (19) we see that a steady state solution is D ( x ) = 0, that is, φ f ( x ) = φ m ( x ). Then the de-terministic equations for φ f and φ m collapse into eachother, and agree with the deterministic equation in thehermaphroditic case, apart from the factor of 1 / φ ( x ) st = Ce − x / (2 σ ) , with σ and C satisfying the same complicated algebraic equations,except that C takes on a value of one quarter the valuefound in Sec. III. Note that the total population is nowhalf that obtained in the hermaphroditic case becausenow we assume that only females can initiate reproduc-tion events.To examine the instability leading to the appearanceof patterns, we can again assume a finite interval − π 12 ˜ D ( x, t ) ⇒ ˜ D ( x, t ) = ˜ D ( x, e − t/ , (20)which shows that the solution D ( x, t ) = 0 is always sta-ble. For ˜ S ( x, t ), it is more convenient to work in Fourierspace (see Appendix A). We find that d ˜ S k ( t ) dt = 12 (cid:2) r k m − k/ − r k m k m − k/ − − g k (cid:3) ˜ S k ( t )+ F k ˜ D k ( t ) , (21)where F k is a function of k only. Equation (20) showsthat ˜ D k ( t ) decays exponentially with t , so Eq. (21) im-plies that one obtains the same stability condition as theone found for φ in the hermaphroditic case, Eq. (10).We can therefore conclude that the stability boundariesin the case of two sexes are identical to those found inthe hermaphroditic case. This is confirmed by numericalsimulations of the stochastic version of the model. Wealso find that stochastic pattern formation takes place asin the hermaphroditic case, with results from the case oftwo sexes being equivalent to those of the hermaphroditiccase, but with a factor of 1 / VI. CONCLUSION The precise definition of exactly what constitutes aspecies is still open to debate [21, 22]. One of several dif-ferent alternatives is the ‘phenotypic clustering speciesconcept’, in which species correspond to distinct pheno-typic clusters, analogous to Mallet’s ‘genotypic clusteringspecies concept’ [21]. In this view, the clustering of indi-viduals in trait or gene space that we recognise as speciesis a pattern that emerges from underlying ecological andevolutionary mechanisms. Just as mixtures of chemicalconstituents which react and diffuse may create patterns(for instance, spots and stripes), so individuals which react (for example, compete) and diffuse (for example,mutate in trait or gene space) may create patterns (clus-ters). The traditional approach of the theoretical physi-cist would then be to construct a simple model to see ifthe effect appears, and if so, then see if a deeper under-standing of the effect can be gained from an analysis ofthe model.This is the approach that we have adopted here. Wehave started from a simple model that only containedbirth, death (through competition) and mutation, andasked under which conditions clusters of individuals inphenotype space were formed. As discussed in the Intro-duction, this is a question which has been investigated byseveral authors, however our study focused on individu-als who gave birth only after mating with another in-dividual, whereas most previous investigations assumedasexual reproduction.We also investigated stochastic pattern formation [23,24]. Most previous work was carried out in the case of in-finite carrying capacity (in our notation, K → ∞ ) wherethe governing equations are deterministic. The standardway of proceeding in this case is to look under what con-ditions the constant density (homogeneous) solution ofthis equation is unstable. This is carried out by perform-ing a linear stability analysis about this homogeneoussolution. However it has been found that frequently pat-terns are still found in regions of parameter space wherethe homogeneous solution of the deterministic equationis stable. These patterns typically are stochastic, andcan be found by analysing the governing equations at fi-nite carrying capacity, K . Interestingly, the scaling ofthese patterns with K has particular biological implica-tions [17].The search for stochastic patterns in models of asexualreproduction has been carried out previously [15, 16]. Itwas found that noise originating in the discrete natureof individuals can lead to the spontaneous formation ofspecies in situations where this would not happen de-terministically. The main purpose of this paper was toextend this to the case of sexual reproduction. We firstsupposed that individuals played both the male and fe-male role (i.e. were hermaphrodite). We found that whenmating is assortative (i.e. organisms show preference forlike individuals) stochastic patterns can appear in the re-gion of stability of the deterministic homogeneous solu-tion. These patterns were, however, somewhat restrictedto parameter values not too far from the deterministic in-stability boundary. When mating is random the stochas-tic effects are stronger, and the phenotype distributionfor moderate K is always relatively narrow, in contrastwith the deterministic predictions. The case where theindividuals were either male or female led to very similarresults, which in some cases could be mapped directlyonto the the hermaphroditic case.There are several extensions of the current work whichcould be carried out. One could distinguish between “ge-netic noise” (arising from recombination and mutations),which affects the inheritable traits, and “environmental”or “developmental” noise, which leads to two individu-als with same genotype to have different phenotypes andwhich is not inherited. These two types of noise are likelyto have rather different effects in the phenotype distribu-tion. Also, if a dominant part of the genetic noise is dueto recombination, that is, the trait considered is deter-mined by the effect of many genes in a diploid organism,and the effect of the different genes is additive (no epis-tasis or dominance), then the noise should depend on theposition in genotype space (i.e. it would be multiplica-tive noise). Multiplicative noise would break the spa-tial symmetry and could lead to interesting effects, butwould be more difficult to study analytically. Anotherpossible extension is to include the Allee effect, which wehave here ignored for simplicity. We believe, however,that the work discussed here shows that the inclusion ofstochastic effects is vital if we wish to predict the rangeof parameters for which patterns, and therefore possiblyspecies, occur. Acknowledgments We wish to thank Tim Rogers for useful discussions.This work was supported in part by the EPSRC (UK)under grant number EP/H02171X. Appendix A: The mesoscopic evolution equation In the limit of large carrying capacity, the master equa-tion (5) can be expanded in powers of K − to give thefollowing functional Fokker-Planck equation (analogousto the derivation in Ref. [15] for the asexual case): ∂∂t P ( φ, t ) = − Z Z δδφ ( x ) [ A ( φ, x, y ) P ( φ )] dxdy (A1)+ 12 K Z Z δ δφ ( x ) [ B ( φ, x, y ) P ( φ )] dxdy, where terms of K − and higher in the expansion havebeen neglected. Here A ( φ, x, y ) = Ψ( φ, x, y ) − Ξ( φ, x, y ) , B ( φ, x, y ) = Ψ( φ, x, y ) + Ξ( φ, x, y ) , (A2)whereΨ( x, y ) = Z φ ( y ) φ ( z ) φ ∗ m ( y ) m ( y − z ) r ( x − ( y + z ) / dz, Ξ( x, y ) = φ ( x ) φ ( y ) g ( x − y ) f ( x ) . (A3)A completely equivalent way of expressing the stochas-tic dynamics of the system, is to write down the equiv-alent stochastic differential equation. This takes theform [25, 26] ∂φ ( x, t ) ∂t = Z A ( φ, x, y ) dy + η ( x, t ) √ K , (A4) where η ( x, t ) is a Gaussian white noise with zero meanand with correlator h η ( x, t ) η ( x ′ , t ′ ) i = Z B ( φ, x, y ) dyδ ( x − x ′ ) δ ( t − t ′ ) , (A5)understood in the sense of It¯o.As described in Sec. III of the main text, we can ob-tain some insight into the transition into a multi-modaldistribution by going over to a finite phenotypic space;specifically we assume that − π < x ≤ π . We begin byanalysing the deterministic dynamics of the model, whichis found by letting K → ∞ in Eq. (6). That is, ∂φ ( x, t ) ∂t = − Z φ ( x, t ) φ ( y, t ) g ( x − y ) f ( x ) dy (A6)+ Z φ ( y, t ) φ ( z, t ) φ ∗ m ( y ) m ( y − z ) r ( x − ( y + z ) / dy dz. We substitute φ ( x, t ) = 1 + ˜ φ ( x, t ) , (A7)into Eq. (A6) and only keep linear terms in ˜ φ ( x, t ). As-suming that g, m and r are periodic with period 2 π andnormalised to unity in the interval ( − π, π ), and f ( x ) = 1,one finds that ddt ˜ φ k ( t ) = (cid:2) − − g k + 2 r k m − k/ − r k m k m − k/ (cid:3) ˜ φ k ( t ) . (A8)Here we have gone over to Fourier space, since the linearnature of the problem, and the translational invariance,make this a natural choice. The Fourier modes are de-fined by h k = Z π − π h ( x ) e − ikx dx, h ( x ) = 12 π X k h k e ikx . (A9)From Eq. (A8) we see that if the condition given inEq. (10) of the main text holds, then the homogeneoussolution φ = 1 is unstable.If we wish to carry out a system-size expansion, asdiscussed in Sec. IV A, then we write φ ( x, t ) = 1 + K − / ξ ( x, t ) and expand in terms of K − / . As discussedin the main text, this leads to Eq. (11): ddt ξ k ( t ) = − ρ k ξ k ( t ) + η k ( t ) , (A10)where ρ k ≡ g k − r k m − k/ + r k m k m − k/ . Multiplyingby e ρ k t , this can be integrated to yield ξ k ( t ) = ξ k (0) e − ρ k t + e − ρ k t Z t dt ′ e ρ k t ′ η k ( t ′ ) . (A11)Assuming that we begin with zero noise, ξ k (0) = 0,Eq. (A11) implies that h ξ k ( t ) ξ − k ( t ) i = exp − ( ρ k + ρ − k ) t × Z t Z t dt ′ dt ′′ exp ( ρ k t ′ + ρ − k t ′′ ) h η k ( t ′ ) η − k ( t ′′ ) i . (A12)Using Eq. (12) and letting t → ∞ , to obtain the resultin the stationary state, one finds that if ( ρ k + ρ − k ) > t →∞ h ξ k ( t ) ξ − k ( t ) i = 2(2 π )( ρ k + ρ − k ) . (A13) Appendix B: Numerical simulation of theindividual-based model The numerical simulations of the individual-basedmodel are performed using the Gillespie algorithm [27].Before providing the details of our algorithm, we recallsome basic elements of the process. The state of thesystem at time t is described by N ( t ) real numbers, x i , i = 1 , . . . , N ( t ), corresponding to the phenotypes ofthe N ( t ) individuals present. The probability that indi-vidual i initiates a reproduction event, given that a birthtakes place, is 1 /N ( t ) (since we consider no Allee effect,the reproduction probability is independent of the den-sity of suitable mates); the probability that individual i chooses individual j to mate, given that individual i isreproducing, is proportional to m ( x i − x j ) (see Sec. II);finally, the probability density that the newborn individ-ual has phenotype x , given that individuals i and j arereproducing, is given by r ( x − ( x i + x j ) / 2) (again, seeSec. II).There are two possible events (assuming N ( t ) ≥ r = R γ ( x, φ ) dx = P N ( t ) i,j =1 g ( x i − x j ) /Kf ( x i ) (see Eqs. (1) and (3)).(ii) Birth of a new individual, with a rate r = R β ( x, φ ) dx = P N ( t ) i,j =1 m ( x i − x j ) / P l m ( x i − x l )(= N ( t )) (see Eq. (4)).Therefore the probability that individual with phenotype x i (that we will denote as individual i ) dies, given that adeath event takes place, is proportional to P N ( t ) j =1 g ( x i − x j ) /Kf ( x i ). The numerical simulations are, then, based on the fol-lowing algorithm:1. Set the initial state of the system, that is, the ini-tial number of individuals, N ( t ), their correspond-ing phenotypes, x i , i = 1 , . . . , N ( t ), and the initialtime, t = t .2. Compute r and r , as described earlier. Computethe time increment after which the next event takesplace (∆ t ), which is an exponential random variablewith average ( r + r ) − , so it can be computed as∆ t = − ln( u ) / ( r + r ), with u a pseudo-randomnumber with a uniform distribution in the interval(0 , t = t + ∆ t .3. Establish what type of event takes place. Withprobability r / ( r + r ) a death event takes place;go to 4a. Otherwise a birth event takes place; goto 4b.4a. Establish which individual dies. Choose individ-ual i at random, with probability proportional to P Nj =1 g ( x i − x j ) /Kf ( x i ). Eliminate individual i .Update N, N = N − 1. If N = 0 the populationbecomes extinct and the simulation ends.4b. Set the phenotype of the new individual. Chooseindividual i uniformly at random to initiate re-production. Choose individual j to mate, at ran-dom with probability proportional to m ( x i − x j ).Set the phenotype of the new individual as x =( x i + x j ) / ζ , where ζ is a random variable withprobability density function given by r ( x ). Update N, N = N + 15. Go to 2 or finish.When considering periodic boundary conditions, this hasto be taken into account when computing the competi-tion and mating functions as well as the phenotype of thenew individual. [1] S. Gavrilets, Fitness Landscapes and the Origin of Species (Princeton University Press, Princeton, 2004).[2] J. Maynard Smith and E. Szathmary, The Major Tran-sitions in Evolution (Oxford University Press, Oxford,1997).[3] M. Doebeli and I. Ispolatov, Science , 494 (2010).[4] J. A. Coyne, Curr. Biol. , R787 (2007).[5] D. I. Bolnick and B. M. Fitzpatrick, Ann. Rev. Ecol.Evol. Syst. , 459 (2007).[6] R. Mac Arthur and R. Levins, Am. Nat. , 377 (1967).[7] A. Sasaki, J. Theor. Biol. , 415 (1997).[8] A. J. Noest, Proc. Roy. Soc. B , 1389 (1997).[9] U. Dieckmann and M. Doebeli, Nature , 354 (1999).[10] M. A. Fuentes, M. N. Kuperman and V. M. Kenkre, Phys. Rev. Lett. , 158104 (2003).[11] M. Scheffer and E. H. van Nes, Proc. Natl. Acad. Sci.(USA) , 6230 (2006).[12] S. Pigolotti, C. L´opez and E. Hern´andez-Garc´ıa, Phys.Rev. Lett. , 258101 (2007).[13] S. Pigolotti, C. L´opez, E. Hern´andez-Garc´ıa and K. H.Andersen, Theor. Ecol. , 89 (2010).[14] M. Doebeli, H. J. Blok, O. Leimar and U. Dieckmann,Proc. Roy. Soc. B , 347 (2007).[15] T. Rogers, A. J. McKane and A. G. Rossberg, Europhys.Lett. , 40008 (2012).[16] T. Rogers, A. J. McKane and A. G. Rossberg, Phys. Biol. , 066002 (2012).[17] A. G. Rossberg, T. Rogers and A. J. McKane, Proc. Roy. Soc. B , 20131248 (2013).[18] J. Roughgarden, Theory of Population Genetics and Evo-lutionary Ecology (Macmillan, London, 1979).[19] M. G. Bulmer, The Mathematical Theory of QuantitativeGenetics (Clarendon Press, Oxford, 1980).[20] N. G. van Kampen, Stochastic Processes in Physics andChemistry (Elsevier Science, Amsterdam, 2007), 3rd ed.[21] J. Mallet, Trends Ecol. Evol. , 294 (1995).[22] J. A. Coyne and H. A. Orr, Speciation (Sinauer Asso-ciates, Suderland, MA, 2004). [23] A. J. Black and A. J. McKane, Trends Ecol. Evol. ,337 (2012).[24] A. J. McKane, T. Biancalani and T. Rogers, Bull. Math.Biol. , 895 (2014).[25] C. W. Gardiner, Handbook of Stochastic Methods (Springer, Berlin, 2009), Fourth edition.[26] H. Risken, The Fokker-Planck Equation - Methods of So-lution and Applications (Springer, Berlin, 1989), 2nd ed.[27] D. T. Gillespie, J. Phys. Chem.81