Allele-specific collateral and fitness effects determine the dynamics of fluoroquinolone resistance evolution.

Collateral sensitivity (CS), which arises when resistance to one antibiotic increases sensitivity toward other antibiotics, offers treatment opportunities to constrain or reverse the evolution of antibiotic resistance. The applicability of CS-informed treatments remains uncertain, in part because we lack an understanding of the generality of CS effects for different resistance mutations, singly or in combination. Here, we address this issue in the gram-positive pathogen Streptococcus pneumoniae by measuring collateral and fitness effects of clinically relevant gyrA and parC alleles and their combinations that confer resistance to fluoroquinolones. We integrated these results in a mathematical model that allowed us to evaluate how different in silico combination treatments impact the dynamics of resistance evolution. We identified common and conserved CS effects of different gyrA and parC alleles; however, the spectrum of collateral effects was unique for each allele or allelic pair. This indicated that allelic identity can impact the evolutionary dynamics of resistance evolution during monotreatment and combination treatment. Our model simulations, which included the experimentally derived antibiotic susceptibilities and fitness effects, and antibiotic-specific pharmacodynamics revealed that both collateral and fitness effects impact the population dynamics of resistance evolution. Overall, we provide evidence that allelic identity and interactions can have a pronounced impact on collateral effects to different antibiotics and suggest that these need to be considered in models examining CS-based therapies.

Antibiotics are a cornerstone of the prevention and treatment of bacterial infections; however, their efficacy is rapidly declining due to the emergence and spread of antibiotic resistance (1). Advances in antibiotic discovery and design have not kept pace with resistance evolution (2), necessitating experimentally validated treatment strategies (3). Selection inversion, a strategy that aims to circumvent the emergence and dissemination of antibiotic resistance by combining existing antibiotics based on their physiological and/or evolutionary interactions, has recently gained prominence (4). Among these, collateral sensitivity (CS)-informed strategies are particularly promising (5).

CS occurs when mutations conferring resistance to one antibiotic increase sensitivity toward other antibiotics in the same or different functional class (6). Because of this effective trade-off, bacteria treated with a pair of drugs exhibiting reciprocal CS are unable to simultaneously evolve resistance to both agents (4). Recent studies have examined the frequency and mechanisms of CS, although these have been largely restricted to laboratory strains of a small number of bacterial species (5, 7–14). In addition, none of these studies have examined if collateral effects are conserved for different clinically circulating resistance alleles to the same antibiotics. This is especially relevant given that resistance to any given antibiotic can arise from mutations in different genes or at different sites within a gene (15). Because these resistance alleles cause distinct phenotypic effects, both with respect to changes in minimum inhibitory concentration (MIC) and bacterial fitness (16), they may also underlie a distinct spectrum of collateral effects, only some of which would be suitable for CS-informed therapies. To address this limitation, we examine collateral responses to distinct genes and alleles that confer resistance to fluoroquinolones (FQs) in the Gram-positive pathogen Streptococcus pneumoniae.

S. pneumoniae invasive infections are responsible for the most deaths among vaccine-preventable diseases globally (17). FQs are a mainstay of treatment against invasive pneumococcal disease (18), but successful FQ treatment is threatened by the emergence of FQ-resistant strains, which have been reported to be as high as 10.5% (19). De novo FQ resistance in S. pneumoniae arises predominantly via the stepwise accumulation of chromosomal mutations in the quinolone resistance-determining regions (QRDRs) of the DNA gyrase (gyrA) and/or topoisomerase IV (parC) genes (19). Mutations in either gyrA or parC result in low-level FQ resistance, whereas the combination of mutations in both genes results in high-level resistance, often affecting multiple agents within the class (18, 20). Alleles carrying different mutations in QRDR regions, in pneumococci and other species, have varying impacts on FQ resistance and on the fitness of strains carrying them (21, 22); these effects are believed to influence the population frequencies of different alleles during monotherapy (23). If similar phenotypic heterogeneity exists for collateral responses, this could impact the efficacy of CS-based combination therapies and the persistence and spread of FQ resistance.

Here, we studied collateral effects of FQ resistance by generating S. pneumoniae mutant strains via allelic replacement that encode clinically relevant FQ resistance mutations in gyrA and parC (19). These mutants were used to assess how different alleles conferring resistance to the same antibiotic influence the susceptibility to other antibiotics and whether interactions between these alleles modulate their collateral effects. We then developed and applied a mathematical model to quantitatively study the population dynamics of resistance evolution during different antibiotic combination treatment regimes.

Results

Extensive and Conserved Collateral Effects to FQ Resistance.

To investigate collateral effects to FQ resistance in S. pneumoniae, we generated a panel of 16 D39 strains harboring mutations in the QRDRs of gyrA and/or parC by natural transformation (Fig. 1). Specifically, we generated eight strains encoding alleles with FQ resistance-associated amino acid mutations in either gyrA (gx, with x = S81F, S81Y, E85G, or E85K) or parC (py, with y = S79F, S79Y, D83N, or D83Y) and eight strains with combinations of both alleles (gxpy). Each strain is hereafter denoted based on its FQ resistance-related amino acid mutations. All mutant strains had decreased susceptibility to ciprofloxacin, where the gxpy showed the most pronounced MIC increase (32 to 64 mg/L), followed by py (16 to 32 mg/L), and lastly gx (8 mg/L) (Fig. 1). In addition, we observed significant fitness reductions in transformed strains (Fig. 1), with gxpy strains exhibiting the lowest mean relative growth rate (0.56; with a range from 0.29 to 0.74), followed by gx (0.60; 0.41 to 0.77), and lastly py (0.85; 0.39 to 1.47) strains.

Fig. 1.

