Phenotypic Heterogeneity Facilitates Survival While Hindering the Evolution of Drug Resistance

Rising rates of resistance to antimicrobial drugs threatens the effective treatment of infections across the globe. Recently, it has been shown that drug resistance can emerge from non-genetic mechanisms, such as fluctuations in gene expression, as well as from genetic mutations. However, it is unclear how non-genetic drug resistance affects the evolution of genetic drug resistance. We develop deterministic and stochastic population models to quantitatively investigate the transition from non-genetic to genetic resistance during the exposure to static and cidal drugs. We find that non-genetic resistance facilitates the survival of cell populations during drug treatment, but that it hinders the development of genetic resistance due to competition between the non-genetically and genetically resistant subpopulations. The presence of non-genetic drug resistance is found to increase the first-appearance, establishment, and fixation times of drug resistance mutations, while increasing the probability of mutation before population extinction during cidal drug treatment. These findings advance our fundamental understanding of the evolution of drug resistance and may guide novel treatment strategies for patients with drug-resistant infections. SIGNIFICANCE Antimicrobial (drug) resistance is predicted to kill as many as 10 million people per year and cost over 100 trillion USD in cumulative lost production globally by 2050. To mitigate these socio-economic costs, we need to fundamentally understand the drug resistance process. We investigate the effect that different forms of drug resistance have on the evolution of drug resistance using mathematical modeling and computer simulations. We find that the presence of non-genetically drug-resistant cells (whose resistance is temporary and not encoded in a genetic mutation) allows the population to survive drug treatment, while at the same time slowing down the evolution of permanent genetic drug resistance. These findings have important implications for advancing evolutionary theory and for developing effective “resistance-proof” treatments.


INTRODUCTION
Antimicrobial (drug) resistance occurs when bacteria, viruses, fungi, and parasites no longer respond to drug therapy, making infections difficult or impossible to treat, which increases the risk of disease transmission, severe illness, and death (1). The evolution of drug resistance is well known to arise from the natural selection of genetic (DNA or RNA) mutations that provide microbes with the ability to survive and proliferate during treatment (2). More recently, it has been shown that non-genetic mechanisms promote microbial phenotypic diversification and survival strategies in selective drug environments (3).
Phenotypic heterogeneity has important implications for drug resistance (4), with heritable resistance potentially arising independently of genetic mechanisms (5). The stochastic or "noisy" expression of genes (6,7) introduces phenotypic variation among genetically identical cells in the same drug environment, which can result in the fractional killing of clonal microbial populations (3,8). This stochasticity is due in part to the inherently random nature of the biochemical reactions involved in the transcription and translation of genetic material, and can lead to the emergence of phenotypically distinct subpopulations within a clonal cell population (7,9). Non-genetic drug resistance (resistance that is not encoded in a genetic mutation) has been proposed to promote the development of genetic drug resistance (resistance that is acquired through genetic mutation) (4,10). This process may be enhanced by the interaction between non-genetic and genetic mechanisms inside the cell (4). For instance, promoter mutations can alter the expression noise levels of drug resistance genes (8), genetic networks can modulate gene expression noise to enhance drug resistance (11)(12)(13), and stress response genes can evolve elevated transcriptional variability through natural selection (14,15). Non-genetic mechanisms can facilitate the generation of genetic diversity by increasing the expression of key regulators involved in DNA replication, recombination, or repair (16,17), as well as by enhancing the adaptive value of beneficial mutations during drug treatment (18) and promoting the fixation of favorable expression level altering mutations (19). However, there are conflicting views on how phenotypic heterogeneity may facilitate adaptive evolution (20). Overall, the interactions between non-genetic and genetic mechanisms in drug resistance are largely unknown and remain to be quantified (4).
Non-genetic phenotypic variability can impact cellular population dynamics by providing a link between micro-scale dynamics (such as stochasticity at the molecular level) and macro-scale biological phenomena (including the fate of interacting cell populations) (21). Noise in biological systems may facilitate the adaption to environmental stress by allowing distinct, co-existing cellular states in a population to find the best adaptive solution from multiple starting points (22). Phenotypic heterogeneity can also promote interactions among subpopulations as well as the division of labor between individual cells, providing clonal microbial populations with new functionalities (23).
In this study, we investigate the transition from non-genetic to genetic resistance during static drug (drugs that slow or stop cell growth) and cidal drug (drugs that kill cells) treatment using deterministic and stochastic population models (24). Overall, we find that non-genetic resistance facilitates the initial survival of cell populations undergoing drug treatment, while hindering the establishment and fixation of genetic mutations due to competition effects between non-genetically and genetically resistant subpopulations. Quantitatively understanding the interaction between non-genetic and genetic mechanisms will improve our fundamental understanding of drug resistance evolution and may lead to more effective treatments for patients with drug-resistant infections.

