Phenotypic delay: mechanistic models and their implications for the evolution of resistance to antibiotics

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-Delbruck 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-Delbruck distribution of mutant numbers, and we show this modified distribution provides an improved explanation of previously published experimental data.


Introduction
The emergence of resistance to drugs is a significant problem in the treatment of diseases such as cancer [1], and viral [2] and bacterial infections [3]. In the case of infectious diseases, treatment failure can often be attributed to the infecting organism being resistant [4]. However, resistance can also be caused by new genetic mutations in infections with sufficiently high pathogen load [5]. Such de novo resistance has been observed in infections caused by Staphylococcus aureus during endocarditis [6,7], Burkholderia dolosa [8,5] and Pseudomonas aeruginosa [9,10] during cystic fibrosis, Escherichia coli during asymptomatic bacteriuria [11], and Helicobacter pylori [12,13].
The emergence and spreading of resistant variants in populations of pathogenic cells has received much experimental [14,15,16,17,18] and theoretical attention [19,20,21,22]. However, most mathematical models assume that a genetic mutation immediately transforms a sensitive cell into a resistant cell [23,24,25,26,27,28]. In reality, a new allele (genetic variant) must be expressed to a sufficient level before the cell becomes phenotypically resistant. The time between the occurrence of a genetic mutation and its phenotypic expression is called phenotypic delay and is also referred to as delayed phenotypic expression, phenotypic lag, cytoplasmic lag or phenomic lag.
Phenotypic delay was first observed in 1934 by Sonnenborn and Lynch [29] when studying the effect of conjugation on the fission rate of Paramecium aurelia. It was observed that, even though the fission rate is determined genetically, the new hybrids initially assumed the fission rates of the parent with whom they shared their cytoplasm, and the new genotype manifested phenotypically only after a few generations. Phenotypic delays were studied during the 1940s and 1950s, both theoretically [30] and experimentally [31,32]. Interestingly, in their hallmark work on the randomness of mutations in bacteria [33], Luria and Delbrück discussed the possible effect of a phenotypic delay on the estimation of mutation rates. The interest in phenotypic delays waned for the next fifty years, mostly because experimental data failed to reveal their presence [34,33]. However, Sun et al. [35] recently demonstrated the existence of a phenotypic delay of 3-4 generations in the evolution of resistance of Escherichia coli to antibiotics rifampicin, nalidixic acid and streptomycin. Sun et al. attributed this delay to effective polyploidy (see below) and argued that not accounting for a phenotypic delay can cause an underestimation of the mutation rate in fluctuation tests [33,36].
Here, we generalize these observations and investigate other reasons that may lead to phenotypic delay. We consider three mechanisms: (i) effective polyploidy as in [35], (ii) the dilution of sensitive molecules targeted by the drug, and (iii) the accumulation of resistance-enhancing molecules.
Effective polyploidy refers to the fact that a single cell can contain multiple copies of a given gene. This can be due to gene duplication events or carriage of multicopy plasmids; it also occurs in fast-growing bacteria, which initiate new rounds of DNA replication before the previous round has finished, allowing for a shorter generation time than the time needed to replicate the chromosome [37,38,39]. Since a de novo resistance mutation happens in only one of the multiple gene copies, it may take several generations before a cell emerges in which all gene copies contain the mutated allele. Until then, sensitive and resistant variants of the target protein coexist in the cell. A phenotypic delay occurs when the resistance mutation is recessive, i.e., the sensitive variant must be replaced by the resistant variant for the cell to become resistant. This is the case for antibiotics which form toxic adducts with their targets [40,41]. Examples are quinolones that lock the enzyme DNA gyrase onto the DNA and prevent DNA replication [42], and polymixins that bind to lipids in the outer membrane which causes membrane perforation [43,44].
The dilution mechanism also assumes the mutation to be recessive, but in contrast to the polyploidy mechanism it focuses on the removal of the sensitive target protein through the process of cell growth and division. As a mutated cell grows, resistant version of the protein accumulates; a subsequent division creates two cells in which the fraction of the sensitive variant is less than in the parent cell. A single copy of the gene may therefore still cause a considerable delay if the number of sensitive proteins to dilute out is large before resistance can be established.
The accumulation mechanism posits that sufficiently many resistant variants of the protein must be produced to cause resistance. This applies to mutations that enhance the expression of drug efflux pumps [45], β-lactamase enzymes that hydrolyze β-lactam antibiotics [46,47], or mutations that protect ribosomes from tetracycline [48] which restores the active ribosome pool [49]. In these cases, a phenotypic delay could emerge due to the time required for the resistance-enhancing molecule to accumulate in the cell to a level high enough to cause resistance.
We first analyse the three mechanisms using computer simulations and analytic calculations. We find that the accumulation of resistant-enhancing molecules only leads to phenotypic delay within a limited parameter range, while effective polyploidy and the dilution of sensitive molecules lead to phenotypic delay for a broad range of parameters. We then show that while the effective polyploidy mechanism does not affect the probability that a population survives antibiotic challenge, dilution of sensitive protein leads to decreased probability of survival under drug treatment.
We then investigate the possibility of detecting a phenotypic delay experimentally. We first show that the dilution mechanism would cause an underestimation of mutation rates in fluctuation tests compared to the true genetic rate of mutations. In a fluctuation test, one measures the distribution of mutant numbers in replicate populations that have been allowed to grow and evolve for a fixed number of generations. The mutation rate is then estimated by fitting a population dynamics model to the experimental distribution [33,50]. In agreement with our prediction, the mutation rate of Escherichia coli obtained in fluctuation tests has been found to be an order of magnitude smaller than the rate obtained by DNA sequencing [51].
We also show that the dilution mechanism subtly alters the shape of the Luria-Delbrück distri-bution of mutant numbers. Discrepancies between the shapes of the experimental and theoretically predicted mutant number distributions have been observed since the original experiments of Luria and Delbrück [33,34,30,52], but have never been satisfactorily explained. Using the experimental data set from Boe et al. [52] for fluoroquinolone antibiotics, we show that a mathematical model that includes the dilution mechanism fits the data better than the no-delay Luria-Delbrück model, thus providing indirect evidence for the existence of this type of phenotypic delay in de novo evolution of resistance to fluoroquinolones.

