Multi-step vs. single-step resistance evolution under different drugs, pharmacokinetics, and treatment regimens

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

Version of Record published: June 7, 2021 (This version)
Accepted Manuscript published: May 18, 2021 (Go to version)
Accepted: May 4, 2021
Received: October 18, 2020

1. Part of Collection
Evolutionary Medicine: A Special Issue

Edited by George H Perry et al.
Further reading

Article
Figures and data
Abstract
eLife digest
Introduction
Results
Discussion
Materials and methods
Data availability
References
Decision letter
Author response
Article and author information
Metrics

Abstract

The success of antimicrobial treatment is threatened by the evolution of drug resistance. Population genetic models are an important tool in mitigating that threat. However, most such models consider resistance emergence via a single mutational step. Here, we assembled experimental evidence that drug resistance evolution follows two patterns: (i) a single mutation, which provides a large resistance benefit, or (ii) multiple mutations, each conferring a small benefit, which combine to yield high-level resistance. Using stochastic modeling, we then investigated the consequences of these two patterns for treatment failure and population diversity under various treatments. We find that resistance evolution is substantially limited if more than two mutations are required and that the extent of this limitation depends on the combination of drug type and pharmacokinetic profile. Further, if multiple mutations are necessary, adaptive treatment, which only suppresses the bacterial population, delays treatment failure due to resistance for a longer time than aggressive treatment, which aims at eradication.

eLife digest

The rise in antibiotic resistance is threatening our ability to treat bacterial infections. Bacteria often evolve resistance by acquiring new genetic mutations during the treatment period. Understanding how resistance emerges and spreads through a bacterial population is crucial to prevent antibiotic drugs from failing.

Mathematical models are a useful tool for exploring how bacteria will respond to antibiotics and assessing the risk of resistance. Usually, these models only consider instances where bacteria acquire one genetic mutation that makes them virtually impervious to treatment. But, in nature, this is not the only possibility. Although some mutations do give bacteria a high level of resistance, numerous others only provide small amounts of protection against the drug. If these mutations accumulate in the same bacterial cell, their effects can combine to make the strain highly resistant to treatment. But it was unclear how the emergence of multiple mutations affects the risk of treatment failure and the diversity of the bacterial population.

To answer this question, Igler et al. devised a mathematical model in which each bacterium is able to mutate multiple times during the treatment period. The model revealed that if one mutation provides a high level of resistance on its own, the risk of bacteria surviving treatment is very high. But, if it takes more than two mutations to achieve a high level of resistance, the risk drops to almost nothing.

Igler et al. also found that the chance of bacteria evolving high enough resistance is affected by the type of antibiotics used and how fast the drug decays. With low-level resistance mutations, adapting treatment to maintain an acceptable number of sensitive bacteria as competitors for (a small number of) resistant bacteria was more effective at delaying treatment failure than trying to kill all the bacteria at once.

These findings suggest that adjusting the treatment strategy used for bacterial infections according to the proportion of low- and high-level resistance mutations could slow down the evolution of resistance. To apply these models in the real world, it will be important to measure the level of resistance conferred by single mutations. The type of models used here could also predict the response of other diseases that resist treatment, such as cancer.

Introduction

The rapid rise and spread of antimicrobial resistance severely curb the efficacy of drug treatments against pathogen infections. Treatment strategies are designed to maximize efficacy and limit toxicity, but their long-term applicability depends on the risk of resistance evolution (Nielsen and Friberg, 2013; Foo and Michor, 2009; Bonhoeffer et al., 1997). This highlights the importance of careful consideration of drug type, dose, and duration to guarantee the desired patient outcome whilst also reducing the risk of resistance evolution (Nielsen and Friberg, 2013; Drusano, 2004). In order to prevent drug resistance and preserve drug efficacy, treatment strategies should also be guided by an understanding of resistance evolution and the ability to assess its risk (Andersson et al., 2020; Read et al., 2011) – a task that is substantially facilitated through mathematical modeling (Nielsen and Friberg, 2013; Andersson et al., 2020; Read et al., 2011; Clarelli et al., 2020).

The main class of models used to predict drug action and treatment outcome are pharmacokinetic and pharmacodynamic (PKPD) models (Clarelli et al., 2020; Drlica, 2003; Abel Zur Wiesch et al., 2017; Chakrabarti and Michor, 2017), which describe the change in drug concentration over time (PKs) and the corresponding effect on a pathogen population (PDs). PKPD approaches are most commonly employed to study the efficacy of treatment without considering the possibility of resistance evolution, but coupled with bacterial population models, they can be used to investigate drug resistance evolution over time (Yu et al., 2018).

One severely understudied aspect in such approaches is that there are two fundamentally different patterns of de novo antibiotic (AB) resistance evolution: (i) ‘single-step’ resistance: a single mutation provides high drug resistance (Nielsen and Friberg, 2013; Drlica, 2003; Yu et al., 2018); or (ii) ‘multi-step’ resistance: the accumulation of several mutations of low individual benefit is necessary for high-level resistance (where high resistance here means higher than a given treatment dose). The availability of either pattern to a pathogen population under drug selection will affect the potential for resistance evolution and therefore the evolutionary dynamics in response to various treatment strategies.

We focus on resistance by de novo mutations as long-lasting infections such as those caused by Pseudomonas aeruginosa become hard to treat due to resistance evolving via mutations within the host during the course of the treatment (Oliver et al., 2000). Another example is tuberculosis (TB), arguably the infectious disease that has caused the highest number of deaths globally (Castro et al., 2020). During persistent TB infections, drug resistance evolves by chromosomal mutations while resistance by horizontal gene transfer (HGT) has not been observed (Castro et al., 2020). HGT is a common path to resistance in hospital-acquired infections and in cases of shorter treatment durations, as exemplified by Staphylococcus epidermidis infections that became resistant by acquiring plasmids carrying genes for linezolid resistance (Dortet et al., 2018).

In this study, we will comprehensively study the influence of the mechanistic pattern of resistance evolution itself (namely the benefits and costs of mutations) by considering ‘single-step’ resistance vs. ‘multi-step’ resistance. The emergence of mutations and their selection depend on an interplay between various treatment factors like drug type, dose, and treatment duration. These factors have been studied before to various extent in isolation (Nielsen and Friberg, 2013; Drusano, 2004), although rarely how their interactions shape resistance evolution (Martinez et al., 2012; Olofsson and Cars, 2007). We will first establish the existence of single-step and multi-step resistance patterns by reviewing evidence in the experimental literature, and then use the obtained parameter values to inform a stochastic PKPD model of multi-step resistance evolution, which we will explore under various treatment regimens.

We will establish the fundamental differences between evolutionary dynamics emerging from these two patterns in one specific treatment setting, but also explore the impact of various clinically relevant treatment strategies. First, we will compare two types of drugs, ABs and antimicrobial peptides (AMPs). AMPs are key components of innate defenses but also important new antimicrobial drugs, which work by disrupting the bacterial membrane (Zasloff, 2002; Mookherjee et al., 2020) – as opposed to ABs, which usually target intracellular structures. AMPs have been found previously to significantly reduce the risk of resistance evolution compared to conventional ABs (Yu et al., 2018; Lazzaro et al., 2020), partly explained by their distinct PDs like higher killing rates (Yu et al., 2018). Second, we will consider three different shapes of drug PKs, which are all clinically relevant (Nielsen and Friberg, 2013), but have rarely been compared in a systematic manner (Chakrabarti and Michor, 2017; Foo et al., 2012). These comprise fluctuating drug concentrations, increasing concentrations (which are then maintained at the highest level), and finally constant (which can be achieved in high-dose IV (intravenous) interventions). Third, as a number of recent studies have questioned the practice of ‘radical pathogen elimination’ (Read et al., 2011; Hansen et al., 2017; Hansen et al., 2020), we will compare aggressive elimination treatment with adaptive suppression – a strategy where the drug concentration is regularly adapted to the pathogen load – in a multi-step mutational framework (Hansen et al., 2017; Gatenby et al., 2009). Lastly, depending on the drug type, resistance evolution can be shaped either by chromosomal mutations or HGT, or both (van Hoek et al., 2011; Woodford and Ellington, 2007). Assuming a scenario where both options are available, we will study the relative importance of resistance resulting from de novo mutations as compared to HGT, which plays an important role in AB resistance evolution (van Hoek et al., 2011), although likely not as much in AMP resistance (Kintses et al., 2019). Taken together, this allows us to obtain an empirically informed modeling framework, which predicts evolutionary dynamics of ‘single-step’ resistance vs. ‘multi-step’ resistance in the context of drug type, PKs, and treatment strategies. We show how this framework provides critical insights into drug resistance emergence in clinically relevant treatment settings.

Results

Antibiotic resistance evolves via multiple low- or single high-benefit mutation(s)

Experimental studies document single target mutations as well as a sequence of mutational steps to drug resistance evolution in bacterial populations (Spohn et al., 2019; Chevereau et al., 2015; Melnyk et al., 2015; Lofton et al., 2013; Makarova et al., 2018; Kubicek-Sutherland et al., 2017), but no systematic review of these patterns has been conducted so far. Here, we only selected studies that report on both parameters, benefit and costs of resistance (see Materials and methods) (Spohn et al., 2019; Chevereau et al., 2015; Melnyk et al., 2015; Lofton et al., 2013; Makarova et al., 2018; Kubicek-Sutherland et al., 2017), in order to obtain a complete picture of the mutational effects. We define the benefit and cost of a mutation as an increase in the minimum inhibitory concentration (MIC) and as a reduction in growth (in the absence of drug), respectively. Despite differences in study setup and type of resistance mutations, we clearly found a wide range of effects, with a large number of benefits below typical clinical MIC breakpoint values (defining whether a bacterial strain is resistant), which are often 10xMIC or higher (EUCAST, 2020; Table 1, Figure 1) – hence likely necessitating multiple mutations for high resistance. The corresponding fitness costs range from almost none to a 25% reduction of the population growth rate and show a very weak positive correlation (R² = 0.07, p=0.09) with (log) benefit over all studies taken together (Figure 1B, Figure 1—figure supplement 1). In general, mutations seem likely to incur more costs than benefits. Notably, our literature search suggests a difference in mutational benefit available for two different antimicrobials: the average benefit of resistance mutations to AMPs is substantially lower than for commonly used ABs (Table 1, Figure 1—figure supplement 1). In the following, we use the correlation observed with these assembled benefit and cost values to inform a PKPD model that reflects the two patterns of resistance evolution.

Figure 1 with 2 supplements see all

Download asset Open asset

Pharmacodynamic (PD) model of single- and multi-step resistance.

(A) The PD curve relating bacterial net growth $ψ$ (which is between the maximal growth rate $ψ_{m a x}$ and the maximal killing rate $ψ_{m i n}$ ) to antimicrobial drug concentration (given in fold minimum inhibitory concentration (MIC)) illustrating a sensitive wildtype (black) and mutants with either small (blue) or large (orange) MIC increases per mutation (benefit), assuming a typical clinical drug dose of 10xMIC. Characteristic PD curves for antibiotic (AB) (solid lines) and antimicrobial peptide (AMP) (dashed lines) single-step resistance are compared (orange), demonstrating the typically steeper decrease ( $κ$ ) and lower $ψ_{m i n}$ observed with AMPs. (B) Shown are costs and benefits from various empirical studies (Table 1), each dot representing resistance mutations to a specific AB or AMP class. The cost of a mutation shows a very weak positive correlation with the log (benefit) (R²= 0.07, p=0.09). Blue and orange colors show multi- or single-step resistance benefits given the drug dose in (A). (C) Schematic of the PD model with several mutated subpopulations ( $M_{i}$ ), which grow with a cost *r(1-c)ⁱ*, determined by the number of mutations i, mutate with rate µ, and die at a constant rate $γ$ and a drug-specific rate $γ_{i} (b_{i})$ , which is dependent on the benefit $b_{i}$ conferred per mutation.

Figure 1—source data 1 Empirical data used to obtain Figure 1B. Calculated mean and standard deviation for mutational benefits and costs taken from Spohn et al., 2019 or Melnyk et al., 2015 for individual antimicrobial classes are given.: https://cdn.elifesciences.org/articles/64116/elife-64116-fig1-data1-v2.xlsx
Download elife-64116-fig1-data1-v2.xlsx

Table 1

Benefits and costs of drug resistance mutations from experimental studies reported for antibiotics (ABs) and antimicrobial peptides (AMPs), with small mutational benefits (likely giving rise to multi-step resistance patterns) given in blue and large ones (likely giving rise to single-step resistance patterns) in red, assuming a typical clinical drug dose of about 10× minimum inhibitory concentration (MIC) (Figure 1A, B).

Source	Drug type	Organism and evolution environment	Number and type of mutations	Benefit per mutational event	Cost per mutation event	Benefit measurement	Fitness measurement
Spohn et al., 2019 (benefits and costs calculated for individual AB classes are given in Figure 1—source data 1)	AMP	Escherichia coli K-12 BW25113 populations were evolved in minimal salt (MS) medium over 20 transfers every 72 hr at 30°C (~120 generations) with successively increasing dosages of the antimicrobial	197 independent mutational events (deletions, insertions, SNPs (Single nucleotide polymorphisms), and intergenic mutations) ~5.2 (±0.8) mutational events/genome	All MIC samples: 5.1 (±7.2) xMIC ~=1 xMIC/mutation Only the ones where costs were measured as well: 31.4 (±8.5) xMIC ~=6.0 xMIC/mutation	0.2 (±0.16) ~=0.04/mutation	Serial broth (MS) dilution medium; MIC was defined as OD600 < 0.05 fold-increase compared to ancestor	Continuous monitoring of optical density of liquid cultures (in MS); area under the growth curve from 1 to 24 hr, normalized by the wildtype (WT)
	AB		N.D. Assumed: ~5	All MIC samples: 72.0 (±2.8) xMIC ~=14.4 xMIC/mutation Only the ones where costs were measured as well: 120.0 (±2.8) xMIC ~=24.0 xMIC/mutation	0.47 (±0.29) ~=0.1/mutation
