Non-Selective Evolution of Growing Populations

Non-selective effects, like genetic drift, are an important factor in modern conceptions of evolution, and have been extensively studied for constant population sizes (Kimura, 1955; Otto and Whitlock, 1997). Here, we consider non-selective evolution in the case of growing populations that are of small size and have varying trait compositions (e.g. after a population bottleneck). We find that, in these conditions, populations never fixate to a trait, but tend to a random limit composition, and that the distribution of compositions “freezes” to a steady state. This final state is crucially influenced by the initial conditions. We obtain these findings from a combined theoretical and experimental approach, using multiple mixed subpopulations of two Pseudomonas putida strains in non-selective growth conditions (Matthijs et al, 2009) as model system. The experimental results for the population dynamics match the theoretical predictions based on the Pólya urn model (Eggenberger and Pólya, 1923) for all analyzed parameter regimes. In summary, we show that exponential growth stops genetic drift. This result contrasts with previous theoretical analyses of non-selective evolution (e.g. genetic drift), which investigated how traits spread and eventually take over populations (fixate) (Kimura, 1955; Otto and Whitlock, 1997). Moreover, our work highlights how deeply growth influences non-selective evolution, and how it plays a key role in maintaining genetic variability. Consequently, it is of particular importance in life-cycles models (Melbinger et al, 2010; Cremer et al, 2011; Cremer et al, 2012) of periodically shrinking and expanding populations.


Introduction
Stochastic effects play an important role in population dynamics [8][9][10][11], particularly when competition between individuals is non-selective. Most previous theoretical analyses have studied how a non-selectively evolving trait can spread and eventually replace all other variants (fixate) under conditions in which the population size remains constant [2,12,13]. However, both natural and laboratory populations frequently experience exponential growth. Here we show that genetic diversity in growing populations is maintained despite demographic noise, and reaches a stationary but random limit. We used a well-controlled model system in which well-mixed co-cultures of a wild-type Pseudomonas putida strain and an isogenic mutant were grown under non-selective conditions. Multiple subpopulations were generated, each containing a random number of individuals of each strain. Depending on the average initial population size and the strain ratio, we observed distinct stationary probability distributions for their genetic composition. Moreover, we showed that the dynamics of growing populations can be mapped to Pólya urn models [4], permitting the observed maintenance of genetic diversity to be understood as the random limit property of a fair game between individual strains. Generalizing the Pólya urn model to include the effects of random initial sampling and exponential growth allowed us to predict the evolution of the composition distribution. Using numerical and analytical methods we found that the distribution broadens at first but quickly "freezes" to a stationary distribution, which agrees with the experimental findings. Our results provide new insights into the role of demographic noise in growing populations.