Modeling the emergence of phenotypic delay
To explore the characteristic features of the three different phenotypic delay mechanisms -dilution of sensitive molecules, effective polyploidy, and accumulation of the resistant variant -we first simulate an idealised mutagenesis experiment (figure 1a). A population of sensitive bacteria is exposed to a mutagen (e.g., UV radiation [53,54]) which instantaneously induces mutations in a small fraction of the cells [55,56], thus providing a well-defined reference point from which we count time. Cells immediately begin to express the mutated allele, but because of the existence of phenotypic delay, they remain sensitive to the antibiotic for some time; phenotypically resistant cells emerge only after a few generations. We shall investigate the emergence of resistance in two different ways. The first way is to follow a random lineage, starting from a single mutant bacterium (i.e., at each division we follow one of the randomly selected daughter cells) and we measure the waiting time before a phenotypically resistant cell emerges in that lineage (SI figure S1). The second way is to track the entire population post mutagenesis, and examine the waiting time before the first phenotypically resistant cell emerges in the whole population (figure 1).
Dilution of antibiotic-sensitive molecules: If the resistance mutation is recessive, such that a small number of sensitive target molecules are enough to cause antibiotic sensitivity, a phenotypic delay can arise from the time taken to replace sensitive target molecules by resistant ones. To model this, we assume that each cell has a fixed number n of target molecules, that can be either sensitive or resistant. Initially, all n molecules are sensitive but once a mutation has happened, the cell starts to produce resistant molecules. We suppose that, upon cell division, the n molecules are partitioned stochastically without bias between the two daughter cells, and new molecules are produced to bring the total back to n per cell (figure 1b). For simplicity, cells are considered phenotypically resistant only when they contain no sensitive molecules.
For this mechanism, the length of the phenotypic delay increases approximately logarithmically with the number n of sensitive molecules that need to be diluted for resistance to emerge (figure 1c). To understand this, suppose momentarily that n is a power of 2 and stochasticity can be neglected so that each daughter cell receives exactly half the number of molecules of the parent cell. Then for any lineage stemming from a genotypically mutant cell, the number of inherited sensitive molcules will be 2 n−1 , 2 n−2 , . . ., as the generations progress. After log 2 n generations all cells will have a single sensitive molecule and hence the first phenotypically resistant cell will then emerge after 1 + log 2 n generations. In this deterministic setting this will also be the time for the population to become resistant.
In the more realistic case of stochastic segregation of molecules, the probability of resistance along a random lineage after g generations is approximately exp(−n2 −g ) (see SI Section 1.1). Hence the probability of resistance emerging in a lineage is negligible until generation g set by 2 g ≈ n, when the probability rapidly rises to 1. Therefore, in line with our deterministic reasoning, resistance along a random lineage will emerge after g ≈ log 2 n generations. Interestingly, however, we obtain a different result for the probability that the population produces at least one resistant cell. If we start from x genotypically mutant cells in the population, the first phenotypically resistant cell in the population emerges, on average, after an approximate time 1 + log 2 (n/ log(xn)) (SI Section 1.1). We can also calculate the resistance probability through a recursion relation (SI Section 1.1); the results fully reproduce the simulation (figure 1c.). The emergence of resistance at the population level is thus accelerated compared to what one would obtain based on deterministic dilution. An extended model in which molecules are distributed in a biased way between the two The copyright holder for this preprint (which was not peer-reviewed) is the . daughter cells, inspired by recent evidence on accumulation of membrane proteins in the daughter cell with the older pole [57,58,59,60], leads to a similar result (SI Section 3). However, the bias decreases slightly the phenotypic delay at a population level (SI figure S3). This is because the bias creates lineages which will be low in the number of resistant molecules.
Effective polyploidy: Rapidly dividing bacteria can become effectively polyploid when they initiate a round of DNA replication before the previous round has finished; this leads to the presence of multiple copies of at least some parts of the chromosome [37] (figure 1d). Crucially, the degree of polyploidy (number of gene copies) depends on the bacterial growth rate, as well as on other factors such as the genetic locus. To model phenotypic delay caused by effective polyploidy, we assume that each cell has a given number c of chromosome copies that is growth-rate dependent according to the well-established Cooper-Helmstetter model of E.coli chromosome replication [37] (Methods). Each chromosome copy contains a single allele, encoding the antibiotic target, that can be either sensitive or resistant. Initially all chromosomes have the sensitive allele but a mutation changes one allele from sensitive to resistant. We then simulate the process of DNA replication and cell division, taking account of the fact that duplicated resistant alleles are co-inherited -for example, if a cell has two chromosome copies, one with a resistant allele and the other with a sensitive allele, then upon replication and division, one daughter cell will have two sensitive alleles and the other daughter cell will have two resistant alleles [38] (Methods). We assume that a cell becomes phenotypically resistant when none of its chromosomes contain the sensitive allele (i.e., the resistant allele is assumed to be recessive). In this model, the waiting time until a cell acquires a full suite of resistant chromosomes, i.e. the phenotypic delay, is log 2 c generations (figure 1e). This delay time is the same whether we track a given lineage or the entire population (since it is deterministic). However, resistance will not occur in all lineages; of the c lineages descended from the original mutant cell, resistance will eventually occur in only one of them [35] (SI figure S1).
We note that effective polyploidy genrally causes a shorter delay than dilution of sensitive molecules: 2 to 3 generations for rapidly growing bacteria (c = 4 or 8 [37,35]), versus 5 generations for the dilution mechanism (assuming n ≈ 500 typical for the gyrase enzyme targeted by fluoroquinolones [61,62]). Transition in the probability of resistance is also sharper than for the dilution mechanism in which stochasticity of the segregation process smoothens out the transition (compare figures 1a and e). Finally, for effective polyploidy, we expect only one in every c lineages to become resistant, while for dilution of sensitive molecules, effectively all lineages will eventually become resistant.
Accumulation of resistance-enhancing molecules: Phenotypic delay can also emerge due to the time needed to accumulate resistance-enhancing molecules to a sufficiently high level (figure 1f). To model this mechanism, we suppose that during each cell cycle a genotypically resistant cell produces M p resistance-enhancing molecules, which are randomly distributed between daughter cells at division. A cell becomes resistant when it has M r or more resistance-enhancing molecules. Interestingly, considering either a single lineage (SI figure S1) or the entire population (figure 1g), we find that phenotypic delay emerges only within a limited parameter range: 1 m 2, where m =