Melnyk et al., 2015 (benefits and costs calculated for individual AB classes are given in Figure 1—source data 1) (and references therein: synthesis of 24 studies)	AB	Borrelia burgdorferi, Campylobacter jejuni, E. coli, Enterococcus faecium, Mycobacterium smegmatis, Mycobacterium tuberculosis, Staphylococcus aureus, Streptococcus pneumonia (all pathogenic) in various environments	Single mutational events (in total 128 mutations)	96.4 (±19.8) xMIC	0.13 (±0.24)	MIC fold-increase to ancestor	Competitive fitness (via in vitro growth assays with WT)
Chevereau et al., 2015	AB	E. coli K-12 BW25113 knockout strains (Keio collection) were incubated at 30°C for 20 hr in rich media (LB) with various ABs	Single-gene deletions (3913 mutant strains) ~4 mutations, all types (day 10)	Mecillinam: 1.2 (±1.1) xIC₅₀ Trimethoprim: 1.2 (±1) xIC₅₀ Ciprofloxacin: 1.3 (±1.2) xIC₅₀ ~18.7 (±1.1) xIC₅₀	0.13 (±0.07) N.D.	IC₅₀ (in LB)	Reduction in growth rate (in LB); calculated from a linear fit of log(OD) in the range 0.022 < OD < 0.22; given relative to WT
Lofton et al., 2013	AMP	Salmonella typhimurium LT2 was passaged daily in rich media (refined LB) at 37°C with successively increasing AMP concentrations for 400–500 generations	SNPs and deletions 2–3 mutational events	LL-37: 1.5–6 xMIC ~=0.5–3 xMIC/mutation WGH: 12–48 xMIC ~=4–24 xMIC/mutation CNY100HL: 2–6 xMIC ~=1–3 xMIC/mutation	0.076 (±1.1) ~=0.025–0.038/mutation 0.11 ~=0.037–0.055/mutation 0.17 ~=0.057–0.085/mutation	MIC through serial broth (refined LB) dilution	Growth rate from OD600 measurements (in refined LB) in the range 0.02 < OD < 0.2; given relative to WT
Kubicek-Sutherland et al., 2017	AMP	S. aureus (MRSA) WT JE2 (DA28823; clinical isolate) was passaged daily in minimal media (MIEM) at 37°C with successively increasing AMP concentrations	1–3 (adaptive) amino acid substitutions	LL-37: 6.5 (±6.2) × survival (for this AMP benefit was measured not as fold MIC increase but fold bacterial survival) ~=2.2–6.5 xMIC/mutation WGH: 13.9 (±1.2) xMIC ~=4.6–14 xMIC/mutation PR-39: 2 (±1) xMIC ~=0.7–2 xMIC/mutation	<= 0 (no detectable cost, rather an advantage, was found in the media used in the evolution experiment)	MIC through serial broth (MIEM) dilution, except for LL-37 where kill curves were used (in tryptic soy broth, rich media)	Maximum growth rate in MIEM based on (exponential phase) OD600 measurements; given relative to WT
Makarova et al., 2018	AMP	S. aureus SH1000 and E. coli MG1655 were transferred for seven daily passages in rich media (Mueller Hinton Broth [MHB]) at 37°C with successively increasing AMP dosages	All types ~2.4 (±0.9) mutational events	Tenecin 1: 3.1 (±1.1) xMIC of wt ~=1.3 xMIC/mutation	0.25 (±0.13) ~=0.1/mutation	MIC through serial broth (MHB) dilution; defined as inhibition of visible growth after 24 hr	Maximum growth rate in MHB based on OD600 measurements
	AMP			Mean: 4.5 (±5.3) xMIC/mutation	Mean: 0.04 (±0.03)/mutation
	AB			Mean: 28.0 (±7.1) xMIC/mutation	Mean: 0.10 (±0.07)/mutation

The PD model

We mainly investigated the rise of de novo resistance in a clonal, susceptible pathogen population, which is a common starting point for many clinically relevant infections (Balmer and Tanner, 2011), by extending a previously described stochastic PKPD model (Materials and methods) (Yu et al., 2018). Specifically, we considered not only a single resistance mutation, but the potential emergence of a sequence of mutations, with each mutation conferring a certain (additional) benefit and cost (Figure 1). The number of mutations needed for ‘full’ resistance depends on the applied drug dose, but generally low mutational benefits are more likely to necessitate multi-step resistance evolution. To compare scenarios where a single mutation is sufficient to scenarios where several mutations have to arise in one cell, we ran the simulations over a range of mutational benefits (2–100 xMIC) – and their correlated fitness costs (Table 1, Figure 1) – in combination with various drug doses (0.5–100 xMIC). Hence, the minimum number of mutations necessary for resistance was predetermined (Figure 1—figure supplement 2), and we investigated how this affects the potential for pathogen survival and mutational diversity under various treatment strategies (PKs) and for two different antimicrobials (PDs) as described below. Competition between the various mutant subpopulations was modeled by imposing a carrying capacity for bacterial growth and very low turnover as soon as this capacity is reached.

Multi-step resistance patterns show lower risk of treatment failure and lower genetic diversity

First, we determined the probability of treatment failure by simulating change of the pathogen population over 200 hr under treatment with drugs (PD) parameters typical for bactericidal ABs (Yu et al., 2018; Supplementary file 1) being applied once every 24 hr (PK). We assumed that the pathogen population initially consists of completely susceptible bacteria and defined a treatment as failed if the pathogen population was not eradicated after 200 hr. We found that the probability of treatment failure was always close to 1 for single-step resistance evolution, but decreased rapidly if multiple mutations were required. Notably, already if three mutations were necessary to overcome the applied dose, the probability of pathogen survival approached 0 (Figure 2A, Figure 1—figure supplement 2). The qualitative picture of these results was not dependent on the specific cost-benefit correlation that we are assuming for most of our simulations (Figure 2—figure supplement 1).

Figure 2 with 6 supplements see all

Download asset Open asset

Resistance evolution with single- and multi-step patterns for peak pharmacokinetic (PK).

(A) Treatment failure probability (measured at the end of the treatment period) and (B) mutational diversity (calculated over the whole treatment period) are shown for two different antimicrobial classes (antibiotics (ABs), left; antimicrobial peptides (AMPs), right) for different combinations of mutational benefits (xMIC) and drug doses (xMIC). The diagonal line shows where the benefit per mutation is exactly equal to the given drug dose and separates single-step (ss, lower orange triangle), where one mutation gives a benefit higher than the applied dose, from multi-step (ms, upper blue triangle) resistance, where more than one mutation is needed for the accumulated benefit to match the drug dose. The arrows indicate the decrease in diversity either due to increasing extinction (too many mutations are needed for survival) or due to increasing cost per mutation (costs are too high compared to the additional benefit). A representative example of subMIC mutational diversity is shown magnified below the plots in (B). MIC: minimum inhibitory concentration.

One aspect of resistance evolution that is especially important when considering multiple mutations is the mutational diversity that arises in the pathogen population: high genetic diversity (here meaning diversity in the resistance phenotype) increases the probability that some individuals will be able to survive a given environment – such as treatment with other drugs (Castro et al., 2020) – and increases the adaptive potential overall (Van Egeren et al., 2018). Using the Shannon index to determine the highest mutational diversity obtained in the population over the treatment period, we clearly observed higher diversity with single-step than multi-step resistance evolution (Figure 2B), even if we increased the mutation rate proportionally to the number of mutations required (Figure 2—figure supplement 2). It can be shown analytically that a mutant strain can invade at the mutant-free equilibrium if the death rate of the sensitive strain is higher than the death rate of the mutant, where the mutant death rate is a combination of intrinsic and drug-induced death as well as the mutational cost (Materials and methods). The observed higher diversity with single-step patterns seems counterintuitive as the need for multiple mutations should increase diversity (Figure 2B), but can be explained as follows: at high drug doses and low benefits, this effect is due to extinction that effectively reduces genetic diversity, while at low doses and high benefits, high mutational costs inhibit the build-up of diversity. These findings agree with an experimental study showing that resistance alleles with low costs are favored (Wichelhaus et al., 2002).

Consistently lower treatment failure with multi-step resistance for various PKs and PDs

Our results clearly show less resistance if multiple mutations are necessary, but the relative importance of the number of resistance mutations compared to other treatment considerations like the dose-response profile of a drug (PD) (Yu et al., 2018; Spohn et al., 2019) or the administration mode (PK) required further investigation. Hence, we compared three different PKs: ‘peak’ (fast absorption and exponential decay), ‘ramp’ (slow, linear absorption and no decay), and ‘constant’ (immediate absorption and no decay) (Figure 3A). Whereas constant PKs distinctly lowered the probability of treatment failure and the emergence of mutational diversity, peak and ramp PKs showed similar magnitudes of resistance evolution (Figure 3B, C, Figure 3—figure supplements 1 and 2). However, ramp PKs lead to more than twice the mutational diversity with multi-step resistance patterns (Figure 3—figure supplement 1), which suggests that treatment failure and pathogen diversity are connected in a non-trivial manner: while higher mutational diversity increases the risk of resistance evolution, neither its presence nor absence is obviously predictive of the treatment outcome (Figure 2, Figure 3—figure supplements 1 and 2).

Figure 3 with 7 supplements see all

Download asset Open asset

Resistance evolution patterns with different pharmacokinetics (PKs).

(A) The three PKs used in the model are shown over time (in hours) for the same peak drug concentration (xMIC). (B) The treatment failure probabilities and (C) mutational diversities are given for the three PKs from (A) and two antimicrobial drug classes (antibiotics (ABs) and antimicrobial peptides (AMPs)). Blue (orange) bars show averages calculated over the blue (orange) framed triangular areas of multi-step (ms) (single-step (ss)) resistance evolution shown in Figure 2, Figure 3—figure supplements 1 and 2. MIC: minimum inhibitory concentration.

The evolutionary dynamics can also be contrasted for different antimicrobial drugs, AMPs and ABs, by using two different PD parameter sets (Materials and methods, Figure 1—source data 1). Briefly, AMP treatments are characterized by higher killing rates, steeper dose-response curves (Figure 1A), and lower mutation rates than AB ones (Yu et al., 2018). Consistent with previous findings that AMPs lead to a lower risk of resistance evolution and a narrower mutant selection window (MSW) than ABs (Yu et al., 2018), treatment failure and mutational diversity was lower for AMPs with peak and constant PK treatments (Figures 2 and 3, Figure 3—figure supplement 2). Notably, in accordance with empirical studies (Andersson and Hughes, 2014), we generally see mutations accumulating at sublethal drug doses, but the maximal diversity is substantially lower in AMP treatments (Figure 2B, Figure 2—figure supplement 1).

Interestingly, the steeper dose-response curve of AMPs seems to make their resistance dynamics more sensitive to the shape of the PK than those of ABs (Figures 2 and 3, Figure 3—figure supplements 1, 2, and 4): in contrast to the other two PK profiles, ramp PKs lead to a drastic increase in treatment failure with AMPs, especially in multi-step scenarios (Figure 3, Figure 3—figure supplement 1). Accordingly, for ramp PKs, AMPs did not perform better and under some conditions even worse than ABs (Figure 3—figure supplement 4). By varying the ramp duration (or equivalently the rate of drug uptake), we found that there is an intermediate range (48–84 hr), which showed increased treatment failure with AMPs over ABs (Figure 3—figure supplement 5A). Paradoxically, while a narrow MSW generally hinders the emergence of numerous mutations in the population, for ramp PKs it can lead to optimal selection conditions for the sequential emergence of increasingly higher resistance mutations due to the strong selection for the next mutation combined with sufficient time for its emergence. Hence, especially the risk of multi-step resistance is increased if AMPs are used with ramp treatments as compared to the other PKs (Figure 3B, C). The broader selection windows in the presence of ABs, on the other hand, overlap and resistance mutations are less strongly favored (Figure 3—figure supplement 5B). Overall, the number of resistance mutations was the main determinant of treatment outcome, but we also found a complex dependence on PK and PD characteristics.

This complexity in resistance determinants raises the question in how far the type of drug action influences treatment outcome. Specifically, antimicrobials can have bactericidal action (which we were modeling so far, through a drug-dependent death rate), but they can also act bacteriostatically, that is, decreasing bacterial growth. We would expect bacteriostatic antimicrobials to slow down the rise of mutations in comparison to bactericidal ones as the acquisition of mutations is also coupled to bacterial growth. However, we find that this is only true for ramp and constant PK treatments (Figure 3—figure supplement 6A). Peak PKs allow for regrowth of bacterial cells due to drug decay, which increases bacterial survival and treatment failure, especially with multi-step resistance. Interestingly, mutational diversity only increased for AMP treatments (Figure 3—figure supplement 6B).

Multi-step resistance can lower the threshold for adaptive treatment application

The conventional treatment goal is to ‘eradicate’ the pathogen population, but it has been suggested that under certain conditions ‘mitigation’ could be a superior strategy (Hansen et al., 2017; Hansen et al., 2020; Gatenby et al., 2009), for example, if it is likely that a resistant subpopulation already exists at the beginning of the treatment. This strategy is called adaptive treatment as drug doses are adapted to keep the sensitive population as big as possible and the total pathogen burden below a given limit. (In practice, this is challenging as it requires monitoring of the pathogen burden and adjusting drug doses accordingly, which is difficult to implement even for measurements of total within-patient loads.) In adaptive treatment, the sensitive population provides a benefit by competitively inhibiting the resistant subpopulation, but also a risk by supplying mutational input (Figure 4). This trade-off creates a threshold for the size of the pre-existent resistant subpopulation above which adaptive treatment is more effective than aggressive ‘eradication’ in containing the infection (Hansen et al., 2017).

Figure 4 with 2 supplements see all

Download asset Open asset

Adaptive treatment threshold.