Results and discussion
Evolutionary dynamics is driven by the complex interplay between selective and non-selective (or neutral) effects. The paradigm of non-selective evolution originates from the seminal work of Kimura [1], in which he solved the Wright-Fisher model, thus showing that non-selective effects -and specifically genetic drift -can have a determinant role in evolution. His results sparked an ongoing debate about the nature and potency of randomness as a fundamental evolutionary force [13,18,19]. For very small populations genetic drift is generally considered an important factor [13], as the theory successfully predicts the outcomes of neutral evolution experiments [9,20].
In most theoretical analyses, constant (or effectively constant) population sizes are assumed, and the role of population growth is neglected. Bacterial populations, however, often undergo rapid growthespecially when they are small. For example, as few as 10 individuals of some highly virulent pathogens (e.g. enterohemorrhagic Escherichia coli or Shigella dysenteriae) suffice to initiate a deadly infection in a human host [21,22]. Another case of small, growing populations are water-borne bacteria that feed on phytoplankton products. Due to nutrient limitation in open water, these bacteria typically live in small populations in close proximity to the planktonic organism [23]. During spring blooms, the phytoplankton releases more organic material, boosting the bacterial growth rate [23][24][25]. In nature, such small populations often form by adventitious dispersal from a larger reservoir population [26]. A typical example is the spreading of pathogens from host to host. This random "sampling" from a reservoir yields small populations whose genetic compositions differ from that of the reservoir (a phenomenon known as the founder effect [27]). Recent studies also showed that the combination of population growth and stochastic fluctuations can have a major impact on the evolution [5-7, 28] and genetics [29] of small populations.
To probe how population growth shapes genetic diversity, we used a well-characterized microbial model system, namely the soil bacterium Pseudomonas putida KT2440 [3,30,31]. The wildtype strain KT2440 produces pyoverdine, an ironscavenging molecule that supports growth when iron becomes scarce in the environment. Here we consider co-cultures of two genetically distinct strains: the wild-type, pyoverdine-producing strain KT2440 (strain A) and the mutant non-producer strain 3E2 (strain B). We set up conditions of non-selective competition between these strains by using an ironreplete medium (casamino acids, supplemented with 200 µM FeCl 3 ). In this medium, production of pyoverdine is effectively repressed [31], such that both strains have the same growth rate and neither has an advantage (see SI). Producer (KT2440 wild type) and non-producer (3E2) strains were first mixed and diluted to yield Poisson dilution conditions. Then we initiated a large number of subpopulations from this reservoir by pipetting aliquots of the cell suspension into the wells of a 96-well plate, thereby generating a large ensemble of subpopulations with a random distribution of initial cell number N 0 and producer fraction x 0 ( Figure 1). Use of shaken liquid cultures ensured homogeneous well-mixed conditions for all cells in the same well (access to nutrients, oxygen, etc.), and exponential growth was observed in all cases.
This experimental setting is well described within the mathematical framework of a Pólya urn model. Consider each bacterium in the population as a marble in an urn, and its genotype as the color of the marble (e.g. red for strain A, and blue for strain B). Population growth results from single reproduction events in which an individual randomly divides. This is mathematically equivalent to a stochastic event in which a marble is chosen at random from the urn and put back, together with another one of the same color. This random process, introduced by Eggenberger and Pólya [4], exhibits several important properties [32][33][34][35]. It is self-reinforcing: each time a marble is extracted, another one of the same color is added, increasing the likelihood of extracting a marble of that color again. In the context of bacterial populations, this means that every birth event for one strain makes it more likely that further birth events of that same strain will occur in the future. Note, however, that fixation, i.e., complete loss of one type of marble from the population, cannot occur, simply because in the Pólya urn model marbles are neither removed nor do they change their color. This fully reflects the experimental conditions: During exponential growth, rates of cell death are negligible, and within the observation time mutations will be extremely rare, given the population sizes considered. The bacteria in each well reproduce randomly at a per-capita (average) rate µ. To translate this conditions. An infinite reservoir contains a diluted mixture of bacteria, a fractionx 0 of which are of strain A. We draw small volumes of liquid from the reservoir containing small, random numbers of individuals, which conform to a Poisson distribution with mean (determined by the dilution of the reservoir population). A certain fraction of this initial population is of strain A. The mean value of this fraction is equal tox 0 . We use these individuals to initiate populations in the wells of a microtiter plate, so that each population starts with a random size N 0 and a random fraction of A-individuals x 0 . (b) Illustration of the Pólya urn model. If a bacterial population is represented as an urn, each individual as a marble and each bacterial strain as a color, this urn model captures the essentials of bacterial reproduction in our populations. At each iteration, a marble is drawn at random and returned to the urn, together with another one of the same color. The probability of extracting a marble of either color is determined solely by its relative abundance, making the process non-selective (since no strain has inherent advantages). The rate of growth in population size can be rendered exponential by letting the waiting time between successive iterations be exponentially distributed (also known as Poissonization). to the urn model, drawing of a marble is assumed to be a stochastic Poisson process, with a "per-marble" rate µ (a procedure known as Poissonization or embedding [36,37]). Mathematically, the growth process in then described by a Master Equation: The time evolution for the probability P (N A , t) of finding N A individuals of strain A at time t reads where we have set the growth rate to µ = 1 in order to fix the time scale (for an introduction to the mathematical concepts see, e.g., [38]); the corresponding Master equation for individuals of strain B is of identical form. To study the composition of the populations, we use the more convenient quantities N = N A + N B (total size) and x = N A /N (fraction of individuals of strain A).
To start the experiment, we inoculated the wells of 96-well-plates by drawing small volumes of diluted liquid bacterial culture from a large reservoir ( Fig. 1(a)). Each volume contains a random number of bacteria whose mean value is controlled by the dilution of the reservoir. The fraction of bacteria of strain A (wild type) in that volume is also random, with its mean valuex 0 given by the fraction of strain A in the reservoir. In the mathematical formulation, this setup corresponds to stochastic initial conditions for the Pólya urn model: the initial population size N 0 for each well is given by a Poisson distribution with meanN 0 , and each individual is assigned to strain A or B with probabilityx 0 and 1 −x 0 ), respectively. This procedure is also equivalent to treating the initial numbers of A-and B-individuals as independent, Poisson-distributed random variables with mean valuesN 0x0 andN 0 (1 −x 0 )), respectively [6]. Figure 2 shows a time series of the histogram for the composition x of all subpopulations considered, as obtained from a stochastic simulation of the Master Equation (1) for a given random initial condition (withN 0 = 10 andx 0 = 0.33). Surprisingly, the distribution first broadens, but then quickly "freezes" to a steady state (see supporting video). This is genuinely different from Kimura's result for populations with constant size [1] (or similar results with effectively constant size [2]) where the balance between stochastic birth and death events leads to genetic drift, and thereby eventually to the extinction of one of the two strains. In contrast, for a growing population, death events are negligible, and therefore there is no fixation of the population during growth. Instead, fixation arises as a direct consequence of the initial sampling process, as can be seen from the heights of the black bins in the histogram (at x = 0 and x = 1), which remain constant over time (Fig. 2). During growth, the composition of each subpopulation, instead of drifting to fixation at either x = 0 or x = 1, reaches a stationary limit value x * , where it remains thereafter [39]. This limit value is random: starting several subpopulations (urns) from the exact same initial composition of strain A and B (blue and red marbles), each reaches a limit, but in general these limits differ from one another. Once all of the subpopulations in an ensemble reach their limit, the distribution of the population composition freezes to a steady state, which is equal to the probability distribution of x * . Similar random limit properties appear in other fields, with lock-in in economics as the best-known example [34].
The inset in Fig. 2 shows approximate solutions for the time evolution of mean, standard deviation, and skewness of the composition x, which we obtained by analytically solving the Master Equation (1) (see SI). The analytical results agree well with their numerical counterparts. In particular, the mean value remains constant over time, as it must for a non-selective process. For the time evolution of the variance, which is a measure for the spread of a distribution, we obtain to leading order in population size The broadening and freezing of the distribution is reflected in the exponential decay term of the variance. Note that the skewness increases as well, because growth is self-reinforcing (see inset in Fig. 2).
To further test the validity of the stochastic simulations, we also calculated the limit values of the average and variance after extended periods of evolution exactly, and found that they match the numerical solutions of the Master Equation perfectly (see SI). We tested these theoretical predictions using P. putida as a bacterial model system. We mixed the wild-type and mutant strains in order to obtain different initial fractionsx 0 . The degree of dilution of the mixture determines the average initial cell num-berN 0 , with which we inoculated 120 wells per ex- In each well the population follows a stochastic path and reaches a (random) limit composition, and the distribution freezes only when all populations reached their limit. The parameter values used in the simulation areN 0 = 10 andx 0 = 0.33 The inset shows the mean, standard deviation and skewness as a function of the number of generations, with symbols denoting numerical simulations, and the solid lines corresponding to the theoretical prediction of Eq (2) and the analogous ones for the other moments (see SI). Analytical and numerical values agree. The mean x remains constant throughout the evolution, as expected for a non-selective process; standard deviation and skewness saturate to limit values, confirming the freezing of the distribution. periment (96-well plate format). In order to compare the experimental data with our model, we set up simulations that matched the experimental configuration by initializingN 0 andx 0 with the same values as measured in the experiments. We simulated the time evolution of about 10 4 populations, grouped in "virtual plates" of 120 wells each. Every virtual plate produced a histogram like the one we obtained from experiments. We then generated an average histogram of the virtual plates and used its values to compute the binomial confidence intervals [40] for the count in each bin, and compared those with the experimental distribution. Figure 3(a) shows a representative experimental histogram of the initial population sizes N 0 for strong dilution withN 0 = 2.55. It is well approximated by a Poisson distribution, and agrees with the simulation results within statistical errors (blue line and shaded gray areas in Fig. 3(a)). Figure 3(b) shows the probability distribution of the corresponding initial compositions x 0 of the populations, where again theoretical and experimental values agree well within statistical error. Note also that in every well the composition x 0 must be a simple fraction; this means that only a few numerical values are possible for small initial population sizes N 0 . This smallnumber effect explains why the distribution of x 0 in Fig. 3(b) is so ragged. The distribution becomes much smoother for larger initial population sizes (see SI). Taken together, these results for the distribution of initial population size and composition confirm that the inoculation of the individual wells is a stochastic sampling process with Poissonian statistics.
Next, we were interested in how the composition of the bacterial population would evolve under nonselective (neutral) growth conditions. To this end we let the 120 populations grow for an 11-hour period, during which they remained in exponential growth phase. Then we measured the population size N (t) in each well by counting colony-forming units, and x(t) by counting the pyoverdine-producing colonies (see Materials and Methods). Figure 4 shows the final outcome for four different initial conditions, i.e. combinations of the initial average population sizē N 0 and compositionx 0 . We first wanted to know what determines the number of wells that contain only individuals of either strain A or strain B, i.e. that are fixated. To this end we compared the experimentally observed values with the correspond-  N 0 (panel (a)) and x 0 (panel (a)) are measured from 120-well ensembles. The averageN 0 andx 0 calculated from the measured values determine the parameters for the simulated distributions. The theoretical average distribution (solid blue line) is the average of the same distributions generated for 84 sets of 120 wells. Using that average we calculate the Wilson binomial confidence intervals (gray areas) for 68% (between dashed lines), 95% (between dotted lines) and 99.73% confidence. The measured and simulated distributions agree well within statistical error, confirming our assumption that individuals of strain A and B in the experiments start Poisson-distributed with meanN 0x0 andN 0 (1 −x 0 ), respectively. The ragged distribution of x 0 derives from a small-number effect, and disappears at larger values of N 0 .
ing predictions from the numerical simulations of the Pólya urn model (Fig. 4). Since both results agree within statistical error, we conclude that fixation of a population is a consequence of the initial sampling process and is not due to fixation during population growth. This is especially obvious for small average initial population size or compositions close to x = 0 or x = 1, where a large fraction of the wells contains cells of strain A or B only ( Fig. 4(a) and Fig. 4(d)). Next we wished to learn how the final distribution of the population composition (i.e. the random limits, x * ) depends on the initial average composition x 0 . Forx 0 = 0.5, we observed both by experiment and theoretically that the initial distribution significantly broadened (by a factor of √ 2) but remained symmetrical (Fig. 4(c)). In contrast, starting from distributions with average values below or above 0.5 caused the final distribution to broaden and also become skewed towards smaller or larger values of x, respectively ( Fig. 4(b) and Fig. 4(d)). Moreover, we found quantitative agreement between experiment and numerical simulations within statistical errors in all analyzed parameter regimes (see blue lines and shaded areas in Fig. 4): most experimental histograms fall within the first confidence interval of the prediction (darkest gray areas, between dashed lines), and almost all fall within the 99.73% confidence interval.
Figure 4: Steady-state distributions of population composition x for different initial conditions. The experimental distribution (bars) is the result of growth on 120 independent wells. We use the measured average x 0 and N 0 from the experiments to initialize the simulations of several 120-well ensembles. After growth, we compute the histogram for each of these ensembles and obtain the average theoretical distribution (blue line). Using the values from this distribution, we compute the three confidence intervals (shaded gray areas) for each bin for 68% (between dashed lines), 95% (between dotted lines) and 99.73% confidence. The two sets of data match: most experimental data agree with the first prediction confidence region, practically all with the second one. The limit distributions are clearly different from the initial ones (see SI). The importance of growth in changing the distributions depends on the initial size N 0 . Parameter values:N 0 = 2.9,x 0 = 0.32 (panel (a));N 0 = 18.4, Taken together, our combined theoretical and experimental analysis gives a coherent picture of evolution during non-selective (exponential) growth. We have shown, experimentally and by analogy with the Pólya urn model, that for each well-mixed population the composition of the population reaches a random stationary limit , and, unlike populations with constant size, generally does not fixate. For a large ensemble of populations, this implies that the probability distribution for the population composition converges to limit distributions ( Fig. 2 and Fig. 4), which are nothing like Kimura's result for constant-sized populations. Our result is also quite different from that obtained in range expansion experiments [41] or other settings featuring population growth without death on two-dimensional substrates. There, monoclonal sectoring patterns arise as a consequence of random genetic drift, which drives population differentiation along the expanding fronts of bacterial colonies, unlike our well-mixed populations that freeze to coexistence.
Our study also shows that, in a growing population with stochastic initial conditions, demographic noise has two possible sources: the initial sampling process by which subpopulations are formed, and the subsequent growth process. The initial average population sizeN 0 sets their relative weight (see SI). For very smallN 0 , of the order of one or two individuals, the formation process already determines the final composition distribution: most populations start off fixated, many with just a single founder individual, and the composition of each well remains the same during growth. For very largeN 0 , of the order of a few hundreds, the sampling process is again central: the composition distribution changes very little before freezing, and growth generates only a very limited amount of variation. In these two limiting regimes, neglecting stochastic effects during growth leaves the evolutionary outcome practically unchanged. In contrast, for small founder colonies such as those typically found during population bottlenecks [22,23,42] (N 0 ∼ 10), population growth is responsible for the major part of the variation observed in the final distribution.
Moreover, our results reveal that a growing population reaches a random limit composition much faster than genetic drift leads to fixation in populations of constant size. Typical fixation times for genetic drift increase logarithmically with the population size [11], while the time scale for freezing is independent of population size. This has important consequences for the role of stochastic effects when a population passes from exponential growth phase to stationary phase, in which growth rate and death rate are equal. Then, the composition of the population shows both freezing and fixation, albeit at quite distinct times because the relevant time scales differ markedly. During growth the composition distribution quickly freezes, as described above. Once the population reaches its stationary size, it slowly drifts to fixation, following Kimura-like dynamics.
We also believe that our results have a broad range of applications since growing populations are ubiquitous in nature. For example, experimental studies of P. aeruginosa [26,43] have shown that typical life cycles pass through different steps with regularly occurring dispersal events being followed by the formation of new colonies. As initial colony sizes are typically small, such dispersal events coincide with population bottlenecks and subsequent exponential growth. During these phases of the life cycle, population dynamics is often selectively neutral and hence falls within the framework of the present work. The degree of diversity generated during these population bottlenecks has been shown to be crucial for some proposed mechanisms for the evolution of cooperation under selective pressure [5-7, [44][45][46]. Our analysis quantifies the ensuing degree of diversity and points to the relative importance of sampling versus growth for long-term behavior of the reservoirs. This may have important consequences for the degree of genetic diversity observed in natural populations with life-cycle structures [42].

Strains and cultivation conditions
The P. putida strains KT2440 (wild type) and 3E2 (mutant with defective pyoverdine synthesis) [3] were used as pyoverdine producers and non-producers, respectively. Cells were grown in casamino acid medium (CAA) containing per liter: 5 g casamino acids, 0.8445 g K 2 HPO 4 , 0.1404g MgSO 4 •(H 2 O) [3]. The CAA medium was supplemented with 200 µM FeCl 3 (CAA-Fe) to suppress pyoverdine production (see Table S9). Overnight cultures of the individual strains in CAA-Fe medium were adjusted to an OD 600 of 1, diluted 10 -2 fold, mixed to yield the desired producer fraction, and further diluted to create Poisson distribution conditions. Producer/nonproducer co-cultures were started by inoculating the central 60 wells of two 96-well plates thereby adjusting the average initial cell numberN 0 to values between 2 and 25 cells/150 µL. Wells at the border of the plates were filled with water to minimize evaporation from central wells. For non-selective growth, co-cultures were grown in CAA-Fe medium shaking at 30°C for given periods of time. Due to the random distribution of initial cell number N 0 and producer faction x 0 in the 120 wells, each experiment was unique. An experiment was limited to 120 wells to allow initiation of the analysis of the subpopulations in the individual wells without uncontrolled changes of growth parameters during analysis. The experiment duration was set to 11h to allow evolution to act for a significant number of generations, while leaving bacteria in exponential growth phase (see SI)

Determination of growth parameters
Cell numbers N 0 and N (t) were determined by counting the colony forming units (cfu) of individual wells. For this purpose 100µL aliquots of the individual wells were plated on cetrimide [14] or King's B agar (contains per liter: 20 g peptone, 10 g glycerol, 1.

Simulation of growing populations
We performed simulations of 10080 wells using a Gillespie algorithm [17]. The initial numbers of "cells" per well were drawn at random from a Poisson distribution with a mean value ofN 0 measured in the corresponding experiment. The strain assigned to every individual in each well was determined by the outcome of a Bernoulli trial (i.e., coin-flip-like process) and the probability of assignment to strain A was set to the value ofx 0 measured in the experiment. After initialization, wells were grouped into 84 virtual 120-well "plates", and a random waiting time was selected for each well, drawn from an exponential distribution with the population size as parameter. The Gillespie algorithm was run until the average size across all wells matched the average size measured at the end of the growth experiments, or until a specified time had elapsed. pseudomonas entomophila l48 and its close relative pseudomonas putida KT2440. The distribution of compositions x first broadens due to demographic noise, but soon "freezes" to a steady state. The steady state form is maintained as long as the populations grow. Parameter values areN 0 = 10, x 0 = 0.33 (as for Fig. 2).

Exact calculations for steady-state composition distribution and moments.
Calculation of probability distribution.
Each population in the ensemble is initialized with A 0 individuals of type A and B 0 of type B. In the general case, A 0 and B 0 are independent random variables for each population with distributions P (A 0 ) and P (B 0 ). All populations evolve for ∆N reproduction events, of which a random amount ∆A generate new A-individuals. From the mathematical literature [35], it is well-known that ∆A follows a beta-binomial, with A 0 , B 0 and ∆N as parameters. The fraction of A-individuals x, then follows the probability where the sums run over all allowed values of their respective indices. P (∆A = k|A 0 , B 0 , ∆N ) is the probability of ∆A being equal to k, given the values of A 0 , B 0 and ∆N . The sum may easily be performed numerically. For the moments of the distribution there are, however, also closed-form analytic expressions.

Exact calculation of asymptotic moment values
Let · 0 be the average over the initial conditions, · ∆A be an average over ∆A, and · be an average over both quantities. From the properties of the beta-binomial distribution we know that, for a given initial condition, we have For the mean of x , one obtains Hence, the average composition is exactly conserved throughout the time evolution of the populations. In other words, the stochastic process is a martingale.
For the variance we obtain For long times (i.e., ∆N 1), ∆N + N 0 ∆N and (10) reduces to The argument up to here is completely independent of the particular choice of initial conditions. If the initial distribution is known, we may even make the value of the variance more explicit. In particular, consider the distribution we obtain from experiments: in each well, N 0 is Poisson-distributed with meanN 0 . Then one gets Within each well of (random) size N 0 there is an initial random number A 0 of A-individuals, which follows a Binomial distribution with parameters N 0 andx 0 . For this choice of distribution, it is possible that N 0 = 0, which would lead to an undetermined value of x 0 = A 0 /N 0 , and therefore also for the average x 0 . We can solve this problem by redefining x 0 : so that x 0 and its average have definite values, and x 0 0 =x 0 . With this we can compute the second moment of x 0 : The sum inside the braces can be solved using exponential and binomial series and yields The remaining series is an exponential integral (Ei), and therefore where we defined ϕ(N 0 ) := e −N 0 Ei(N 0 ) − γ − ln(N 0 ) . Then the variance of x reads For largeN 0 , through an asymptotic expansion [47], ϕ(N 0 ) can be approximated by To leading order inN 0 , then, the variance of x becomes in perfect agreement with our approximate results based on Master equations (Eq. (2) in main text, see also below).

Approximate calculations for the time evolution of the distribution moments
Using the Master Equation for the number of individuals of each strain (1), we are able to obtain the time evolution of the first three moments of the distribution of x. Equation (1) is sometimes called "Simple Growth Equation" and can be exactly solved (see, for example, [38]) using generating functions like To approximate the time evolution of the first three moments of x, however, we do not need the full solution, but only the first three moments of N A and N B . To this end, we insert the Master Equation (Eq. (1) in main text) in the definition of the generating function to get the time derivative for F (a, t): To obtain the time evolution of the nth moment, we apply the nth derivative with respect to a on both sides of equation (22), and solve for the corresponding moment. For the first three moments, the solution is K 1 , K 2 , K 3 are integration constants, which depend on the initial conditions. We consider the case of Poisson initial conditions. This means that the initial number of A is Poisson-distributed with mean valueN A,0 , and, since for the Poisson distribution the variance equals the mean, we get Employing these conditions in the solutions of the differential equations we found in Eq. (23) and (24), we get By the known properties of the Poisson distribution, the skewness of our initial distribution equals to For N B , the calculations are analogous. Note also that all calculations were exact so far. With the moments of N A and N B we can find the (approximate) time evolution of variance and skewness of x = N A /(N A + N B ). For the mean of x we have already seen in the exact calculation (see Eq. (6)) that it does not change with time, and hence its time evolution is already known.
To calculate the time evolution of the variance of x, we consider x as a function of N A and N B : From this we obtained Eq. (2) in main text. For infinite times the approximate result for the variance matches the exact one of Eq. (11). The skewness of the x distribution is calculated analogously: Comparison of initial and steady-state distributions of x, and entropy of the steady state distribution conditioned on the initial one. , the evolutionary fate of the population is largely determined by the initial population sampling. Therefore, the initial distribution (red bars) and the steady-state one (green bars) look qualitatively very similar. For intermediate values ofN 0 , however, population growth becomes more important, and the distributions look very different. The amount of composition values the population can access through growth can be quantified looking at the "unpredictability" of the steady-state composition, once the initial one is known: the more unpredictable, the more are made accessible by growth. Mathematically, the measure for this is called conditional entropy: the higher the entropy, the more unpredictable the outcome. Panel We simulate an ensemble of populations starting from Poisson initial conditions, and track their time evolution until the x distribution freezes. Once it freezes, we can build a joint histogram of initial and final compositions, which approximates the joint distribution P joint (x 0 , x f ). From P joint we can obtain the initial and final distributions as its marginal distributions, integrating over all values of x f and x 0 , respectively. The joint information (Shannon) entropy is defined as [48] The marginal entropies H(x 0 ) and H(x f ) are defined, analogously, through integrals only of P (x 0 ) over x 0 , and P (x f ) over x f , respectively. The conditional entropy of the final distribution given the initial is defined as It measures the amount of information necessary to describe the final distribution, once all information about the distribution of x 0 is known. Therefore, H(x f |x 0 ) provides a measure of how entropic (or "noisy") growth itself is [49]-or, in other words, how many different final compositions are possible given the initial condition. Figure 5(d) shows H(x f |x 0 ) from repeated simulations, all with the same initial distribution form, the samex 0 , but differentN 0 . For very smallN 0 (of the order of one or two individuals) the group formation almost completely determines the fate of populations: most populations start fixated, many with just a single founder individual, and the composition of each well remains the same during growth. The path followed by x in each population during time is a straight line, as the compositions stay constant. Therefore, x for different populations follow in time paths that do not cross or "mix". Growth produces very little demographic noise, and its conditional entropy tends to zero. For very largeN 0 (of the order of a few hundreds), the group sampling is again central to determine the final distribution. Very large populations, in fact, all start with similar compositions (according to the Law of Large Numbers), and their compositions are difficult to change, as each individual event has little impact. The composition distribution changes very little before freezing; time evolution paths of different populations "mix" very little. Entropy in this regime saturates for increasing initial sizes, and is rather low. Between the small size regime (where paths do not "mix") and the large size regime (where size limits "mixing"), we find a window where populations are small enough to significantly change their composition, but also large enough to not start fixated. This is the region where the conditional entropy peaks, and growth is the most important in determining the final distribution. Intuitively, the difference in variance between initial and final distribution could provide an alternative measure of the noise introduced by growth. However, of all x distributions between 0 and 1 with fixed x 0 , the one with maximal variance is the one for which x is only 0 or 1, i.e., when all populations start off fixated. In this case, the compositions never change during growth and the variance stays constant. Moreover, independently on the choice of initial distribution, the difference between initial and steady-state variance decreases for increasingN 0 (see Eq.(11)). Therefore, all considerations on noise sources based on variance would indicate that growth matters more when initial populations are smaller, in contrast with our observations. .While experiments for constant-sized populations of Drosophila observe significant fixations within the first tens of generations, we instead observe freezing of the probability distribution for the population composition, without any fixation. Figure 6: Growth curve of a mixed population. The population consists of pyoverdine producer (P. putida KT2440) and non-producer (P. putida 3E2) under non-selective (iron replete) conditions. Individual precultures of the strains were mixed and diluted in iron replete medium to yieldN 0 = 4 (in 150 µL), and x 0 = 0.5. Cells were grown aerobically at 30°C for 24 hours. The dots represent the mean N (t) of three independent replications, the bars the corresponding standard deviation. After a lag phase of about 2 hours, the cells start to grow exponentially and reach the stationary phase after about 14 h of growth. For the non-selective growth experiments used to test the predictions of the Pólya urn model, cells were grown for 11.5 h to ensure exponential growth conditions.  Table 2: Comparison of growth and pyoverdine production per cell of P. putida KT2440 and 3E2. Separate cultures of producer (P. putida KT2440) and non-producer (P. putida 3E2) were grown in ironlimiting (no addition of FeCl 3 ) and iron-replete medium (addition of 200 µM FeCl 3 ) at 30°C. The cell density was measured at 600 nm, and specific growth rates were calculated from density values of the exponential phase. The pyoverdine production was determined by fluorescence emission measurements (excitation 400 nm, emission at 460 nm). The pyoverdine production per cell represents the ratio of pyoverdine fluorescence and optical density measured after 24 h of growth. The values in the table are averages over a minimum of five experiments, with the corresponding standard deviation. The fluorescence value for the non-producing mutant in iron-limiting medium is 0 because the culture failed to grow.