Mp
Mr is the ratio of the number of molecules produced during a cell cycle and the number of molecules needed for resistance. The origin of this limited parameter range is most easily explained by considering a single lineage. Tracking a lineage arising from a single mutant cell, the cell in the gth generation will be born with an average of M p (1 − 2 −g ) molecules (SI Section 1.2). The steady-state number of molecules (found by taking g → ∞) is M p . Thus if m < 1, the steady state number of molecules will be always smaller than the minimum required number M r , and the lineage will never become phenotypically resistant. Conversely, if m > 1, phenotypic resistance will emerge after approximately τ = − log 2 (1 − 1/m) generations when the average number of resistance-enhancing molecules exceeds M r . But for the delay to be detectable -at least one generation (τ ≥ 1) -we require m ≤ 2. Considering now the scenario where we track the entire population, we again expect the steady-state molecule number M p to be rapidly reached for all cells, so that there will be no phenotypic resistance for m < 1. Further, if resistance does emerge (for m > 1), it will do so more quickly in the entire population than along the random lineage (as resistance may be acquired in any lineage). We thus expect an even tighter upper bound on the value of m for phenotypic delay to manifest itself on the population level.
Since our analysis shows that, for this mechanism, phenotypic delay only emerges in a narrow parameter range, we conclude that the accumulation of resistance-enhancing molecules is unlikely to be biologically relevant in causing phenotypic delay. Therefore we do not explore this mechanism further.