(A) The dependence of the threshold number of resistant bacteria cells M_m is given for different proportions of the risk population to the whole competitor population for single-step (orange stars), where the risk population is always equal to 1, and multi-step (blue lines) resistance patterns. Different benefits (and correlated costs) per mutation are shown as different color shades. (B) The minimum inhibitory concentration (MIC) increases with every mutation (given by *b_i*), but only an MIC above the given drug dose will lead to a reproductive number R₀ >1, that is, growth of the population (resistant cells). All other subpopulations serve as competitors, and the subpopulation one mutation away from resistance is the risk population.

Previously, the threshold for adaptive treatment was derived in a single-step resistance scenario (Hansen et al., 2017). When we incorporated adaptive treatment in our multi-step resistance framework (Figure 4—figure supplement 1), we found that the resistant subpopulation threshold above which adaptive treatment is more beneficial can be much lower in the multi-step scenario than in the single-step one (Figure 4A). This can be intuitively explained by the fact that all (partially) sensitive bacteria serve as competitors for fully resistant cells, but only the subpopulation one mutation away from being fully resistant constitutes the risk population (Figure 4B). Thus, with multi-step resistance there is a smaller population to supply resistant bacteria than with single-step resistance, changing the trade-off towards adaptive treatment. Additionally (in scenarios where adaptive treatment is favorable), the difference between adaptive and aggressive treatment in the duration until treatment failure can be several-fold larger for multi-step than single-step resistance patterns (Figure 4—figure supplement 2). Hence, assuming either single- or multi-step evolution could lead to considerably different treatment strategy assessments with regard to treatment failure through resistant pathogens.

HGT does not change the treatment failure probability

In addition to chromosomal mutations (Woodford and Ellington, 2007), antimicrobial resistance can be conferred through HGT (van Hoek et al., 2011), which could facilitate resistance in multi-step scenarios. To account for this possibility, we extended the model to allow for acquisition of a gene conferring full resistance, initially only at a low rate from the environment, and then at a density-dependent rate from other cells carrying the HGT gene (for assumption and implementation details, see Materials and methods). The HGT gene always provided immediate resistance to the applied maximal dose, regardless of the benefit or costs conferred by mutations. In order to compare the population dynamics of these two main antimicrobial resistance acquisition mechanisms, we assumed that resistance through mutations or HGT can be acquired independently of each other and that their effects are multiplicative.

Even though HGT carriers dominated the remaining pathogen population at the end of the treatment (Figure 2—figure supplement 3), the addition of HGT did not change the probability of treatment failure (Figure 2—figure supplement 4). This result holds true as long as the acquisition rate from the environment is lower than the mutation rate (this constraint is examined further in the Discussion). Consequently, initial rescue of the population is due to mutations – and therefore dependent on the magnitude of the mutational benefit – whereas HGT resistance is acquired later during the infection, after which it spreads rapidly.

Discussion

In this study, we compared the risk of drug resistance evolution patterns that either feature single resistance mutations with large costs and benefits or multiple steps involving mutations with smaller costs and benefits. We extended this comparison across a wide range of PD and PK profiles, which cover a multitude of antimicrobials and treatment strategies. We first showed that the single- and multi-step resistance patterns are relevant by gathering evidence of multi-step resistance patterns in the experimental literature (Table 1, Figure 1). While it is intuitive that drug resistance requiring more than one mutation will arise more slowly, we find that it can be a surprisingly strong inhibitor of resistance evolution and mutational diversity, depending on the drug class and administration route (Figure 3). We demonstrated that the number of mutations necessary for resistance strongly affects predictions of treatment outcome and optimality with regard to antimicrobial resistance – in a manner that is robust to variations in mutation rates and in the cost per mutation (Figure 2—figure supplements 2 and 5). Experimental support for our simulation results comes from studies reporting that mutational input limited to low benefits (Drlica, 2003) leads to decreased drug resistance evolution as compared to systems, in which high-benefit mutations are available (Allen et al., 2004). Moreover, limited access to high-benefit mutations seems to curtail MIC increase beyond a certain threshold (Chevereau et al., 2015).

The pattern of resistance evolution (single- and multi-step) is likely to be associated with the molecular mechanisms of resistance for a given antimicrobial: as an overall rule, the magnitude of the resistance benefit correlates with the mechanism of resistance, for example, efflux pumps yield low benefits, whereas specific drug target mutations yield high benefits (Hughes and Andersson, 2017).

Unfortunately, the specific mutations linked to the benefit and cost of mutations in our literature analysis (Table 1) are generally not known. Overall, however, MIC increase was low for drugs, which typically show unspecific resistance mechanisms via two-component systems or lipopolysaccharide modifications – as generally seen for AMP resistance (Lofton et al., 2013; Makarova et al., 2018; Kubicek-Sutherland et al., 2017) – and high for drugs with typical resistance via specific target modifications, as seen for some AB classes (e.g., rifampicin resistance via RNAP subunit mutations) (Goldstein, 2014). For most drugs, the prevalent resistance mechanisms are known (van Hoek et al., 2011); hence, this information can be used to determine drug and dosing regimens that minimize resistance evolution based on the inferred pattern of resistance evolution (i.e., using the probability that a single- or multi-step pattern is underlying resistance evolution). A recent study also suggests that resistance evolution in biofilms, which are often associated with clinical infections, is prone to occur through unspecific mechanisms, even if specific mechanisms are favored in planktonic cultures (Santos-Lopez et al., 2019).

Interestingly, the risk of resistance evolution does not seem to be related to the emerging mutational diversity in the population in a trivial manner as it is either limited by fast extinction or high mutational cost (Figure 2, Figure 3—figure supplements 1 and 2). Reducing mutational diversity is however a worthwhile goal in its own right as mutational diversity can increase adaptation by fixing more mildly deleterious mutations, which can then act as stepping stones for multi-drug resistance evolution (Van Egeren et al., 2018). Further, the diversity arising during the treatment period will help to determine if escalating the drug dose is expected to be beneficial or if, conversely, it would be detrimental because higher-resistance mutations are already present in the population and would be selected. The strength of this selection is determined by the MSW of the antimicrobial. Hence, the number of resistance mutations emerging over the treatment can be useful in estimating the width of the MSW – even though the diversity remaining at the end of the treatment will likely be lower (Figure 3—figure supplement 3).

We find that mutational diversity arises from a combination of selection pressure, bacterial growth, and fitness costs and cannot be predicted from the mutational benefit or the probability of treatment failure alone. Diversity is also shaped in unexpected ways by interactions between the drug type and drug concentration changes, making drug choice not only dependent on the PD characteristics, but also the specific drug PK in the target body compartment. Notably, this can lead to more favorable assessment of a specific drug application mode for one type of drug (e.g., AMPs for bolus drug application), but a different mode for another drug (e.g., ABs for drug infusions). While we mostly focused on the action of bactericidal drugs in this study, we note that purely bacteriostatic effects can lead to different trends for PK and PD influence on treatment failure and mutational diversity, for example, making peak PKs the least favorable drug administration route (Figure 3—figure supplement 6). The unexpected complexity in predicting which treatment strategies will minimize resistance evolution highlights the need of critically evaluating assumptions such as single-step resistance made in current PKPD models. The role of specific drug characteristics in resistance evolution is exemplified by the steepness of the PD curve, κ. By analyzing the selection coefficients for various treatments, we find that κ governs not only the size of the MSW (Yu et al., 2018; Chevereau et al., 2015), but generally shapes the selection pressure for resistance evolution in a qualitative manner. $ψ_{m i n}$ , the minimal bacterial growth rate, on the other hand, leads to substantial quantitative changes in selection pressure, meaning that κ and $ψ_{m i n}$ shape the form and strength of drug selection independently (Figure 3—figure supplement 5C). Ultimately, the interactions between PD and PK characteristics give rise to complex, and dynamic, fitness landscapes that are navigated by mutations of various benefit and cost sizes.

Interestingly, AMP-like drugs show considerably more resistance evolution with ramp PKs than in the other PK scenarios. This is noteworthy as natural AMP expression patterns in the producing organisms resemble ramp PKs (Johnston et al., 2014; Haine et al., 2008). This finding could suggest another reason why natural AMP production in cocktails is favorable (Zanchi et al., 2017) as AMP cocktails will limit the selection pressure and potential for resistance evolution to individual components. Intuitively, we would expect that a gradual increase in drug concentration would facilitate the rise of multiple mutations and indeed we find that ramp PKs lead to the highest probability of treatment failure and mutational diversity (Figure 3). However, a high probability of treatment failure is still mostly observed with high mutational benefits (Figure 3—figure supplement 1), that is, limited with the small-benefit mutations likely associated with multi-step resistance (Jochumsen et al., 2016). For clinical settings, our simulations caution that attention should be paid to the drug application mode when using AMPs. AMP-like colistin and daptomycin, for example, are typically applied as (short) IV treatments (Liu et al., 2011; Tsuji et al., 2019), which resemble peak PKs, and they are still active as last-resort drugs for multi-drug-resistant bacterial pathogens (Liu et al., 2011; Tsuji et al., 2019). Overall, our results agree with Yu et al., 2018 in that AMP treatment lowers resistance evolution and mutational diversity. This is particularly notable as multi-step patterns seem to be the common mechanism by which AMP resistance evolves (Table 1; Spohn et al., 2019; Jochumsen et al., 2016; Joo et al., 2016) – thereby suggesting another advantage over ABs, for which single- and multi-step evolution is common (Drlica, 2003; Weinreich et al., 2006; Wistrand-Yuen et al., 2018; Marcusson et al., 2009; Jin and Gross, 1988).

Unfortunately, distributions of mutational effects have rarely been characterized experimentally for drug resistance, and even then only for a single mutational step (Chevereau et al., 2015). We show, however, that this information is crucial as input for PKPD models to accurately predict resistance evolution and population diversity in response to drug treatment. Even between mutations involved in multi-step resistance to a single drug, benefit and costs of individual mutations are likely to vary (Figure 1B). In addition, epistasis in benefit and or cost magnitude can facilitate or preclude certain evolutionary pathways (Jochumsen et al., 2016). Both options can be easily included in our model, but empirical data in this regard is sparse, and we expect our main results with regard to PD and PK influence on single- and multi-step resistance to be robust to such changes. The empirical data that we used to inform our simulations did also not provide explicit information about potential compensatory mutations, which arguably can influence the dynamics of resistance evolution (Andersson and Hughes, 2010) – although likely in a very complex manner, as recent studies suggest (Dunai et al., 2019). According to our results, these mutations might even be a necessary means to allow multi-step resistance patterns to arise. If they emerge fast enough to compensate for the cost of the first mutation, they would increase the selection coefficient of this mutational subpopulation and thereby provide a stepping stone to high-level resistance. This might either be akin to crossing a fitness valley, if the first mutation does not provide a benefit, or it might facilitate climbing a fitness peak by making low-benefit mutations more favorable.

For many antimicrobial drugs, resistance evolution can not only arise through chromosomal mutations, but also by acquisition of resistance genes through HGT (van Hoek et al., 2011). Notably, our results highlight the importance of transfer rates as we find rescue of the pathogen population through HGT resistance only if the initial acquisition rate is higher than the mutation rate. HGT is not only dependent on the recipient population size but also on the donor population size, hence using typical experimentally measured conjugation rates of 10⁻¹¹–10⁻¹³ ml cell⁻¹h^-1 (Licht et al., 1999), environmental donors have to be more abundant than 10⁵ cells ml⁻¹ to be faster than chromosomal mutation rates (Rodríguez-Rojas et al., 2014) of 10⁻⁶, which might not always be the case at bacterial infection sites (Stecher et al., 2012). This implies either (i) that HGT resistance is acquired after chromosomal mutations, (ii) that HGT spreads mostly at sublethal drug doses, or (iii) that acquisition rates from a pre-existent pool of HGT carriers are high. Plasmid transfer rates are likely increased at low AB doses (Cairns et al., 2018), but generally they are highly variable, and even though they are biased towards spread between clone-mates, there seems to be no obvious correlation between transfer rates and genetic distance of donors and recipients (Dimitriu et al., 2019). Hence, determining the relative importance of resistance evolution through HGT or chromosomal mutations is difficult, but for specific drugs like AMPs, for which spread of HGT resistance from the gut microbiota seems to be low (Kintses et al., 2019), the risk of treatment failure is mainly shaped by the beneficial mutations available to the population.

Most of our results assumed a completely susceptible pathogen population at onset of treatment, as seen in many bacterial infections (Balmer and Tanner, 2011). However, the fast growth and high mutation rates can lead to significant heterogeneity in bacterial populations and we would expect this (neutral) heterogeneity to increase treatment failure, even with multi-step resistance patterns, by giving the population a ‘head-start’ in the accumulation of mutations. This is indeed what we see with our model when we start from a heterogeneous population, but we still find on average less than 50% treatment failure in each multi-step resistance scenario (considering various PKs and PDs) (Materials and methods, Figure 2—figure supplement 6, Figure 3—figure supplement 7).

When starting from populations that likely already contain resistance mutations, aiming for ‘mitigation’ (adaptive treatment) can be more effective in reducing resistance spread than trying to completely ‘eradicate’ the pathogen population (aggressive treatment). If multiple steps are necessary to obtain full resistance to the highest possible drug dose, we find that the threshold for choosing adaptive over aggressive treatment can be much lower than if only a single mutation were necessary (Figure 4). Additionally, in drug-free environments, we expect a lower frequency of resistant cells for multi-step patterns as it is less likely that neutral heterogeneity produces cells carrying all resistance mutations. Hence, the high competitive benefit is paired with a low risk for resistance evolution. Even though determination of the number and size of resistant subpopulations is very difficult in practice, this suggests that adaptive treatment is likely to be superior in containment of resistant infections for many drugs, for which multi-step patterns are the most common resistance mechanism. Further, the assumptions in our model are not specific to bacterial populations or antimicrobials, which makes them more broadly applicable to other drug treatments, like cancer therapy (Chakrabarti and Michor, 2017). Our results suggest a way forward to develop treatment strategies that – in addition to all other important considerations – explicitly account for the risk of drug resistance evolution.