Experimental design and determined MIC and growth rates. (A) Overview of experimental workflow. Mutants of S. pneumoniae D39 strain encoding gyrA (gx, yellow), parC (py, light blue), or double gyrA::parC (gxpy, red) alleles that confer FQ resistance were generated via transformation followed by allelic replacement. Mutants were subjected to antimicrobial susceptibility testing for 12 clinically relevant antibiotics and ciprofloxacin (CIP). A decrease in MIC of each mutant relatively to that of the WT D39 strain (WT; black) was determined as CS and an increase as CR. The growth rate of the WT and each mutant was determined to estimate in vitro fitness. AB, antibiotic. (B) CIP MICs for the WT S. pneumoniae D39 (WT, black), gyrA mutants (gx, yellow), parC mutants (py, light blue), and their corresponding gyrA::parC double-allele mutant (gxpy, red). (C) Relative growth rate of gyrA (gx, yellow), parC (py, light blue), and gyrA::parC (gxpy, red) mutants compared to WT (black) using a two-sided t test. Significant differences are indicated by asterisks, as follows: *, P < 0.05; **, P < 0.01; ***, P < 0.001; and ****, P < 0.0001.

The collateral effects of these mutants were assessed against 12 antibiotics belonging to a broad range of classes, including commonly used antipneumococcal agents. Collateral effects were observed in ∼52% of the possible instances (100 out of 192) against all tested antibiotics (Fig. 2). Among these, 87% were CS and the remaining 13% were CR. For 7 of 12 antibiotics, collateral effects were in the same direction for all mutants (CS), whereas for the remaining antibiotics collateral effects were mixed between CS and collateral resistance (CR). All mutants (100%) exhibited CS to gentamicin, 81.3% to clindamycin, 75.0% to tetracycline, 56.3% to trimethoprim/sulfamethoxazole, 43.8% to penicillin and fusidic acid, and only 6.3% to vancomycin.

Fig. 2.

Collateral effect (CE) profiles for FQ-resistant S. pneumoniae, including gyrA (gx) and parC (py) mutants, and gyrA::parC double mutants (gxpy). Antibiotics are ranked based on the hierarchical clustering of their associated CEs. Color indicates the magnitude of the CE quantified as the mean log2 relative change of MIC compared to the WT. Red indicates CR and blue CS. CHL, chloramphenicol; CIP, ciprofloxacin; CLI, clindamycin; DAP, daptomycin; ERY, erythromycin; FUS, fusidic acid; GEN, gentamicin; LNZ, linezolid; PEN, penicillin; RIF, rifampicin; SXT, trimethoprim/sulfamethoxazole; TET, tetracycline; VAN, vancomycin.

Single-allele strains, namely, gx and py, exhibited CS to a median of four antibiotics with log-scale fold changes relative to the isogenic wild type (WT) varying from −0.42 to −1.58, and CR to a median of one antibiotic with a fold change of 0.42. Double gxpy strains exhibited CS to a median of seven antibiotics with log-scale fold changes relative to the isogenic WT varying from −0.42 to −2 and CR to a median of 0.5 antibiotic with fold changes of 0.42 to 1. Conserved CS effects to gentamicin and clindamycin were pronounced, with a median log-fold change in MICs of −1.2 and −1, respectively.

Hierarchical clustering revealed four antibiotic groups for which resistant mutants exhibited similar collateral effects (); only one of these groups consisted exclusively of functionally similar antibiotics, including the protein synthesis-inhibiting antibiotics gentamicin and clindamycin, confirming that collateral effects cannot be simply predicted based on the antibiotic target. Single- and double-allele mutants did not cluster distinctively based on their effects within the five strain clusters that we observed (); this suggests that many of the collateral effects in double mutants are potentially caused by the mutations of each single allele, although the magnitude of these effects is altered in the double mutants, as discussed below.

Heterogeneity in Collateral Effects to FQ Resistance.

Despite the extensive and conserved collateral effects that we observed, we found that the single-allele gyrA and parC strains exhibited heterogeneity in their collateral effects even between strains carrying different resistant alleles of the same gene (Fig. 2). For instance, strains carrying the gyrA alleles coding for S81F and S81Y exhibited CS to three and eight antibiotics, respectively, and a median fold decrease of −0.58 in MICs. Similarly, the strain carrying the parC allele coding for S79Y exhibited CS to three antibiotics with a median fold decrease of −1 in MICs, whereas the strain carrying the S79F allele showed CS to five antibiotics with a median fold decrease of −0.58. Interestingly, the heterogeneity was not limited to the number or intensity but also the direction of the effects. For example, strains carrying the gyrA alleles coding for S81F and S81Y exhibited CR and CS to linezolid, respectively; a similar difference was observed with daptomycin for strains carrying the gyrA alleles coding for E85G and E85K.

Collateral effects among double-allele gxpy strains (Fig. 2) varied in their number, intensity, and/or direction as well. For instance, the double-allele strain gS81FpS79F exhibited CS to seven antibiotics with a median fold decrease of −0.58 in MICs, whereas the strain gS81FpS79Y showed CS to three antibiotics with a median fold decrease of −0.57. Heterogeneity in the direction of the collateral effects was observed for erythromycin where several of the double-allele strains exhibited either CS (gS81YpS79Y, gS81FpS79F, and gS81FpD83N) or CR (gS81YpS79F, gS81FpD83Y, and gE85GpS79F). In addition, collateral effects in double-allele strains appeared to differ from those expected based on the strains carrying the respective single alleles, resulting in either lower (i.e., gS81YpS79Y for fusidic acid) or higher (i.e., gS81YpS79F for gentamycin) MIC for the collaterally affected antibiotic.

Antibiotic Resistance-Linked Fitness Does Not Predict Collateral Effects.

All gyrA and parC alleles, aside from the parC coding for S79F, caused significantly reduced growth rates compared to the WT strain (Fig. 1 and ). In addition, significant differences (ANOVA, P < 0.05) were observed between strains carrying different FQ-resistant alleles of either the gyrA (gS81F and gE85K absolute difference in growth rate (|ΔKG|) = 0.0134 h−1) or the parC (pS79F and pD83N |ΔKG| = 0.0663 h−1, pS79F and pD83Y |ΔKG| = 0.0431 h−1, pS79Y and pS79F |ΔKG| = 0.0461 h−1), but not between the strains carrying resistant alleles of both genes (). Our analysis indicated that there was no clear overall correlation between the fitness of the generated strains and their collateral effects (). This suggests that the fitness costs of resistance have only a limited impact on observed CS effects overall, with this impact being heterogeneous and antibiotic dependent.