Combining effective polyploidy and dilution
In reality, for a recessive resistance mutation, we expect both the effective polyploidy and dilution mechanisms to contribute to the phenotypic delay. To see this, we simulated a model combining the two mechanisms, tracking the emergence of resistance at a single-cell and population level. Our simulations predict a phenotypic delay with characteristics of both mechanisms (figure 2).
Focusing first on a single lineage ( figure 2a,b), we observe that the long-term probability of phenotypic resistance depends on the ploidy c, tending approximately to 1/c, as expected for the effective polyploidy mechanism, while the approach to this value is gradual as expected for the dilution mechanism. Combining both mechanisms increases the length of the delay compared to either mechanism acting in isolation. Following Sun et al. [35], we also calculate the phenotypic penetrance defined as the proportion of genetic mutants which are phenotypically resistant, in the entire population. The expected phenotypic penetrance is (see SI for derivation Section 1.3): Note that n = 0 corresponds to only the effective polyploidy mechanism, while c = 1 corresponds to only the dilution mechanism being present. The piecewise form of Eq. (1) arises because no cell can become phenotypically resistant until all its chromosomes have the resistant allele. Figure  2d shows that the phenotypic penetrance predicted by Eq. (1) increases gradually with time (characteristic of the dilution mechanism) but with a delay determined by effective polyploidy. We now return to computer simulations to study the emergence of resistance on the population level following mutagenesis (figure 1). In general, both the ploidy c and the number of antibiotic target molecules per cell n will depend on the doubling time t d (or growth rate) of cells. To be more specific, we consider resistance of E. coli to fluoroquinolone antibiotics that arises through mutations in the DNA gyrase (protein targeted by the antibiotic). Gyrase abundance as a fraction of the proteome has been found to be independent of the growth rate [63]. We therefore assume that the number n of gyrases per cell is proportional to the cell volume V . We model the volume as V ∝ 2 λ/λ0 , where λ = (ln 2)/t d is the growth rate and λ 0 = 1h −1 [64,65,66,67], and polyploidy using the Cooper-Helmstetter model [37] (see Methods and Model for details). We find for slowgrowing cells (t d = 60 min) that c = 2 and n = 20, while for fast-growing cells (t d = 30 min), c = 4 and n = 40. Note that here we do not assume realistic values of n because the minimum number n r of poisoned sensitive gyrase molecules required to inhibit growth is probably much higher than n r = 1 assumed in the model. n should be therefore interpreted more correctly as the number of "units" of gyrase, with one unit equivalent to n r molecules. Figure 2c shows that the phenotypic delay is longer for the fast-growing population, and that this is mostly caused by the increase in the number of molecules n (SI figure S4). We also observe that protein dilution leads to a smoother transition between sensitivity and resistance than the transition due to effective polyploidy alone.

Dilution mechanism, but not effective polyploidy, affects the probability of clearing an infection
To understand better the practical significance of phenotypic delay, we simulated antibiotic treatment of an idealised bacterial infection (figure 3). We assume that, before treatment, the population of bacteria grows exponentially in discrete generations, and cells mutate with probability µ = 10 −7 per cell per replication. When the population size reaches 10 7 , an antibiotic is introduced; this causes all phenotypically sensitive bacteria to die, leaving only phenotypically resistant cells (figure 3b). We are interested in the probability that the bacterial infection survives the antibiotic treatment, a concept closely related to evolutionary rescue probability, i.e., the probability that cells can survive a sudden environmental change thanks to an adaptive mutation [68,35,69]. Since sensitive cells do not reproduce in our simulations in the presence of the antibiotic, survival can only be due to pre-existing mutations (standing genetic variation).
6 . CC-BY-NC-ND 4.0 International license author/funder. It is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the . Figure 2: Joint effect of the dilution and effective polyploidy mechanisms. (a-b) Probability of resistance of a single mutated cell. While the long-term probability is defined by the effective polyploidy, short-term behaviour is determined by the dilution mechanism, leading to longer phenotypic delays than the effective polyploidy mechanism would do. (c) Population-level probability of resistance versus the number of generations from the mutation event, for n = 20. The combined mechanism leads to smoother curves than the effective polyploidy mechanism and longer delays than for both mechanism individually. (d) Phenotypic penetrance (ratio of phenotypically resistant to genotypically resistant cells, obtained from Eq. (1)) for the different mechanisms, for c = 8, n = 8. The dashed red line indicates when the phenotypic penetrance surpasses 1/2, which is the threshold used by Sun et al. [35] to define the emergence of phenotypic resistance. With this definition, the dilution mechanism plus effective polyploidy doubles the delay (generation 6 as opposed to generation 3).   We first consider the effective polyploidy model, with ploidy c controlled by the doubling time t d . In agreement with Sun et al. [35], we find that t d has no effect on the survival probability ( figure  3c). This is due to a cancellation of two effects: the increased number of gene copies increases the per-cell chance of genetic mutation, but also increases the length of the phenotypic delay. In contrast, phenotypic delay caused by the dilution of sensitive molecules does affect the survival probability (figure 3d). The survival probability strongly depends on n, and decreases significantly from 0.69 for n = 0 to 0.06 for n = 100.
We also simulated a mixed case where both the effective polyploidy and dilution mechanisms combined, with ploidy c and molecule number n determined by the doubling time t d as described in Sec. 2.2. In this case the survival probability does depend on the doubling time (figure 3e; blue line). This is mostly caused by the change in the molecular number n as a function of doubling time. If we neglect the dependence of n on t d , the effect is much smaller, although there is still some dependence on t d because the rate of resistant protein production depends on the resistant gene copy number, which increases en route to the full suite of resistant chromosomes (SI figure  S4).