Materials and methods

Key resources table

Reagent type (species) or resource	Designation	Source or reference	Identifiers	Additional information
Software, algorithm	R package adaptivetau	Johnson, P. Tau-Leaping Stochastic Simulation. R package version 2.2–3 (2019) (Johnson, 2019)	RRID:SCR_001905	The epsilon parameter (i.e., tolerance of relative rate changes in step size selection) was changed to 0.01 to increase accuracy
Other	Previously published datasets	Melnyk A, Wong A, Kassen R. The fitness costs of antibiotic resistance mutations. (2015) (Melnyk et al., 2015)	Dryad repository: http://doi.org/ 10.5061/dryad.5rc47

Literature review of costs and benefits of antimicrobial drug resistance mutations

Request a detailed protocol

We compounded a comprehensive set of experimental evolution studies (or reviews thereof) that measured both fitness costs (usually growth rate reductions in the absence of drugs) and benefits (usually increases in MIC) of AB or AMP resistance mutations within the same set of experiments. The studies used various bacterial species, including pathogenic isolates (see Table 1). From empirically measured data of sample replicates, we calculated costs as the arithmetic mean of 1-(relative fitness to wildtype) and the benefit as the geometric mean (due to the logarithmic scale of MIC/IC₅₀ measures) of MIC or IC₅₀ increase relative to the wildtype. (Note: Chevereau et al., 2015 used IC₅₀ instead of MIC but our calculation of IC₅₀ and IC₉₀ – which is likely very close to MIC – in their data gave a good correlation [R² = 0.45, p<0.001], which indicates that the benefits obtained from IC₅₀ measurements are comparable to ones obtained from MIC measurements.) As fitness measure, we considered only the measurements done in the same conditions (media and temperature) that was also used for experimental evolution, even if growth was also measured in different environments. We list the conditions of the various evolution experiments, MIC and fitness measurements in Table 1, with the exception of Melnyk et al., 2015, where we only list the eight different pathogenic strains used, as this paper synthesizes 24 different studies, grown under various conditions.

We obtained the type and average number of mutational events observed from supplemental data in most studies, but there was generally no possibility to link any individual resistance mutation with a specific cost and benefit. Hence, we divided the overall costs and benefits by the average number of observed (adaptive) mutations (i.e., mutations that were not observed in control lines), assuming that each mutation provides a similar share to the overall magnitude. As most studies have a very low number of mutational events linked to resistance, this assumption is not expected to lead to strong biases. Overall, the results from all of the studies gave only a very weak positive linear correlation between the log(benefit) and cost of a mutational event (Figure 1B). Mutations seem to be more likely to incur costs than benefits. This result is largely determined by the large data set from Spohn et al., 2019, which gives a very weak correlation between cost and benefit (Figure 1—figure supplement 1), similar to the data points from Melnyk et al., 2015. The dataset from Spohn et al., 2019 is the only one that fulfilled our criteria and directly compared AB and AMP mutational effects, which we summarize in Figure 1—figure supplement 1. The calculated benefit and cost per mutation for each individual AB and AMP in the Spohn et al., 2019 and Melnyk et al., 2015 data is given in Figure 1—source data 1.

PD model

Request a detailed protocol

We combined a PD model, which connects the growth of bacterial (mutant) subpopulations to antimicrobial drug concentration (Figure 1A; Nielsen and Friberg, 2013; Andersson et al., 2020; Read et al., 2011; Clarelli et al., 2020; Yu et al., 2018; Regoes et al., 2004), with a population model to predict the emergence of resistance mutations in individual bacterial cells.

In the population model, bacteria can grow up to a certain carrying capacity and can accumulate mutations during replication at a certain rate. In order to simulate de novo mutation emergence, we started most of our simulations from a completely susceptible population M₀, but we also ran simulations starting from neutral diversity (meaning that we ran the simulation for 50 hr without AB treatment and then started from the observed neutral heterogeneity) (Figure 2—figure supplement 6, Figure 3—figure supplement 7). We do not allow for reversion of resistance mutations, which has been found to be rare (Dunai et al., 2019) and likely does not play a role in multi-step resistance networks (over the time frame of a single treatment period). The population dynamics is captured by the following deterministic equations (which were implemented in a stochastic manner):

\frac{d M_{i}}{d t} = r ∙ {(1 - c)}^{i - 1} ∙ u ∙ M_{i - 1} ∙ (1 - \frac{M}{K}) + r ∙ {(1 - c)}^{i} ∙ (1 - u) ∙ M_{i} ∙ (1 - \frac{M}{K}) - (γ + γ_{i}) ∙ M_{i}

w i t h i = 0, 1, 2, \dots

M = \sum_{i} M_{i}

Here, M_i is the bacterial subpopulation carrying i mutations, r the wildtype growth rate (set to 1 in our simulations), c the cost of each mutation, u the mutation rate, K the carrying capacity of the system, $γ$ the natural death rate, and $γ_{i}$ the death rate caused by drugs (which captures the PD properties of a drug and the resistance level of the bacterial population via the mutational benefit).

The PD function

Request a detailed protocol

In our population model, cells die at a low intrinsic rate $γ$ , whereas death due to antimicrobials ( $γ_{i}$ ) is dependent on the properties of the antimicrobial applied, the benefit conferred by each mutation, and the PK profile. Specifically, $γ_{i}$ is calculated from the maximal and minimal growth rates $ψ_{m a x}$ and $ψ_{m i n}$ (note that $ψ_{m i n}$ can be negative in the presence of drugs, meaning that we generally consider bactericidal AB action), the (time-dependent) concentration of the drug a, the MIC of the mutation-free population (set to 1 in our simulations), the benefit b_i conferred by each mutation, and the sensitivity of the dose–growth relationship κ (the Hill coefficient or steepness of the curve):

γ_{i} = \frac{{(ψ}_{m a x} - ψ_{m i n}) ∙ {(a / (M I C ∙ b_{i}))}^{κ}}{{(a / (M I C ∙ b_{i}))}^{κ} - ψ_{m i n} / ψ_{m a x}}

ψ_{m a x} = r ∙ {(1 - c)}^{i} - γ

Considering bacteriostatic antimicrobial action can be achieved in our model by using a small $ψ_{m i n}$ value and incorporating antimicrobial effect into the growth, not the death term. Note that introducing antimicrobial action into the birth term here leads to density-dependent antimicrobial effects. This is not entirely unrealistic, considering persister bacteria, whose dormant state protects them from killing by ABs (Kussell et al., 2005). However, bacteriostatic action in itself would result in a high bacterial presence at the end of the treatment – even if bacteria are fully susceptible to the antimicrobial – as intrinsic bacterial death is very low. Hence, we incorporated an extrinsic removal rate $- γ_{c l} ∙ M_{i}$ , akin to immune system clearance of inert bacterial cells, with $γ_{c l} = 0.1 h^{- 1}$ being in a realistic range (Roach et al., 2017). The model for bacteriostatic drug action is then given by

\frac{d M_{i}}{d t} = (r ∙ {(1 - c)}^{i - 1} - γ_{i - 1}) ∙ u ∙ M_{i - 1} ∙ (1 - \frac{M}{K}) + (r ∙ {(1 - c)}^{i} - γ_{i}) ∙ (1 - u) ∙ M_{i} ∙ (1 - \frac{M}{K}) - (γ + γ_{c l}) ∙ M_{i}

w i t h i = 0, 1, 2, \dots

M = \sum_{i} M_{i}

Benefits and costs of mutations

Request a detailed protocol

The main interest of our study is the comparison of AB resistance evolution via ‘typical’ single mutations with complex, multi-step processes as shown in Jochumsen et al., 2016. The latter are characterized by a network of mutations of small benefits in multiple genomic resistance loci that create evolutionary pathways to high-level AB resistance (Jochumsen et al., 2016). We model this mutation accumulation via sequential acquisition of mutations with a certain benefit and cost, that is, decreases in drug-induced death and decreases in the maximum growth rate. Benefits and costs of each mutation were taken from the positive correlation that was observed with literature values (slope = 0.0087) – except for simulations testing the dependence of our results on this relationship, where we took a steeper correlation (slope = 0.0467) (Figure 2—figure supplement 1).

As benefits and costs are likely to vary, we also confirmed that our results are robust with regard to drawing benefits and costs of each mutation from a normal distribution. Similarly, we ran simulations with ‘peak PK’, where only the first mutational benefit/cost was fixed (i.e., deciding if a single- or multi-step pattern was necessary) and the other mutations were sampled from the whole range of benefits and costs obtained from the literature, independently of each other (Figure 2—figure supplement 5). We ignore the possibility of positive epistasis between these mutations (which would speed up resistance evolution as fewer mutations would be required for higher levels of resistance), as well as the possibility of negative epistasis, which would limit access to some mutations and the available pathways to resistance (thereby slowing down resistance evolution as ‘effective’ mutation rates might be lower than we assume in our model). Both of these processes are complex and not well understood, hence by ignoring these possibilities we aim to provide a more fundamental and intuitive comparison between single- and multi-step resistance evolution.

Resistance mutation rates were generally kept the same for each simulation run (i.e., regardless of the benefit magnitude). In reality, there might be more mutations available that provide a low benefit – which are likely to be less specific and therefore have a larger genomic target size, but using higher mutation rates for mutations with lower benefits and costs – which was done proportional to the number of steps needed to obtain resistance – did not change our results noticeably (Figure 2—figure supplement 2).

PK functions

Request a detailed protocol

In our simulations, we used three different PK functions to evaluate resistance evolution dynamics. ‘Peak PK’ describes the intake of a drug with a certain period τ, which is absorbed instantaneously and then decays exponentially at rate k (Yu et al., 2018):

a (t) = \sum_{n} d ∙ (e^{- k [t - (n - 1) τ]})

with n = 1,2,… the number of times the treatment dose d is applied.

For ‘constant PK’, the drug concentration is independent of time and simplifies to a = d, whereas for ‘ramp PK’, the drug concentration increases linearly over a time k2cmax (hence the rate of drug concentration increase is given as d/k2cmax) and then stays constant for the rest of the treatment period. The value for k2cmax used for most simulations (48 hr) was taken from literature and describes an example of AMP production timing during a natural immune response (Haine et al., 2008).

Implementation and simulation

Request a detailed protocol

The model was implemented in R using the package adaptivetau (Johnson, 2019) for stochastic implementations via the Gillespie algorithm. We focused on stochastic simulations as we were particularly interested in the timing and probability of the de novo rise and fixation of multiple mutations. We used the package adaptivetau because it allows for time-varying reaction rates, which was necessary in order to incorporate drug-concentration-dependent death rates. It also allows for deterministic simulation of a subset of rates as we did not want the AB concentration to be stochastic. For increased accuracy, we changed the epsilon parameter (which describes the tolerance of relative rate changes in step size selection) to 0.01, which we found gave the same results as exact simulations.

We calculated treatment failure probability as the frequency of runs, in which bacteria were not eradicated at the end of the treatment period (200 hr). Mutational diversity was calculated using the Shannon index, which takes into account the richness and evenness of the distribution of mutant subpopulations, and either averaging the maximum per treatment period (for most results) or the end diversity (Figure 3—figure supplement 3) over all simulation runs. Treatment failure probability and mutational diversity were plotted using the R function filled.contour, which, as far as we could ascertain, interpolates linearly between (potentially irregularly) spaced grid points. To increase the appeal of our figures, we increased the option nlevels from the default value of 20 to 50. The values for treatment failure and diversity shown in the contour plots were then averaged over the whole multi- or single-step area (colored triangles shown in Figure 2) in order to compare different treatment strategies.

The difference in treatment failure and mutational diversity between the two antimicrobial classes (PK profiles) was obtained by subtracting the corresponding values after every simulation of an AMP treatment from the one obtained in a simulation for an AB treatment and plotting the individual resulting differences (for 500 simulations) as well as the density via violin plots.

Model parameters other than benefit, cost, and drug dose are taken from Yu et al., 2018 (Supplementary file 1). The two different antimicrobial drug classes were defined based on previous experimental and theoretical work (Yu et al., 2018; Rodríguez-Rojas et al., 2014; Yu et al., 2016) by using two parameter sets: for the AB class, the mutation rate was 3 * 10⁻⁶ per division, κ was 1.5 and $ψ_{m i n}$ was –5 h⁻¹; whereas for the AMP class, the mutation rate was 10⁻⁶, κ was 5 and $ψ_{m i n}$ was –50 h⁻¹.

The described code has been made available as an R package (Source code 1).

Selection coefficient analysis

Request a detailed protocol

Selection coefficients for our PD model were calculated under the assumption that the sensitive population is very small compared to the carrying capacity, which is a good approximation to the selection pressure at the start of an infection. This means that we can neglect the logistic growth term in our calculations. As the results were very similar to assuming a population at the carrying capacity (which is an approximation for an infection that has had time to establish itself), we will focus on the selection coefficients with a small starting population. Selection coefficients were determined through eigenvalues obtained from the Jacobi matrix given by $\frac{d M}{d t} = [\begin{matrix} \frac{d M_{0}}{d t} \\ ⋮ \\ \frac{d M_{n}}{d t} \end{matrix}] = N (i) * M (t)$ :

N (i) = (\begin{matrix} \begin{matrix} r ∙ (1 - u) - (γ + γ_{0}) & 0 \\ r ∙ u & r ∙ {(1 - c)}^{1} ∙ (1 - u) - (γ + γ_{1}) \end{matrix} & \begin{matrix} \dots \\ \dots \end{matrix} & \begin{matrix} \dots & 0 \\ \dots & ⋮ \end{matrix} \\ \begin{matrix} 0 & r ∙ {(1 - c)}^{1} ∙ u \end{matrix} & ⋱ & ⋮ \\ \begin{matrix} ⋮ & ⋱ \\ 0 & \dots \end{matrix} & ⋱ & \begin{matrix} r ∙ {(1 - c)}^{n - 1} ∙ (1 - u) - (γ + γ_{n - 1}) & 0 \\ r ∙ {(1 - c)}^{n - 1} ∙ u & r ∙ {(1 - c)}^{n} ∙ (1 - u) - (γ + γ_{n}) \end{matrix} \end{matrix})