Combination Treatment with Clinical Dosing Regimens Suppresses De Novo Antibiotic Resistance Development.

Given the distinct collateral and fitness effects of FQ resistance alleles, we used simulations to assess the emergence and fixation of resistant strains based on eight evolutionary trajectories (ETs) leading to high-level FQ-resistance (Fig. 3), where each ET included one double-allele mutant and its corresponding single-allele mutants. To this end, we developed a mathematical model and studied the impact of these differences on treatment outcomes (Fig. 3) under a clinical dosing schedule for ciprofloxacin monotreatment (500 mg twice daily [b.i.d.], average unbound steady state concentration [Css] 1.39 mg/L) and in combination with erythromycin (600 mg b.i.d., Css 0.48 mg/L), linezolid (600 mg b.i.d., Css 7.33 mg/L), and penicillin (3000 mg b.i.d., Css 6.95 mg/L). The model incorporated our experimentally derived antibiotic MICs, collateral and fitness effects, and antibiotic-specific pharmacodynamic relationships (details and estimated parameters can be found in ). For the different simulation scenarios, we then computed the probability of resistance evolution for each of the eight ET leading to high-level FQ resistance shown in Fig. 3. Resistance was defined as the end of treatment fixation of the gx, py, and/or gxpy subpopulation, which occurred when the subpopulation-specific bacterial density exceeded the resistance cutoff. The resistance cutoff was set to 104 colony-forming unit (CFU)/mL and represents an established infection (24).

Fig. 3.

Overview of treatment simulation workflow. (A) Trajectories leading to high-level FQ resistance (ETs) used for simulations, with each trajectory consisting of four subpopulations representing the WT S. pneumoniae D39 (WT, black), a gyrA mutant (gx, yellow), a parC mutant (py, light blue), and their corresponding gyrA::parC double-allele mutant (gxpy, red). The experimentally determined MICs and relative growth rate (KG) for each specific strain were used to inform the model. (B) Schematic of the stochastic hybrid ODE model use for simulations. The bacterial model comprised of ODEs describing the change of bacterial density over time for each subpopulation. The subpopulation-specific net bacterial growth (knet,z) was modeled considering capacity limitation (G), subpopulation-specific growth rate (k), and total antibiotic effect (E). The capacity limitation (G) was implemented using a maximal carrying capacity (B). Unidirectional mutations were implemented as a stochastic event using random sampling from a binomial distribution with a probability representing the mutation frequency (μ) and sample size equal to the density of each bacterial subpopulation. The bacterial model was linked to a pharmacokinetic model of CIP, i.e., AB1; and a pharmacokinetic model of erythromycin (ERY), linezolid (LNZ), or penicillin (PEN), i.e., AB2. k, k, and k are rate constants and V and V are the volumes of the central and peripheral compartment (only relevant for ERY), respectively. Antibiotic effect for CIP (E) and the second antibiotic (E) were included on each subpopulation and were driven by the simulated plasma concentrations of the respective drugs. Pharmacodynamic drug interactions (αAB1*AB2) were included for drug combinations for which significant nonadditivity was identified. (C) Treatments were simulated 500 times per ET. Pharmacokinetic profiles of an example combination treatment are depicted in the lower graph, while the nth realization of the bacterial dynamics corresponding to the example treatment are shown in the top graph. The probability of resistance (PR,Mz) (mutants) or failure of eradication (WT) was calculated for each trajectory and treatment by summarizing the number of simulations for which bacterial concentrations exceeded the resistance cutoff (dashed line) at end of the treatment for each individual subpopulation.

Our simulations of clinical ciprofloxacin dosing schedules, where Css equates to 0.9-fold of WT MIC (1.5 mg/L), revealed that monotreatment failed to eradicate the WT and promoted resistance evolution (Fig. 4 ) for all ETs due to the emergence of low-level FQ-resistant mutants (gx and/or py). However, adding a second drug effectively eradicated the WT (Fig. 4). Treatments using linezolid (1.5 mg/L, 4.9× WT MIC) or penicillin (0.012 mg/L, 594× WT MIC) resulted in a 0% probability of resistance, whereas erythromycin treatment (0.19 mg/L, 2.6× WT MIC) promoted the emergence of either exclusively py (75.6 to 80.4%) or both py (75.6%) and gx (77%) strains for all but one ET (ET 1). The emergence of gxpy double-allele strains was only observed for a single trajectory (ET 2), with a very low probability of resistance (0.2%). Overall, our findings support the idea that antibiotic combinations can minimize the risk of de novo antibiotic resistance evolution, but the benefit depends on the specific antibiotic combination.

Fig. 4.

Treatment outcomes of simulated clinical CIP mono- or combination therapies. Each simulation includes WT (black), a gyrA mutant gx (yellow), a parC mutant py (light blue), and the double-allele mutant gxpy (red). (A) Treatment outcomes of simulated clinical CIP monotreatment (CIP, 500 mg b.i.d., Css 1.39 mg/L) vary between different trajectories (ETs). Solid line represents the median bacterial growth dynamics of 500 simulations where bacterial subpopulations are indicated by the respective color. The shaded areas represent the 5th and 95th simulated percentiles. The black dashed line represents the resistance cutoff, which is equal to the WT bacterial density at the start of treatment. Although the bacterial dynamics and resistance development differed between the ETs, the simulated CIP treatment was ineffective for all scenarios. (B) End-of-treatment probability of resistance (mutants) and WT eradication failure during CIP monotreatment were calculated for different ETs (panels) and subpopulations (x-axis) separately. (C) Treatment outcome for clinical dosing regimens of CIP combination therapies. Simulation outcomes with clinically relevant dosing regimens of CIP in combination with ERY (600 mg b.i.d., Css 0.48 mg/L), LNZ (600 mg b.i.d., Css 7.33 mg/L), or PEN (3 g b.i.d., Css 6.95 mg/L). CIP as monotreatment had a high probability of resistance, foremost selecting for strains encoding FQ-resistant parC alleles (py, light blue) and secondly gyrA alleles (gx, yellow). The addition either LNZ or PEN suppressed the probability of resistance for all ETs, whereas that of ERY resulted in high probability of resistance (>75%) for at least one of the single-allele mutants for seven out of the eight trajectories. All treatment combinations eradicated the WT (black) and suppressed the double-allele (gxpy, red) subpopulations.

