Fixation and fluctuations in two-species cooperation

Journal of Physics: Complexity PAPER • OPEN ACCESS Fixation and fluctuations in two-species cooperation To cite this article: Jordi Piñero et al 2022 J. Phys. Complex. 3 015011 View the article online for updates and enhancements. You may also like An ALMA and MagAO Study of the Substellar Companion GQ Lup B Ya-Lin Wu, Patrick D. Sheehan, Jared R. Males et al. - SCExAO/CHARIS Direct Imaging of A Low-mass Companion At A Saturn-like Separation from an Accelerating Young A7 Star Jeffrey Chilcote, Taylor Tobin, Thayne Currie et al. - Fixation in fluctuating populations Deepak Bhat, Jordi Piñero and S Redner - This content was downloaded from IP address 161.111.10.228 on 19/09/2022 at 10:27 J.Phys.Complex. 3(2022) 015011 (15pp) https://doi.org/10.1088/2632-072X/ac52e7 OPEN ACCESS RECEIVED 19 September 2021 REVISED 1 February 2022 ACCEPTED FOR PUBLICATION 8 February 2022 PUBLISHED 28 February 2022 Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. PAPER Fixation and fluctuations in two-species cooperation Jordi Piñero1,SRedner 2,∗and Ricard Sol´ e3,∗ 1ICREA-Complex Systems Lab, Universitat Pompeu Fabra, 08003 Barcelona, Spain 2Institut de Biologia Evolutiva (CSIC-UPF), Psg Maritim Barceloneta, 37, 08003 Barcelona, Spain 3Santa Fe Institute, 1399 Hyde Park Road, Santa Fe NM 87501, United States of America ∗Authors to whom any correspondence should be addressed. E-mail: [email protected],[email protected] and ricard.sol[email protected] Keywords: fixation, extinction, first passage, stochastic processes, cooperation Abstract Cooperative interactions pervade in a broad range of many-body populations, such as ecological communities, social organizations, and economic webs. We investigate the dynamics of a population of two equivalent species A and B that are driven by cooperative and symmetric interactions between these species. For an isolated population, we determine the probability to reach fixation, where only one species remains, as a function of the initial concentrations of the two species, as well as the time to reach fixation. The latter scales exponentially with the population size. When members of each species migrate into the population at rate λand replace a randomly selected individual, surprisingly rich dynamics ensues. Ostensibly, the population reaches a steady state, but the steady-state population distribution undergoes a unimodal to trimodal transition as the migration rate decreases below a critical value λc.Inthelow-migrationregime,λ<λ c,the steady state is not truly steady, but instead strongly fluctuates between near-fixation states, where the population consists of mostly A’s or of mostly B’s. The characteristic time scale of these fluctuations diverges as λ−1. Thus in spite of the cooperative interaction, a typical snapshot of the population will contain almost all A’s or almost all B’s. 1. Introduction Competitive interactions have played a prominent role in the literature of ecological and evolutionary dynamics, as well as in economics and sociology [1–3]. Resource limitations and their impact in defining the outcome of competition among species has shaped a large part of evolutionary thinking. A counterpoint to competition is cooperativity in which there are positive interactions and feedback loops between species. These mechanisms have received increasing attention recently [4,5]. In fact, it is the presence of cooperative interactions, where positive reciprocal exchanges are at work, that appear to drive innovations in evolution and also maintain biodiversity in Nature [5]. Cooperation, or mutualism, has been part of mathematical models of populations since the formulation of Lotka–Volterra equations [2] and a variety of statistical physics models of human cooperation [6,7]. In its most abstract form, two species (for example, Aand Bin figure 1(a)) ‘help’ each other by means of a mutual positive interaction; in some cases, both partners completely rely on one another for survival. This feature underlies the two-species system in figure 1(b), where a given species requires the other to replicate because each species needs a molecule that is produced by the partner species. Recent experimental studies have shown that such cooperative populations can in fact be engineered. By following the scheme in figure 1(b), it is possible to create a completely symmetric pairwise dependence and make these mixed populations grow on a Petri dish [8–10]. Figure 1(c) shows the outcome of symmetric competition when each strain is marked with a different fluorescent protein: each strain locally out competes the other, thereby generating stripes of segregated domains. The cooperatively interacting population, on the other hand, constrains both species to remain in proximity, leading to a well-mixed population (figure 1(d)). These simple engineered, or synthetic, bacterial populations, which can be tuned so that they become virtually © 2022 The Author(s). Published by IOP Publishing Ltd J.Phys.Complex. 3(2022) 015011 (15pp) JPiñeroet al Figure 1. In pairwise cooperation (a) the replication (closed arrows) of a given species A requires the help of B and vice versa. This feedback occurs, for example, when two types of bacteria A and B each lack a metabolite for reproduction that is supplied by the other species (small circles in (b)). Such cooperative feedbacks are commonplace and help maintaining diversity. Experimental set-ups using engineered bacteria (c) and (d) reveal the marked difference between competitive and cooperative interactions. In (c), two equal competitors grow on a Petri dish and locally exclude each other, as shown by the growing bands that indicate the presence of only one strain. With cooperativity (d), the mutual dependency drives both strains to persist and mix. symmetric, allow one to explore the fundamental dynamical features of interacting living consortia and also study the impact of stochasticity [11–13]. In this paper, we present an analytic approach to understand the role of stochasticity in a simple twospecies stochastic model of cooperation. This model represents a special case of evolutionary game dynamics [14–17], with a specific and particularly simple form of the payoff matrix. We emphasize that we are not treating a general ecological model, but rather a simplified system in which only cooperative interactions occur. Thisscenarioappearstobemorerelevantforthemicrobiome[18–20]. Current advances in microbial ecology involve experimental setups with a small number of interacting species [21,22]. In this model, we treat both closed and open populations, in which there is either no migration or a finite rate of migration into the population, respectively. We define microscopic rules that incorporate both cooperativity, in which each species helps the other, as well as neutrality, in which neither species is preferred over the other. We first determine the steady state of the population in the absence of stochastic fluctuations. For a finite population, we then incorporate stochasticity and determine the time until fixation is reached for the situation where no migration can occur. When migration is allowed (with compensatory removal), the population now reaches a steady state; however, the character of this steady state dramatically changes as function of the migration rate. For large migration rate, both species are present in roughly equal abundances. However, for a sufficiently small migration rate, the population strongly fluctuates between consisting of nearly all A or all B. Thus a typical realization of the population has a completely different character that the average population. This change in behavior is mirrored by a bimodal to trimodal transition in the shape of the steady-state probability distribution of species abundances. In section 2, we outline basic features of our two-species cooperation model in the absence of migration. We solve the model in the mean-field approximation and then include the role of stochasticity due to the finiteness of the population. For a finite population, we determine the fixation probability as a function of the initial population composition and the time until fixation, where only a single species remains. In section 5, we incorporate migratory inflow, with compensatory removal, so that the population size remains fixed and reaches a steady state. We discuss basic features of this steady state, including the intriguing feature that the species abundances can exhibit huge fluctuations, even though time-averaged properties are constant. We give some concluding remarks in section 6. 2. Two-species cooperation We investigate a finite population of Nparticles, with nof species A and N−nof species B. The population undergoes repeated reaction events in which each event consists of the following steps (figure 2): (a) Pick a random pair of particles. (b) If the pair is AB, one member of this pair reproduces; if the pair is AA or BB, nothing happens. (c) The newly reproduced offspring replaces one randomly selected particle in the remainder of the population. 2 J.Phys.Complex. 3(2022) 015011 (15pp) JPiñeroet al Figure 2. The reaction step in two-species cooperation. Two randomly selected particles happen to be from different species, namely A and B (red and blue). One of them reproduces (B, blue), an event that is aided by the presence of the other (A, red). The offspring replaces another randomly selected particle from the remainder of the population. In the example shown, the newly generated B replaces an A. Thus interactions between members of different species are cooperative in nature, while members of the same species are non interacting. The replacement step (iii) ensures that the total population remains fixed. The lack of interactions between AA and BB pairs follows from the assumption of strict mutualism, i.e., replication occurs if and only if both species are present (as sketched in figure 2). After each update, time is incremented by 1 N. This time increment corresponds to each particle undergoing, on average, steps (i)–(iii) in a single time unit. While this dynamics manifestly conserves the total number of particles, the composition of the population can change. When the population consists entirely of a single species—either all A’s or all B’s—there is no further dynamics and the population fixates. If an A reproduces in a single interaction, then with probability 1 −n N≡1−x, the A offspring replaces a Bandn→n+1. Conversely, with probability x=n N, the A offspring replaces an existing A so that ndoes not change. The probability anat which n→n+1 therefore is an=2x(1 −x)×1 2×(1 −x)=x(1 −x)2.(1a) Throughout, we use the variables nand x=n Ninterchangeably. In (1a), the factor 2x(1 −x) gives the probability a randomly selected pair is AB, the factor 1 2gives the probability that the A in this pair reproduces, and the factor 1 −xgives the probability that the A offspring replaces a B. By the same reasoning, when a B reproduces in an interaction, the probability at which n→n−1is bn=2x(1 −x)×1 2×x=x2(1 −x).(1b) Here the trailing factor xin (1b) accounts for the probability that the B offspring replaces an A. Finally, the probability that the number of A’s and B’s does not change is given by x2+(1 −x)2+2x(1 −x)1 2x+1 2(1 −x)=1−x(1 −x)=1−an−bn.(1c) The terms x2+(1 −x)2give the probability to pick either an AA or BB pair, for which no change in noccurs. For the last term, 2x(1 −x) is again the probability of picking an AB pair, while the factor in the square brackets is the probability that the offspring (either A or B with probability 1 2)replacesitsownkindsothatndoes not change. 3. Mean-field approaches 3.1. Rate equation Using the probabilities in equation (1),therateequationfortheaveragenumberofA’sis ˙ n=N(an−bn)=Nx(1 −x)(1 −2x), (2a) 3 J.Phys.Complex. 3(2022) 015011 (15pp) JPiñeroet al or equivalently, ˙ x=x(1 −x)(1 −2x).(2b) To keep the notation simple, nand xrefer to average values in this section; that is, we do not write angle brackets. This rate equation has a stable fixed point at x=1 2and unstable fixed points at x=0andx=1. The stability of the fixed point at x=1 2arises because the transition probabilities (1) tend to reduce population imbalances. Thus the steady state in this continuum description is a static population that consists of equal densities of A’s and B’s. That is, cooperativity promotes diversity in the mean-field description. The solution to the rate equation (2b) may be straightforwardly obtained by first performing a partial fraction decomposition: dt=dx x(1 −x)(1 −2x)=dx1 x−1 1−x+4 1−2x, from which t=x x0 dy1 y−1 1−y+4 1−2y=−4ln x(1 −x)(1 −2x) x0(1 −x0)(1 −2x0). We then obtain x(t) by solving the resulting cubic equation. For t→∞, the limiting behavior is x(t)≃1 2−2x0(1 −x0)(1 −2x0)e−t/4,(3) so that the stable fixed point x∗=1 2is approached exponentially quickly in time. 3.2. Master equation and its moments In the stochastic dynamics where nis a discrete variable, the true fixed points are at x=0andx=1. Even though the fixed point at x∗=1 2is stable in the continuum limit, stochastic fluctuations allow the population to explore the full state space and eventually get trapped at either x=0orx=1. This behavior is analogous to the extinction phenomena that arise, for example, in the logistic birth-death process, A→2Aand 2A→0, as well as other reactions of this genre [23–26]. In these reactions, the rate equation predicts a steady population, Ns, whichis determined by the balance between the birth and deathrates. However, in the true stochastic dynamics, the population fluctuates around Ns, which actually is the quasi steady-state value. Ultimately, a sufficiently large fluctuation occurs that leads to extinction, from which there can be no escape, with an extinction time that scales exponentially in Ns[23–27]. To understand the stochastic dynamics for two-species cooperation, we study Pn(t), the probability that the population consists of nA’s and (N−n)B’sattimet. The time dependence of this probability distribution is given by ˙ Pn=Nan−1Pn−1+bn+1Pn+1−(an+bn)Pn.(4) When the number of particles Nis small, the set of equation (4) can be straightforwardly solved. For the initial condition of equal numbers of A’s and B’s, both P0(t)andPN(t) approach 1 2as t→∞, while all the other Pn(t) vanish exponentially quickly in time. This direct approach quickly becomes tedious as Nincreases, however, and to gain insight into the long-time dynamics for general N, it is useful to study low-order moments of Pn. From equation (4)andusinganand bnfrom equation (1), the first moment obeys ˙ x=1 N n n˙ Pn= 1⩽n⩽N{nan−1Pn−1+nbn+1Pn+1−n(an+bn)Pn} = 1⩽n⩽N{(n+1)anPn+(n−1)bnPn−n(an+bn)Pn} = 1⩽n⩽N (an−bn)Pn=x(1 −x)(1 −2x).(5a) Here we now explicitly write angle brackets to denote average values. Under the assumption of no correlations, that is, xk=xk,(5a) reproduces the rate equation (2b). Similarly, the equation of motion for the second moment is 4 J.Phys.Complex. 3(2022) 015011 (15pp) JPiñeroet al ˙ x2=1 N2 n n2˙ Pn=1 N 1⩽n⩽Nn2an−1Pn−1+n2bn+1Pn+1−n2(an+bn)Pn =1 N 1⩽n⩽N(n+1)2anPn+(n−1)2bnPn−n2(an+bn)Pn =1 Nx(1 −x)+2x2(1 −x)(1 −2x).(5b) It is more convenient to express (5a)and(5b)intermsofz≡2x−1, which lies in the range [−1, 1]. Doing so, we obtain ˙ z=−1 2z(1 −z2) ˙ z2=(1 −z2)1 N−z2, (6) which are both symmetric about z=0. If we make the assumption of no correlations, that is, zk=zk,then the first equation reproduces the result that z=0 is a stable fixed point. The second equation then predicts that the width of the distribution initially grows and eventually ‘sticks’ at the value √N. To check this point, we numerically integrated equation (4) for small values of N. From the resulting solution, we find that the width of the probability distribution initially grows with time and later approaches a nearly fixed value that is proportional to √N. However, at very long times, there is slow leakage of the probability distribution to the true stochastic fixed points at z=±1. Thus the probability distribution eventually approaches two deltafunction peaks at these fixed points. This behavior cannot be captured by low-order moment equations, such as (6). Instead, we need to study the full stochastic dynamics; this is done in the following section. 4. Fixation probability and fixation time We now turn to two quantities of primary interest in the stochastic dynamics, namely, (i) the fixation (or exit) probability En, and (ii) the fixation time Tn. The fixation probability Enis defined as the probability that a population of size Nthat initially contains nparticles of type A reaches the static fixation state of all A’s. We use the backward Kolmogorov equation [28,29] to compute the fixation probability. In this approach, Ensatisfies the recursion En=anEn+1+bnEn−1+(1−an−bn)En.(7) Since the process renews itself after each event, we can express the fixation probability from the state that contains nA’s in terms of the appropriately weighted average of the fixation probabilities after a single step to the states n−1, n,andn+1. The weights are merely the hopping probabilities to these respective states. Equation (7) is subject to the boundary conditions E0=0andEN=1. The first condition corresponds to the impossibility of reaching a population of all A’s if the initial state contains no A’s, while the second condition corresponds to the initial state coinciding with the desired final state of all A’s. The solution to (7) is (see appendix Afor the calculational details) En= n−1  m=0N−1 m−1N−1  m=0N−1 m−1 .(8) Neither of these sums has a closed form, but for N→∞the denominator approaches 2 [30]. For N1, En and its continuum counterpart E(x) (see also appendix A) are nearly independent of x=n Nwhen xis not close to 0 or 1. Figure 3showsthisdependenceofEnon n. Also shown are the corresponding results from discrete simulations of the fixation process. Simulations are carried out by setting a finite size (N) array of particles with two possible states. Each iteration the particles interact following rules (ii) and (iii) with the interaction rates given by expressions (1a)and(1b). Time is updated by Δt=[N(an+bn)]−1after each iteration. The anti-sigmoidal shape of Enarises from the underlying drift that tends to drive any initial population towards x=1 2. Eventually, a large and rare stochastic fluctuation causes the population to escape this effective potential well and reach fixation. This anti-sigmoidal shape also strongly contrasts with the Moran process [31], which is symmetric (neutral), but non-cooperative. Here an AB pair equiprobably converts to either AA or to BB. As a result of this lack of cooperativity, the fixation probability in the strictly neutral Moran process is simply the linear function E(x)=x[29,31–33]. We now investigate the fixation time Tn, which is defined as the average time for the population of N particles to first reach either of the two fixation states, n=0orn=N, when the population initially contains 5 J.Phys.Complex. 3(2022) 015011 (15pp) JPiñeroet al Figure 3. Dependence of the discrete and continuum fixation probabilities, Enand E(x), versus xfor the cases N=8 and 16. The smooth curves represent E(x)fromequation(A.5) and the dots represent simulation results. nA’s. Within the same backward Kolmogorov framework as that used for the fixation probability, the fixation time satisfies [28,29] Tn=anTn+1+bnTn−1+(1−an−bn)Tn+δt.(9) Again Tnmay be expressed as the appropriately weighted average of the fixation time after a single hopping event to the states n−1, n,andn+1, plus the time δt=1 Nrequired for this single step. The latter corresponds to each particle being updated once, on average, in a single time unit. The equation for Tnis subject to the boundary conditions T0=TN=0; namely, if the population starts in a fixation state, the time to reach this state is 0. The result for the fixation time is (see appendix B) Tn=En N−1  m=1 Qm− n−1  m=1 Qm, (10) where Qn≡αn+rnαn−1+rnrn−1αn−2+···+rnrn−1...r2α1, with Q0=0, rn=bn/an,andαn≡δt/an. It does not seem possible to reduce (10) to a compact form, but the main feature of this exact expression is that the fixation time scales exponentially in Nand is nearly independent of n(or, equivalently, x), except for nclose to 0 or to N(figure 4). The exponential dependence on Nagain arises because of the existence of an effective potential well, whose depth grows linearly with N, which draws the population toward the state x=1 2. The near independence of the fixation time on the initial condition is a consequence of the population being drawn toward the bottom of this potential well, where the concentrations of A and B are equal. As a result, the value of the fixation time for any initial value of xis close to the fixation time when the population starts from the bottom of the potential well at x=1 2. It is possible, however, to obtain an analytical expression for the average fixation time by the WKB method [23–26]. The idea of this approach is that the probability distribution settles into a quasi-steady state that assumes an exponential large-deviation form. From equation (4), there is a slow leakage from this quasi-steady state to the fixation state whose rate, Γ(N), is given by Γ(N)δt=b1 P1+aN−1 PN−1=2b1 P1, (11) i.e., the flux from states that are one step away from fixation to the fixation states. Here the tilde denotes the steady-state distribution and we also use the symmetry n↔N−n. We then identify the inverse of this leakage rate with the fixation time. We obtain an approximate equation for the continuum probability distribution  Pn→ P(x) by setting the time derivative in the master equation (4)tozero,andwritingn±1asx±δxto give a(x−δx) P(x−δx)+b(x+δx) P(x+δx)=[a(x)+b(x)] P(x).(12) We now assume that  P(x) has the exponential form  P(x)∼eS(x)/δx=eNS0(x)+S1(x)+··· and substitute this form into (12)togive(uptoO(1)) S0(x)=x dzlog a(z) b(z),S1(x)=−1 2log [a(x)b(x)].(13) 6 J.Phys.Complex. 3(2022) 015011 (15pp) JPiñeroet al Figure 4. (a) Dependence of the fixation time Tnversus x=n Nusing a data-collapse scheme by resetting the scale Tn→Ne−Nlog2Tnfor each Nvalue. This scaling factor corresponds to the prediction made in equation (16). The solid lines follow the predictions obtained from equation (10). The plot marks {×,+and ∗} correspond to simulated average fixation times for N=8, 16 and 24 in blue, green and red colors, respectively. Each data point is obtained by averaging 103simulated fixation processes at corresponding values of nand N. (b) Dependence of the fixation time from the symmetric initial state, TN/2(red dots) computed by (10) and the WKB prediction for the inverse leakage rate from the quasi-steady state (blue line), following (16). Now using (1a)and(1b)fora(x)andb(x), we have ˜ P(x)∼eNS0(x) x3/2(1 −x)3/2, (14) with S0(x)=−xlog x−(1 −x)log(1 −x). Note that the action S0(x) is peaked at the quasi-state state x=1 2. We normalize  P(x) by using the Laplace method for N→∞[34], 1 0 eNS0(x) x3/2(1 −x)3/2dx≈32 NeNlog 2, so that  P(x)≃N √32π eN[−xlog x−(1−x)log(1−x)−log 2] x3/2(1 −x)3/2.(15) We now compute the fixation rate Γby substituting  P1 N≃N3 √32πe1−Nlog 2 and b1 N≃1 N2, into equation (11)togive Γ(N)=Ne √8πe−Nlog 2.(16) We now identify the inverse of this rate with the average fixation time. As shown in figure 4(b), this inverse rate accurately matches the simulation data for the fixation time. 5. Two-species cooperation with migration We now incorporate migration into the dynamics, in which particles of either species migrate into the population at the same fixed rate λ, and each new particle replaces a randomly selected existing particle. Because migration is accompanied by replacement, the population size remains fixed, which is the physically most relevant case. Now the population is driven to a steady state rather than to fixation and we want to understand thenatureofthissteadystate. 5.1. Probability distribution For a population that consists of nA’s and (N−n) B’s, suppose that the migrant is an A. With probability 1 2(1 −x), the A migrant replaces a B and n→n+1, while with probability 1 2x, the A migrant replaces an A, and the composition of the population remains the same. Similar reasoning applies when the migrant is a B. As a result of a migration event, the average change in the number of A’s is 1 2(1 −x)−1 2x.Therateequation for nnow is (compare with equation (2)) ˙ n=N(1 −λ)[x(1 −x)(1 −2x)]+1 2Nλ(1 −2x).(17) 7 J.Phys.Complex. 3(2022) 015011 (15pp) JPiñeroet al For λ>0, x=0andx=1 are no longer fixed points and only the remaining fixed point at x=1 2is stable. In the absence of fluctuations, the population is thus driven to a steady-state distribution, Pn(t→∞), that is peaked about x=1 2. Because there is no absorbing state in the stochastic dynamics, we might anticipate a similarbehaviorforPn(t→∞) when stochasticity is accounted for. We will show, however, that within the Fokker–Planck approximation the steady-state distribution can either be unimodal or trimodal in shape and the latter case corresponds to a steady state that is not truly steady. The probability distribution Pnis now governed by the master equation ˙ Pn=N(1 −λ)an−1Pn−1+bn+1Pn+1−(an+bn)Pn +Nλcn−1Pn−1+dn+1Pn+1−(cn+dn)Pn, (18) with hopping probabilities due to migration that are given by cn=1 21−n Ndn=1 2 n N.(19) We now determine the continuum probability distribution in the Fokker–Planck approximation. As we shall see, this continuum expression for the probability distribution matches simulation data quite well, thus justifying the Fokker–Planck approximation ex post facto as a way to probe steady-state properties. In terms of x=n N,dx=1 N,Pn→P(x), we expand (18) in a Taylor series up to second order. This gives the Fokker–Planck equation [28,35] Pt=−(1−2x)(1−λ)x(1−x)+λ 2P(x,t) x +1 2N(1−λ)x(1−x)+λ 2P(x,t) xx ≡−{v(x)P(x,t)}x+{D(x)P(x,t)}xx, (20) where the subscripts denote partial derivatives. The steady state is defined by solving this equation with the left-hand side set to zero. Integrating once gives (DP)x−vP =B,whereBis a constant. We determine the constant by evaluating this equation at the symmetry point x=1 2. Because the probability distribution is symmetric about x=1 2,Px(x=1 2)=0. Moreover, at x= 1 2,v=0andDx=0, which implies that B=0. Thus we only need to solve (DP)x−vP =0, whose solution is P(x)=Cexp x dyv(y)−Dy(y) D(y)=Cexp −log D(x)+x dyv(y) D(y) =C D(x)exp x dy2N(1 −2y) =C1 (1−λ)x(1−x)+λ 2e2Nx(1−x), (21) where the constant Cis determined by normalization. For λ→0, P(x)isconcentratednearx=0 and near x=1; these peaks correspond to the near-fixation states. Because of rare fluctuations, however, the population stochastically switches between states where almost all particles are of type A to states where almost all particles are of type B. Naively, one therefore anticipates that the steady-state distribution should be bimodal, with a peak at each of the two near-fixation states. Unexpectedly, there always remains a peak at x=1 2(which may be vanishingly small), so that the this distribution is trimodal in the small-λregime. As λincreases beyond a critical value, the steady-state distribution undergoes a trimodal to unimodal transition (figure 5). For fixed N, we determine the transition between trimodality and unimodality by finding the point(s) where P(x)=0. This calculation gives, after straightforward algebra, P(x)∝(1 −2x)e 2Nx(1−x)2N−(λ−1) D(x)2, where again D(x) is the diffusion coefficient defined by equation (20). The leading factor of 1 −2xequals 0 at x=1 2and corresponds to the extremum in the distribution at x=1 2. However, there are additional extrema at the points where the factor in the square brackets equals 0. To determine these extrema, we first determine the zero of this factor at x=0, 1. Thus we have the condition 2N=(1 −λ)/(λ/2)2.Sincewewillfindthatλ1, we also neglect λcompared to 1 to give λc=2 N.(22) 8 J.Phys.Complex. 3(2022) 015011 (15pp) JPiñeroet al [4] May R M, Levin S A and Sugihara G 2008 Nature 451 893–4 [5] Bronstein J L 2015 Mutualism (Oxford: Oxford University Press) [6] Perc M 2016 Phys. Lett. A380 2803–8 [7] Perc M, Jordan J J, Rand D G, Wang Z, Boccaletti S and Szolnoki A 2017 Phys. Rep. 687 1–51 [8] Mitri S, Clarke E and Foster K R 2016 ISME J. 10 1471–82 [9] Nadell C D, Drescher K and Foster K R 2016 Nat. Rev. Microbiol. 14 589–600 [10] Müller M J I, Neugeboren B I, Nelson D R and Murray A W 2014 Proc. Natl Acad. Sci. 111 1037–42 [11] Shou W, Ram S and Vilar J M G 2007 Proc. Natl Acad. Sci. 104 1877–82 [12] Amor D R, Montañez R, Duran-Nebreda S and Sol´ e R 2017 PLoS Comput. Biol. 13 e1005689 [13] Amor D R and Dal Bello M 2019 Life 922 [14] Nowak M A 2006 Evolutionary Dynamics: Exploring the Equations of Life (Cambridge, MA: Harvard University Press) [15] Antal T and Scheuring I 2006 Bull. Math. Biol. 68 1923–44 [16] Altrock P M and Traulsen A 2009 New J. Phys. 11 013012 [17] Black A J, Traulsen A and Galla T 2012 Phys.Rev.Lett.109 028101 [18] Nadell C D, Foster K R and Xavier J B 2010 PLoS Comput. Biol. 6e1000716 [19] Rakoff-Nahoum S, Foster K R and Comstock L E 2016 Nature 533 255–9 [20] Foster K R, Schluter J, Coyte K Z and Rakoff-Nahoum S 2017 Nature 548 43–51 [21] Friedman J and Gore J 2017 Curr. Opin. Syst. Biol. 1114–21 [22] Vega N M and Gore J 2018 Curr. Opin. Syst. Biol. 45 195–202 [23] Elgart V and Kamenev A 2004 Phys. Rev. E70 041106 [24] Kessler D A and Shnerb N M 2007 J. Stat. Phys. 127 861–86 [25] Assaf M and Meerson B 2010 Phys. Rev. E81 021116 [26] Assaf M and Meerson B 2017 J. Phys. A: Math. Theor. 50 263001 [27] Krapivsky P L, Redner S and Ben-Naim E 2010 A Kinetic View of Statistical Physics (New York: Cambridge University Press) [28] Van Kampen N G 1992 Stochastic Processes in Physics and Chemistry vol 1 (Amsterdam: Elsevier) [29] Redner S 2001 A Guide to First-Passage Processes (New York: Cambridge University Press) [30] Lee S et al 2012 Mathematics Stack Exchange https://math.stackexchange.com/questions/151441/calculate-sums-of-inversesof-binomial-coefficients [31] Moran P A P et al 1962 The Statistical Processes of Evolutionary Theory (Oxford: Oxford University Press) [32] Kimura M 1983 The Neutral Theory of Molecular Evolution (Cambridge: Cambridge University Press) [33] Ewens W J 2012 Mathematical Population Genetics 1: Theoretical Introduction (New York: Springer) [34] Bender C M and Orszag S A 1999 Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory (New York: Springer) [35] Gardiner C W et al 1985 Handbook of Stochastic Methods 3rd edn (Berlin: Springer) [36] Fichthorn K, Gulari E and Ziff R 1989 Phys.Rev.Lett.63 1527 [37] Considine D, Redner S and Takayasu H 1989 Phys. Rev. Lett. 63 2857 [38] Carro A, Toral R and Miguel M S 2016 Sci. Rep. 624775 [39] Herrerías-Azcu´ e F and Galla T 2019 Phys. Rev. E100 022304 [40] Lambiotte R and Redner S 2007 J. Stat. Mech. L10001 [41] Bunin G 2017 Phys. Rev. E95 042414 [42] Pearce M T, Agarwala A and Fisher D S 2020 Proc. Natl Acad. Sci. USA 117 14572–83 [43] Ferry M S, Razinkov I A and Hasty J 2011 Methods Enzymol. 497 295–372 [44] Luke C S, Selimkhanov J, Baumgart L, Cohen S E, Golden S S, Cookson N A and Hasty J 2016 ACS Synth. Biol. 58–14 [45] Karlin S and Taylor H M 1975 A First Course in Stochastic Processes (San Diego, CA: Academic) 15