Deterministic Population Model
The deterministic population model describes changes in cellular subpopulation concentrations over time under the influence of a drug. Three different subpopulations comprising the total population are described in this model: a susceptible subpopulation , a non-genetically resistant subpopulation , and a genetically resistant subpopulation . Cells may switch between the and subpopulations, and cells in the subpopulation can mutate into the subpopulation (Fig. 1). The mathematical model is described by a set of coupled ordinary differential equations (ODEs): where , is the switching rate from to , , is the switching rate from to , , is the mutation rate from to , and , , and are the death rates of , , and , respectively. There is no mutational pathway from to , as we are considering drugs that completely arrest growth and division of cells ( = 0) and therefore genetic mutation due to DNA replication errors do not occur. We also assume that genetic mutation provides complete and permanent resistance to the drug ( = 0).
To model cidal drug treatments of varying strength we scaled the deaths rates in Eqs. 1-2 by a factor . Summing Eqs. 1-3 yields the following equations: Resource competition between the subpopulations was modelled by scaling and by a Baranyi-Hill type function, which depends on and results in logistic growth (24). For subpopulation (where, ∈ { , }), this is given by: where is the growth rate for subpopulation (which leads to exponential growth in the absence of competition for limited resources), is the Hill coefficient, and ℎ is the point at which the competition function ( ) is half of its maximum value.
The growth dynamics of , , and were obtained numerically by solving the deterministic model, starting from initial population sizes , , and and integrating Eqs. 1-3 over a total time using time steps Δ . This numerical integration was performed using the ode45 ODE solver, which is based on an explicit Runge-Kutta method, in MATLAB (25).

Stochastic Population Model
Next, we developed a stochastic population model corresponding to the deterministic population model to study the effect of non-genetic resistance on the evolution of genetic resistance in low cell number regimes. Low numbers of infectious cells can occur during drug treatment and is the regime where stochastic fluctuations are expected to have a significant effect on population dynamics. Accordingly, Eqs. 1-3 were translated into the following set of reactions: Eqs. (7)-(13) were simulated using the Gillespie stochastic simulation algorithm (26, 27).
While the deterministic population model (Eqs. 1-3) was suitable for investigating large population dynamics under drug treatment, the corresponding stochastic population model (Eqs. 7-13) was necessary to accurately quantify establishment time, fixation time, and mutation first-appearance time distributions and extinction events in small populations undergoing cidal drug treatment.

Quantitative Measures
We used several quantitative measures to determine how non-genetic and genetic drug resistance subpopulations interact to produce a drug-resistant population.
The establishment time provides a quantitative measure of how long it takes for to make up a sufficiently large fraction of the total population, such that is unlikely to be affected due to stochastic fluctuations (i.e., the point at which is sufficiently large that genetic drift is weak and these mutants increases or decreases in frequency deterministically according to their relative fitness) (28).
was defined as the time it takes for to comprise 5% of the total population: The fixation time provides a quantitative measure of how long it takes for to become dominant in the population (28). was defined as the time it takes for to comprise 95% of the total population: To quantify the effect of non-genetic drug resistance on the evolution of genetic drug resistance in the stochastic population model, we obtained the first-appearance time ( ), establishment time ( ), and fixation time ( ) distributions of . For parameter regimes where population extinction could occur during the cidal drug treatment simulations (i.e., when and go extinct before appears), we determined the effect that the fitness of had on the probability of emerging ( ) before the extinction of the population. This was calculated from the number of population extinction events that occurred over a large number of simulations for different values of .