The eigenvalues of N(i) are its diagonal entries, which correspond to the net growth of each population. We calculated selection coefficients for each of the mutational subpopulations in our model as the difference in growth rates between bacteria with i mutations and bacteria with i–1 mutations (i.e., the difference between their eigenvalues):

s_{i} = g r o w t h (M_{i}) - g r o w t h (M_{i - 1}) = (r ∙ {(1 - c)}^{i} ∙ (1 - u) - (γ + γ_{i})) - (r ∙ {(1 - c)}^{i - 1} ∙ (1 - u) - (γ + γ_{i - 1}))

The AB concentration over time was calculated deterministically using the R package deSolve (Soetaert et al., 2010) in order to calculate the death rates due to antimicrobial treatment.

The difference between the parameter sets for the two antimicrobial classes used here lies in the higher mutation rate u, lower κ, and higher $ψ_{m i n}$ for AB treatments (Yu et al., 2018). Hence, we investigated the importance of the two PD parameters κ and $ψ_{m i n}$ by calculating the selection coefficients using the AMP parameter set and swapping either κ or $ψ_{m i n}$ with that of the AB parameter set.

More generally, we can consider the Jacobi matrix for the resistant populations invading at the mutant-free equilibrium:

N (i) = (\begin{array}{ccc} \begin{array}{cc} r \cdot {(1 - c)}^{1} \cdot (1 - u) \cdot (1 - \frac{M_{0}^{*}}{K}) - (γ + γ_{1}) & 0 \\ r \cdot {(1 - c)}^{1} \cdot u \cdot (1 - \frac{M_{0}^{*}}{K}) & r \cdot {(1 - c)}^{2} \cdot (1 - u) \cdot (1 - \frac{M_{0}^{*}}{K}) - (γ + γ_{2}) \end{array} & \begin{matrix} \dots \\ 0 \end{matrix} & \begin{array}{cc} \dots & 0 \\ \dots & ⋮ \end{array} \\ \begin{array}{cc} 0 & \dots \end{array} & ⋱ & \begin{array}{cc} ⋱ & ⋮ \end{array} \\ \begin{array}{cc} ⋮ & ⋱ \\ 0 & \dots \end{array} & \dots & \begin{array}{cc} r \cdot {(1 - c)}^{n - 1} \cdot (1 - u) \cdot (1 - \frac{M_{0}^{*}}{K}) - (γ + γ_{n - 1}) & 0 \\ r \cdot {(1 - c)}^{n} \cdot u \cdot (1 - \frac{M_{0}^{*}}{K}) & r \cdot {(1 - c)}^{n} \cdot (1 - u) \cdot (1 - \frac{M^{*}}{K}) - (γ + γ_{n}) \end{array} \end{array})

The criteria for invasion of a mutant into the susceptible population is then that the eigenvalue of the mutant has to be bigger than zero, that is,

λ_{n} = r \cdot {(1 - c)}^{n} \cdot (1 - u) \cdot (1 - \frac{M_{0}^{*}}{K}) - (γ + γ_{n}) > 0

Inserting the mutant-free equilibrium ${M_{0}}^{*} = K ∙ (1 - \frac{γ + γ_{0}}{r ∙ (1 - u)})$ yields

γ + γ_{0} > \frac{γ + γ_{n}}{{(1 - c)}^{n}}

which means that bacterial cells with n mutations can invade if the death rate of the sensitive strain is higher than the death rate of the mutant normalized by the cost of the mutation(s).

Horizontal gene transfer

Request a detailed protocol

We added HGT to the model by allowing for an additional resistance gene (with benefit b_p and cost c_p) to be acquired, which gives resistance in a single step. Hence, the benefit b_p and the corresponding cost c_p were adjusted with respect to the drug dose applied by using a benefit that would increase the MIC 20% above the applied drug dose and calculating the cost accordingly through the linear correlation obtained from Table 1. We assume that the bacterial population under investigation has not yet acquired the HGT element, and initial transfer has to come from the environment, that is, initial conditions were the same as for simulations without HGT and $M_{p} (0)$ = 0. This gene can be acquired at a low rate $α$ from the environment or at a density-dependent rate $β$ , which we assumed to be on the same order of magnitude as the mutation rate (Bakkeren et al., 2019).

The horizontally transferred gene can be acquired by sensitive or mutant bacterial populations, and cells containing HGT resistance can still acquire further mutations (but not further HGT resistance). Hence, we assume that HGT resistance is, for example, acquired via a specific resistance gene on a plasmid (typically a plasmid can only be acquired once per cell) and that the resistance gene from this plasmid (e.g., using enzymatic drug inactivation) acts through a different mechanism than resistance by chromosomal mutation (e.g., modification of the drug target) (van Hoek et al., 2011).

The equations were modified as follows:

\frac{d M_{i}}{d t} = r ∙ {(1 - c)}^{i - 1} ∙ u ∙ M_{i - 1} ∙ (1 - \frac{M}{K}) + r ∙ {(1 - c)}^{i} ∙ (1 - u) ∙ M_{i} ∙ (1 - \frac{M}{K}) - (γ + γ_{i}) ∙ M_{i} - (α + β ∙ M_{p}) ∙ M_{i}

\begin{array}{ll} \frac{d M_{p i}}{d t} = r \cdot ({(1 - c)}^{i - 1} * c_{p}) \cdot u \cdot M_{p i - 1} \cdot (1 - \frac{M}{K}) + r \cdot ({(1 - c)}^{i} * c_{p}) \cdot (1 - u) \cdot M_{p i} \cdot (1 - \frac{M}{K}) - \\ (γ + γ_{p i}) \cdot M_{p i} + (α + β \cdot M_{p}) \cdot M_{i} \end{array}

M = \sum_{i} M_{i} + \sum_{i} M_{p i}

M_{p} = \sum_{i} M_{p i}

Here, M_pi is the bacterial subpopulation carrying the HGT gene and i mutations, M_p the total number of HGT subpopulations, and M the total number of all bacterial populations. Relative population frequencies were calculated at the end of the treatment period by dividing the cell number of each subpopulation through the whole population size.

Adaptive treatment

Request a detailed protocol

In adaptive treatment, the goal is not to eradicate the bacterial population entirely but to adjust the treatment dose continuously in order to keep the pathogen level below a certain upper limit. Hansen et al., 2017 calculated the threshold of resistant cells that are necessary at the beginning of the treatment for adaptive treatment to outperform aggressive treatment (i.e., giving the full dose right away), which is based on the idea that sensitive cells provide a risk for becoming resistant through mutation and a benefit through growth competition with the resistant cells at the same time. Hansen et al., 2017 only considered one mutation to resistance, which means that their risk subpopulation and competitor subpopulation was the same. If we consider however sequential mutational steps, then the risk population only consists of the subpopulation one mutation away from full resistance (which will be the mth mutation), whereas the competitor population for the fully resistant strain contains all (partially) sensitive bacteria (i.e., including mutant strains, which are not fully resistant to the highest possible treatment dose). Therefore, the threshold of resistant bacteria is given by (compare to [4] in Hansen et al., 2017):

C o m p e t i t i v e b e n e f i t o f (p a r t i a l l y) s e n s i t i v e s = M u t a t i o n a l r i s k f r o m M_{m - 1}

r ∙ {(1 - c)}^{m} ∙ M_{m} ∙ δ ∙ (P_{m a x} - M_{m}) = u ∙ r ∙ {(1 - c)}^{m - 1} ∙ (1 - δ P_{m a x}) ∙ M_{m - 1}

M_{m} = \frac{u \cdot (1 - δ P_{m a x}) \cdot M_{m - 1}}{(1 - c) \cdot δ \cdot \sum_{j = 0}^{m - 1} M_{j}} w i t h (P_{m a x} - M_{m}) = \sum_{j = 0}^{m - 1} M_{j}

Here, $δ$ describes the strength of competition and $P_{m a x}$ the upper limit of acceptable pathogen burden. Note that here it is assumed that all pathogens (regardless of drug sensitivity) contribute equally to competition (Hansen et al., 2017). This leads to a quadratic equation for the subpopulation with m mutations, $M_{m}$ , which we used to calculate how the resistant population threshold for adaptive treatment (i.e., the initial density $M_{m 0}$ above which adaptive treatment is more favorable) differs between multi- and single-step resistance patterns (Figure 4).

We implemented adaptive treatment in our model by setting a defined upper bound of acceptable pathogen cells, which was equal to the starting density in these simulations (Hansen et al., 2017) (i.e., assuming that the bacterial infection already progressed to a level at which treatment becomes necessary). We used a relatively low acceptable burden of 10⁵ CFU, which is supported by bacterial loads in, for example, urinary tract infections (Schmiemann et al., 2010). Note that defining an acceptable limit of pathogen burden in clinical settings is far more intricate as a patient’s individual biology will play a significant role and is beyond the scope of this paper. We adjusted the treatment dose in order to keep the pathogen load at or below this threshold value but the subpopulations of at least partially sensitive cells as big as possible (Figure 4—figure supplement 1): specifically, we increased the treatment dose to the MIC of the highest resistant subpopulation when its frequency exceeded 1% of the total population and the total pathogen load was higher than our set acceptable burden – until the maximum dose set for a specific treatment simulation was reached; after which the maximum dose was applied continuously. The (partially) sensitive cells serve as competitors for the resistant strain that carries a mutational growth cost and can be outcompeted at low drug doses (Hansen et al., 2017). At the same time, subpopulations that are one step away from the resistant population provide a risk population as they are likely to gain resistance.

For simulations of adaptive and aggressive treatment, we started from a population with neutral heterogeneity, meaning that we calculated the steady-state number of cells with a specific number of mutations given a certain cost (and benefit) in the absence of drug selection. As we want to compare the time difference to treatment failure between adaptive and aggressive treatment for single- and multi-step patterns, we initially add to this ‘neutral population’ the predicted number of resistant cells necessary to make adaptive treatment superior to aggressive treatment. The drug dose in adaptive treatments was then adjusted to keep the number of pathogens below the acceptable burden as described above. The time of treatment failure was determined as the time where the total pathogen population crossed 10⁸ CFUs. We compared adaptive and aggressive treatment by dividing the time to treatment failure obtained from the adaptive strategy by the one obtained with the aggressive strategy, yielding the fold difference in treatment success duration.

Data availability

All data and code generated or analysed during this study are included in the manuscript and supporting files. Source code has been provided for Figures 2-4, as well as Figure 1—figure supplement 2, Figure 2—figure supplements 1–6, Figure 3—figure supplements 1–7, and Figure 4—figure supplements 1–2 in the form of an R package. Source data has been provided for Table 1, Figure 1B and Figure 1—figure supplement 1.

The following previously published data sets were used

1. Melnyk A
2. Wong A
3. Kassen R
(2015) Dryad Digital Repository
The fitness costs of antibiotic resistance mutations.
https://doi.org/10.5061/dryad.5rc47

References