Probability of De Novo Resistance Evolution Varies between Trajectories.

Next, to examine the effects of CS-informed treatments using lower antibiotic concentrations due to collateral effects, we simulated combination treatments with Css ranging from 0.25 to 4× MIC (i.e., −2 to 2 log2 MIC) of the WT population. Our simulations indicated that eradication of WT bacteria was primarily observed for dosing schedules with Css of at least one of the antibiotics being equal to or larger than the MIC of the WT population (Fig. 5 and ). At these concentrations, a clear impact of collateral effects was observed. Specifically, subpopulations with CS were suppressed while subpopulations with CR were associated with a high probability of resistance evolution.

Fig. 5.

Treatment outcome of CIP in combination with LNZ. The treatments were evaluated on eight different trajectories leading to high-level FQ resistance (see Fig. 4). (A) Simulation outcomes of CIP in combination with LNZ with Css around the MIC of the WT population (Css 0.25 to 4× MIC). Concentration dependency related to WT strain eradication and resistance evolution. Each point represents a bacterial subpopulation from a specific ET (indicated by color). CEs are indicated by shape and line type, where solid lines with round points represent no CE, dashed lines with triangles pointing up represent CR, and dotted lines with triangles pointing down represent CS. (B) End-of-treatment probability of resistance (mutants) and treatment failure (WT) for different ETs treated with antibiotic doses resulting in an average steady state plasma concentration of CIP and LNZ equal to the 1× MIC of the WT. The presence or not of CEs are indicated above each specific mutant. The ET, i.e., the multistep alternative pathways for the de novo emergence of a particular high-level FQ-resistant double-mutant (gxpy), had an impact on treatment outcome, although it was not the sole determinant.

In addition, our simulations showed that the population dynamics and treatment outcomes for ciprofloxacin mono- and combination treatments with Css around the MIC of the WT population varied both within and between antibiotic treatments and ETs (). In particular, the fixation of each specific gx and py strain in the population varied between trajectories due to competition between these two arising clonal lineages. For example, when comparing our eight different ETs treated with 1× MIC for both ciprofloxacin and linezolid (Fig. 5), the fixation of specific gx, and therefore the emergence of resistance, largely differed depending on the concurrent py strain. The specific treatment schedule resulted in high probabilities of fixation for at least one FQ-resistant subpopulation (>95%) for six out the eight ETs. However, the probability of fixation for different strains varied across these trajectories, e.g., whereas the probability of fixation of pS794 was 0.032% when gS81F was the concurrent strain (ET 4), it rose to 98.5% when gE85G was the concurrent strain (ET 3). Overall, these findings suggest that strain-specific collateral and fitness effects as well as clonal competition have an impact on the dynamics of FQ resistance evolution.

Collateral and Fitness Effects Determine the Outcome of Antibiotic Combination Treatments.

The heterogeneous fixation of individual strains observed (Figs. 4 and 5 and ) prompted us to assess and disentangle the impact of collateral and fitness effects on this phenomenon. To this end, we simulated four different scenarios of treatment combinations with ciprofloxacin (1× MIC) and linezolid (1× MIC), ciprofloxacin (1× MIC) and erythromycin (1× MIC), as well as ciprofloxacin (0.5× MIC) and penicillin (1× MIC), corresponding to the combinations associated with large between-trajectory variability in our previous set of simulations. Specifically, the first scenario (scenario 1) did not include the experimentally derived collateral or fitness effects, while the remaining simulations included either the fitness (scenario 2) or the collateral effect (scenario 3), or both (scenario 4). The following two different model structures were used for the simulations of each scenario: one was a simplified structure where each ET was evaluated separately as described in Fig. 3, allowing a low degree of clonal competition (approach A) but high power to disentangle impact of collateral and fitness effects, and the other was an adapted structure where all strains used in our study were included simultaneously (free-for-all simulation; approach B), allowing for a higher degree of clonal competition and capturing a more real-life scenario compared to our eight-defined-trajectory approach.

When simulating each ET separately and in the absence of both collateral and fitness effects (scenario 1A), the ciprofloxacin-linezolid treatment (Fig. 6) resulted in a fixation probability of 47.4% for py strains for all ETs but ET 1, with a probability of resistance of 89.0%. Notably, ET 1 is the only trajectory including pD83N, which is the single-allele strain associated with the highest ciprofloxacin MIC (32 mg/L) in our experiments. Both the gx and gxpy strains were completely suppressed in these simulations. For the two other simulated antibiotic combinations (), gx and py strains had a similar probability of resistance (∼80%), whereas the gxpy strains had a probability lower than 1%. The latter suggests that due to the low antibiotic pressure in our simulations by the time a double mutant emerges via a two-step mutation, the single-mutant populations are large enough to suppress its fixation. The outcome of these simulations suggests that the majority of the heterogeneity observed in the previous simulations are due to the collateral and/or fitness effects, and not differences in sensitivity toward ciprofloxacin. In the free-for-all simulation (scenario 1B), we observed a similar high probability of resistance development during both ciprofloxacin-erythromycin (>80%) and ciprofloxacin-penicillin combination treatments (>70%) for all the single-allele strains tested (Fig. 6 and ). During ciprofloxacin-linezolid combination treatment however, the highest probability of resistance was due to the fixation of the pD83N strain (89.2%), followed by that of other py strains (pS79F, pS79Y, and pD83Y) with resistance probabilities close to 50%, while the other FQ-resistant strains had a 0% probability of fixation (Fig. 6).

Fig. 6.