Phenotypic delay due to dilution changes the Luria-Delbrück distribution and biases mutation rate estimates
The scenario discussed in the previous section is equivalent to the Luria-Delbrück fluctuation test [33,70], which has been extensively studied theoretically [71,72,73,50,74,75,76,77,78]. A small number of sensitive bacteria are allowed to grow until the population reaches a certain size. The cells are then plated on a selective medium (often an antibiotic) to reveal the number of mutated bacteria in the population. The distribution of the number of mutants (measured over replicate experiments) is termed the Luria-Delbrück distribution. This distribution has a power-law tail caused by mutational "jackpot" events [33,70,77] in which rare, early-occurring mutants produce many descendants in the population. The fluctuation test and the corresponding mathematical model have been used to estimate mutation rates in bacteria. Here, we discuss the effect of phenotypic delay on the Luria-Delbrück distribution and on the resulting mutation rate estimate. First, we note that phenotypic delay caused by effective polyploidy alone does not affect the Luria-Delbrück distribution. As discussed in the previous section, this is due to an exact cancellation of two effects: increased ploidy leads to more mutations per bacterium but also a longer phenotypic delay [35]. In contrast, the dilution model does alter the Luria-Delbrück distribution. Figure 4a shows that, for a fixed mutation probability µ and initial and final population sizes, phe- The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2019.12.19.883132 doi: bioRxiv preprint due to jackpot events, is unaffected by the delay -this is because mutants that arise early will have had sufficient time to dilute out all sensitive molecules before being exposed to the antibiotic.
From a practical point of view, the mutation probability is often unknown and the fluctuation test is used to estimate it. To investigate the effect of phenotypic delay on the measured mutation probability, we simulated the fluctuation test for the dilution model with n = 16, for a range of mutation probabilities. We compared the resulting mutant number distributions to that obtained in an equivalent simulation without phenotypic delay, with mutation probability µ = 10 −7 . Using a genetic algorithm [79] to minimize the L 2 norm between the distributions with and without phenotypic delay, we found that the phenotypic delay model required a much larger mutation probability (µ = 8 × 10 −7 ) to reproduce the distribution of the no-delay model. This suggests that neglecting phenotypic delay when fitting theory to fluctuation test data could significantly underestimate the true mutation probability. We also note that the "closest match" distributions with and without phenotypic delay are not exactly identical (figure 4b). The model with phenotypic delay leads to a larger number of jackpot events (as might be expected since the mutation probability is higher) and a reduced number of replicates with few mutants, consistent with suppression of late-occuring mutants by the phenotypic delay.
Our result explains the apparent discrepancy between mutations probabilities estimated through different methods. Lee et al. measured the mutation probability of E. coli using both fluctuation tests (fluoroquinolone nalidixic acid as selective agent) and whole-genome sequencing [51]. The fluctuation test underestimated the mutation probability by a factor of 9.5; Lee et al. suggested that this could be caused by phenotypic delay [51]. To see whether our dilution model could explain this, we simulated the 40-replicate, 20 generation fluctuation test experiment of Lee et al. [51], using the mutation probability as estimated by whole-genome sequencing (µ = 3.98 × 10 −9 , total for all mutations givins sufficient resistance to nalidixic acid), for differing values of the number n of target "units" (gyrase molecules). For each n we simulated 1000 realisations of the 40-replicate experiment, and for each realisation we estimated the mutation probability under the no-delay model using the maximum likelihood method [50] (the same as used by Lee et al.) implemented in the package flan [80]. This procedure correctly reproduced the mutation probability for data from simulations without delay (n = 0; SI Fig. S6). For the model with delay, the maximum likelihood fit returned a mutation probability that was lower than the true one (figure 5a); the discrepancy increased with the phenotypic delay. To obtain an apparent mutation probability that is underestimated by a factor of 9.5, as observed by Lee et al. [51], we require n ≈ 30; i.e. roughly 30 sensitive 'units' of the antibiotic target must be diluted out before a cell becomes phenotypically resistant. Thus, while our simulations do not prove that phenotypic delay is responsible for the discrepancy observed by Lee et al., they suggest that it is a plausible explanation.

