Phenotypic delay in the evolution of bacterial antibiotic resistance: Mechanistic models and their implications

Phenotypic delay—the time delay between genetic mutation and expression of the corresponding phenotype—is generally neglected in evolutionary models, yet recent work suggests that it may be more common than previously assumed. Here, we use computer simulations and theory to investigate the significance of phenotypic delay for the evolution of bacterial resistance to antibiotics. We consider three mechanisms which could potentially cause phenotypic delay: effective polyploidy, dilution of antibiotic-sensitive molecules and accumulation of resistance-enhancing molecules. We find that the accumulation of resistant molecules is relevant only within a narrow parameter range, but both the dilution of sensitive molecules and effective polyploidy can cause phenotypic delay over a wide range of parameters. We further investigate whether these mechanisms could affect population survival under drug treatment and thereby explain observed discrepancies in mutation rates estimated by Luria-Delbrück fluctuation tests. While the effective polyploidy mechanism does not affect population survival, the dilution of sensitive molecules leads both to decreased probability of survival under drug treatment and underestimation of mutation rates in fluctuation tests. The dilution mechanism also changes the shape of the Luria-Delbrück distribution of mutant numbers, and we show that this modified distribution provides an improved fit to previously published experimental data.


1
Mutagenesis experiment: mathematical derivations Here we present derivations of various results discussed in the main text, for the dilution and accumulation mechanisms. The polyploidy mechanism is not discussed here as all results are presented in the main text.

Dilution of antibiotic-sensitive molecules
In this subsection all cells are assumed monoploid (c = 1). Suppose we start from a single genotypically resistant cell containing n sensitive molecules. We follow a single lineage, as in Fig 1a, and consider the cell in that lineage after g generations. Let the number of sensitive molecules in that cell be z n (g). A cell is to be considered phenotypically resistant if z n (g) = 0. During cell division, each of the n sensitive molecules may be lost (from the lineage we track) with probability 1/2. Therefore, the probability that any of the original sensitive molecules present in the initial cell remains in our chosen cell is 2 −g . Hence, z n (g) is a binomially distributed random variable with n trials and probability parameter 2 −g . The probability of phenotypic resistance in generation g is therefore The exponential approximation in Eq (S1) holds when n is large and 2 g ∝ n. Resistance thus emerges when n ≈ 2 g or for g ≈ log 2 (n), in agreement with the simple argument presented in the main text. We now turn to the probability of resistance occurring in the whole population (for which at least one cell must be resistant), for a population that is initiated with a single mutant cell with n sensitive molecules. Let p n (g) be the probability that no cell exists with zero sensitive molecules after g generations. Considering the distribution of molecules after the rst division we have This recursion allows us to calculate the probability of no resistance after g generations. If we start from x mutant cells rather than just a single mutant, the probability of no resistance after g generations is (p n (g)) x because mutant lineages are assumed to evolve independently of each other.
To obtain an approximate time for resistance to emerge, we consider the mean number of phenotypically resistant cells. Starting from x cells, after g generations the number of resistant cells can be written as where η i is a random variable equal to 1 if the ith cell has zero sensitive molecules, and 0 otherwise. The random variables (η i ) x2 g i=1 are dependent but have identical marginal distributions. Taking expectations and using E[η i ] = (1 − 2 −g ) n , we obtain the expected number of phenotypically resistant cells: Let τ denote the expected time to resistance in the population. From Eq (S2) the expected time until a resistant cell emerges can be expressed exactly as where p n (k) can be computed recursively from Eq (S2). In order to obtain an intuitively meaningful result, we reason that phenotypic delay corresponds to the time period during which the expected number of phenotypic mutants is less than 1. Thus a rough approximation for the expected value of τ is equal to one plus the number of the last generation in which the expected number of resistant cells is less than 1. We use the approximation of Eq (S4) and thus we wish to nd max{g : x2 g exp(−n2 −g ) < 1}. It can be veried by direct substitution that g = W (nx)−log(x) log (2) is the solution, where W is Lambert's function. Approximating W (z) ≈ log z − log log z for z → ∞, we obtain that E[τ ] ≈ 1 + log n − log log xn log(2) = 1 + log 2 (n/ log(nx)).