Relationship between CE, relative growth rate, and probability of resistance under treatment with antibiotic combinations. (A) Simulation of the eight defined S. pneumoniae ETs leading to high-level FQ resistance (indicated by number). (B) Free-for-all simulation of resistance evolution in S. pneumoniae allowing for the simultaneous emergence of 16 different FQ-resistant subpopulations (4 gx, 4 py, and 8 gxpy). In both simulations, the treatment combinations consisted of CIP and LNZ with average steady state plasma concentrations equal to 1× MIC of the WT. Each unique strain is represented on the x-axis. The CE of LNZ is depicted using shapes, where a triangle pointing up is CR, a triangle pointing down is CS, and a circle represents no CE. The probability of resistance at the end of treatment is shown on the y axis. Colors indicate the relative growth rate of each mutant compared to the WT. Each panel column denotes an FQ-resistant subpopulation (gx, py, or gxpy) with the Top panels representing simulations where both collateral and fitness effects were included, Top Middle panels representing simulations based solely on the experimentally determined CEs, Lower Middle panels based solely on the experimentally determined relative growth rate, and the Lowest panel simulations where neither collateral nor fitness effects were included.

Reintroducing fitness effects in the absence of collateral effects and simulating each trajectory separately (Scenario 2A) revealed that only the strain with a fitness advantage (pS79Y; ET 3, 4, and 5) was fixed during the simulated ciprofloxacin-linezolid treatment (probability of resistance 95.8%) (Fig. 6). This differed from the other two treatments simulated for this scenario, where there was large variation between trajectories and specific strains, e.g., under ciprofloxacin-penicillin treatment, gS81F paired with pS79F (ET 4) was associated with a probability of resistance of 16.2%, while when paired with pD83Y (ET 2) or pS79Y (ET 7), the probability was 59.2% and 76.6% (), respectively. These intertrajectory differences relating to specific strains again highlight the importance of clonal competition in resistance evolution. During the free-for-all simulation (scenario 2B), only the pS79F FQ-resistant strain, which exhibited a fitness advantage, was associated with a high probability of fixation (>75%). When treated with the ciprofloxacin-linezolid combination (Fig. 6), all other subpopulations were fully suppressed, while the two other antibiotic combinations resulted in low levels (<25%) of FQ resistance due to the fixation of a few specific subpopulations (pS79Y, pD83Y, gS81F, and gS81Y) ().

The ciprofloxacin-linezolid treatment simulations informed solely by the collateral effects (scenario 3A) showed a decrease in the emergence and fixation of pD83Y (ET 2) (Fig. 6), the only py strain exhibiting CS toward linezolid, compared to the simulations without any effects (scenario 1A). All gx and py strains exhibiting CR were associated with a high probability of resistance (86 to 97%). When such a gx strain was combined with a py strain not showing any collateral effect (e.g., ET 4, gS81F and pS79F), the py-associated probability of resistance was decreased compared to when the same stain was paired with a gx showing no CR (e.g., ET 3, gE85G and pS79F). Simulations of the same scenario using a combination of ciprofloxacin and penicillin showed that for all single-allele strains exhibiting CS to penicillin, the probability of resistance was reduced compared to scenario 1A. However, no ETs included a pair of single-allele strains where both show CS. Therefore, we are unable to elucidate if reciprocal CS could suppress all resistance development. During the free-for-all simulation, the fixation of specific FQ-resistant subpopulations differed between the treatments (Fig. 6 and ). However, the de novo emergence of single-allele subpopulations exhibiting CR was associated with a high probability of resistance (>85%), while the emergence of subpopulations exhibiting CS was associated with resistance suppression for all three antibiotic combinations.

Finally, when we compared the treatment outcome of scenarios 1 to 3 to scenario 4, with the latter taking into consideration the fitness and the collateral effects, we showed that both effects impact the probability of resistance (Fig. 6 and ). However, for the single-allele strains, the collateral effects seem to have a larger impact if present, as low relative growth rate combined with CR would result in a high probability of resistance. The probability of the double-allele strains fixing during ciprofloxacin and erythromycin combination treatment was suppressed in the presence of CS while clearly affected by both collateral and fitness effects in the presence of CR (). Overall, our findings suggest that both collateral and fitness effects can impact resistance evolution. In contrast, the simulation using the free-for-all model revealed that the fitness effect was the main driver of the resistance development during treatment with antibiotic combinations (Fig. 6 and ). However, due to the strong fitness advantage of pS79F observed here, which outcompete the other subpopulations, the impact of the fitness cost might be masked by the clonal suppression elicited by this subpopulation.

Discussion

Our aim in this study was to integrate experiments and simulations to evaluate CS-informed multidrug treatment strategies to prevent FQ resistance in S. pneumoniae. First, we characterized the collateral and fitness effects of clinically relevant FQ resistance alleles on S. pneumoniae susceptibility to a wide range of antibiotics (19). We observed extensive collateral effects toward the majority of the antibiotics tested, especially antibiotics inhibiting protein synthesis, consistent with earlier results derived from laboratory-evolved strains of Enterococcus faecalis (7), Pseudomonas aeruginosa (8, 9), Escherichia coli (5, 13), and Acinetobacter baumannii (14). CS was much more predominant than CR and the magnitude of both effects was moderate, with the exception of clindamycin and gentamicin, which is in line with previous findings in E. coli (13). Most FQ-resistant mutants also exhibited CS toward tetracycline and chloramphenicol, consistent with results in E. faecalis (7). FQ resistance mutations in gyrA and parC result in three-dimensional structural changes of the DNA gyrase and the topoisomerase IV that are associated with the modification of DNA topology altering the global supercoiling of the bacterial genome (25). In Salmonella enterica serovar Typhimurium (26), this leads to an altered transcriptome and the subsequent global reprogramming of gene expression (27, 28), which in turn leads to collateral effects. We hypothesize that similar gene expression reprogramming, involving, i.e., the altered expression of stress response pathways, may also be responsible for the extensive collateral effects we observed in S. pneumoniae, although this remains to be elucidated in future studies.