Mutant number distributions may support the existence of phenotypic delay
Our results suggest that a phenotypic delay caused by dilution produces a characteristic (though small) change in the shape of the observed mutant number distribution (figure 4b Qualitatively, this seems to be consistent with our expectations for the dilution model ( figure 4a). To see if the dilution model of phenotypic delay indeed provides a superior fit to Boe et al's data, we used an approximate Bayesian computation (ABC) approach [81] (Methods). We simulated a 1104-replicate fluctuation experiment 10 4 times, for both the model with and without delay, with initial and final population sizes of 1.2 × 10 4 and 1.2 × 10 9 matching those of Boe et al. [52]. We then determined the posterior Bayesian probability that the experimental data is generated by the We simulated the fluctuation experiment of Ref. [51], where the authors report a factor of 9.5 difference between µ obtained by DNA sequencing and fluctuation tests. For each n we simulated 1000 experiments with the sequencing-derived mutation probability µ = 3.98 × 10 −9 and then used the same estimation procedure as Ref. [51] to infer µ assuming no delay exists. n = 30 sensitive molecules are required to account for the discrepancy observed. Error bars are 1.96 × standard error. (b) The experimental cumulative mutant frequency distribution reported by Boe et al. [52] (black points) and the best-fit simulated distribution (green line) for the dilution phenotypic delay model. (c) Histograms of the probability of the delay model obtained by applying the approximate Bayesian computation scheme to simulated data. Our classification algorithm correctly discriminates between the models.

11
. CC-BY-NC-ND 4.0 International license author/funder. It is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2019.12.19.883132 doi: bioRxiv preprint delay model as opposed to the no-delay model. We also examined the validity of our approach on synthetic data, see figure 5c and Methods and models for details. We obtain that the probability of the experimental data coming from the model with phenotypic delay is 0.97. We thus conclude that the Boe et al. data supports the existence of phenotypic delay caused by the dilution mechanism.

Discussion
Quantitative models for de novo evolution of drug resistance are an important tool in tackling bacterial antimicrobial resistance, as well as viral infections and cancer. However, our quantitative understanding of how resistance emerges is still limited. The possibility of a phenotypic delay between the occurrence of a genetic mutation and its phenotypic expression has long been discussed [29,30,31,32,33], but its relevance for bacteria evolution has been questioned until recently [35]. Here, we have used computer simulations and theory to study the effects of phenotypic delay on the emergence of bacterial resistance to antibiotics. We investigated three different mechanisms that could lead to phenotypic delay: (i) dilution of antibiotic-sensitive molecules, (ii) effective polyploidy, and (iii) accumulation of resistance-enhancing molecules. We observe that the third mechanism only leads to phenotypic delay under a limited range of parameters, which makes it unlikely to be biologically relevant. The other two mechanisms have different "control parameters" (the degree of ploidy c versus the number of target molecules n) and different effects on the population dynamics. In particular, we show that protein dilution, but not effective polyploidy, can affect the probability that a growing population survives antibiotic treatment. The same mechanism is responsible for biasing the estimated mutation rate in a Luria-Delbrück fluctuation test. Effective polyploidy does not play a role here because of two cancelling effects: increased ploidy increases the number of mutations per cell in the growing population, but also increases the length of the phenotypic delay. These effects counterbalance such that the Luria-Delbrück distribution remains unaffected [35].
Effect of the dilution mechanism on fluctuation test data. Luria-Delbrück fluctuation tests are the standard microbiological method still being used for estimating mutation rates, yet it has often been noted that the measured distributions of mutant numbers do not quite fit the theoretical distribution [33,30,34,52]. A comparison with a more direct approach (DNA sequencing) suggests that fluctuation tests can significantly underestimate mutation rates [51]. Although phenotypic delay has been suggested as a possible explanation for these effects [33,51], our study is the first to investigate in detail how specific mechanisms of phenotypic delay alter the shape of the Luria-Delbrück distribution, and to demonstrate that it can indeed produce a biased mutation rate estimate of the same order of magnitude as that observed experimentally [51]. We also show that the simulated distribution from the dilution model fits the experimental fluctuation test data of Boe et al. [52] better than the standard model without phenotypic delay. We note that this result should however be taken cautiously. Boe et al.'s experimental protocol is not ideal for detecting phenotypic delay: for example, their bacterial cultures were allowed to reach stationary phase before plating. Moreover, our work shows that while phenotypic delay due to dilution affects the mutant number distribution, the change is subtle, requiring many replicate experiments to produce statistically significant results. While the usual number of replicates in a fluctuation test is less than 100, recent developments in automated culture methods should make it possible to run fluctuation tests with many more replicates, which may provide a way to probe the effects of phenotypic delay on the Luria-Delbrück distribution in more detail.
From molecular detail to evolutionary population dynamics. Our work presents an example of how molecular details at the intracellular level (here, protein dilution and the details of DNA replication) can have a direct effect on evolution at the population level [82,83,84]. This observation complements other works showing that, for example, molecular processes such as transcription and translation affect population-level distributions of protein numbers [85,86] and that noise in gene expression can directly affect the survival of populations in a fluctuating environment [87].
Importantly, both the effective polyploidy mechanism and the dilution mechanism cause a phenotypic delay only if the resistance mutation is recessive. For effective polyploidy this implies that a cell must contain only resistant alleles in order to be phenotypically resistant, while for the dilution mechanism we have assumed that sensitive target molecules need to be diluted out (or otherwise removed). This implies that we would expect to see phenotypic delay in the evolution of resistance to some antibiotics, but not to others. In particular, we would expect phenotypic delay due to dilution if the antibiotic acts by binding to its molecular target to make a toxic adduct, and resistance involves production of a resistant target. This is the case for fluoroquinolone antibiotics, which bind to DNA gyrase, causing DNA double-strand breaks; resistance is caused by production of mutant gyrase that no longer binds the antibiotic [88]. The fact that both Boe et al. [52] and Lee et al. [51] observed discrepancies in fluctuation test data for resistance to the fluoroquinolone nalidixic acid, is consistent with this expectation.
Assumptions of the model. Our simulations and theoretical calculations have involved a number of simplifying assumptions. Firstly, we ignore any possible fitness costs of mutations, assuming equal growth rates for wild-type and mutant cells in the absence of the antibiotic. While many resistance mutations incur a fitness costs [89,90], many clinically-relevant mutations have either no cost or even a small growth advantage [90,91].
For the molecular dilution mechanism, we have assumed that the degradation rate of target molecules is negligible, so that sensitive molecules can only be removed through cell division and dilution. While this seems to be (mostly) the case for bacterial enzymes targeted by antibiotics [92,93], it does not have to be the case for mammalian cells in which degradation plays a bigger role than dilution [94].
We have also assumed that in the dilution mechanism, all sensitive molecules need to be removed for the cell to become phenotypically resistant. In reality, resistance is likely to gradually increase as the number of sensitive molecules decreases. Our general conclusions remain valid in this case, however the mutant distribution may change. To construct more accurate models, we need measurements of the degree of antibiotic sensitivity as a function of the intracellular numbers of resistant and sensitive antibiotic targets. While technically challenging, such measurements could be carried out e.g. by fluorescently labelling target molecules [62].
Experimental tests for phenotypic delay Our work suggests that, at least in principle, the mutant number distribution could be used to detect the existence of a phenotypic delay caused by molecular dilution, although this would require many replicate experiments. Sun et al. have demonstrated phenotypic delay due to effective polyploidy more directly by tracking expression of a genetically engineered fluorescent marker in bacterial lineages [35]. Another approach would be an experiment similar to that from figure 1, in which a mutagen such as UV irradiation creates a burst of mutants. Other signatures of phenotypic delay may be detected in experiments where the timing of antibiotic exposure, and of resistance evolution, can be precisely controlled in turbidostat-like continuous culture devices [95].
Broader significance of phenotypic delay We have shown here that phenotypic delay (caused by molecular dilution) can affect mutation rate estimates from fluctuation tests, as well as the probability that a bacterial infection survives antibiotic treatment. Phenotypic delay may also affect other processes. For example, it was recently shown that a delay in evolutionary adaptation can lead to coexistence of spatial populations, in cases where immediate adaptation would eradicate coexistence [96,97]. In another example, it was suggested that in order to explain the effect of antibiotic pulses of different lengths on the probability of emergence of antibiotic resistance, a delay in evolutionary adaptation (called physiological memory by the authors) must be taken into account [98].

Methods and models
In all our simulations we use an agent-based model to simulate how mutated cells gain phenotypic resistance. Each cell has a number of attributes depending on the studied mechanism, such as the numbers of sensitive and mutated DNA copies, and the numbers of sensitive and resistant proteins, as specified below. Cells divide after time t d since last division.
For the population-level simulations (section 2.1), we simulate 100 cells which have just become genotypically resistant. Population-level simulations are repeated 1,000 times and single-cell simulations are repeated 10,000 times.

Modelling effective polyploidy
To describe how the copy number (ploidy) c changes during cell growth and division we use the Cooper-Helmstetter model [37]. We assume that it takes t 1 = 40 min for a DNA replication fork to travel from the origin of replication to the replication terminus, and that the cell divides t 2 = 20 min after DNA replication terminated (t 1 = C and t 2 = D in the original nomenclature of Ref. [37]; values representative for Escherichia coli strain B/r). During balanced ("steady state") growth assumed in this work, the number of chromosomes must double during the time t d between cell divisions (population doubling time). This means that for any t d < t 1 + t 2 = 60min, the cell must have multiple replication forks and more than one copy of the chromosome. The number of chromosomes will fluctuate during cell growth: it will double some time before the division, and decrease by a half just after the division. If t ini is the time since last division at which new replication forks are initiated, we must have (t ini + t 1 ) mod t d = t d − t 2 . This equation states that the time when a replication round initiated in the parent cell finishes in the offspring cell ((t ini + t 1 ) mod t d ) must be the same as the time t d − t 2 when the cell division process (lasting t 2 min) is initiated. It can be shown that this gives t ini = t d − (t 1 + t 2 ) mod t d . We proceed in a similar way to determine the time t rep at which a gene which confers resistance is replicated. If the gene is located in the middle of the genome, as is the case for the gyrA gene relevant for fluoroquinolone resistance, it will be copied t 1 /2 minutes after chromosome replication initiation. This implies that At this time point during the cell cycle the copy number of the gene of interest will double. The effective polyploidy immediately after this event is maximal and equal to where . . . denotes the ceiling function. We use c from Eq. (3) as the control variable in simulations of the polyploidy model. To simulate a cell or a population of cells with effective polyploidy we use the following algorithm. We initialize the simulation with all cells having c/2 sensitive alleles. Cells replicate in discrete generations every t d minutes. The number of allele copies double at t rep (Eq. (2)) since the last division in such a way that a sensitive/resistant allele gives rise to a sensitive/resistant copy, respectively. Sensitive alleles have a probability µ of mutating to a resistant allele. When a cell divides, the copies are split between the two daughter cells, with those linked by the most recent replication fork ending up in the same cell. We assume that the resistant mutation is recessive, which implies that a cell becomes resistant when all of its gene copies are resistant.

Modelling the dilution of sensitive molecules
For the dilution mechanism, we track the number of sensitive target molecules in each cell. We assume that at time zero, all cells have n sensitive target molecules and no resistant ones. When a mutation happens, the mutated cell begins to produce resistant target molecules and ceases to produce new sensitive molecules. At cell division, the sensitive molecules are partitioned between the two daughter cells following a binomial distribution with probability 0.5. We consider that a cell becomes phenotypically resistant when it contains no sensitive molecules. In the supplementary information we relax this assumption and study the case where a cell is considered resistant when the number of sensitive molecules falls below a (non-zero) threshold value (SI figure S5).

Modelling the accumulation of resistance-enhancing molecules
To model the accumulation of resistance-enhancing molecules, we explicitly simulate the production of M p resistance-enhancing molecules per cell cycle and their stochastic division between daughter cells, via a binomial distribution, at cell division. A cell is considered resistant when it contains more than M r resistance-enhancing molecules. In all simulations we fix M r = 1000 and vary M p to explore a range of m =

Combining effective polyploidy and molecular dilution
To include both effective polyploidy and molecular dilution, we track explicitly the total gene copies, the resistant gene copies and the number of sensitive proteins, as explained in Secs. 4.1 and 4.2. We assume that the number of resistant proteins produced in one cell cycle is proportional to the ratio of resistant to total gene copies. Both types of proteins (sensitive and resistant) are partitioned at cell division following a binomial distribution with probability 0.5. We consider that a cell becomes phenotypically resistant when it contains no sensitive molecules.

Simulating a growing infection
We start our simulations with 100 sensitive bacteria. Bacteria reproduce in discrete generations with doubling time t d . Upon reproduction, each bacterium can mutate with probability µ = 10 −7 . When the population reaches 10 7 cells, all phenotypically sensitive cells are removed (killed); this represent antimicrobial therapy. We repeat the simulation 1000 times to obtain the survival probability as a fraction of simulations in which phenotypically resistant cells emerge before the population dies out.

Simulating Luria-Delbrück fluctuation tests
To generate mutant size distributions for realistically large population sizes of sensitive cells required for comparing the model with experimental data, we use an algorithm based on the Çinlar's method [99,100]. The algorithm does not simulate the sensitive population explicitly, but it generates a set of times {t i } at which mutants emerge from the exponentially growing sensitive population (formally, t i is the time generated from a Poisson process with probability µλ s e λs(1−µ)ti ): is the final time, N f is the final population size, N i is the initial population size, λ s = ln(2) t d (1 − µ) is the growth rate of the sensitive bacteria, t d is the doubling time, µ is the mutation probability and U (0, 1) is a random variable uniformly distributed between 0 and 1.
For each t i , we then calculate the number of generations until the final time t f as For all of the simulations in sections 2.4 and 2.5, we assume t d = 60 min. We then simulate each clone for g i generations, including dilution of sensitive target molecules, and measure the number of resistant cells for each clone. Finally, we measure the number of resistant cells for each replicate by summing up over all clones.

Approximate Bayesian computation
We use an approximate Bayesian computation method to determine posterior probabilities of the non-delay and the dilution model. Briefly, the method relies on generating many (here: 10 4 ) independent samples of the simulated experiment mimicking Boe et al. [52] for both models. Model parameters are sampled from suitable prior distributions, we then select samples that approximate well the real data, and calculate the fraction of best-fit samples corresponding to each model. A single sample corresponds to 1104 simulated replicates of the fluctuation experiment at fixed parameters, for a given model. For each sample, parameters are randomly chosen from the . CC-BY-NC-ND 4.0 International license author/funder. It is made available under a The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/2019.12.19.883132 doi: bioRxiv preprint following prior distributions: log 10 (µ) uniform on [−10, −8], and log 2 (n) uniform on [0,8] (for the delay model). The tail cumulative mutation function 1 F (k) = Number of experiments yielding ≥ k mutants, 0 ≤ k ≤ 513, is calculated for each sample i (F i ), and also for the experimental data from Boe et al. [52] (F obs ). We then select 100 out of the 2 × 10 4 (10 4 from each model) generated samples with the smallest Euclidean distance ||F i −F obs || 2 (simulated distributions closest to the experimental data). The proportion of these which come from the phenotypic delay model is an approximation of the posterior probability that the experimental data was generated by the delay model (under the assumption that the experimental data was generated by one of the models). In reality the data generation process is likely to be far more complex than our idealised models, but the posterior probability of 0.97 implies the delay model provides a superior explanation to the model with no delay.
To examine the validity of our approach, we performed cross validation. For each model we randomly chose one sample corresponding to that model. We then computed the probability the simulated data was generated by the model with phenotypic delay, via the approximate method detailed above. This was carried out 500 times for each model. The proportion of simulations that were misclassified (as being with delay when they were not, or vice versa) was low (0.006, figure  5c), showing that our model selection framework is able to discriminate between the two models. We provide a further sensitivity analysis of this inference method in Section 7 of the SI.