Accumulation of resistance-enhancing proteins
We consider the scenario that prior to cell division each genotypically resistant cells create M p resistance enhancing molecules. These are binomially distributed between the daughter cells at division. A cell is assumed to be phenotypically resistant once it has acquired M r sensitive molecules.
When tracking a single cell (or random lineage), let r g be the number of resistant molecules the cell possesses at generation g. From the model specication we have the stochastic recursion where d = denotes equality in distribution and we have abused notation using Bin(n, p) to denote independent binomial random variables with n trials with success probability p. Taking expectations over Eq (S7), and solving the resulting recursion, leads to E[r g ] = M p (1−2 −g ). In fact the full distribution of r g can be obtained by observing that it is distributionally equal to corresponding to the molecular contributions to our chosen cell from each of the previous generations. This sum has a Poisson-binomial distribution, however the cumulative distribution function (of interest as we care about Pr(r g ≥ M r )) is relatively uninformative and numerically unstable. Therefore we simply note that the number of molecules after many generations is M p on average, with uctuations around this value. The variance Var(r g ) = M p (2/3 − 2 −g + 4 −g /3) tends to 2M p /3 in the limit g → ∞, hence uctuations become less important for large M p (the coecient of variation tends to zero).
As a rst approximation, since our condition for resistance is that E[r g ] ≥ M r , then resistance will occur at τ ≈ − log 2 (1 − 1/m) as long as m = M p /M r > 1. Note that if m < 1, r g may stray above M r due to uctuations, but this cannot produce sustained resistance as such high r g values will be transient. We therefore do not consider the case of m < 1 in detail. We also omit the special case m = 1 for lack of biological realism.
Turning to the population as whole, but still starting with a single cell, we rstly note that, as with Eq (S2), the following recursion holds for the probability of no resistance after g generations starting with n resistance molecules: Initial conditions for the recursion are p n (0) = 1 for n < M r , and 0 otherwise. Again, for x > 1, the probability of no resistance is (p n (g)) x . In summary, for m < 1 resistance will not be stably achieved (steady state less than M r ), while if we seek a delay of at least a generation (τ ≥ 1) we require m ≤ 2. A noticeable phenotypic delay occurs thus only in the narrow parameter range 1 ≤ m ≤ 2. As mentioned in the main text, resistance will occur faster in the whole population than down any lineage, further narrowing the parameter regime in that setting. 1.3 Combining the dilution and polyploidy mechanism: phenotypic penetrance We now consider bacteria with ploidy c. Starting from a cell possessing a mutated allele on a single chromosome, g c = log 2 (c) divisions are required to generate a cell with all c copies having the resistant allele (henceforth termed chromosomal resistant). To include the dilution mechanism, suppose n sensitive molecules exist in the initial cell and that each wild-type chromosome produces n/c sensitive molecules between cell divisions. A descendant of the initial mutated cell that emerges with a full suite of resistant chromosomes will eventually have phenotypically resistant cells amongst their progeny (after dilution of any sensitive molecules). This cell will initially have n c sensitive molecules. Binomial partitioning of the n original sensitive molecules and those created between the appearance of the rst mutant chromosome and the genotypically resistant cell then allows us to write where Bin(. . . ) denotes, as before, independent binomial random numbers. Note that here we restrict ourselves to n such that (1 − 2 i /c)n is an integer, for each 0 ≤ i ≤ g c − 1. This is satised if each sensitive chromosome produces an integer number of sensitive molecules (n/c). Any phenotypically resistant cells will be descendants of the initial genotypically resistant cell. The question of whether there is any such resistant cell by generation g is therefore equivalent to asking whether there is any resistant cell by generation g − g c but initiating the process with the initial genotypically resistant bacteria with n c sensitive molecules. To nd the expected number of phenotypically resistant cells, we average over (S4) with respect to the distribution of n c , that is we compute Here, the generating function for the binomial distribution has been used. If, as before, we dene genotypically resistant cells as having at least one resistant chromosome, then for any g ≤ g c only a single genotypically resistant bacterium will exist because we have assumed co-inheritance of recently linked chromosomes. After generation g c all genotypically resistant cells are the descendants of the initial genotypically resistant cell, and their number is 2 g−gc . Hence, of the genotypically resistant cells, the expected fraction of phenotypically resistant cells is given by which gives Eq (1) from the main text. Following [1] we call this quantity the phenotypic penetrance. The limiting case of n = 0 gives a Heaviside step function θ(g − g c ) as expected. Note that while this equation was derived by considering the scenario where we start with a single mutated cell, it holds also when there are initially many mutated cells (x > 1). This is because phenotypic penetrance is the ratio of phenotypically resistant cells to genotypically resistant cells.