(2017) Using chemical reaction kinetics to predict optimal antibiotic treatment strategies
PLOS Computational Biology 13:e1005321.
https://doi.org/10.1371/journal.pcbi.1005321
- PubMed
- Google Scholar
(2004) In vitro activities of mutant prevention concentration-targeted concentrations of fluoroquinolones against Staphylococcus aureus in a pharmacodynamic model
International Journal of Antimicrobial Agents 24:150–160.
https://doi.org/10.1016/j.ijantimicag.2004.03.011
- PubMed
- Google Scholar
1. Andersson DI
2. Balaban NQ
3. Baquero F
4. Courvalin P
5. Glaser P
6. Gophna U
7. Kishony R
8. Molin S
9. Tønjum T
(2020) Antibiotic resistance: turning evolutionary principles into clinical reality
FEMS Microbiology Reviews 44:171–188.
https://doi.org/10.1093/femsre/fuaa001
- PubMed
- Google Scholar
1. Andersson DI
2. Hughes D
(2010) Antibiotic resistance and its cost: is it possible to reverse resistance?
Nature Reviews Microbiology 8:260–271.
https://doi.org/10.1038/nrmicro2319
- PubMed
- Google Scholar
1. Andersson DI
2. Hughes D
(2014) Microbiological effects of sublethal levels of antibiotics
Nature Reviews Microbiology 12:465–478.
https://doi.org/10.1038/nrmicro3270
- PubMed
- Google Scholar
1. Bakkeren E
2. Huisman JS
3. Fattinger SA
4. Hausmann A
5. Furter M
6. Egli A
7. Slack E
8. Sellin ME
9. Bonhoeffer S
10. Regoes RR
11. Diard M
12. Hardt WD
(2019) Salmonella persisters promote the spread of antibiotic resistance plasmids in the gut
Nature 573:276–280.
https://doi.org/10.1038/s41586-019-1521-8
- PubMed
- Google Scholar
1. Balmer O
2. Tanner M
(2011) Prevalence and implications of multiple-strain infections
The Lancet Infectious Diseases 11:868–878.
https://doi.org/10.1016/S1473-3099(11)70241-9
- PubMed
- Google Scholar
(1997) Evaluating treatment protocols to prevent antibiotic resistance
PNAS 94:12106–12111.
https://doi.org/10.1073/pnas.94.22.12106
- PubMed
- Google Scholar
(2018) Ecology determines how low antibiotic concentration impacts community composition and horizontal transfer of resistance genes
Communications Biology 1:1–8.
https://doi.org/10.1038/s42003-018-0041-7
- PubMed
- Google Scholar
(2020) The within-host evolution of antimicrobial resistance in Mycobacterium tuberculosis
FEMS Microbiology Reviews 1:fuaa0711.
https://doi.org/10.1093/femsre/fuaa071
- Google Scholar
1. Chakrabarti S
2. Michor F
(2017) Pharmacokinetics and drug interactions determine optimum combination strategies in computational models of Cancer evolution
Cancer Research 77:3908–3921.
https://doi.org/10.1158/0008-5472.CAN-16-2871
- PubMed
- Google Scholar
1. Chevereau G
2. Dravecká M
3. Batur T
4. Guvenek A
5. Ayhan DH
6. Toprak E
7. Bollenbach T
(2015) Quantifying the determinants of evolutionary dynamics leading to drug resistance
PLOS Biology 13:e1002299-18.
https://doi.org/10.1371/journal.pbio.1002299
- Google Scholar
(2020) Multi-scale modeling of drug binding kinetics to predict drug efficacy
Cellular and Molecular Life Sciences 77:381–394.
https://doi.org/10.1007/s00018-019-03376-y
- PubMed
- Google Scholar
(2019) Bacteria from natural populations transfer plasmids mostly towards their kin
PNAS 286:20191110.
https://doi.org/10.1098/rspb.2019.1110
- Google Scholar
1. Dortet L
2. Glaser P
3. Kassis-Chikhani N
4. Girlich D
5. Ichai P
6. Boudon M
7. Samuel D
8. Creton E
9. Imanci D
10. Bonnin R
11. Fortineau N
12. Naas T
(2018) Long-lasting successful dissemination of resistance to oxazolidinones in MDR Staphylococcus epidermidis clinical isolates in a tertiary care hospital in France
Journal of Antimicrobial Chemotherapy 73:41–51.
https://doi.org/10.1093/jac/dkx370
- PubMed
- Google Scholar
1. Drlica K
(2003) The mutant selection window and antimicrobial resistance
Journal of Antimicrobial Chemotherapy 52:11–17.
https://doi.org/10.1093/jac/dkg269
- PubMed
- Google Scholar
1. Drusano GL
(2004) Antimicrobial pharmacodynamics: critical interactions of 'bug and drug'
Nature Reviews Microbiology 2:289–300.
https://doi.org/10.1038/nrmicro862
- PubMed
- Google Scholar
1. Dunai A
2. Spohn R
3. Farkas Z
4. Lázár V
5. Györkei Á
6. Apjok G
7. Boross G
8. Szappanos B
9. Grézal G
10. Faragó A
11. Bodai L
12. Papp B
13. Pál C
(2019) Rapid decline of bacterial drug-resistance in an antibiotic-free environment through phenotypic reversion
eLife 8:e47088.
https://doi.org/10.7554/eLife.47088
- PubMed
- Google Scholar
Conference
1. EUCAST
(2020)
European committee on antimicrobial susceptibility testing

Testing Breakpoint Tables for Interpretation of MICs and Zone Diameters.
- Google Scholar
1. Foo J
2. Chmielecki J
3. Pao W
4. Michor F
(2012) Effects of pharmacokinetic processes and varied dosing schedules on the dynamics of acquired resistance to erlotinib in EGFR-Mutant lung Cancer
Journal of Thoracic Oncology 7:1583–1593.
https://doi.org/10.1097/JTO.0b013e31826146ee
- Google Scholar
1. Foo J
2. Michor F
(2009) Evolution of resistance to targeted anti-cancer therapies during continuous and pulsed administration strategies
PLOS Computational Biology 5:e1000557.
https://doi.org/10.1371/journal.pcbi.1000557
- PubMed
- Google Scholar
(2009) Adaptive therapy
Cancer Research 69:4894–4903.
https://doi.org/10.1158/0008-5472.CAN-08-3658
- PubMed
- Google Scholar
1. Goldstein BP
(2014) Resistance to rifampicin: a review
The Journal of Antibiotics 67:625–630.
https://doi.org/10.1038/ja.2014.107
- Google Scholar
1. Haine ER
2. Pollitt LC
3. Moret Y
4. Siva-Jothy MT
5. Rolff J
(2008) Temporal patterns in immune responses to a range of microbial insults (Tenebrio molitor)
Journal of Insect Physiology 54:1090–1097.
https://doi.org/10.1016/j.jinsphys.2008.04.013
- PubMed
- Google Scholar
(2017) How to use a chemotherapeutic agent when resistance to it threatens the patient
PLOS Biology 15:e2001110.
https://doi.org/10.1371/journal.pbio.2001110
- PubMed
- Google Scholar
1. Hansen E
2. Karslake J
3. Woods RJ
4. Read AF
5. Wood KB
(2020) Antibiotics can be used to contain drug-resistant Bacteria by maintaining sufficiently large sensitive populations
PLOS Biology 18:e3000713.
https://doi.org/10.1371/journal.pbio.3000713
- PubMed
- Google Scholar
1. Hughes D
2. Andersson DI
(2017) Evolutionary trajectories to antibiotic resistance
Annual Review of Microbiology 71:579–596.
https://doi.org/10.1146/annurev-micro-090816-093813
- PubMed
- Google Scholar
1. Jin DJ
2. Gross CA
(1988) Mapping and sequencing of mutations in the Escherichia coli rpoB gene that lead to rifampicin resistance
Journal of Molecular Biology 202:45–58.
https://doi.org/10.1016/0022-2836(88)90517-7
- PubMed
- Google Scholar
1. Jochumsen N
2. Marvig RL
3. Damkiær S
4. Jensen RL
5. Paulander W
6. Molin S
7. Jelsbak L
8. Folkesson A
(2016) The evolution of antimicrobial peptide resistance in Pseudomonas aeruginosa is shaped by strong epistatic interactions
Nature Communications 7:13002.
https://doi.org/10.1038/ncomms13002
- PubMed
- Google Scholar
Software
1. Johnson P
(2019)
Tau-Leaping Stochastic Simulation, version R package version 2.2-3

Tau-Leaping Stochastic Simulation.
(2014) Inducible Defenses Stay Up Late: Temporal Patterns of Immune Gene Expression in Tenebrio molitor
G3: Genes, Genomes, Genetics 4:947–955.
https://doi.org/10.1534/g3.113.008516
- Google Scholar
1. Joo H-S
2. Fu C
3. Otto M
(2016) Bacterial strategies of resistance to antimicrobial peptides
Philosophical Transactions of the Royal Society B: Biological Sciences 371:20150292.
https://doi.org/10.1098/rstb.2015.0292
- Google Scholar
1. Kintses B
2. Méhi O
3. Ari E
4. Számel M
5. Györkei Á
6. Jangir PK
7. Nagy I
8. Pál F
9. Fekete G
10. Tengölics R
11. Nyerges Á
12. Likó I
13. Bálint A
14. Molnár T
15. Bálint B
16. Vásárhelyi BM
17. Bustamante M
18. Papp B
19. Pál C
(2019) Phylogenetic barriers to horizontal transfer of antimicrobial peptide resistance genes in the human gut Microbiota
Nature Microbiology 4:447–458.
https://doi.org/10.1038/s41564-018-0313-5
- PubMed
- Google Scholar
(2017) Antimicrobial peptide exposure selects for Staphylococcus aureus resistance to human defence peptides
Journal of Antimicrobial Chemotherapy 72:115–127.
https://doi.org/10.1093/jac/dkw381
- PubMed
- Google Scholar
(2005) Bacterial persistence: a model of survival in changing environments
Genetics 169:1807–1814.
https://doi.org/10.1534/genetics.104.035352
- PubMed
- Google Scholar
(2020) Antimicrobial peptides: application informed by evolution
Science 368:eaau5480.
https://doi.org/10.1126/science.aau5480
- PubMed
- Google Scholar
(1999) Plasmid transfer in the animal intestine and other dynamic bacterial populations: the role of community structure and environment
Microbiology 145:2615–2622.
https://doi.org/10.1099/00221287-145-9-2615
- PubMed
- Google Scholar
1. Liu C
2. Bayer A
3. Cosgrove SE
4. Daum RS
5. Fridkin SK
6. Gorwitz RJ
7. Kaplan SL
8. Karchmer AW
9. Levine DP
10. Murray BE
11. J Rybak M
12. Talan DA
13. Chambers HF
(2011) Clinical practice guidelines by the infectious diseases society of america for the treatment of methicillin-resistant Staphylococcus aureus infections in adults and children: executive summary
Clinical Infectious Diseases 52:285–292.
https://doi.org/10.1093/cid/cir034
- PubMed
- Google Scholar
(2013) Mechanisms and fitness costs of resistance to antimicrobial peptides LL-37, CNY100HL and wheat germ histones
PLOS ONE 8:e68875.
https://doi.org/10.1371/journal.pone.0068875
- PubMed
- Google Scholar
(2018) Genomics of experimental adaptation of Staphylococcus aureus to a natural combination of insect antimicrobial peptides
Scientific Reports 8:15359.
https://doi.org/10.1038/s41598-018-33593-7
- PubMed
- Google Scholar
(2009) Interplay in the selection of fluoroquinolone resistance and bacterial fitness
PLOS Pathogens 5:e1000541.
https://doi.org/10.1371/journal.ppat.1000541
- PubMed
- Google Scholar
(2012) Dosing regimen matters: the importance of early intervention and rapid attainment of the pharmacokinetic/pharmacodynamic target
Antimicrobial Agents and Chemotherapy 56:2795–2805.
https://doi.org/10.1128/AAC.05360-11
- PubMed
- Google Scholar
(2015) The fitness costs of antibiotic resistance mutations
Evolutionary Applications 8:273–283.
https://doi.org/10.1111/eva.12196
- PubMed
- Google Scholar
(2020) Antimicrobial host defence peptides: functions and clinical potential
Nature Reviews Drug Discovery 19:311–332.
https://doi.org/10.1038/s41573-019-0058-8
- PubMed
- Google Scholar
1. Nielsen EI
2. Friberg LE
(2013) Pharmacokinetic-pharmacodynamic modeling of antibacterial drugs
Pharmacological Reviews 65:1053–1090.
https://doi.org/10.1124/pr.111.005769
- PubMed
- Google Scholar
1. Oliver A
2. Cantón R
3. Campo P
4. Baquero F
5. Blázquez J
(2000) High frequency of hypermutable Pseudomonas aeruginosa in cystic fibrosis lung infection
Science 288:1251–1253.
https://doi.org/10.1126/science.288.5469.1251
- PubMed
- Google Scholar
1. Olofsson SK
2. Cars O
(2007) Optimizing drug exposure to minimize selection of antibiotic resistance
Clinical Infectious Diseases 45:S129–S136.
https://doi.org/10.1086/519256
- PubMed
- Google Scholar
1. Read AF
2. Day T
3. Huijben S
(2011) The evolution of drug resistance and the curious orthodoxy of aggressive chemotherapy
PNAS 108:10871–10877.
https://doi.org/10.1073/pnas.1100299108
- PubMed
- Google Scholar
1. Regoes RR
2. Wiuff C
3. Zappala RM
4. Garner KN
5. Baquero F
6. Levin BR
(2004) Pharmacodynamic functions: a multiparameter approach to the design of antibiotic treatment regimens
Antimicrobial Agents and Chemotherapy 48:3670–3676.
https://doi.org/10.1128/AAC.48.10.3670-3676.2004
- PubMed
- Google Scholar
1. Roach DR
2. Leung CY
3. Henry M
4. Morello E
5. Singh D
6. Di Santo JP
7. Weitz JS
8. Debarbieux L
(2017) Synergy between the host immune system and bacteriophage is essential for successful phage therapy against an acute respiratory pathogen
Cell Host & Microbe 22:38–47.
https://doi.org/10.1016/j.chom.2017.06.018
- PubMed
- Google Scholar
(2014) Antimicrobials, stress and mutagenesis
PLOS Pathogens 10:e1004445.
https://doi.org/10.1371/journal.ppat.1004445
- PubMed
- Google Scholar
(2019) Evolutionary pathways to antibiotic resistance are dependent upon environmental structure and bacterial lifestyle
eLife 8:e47612.
https://doi.org/10.7554/eLife.47612
- PubMed
- Google Scholar
(2010) The diagnosis of urinary tract infection
Deutsches Aerzteblatt Online 107:361–367.
https://doi.org/10.3238/arztebl.2010.0361
- Google Scholar
(2010) Solving differential equations in R: package deSolve
Journal of Statistical Software 33:1–25.
https://doi.org/10.18637/jss.v033.i09
- Google Scholar
1. Spohn R
2. Daruka L
3. Lázár V
4. Martins A
5. Vidovics F
6. Grézal G
7. Méhi O
8. Kintses B
9. Számel M
10. Jangir PK
11. Csörgő B
12. Györkei Á
13. Bódi Z
14. Faragó A
15. Bodai L
16. Földesi I
17. Kata D
18. Maróti G
19. Pap B
20. Wirth R
21. Papp B
22. Pál C
(2019) Integrated evolutionary analysis reveals antimicrobial peptides with limited resistance
Nature Communications 10:4538.
https://doi.org/10.1038/s41467-019-12364-6
- PubMed
- Google Scholar
1. Stecher B
2. Denzler R
3. Maier L
4. Bernet F
5. Sanders MJ
6. Pickard DJ
7. Barthel M
8. Westendorf AM
9. Krogfelt KA
10. Walker AW
11. Ackermann M
12. Dobrindt U
13. Thomson NR
14. Hardt WD
(2012) Gut inflammation can boost horizontal gene transfer between pathogenic and commensal Enterobacteriaceae
PNAS 109:1269–1274.
https://doi.org/10.1073/pnas.1113246109
- PubMed
- Google Scholar
1. Tsuji BT
2. Pogue JM
3. Zavascki AP
4. Paul M
5. Daikos GL
6. Forrest A
7. Giacobbe DR
8. Viscoli C
9. Giamarellou H
10. Karaiskos I
11. Kaye D
12. Mouton JW
13. Tam VH
14. Thamlikitkul V
15. Wunderink RG
16. Li J
17. Nation RL
18. Kaye KS
(2019) International consensus guidelines for the optimal use of the polymyxins: endorsed by the american college of clinical pharmacy (ACCP), European society of clinical microbiology and infectious diseases (ESCMID), Infectious diseases society of america (IDSA), International society for Anti-infective pharmacology (ISAP), Society of critical care medicine (SCCM), and society of infectious diseases pharmacists (SIDP)
Pharmacotherapy: The Journal of Human Pharmacology and Drug Therapy 39:10–39.
https://doi.org/10.1002/phar.2209
- PubMed
- Google Scholar
(2018) Fitness variation in Isogenic populations leads to a novel evolutionary mechanism for crossing fitness valleys
Communications Biology 1:151.
https://doi.org/10.1038/s42003-018-0160-1
- PubMed
- Google Scholar
1. van Hoek AH
2. Mevius D
3. Guerra B
4. Mullany P
5. Roberts AP
6. Aarts HJ
(2011) Acquired antibiotic resistance genes: an overview
Frontiers in Microbiology 2:1–27.
https://doi.org/10.3389/fmicb.2011.00203
- PubMed
- Google Scholar
(2006) Darwinian evolution can follow only very few mutational paths to fitter proteins
Science 312:111–114.
https://doi.org/10.1126/science.1123539
- PubMed
- Google Scholar
(2002) Biological cost of rifampin resistance from the perspective of Staphylococcus aureus
Antimicrobial Agents and Chemotherapy 46:3381–3385.
https://doi.org/10.1128/AAC.46.11.3381-3385.2002
- PubMed
- Google Scholar
(2018) Evolution of high-level resistance during low-level antibiotic exposure
Nature Communications 9:1599.
https://doi.org/10.1038/s41467-018-04059-1
- PubMed
- Google Scholar
1. Woodford N
2. Ellington MJ
(2007) The emergence of antibiotic resistance by mutation
Clinical Microbiology and Infection 13:5–18.
https://doi.org/10.1111/j.1469-0691.2006.01492.x
- PubMed
- Google Scholar
1. Yu G
2. Baeder DY
3. Regoes RR
4. Rolff J
(2016) Combination effects of antimicrobial peptides
Antimicrobial Agents and Chemotherapy 60:1717–1724.
https://doi.org/10.1128/AAC.02434-15
- PubMed
- Google Scholar
1. Yu G
2. Baeder DY
3. Regoes RR
4. Rolff J
(2018) Predicting drug resistance evolution: insights from antimicrobial peptides and antibiotics
PNAS 285:20172687.
https://doi.org/10.1098/rspb.2017.2687
- Google Scholar
(2017) Evolution of defence cocktails: antimicrobial peptide combinations reduce mortality and persistent infection
Molecular Ecology 26:5334–5343.
https://doi.org/10.1111/mec.14267
- PubMed
- Google Scholar
1. Zasloff M
(2002) Antimicrobial peptides of multicellular organisms
Nature 415:389–395.
https://doi.org/10.1038/415389a
- PubMed
- Google Scholar