In addition to consistent CS responses, we also observed significant heterogeneity in the CS responses of different gyrA and parC genes and alleles, even for FQ-related substitutions at single amino acids in the same gene. Most of these heterogenous effects resulted in changes to the magnitude of CS effects, but in a few cases, they also influenced the direction of the response. As yet, the cause of these differences, and of collateral effects overall, remain unclear. Because our design used the transformation of gyrA and parC alleles from clinically resistant isolates of S. pneumoniae (19), our transformants include known QRDR mutations plus a small number of secondary mutations that are not known to be associated with FQ resistance (). Although this allelic diversity reflects the natural diversity at these loci, we cannot fully exclude the influence of these non-QRDR mutations; subsequent studies using single-site mutagenesis, together with more detailed studies of transcriptional responses, will be valuable to more fully disentangle these possible influences. Overall, although our experimental results suggest that CS-informed combination treatments may be promising to limit the evolution of FQ resistance, heterogeneous collateral effects may necessitate systematic screening of allelic identity to establish appropriate antibiotic combinations for the successful selection inversion strategies.

Next, using a modeling approach, we combined our experimental findings on collateral and fitness effects, with antibiotic pharmacokinetics and pharmacodynamics in order to understand multistep resistance evolution in S. pneumoniae under different treatment conditions. Despite heterogeneous collateral responses for the antibiotics tested, our results simulating clinical dosing schedules for ciprofloxacin monotreatment and in combination with erythromycin, linezolid, and penicillin showed a clear benefit of combination treatment over ciprofloxacin monotreatment, often successfully suppressing the emergence of FQ resistance and clearing the infection. However, given that the clinical dosing regimens for the three antibiotics that were coadministered here with ciprofloxacin generally resulted in a plasma concentration greater than the MIC of the WT (2.6 to 594× WT MIC) and that the magnitude of the CS was generally moderate (−0.58 median fold decrease in the MIC) for these specific antibiotics, the beneficial outcome is likely driven by not only the collateral effects but also the intrinsically increased bacterial killing associated with the addition of the second antibiotic.

In silico treatment combinations with Css for each of the coadministered antibiotics around the MIC of the WT strain (Css 0.25 to 2× MIC; far below the clinical dosing concentrations) inhibited the fixation of FQ resistant subpopulations, and further eradicated the infection, when FQ resistance conferred CS to the one of the coadministered antibiotics. These results indicate the potential effectiveness of CS-informed treatment combinations even with lower dosages for each of the coadministered antibiotics in eradicating infections while minimizing the risk for resistance development. This, in turn, suggests the exploitability of CS-informed treatment combinations especially for antibiotics having a narrow window between their effective doses and those doses at which they give rise to adverse toxic effects, which is in line with previous research (29).

In addition to the collateral effects, our simulations using subclinical dosing revealed that fitness effects had an impact on resistance development during treatment. In fact, antibiotic-resistant clonal variants encoding distinct mutations and exhibiting varying collateral and fitness effects emerge during antibiotic treatment. The clonal competition among these variants determined by their collateral and fitness effects, in turn, shapes the population dynamics and therefore the overall treatment outcome. This highlights the need to consider fitness, in addition to collateral effects, when attempting to design effective CS-informed treatments.

Translating a complex system into a simplified model requires a number of assumptions, which should be taken into consideration when interpreting our simulation results. Our pharmacodynamic model includes parameters derived from early phase data from several in vitro time-kill literature studies, using different S. pneumoniae strains (30–35). This approach allowed us to obtain typical effect parameters of S. pneumoniae rather than strain-specific effects. We assumed that the early phase data represent the antibiotic-mediated killing of one homogenous population, thus ignoring the possibility of the observed killing rate being affected by the growth of less susceptible subpopulations. Furthermore, we assumed that there were no pharmacokinetic interactions between the antibiotics tested. Only FQ resistance was considered in the model while resistance evolution and possible collateral effects of the second antibiotic were ignored. The model could be further expanded to include such information, thus allowing us to investigate the impact of reciprocal versus nonreciprocal collateral effects, as well as sequential treatment regimens to exploit these reciprocal effects.

In conclusion, we showed that the FQ resistance allelic identity can have a pronounced impact on the fitness and collateral effects to different antibiotics in S. pneumoniae and that they can strongly influence the multistep evolution of high-level FQ resistance under antibiotic mixing regimens. Our results highlight the importance of estimating CS for different resistance mutations to the same antibiotic; just as different mutations impart different fitness costs, so too can they drive distinct collateral responses, in turn affecting their probability of fixation during combination treatments. Despite heterogeneity in collateral effects, our results suggest that conserved CS can be exploited to develop CS-informed antibiotic treatment combinations to successfully eradicate infections while suppressing the de novo emergence of resistance during treatment.

Materials and Methods

Bacterial Strains, Growth Conditions, and Media.

Six S. pneumoniae strains encoding gyrA and/or parC alleles that confer FQ resistance were provided by the CDC Streptococcus Laboratory and used as donors for all subsequent transformations with S. pneumoniae D39 strain as the WT recipient strain. S. pneumoniae ATCC 49619 was used as a quality control for antimicrobial susceptibility testing. All strains used in this study and their relevant characteristics are listed in . Strains were routinely grown either on tryptic soy agar (BD) supplemented with 0.5% wt/vol yeast extract (BD) and 5% vol/vol sheep blood (Sanbio) (TSYA) or in tryptic soy broth (BD) supplemented with 0.5% wt/vol yeast extract (TSYB) at 37 °C under 5% CO2 for 18 to 24 h, unless otherwise mentioned. The 14 antibiotics used in this study were prepared from powder stock according to manufacturers’ recommendations, stored at −20 °C or −80 °C, and are listed in .

Construction of FQ-Resistant Mutants.