Single-cell simulations
One of the key arguments used by Sun et al. [1] to arm that phenotypic delay is caused mostly by the eective polyploidy mechanism was the asymmetric inheritance of resistance, i.e., the observation that some lineages become resistant while others do not. This eect can only be observed by studying the behaviour of individual cells. Hence, to understand the eect of the dierent phenotypic delay mechanisms, we simulate a single cell immediately after a mutation (Fig 1a). When the cell divides, we follow one randomly selected daughter and continue tracking that cell until phenotypic resistance emerges or the cell loses all resistant alleles. We repeat the simulation 10000 times and calculate the probability of the cell becoming phenotypically resistant as a function of the number of generations from the mutation.
In the case of the dilution mechanism, the phenotypic delay is controlled by the number of molecules n that need to be diluted out (Fig 1b). More molecules lead to a longer delay because more generations must pass before random segregation produces a cell without any sensitive molecules. In contrast to the eective polyploidy model, the probability of phenotypic resistance always approaches 1 for long times, and the approach to this end point is also more gradual. This is because, unlike in the polyploidy model, molecules segregate independently of each other, and the probability of producing a cell with only resistant molecules is non zero (albeit very small) already after the rst division. Since genetic resistance cannot be lost in this model (there is no chromosomal segregation), all cells eventually become resistant. Analytically, we showed (Eq S1) that the probability of resistance is equal to : where g is the number of generations from the mutation. In the case of the eective polyploidy mechanism, the probability of resistance is exactly zero until generation 1 + log 2 c, after which it assumes a constant value 1/c (Fig 1c), where c is the eective polyploidy. This is because if only one of the c = 1, 2, 4, . . . chromosomes has the resistant allele, exactly one cell out of c progeny of the initial cells ends up with resistant chromosomes, and the rest c − 1 cells have only sensitive chromosomes. Therefore, (S14) In Fig 1c, we observe close agreement between the theory and our stochastic simulations. The third mechanism that could lead to a phenotypic delay is the accumulation of resistanceenhancing molecules (Fig 1d). As explained in the main text and in Section 1.2, a substantial phenotypic delay is only observed within a narrow parameter range for this mechanism: 1 m 2.