Decision letter

Ben S Cooper

Reviewing Editor; Mahidol University, Thailand
Patricia J Wittkopp

Senior Editor; University of Michigan, United States

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

Previously mathematical models have explored how different antibiotic treatment strategies might affect the emergence of antibiotic resistance and subsequent treatment failure, but these have all assumed that resistance arises in a single mutational step. In practice, two or more mutational steps are often needed before clinically significant resistance arises. This paper address this gap, and explores through a series of computer simulations the potential impact of multiple mutational steps on the likelihood of resistance emergence, and the impact of different treatment strategies on such resistance emerging and reveals unexpected complexity in predicting which treatment strategies will minimize resistance evolution.

Decision letter after peer review:

Thank you for submitting your article "Multi-step vs. single-step resistance evolution under different drugs, pharmacokinetics and treatment regimens" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Patricia Wittkopp as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

We would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). Specifically, we are asking editors to accept without delay manuscripts, like yours, that they judge can stand as eLife papers without additional data, even if they feel that they would make the manuscript stronger. Thus the revisions requested below only address clarity and presentation.

Summary:

The emergence of antibiotic resistance through mutations that arise during treatment is of clinical importance for some bacterial species and certain antibiotics. Previously mathematical models have explored how different antibiotic treatment strategies might affect the emergence of such resistance and subsequent treatment failure. However, these models have all assumed that resistance arises in a single mutational step. In practice, two or more mutational steps are often needed before clinically significant resistance arises. This paper address this gap, and explores through a series of computer simulations the potential impact of multiple mutational steps on the likelihood of resistance emergence, and the impact of different treatment strategies on such resistance emerging. The paper also includes a literature review of the costs and benefits of mutations conferring resistance.

Essential revisions:

There are a number of areas where there is a need to clarify aspects of the work, provide additional information, or add to the discussion. Full reviews are appended and all points should be addressed with the exception that, following the discussion, the reviewers and editor agreed there was not a need to reduce the amount of material in the paper.

Specific points that need addressing (not all of which are mentioned in the appended reviews) are:

1. There is a need to improve the framing of the work. There is a danger that the work gives the impression that emergence of resistance through mutations is the primary driver of clinically important antimicrobial resistance. As reference 22 makes clear, this is not the case, though this route to resistance is important for certain bacteria (such as TB) and certain antimicrobial agents. The introduction should highlight this point, and give examples of organisms and antibiotics for which this work is relevant. It also needs to be clearer that the treatment of HGT is only relevant for the limited range of bug-drug combinations where such emergence of resistance during treatment is important.

2. There is also a danger that readers may get the mistaken impression that emergence of resistance is the only consideration when considering dosing. This is clearly not the case, and the best treatment strategies for minimising risk of treatment failure due to resistance are not necessarily the best strategies for patient outcomes overall. Again, the authors should make this clear in the introduction, and make it clear that any findings in this paper should not be used to directly influence patient treatment decisions, and highlight how they do envisage these findings being used. This is important, as modelling studies are often misinterpreted and taken out of context, and without such explicit caveats there is a danger that they are used to support dosing decisions that might cause real harm to patients.

3. Table 1 needs more detail about the nature of the organisms and the environments they are growing in.

4. Clarification of rationale for choice of initial conditions in simulations is needed

5. Need to explain more clearly why diversity is being looked at and why the measure used to measure it is appropriate.

6. More detailed discussion that makes it clearer how the modelling has added to our understanding is needed.

7. There is a need to highlight that antimicrobials are assumed to affect the death rate of bacteria and not their birth rate (which has implications for the type of antibiotics this modelling is relevant for). Can the discussion consider what changes would be expected if this assumption were to be relaxed?

8. Methods need to be described with greater precision – which methods were used when, what approximations were made, and why methods were chosen.

9. Figure 2 needs more text to help readers to interpret it.

Reviewer #1:

This manuscript considers how the details of resistance evolution impact genetic diversity and the probability of treatment failure. Specifically it uses an in silico analysis to compare resistance evolving from single step mutation versus multi-step mutation. Some of the basic patterns concerning costs and benefits of mutation are motivated by a literature search. This study investigates an important and interesting question and is quite ambitious in its undertaking. In addition to the main topic described above it also considers the additional complications of PK/PD dynamics (by considering three different dosing profiles and two types of antimicrobial), horizontal gene transfer and adaptive therapy. My sense is that the manuscript would be much stronger if it covered fewer topics but covered these topics more comprehensively (for example, a more detailed discussion on diversity would be more informative to the reader). In general it was very hard for me to understand how heavily the results depended on the model structure/assumptions.

1. Table 1 is very nice and contains a lot of interesting information. The cost-benefit relationship of mutations is probably heavily influenced by (i) the environment that the mutations evolved in and (ii) the environment the costs and benefits are measured in. Can you indicate in the table when these environments differ from each other? Also, do any of the strains originate from clinical isolates? I wasn't sure if the data in the Table reflected trends I would expect to see in a patient.

2. Simulations assumed that the initial population was completely susceptible (except for the adaptive therapy analysis). I don't understand why this was done. Why wouldn't you start the simulations with a single sensitive bacteria and then use the dynamical equations to grow the population to the "starting population size". That would give you an initial population with a diversity that reflects the assumed dynamics.

3. Can you clarify the reason you are looking at diversity and why the measure you use is appropriate to investigate the questions you have? I was confused about the decision to average diversity over 200 hours, especially since some populations must go extinct well before 200 hours? For populations that don't go extinct I'd be interested to know how the process of treatment modified the diversity (i.e., what is the composition of the population at the end of 200 hours?) In general I found the discussion around diversity confusing. (I suspect the conclusion on L313-314 is heavily influenced by the way diversity is measured in this study.)

4. I would find the paper much more impactful if the main findings listed in the abstract (L21-25) were discussed in more detail in the discussion. For example, the first result that "resistance evolution is substantially limited if two or more mutations are required": before the modeling analysis, under what situations would you not expect this to be the case? How did the modeling help rule these situations out? Etc.

Reviewer #2:

In their manuscript titled "Multi-step vs. single-step resistance evolution under different drugs, pharmacokinetics and treatment regimens", the authors present a very comprehensive study of antimicrobial resistance evolution, based on stochastic simulations of microbial population growth incorporating pharmacokinetic and pharmacodynamic models of various antimicrobial treatments, as well as experimentally-motivated properties of the mutational pathways to resistance. They compare resistance evolution via a single large-effect mutation to resistance evolution via a sequence of mutations, antibiotics to antimicrobial peptides, and different treatment regimens, including adaptive treatment. They mainly focus on de novo mutations but also include horizontal gene transfer. I found this paper extremely interesting and well-written, incorporating many crucial aspects of resistance evolution in an attractive unified framework, and providing multiple important and nontrivial conclusions on how to optimize antimicrobial treatments. It is also very timely, given the pressing public health issue raised by antimicrobial resistance.

In the model developed here (line 443), it is assumed that the antimicrobial affects the death rate of bacteria and not their birth rate, since its effect is in the death rate γ_i caused by drugs. In other words, the focus is on bactericidal drugs. In practice, some antibiotics are bacteriostatic, and many possess combined bacteriostatic and bactericidal effects. It would be good to point out that this assumption is made, and also to discuss what changes would be expected if bacteriostatic or mixed effects were considered.

The description of the model details in "Implementation and simulation" is rather brief. In particular, it is mentioned that "deterministic analysis" and "stochastic implementations via the Gillespie algorithm" were both carried out. It would be important to specify where each of these approaches was employed. More generally, it would also be helpful to discuss why including stochastic effects is important in the analysis, and for which specific points it is crucial. From a technical point of view, it would also be valuable to explain what approximations are made in the specific tau-leaping R package employed for the Gillespie simulations, and how time-varying reaction rate constants (here, death rates impacted by concentration changes) are dealt with in this Gillespie algorithm.

The presentation of Figure 2 (and Figures S2 to S6 and Figure S13) is a bit unusual, with two triangular regions corresponding to single-step and multiple-step mutational pathways. I think that adding a few extra words about this would help the reader. In particular, it would be good to explain why one transitions from single-step to multiple-step pathways on the diagonal – it makes sense, but it is not entirely straightforward to figure this out when first seeing such plots. In addition, it would be good to clarify how (discrete) simulation results are interpolated to obtain these plots.

Apart from my three main remarks and suggestions above, here are some additional more detailed points.

It would also be good to discuss some assumptions in more detail:

l.276 : "We assumed that resistance through mutation or HGT can be acquired independently of each other and that their effects are multiplicative": it's not obvious to me why this should be the case, as one might imagine that a mutation giving full resistance can also be acquired by HGT, and then the effect of having 1 or 2 copies would depend on the resistance mechanism.

l. 409: "the benefit as the geometric mean of MIC or IC 50 increase relative to the wildtype": Please explain this choice.

l.439: "we don't allow for loss of mutations": Please justify this.

l.462: "Mutation rates were generally kept constant but considering higher mutation rates for mutations with lower benefits and costs did not change our results noticeably (Figure S4)": It would be good to explain why higher mutation rates may be associated to mutations with lower benefits and costs. Is it just that there should be more of these?

l. 550: "Hence, the benefit b_p and the corresponding cost c_p were adjusted according to the (maximum) drug dose applied.": please explain how this was done.

l. 560: Please explain what initial conditions are taken, especially for the number Mp of bact that can transfer the resistance by HGT.

l.597: "setting an upper bound of acceptable pathogen cells": please explain how it should be chosen.

l. 597: "We implemented adaptive treatment in our model by setting an upper bound of acceptable pathogen cells, and adjusting the treatment dose in order to keep the pathogen load at or below this threshold but the subpopulations of at least partially sensitive cells as big as possible (Figure S9)": please explain how this was done.

https://doi.org/10.7554/eLife.64116.sa1

Author response

Essential revisions:

There are a number of areas where there is a need to clarify aspects of the work, provide additional information, or add to the discussion. Full reviews are appended and all points should be addressed with the exception that, following the discussion, the reviewers and editor agreed there was not a need to reduce the amount of material in the paper.

Specific points that need addressing (not all of which are mentioned in the appended reviews) are:

1. There is a need to improve the framing of the work. There is a danger that the work gives the impression that emergence of resistance through mutations is the primary driver of clinically important antimicrobial resistance. As reference 22 makes clear, this is not the case, though this route to resistance is important for certain bacteria (such as TB) and certain antimicrobial agents. The introduction should highlight this point, and give examples of organisms and antibiotics for which this work is relevant. It also needs to be clearer that the treatment of HGT is only relevant for the limited range of bug-drug combinations where such emergence of resistance during treatment is important.

We agree that we did not make it clear enough for which types of infections and treatments our work will be relevant. We included a paragraph discussing clinical examples in the introduction (L55-63).

