Theoretical considerations and empirical predictions of the pharmaco- and population dynamics of heteroresistance

We open this consideration of the pharmaco- and population dynamics of HR with a description of the two mathematical models employed. For both models of HR, we assume a Hill function for the relationship between the concentration of the antibiotic, the concentration of the limiting resource, and the rates of growth and death of the bacteria, known as the pharmacodynamics (13–15).

The two models of HR used here are depicted in Fig. 1. In the progressive model (Fig. 1A), the increasingly resistant states are generated by a transition from a less resistant state to a more resistant state, and the more resistant states generate the less resistant states sequentially. In the non-progressive model (Fig. 1B), the different resistant states are generated directly by a transition from the susceptible state, and the more resistant states transition directly back to the most sensitive state.

In this model (Fig. 1A), the bacteria transition between four different states: sensitive, S, and increasingly resistant, R1, R2, and R3, which are the designations and densities in cells per mL of bacteria of these different states. The total number of cells in a state is given by the product of the density of the state and the total volume, Vol. Cells of the S state transition to R1, R1 transitions to R2, and R2 transitions to R3 at rates µ_S1, µ_R12, and µ_R23 per cell per hour, respectively. Cells of resistant states progressively transition to the less resistant states, R3 to R2, R2 to R1, and R1 to S, with rates µ_R32, µ_R21, and µ_R1S per cell per hour. We simulate these transitions with a Monte Carlo process to account for the stochasticity that occurs during cell division related to either gene amplification or a point mutation at a specific locus due to the polymerase error rate (16). In this Monte Carlo process, a pseudorandom number x (0 ≤ x ≤1) is generated from a rectangular distribution (17). If x is less than the product of the total number of cells in the generating state, the transition rate (µ), the Monod function ($\psi \left(r\right)$), and the step size (dt) of the Euler method (18) employed for solving differential equation, then 1/(dt*Vol) cells are added to the recipient population and removed from the generating population. In the equations below, the transition of cells due to this process is expressed as MIJ, e.g., MR1S for transitions from R1 to S and MR21 for transitions from R2 to R1. For example, if S is 10⁵, µ_SR1 is 10⁻⁷, dt is 10⁻⁴, Vol is 1, and x < S*µ_SR1*dt*Vol* $\psi \left(r\right)$, which in this case is 10⁻⁶* $\psi \left(r\right)$, then 1/(dt*Vol) cells are removed from the S population and added to the R1 population, which in this case would mean that MSR1 is 10⁴ cells. Practically, when MSR1 is plugged into the differential equations below, it means that for any given step of the Euler method from t to t + dt we add (1/(dt*Vol))*dt cells, so here, there would be one cell added. In this model the assumption that there are 1 sensitive and 3 resistant states is arbitrary, but in reality, there could be more or fewer states depending on the mechanism of HR. In the simulations presented in the main text, we assume the transition rates between states are equal in both directions and between states. We also assume that the transitions between states slow down with the decline in the physiological state of the cells in direct proportion to $\psi \left(r\right)$. With these definitions, assumptions, and the parameters defined and presented in SI Appendix, Table S1, the rates of change in the densities of the different populations are given by,

In this model (Fig. 1B), all the resistant states, R1, R2, and R3 are derived from the susceptible state and, by transition, return directly to the susceptible state, S. The rates of transition from state S are, respectively, µ_SR1, µ_SR2, and µ_SR3 per cell per hour. The rates of return to the susceptible state are µ_R1S, µ_R2S, and µ_R3S per cell per hour. The transitions between states are via a Monte Carlo process (16), using a routine like that for the progressive model. When transients from S to the different R states are generated (MSR1, MSR2, and MSR3), 1/(dt*Vol) are added to the R1, R2, and R3 populations and are removed from the S population. When transients from the R1, R2, and R3 populations are generated (MR1S, MR2S, and MR3S), 1/(dt*Vol) are added to the S population and removed from the R1, R2, and R3 populations, respectively. Here, the primary assumption is that all resistant states are derived from and transition back to the sensitive state and we continue to assume there are four states with equal transition rates between them. With these definitions, assumptions, and the parameters defined and presented in SI Appendix, Table S1, the rates of change in the densities of the different populations are given by:

Here, we consider the population dynamics of HR with the distributions of the different resistant states generated from single cells grown up to full densities for the progressive and non-progressive models with four transition rates (Fig. 2). We further consider a greater range of transition rates for the non-progressive model to determine the minimum rate for which we generate sufficiently large minority populations in SI Appendix, Fig. S2 A–C.

For the progressive model, only in the runs with the highest transition rates, µ = 10⁻² and µ = 10⁻³ per cell per hour, is the subpopulation with the highest resistance level, R3, present. A very different situation obtains for the non-progressive model, as at every transition rate the R3 population is present. We also consider the effect that the relative fitness cost of each state has on these stationary phase densities (SI Appendix, Fig. S3) and find modest differences in these distributions.

Using these same parameters for both models, another difference can be seen between the progressive and non-progressive models in the PAP tests of each (Fig. 3). For these PAP tests, we calculate the ratio of the number of cells generated at a particular antibiotic concentration compared to the number of cells present when there is no antibiotic [N(A)/N(0)] for 0, 1, 2, 4, 8, and 16 times the MIC of the susceptible population.

The PAP test results anticipated from the progressive model are very different than those anticipated from the non-progressive model, two extreme HR cases. The presence of four subpopulations with different MICs is apparent from the PAP test of the progressive model with the parameters considered. For the non-progressive model, the differences in the relative densities of the subpopulations are too low to be detected by a PAP test performed in the lab. In general, the plateaus shown in Fig. 3 A and B are sharper and more dramatic than would be seen in the lab. This is a consequence of having only three resistant states. Moreover, using the standard HR criteria of having a subpopulation with an MIC of >8 times at a frequency of at least 10⁻⁷, the progressive model only meets these criteria at high transition rates (exceeding 10⁻⁴). On the other hand, the non-progressive model meets these criteria at transition rates as low as 10⁻⁷ (SI Appendix, Fig. S2D).

To explore how these models differ in their response to antibiotic treatment, we follow the changes in the densities and MICs of heteroresistant populations exposed to two antibiotic concentrations (5 µg/mL and 10 µg/mL corresponding to 5× and 10× the MIC of the susceptible population) when µ = 10⁻² and µ = 10⁻⁵ per cell per hour (SI Appendix, Fig. S4). We initiate these simulations with 1/100 of the stationary phase densities of the different states anticipated from the heteroresistant populations depicted in Fig. 2 A,D, H, and E for the progressive and non-progressive models of heteroresistance, respectively. In SI Appendix, Fig. S5, we consider the effect that the fitness cost of resistance has on these dynamics and find the effect modest at best, just slowing the response time to the antibiotic.

There are apparent differences in the bacterial response to antibiotics between the progressive and non-progressive models of HR. For both models, when µ = 10⁻² per cell per hour, the R3 population comes to dominate and the MIC increases to the maximum (15 µg/mL), though with the higher concentration of the drug, it takes longer for the R3 population to become dominant. With the lower transition rate of µ = 10⁻⁵ per cell per hour, at 5 µg/mL of antibiotic, the R2 population comes to dominate in the progressive model and the R3 population remains minor. At this same transition rate and at 10 µg/mL, resistance does not evolve, and the bacterial populations are lost. In both cases, the MIC does increase but does not go to the maximum value. For these conditions, in the non-progressive model, subpopulations are always able to respond to the antibiotic and are never eliminated.

Upon removal of the antibiotic, the heteroresistant bacterial population reverts to the sensitive state. This reversion is the case for both the progressive and non-progressive models of HR considered here. To illustrate this and elucidate the relative contributions of the rates of transition between states and the fitness cost of resistance (as measured by the growth rates) to the dynamics and the time needed to restore susceptibility, we use serial transfer forms of the progressive and non-progressive versions of the HR models. In these simulations, the populations are grown for 24 h, diluted by a factor of 100, and fresh resources added. In Fig. 4, we present the results of simulations of the changes in the densities of the susceptible and resistant populations as well as the change in average MIC in serial transfer following the removal of the antibiotics. These serial transfer simulations were initiated with 10⁷ bacteria per mL of the highest resistance level, R3. We consider two major conditions: one where the fitness cost of resistance is high and another where the fitness cost of resistance is low. In SI Appendix, we consider the dynamics of reversion when a set of even higher fitness costs are used (SI Appendix, Fig. S6).