Deterministic Population and Evolutionary Dynamics Under Static Drug Exposure
To investigate the relative fitness effects of the non-genetically resistant and genetically resistant subpopulations on the evolution of drug resistance, we numerically solve the deterministic population model and simulate the stochastic population model. Specifically, we obtain the time series subpopulations concentrations/numbers, as well as the establishment, fixation, and first-appearance time distributions for drug resistance mutations that emerge during drug treatment. The parameters for the deterministic and stochastic population models are provided in the Supporting Material (Table S1).
The concentration of initially decreases after the application of the static drug as a result of cells switching from to , and then increases logistically due to switching from to ( Fig. 2A). The growth of follows a logistic-type curve (Fig. 2B). Overall, the growth of and increase as the fitness of increases ( Fig. 2A,B).
The concentration of (Fig. 2C) and the fraction of in the total population ( Fig. 2D) reveal that an increase in the fitness of decreases the expansion of . This can be attributed to resource competition between the subpopulations, as a higher total population size reduces the growth rate of (Eq. (6)).
The growth rate of the total population over time increases before sharply decreasing after it reaches a maximum (Fig. 2E). This is a result of the growth of the population beginning to slow down as it increases in size ( / → 0 as → ∞), which is expected for logistic-type growth (Eq. (6)) (24). Despite the decrease in the expansion of , the total population growth rate increases as the fitness of increases (Fig. 2E). Thus, increasing the fitness of in the static drug environment enhances the growth of the population, while at the same time hindering the expansion of .
We quantified how the fitness of the and subpopulations affects the establishment and fixation of the mutated subpopulation. As expected, an increase in relative to shortens the establishment time ( Fig. 3A) and the fixation time ( Fig. 3B) of in the population. Importantly, increasing relative to lengthens the establishment time ( Fig. 3A) and the fixation time ( Fig. 3B) of , due to competition decreasing the growth of ( Fig. 2C,D).