2. There is also a danger that readers may get the mistaken impression that emergence of resistance is the only consideration when considering dosing. This is clearly not the case, and the best treatment strategies for minimising risk of treatment failure due to resistance are not necessarily the best strategies for patient outcomes overall. Again, the authors should make this clear in the introduction, and make it clear that any findings in this paper should not be used to directly influence patient treatment decisions, and highlight how they do envisage these findings being used. This is important, as modelling studies are often misinterpreted and taken out of context, and without such explicit caveats there is a danger that they are used to support dosing decisions that might cause real harm to patients.

We thank the reviewers for pointing this out as we did not intend this work to be a treatment guideline, but to provide a more comprehensive view of antimicrobial resistance evolution via mutations. This has been addressed at several points in the introduction and discussion now (e.g. L29-33, L442ff).

3. Table 1 needs more detail about the nature of the organisms and the environments they are growing in.

This information was added to Table 1 by adding the column “Organism and evolution environment”, and extending the description of benefit and fitness measurements – with the one exception of the Melnyk, Wong and Kassen, 2015 paper (which we explain in the Methods section L449-453), as they themselves review and synthesize 24 different experimental studies, performed under various conditions and using different species (where we only list the species used but not all of the conditions).

4. Clarification of rationale for choice of initial conditions in simulations is needed

We explain now in more detail what our initial conditions in each of the simulations, which was added in the main text (L146f) and in more detail in the Methods section (L489-496). Further, we ran simulations with different initial conditions (starting from neutral heterogeneity), which is discussed in L418-426 and shown in Figure 2—figure supplement 6, Figure 3—figure supplement 7.

5. Need to explain more clearly why diversity is being looked at and why the measure used to measure it is appropriate.

We agree that this measure is not self-explanatory and needed more discussion. This was done in L158-166 and the Methods section (L588-592). We also compare our measure of diversity to a potentially more intuitive one, by calculating diversity only at treatment endpoints, which reduces diversity in comparison to our measure and show the results in Figure 3—figure supplement 3. We also highlight that we are interested in the maximal diversity over the treatment period because it provides information about the usefulness of increasing drug doses during the treatment, or if higher-level resistance mutations are likely to exist already in the population and will then be selected for (L326-333).

6. More detailed discussion that makes it clearer how the modelling has added to our understanding is needed.

We added more detail to the discussion explaining what the expectations prior to this modelling analysis would be and how the outcomes either contradict or confirm them. For example, while it is intuitive that the necessity for more mutations will limit resistance evolution, it is not a priori clear how strong that limitation is and how that plays out for different drug PKs, like ramp PKs, where acquisition of mutations should be facilitated. These points were added to the discussion in various places, for example L291-294, L363-369, L440-444.

7. There is a need to highlight that antimicrobials are assumed to affect the death rate of bacteria and not their birth rate (which has implications for the type of antibiotics this modelling is relevant for). Can the discussion consider what changes would be expected if this assumption were to be relaxed?

We thank the reviewers for drawing our attention to this point. We explain now in the main part and the Methods section more clearly that this assumption was made. But to get more insight into the implications of this assumption, we also ran simulations of bacteriostatic drug action and included its discussion in the Results section (L223-232), the discussion (L342ff) and the Methods section (L516-529). Even though we would have expected a slow-down of resistance evolution overall, due to the coupling of mutation emergence to replication (which is slowed down by bacteriostatic drugs), we found that this was not true for peak PKs, which showed an even higher probability of treatment failure.

8. Methods need to be described with greater precision – which methods were used when, what approximations were made, and why methods were chosen.

This was added to some extent in the main text but in more detail in every section of the Methods section describing a set of simulation runs.

9. Figure 2 needs more text to help readers to interpret it.

We added more details on the presentation and interpretation of our results to the legend of Figure 2:

“The diagonal line shows where the benefit per mutation is exactly equal to the given drug dose and separates single-step (ss, lower orange triangle), where one mutation gives a benefit higher than the applied dose, from multi-step (ms, upper blue triangle) resistance, where more than one mutation is needed for the accumulated benefit to match the drug dose.”

Reviewer #1:

This manuscript considers how the details of resistance evolution impact genetic diversity and the probability of treatment failure. Specifically it uses an in silico analysis to compare resistance evolving from single step mutation versus multi-step mutation. Some of the basic patterns concerning costs and benefits of mutation are motivated by a literature search. This study investigates an important and interesting question and is quite ambitious in its undertaking. In addition to the main topic described above it also considers the additional complications of PK/PD dynamics (by considering three different dosing profiles and two types of antimicrobial), horizontal gene transfer and adaptive therapy. My sense is that the manuscript would be much stronger if it covered fewer topics but covered these topics more comprehensively (for example, a more detailed discussion on diversity would be more informative to the reader). In general it was very hard for me to understand how heavily the results depended on the model structure/assumptions.

1. Table 1 is very nice and contains a lot of interesting information. The cost-benefit relationship of mutations is probably heavily influenced by (i) the environment that the mutations evolved in and (ii) the environment the costs and benefits are measured in. Can you indicate in the table when these environments differ from each other? Also, do any of the strains originate from clinical isolates? I wasn't sure if the data in the Table reflected trends I would expect to see in a patient.

We added the column “Organism and evolution environment” to Table 1 and gave more detailed descriptions in the columns for benefit and fitness measurements. However, for Melnyk, Wong and Kassen, 2015, we only state overall the various species used in the experiments (all of which were pathogenic), as they themselves review and synthesize 24 different experimental studies, which were performed under various conditions (we also explain this in the Methods section L449-453). (See also reply to Editor)

2. Simulations assumed that the initial population was completely susceptible (except for the adaptive therapy analysis). I don't understand why this was done. Why wouldn't you start the simulations with a single sensitive bacteria and then use the dynamical equations to grow the population to the "starting population size". That would give you an initial population with a diversity that reflects the assumed dynamics.

The reviewer brings up a very good point. While we were generally interested in de novo acquisition of resistance mutations, we ran further simulations starting from a neutral population structure, which gave very similar results but an (expectedly) higher treatment failure probability overall, and added this discussion in L418-426 and Figure 2—figure supplement 6, Figure 3—figure supplement 7. (See also reply to Editor)

3. Can you clarify the reason you are looking at diversity and why the measure you use is appropriate to investigate the questions you have? I was confused about the decision to average diversity over 200 hours, especially since some populations must go extinct well before 200 hours? For populations that don't go extinct I'd be interested to know how the process of treatment modified the diversity (i.e., what is the composition of the population at the end of 200 hours?) In general I found the discussion around diversity confusing. (I suspect the conclusion on L313-314 is heavily influenced by the way diversity is measured in this study.)

We agree with the reviewer that our measure of diversity is not immediately intuitive and needs more explanation and discussion. We added this in the main text in L158-166, L326-333 and the Methods section (L588-592). Further, we compare our measure with diversity obtained at the end of treatment, which indeed is much lower than the diversity that can be attained during the treatment, and added this Figure to the SI (Figure 3—figure supplement 3). (See also reply to Editor).

4. I would find the paper much more impactful if the main findings listed in the abstract (L21-25) were discussed in more detail in the discussion. For example, the first result that "resistance evolution is substantially limited if two or more mutations are required": before the modeling analysis, under what situations would you not expect this to be the case? How did the modeling help rule these situations out? Etc.

We discuss now in more detail when we would expect multi-step resistance to deter resistance evolution or when we would expect it to be facilitated and what our model shows in that regard (e.g. L291-294, L363-369, L440-444). (See also reply to Editor).

Reviewer #2:

[…] In the model developed here (line 443), it is assumed that the antimicrobial affects the death rate of bacteria and not their birth rate, since its effect is in the death rate γ_i caused by drugs. In other words, the focus is on bactericidal drugs. In practice, some antibiotics are bacteriostatic, and many possess combined bacteriostatic and bactericidal effects. It would be good to point out that this assumption is made, and also to discuss what changes would be expected if bacteriostatic or mixed effects were considered.

We thank the reviewer for bringing our attention to this important distinction. We now clarify that we are looking at bactericidal drug effects in the main text (L148) and the Methods section (L542f). Further, we simulated bacteriostatic drug action (together with ‘immune system’ removal) and included a discussion on the comparison between bactericidal and bacteriostatic effects in the Results section (L223-232), the discussion (L342ff) and the Methods section (L516-529). (See also reply to Editor)

The description of the model details in "Implementation and simulation" is rather brief. In particular, it is mentioned that "deterministic analysis" and "stochastic implementations via the Gillespie algorithm" were both carried out. It would be important to specify where each of these approaches was employed. More generally, it would also be helpful to discuss why including stochastic effects is important in the analysis, and for which specific points it is crucial. From a technical point of view, it would also be valuable to explain what approximations are made in the specific tau-leaping R package employed for the Gillespie simulations, and how time-varying reaction rate constants (here, death rates impacted by concentration changes) are dealt with in this Gillespie algorithm.

We agree that this was written in a confusing and insufficient manner. In fact, deterministic simulations were only used for selection coefficient analysis, all other simulations were done stochastically. We explain the rationale behind using stochastic simulations (‘to investigate the timing and probability of the de novo resistance’) and give more details on the package used in the Methods section (L578-586).

The presentation of Figure 2 (and Figures S2 to S6 and Figure S13) is a bit unusual, with two triangular regions corresponding to single-step and multiple-step mutational pathways. I think that adding a few extra words about this would help the reader. In particular, it would be good to explain why one transitions from single-step to multiple-step pathways on the diagonal – it makes sense, but it is not entirely straightforward to figure this out when first seeing such plots. In addition, it would be good to clarify how (discrete) simulation results are interpolated to obtain these plots.

We thank the reviewer for the suggestions and we included more details on the presentation as well as the intuition behind the plots in the figure legend:

Further, we added a sentence explaining that the R function filled.contour seems to interpolate discrete grid points linearly to obtain contour plots as shown (L592-595).

Apart from my three main remarks and suggestions above, here are some additional more detailed points.

It would also be good to discuss some assumptions in more detail:

l.276 : "We assumed that resistance through mutation or HGT can be acquired independently of each other and that their effects are multiplicative": it's not obvious to me why this should be the case, as one might imagine that a mutation giving full resistance can also be acquired by HGT, and then the effect of having 1 or 2 copies would depend on the resistance mechanism.

We explain now in the main text (L272ff) as well as in the Methods section (L664-668) that this was done because we are assuming that HGT resistance arises through a specific gene, while mutational resistance occurs due to changes in drug targets. Hence, we assume that these pathways are independent of each other and that HGT resistance is already ‘maximal’ and cannot be acquired twice.

l. 409: "the benefit as the geometric mean of MIC or IC 50 increase relative to the wildtype": Please explain this choice.

We added that this was done because MIC and IC50 are both measures with logarithmic scales (L455f).

l.439: "we don't allow for loss of mutations": Please justify this.

We added a paragraph explaining that we chose to specifically model sequential mutations in order to mimic mutational resistance networks as found in Jochumsen et al., 2016, but that we ignore other effects like loss of resistance mutations, which was found to be rare and is likely playing a minor role in evolutionary trajectories in such networks (L492ff, L532-538).

l.462: "Mutation rates were generally kept constant but considering higher mutation rates for mutations with lower benefits and costs did not change our results noticeably (Figure S4)": It would be good to explain why higher mutation rates may be associated to mutations with lower benefits and costs. Is it just that there should be more of these?

Yes, we would assume that mutations that provide only a small benefit would be more numerous and therefore provide a larger mutational target. We clarified this in L556-560.

l. 550: "Hence, the benefit b_p and the corresponding cost c_p were adjusted according to the (maximum) drug dose applied.": please explain how this was done.

We added that the benefit was chosen to provide an MIC 20% above the given drug dose and that cost was calculated via the correlation with benefit (L653-656).

l. 560: Please explain what initial conditions are taken, especially for the number Mp of bact that can transfer the resistance by HGT.

We clarified that we are starting from a population that has to acquire HGT resistance first from the environment, i.e. Mp=0, in the main text (L268ff) and in the Methods (L656-661).

l.597: "setting an upper bound of acceptable pathogen cells": please explain how it should be chosen.

We added the threshold value that we were using in the simulations and an explanation that this problem is much more intricate in clinical settings (L712-718).

l. 597: "We implemented adaptive treatment in our model by setting an upper bound of acceptable pathogen cells, and adjusting the treatment dose in order to keep the pathogen load at or below this threshold but the subpopulations of at least partially sensitive cells as big as possible (Figure S9)": please explain how this was done.

We added an detailed explanation to the Methods section: ‘we increased the treatment dose to the MIC of the highest resistant subpopulation when its frequency exceeded 1% of the total population and the total pathogen load was higher than our set acceptable burden – until the maximum dose set for a specific treatment simulation was reached; after which the maximum dose was applied continuously’ (L718-724).

https://doi.org/10.7554/eLife.64116.sa2

Article and author information

Author details

Claudia Igler

Institute of Integrative Biology, ETH Zurich, Zurich, Switzerland

Contribution
Conceptualization, Data curation, Software, Formal analysis, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing

For correspondence
claudia.igler@env.ethz.ch

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-7777-546X
Jens Rolff

Evolutionary Biology, Institute for Biology, Freie Universität Berlin, Berlin, Germany

Contribution
Conceptualization, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-1529-5409
Roland Regoes

Institute of Integrative Biology, ETH Zurich, Zurich, Switzerland

Contribution
Conceptualization, Supervision, Funding acquisition, Methodology, Writing - review and editing

For correspondence
roland.regoes@env.ethz.ch

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-8319-5293

Funding

Volkswagen Foundation (96517)

Claudia Igler
Jens Rolff
Roland Regoes

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank D Baeder, S Bonhoeffer, S Lehtinen, and H Alexander for useful discussions and comments on the manuscript. This work was supported by a grant from the Volkswagen Foundation (grant no. 96517).

Senior Editor

Patricia J Wittkopp, University of Michigan, United States

Reviewing Editor

Ben S Cooper, Mahidol University, Thailand

Publication history

Received: October 18, 2020
Accepted: May 4, 2021
Accepted Manuscript published: May 18, 2021 (version 1)
Version of Record published: June 7, 2021 (version 2)

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.