Amplified fragments from genomic DNA of the FQ-resistant S. pneumoniae strains encoding each of the non-WT gyrA or parC gene allele and 3 Kb of its flanking regions necessary for the integration by double crossover were used to transform the S. pneumoniae D39 strain. All the primer sequences used in this study are listed in . Transformation was performed using a saturating concentration of 1 μg/mL of amplified DNA and 0.1 μg/mL synthetic competence stimulating peptide 1 (CSP-1; GenScript), as previously described (36). Putative single-allele gyrA (gx) and parC (py) mutants were selected on TSYA plates supplemented with 0.5 mg/L sparfloxacin (Santa Cruz Biotechnology) or 4 mg/L ciprofloxacin (Acros Organics), respectively. Double-allele (gxpy) mutants were generated by transforming the generated parC mutants with the respective gyrA non-WT alleles and selecting on TSYA plates supplemented with 12 mg/L ciprofloxacin (Acros Organics). GyrA and parC allelic replacements in the transformants were confirmed by PCR and Sanger sequencing.

Whole-Genome Sequencing.

Genomic DNA was extracted from the S. pneumoniae D39 WT and all FQ-resistant mutants using the DNeasy blood and tissue kit (QIAGEN) according to the manufacturer’s recommendations for Gram-positive bacteria. Quality control of purified genomic DNA was performed using the Quant-iT double-stranded DNA broad-range assay kit (Thermo Scientific) and 1% agarose gel electrophoresis. Genomic DNA was fragmented by ultrasound on Covaris S/E210 according to the manufacturer’s recommendations and was subsequently used for library preparation as previously described (37). Whole-genome sequencing was commercially performed using 100-bp paired-end libraries on a BGISEQ-500 platform at BGI Tech Solutions or on a NextSEq 2000 platform at MiGS. High-quality filtered reads were mapped to the reference genome of S. pneumoniae D39V available in GenBank (accession number NZ_CP027540.1) in order to identify and annotate genetic differences found between our FQ-resistant mutants and their parental D39 WT strain using the open-source computational pipeline breseq with the default parameter settings (38).

Antimicrobial Susceptibility Testing.