In the absence of antibiotics, the populations become increasingly dominated by more susceptible populations for both the progressive and non-progressive models of HR. This change in the composition of the populations is also reflected in a decline in the average MIC, approaching the level of the susceptible population. With the same fitness parameter and transition rates between states, µ, the rate of return to the susceptible state is greater for the non-progressive model than the progressive model. For both models, the rate of return to the sensitive state is proportional to the transition rate between states, µ, and the relative fitness cost of resistance. Notably, the new apparent equilibria obtained for both models differ substantially. In the progressive model, the most resistant populations are in continuous decline and will ultimately be lost or nearly so, while in the non-progressive model, all resistant populations are present at roughly equal frequency and appear to be in equilibrium. Of note is the vast difference in the time needed for the susceptible population to come to dominate; we list these times in SI Appendix, Table S2.

For our numerical analysis of the coupled, ordered differential equations presented (Eqs. 2–11), we used Berkeley Madonna with the parameters presented in SI Appendix, Table S1. Copies of the Berkeley Madonna programs used for these simulations are available at https://www.eclf.net. In the analysis of our simulations, to calculate average MIC, we take a weighted average of the MIC of each population.

Enterobacter cloacae Mu208 is a carbapenem-resistant isolate collected by the Georgia Emerging Infections Program Multi-site Gram-negative Surveillance Initiative and described previously (23). Burkholderia cepacia complex isolate JC8 is a cystic fibrosis patient isolate collected by the Georgia Emerging Infections Program Multi-site Gram-negative Surveillance Initiative. E. coli MG1655 was obtained from the Levin Lab’s bacterial collection.

Single colonies of E. coli MG1655 were inoculated into 10 mL lysogeny broth (BD, Product #244610) and grown overnight at 37 °C with shaking. Cultures were serially diluted in saline and all dilutions (10⁰ to 10⁻⁷) plated on LB agar plates (BD, Product #244510) containing 0, 1, 2, 4, 8, and 16 times the MIC of rifampin (Thermo Fisher, Product #J60836.03). Plates were grown at 37 °C for 48 h before the density of surviving colonies was estimated.

Single colonies of B. cepacia complex isolate JC8 and E. cloacae Mu208 were inoculated into 1.5 mL Mueller–Hinton broth (BD, Product #275730), and cultures were grown overnight at 37 °C with shaking. Cultures were serially diluted in phosphate-buffered saline, and 10 µL of each dilution was spotted on Mueller–Hinton agar (BD, Product #225250) plates containing 0, 0·125, 0·25, 0·5, 1, 2, and 4 times the breakpoint concentration of each antibiotic. Antibiotics used were ticarcillin disodium (BioVision, Product #B1536) with clavulanate potassium salt (Cayman Chemical Company, Procut #19456), amikacin sulfate (AstaTech, Product # 40003), colistin sulfate salt (Sigma-Aldrich, Product # C4461), and fosfomycin disodium salt (TCI America, Product # F0889). For Mu208 on fosfomycin, broth and agar included 25 µg/mL glucose-6-phosphate (Sigma-Aldrich, Product #10127647001). Plates were maintained at 37 °C overnight for Mu208 and for 36 to 60 h for JC8.

B.R.L., B.A.B., T.G.-G., and F.B. designed research; B.R.L., B.A.B., T.G.-G., J.A.M., J.E.C., and D.I.A. performed research; B.R.L., B.A.B., T.G.-G., and J.A.M. contributed new reagents/analytic tools; B.R.L., B.A.B., T.G.-G., J.A.M., D.I.A., D.S.W., and F.B. analyzed data; and B.R.L., B.A.B., T.G.-G., J.A.M., A.P.S., J.E.C., D.I.A., D.S.W., and F.B. wrote the paper.

The authors declare no competing interest.