Deterministic Population and Evolutionary Dynamics Under Cidal Drug Exposure
Next, we investigated how the relative fitness of the non-genetically resistant and genetically resistant subpopulations affected the evolutionary dynamics of the population under cidal drug treatment.
The concentration of quickly declines after exposure to the cidal drug, with switching from providing temporary survival before dying off (Fig. 4A). The subpopulation shows slight temporary growth for higher values of before dying off, while lower values of produce flat growth curves before going extinct (Fig. 4B). Higher values prolong the temporary survival of and compared to lower values (Figs. 4A,B). The logistic growth of was unaffected by the fitness of ( Fig. 4C). Increasing the fitness of decreases the fraction of in the total population, as higher values result in more cells and consequently lower / values ( Figure 4D).
The number of cells in the total population initially declines (negative growth rate) before stabilizing at zero population growth (Fig. 2E). Then the population growth rate increases as the genetically drug resistant subpopulation expands to take over the population. This is followed by a decrease in the population growth rate, as the total population size moves towards saturation after the maximum population growth rate is reached.
Similar to static drug treatment, increasing the fitness of hinders the establishment time (

Stochastic Evolutionary Dynamics Under Cidal Drug Treatment
We simulated the stochastic population model corresponding to the deterministic population model to investigate the transition from non-genetic to genetic drug resistance in cell populations moving towards extinction during cidal drug treatment. This is important as fluctuations in small subpopulation sizes may impact the evolutionary dynamics of the population. Given that and are killed by differing degrees by the cidal drug, and that the evolution of depends directly on , we anticipated that stochastic fluctuations in the size of would in some cases lead to the extinction of the total population. To quantify the population and evolutionary drug resistance dynamics in this regime, we determined how the fitness of affects the probability of population extinction (which occurs when and go extinct before emerges), along with the first-appearance time distribution of a genetic mutation, and the establishment and fixation times which are bound to occur in our model once is present in the population.
Increasing the fitness of increased the likelihood of appearing and rescuing the cell population during intermediate strength cidal drug treatment ( = 2; Fig. 6A). This is in agreement with previous work that found that increasing the fluctuation relaxation time of a drug resistance gene increased the probability of acquiring a drug resistance mutation (5). When had low fitness in the cidal drug ( = 0.1733/hr) the population went extinct 83% of the time. The extinction probability decreased exponentially as the fitness of increased, with the population going extinct 69% of the time for moderate fitness ( = 0.2600/hr) and 0.00163% of the time for high fitness ( = 0.3466/hr). These results show that the presence of non-genetic resistance enhances the survival of the population when there are no pre-existing drug resistance mutations prior to drug exposure, and that non-genetic resistance increases the chance that a mutation will occur by providing the population with more time before extinction due to drug treatment. As expected, always appeared in the population when the strength of the cidal drug was low ( = 1), and conversely, the population almost always went extinct before could emerge when the strength of the cidal drug was high ( = 4). Decreasing the fitness of resulted in an earlier but more gradual transition from population survival to population extinct as the strength of the cidal drug increased (Fig. 6B).
When the the fitness of increased in the cidal drug environment, so did the means of the mutation first-appearance time, establishment time, and fixation time distributions (Fig. 7). This indicates that the presence of non-genetic drug resistance slows the evolution of genetic drug resistance, in agreement with the results obtained from the deterministic population model ( Fig. 3 and Fig. 5). While the coefficient of variation (CV; defined as the standard deviation divided by the mean) of the first-appearance time distributions (Fig. 7A,D,G) was only marginally dependent on the fitness of , the CV of the establishment time distributions (Fig. 7B,E,H) and fixation time distributions (Fig. 7C,F,I) increased approximately three fold as the fitness of increased from low to high. Therefore, the presence of increased non-genetic drug resistance is predicted to not only slow down genetic drug resistance, but also to increase the uncertainty in its evolution.

CONCLUSION
We found using deterministic and stochastic population models that non-genetic resistance enhances population survival during drug treatment, while hindering the evolution of genetic resistance. More specifically, increasing the fitness of the non-genetically resistant subpopulation (which allows the population to survive initial drug exposure) exponentially increased the chance of a genetically resistant subpopulation appearing and rescuing the population from extinction (29) during cidal drug treatment. However, increasing the fitness of the non-genetically resistant subpopulation slowed down the establishment and fixation of genetic drug resistance mutations due to competition effects, when no pre-existing mutations were present in the population. Incorporating pre-existing mutations into the model is expected to reduce the timescales of the establishment and fixation, as it would remove the growth delay of the genetically resistant subpopulation which results from the time it takes for non-genetically resistant cells to mutate into genetically resistant cells.
It will be important to study the effects of non-genetic resistance on the development of genetic resistance in the context of fluctuating environments, which may be governed by environment-sensing genetic networks (30), along with cellular trade-offs that may occur in drug environments (31). Fluctuating environmental stressors have been shown to facilitate "bet-hedging" in cell populations (32, 33), whereby some cells adopt a non-growing, stress-resistant phenotype to increase the long-term fitness of the population. This could be modeled using stochastic hybrid processes (34), for instance by using an stochastic ON-OFF switch coupled to a system of ordinary differential equations describing subpopulation dynamics in the presence of a drug. Furthermore, the first-appearance time, establishment time, fixation time, and extinction events could be described analytically in future studies using a first-passage time framework (5,35).
A complete understanding of the drug resistance process, including the interplay between non-genetic and genetic forms of drug resistance, will be important for mitigating the socio-economic costs of antimicrobial resistance (4). The drug resistance mutation appearance probabilities and first-appearance time distributions determined using stochastic population models may prove useful for guiding patient treatment. During the course of treatment, a drug could be substituted or combined with another drug (36, 37) before permanent, genetic drug resistance is predicted to emerge. Fluctuations in mutation first-appearance times are also important, as they can be the difference between the eradication or the establishment of drug-resistant infection during drug therapy. Finally, as drug exposure is more generally a form of selective pressure, the results of this study are anticipated to be applicable to evolutionary theory.

SUPPLEMENTARY MATERIAL
An online supplement to this article can be found by visiting BJ Online at http://www.biophysj.org.     (B) in hours. As is less than or equal to , numerical simulation data does not appear in the upper diagonal of the heat maps. The cidal death rate scaling factor was set to = 2 for these simulations.

Parameters
The parameter values used in our study (Table 1) were based on the literature or were scanned over a range of values. Drug susceptible cells S cells formed the majority of the total initial population T i , and non-genetically drug resistant cells N formed 1-10% of T i [2]. We assumed there were no pre-existing genetic drug resistance mutations (G i = 0) for all numerical and stochastic simulations. Numerical simulation of the deterministic population started with S i = 5.5x10 5 cells/mL and N i = 5.5x10 4 cells/mL, which are common initial concentrations for 'log-phase' laboratory experiments [1].
Depending on the concentration of the drug being considered, constant growth rates k N and k G (which here determine the fitness of the given subpopulation in the presence of a drug) were assigned values between 0.1733/hr and 0.3466/hr, based on growth rates measured in standard yeast cell culture experiments [4]. These parameter ranges were the basis of parameter scans that were used to predict how the relative fitness between N and G will affect population dynamics and the evolution of genetic drug resistance in our study.
To model partial drug resistance due to gene-expression noise, it was assumed that k N ≤ k G , with genetic drug resistance mutations providing the greatest level of resistance. Genetic mutations were also assumed to be permanent in our simulations.
The switching rates between S and N are based on experimental estimates [3] and ranged between r S,N = 0.0035/hr and r N,S = 0.0625/hr. The mutation rate from N to G was based on a previous modeling-experimental study [2], which ranged from 10 −6 to 10 −7 per cell division, and was assigned a value of r G,N = 0.3x10 −6 /hr in our simulations. The small magnitude of r G,N relative to the growth and switching rate parameters was found to cause negligible differences in our simulation results.
The death rates δ S and δ N in cidal drug environments were scanned over similar order of magnitudes as the growth rates. To satisfy the assumption that the N subpopulation is more fit than the S subpopulation, we considered constant death rates given by: where f δ is a scaling constant. To model low-strength, intermediate-strength, and high-strength cidal drug treatments we performed simulations for f δ = 1, f δ = 2, and f δ = 4, respectively. The deterministic and stochastic population models provide a means to simulate the effects of non-genetic resistance on the evolution of drug resistance. The effects of static drugs (which reduce cell growth but do no kill cells) were modelled by arresting the growth of S (k S = 0) and reducing the growth rates of N and G (k N and k G , respectively) and by setting the death rates to zero (δ S = δ N = δ G = 0). Cidal drugs (which eventually kill all non-genetically resistant cells) were modelled by non-zero death rates for S and N (δ S and δ N , respectively). To model the partial drug resistance of N , δ N < δ S .
The parameter values given in Table S1 were also used for the stochastic simulations, as all the corresponding reactions in the stochastic population model were of zeroth-order or first-order [5]. Note that in the stochastic population model the values of S i , N i , and G i are exact numbers of cells (S i = 5.5x10 5 cells and N i = 5.5x10 4 cells), which we used as the initial conditions to investigate the stochastic transition from non-genetic to genetic resistance as the number of cells in the population approached zero (extinction) due to cidal drug treatment. The S to N (and vice-versa) phenotype switching rates (r N,S and r S,N , respectively), the N to G mutation rate (r G,N ), and the growth and death rates (k i and δ i , respectively, where i ∈ {S, N, G}) are all given as probability per unit time in the stochastic population model.

Low-Strength Cidal Drug Condition
For the low-strength cidal drug environment, we considered a cidal death rate scaling factor of f δ = 1.
Subpopulation trajectories, the fraction of genetically drug resistant cells in the population G/T , and the population rate of change dT /dt for the lowstrength cidal drug scenario are shown in Fig. S1.
Establishment time and fixation time heat maps for this case are shown in Fig. S2.

High-Strength Cidal Drug Condition
For the high-strength cidal drug environment, we considered a cidal death rate scaling factor of f δ = 4.
Subpopulation trajectories, the fraction of genetically drug resistant cells in the population G/T , and the population rate of change dT /dt for the lowstrength cidal drug scenario are shown in Fig. S3.
Establishment time and fixation time heat maps for this case are shown in Fig. S4.

Survival and Extinction Scenarios
Representative survival and extinction trajectories from the stochastic simulations are shown in Fig. S5. Fig. S5A shows a case where the cell population survives cidal drug treatment due to G emerging before S and N go extinct. G subsequently takes over the population in this scenario. In contrast, Fig. S5B shows a case where the cell population goes extinct due to S and N reaching zero cells before G appears.