3
Biasing the partitioning of molecules at divisions shortens phenotypic delay In the main text, for the dilution mechanism we assumed that target molecules are divided between the two daughters without any bias. While this is a reasonable assumption for cytoplasmic proteins, there is evidence that membrane-associated proteins are segregated in a biased way: a fraction p > 1/2 ends in a daughter cell with the older pole [2,3,4]. In particular, a bias p = 0.62 has been shown for eux pumps [5], whose expression is associated with elevated resistance to antibiotics. Another example is OmpA (outer membrane protein A) which also tends to accumulate in cells with older poles [3]. OmpA is implicated in phage invasion; mutations that decrease or alter the protein create resistance to bacteriophages [6,7].
To investigate the eect of biased protein segregation on the length of phenotypic delay, we repeated our single-cell and population level simulations using the dilution mechanism (n = 1000) with a biased binomial distribution (p = 0.62 instead of p = 0.5 as used in the main text). A very small eect of the bias is present at the single-cell level (Fig 3a), but a much larger eect (delay shortened by one generation) is observed at the population level (Fig 3b). This can be explained as follows. At the population level, it is enough for a single cell to become resistant. Biasing causes cells with the younger pole to be depleted of the sensitive molecules faster and, hence, to become resistant earlier. On the other hand, at the single-cell level we follow a randomly selected daughter cell and not necessary the one with the younger pole. This largely nullies the eect because even though some cells in the lineage get rid of the sensitive molecules faster than in the case of no bias, others take longer. 4 Phenotypic delay when the number of target molecules is independent of the growth rate In the main text, for the combined dilution and eective polyploidy mechanisms we assumed that the number of sensitive target molecules depends on the growth rate. In particular, we assumed that the gyrase enzyme constitutes a roughly constant fraction of the proteome, independently of the growth rate [8]. Assuming that the total protein mass is proportional to cell volume, and that volume ∝ 2 λ/λ0 , where λ is the growth rate and λ 0 = 1h −1 [9,10], we can then predict the change in the number of gyrase molecules per cell associated with a change in growth rate. For example, if we assume that n = 20 for t d = 60 min, then for t d = 30 min we obtain n = 40 (Fig 4a).
In contrast, in this subsection we investigate what happens if the number of target molecules is independent of the doubling time. To be specic, let us take n = 20 for both t d = 60 min and t d = 30 min (Fig 4b). We observe that when the number of target molecules n does not depend on t d , the increase in phenotypic delay decreases from 2 to 1 generation, similar to the results obtained for the eective polyploidy model without the dilution mechanism. Hence, the increase in phenotypic delay in the combined model caused by a decrease of the doubling time t d is caused largely by the increase in the number of target molecules.
In the main text, we also show that the probability of surviving an antibiotic treatment in a simulated infection for the combined model is a function of t d (Fig 4c). In Fig 4d we

5
Phenotypic delay for the dilution mechanism with non-zero molecular threshold for resistance In the main text, for the dilution mechanism, we assumed that a cell becomes resistant only when it loses all n sensitive molecules. However, in reality a cell may be resistant even if a few sensitive molecules are left. Here we repeat the analysis for the dilution mechanism for which the threshold in the number of sensitive molecules for which a cell is phenotypically resistant is n r > 0. The results are summarized in Fig 5. As expected, the phenotypic delay decreases as n r increases, because the number of molecules required to be diluted out for resistance to emerge decreases. Interestingly, we also observe that increasing n r leads to a less gradual appearance of resistance at the single-cell level. 6 Maximum likelihood estimate of the mutation rate In Fig 6 we show the maximum likelihood estimates for the mutation probability, obtained using the package an [11], from 1000 simulations of the uctuation test from Ref. [12]. Simulations were performed using the same method as discussed in the main text, Section 4.6. Each of the 1000 simulated experiments contained 40 independent replicates. Replicates were initiated with N 0 = 100 cells and the nal population size was N f = 2 20 N 0 so that the population of bacteria grew for 20 generations. A relatively narrow range of estimated mutation probabilities was obtained, as shown in Fig 6. This strongly suggests that the observed discrepancy between the mutation rates obtained from sequencing and the uctuation test is not due to sampling variation but is likely to have another (biological) cause. We suggest that this may be phenotypic delay.

7
Sensitivity analysis of the ABC method of model comparison In the main text we report that the probability that the data of Boe et al. [13] is generated by the phenotypic delay model with dilution, compared to the model with no delay, is 0.97.
To check how robust this result is, we performed two sensitivity analyses, altering the size of our simulation bank N sim , and the number of`closest' simulations we compare to, N thresh (N thresh = 100 in the main text). Firstly, we query whether enough simulations have been performed. To do so we sub-sample our bank of simulations. Secondly, we assess how our probability estimate changes as a function of N thresh . Note that as N thresh → N sim the probability estimate tends to 0.5 as an equal proportion of simulations comes from either model. The results are presented in Fig 7, which supports our conclusion that the model with delay via the dilution mechanism is indeed a better generative model for the data than the model without delay.