MICs of all isogenic strains for 12 clinically relevant antibiotics and ciprofloxacin () were determined in triplicate by broth microdilution according to European Committee on Antimicrobial Susceptibility Testing (EUCAST; http://www.eucast.org) and International Organization for Standardization (ISO) 20776-1:2006 guidelines with two minor modifications. A 1.5-fold testing scale was used to include the standard twofold antibiotic concentrations and their median values, and the cation-adjusted Mueller-Hinton broth (CAMHB; BD) was supplemented with 100 U of catalase (Worthington Biochemical Corporation) instead of 5% lysed horse blood. The lowest antibiotic concentration where no turbidity was observed was scored as the MIC.

Static Time-Kill Assays.

Static time-kill assays were performed in triplicate using a starting inoculum of the parental D39 WT strain equating to ∼106 CFU/mL in 15 mL of CAMHB supplemented with 5% lysed blood according to EUCAST and ISO 20776-1:2006 guidelines. Antibiotic-free cultures and cultures with either the addition of antibiotics as monotherapy or in combinations at final concentrations equal to the Css of each antibiotic were incubated at 37 °C under 5% CO2 continuous agitation (100 rpm). At set time points of 0, 3, and 6 h postinoculation, 100-μL samples were collected, serially diluted, and plated on TSYA plates. After 48 h of incubation at 37 °C and 5% CO2, colonies were counted for viable cell titer determination.

Collateral Effect Determination.

Collateral effects were determined as the decrease (CS) or increase (CR) in the MIC of each mutant relative to the isogenic WT D39 strain, and their magnitude was determined as the log2-scaled fold change in the MICs between each mutant and the WT D39 strain (Eq. ), as previously described (7).

Where CE represents the collateral effect of a mutant (gx, py, or gxpy), MIC its corresponding MIC, and MIC the MIC of the D39 WT strain.

The conservation of the collateral effects was assessed based on the CS50 and CR50 thresholds, defined respectively as CS and CR effects occurring in more than 50% of the mutants tested (13).

Growth Rate Measurements.

The in vitro fitness of the WT and each mutant was determined in triplicate using the Malthusian growth model (39). Bacterial inocula of ∼108 CFU/mL were diluted 100-fold into 3 mL of fresh prewarmed TSYB medium. A total of 200 µL of each diluted culture was loaded in wells of a 100-well honeycomb plate and incubated at 37 °C. Subsequently, the optical density at 600 nm was measured hourly using a Bioscreen C Reader (Thermo Scientific), with 5 s of shaking before reads. Growth rates were calculated based on the slope of the line that fitted points displaying log-linear growth. Relative growth rates were calculated by dividing the generation time of each D39 mutant by that of the mean growth rate of the WT D39 strain.

Statistical Analysis.

Statistical analyses were performed using R (version 3.6.3). We performed hierarchical clustering to identify relationships between strains, and between antibiotics, relating to collateral effects using the complete linkage method within the stats R package. We evaluated between- and within-group differences in mean relative growth rate with a one-way ANOVA and a post hoc Tukey pairwise comparison with multiple correction when significant differences were identified (α = 0.05).

Mathematical Modeling.

We performed mathematical modeling of bacterial growth dynamics to evaluate the effect of different FQ mono- and combination treatments (antibiotic mixing). We applied the model to eight different trajectories and evaluated treatment outcomes and the probability of resistance. The model included clinical pharmacokinetics of antibiotics used in the treatment schedule, bacterial growth rates, antibiotic-mediated killing, collateral effects, and stepwise FQ resistance development (Fig. 3). The treatment simulations were conducted for ciprofloxacin monotreatment and in combination with erythromycin, linezolid, and penicillin, which were chosen because of the availability of literature-derived experimental data (30–35) needed for the estimation of pharmacodynamic parameters. We incorporated the experimentally measured relative fitness and collateral effects observed for the strains designed in this study.

The bacterial model was comprised of a four-state stochastic hybrid ordinary differential equation (ODE) model, where each state represents a bacterial subpopulation. Included S. pneumoniae subpopulations were the FQ-sensitive D39 WT (gyrAwt::parCwt), the resistant gyrA mutant (gx; gyrAx::parCwt), the resistant parC mutant (py; gyrAwt::parC), and the resistant gyrA and parC double mutant (gxpy; gyrAx::parC). Evolution of de novo resistance occurred in a stepwise manner with the emergence of either of the low-level FQ-resistant single-allele mutants (gx or py), followed by the emergence of their corresponding high-level FQ-resistant double-allele mutant (gxpy). The mutation events were modeled by a stochastic process calculating the specific number of bacteria mutated for each time interval τ, which was converted to a transition- and interval-specific mutation rate ( using a τ of 1 h. This mutation process was based on a binomial distribution B with a mutation probability equal to the gyrA and parC mutation frequency per hour (μ) (40), as shown in Eq. , and was identically implemented for each state transition.where n is the integer number of bacteria of subpopulation z at start of τ and V is the volume of the infection site.

Each bacterial state (z) included a logistic growth model for which the state-specific maximal growth rate () was a composite parameter of the literature derived growth rate of D39 (41) and our experimental growth rate ratio between mutant and WT. The systems growth capacity limitation parameter (B) and the bacterial concentration at the start of the infection were obtained from an in vivo model of bacteremia (24). The volume available for infection represented a human blood volume (42). All system specific parameters are stated in .

Antibiotic-mediated killing was implemented for each state z () according to Eq. , which includes the subpopulation-specific antibiotic effect of AB1 (), AB2 (), and a term describing the PD interaction between these drugs (), where significant interactions were identified based on our performed static time-kill assays.

The antibiotic-specific effects were described according to a previously described model (43) where the effect of the ith antibiotic on bacterial state z ( related to the unbound concentration () according to Eq. .

Where represents the maximal killing rate, the shape factor of the concentration-effect relationship, the state specific maximal growth rate, and the state specific MIC.

Pharmacodynamic model parameter values were obtained by fitting the model to digitized early phase S. pneumoniae experimental in vitro time-kill data using the nlmixr package (version 1.1.1-7) within the statistical software R (version 3.6.3). Growth in the absence of antibiotic and antibiotic-mediated bacterial killing were modeled sequentially. Initially, bacterial maximal growth rates were obtained by fitting an exponential growth model for all the antibiotic-free growth data. A log-normal random effect was included on maximal growth rate to account for between-experiment variability. The estimated individual maximal growth rates were included as a covariate in the subsequent fitting of the antibiotic-mediated killing. Separate models were fitted for each antibiotic to obtain the drug-specific parameters and . Antibiotic-mediated killing was incorporated separately for each antibiotic on each individual bacterial state with the corresponding experimentally determined MICs ().

Data generated through the interaction study were fitted with nonlinear models including exponential bacterial growth and antibiotic-mediated killing implemented according to Eq. .

Where B represents the bacterial density, B the inoculum, k the bacterial maximal growth rate, the effect size of antibiotic 1 (AB1), the effect size of antibiotic 2 (AB2), and the effect size of the interaction between AB1 and AB2. Here, antibiotic AB1 and AB2 were binary variables indicating the presence (1) or absence (0) of the respective drugs. To facilitate the translation of the integration to the PK-PD model framework, a relative interaction term () was derived according to Eq. .

The PD interaction was subsequently integrated to the framework according to Eq. .

Where is given by Eq. .

The subpopulation-specific rates of change for bacterial density are shown in the Eqs. –. These rates are dependent on the bacterial density of the specific subpopulation, the subpopulation-specific net growth (k), the antibiotic effect (E), and the mutational transition(s) (k,).

Where k is given by Eq. .

The pharmacodynamic model was linked to a previously published pharmacokinetic model of ciprofloxacin (44, 45) and a pharmacokinetic model of the respective antibiotic selected for combination treatment () (46–49). The pharmacokinetic models were used to simulate pharmacokinetic profiles of an individual with a body weight of 70 kg and a creatinine clearance of 75 mL/min/1.73 m2. The estimate of the 95th percentile of parameters related to clearance and volume were used if the interindividual variability was reported. The fraction unbound (fu) was used to obtain the unbound concentration of each antibiotic.

The modeling framework was used to simulate two different sets of treatment scenarios of S. pneumoniae bacteraemia, as follows: 1) using ciprofloxacin as monotreatment (500 mg b.i.d.) or in combination with erythromycin (600 mg b.i.d.), linezolid (600 mg b.i.d.), or penicillin (3 g b.i.d.); and 2) using a series of nonclinical treatment scenarios where the Css of each antibiotic was related to the MIC of the WT (Css 0.25 to 2× MIC). A third set of simulations was performed where we investigated and disentangled the impact of fitness and collateral effect on resistance development by removing one or both of the experimentally derived effects. For these simulations, we focused on the concentrations associated with large between-trajectory variability in the previous set of simulations. The selection of the antipneumococcal antibiotics included in these treatment scenarios was based on the availability of relevant time-kill data, and ciprofloxacin was used as a representative agent of the FQ class. Each of the treatment scenarios was simulated 500 times. The treatment outcomes were evaluated by assessing the probability of failure in eradicating the WT strain and/or the probability of resistance establishment. The failure in WT eradication and resistance establishment were defined as a bacterial density of WT and FQ-resistant mutant strains exceeding 104 CFU/mL at the end of treatment (2 wk), respectively, with the particular density corresponding to both the initial inoculum and the density of an established infection (24). The framework was used to simulate treatments of eight defined trajectories, each depicting the de novo multistep emergence of a particular high-level FQ-resistant double-allele mutant (gxpy), with either of its corresponding low-level FQ-resistant single-allele mutants (gx or py) as an intermediate step (Fig. 2).

The modeling framework was adapted to allow for a free-for-all simulation, where all ETs were included simultaneously. Thus, the adapted model constituted of four gx subpopulations, four py subpopulations, and eight gxpy subpopulations. Appropriate state transitions were included on gx and py, subsequently leading to the eight gxpy studied. This adapted model was used to further investigate the role of the collateral and fitness effects in the context of greater clonal competition.

Software and Model Code.

All model simulations were conducted in R (version 3.6.3), using the ODE solver package RxODE (version 1.0.0). The modeling framework, the adapted free-for-all version, and the associated code are available at Github (https://github.com/vanhasseltlab/PKPD-CSpneumo) (50).