eLife
eLife Sciences Publications, Ltd
Quantifying antibiotic impact on within-patient dynamics of extended-spectrum beta-lactamase resistance
Volume: 9
DOI 10.7554/eLife.49206
•
•
•
• Altmetric

### Notes

Abstract

Bacteria that are resistant to antibiotics are a growing global health crisis. One type of antibiotic resistance arises when certain bacteria that can produce enzymes called extended-spectrum beta-lactamases (or ESBLs for short) become more common in the gut. These enzymes stop important antibiotics, like penicillin, from working. However, exactly which antibiotics and treatment durations contribute to the emergence of this antibiotic resistance remain unknown.

Now, Niehus et al. find certain antibiotics that are associated with an increase in the number of gut bacteria carrying antibiotic resistance genes for ESBL enzymes. First, rectal swabs collected from 133 patients from three European hospitals were analysed to measure the total gut bacteria and the number of genes for ESBL enzymes. These samples had been collected at several time points including when the patient was first admitted to hospital, then every two to three days during their stay, and finally when they were discharged.

Combining the analysis of the samples with details of the patients’ charts showed that treatment with two antibiotics: cefuroxime and ceftriaxone, was linked to an increase in ESBL genes in the gut bacteria. Other antibiotics – namely, meropenem, piperacillin-tazobactam and oral ciprofloxacin – were associated with a decrease in the number of bacteria with ESBL genes. Niehus et al. then performed further analysis to see if different treatment regimens affected how long patients were carrying gut bacteria with ESBL genes. This predicted that a longer course of meropenem, 14 days rather than 5 days, would shorten the length of time patients carried ESBL-resistant bacteria in their guts by 70%, although this effect will likely depend on the location of the hospital and the local prevalence of other types of antibiotic resistance.

This analysis reveals new details about how antibiotic treatment can affect ESBL resistance genes. More studies are needed to understand how antibiotics affect other antibiotic resistance genes and how resistant bacteria spread. This will help scientists understand how much specific antibiotic regimens contribute to antibiotic resistance. It may also help scientists develop new antibiotic treatment strategies that reduce antibiotic resistance.

Keywords
Niehus, van Kleef, Mo, Turlej-Rogacka, Lammens, Carmeli, Goossens, Tacconelli, Carevic, Preotescu, Malhotra-Kumar, Cooper, Garrett, and Colijn: Quantifying antibiotic impact on within-patient dynamics of extended-spectrum beta-lactamase resistance

## Introduction

Antibiotic use can increase resistance prevalence in a host population through multiple pathways (Lipsitch and Samore, 2002). It may: (i) affect the duration of resistance carriage and hence transmission potential; (ii) increase bacterial load of resistant organisms and thus increase transmission; or (iii) selectively suppress host microbial flora where resistance is lacking, which may reduce the potential for transmission of sensitive organisms and also render hosts more susceptible to acquiring resistant bacteria. As well as being important for understanding population dynamics, levels of intestinal resistance are also likely to be important from an individual patient perspective. It has been shown, for instance, that the digestive tract is the primary source of enterobacteria causing bloodstream infections in haematological patients, and a high abundance of beta-lactam resistant enterobacteria in the gut flora is predictive of a high risk of a corresponding drug-resistant bloodstream infection (Woerther et al., 2015). Moreover, colonization with extended-spectrum beta-lactamase (ESBL)-producing Enterobacteriaceae amongst patients receiving cephalosporin-based prophylaxis prior to colorectal surgery is associated with a more than two-fold increase in risk of surgical site infection (Dubinsky-Pertzov et al., 2019). Therefore, quantifying within-host selection dynamics should lead to a better understanding of both individual patient-level and population-level risks and benefits of antibiotic use.

Here we focus on Enterobacteriaceae, a bacterial family that is commonly found in the healthy mammalian gut microbiome (Donnenberg, 1979). Some member genus-species—Klebsiella pneumoniae, Escherichia coli, Enterobacter spp.—are important opportunistic human pathogens that can cause urinary tract, bloodstream, and intra-abdominal infections, as well as hospital-acquired respiratory tract infections. A major concern is the global increase in extended-spectrum beta-lactamase-producing organisms in this family (Coque et al., 2008; Tacconelli et al., 2018; Valverde et al., 2004). ESBL genes – of which the most important and globally widespread is the blaCTX-M gene family – confer resistance to clinically important broad-spectrum antimicrobials, such as third generation cephalosporins (Paterson, 2000). These genes commonly reside on large conjugative plasmids (Bonnet, 2004), and are co-carried with other antibiotic resistance determinants, making them a good marker for multi-drug resistance (MDR) in strains of Enterobacteriaceae (Schwaber et al., 2005). Because Enterobacteriaceae have their main biological niche in the gut microbiome (Masci, 2005), these bacteria are exposed to substantial collateral selection from antibiotics used to treat or prevent infections with other organisms (‘bystander selection’ [Tedijanto et al., 2018]). Quantifying the effects of antibiotic therapy on within-host resistance dynamics will help us to better understand the potential for selection of drug-resistant Enterobacteriaceae associated with different patterns of antibiotic usage.

In this work, we analysed sequential rectal swabs (n = 833) from 133 ESBL positive hospitalised patients from three hospitals (Italy, Romania, Serbia) to study the dynamics of antibiotic resistance gene abundance. Both blaCTX-M gene and, as a proxy for total bacterial load, 16S rRNA gene abundance were determined using quantitative polymerase chain reaction (qPCR). Previously, using a subset of these data, Meletiadis et al. demonstrated a statistical association between exposure to ceftriaxone and increases in blaCTX-M normalised by total bacterial load. Here, we addressed some broader questions. We studied the effects of a range of different antibiotics on the abundance of blaCTX-M and of 16S rRNA, and we aimed to fully characterise the within and between host variation of blaCTX-M and 16S rRNA and their within-host dynamics. For this purpose we developed a novel dynamic model, a state-space model that we fit to fine-grained patient-level measurements and antibiotic exposure data. By incorporating hidden-state dynamics our model allowed us to dissect and quantify different types of data variability, such as noise from qPCR measurement or from the DNA extraction process, and to separate this from the within-host processes. In this way we directly estimated ecologically important parameters such as strength of resistance amplification during antibiotic treatment or the rate of decline of blaCTX-M. We then used our model to make counterfactual predictions about how alternative choices of treatment would impact blaCTX-M carriage duration. The development of this data-driven within-host model and its use in exploring the impact of antibiotic treatment on amplification and loss of resistance is an important step in furthering our quantitative mechanistic understanding of how antibiotic use drives changes in the prevalence of resistance in a population.

## Results

### Patient cohort and treatment

The study enrolled a total of 1102 patients who were screened positive for ESBL producing Enterobacteriaceae at admission, and 133 patients (12%) gave consent to be included in the study: 51 (38%) from Romania; 52 (39%) from Serbia; and 30 (23%) from Italy. The median age was 59 years (range of 23–88), and 46% were female. The median length of hospital stay was 15 days (maximum of 53 days). All patients apart from one had two or more rectal swabs taken, with a median of five swabs per patient (range of 1–15). 114 out of 133 (86%) enrolled patients received antibiotics during their stay and 85 of these 114 (75%) received two or more different antibiotics, which were given both in mono- and combination therapy (see Figure 1). A total of 3993 patient days were observed, of which 2686 (67%) were days with antibiotic therapy (mono- or combination therapy). Table 1 summarises important details of the study. Note that the antibiotics that we considered in this study were exclusively antibacterial drugs, and we ignored treatment with anti-tuberculosis drugs (pyrazinamide and isoniazid), which occurred only in two patients.

Figure 1.
The x-axis scale is identical across panels, the length of one week is given for scale in the top-left corner. Timelines are ordered by length. The y-axis scale differs between panels, with the space between vertical grey lines representing a 10-fold change in the absolute blaCTX-M gene abundance (measured in copy numbers). The left-hand side shows patients who received antibiotic treatment (n = 113), and the two right-hand side columns are patients without antibiotic treatment (n = 19). For clarity, we show only the twelve most frequently used antibiotics in distinct colours and other antibiotics in light grey.Time series plots demonstrating the diverse range of dynamical patterns of blaCTX-M resistance gene abundance across the 132 patients with two or more samples.
Table 1.
Summary of the study.
Number of participating hospitals3 (Serbia, Italy, Romania)
Study duration2 years (Jan 2011 - Dec 2012)
Inclusion criteriaInpatients of medical and surgical wards, adults, non-pregnant, ESBL producting Enterobacteriaceae carriers (at admission)
Number of patients followed up133 (including one with a single swab taken)
Intervals between rectal swabstwo to three days
qPCR targetsblaCTX-M (ESBL resistance gene), 16S rRNA (total bacterial load)
Number of different antibiotics used35
This study (registration number NCT01208519) was conducted by the SATURN consortium, supported by the European Commission under the 7th Framework Program.

The different antibiotic classes, ranked by proportion of treatment days, were cephalosporins 25% ), fluoroquinolones (18%), penicillins (9%), nitroimidazole derivatives (metronidazole) (9%), glycopeptides (8%), carbapenems (5%), and others (26%). Two thirds of antibiotic treatment days were from intravenously administered antibiotics and one third from oral administration. Details on individual antibiotics are given in Table 2.

Table 2.
Intravenous (iv), oral (or), and intramuscular (im) route administration is given in percent of treatment days.
Overview of antibiotic treatments showing the ten most used antibiotics in this patient cohort.
AntibioticNumber of treated patients (total n=133)Route
Ceftriaxone64( 98% iv , 2% im )
Ciprofloxacin34( 67% iv , 33% or )
Metronidazole20( 50% iv , 50% or )
Cefuroxime13( 100% iv , 0% or )
Vancomycin13( 86% iv , 14% or )
Meropenem10( 100% iv , 0% or )
Amikacin9( 100% iv , 0% or )
Amoxicillin-clavulanic acid9( 57% iv , 43% or )
Piperacillin-tazobactam7( 100% iv , 0% or )
Imipenem5( 100% iv , 0% or )

### Resistance dynamics

The time-varying blaCTX-M abundance exhibits a diverse range of dynamic patterns, including monotonic increases and decreases, as well as highly variable non-monotonic behaviour (Figure 1). Qualitatively similar fluctuations in blaCTX-M abundance were seen both in the presence and absence of antibiotic treatment. To determine whether this high level of dynamic variation contained a meaningful biological signal, we first studied temporal autocorrelation. If the observed variability is driven by observation uncertainty – for instance through the swab procedure, DNA extraction, or qPCR process – we expect autocorrelation close to zero in the time series. Conversely, if the observed fluctuations reflect true within-host dynamics in carriage levels, we would generally expect to see positive autocorrelation. We found a clear signal of first-order autocorrelation for both the blaCTX-M and the 16S rRNA gene time series, though autocorrelation was substantially stronger for the blaCTX-M data (Figure 1—figure supplement 1a and b). Using a Bayesian state-space model that decomposes the time series data into an observation component (representing noise due to variability in qPCR runs, and in the procedure of swab taking and sample processing) and a process component (due to the within-host dynamics), we estimated that much of the variability in blaCTX-M and 16S rRNA outcomes was due to measurement error associated with the swab procedure (median estimate of the proportion of total abundance variability attributable to swab error [90% credible interval [CrI]] of 54% [44%, 57%] and 73% [68%, 77%], respectively) (see Figure 1—figure supplement 1c).

However, the blaCTX-M data in particular were found to also contain a strong process component signal, indicating that a median estimate of 36% (90% CrI 30%, 43%) of the variability in the qPCR outcomes was due to underlying within-host dynamics (Figure 1—figure supplement 1c). To further investigate the determinants of blaCTX-M gene variation, we explored how much the blaCTX-M gene load varied between different patients or, over time, within the same patient. Using a Bayesian state-space model (see Methods and Materials) we found 16S rRNA gene abundance to be two orders of magnitude higher than blaCTX-M (median ratio 16S / blaCTX-M [90% CrI] 158 [88, 181]), with an estimated coefficient of variation (ratio of standard deviation to the mean) of 5.5 for 16S rRNA and 32.1 for blaCTX-M. Between-patient abundance of blaCTX-M showed substantially more variability than within-patient abundance (median ratio [90% CrI] 134 [18, 1422]). In contrast, 16S rRNA gene abundance had similar between-patient and within-patient variability (median ratio [90% CrI] 0.8 [0.4,1.7]) (Figure 2). We noted that the rank plots (Figure 1—figure supplement 1a) indicate some convergence problems of ${\sigma }_{bio,16S}$, but several independent runs of the MCMC algorithm with different initial values consistently arrived at the same mean and standard deviation of the posterior estimate.

Figure 2.
In this model, abundance is distributed around individual patient intercepts, which are distributed around a common population intercept. The plot shows individual patient intercepts given as mean posterior estimates (coloured dots) together with posterior predictions for sequence abundance for each patient (grey bars show 80% central quantiles).Variability of 16S abundance and blaCTX-M abundance within and between patient time series using a Bayesian hierarchical model.

### Associating resistance and antibiotic treatment

The change in relative resistance between samples, measured as blaCTX-M abundance divided by 16S rRNA gene abundance, was only slightly elevated in time intervals where antibiotics were given compared to those where they were not (Figure 3a).

Figure 3.
The upper panels show the change in relative resistance between all neighbouring timepoints (black dots), dashed horizontal lines in grey indicate the region of no change. Pairs of violin scatter plots (with the mean values shown as red bars) contrast different treatment that occurred between those timepoints. ‘Yes’ indicates treatment with specified antibiotics and ‘No’ means treatment with other antibiotics or no treatment. The lower three panels show the distribution of mean differences of the change in relative resistance between treatment groups generated through treatment-label permutation (areas in darker grey show 80% central quantiles). The distributions are overlaid with the observed difference (red vertical line). Panel (a) compares treatment with any antibiotic versus no antibiotic (number of intervals without treatment/number of intervals with treatment are 251(N)/449(Y)). Panel (b) compares treatment with antibiotics with activity against Enterobacteriaceae and to which blaCTX-M does not confer resistance (colistin, meropenem, ertapenem, imipenem, amoxicillin-clavulanic acid, ampicillin-sulbactam, piperacillin-tazobactam, gentamicin, amikacin, ciprofloxacin, ofloxacin, levofloxacin, tigecycline, doxicycline) with all other treatment, including no treatment (445(N)/255(Y)). Finally, in panel (c) we consider antibiotics with broad-spectrum activity but to which blaCTX-M does confer resistance (cefepime, ceftazidime, ceftriaxone, cefotaxime, cefuroxime, amoxicillin, ampicillin) (513(N)/187(Y)).Association of antibiotic use with change in relative resistance (abundance of blaCTX-M divided by abundance of 16S rRNA).

However, use of antibiotics with activity against Enterobacteriaceae to which carriage of blaCTX-M does not confer resistance (colistin, meropenem, ertapenem, imipenem, amoxicillin-clavulanic acid, ampicillin-sulbactam, piperacillin-tazobactam, gentamicin, amikacin, ciprofloxacin, ofloxacin, levofloxacin, tigecycline, doxicycline) was associated with a modest decrease in blaCTX-M abundance (Figure 3b). In contrast, the use of antibiotics with broad spectrum killing activity and to which carriage of blaCTX-M does confer resistance (cefepime, ceftazidime, ceftriaxone, cefotaxime, cefuroxime, amoxicillin, ampicillin) was associated with substantially higher increases in relative blaCTX-M abundance (Figure 3c).

### Dynamic antibiotic effect model

Fitting a dynamic model of blaCTX-M abundance and 16S rRNA abundance to the data (133 patients, 833 swabs, 3361 qPCR measurements), we found that cefuroxime and ceftriaxone were associated with increases in both absolute blaCTX-M abundance (mean daily increase [90% CrI] 21% [1%, 42%] and 10% [4%, 17%], respectively) and relative blaCTX-M abundance (14% [-1%, 30%] and 11% [5%, 17%], respectively) (Figure 4). Piperacillin-tazobactam, meropenem and ciprofloxacin (when given orally) were negatively associated with both blaCTX-M (−8% [−18%, 2%], −8% [−17%, 1%], and −8% [−17%, 2%], respectively) and 16S rRNA gene abundance (−3% [−8%, 1%], −3% [−7%, 1%], and −1% [−6%, 3%], respectively), although uncertainty was large (Figure 4). Their effect on relative resistance (blaCTX-M/16S rRNA) also appeared to be negative (−5% [−14%, 5%] for piperacillin-tazobactam, −5% [−14%, 4%] for meropenem, −7% [−15%, 3%] for oral ciprofloxacin). Intravenously administered ciprofloxacin did not show these effects. Imipenem and meropenem had similar effects on blaCTX-M abundance, while no clear effects were evident for amikacin, metronidazole, and amoxicillin-clavulanic acid (Figure 4). The out of sample prediction accuracy using approximated leave-one-out cross-validation (loo-cv) (Watanabe, 2010) (see Materials and methods) for the dynamic model with antibiotic effects is higher than the accuracy of the model without antibiotic effects (Δloo-cv = 4.2 ).

Figure 4.
The bars show estimated daily effects of individual antibiotics on the absolute blaCTX-M abundance (red) and 16S rRNA abundance (light blue) indicating the 80% and 95% highest posterior density intervals (thick and thin horizontal bars, respectively). The model also gives the antibiotic effect on the blaCTX-M/16S relative resistance shown as arrows on the right-hand side. Arrows are in grey for antibiotics with mean effect estimates between −10% and +10%, otherwise they are coloured red (positive selection) and green (negative selection). Route of antibiotic administration (intravenous, iv; oral, or) is indicated in parenthesis.Estimated effects of different antibiotics on within-host dynamics from a multivariable model.

With the dynamic model we are able to make predictions about the time required for the blaCTX-M gene to fall below detection levels. To achieve this, we added to our stochastic model a threshold below which the blaCTX-M gene cannot be detected (see Materials and mthods). Below this threshold the gene may either be lost from the bacterial community, or it may exist in very small reservoirs for example in persister cells (Balaban et al., 2019). The predictions of detectable carriage duration show a high degree of uncertainty, visible as long-tailed predictive distributions (Figure 5). Because of the skew, we report here the median instead of the mean together with 80% credible intervals. We chose the duration of different antimicrobial therapies according to clinical guidelines. Assuming that the estimated antibiotic associations represent causal effects, we find that a single eight day course of cefuroxime or a 14 day course of ceftriaxone substantially prolongs carriage of blaCTX-M, by a median estimate of 147% (80% CrI 13.4%, 577%) for cefuroxime and 120% (80% CrI −8.6%, 492%) for ceftriaxone versus no exposure (Figure 5). Addition of oral ciprofloxacin to a course of amoxicillin-clavulanic acid or ceftriaxone reduces blaCTX-M carriage duration (by approximately 51% [80% CrI −115%, 89%] and 48% [80% CrI −71.1%, 86%]) (Figure 5). A typical 14 day course of meropenem or a 8 day course of piperacillin-tazobactam reduce blaCTX-M carriage duration relative to no treatment (by approximately 42% [80% CrI −25%, 75%] and 41% [80% CrI −45%, 71%], respectively), and each reduces blaCTX-M carriage even more relative to a 7 day course of combined ceftriaxone plus amikacin (by approximately 69% [80% CrI 20%, 89%] and 66% [80% CrI −7%, 88%], respectively) (Figure 5). Finally, a 14 day course of meropenem reduces blaCTX-M resistance carriage relative to a shorter 5 day course (by approximately 69% [80% CrI 20%, 89%]) (Figure 5).

Figure 5.
The distributions on the left-hand side shows model predictions with parameter uncertainty, but assuming deterministic dynamics. The right-hand side shows the predictions with parameter uncertainty as well as Markov process uncertainty. The darker grey areas shows the 50% credible intervals and the white lines show the median predictions. Each density distribution is overlaid with the density line of the no treatment case (dotted line) and its median prediction (dotted vertical line). In the first five rows, we compare predictions for treatment with amoxicillin-clavulanic acid (18 days) and ceftriaxone (14 days), and each in combination treatment with ciprofloxacin. In the next three rows, we compare treatment with ceftriaxone plus amikacin (7 days), meropenem (14 days), and piperacillin-tazobactam (8 days). In the final two rows, we compare a 14 day course of meropenem with a shortened course (5 days). AMC stands for amoxicillin-clavulanic acid.Simulated predictions of blaCTX-M carriage duration under different antibiotic treatments.

## Discussion

By fitting a dynamic model accounting for both observation noise and within-host dynamics to time series data from 133 patients, we quantified the association between antibiotic exposure and changes in rectal swab abundance of gut bacteria and blaCTX-M resistance genes. The largest effects were found for exposures to the second and third generation cephalosporins, cefuroxime and ceftriaxone, both of which were associated with increases in blaCTX-M abundance. Forward simulations indicated that if these associations are causal, exposure to typical courses of these antibiotics would be expected to more than double the carriage duration of blaCTX-M. Both cefuroxime and ceftriaxone have broad-spectrum killing activity (Nahata and Barson, 1985; Neu and Fu, 1978), but have limited activity against ESBL-producing organisms (Livermore and Brown, 2001; Sorlózano et al., 2007). Therefore, a direct selective effect of these two antibiotics is biologically plausible to account for the above finding.

There are are number of advantages to our modelling approach over more classical, associational methods used in related work (Meletiadis et al., 2017). First, because we use a mechanistic model it allows us to directly estimate ecologically important parameters such as strength of resistance amplification under antibiotic selection, which are of inherent interest and can inform further modelling work. Indeed, our predictions of resistance carriage duration (Figure 5) are good examples of the latter. Second, our model uses the variability present in the data to quantify different types of data variability due to the data collection (noise from the qPCR machine, noise from taking the swab and DNA extraction). This allows our model to fully propagate uncertainty to the final estimates. Finally, rather than using only aggregate data (which loses information), our analysis is designed to fully exploit the information available in the time series data. Our work also has a number of important limitations, aside from the obvious risks of confounding present in this observational dataset. Our analysis did not explicitly model changing antibiotic concentrations in the gut, nor did it attempt to explicitly model how antibiotics affect the ecology of the gut bacterial community. While it would have been straightforward to include a pharmacokinetic model of antibiotic concentrations (similar studies have been performed in mice [Jumbe et al., 2003] and pigs [Nguyen et al., 2014]), disentangling direct effects of antibiotic concentrations from indirect effects mediated by other components of the gut flora is far more challenging and beyond the scope of what we considered appropriate with the available data. Instead, our model assumed multiplicative antibiotic effects, which we considered a reasonable simplification of the underlying mechanisms. Multiplicative effects imply that antibiotics alter the daily total bacterial growth rate (16S abundance) and the growth of resistant bacteria relative to the average bacterium (blaCTX-M/16S), but it does not allow for more complex fitness effects due to, for example, synergies between antibiotics (MacLean et al., 2010), or density-dependent effects, whereby antibiotic-mediated killing may depend on bacterial density (Udekwu et al., 2009). Further, many non-antibiotic drugs have been shown to have an impact on human gut bacteria (Maier et al., 2018), but only antibacterials were considered in our analysis. Lastly, all patients in this study were identified (and consenting) ESBL-carriers. Therefore, apart from the potential for selection bias, we assumed that all changes in blaCTX-M abundance were due to within-host dynamics, neglecting the possibility of new acquisitions, which should be the scope for other modelling frameworks that integrate both within- and between-host dynamics.

Antibiotic impact on the human gut microbiome is likely to be an important mediator for the increase of bacterial resistance globally (Donskey, 2004; Relman and Lipsitch, 2018). A large body of theory has been developed that demonstrates the role of within-host processes for understanding population-wide selection of resistance through antibiotic use (Blanquart, 2019; Davies et al., 2019; Webb et al., 2005; Knight et al., 2018; Lipsitch et al., 2000; Lipsitch and Samore, 2002; Webb et al., 2005). However, surprisingly few studies have used data to quantify the effect of antibiotic treatment on resistance abundance within individual gut microbiomes. Two studies involving patients admitted to intensive care units looked at the effect of a preventative antibiotic cocktail (selective digestive decontamination) on gut microbiome resistance in patients, with one study (including n = 13 patients) finding no clear effect (Buelow et al., 2014), and the other (n = 10) showing increases of four different resistance genes associated with treatment (Buelow et al., 2017). Gibson et al. (2016) studied the faecal metagenomics of preterm infants (n = 84) over time and found that treatment with antibiotics was correlated with enrichment of both cognate and noncognate resistance markers. But none of these studies attempted to model resistance dynamics. The modelling framework we have developed enables testable predictions about the impact of different antibiotics and lengths of treatment on the duration of carriage of resistant determinants above detection threshold. Such understanding is an important step toward understanding the spread of antibiotic resistance.

The development of a mechanistic understanding of the relationship between antibiotic use in a population and the proportion of this population in whom resistance can be detected relies on quantifying the antibiotic effects in individual exposed patients, as we do here, but also on quantifying the knock-on effects on transmission to contacts. These indirect effects are likely to be considerable. A recent study in Dutch travellers returning to the Netherlands who had acquired ESBL-producing Enterobacteriaceae carriage overseas found that their new carriage status was associated with a 150% increase in the daily risk of non-carrying household members also becoming ESBL-positive (Arcilla et al., 2017). Developing mechanistic models for the spread of ESBLs and other resistance determinants within host populations accounting for direct and indirect antibiotic effects is an important priority for future research. Such models would help us to understand and predict how changes in antibiotic usage patterns affect the prevalence of antimicrobial resistance in a community and ultimately help to prioritise interventions to reduce the burden of antimicrobial resistance.

## Materials and methods

### Study participants and follow-up

Participants were recruited as part of an observational, prospective, cohort study that included three hospitals (Italy, Serbia and Romania), with known high endemic prevalence of antibiotic resistance in bacterial infections. The hospitals were serving a general urban population. The study was conducted over two years from January 2011 to December 2012 as part of the multi-centre SATURN (‘Impact of Specific Antibiotic Therapies on the prevalence of hUman host ResistaNt bacteria’) project (NO241796; clinicaltrials.gov NTC01208519). The study enrolled adult (>18 y) inpatients of medical and surgical wards, excluding pregnant patients. Enrolled patients were screened at admission for carriage of ESBL-producing Enterobacteriaceae with rectal swabs (E swab, Copan, Italy). Patients who tested positive for ESBL producing Enterobacteriaceae carriage (details below) were included in the follow-up cohort. The target cohort size was calculated to be 400 patients, based on the number of patients that would be required to detect one log difference in resistance abundance with 90% certainty. However, patient recruitment we slower than anticipated. For all follow-up patients (n = 133) rectal swabs were taken every two to three days (as per study protocol) during hospitalisation, which includes one swab at admission and one at discharge. The swabs were stored at −80 degrees Celsius and sent to a central laboratory for processing. Using patient charts, the study also collected information on antibiotics treatment, including antibiotic type, duration, and route of administration. See Table 1 for an overview of the study details. Written informed consent and consent to publish was obtained from all patients before study enrolment. All collected data was entered de-identified into the central study database which sat in Tel Aviv in accordance with the local rules of personal data privacy protection. The study protocol was reviewed and approved by the Catholic University Ethics Commission in Rome, Italy (protocol P/291/CE/2010 approved on 6.4.2010) and the Clinical Center of Serbia Ethic Committee (protocol 451/34 approved on 18.03.2010). At the site in Romania patient screening for multidrug-resistant bacteria was considered, due to the local epidemiology, a quality improvement intervention and did not require institutional ethical approval.

### Identification of ESBL producing organism carriers

Samples taken at admission were cultured on chromogenic agar (Brilliance ESBL, Oxoid, Basingstoke, UK) to test for ESBL producing Enterobacteriaceae. Characteristically coloured colonies for Enterobacteriaceae were isolated (single colony per colour), replated on blood agar and incubated overnight in air. ESBL status was then confirmed with the double disk diffusion method according to CLSI guidelines . These methods were performed in the laboratories on the respective hospitals (all laboratories were ISO accredited). The above methods are not specific to any single bacterial species, but instead identify ESBL producing specimens of the Enterobacteriaceae family, including Escherichia coli, Klebsiella pneumoniae, or Enterobacter cloacae and others. According to the standard definition by the Centers of Disease Control (CDC/NHSN, 2018), samples taken at admission identify community acquired organisms.

### Quantitative PCR

DNA was extracted from the swab samples using QIAampDNA Stool Mini Kit (Qiagen) and a fixed volume (4 µl) of DNA solution was used as a template for quantitative PCR (qPCR) assays. Two singleplex qPCR assays were conducted, one to assess quantity of blaCTX-M gene family with primers CTX-M-A6 (TGGTRAYRTGGMTBAARGGCA) and CTX-M-A8 (TGGGTRAARTARGTSACCAGAA) (product length, 175 bp) and one targeting a conserved bacterial 16S rRNA gene region bacteria using the following primer set, 16S_E939F (GAATTGACGGGGGCCCGCACAAG) and 16S_1492R (TACGGYTACCTTGTTACGACTT) (product length, 597 bp) to assess total bacterial quantity. Each singleplex qPCR run targeted either 16S gene or blaCTX-M, and included reaction tubes with negative controls and tubes containing a standard (16S or blaCTX-M depending on the run) of different concentrations (eight different dilutions), which also served as a positive control. For details on reaction mix and cycling conditions see Lerner et al., 2013. All reactions were carried out in duplicates, sometimes triplicates, representing technical replicates. The qPCR was performed at the Laboratory of Medical Microbiology, University of Antwerp (ISO accredited). While we did not validate the 16S qPCR measurements through, for example, spiked standard DNA preparations, the ability of qPCR to quantify bacterial amounts in faecal samples has been shown previously (Rinttilä et al., 2004).

### Time series autocorrelation

We first transformed all qPCR measurements onto a log-scale. For all patients and each time point we then computed the mean of the qPCR duplicates (or triplicates) for blaCTX-M and 16S rRNA. To get reliable estimates of autocorrelation, we selected only patients with more than five time points. Separately for the blaCTX-M and 16S rRNA gene data, we computed the first-order autocorrelation (disregarding varying spacing between time points) for each patient, and we averaged these values across the patients. We then simulated serially uncorrelated ‘white noise’ time series, again separately for blaCTX-M and 16S rRNA, with the same length as the patient data and with identical time series mean and variance. Similar to the real data, we computed mean autocorrelations for the simulated data and show their distribution for a large number of such simulations (n = 10,000) together with the observed autocorrelation (Figure 1—figure supplement 1a and b). We also computed the proportion of simulated datasets that showed an average autocorrelation equal to, or larger than, the observed data, and we show those numbers on the arrows in Figure 1—figure supplement 1a and b.

### Estimating observation and process noise in time series

To estimate the amount of observation noise and process noise in the time series we constructed a Bayesian state-space model that included qPCR noise, swab noise, and biological noise. This model is given through:

$\begin{array}{ll}{q}_{i,j,k,g}& \sim N\left({s}_{i,j,g},{\sigma }_{qpcr}\right),\\ {s}_{i,j,g}& \sim N\left({x}_{i,j,g},{\sigma }_{swa{b}_{g}}\right),\\ {x}_{i,j+1,g}& \sim N\left({s}_{i,j,g},{\sigma }_{bi{o}_{g}}\right),\end{array}$
where $i$ denotes a given patient, $j$ denotes a swab (one per time point), $k$ denotes a qPCR measurement (multiple repeats per swab), and $g$ denotes the genetic target, either blaCTX-M or 16S rRNA. The term ${q}_{i,j,k,g}$, represents the measured quantity of genetic target $g$, of the kth qPCR replicate (on a log-scale) from patient $i$, at time point $j$. In addition, there are two hidden-state parameter vectors: ${s}_{i,j,g}$ is the underlying, true sequence abundance of genetic target $g$ that a qPCR assay with 100% efficiency could (in theory) measure at time point $j$ for patient $i$, and ${x}_{i,j,g}$ is the actual gene abundance of genetic target $g$, in the patient at time point $j$ for patient $i$, before the added noise through the swab process and gene extraction. The unobserved variables of interest are ${\sigma }_{qpcr}$, the qPCR machine error (assumed to be the same for blaCTX-M and 16S rRNA), ${\sigma }_{swa{b}_{g}}$, the swab variation of the genetic target $g$, and ${\sigma }_{bi{o}_{g}}$, the variation of genetic target $g$ from biological processes. We assigned improper uniform (‘flat’) priors for the hidden-state parameters and generic weakly informative priors (half-normal, N+(0,1)) for the the noise parameters ${\sigma }_{qpcr}$, ${\sigma }_{swa{b}_{g}}$, and ${\sigma }_{bi{o}_{g}}$. We then fitted this model to the blaCTX-M and 16S rRNA measurements. The posteriors of the three noise parameters are shown in Figure 1—figure supplement 1c, where we expressed each type of noise as a fraction of the total noise. The model was fitted using Stan software (v2.19.1) (Carpenter et al., 2017) and with additional analysis in R (R Development Core Team, 2016), and we sampled 80,000 samples from the posterior with four independent chains and a burn-in period of 20,000 samples. We assessed convergence by checking that R̂ (Gelman and Rubin, 1992) was low (<1.01) for all parameters, and visually by looking at the rank plots for all parameters (shown in Figure 2—figure supplement 1a). For rank plots, posterior draws are ranked across all chains and then ranks are plotted as histograms separately for each chain. A uniform shape of all histograms then indicates that all chains target the same posterior (Vehtari et al., 2019).

### Between and within time series variance

For the estimation of between and within time series variation we used a Bayesian hierarchical model, which accounted for unbalanced sampling between patients. This model used the mean posterior estimates of ${x}_{i,j}$ (actual gene abundance in time point $j$ for patient $i$) from the previous model, and it took the form

$\begin{array}{ll}{x}_{i,j,g}& \sim N\left({\mu }_{i,g},{\sigma }_{within,g}\right),\\ {\mu }_{i,g}& \sim N\left({\mu }_{g},{\sigma }_{between,g}\right),\end{array}$
where $m{u}_{i,g}$ is the mean abundance of genetic target $g$ (blaCTX-M or 16S gene) and patient $i$, around which the log-scaled measurements were assumed to be normally distributed with standard deviation ${\sigma }_{within,g}$, the within time series variation. The mean abundances were assumed to follow a normal distribution with a population mean µ and between-patient standard deviation ${\sigma }_{between}$. We assigned improper uniform priors for the population and the patient means, and generic weakly informative priors for the standard deviations (N+(0,1)). We fitted the model using Stan (v2.19.1), with 80,000 posterior draws after a burn-in period of 500 iterations. The R̂ statistic (<1.01) and MCMC rank histogram plots (Figure 2—figure supplement 1b) were used to assess convergence. Model estimates are shown in Figure 2. To calculate the coefficient of variation for the non-log-scaled blaCTX-M and 16S rRNA measurements, we use the transform described by Koopmans et al., 1964:
${c}_{v}=\sqrt{{e}^{{s}_{ln}^{2}}-1},$
where ${s}_{ln}$ is the estimated standard deviation of the log-scaled data.

### Association of antibiotic treatment and changes in resistance

To study the association between antibiotic treatment and resistance we looked at relative abundance of resistance (blaCTX-M abundance/16S rRNA gene abundance) as a marker of natural selection. First, we computed the changes in relative resistance for every pair of adjacent time points and for each antibiotic we used a binary variable indicating whether or not a given antibiotic was administered between these time points. When an antibiotic treatment was on the same day as a swab, this treatment was allocated to the time interval between this day and the next swab. We first looked at how changes in relative resistance are associated with courses of any antibiotics, then with courses of anti-enterbacteriaceae antibiotics to which carriage of blaCTX-M does not confer direct resistance (colistin, meropenem, ertapenem, imipenem, amoxicillin-clavulanic acid, ampicillin-sulbactam, piperacillin-tazobactam, gentamicin, amikacin, ciprofloxacin, ofloxacin, levofloxacin, tigecycline, doxicycline), and finally with antibiotics that have a broad-spectrum activity and to which blaCTX-M does confer resistance (cefepime, ceftazidime, ceftriaxone, cefotaxime, cefuroxime, amoxicillin, ampicillin). Results are shown in Figure 3, upper panel. We evaluated how likely the observed differences between treatments are under the assumption of no association between treatment and resistance. For this we did a permutation or ‘reshuffling’ experiment: we randomly reassigned (without replacement) the antibiotic treatment labels to the data intervals. We compute the distribution of mean differences from 50,000 permutations and compare this to the observed difference (Figure 3, lower panel).

### Dynamic within-host model

We extended previous modelling approaches to extracting ecological parameters from microbial ecosystem dynamics (Faust and Raes, 2012; Stein et al., 2013) by applying a Bayesian hidden-state model, which featured two layers of hidden-state variables: the unobserved mean of the qPCR measurements, and the unobserved true abundances of blaCTX-M or 16S rRNA in the gut. This model structure allowed us to separate process noise (stochastic effects impacting the gene abundance change from one day to the next) from observation noise (stochastic effects impacting the swab efficiency, DNA extraction or the qPCR measurements) and also to account for different spacing between measured time points.

We analysed antibiotic treatment separately by type and route of administration, only including treatments that occurred in five or more patients (amoxicillin-clavulanic acid (iv), piperacillin-tazobactam (iv), cefuroxime (iv), ceftriaxone (iv), meropenem (iv), imipenem (iv), ciprofloxacin (iv), ciprofloxacin (or), amikacin (iv), metronidazole (iv)). Since the blaCTX-M and 16S rRNA measurements are expected to be proportional to the absolute abundance of the resistance gene and bacterial load respectively, the ratio blaCTX-M/16S is a measure of the relative abundance of the blaCTX-M gene in the gut microbiota. Positive or negative selective effects by antibiotics on blaCTX-M mediated resistance are expected to cause shifts in bacteria carrying blaCTX-M versus non-carriers. As a result they affect blaCTX-M/16S, but quantifying their effects on absolute blaCTX-M abundance is important for predicting extinction and persistence of the blaCTX-M gene. Under the assumption that 16S rRNA gene abundance is independent of antibiotic treatment, variation in 16S rRNA would be caused mainly by the swab procedure and DNA extraction (and other steps in the protocol), and it could be used to normalise blaCTX-M abundance. However, as we found in Figure 4, certain antibiotic treatments were associated with changes in 16S rRNA abundance. Thus, we used a dynamic model that explicitly modelled both antibiotic effects on 16S rRNA and on blaCTX-M/16S, from which the effects on blaCTX-M could then be computed.

Studying the standard deviation between qPCR measurement repeats as a function of the mean, we observed that qPCR variation remained relatively stable over five orders of magnitude of the mean measurement (from 1.5 to 6.5 on the log-scale), but it increased quickly for lower magnitudes (Figure 4—figure supplement 1). In the Bayesian model for different sources of variation described above, the parameter ${\sigma }_{qpcr}$ assumed that the qPCR uncertainty is the same across measurements. Here, we aimed to account for the fact that low measurements of gene copy numbers have higher uncertainty. We fitted a smooth spline (choosing five degrees of freedom) to the qPCR measurements (red line in Figure 4—figure supplement 1). This let us assign an estimated qPCR standard deviation to every set of qPCR repeats. We provided those estimates as data to the Bayesian model. This allowed us to use all qPCR measurements, including extremely low values, without removing any data points from the analysis. Our model then took the form:

$\begin{array}{ll}{q}_{i,j,g,k}& \sim N\left({s}_{i,j,g},{\sigma }_{qpc{r}_{i,j,g}}\right),\\ {s}_{i,j,g=16S}& \sim N\left({x}_{i,j,b=16S},{\sigma }_{swab}\right),\\ {x}_{i,j,b=ratio}& ={s}_{i,j,g=CTX-M}-{s}_{i,j,g=16S},\\ {x}_{i,j+1,b}& \sim N\left({x}_{i,j,b}+f\left(abx{\right)}_{i,j,b},\sqrt{{t}_{j+1}-{t}_{j}}{\sigma }_{bio,b}\right),\\ f\left(abx{\right)}_{i,j,b}& =\sum _{{t}_{j}}^{{t}_{j+1}-1}\left({a}_{b}+\sum _{z=1}^{{n}_{z}}{c}_{z,b}{y}_{z,t}\right),\end{array}$
where $g$ denotes either blaCTX-M or 16S rRNA, and $b$ denotes either relative resistance (blaCTX-M/16S rRNA) or 16S rRNA. Then, ${q}_{i,j,g,k}$ is the $k$-th qPCR result (log-scaled) of patient $i$, measured in the sample with index $j$, and genetic target $g$ (blaCTX-M or 16S rRNA). The qPCR standard deviation for sample $j$ and patient $i$ is given through ${\sigma }_{qpc{r}_{i,j,g}}$, and it is estimated as described above. ${s}_{i,j,g}$ is the mean of the qPCR measurements in sample $j$ of patient $i$, and genetic target $g$, and ${x}_{i,j,b=16S}$ is the 16S sequence abundance that is actually present at time point $j$ in the gut of patient $i$. The error introduced by the swab procedure, DNA extraction etc. is given through ${\sigma }_{swab}$, and we assume that this error causes the same perturbations to blaCTX-M and to 16S, such that this error cancels out when computing the ratio of blaCTX-M and 16S rRNA ${x}_{i,j,b=ratio}$, which on a log-scale is computed as the difference (${s}_{i,j,g=CTX-M}-{s}_{i,j,g=16S}$). Further, ${t}_{j}$ denotes the calendar day of sample $j$ in patient $i$, and ${t}_{j+1}$ denotes the calendar day of the following sample in the same patient $i$. The daily biological variability of the blaCTX-M/16S ratio and of 16S are given through ${\sigma }_{bio,b=ratio}$ and ${\sigma }_{bio,b=16S}$, respectively, with $\sqrt{{t}_{j+1}-{t}_{j}}$ adjusting the expected random walk variation by the number of days between observations (Lemons and Langevin, 2002). The ecological dynamics are modelled with the function $f{\left(abx\right)}_{i,j,b}$, which is the change in the expected value of ${x}_{i,j,b}$ between sample $j$ and $j+1$. In the definition of $f{\left(abx\right)}_{i,j,b}$ in line 5 of Equation 4, ${a}_{b}$ denotes the neutral growth or loss of $g$, ${c}_{z,b}$ denotes the effect of antibiotic $z$ on $b$, and ${y}_{z,t}$ is a boolean variable indicating whether or not antibiotic $z$ was given on day $t$. The term inside the bracket takes, for a single calendar day $t$, the neutral growth/loss term (${a}_{b}$) and adds to this the summed effect of all antibiotics given on day $t$. This is computed for each calendar day from ${t}_{j}$ until a day before the subsequent sample (${t}_{j+1}-1$). The effects of all of these days are then summed up. Note, that ${x}_{i,j,b}$ denotes the abundance of 16S rRNA, or the relative abundance of blaCTX-M/16S rRNA, on the log-scale. Exponentiating this variable gives the copy numbers (or copy number ratio) on the real scale. Therefore, summing all effects on the scale of ${x}_{i,j,b}$ is equivalent to multiplying the exponentiated effects on the scale of copy numbers. For example, consider that genetic target $b=16S$ in patient $i=1$ has at the time of sample $j=1$ an abundance of ${10}^{{x}_{i=1,j=1,b=16S}}=100$ copy numbers and a neutral trend of ${a}_{b=16S}=-0.5$. Suppose on the day of this sample (${t}_{i=1,j=1}$) two antibiotics $z=1$ and $z=2$ are given with effects ${c}_{z=1,b=16S}=+0.5$ and ${c}_{z=2,b=16S}=+0.1$, then one day after this sample the genetic target abundance has an expected abundance of $100*{10}^{-0.5}*{10}^{+0.5}*{10}^{+0.1}=126$.

In this model, ${q}_{i,j,g,k}$ and ${y}_{z,t}$ correspond to measured data, ${\sigma }_{qpc{r}_{i,j,g}}$ is computed from this data (see above). All other parameters are estimated: the hidden-state variables are ${s}_{i,j,g}$, ${x}_{i,j,g=16S}$, and ${x}_{i,j,g=ratio}$, the noise variables are ${\sigma }_{swab}$, ${\sigma }_{bio,g=16S}$, and ${\sigma }_{bio,g=ratio}$, and finally the variables describing the ecology are ${a}_{g}$, ${c}_{z,g}$. The model has three likelihood functions. The first (line 1) applies to each single qPCR result and it relates repeat qPCR measurements to their variability and underlying mean. The second (line 2) relates qPCR means of 16S rRNA to their variability and the underlying 16S rRNA gene abundance. The third likelihood (line 4) applies to all sample pairs where a previous and following sample exists from the same patient. This likelihood relates changes in underlying 16S rRNA abundance or underlying blaCTX-M/16S rRNA to the parameters ${a}_{b}$ and ${c}_{z,b}$, and ${\sigma }_{bio}$. All parameters (including all hidden-states) are estimated simultaneously. We can express the posterior distribution over all estimated parameters ($\mathrm{\Theta }$), which is conditional on the set of all data ($D$), and the prior over the parameter space $p\left(\mathrm{\Theta }\right)$. For readability, here we only keep necessary subscripts:

$\begin{array}{ll}p\left(\mathrm{\Theta }|D\right)& =p\left(s,x,a,c,{\sigma }_{swab},{\sigma }_{bio}|q,{\sigma }_{qpcr},y\right),\\ p\left(\mathrm{\Theta }|D\right)& \propto p\left(\mathrm{\Theta }\right)\prod _{i,j,g,k}p\left({q}_{i,j,g,k}|{s}_{i,j,g},{\sigma }_{qpc{r}_{i,j,g}}\right)\prod _{i,j}p\left({s}_{i,j,g=16S}|{x}_{i,j,b=16S},{\sigma }_{swab}\right)\\ & \prod _{i,j,b}p\left({x}_{i,j+1,b}|{c}_{b},{a}_{b},{\sigma }_{bio,b},{x}_{i,j,b},y\right).\end{array}$

On the hidden-state variables we assigned improper uniform priors, on the standard deviations describing swab and biological variability we assigned standard weakly informative priors (N+(0,1)), and on the antibiotic effects ($c$) we assigned conservative priors of the form N(0, 0.1). We fitted the model using Stan software (v2.19.1), and we sampled 80,000 samples from the posterior with four independent chains after a burn-in phase of 10,000 samples. The marginal posterior draws of the ${c}_{z,g}$ parameters are exponentiated to be on the scale of gene counts. Subtracting one allows us to express daily effects as percent of increase or decrease relative to the previous day (Figure 4). We also show marginal posterior distributions together with prior distributions for the ${c}_{z,g}$ parameters (Figure 4—figure supplement 2). Diagnostic plots for the MCMC sampling of the ${c}_{z,g}$ parameters are shown in Figure 4—figure supplement 3.

We compared our model as given through model definition four to a model without antibiotic effects (all $c$ parameters set to zero). The number of patients treated with the same antibiotic is too small to perform cross-validation (iteratively fitting the model to all time series while leaving out data from one patient, which is then used to assess model predictions). Therefore, we used an efficient approximation of Bayesian leave-one-out cross-validation using Pareto-smoothed importance sampling (Vehtari et al., 2017).

We forward simulated blaCTX-M data using the dynamic model above and the posterior distributions from the model fit. We added to the model a threshold below which the blaCTX-M gene becomes extinct or at least undetectable. According to a study of returning European travelers to Southeast Asia, ESBL carriers lose detectable resistant bacteria after a median of 30 days (Arcilla et al., 2017). Accordingly, we simulated blaCTX-M time series without antibiotic treatment and chose an extinction threshold (0.25 blaCTX-M copy numbers) that achieved the same median extinction time. We then used this model to repeatedly (2000 times) simulate blaCTX-M carriage durations, with each simulation using a new draw from the parameter posterior. The resulting distribution of carriage times contains both the uncertainty in the parameter estimates and uncertainty from the Markov process (Figure 5, right-hand side). We also draw a set of parameter values from the posterior to simulate blaCTX-M carriage durations repeatedly (200 times) with the same parameters, and taking the median carriage duration to remove Markov process uncertainty. We repeated this for 300 draws of parameters (Figure 5, left-hand side). We used both of the above methods to simulate carriage time under different alternative antibiotic treatments. The resulting distributions are shown in Figure 5.

## Funding Information

This paper was supported by the following grants:

http://dx.doi.org/10.13039/501100000265Medical Research Council to Rene Niehus, Ben S Cooper.
http://dx.doi.org/10.13039/100004440Wellcome TrustMajor Overseas Programme, 106698/Z/14/Z to Ben S Cooper.
http://dx.doi.org/10.13039/100011102European Union 7th Framework ProgrammeR-GNOSIS Integrated project, grant agreement number 241796 to Rene Niehus, Esther van Kleef, Agata Turlej-Rogacka, Christine Lammens, Yehuda Carmeli, Herman Goossens, Evelina Tacconelli, Biljana Carevic, Surbhi Malhotra-Kumar, Ben S Cooper.
http://dx.doi.org/10.13039/501100000278Department for International DevelopmentMR/K006924/1 to Rene Niehus, Ben S Cooper.
http://dx.doi.org/10.13039/501100001349National Medical Research Council to Yin Mo.

## Acknowledgements

We thank Jonas Schluter, Marc Lipsitch, Thomas Crellen, and Pierre Jacob for valuable feedback along the way. RN as well as the study were supported by funding from the European Community’s R-GNOSIS Integrated project (FP7/2007-2013) under grant agreement number 241796. RN and BSC were also supported by the Medical Research Council and Department for International Development (grant number MR/K006924/1). BSC works within the Wellcome Trust Major Overseas Programme in SE Asia (grant number 106698/Z/14/Z). The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

### Competing interests

Reviewing Editor, eLife.
No competing interests declared.

### Author contributions

Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology.
Conceptualization, Resources, Data curation, Formal analysis, Validation, Investigation, Visualization, Methodology, Project administration.
Conceptualization, Validation, Methodology.
Resources, Validation, Investigation, Methodology.
Resources, Data curation, Validation, Methodology.
Conceptualization, Resources, Supervision, Funding acquisition, Validation, Investigation, Methodology, Project administration.
Conceptualization, Supervision, Funding acquisition, Validation, Investigation, Methodology, Project administration.
Conceptualization, Supervision, Validation, Investigation, Methodology, Project administration.
Resources, Investigation, Methodology.
Resources, Investigation.
Conceptualization, Resources, Supervision, Funding acquisition, Validation, Investigation, Methodology, Project administration.
Conceptualization, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Methodology, Project administration.

### Ethics

Clinical trial registration NCT01208519.
Human subjects: Written informed consent and consent to publish was obtained from all patients before study enrolment. All collected data was entered de-identified into the central study database which sat in Tel Aviv in accordance with the local rules of personal data privacy protection. The study protocol was reviewed and approved by the Catholic University Ethics Commission in Rome, Italy (protocol P/291/CE/2010 approved on 6.4.2010) and the Clinical Center of Serbia Ethic Committee (protocol 451/34 approved on 18.03.2010). At the site in Romania patient screening for multidrug-resistant bacteria was considered, due to the local epidemiology, a quality improvement intervention and did not require institutional ethical approval.

## Data availability

The Bayesian model code in Stan, the R code to fit the model, and the data in Rdata format as well as csv format are deposited in Dryad under https://doi.org/10.5061/dryad.8vf034t.

The following dataset was generated:

Niehus R, van Kleef E, Yin M, Turlej-Rogacka A, Lammens C, Carmeli Y, Goossens H, Tacconelli E, Carevic B, Malhotra-Kumar S, Cooper BS 2020. Quantifying antibiotic impact on within-host dynamics of extended-spectrum beta-lactamase resistance in hospitalized patientsDryad Digital Repository 10.5061/dryad.8vf034t

## References

1

Arcilla MS, van Hattem JM, Haverkate MR, Bootsma MCJ, van Genderen PJJ, Goorhuis A, Grobusch MP, Lashof AMO, Molhoek N, Schultsz C, Stobberingh EE, Verbrugh HA, de Jong MD, Melles DC, Penders J 2017. . Import and spread of extended-spectrum β-lactamase-producing Enterobacteriaceae by international travellers (COMBAT study): a prospective, multicentre cohort study. The Lancet Infectious Diseases 17: , pp.78-85, doi: 10.1016/S1473-3099(16)30319-X

2

Balaban NQ, Helaine S, Lewis K, Ackermann M, Aldridge B, Andersson DI, Brynildsen MP, Bumann D, Camilli A, Collins JJ, Dehio C, Fortune S, Ghigo JM, Hardt WD, Harms A, Heinemann M, Hung DT, Jenal U, Levin BR, Michiels J, Storz G, Tan MW, Tenson T, Van Melderen L, Zinkernagel A 2019. . Definitions and guidelines for research on antibiotic persistence. Nature Reviews Microbiology 17: , pp.441-448, doi: 10.1038/s41579-019-0196-3

3

Blanquart F 2019. . Evolutionary epidemiology models to predict the dynamics of antibiotic resistance. Evolutionary Applications 12: , pp.365-383, doi: 10.1111/eva.12753

4

Bonnet R 2004. . Growing group of Extended-Spectrum -Lactamases: the CTX-M Enzymes. Antimicrobial Agents and Chemotherapy 48: , pp.1-14, doi: 10.1128/AAC.48.1.1-14.2004

5

Buelow E, Gonzalez TB, Versluis D, Oostdijk EA, Ogilvie LA, van Mourik MS, Oosterink E, van Passel MW, Smidt H, D'Andrea MM, de Been M, Jones BV, Willems RJ, Bonten MJ, van Schaik W 2014. . Effects of selective digestive decontamination (SDD) on the gut resistome. Journal of Antimicrobial Chemotherapy 69: , pp.2215-2223, doi: 10.1093/jac/dku092

6

Buelow E, Bello González TDJ, Fuentes S, de Steenhuijsen Piters WAA, Lahti L, Bayjanov JR, Majoor EAM, Braat JC, van Mourik MSM, Oostdijk EAN, Willems RJL, Bonten MJM, van Passel MWJ, Smidt H, van Schaik W 2017. . Comparative gut Microbiota and resistome profiling of intensive care patients receiving selective digestive tract decontamination and healthy subjects. Microbiome 5: 88, doi: 10.1186/s40168-017-0309-z

7

Carpenter B, Gelman A, Hoffman MD, Lee D, Goodrich B, Betancourt M, Brubaker M, Guo J, Li P, Riddell A 2017. . Stan : a probabilistic programming language. Journal of Statistical Software 76: 1, doi: 10.18637/jss.v076.i01

8

CDC/NHSN 2018. Surveillance Definitions for Specific Types of Infections CDC https://www.cdc.gov/nhsn/pdfs/pscmanual/17pscnosinfdef_current.pdf,

9

Coque TM, Baquero F, Canton R 2008. . Increasing prevalence of ESBL-producing Enterobacteriaceae in Europe. Euro Surveillance : Bulletin Europeen Sur Les Maladies Transmissibles = European Communicable Disease Bulletin 13: 19044

10

Davies NG, Flasche S, Jit M, Atkins KE 2019. . Within-host dynamics shape antibiotic resistance in commensal bacteria. Nature Ecology & Evolution 3: , pp.440-449, doi: 10.1038/s41559-018-0786-x

11

Donnenberg MS 1979. EnterobacteriaceaeBennett M. DPrinciples and Practice of Infectious Diseases Elsevier, pp.2503-2517

12

Donskey CJ 2004. . The role of the intestinal tract as a reservoir and source for transmission of nosocomial pathogens. Clinical Infectious Diseases 39: , pp.219-226, doi: 10.1086/422002

13

Dubinsky-Pertzov B, Temkin E, Harbarth S, Fankhauser-Rodriguez C, Carevic B, Radovanovic I, Ris F, Kariv Y, Buchs NC, Schiffer E, Cohen Percia S, Nutman A, Fallach N, Klausner J, Carmeli Y, R-GNOSIS WP4 Study Group 2019. . Carriage of Extended-spectrum Beta-lactamase-producing Enterobacteriaceae and the risk of surgical site infection after colorectal surgery: a prospective cohort study. Clinical Infectious Diseases 68: , pp.1699-1704, doi: 10.1093/cid/ciy768

14

Faust K, Raes J 2012. . Microbial interactions: from networks to models. Nature Reviews Microbiology 10: , pp.538-550, doi: 10.1038/nrmicro2832

15

Garot D, Respaud R, Lanotte P, Simon N, Mercier E, Ehrmann S, Perrotin D, Dequin PF, Le Guellec C 2011. . Population pharmacokinetics of ceftriaxone in critically ill septic patients: a reappraisal. British Journal of Clinical Pharmacology 72: , pp.758-767, doi: 10.1111/j.1365-2125.2011.04005.x

16

Gelman A, Rubin DB 1992. . Inference from iterative simulation using multiple sequences. Statistical Science 7: , pp.457-472, doi: 10.1214/ss/1177011136

17

Gibson MK, Wang B, Ahmadi S, Burnham CA, Tarr PI, Warner BB, Dantas G 2016. . Developmental dynamics of the preterm infant gut Microbiota and antibiotic resistome. Nature Microbiology 1: 16024, doi: 10.1038/nmicrobiol.2016.24

18

Jumbe N, Louie A, Leary R, Liu W, Deziel MR, Tam VH, Bachhawat R, Freeman C, Kahn JB, Bush K, Dudley MN, Miller MH, Drusano GL 2003. . Application of a mathematical model to prevent in vivo amplification of antibiotic-resistant bacterial populations during therapy. Journal of Clinical Investigation 112: , pp.275-285, doi: 10.1172/JCI200316814

19

Knight GM, Costelloe C, Deeny SR, Moore LSP, Hopkins S, Johnson AP, Robotham JV, Holmes AH 2018. . Quantifying where human acquisition of antibiotic resistance occurs: a mathematical modelling study. BMC Medicine 16: 137, doi: 10.1186/s12916-018-1121-8

20

Koopmans LH, Owen DB, Rosenblatt JI 1964. . Confidence intervals for the coefficient of variation for the normal and log normal distributions. Biometrika 51: , pp.25-32, doi: 10.1093/biomet/51.1-2.25

21

Lautenbach E, Strom BL, Bilker WB, Patel JB, Edelstein PH, Fishman NO 2001. . Epidemiological investigation of fluoroquinolone resistance in infections due to extended-spectrum beta-lactamase-producing Escherichia coli and Klebsiella pneumoniae. Clinical Infectious Diseases : An Official Publication of the Infectious Diseases Society of America 33: , pp.1288-1294, doi: 10.1086/322667

22

Lemons DS, Langevin P 2002. An Introduction to Stochastic Processes in Physics JHU Press

23

Lerner A, Romano J, Chmelnitsky I, Navon-Venezia S, Edgar R, Carmeli Y 2013. . Rectal swabs are suitable for quantifying the carriage load of KPC-producing carbapenem-resistant Enterobacteriaceae. Antimicrobial Agents and Chemotherapy 57: , pp.1474-1479, doi: 10.1128/AAC.01275-12

24

Lipsitch M, Bergstrom CT, Levin BR 2000. . The epidemiology of antibiotic resistance in hospitals: Paradoxes and prescriptions. PNAS 97: , pp.1938-1943, doi: 10.1073/pnas.97.4.1938

25

Lipsitch M, Samore MH 2002. . Antimicrobial use and antimicrobial resistance: a population perspective. Emerging Infectious Diseases 8: , pp.347-354, doi: 10.3201/eid0804.010312

26

Livermore DM, Brown DF 2001. . Detection of beta-lactamase-mediated resistance. Journal of Antimicrobial Chemotherapy 48 Suppl 1: , pp.59-64, doi: 10.1093/jac/48.suppl_1.59

27

MacLean RC, Hall AR, Perron GG, Buckling A 2010. . The population genetics of antibiotic resistance: integrating molecular mechanisms and treatment contexts. Nature Reviews Genetics 11: , pp.405-414, doi: 10.1038/nrg2778

28

Maier L, Pruteanu M, Kuhn M, Zeller G, Telzerow A, Anderson EE, Brochado AR, Fernandez KC, Dose H, Mori H, Patil KR, Bork P, Typas A 2018. . Extensive impact of non-antibiotic drugs on human gut Bacteria. Nature 555: , pp.623-628, doi: 10.1038/nature25979

29

Masci JR 2005. Clinical Infectious Diseases Mandell, Douglas, and Bennett’s Principles and Practice of Infectious Diseases An Official Publication of the Infectious Diseases Society of America

30

Masterton R, Drusano G, Paterson DL, Park G 2003. . Appropriate antimicrobial treatment in nosocomial infections-the clinical challenges. Journal of Hospital Infection 55 Suppl 1: , pp.1-12, doi: 10.1016/S0195-6701(03)00294-9

31

Meletiadis J, Turlej-Rogacka A, Lerner A, Adler A, Tacconelli E, Mouton JW 2017. . Amplification of antimicrobial resistance in gut flora of patients treated with ceftriaxone. Antimicrobial Agents and Chemotherapy 61: e00473-17, doi: 10.1128/AAC.00473-17

32

Nahata MC, Barson WJ 1985. . Ceftriaxone: a third-generation cephalosporin. Drug Intelligence & Clinical Pharmacy 19: , pp.900-906, doi: 10.1177/106002808501901203

33

Neu HC, Fu KP 1978. . Cefuroxime, a Beta-Lactamase-Resistant cephalosporin with a broad spectrum of Gram-Positive and -Negative activity. Antimicrobial Agents and Chemotherapy 13: , pp.657-664, doi: 10.1128/AAC.13.4.657

34

Nguyen TT, Guedj J, Chachaty E, de Gunzburg J, Andremont A, Mentré F 2014. . Mathematical modeling of bacterial kinetics to predict the impact of antibiotic colonic exposure and treatment duration on the amount of resistant enterobacteria excreted. PLOS Computational Biology 10: e1003840, doi: 10.1371/journal.pcbi.1003840

35

Paterson DL 2000. . Recommendation for treatment of severe infections caused by Enterobacteriaceae producing extended-spectrum beta-lactamases (ESBLs). Clinical Microbiology and Infection 6: , pp.460-463, doi: 10.1046/j.1469-0691.2000.00107.x

36

Paterson DL, Singh N, Rihs JD, Squier C, Rihs BL, Muder RR 2001. . Control of an outbreak of infection due to extended-spectrum beta-lactamase--producing Escherichia coli in a liver transplantation unit. Clinical Infectious Diseases : An Official Publication of the Infectious Diseases Society of America 33: , pp.126-128, doi: 10.1086/320882

37

Paterson DL 2006. . Resistance in gram-negative Bacteria: enterobacteriaceae. The American Journal of Medicine 119: , pp.S20-S28, doi: 10.1016/j.amjmed.2006.03.013

38

R Development Core Team 2016. R: A Language and Environment for Statistical ComputingVienna, Austria R Foundation for Statistical Computing http://www.r-project.org,

39

Relman DA, Lipsitch M 2018. . Microbiome as a tool and a target in the effort to address antimicrobial resistance. PNAS 115: , pp.12902-12910, doi: 10.1073/pnas.1717163115

40

Rinttilä T, Kassinen A, Malinen E, Krogius L, Palva A 2004. . Development of an extensive set of 16S rDNA-targeted primers for quantification of pathogenic and indigenous Bacteria in faecal samples by real-time PCR. Journal of Applied Microbiology 97: , pp.1166-1177, doi: 10.1111/j.1365-2672.2004.02409.x

41

Schwaber MJ, Navon-Venezia S, Schwartz D, Carmeli Y 2005. . High levels of antimicrobial coresistance among Extended-Spectrum- -Lactamase-Producing Enterobacteriaceae. Antimicrobial Agents and Chemotherapy 49: , pp.2137-2139, doi: 10.1128/AAC.49.5.2137-2139.2005

42

Sorlózano A, Gutiérrez J, Romero JM, de Dios Luna J, Damas M, Piédrola G 2007. . Activity in vitro of twelve antibiotics against clinical isolates of extended-spectrum beta-lactamase producing Escherichia coli. Journal of Basic Microbiology 47: , pp.413-416, doi: 10.1002/jobm.200710318

43

Stein RR, Bucci V, Toussaint NC, Buffie CG, Rätsch G, Pamer EG, Sander C, Xavier JB 2013. . Ecological modeling from time-series inference: insight into dynamics and stability of intestinal Microbiota. PLOS Computational Biology 9: e1003388, doi: 10.1371/journal.pcbi.1003388

44

Tacconelli E, Carrara E, Savoldi A, Harbarth S, Mendelson M, Monnet DL, Pulcini C, Kahlmeter G, Kluytmans J, Carmeli Y, Ouellette M, Outterson K, Patel J, Cavaleri M, Cox EM, Houchens CR, Grayson ML, Hansen P, Singh N, Theuretzbacher U, Magrini N, WHO Pathogens Priority List Working Group 2018. . Discovery, research, and development of new antibiotics: the WHO priority list of antibiotic-resistant Bacteria and tuberculosis. The Lancet Infectious Diseases 18: , pp.318-327, doi: 10.1016/S1473-3099(17)30753-3

45

Tacconelli E, Górska A, De Angelis G, Lammens C, Restuccia G, Schrenzel J, Huson DH, Carević B, Preoţescu L, Carmeli Y, Kazma M, Spanu T, Carrara E, Malhotra-Kumar S, Gladstone BP 2020. . Estimating the association between antibiotic exposure and colonization with extended-spectrum β-lactamase-producing Gram-negative Bacteria using machine learning methods: a multicentre, prospective cohort study. Clinical Microbiology and Infection 26: , pp.87-94, doi: 10.1016/j.cmi.2019.05.013

46

Tedijanto C, Olesen SW, Grad YH, Lipsitch M 2018. . Estimating the proportion of bystander selection for antibiotic resistance among potentially pathogenic bacterial flora. PNAS 115: , pp.E11988-E11995, doi: 10.1073/pnas.1810840115

47

Udekwu KI, Parrish N, Ankomah P, Baquero F, Levin BR 2009. . Functional relationship between bacterial cell density and the efficacy of antibiotics. Journal of Antimicrobial Chemotherapy 63: , pp.745-757, doi: 10.1093/jac/dkn554

48

Valverde A, Coque TM, Sánchez-Moreno MP, Rollán A, Baquero F, Cantón R 2004. . Dramatic increase in prevalence of fecal carriage of extended-spectrum beta-lactamase-producing Enterobacteriaceae during nonoutbreak situations in Spain. Journal of Clinical Microbiology 42: , pp.4769-4775, doi: 10.1128/JCM.42.10.4769-4775.2004

49

Vehtari A, Gelman A, Gabry J 2017. . Practical bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing 27: , pp.1413-1432, doi: 10.1007/s11222-016-9696-4

50

Vehtari A, Gelman A, Simpson D, Carpenter B, Bürkner PC 2019. . Rank-normalization, folding, and localization: an improved R̂ for assessing convergence of MCMC. arXiv https://arxiv.org/pdf/1903.08008.pdf,

51

Watanabe S 2010. . Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research 11: , pp.3571-3594

52

Webb GF, D'Agata EM, Magal P, Ruan S 2005. . A model of antibiotic-resistant bacterial epidemics in hospitals. PNAS 102: , pp.13343-13348, doi: 10.1073/pnas.0504053102

53

Winokur PL, Canton R, Casellas JM, Legakis N 2001. . Variations in the prevalence of strains expressing an extended-spectrum beta-lactamase phenotype and Characterization of isolates from Europe, the Americas, and the western pacific region. Clinical Infectious Diseases 32 Suppl 2: , pp.S94-S103, doi: 10.1086/320182

54

Woerther PL, Micol JB, Angebault C, Pasquier F, Pilorge S, Bourhis JH, de Botton S, Gachot B, Chachaty E 2015. . Monitoring antibiotic-resistant enterobacteria faecal levels is helpful in predicting antibiotic susceptibility of Bacteraemia isolates in patients with haematological malignancies. Journal of Medical Microbiology 64: , pp.676-681, doi: 10.1099/jmm.0.000078

55

Zhang L, Huang Y, Zhou Y, Buckley T, Wang HH 2013. . Antibiotic administration routes significantly influence the levels of antibiotic resistance in gut Microbiota. Antimicrobial Agents and Chemotherapy 57: , pp.3659-3666, doi: 10.1128/AAC.00670-13

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

Acceptance summary:

The dynamics of resistance in response to antibiotic treatment is of central importance in understanding resistance emergence more broadly. This paper develops a data-driven model to describe the within-host dynamics of Enterobacteriaceae that produce extended-spectrum β-lactamase. Supported by the model, the authors compare different antibiotics' effects on resistance.

Decision letter after peer review:

Thank you for submitting your article "Quantifying antibiotic impact on within-patient dynamics of extended-spectrum β-lactamase resistance" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Wendy Garrett as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Amy Mathers (Reviewer #2); Lulla Opatowski (Reviewer #3).

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

Summary:

This work asks what the dynamic consequences of antibiotic perturbations are to the human gut flora, and addresses this question with a novel mathematical model that is fit to data. The authors use the model to estimate the impact of antibiotic treatment on the abundance of resistance (through blaCTX-M) and on the duration of carriage of resistance. Overall, the manuscript represents a potentially important contribution to the field. However, there is some T data overlap in a previous study by the authors. and we have identified a few essential revisions to resolve this issue and that need to be addressed as well.

Essential revisions:

The reviewers brought up the previous publication, Meletiadis, J. et al., 2017, in which a good portion of the cases were reported, and where some of the same conclusions were reached.

The authors should clarify the novel contributions of this work in light of the previous manuscript. It seems that the key claimed results here are (1) variability, Figures 1, 2; (2) effect of treatment on relative resistance; (3) effect of treatment on abundance and (4) the model with fitting, and consequently the estimates re duration of carriage. Each has limitations.

Please set these results, with the novelty and limitations, in context.

The modelling is clearly new and I think it is important. In that direction, the clarity should be improved. I particularly noted the very long paragraph in subsection “Dynamic within-host model”, and the benefit that explaining the model and its assumptions (and limitations) for non-statisticians, non-modellers would likely bring.

Please clarify re interpretation of 16s, abundance as noted by reviewers.

Please improve Figures 1 and 2. While "hieroglyphics" may go a bit far, I agree that Figure 1 is visually striking but hard to see and interpret. Perhaps in Figure 1 you could group the time curves by similarity, or present an expanded view of some representative curves, with the full set in an appendix. Or something like that.

Reviewer #1:

The motivation for this study and its goals are certainly important and timely. The collateral effects of antibiotic treatment on the distribution of bacteria species and strains in the enteric flora (microbiome) of humans and the incidence of antibiotic resistance genes and resistance encoding plasmids are subjects of considerable importance, epidemiologically as well as clinically. In their Introduction, the authors do a fine job of presenting this and particularly so for the extended-spectrum β-lactamase (ESBL)-producing Enterobacteriaceae and the blaCTX family of resistance genes, that are the focus of their study.

This investigation certainly has the virtue of being extensive; fecal samples were taken sequentially from 133 hospitalized patients from three countries (Romania, Serbia and Italy) for a median of five samples from each patient. This virtue does have a downside; these patients were hospitalized for different reasons and treated orally of intravenously with 10 different antibiotics. The results presented suggest that some antibiotics are more likely than others to affect the distribution and abundance of these enteric bacteria and the blaCRX-M genes, their study is correlative rather than mechanistic. From the results presented, it's not clear how the different modes and purpose of treatment the within-patient dynamics ESBL resistance. From what they say in their Discussion, the authors are very much aware of this limitation of this study.

Although their observation some antibiotics are more likely than others to affect the enteric microbiome and the frequency of blaCRX-M genes could have epidemiological and clinical implication, they don't really consider these potential implications of their results. It's also not clear how this information can be used. Presumably the choice of antibiotics, mode and frequency of administration in these hospitalized patients is based on the nature of the infection and information about the susceptibility of the target bacteria to the drugs employed for treatment. While in an ideal world, consideration should to the collateral effects of treatment; we are far from that ideal world. This is particularly so for hospitalized patients. It would be of some interest to do an analogous study for the treatment of community-acquired infections (about 90% of antibiotic use in humans), where information about the collateral effects of treatment with different antibiotics may have broader and implications that could be implemented.

Reviewer #2:

This is an interesting and important piece of work demonstrating the correlation between certain antimicrobial exposure in patients and associated increase in abundance of blaCTX-M. The gene increase is then normalized to 16s abundance for overall internal QC for extraction efficiency and bacterial content. The work appears to be very well done and the methods and approach to analysis appear to support the conclusions of the authors. The figures are easy to interpret and additive to the manuscript (not sure if Figure 1 is critical but it is interesting to look at the variability across the data set and would favor keeping). Also the discussion includes a portion focused on relative biomass and bacterial presence in the gut of patients on antimicrobials but it is not clear they would be able to infer this information from the variability in swab collection and extraction efficiency. Would favor revising the discussion to address this issue.

One concern is the amount this manuscript overlaps with another previously published manuscript (Meletiadis, J. et al., 2017). The authors state in the Introduction section that they used a subset of the data but in reviewing the other manuscript it appears that there are 133 patients in this analysis and 122 from the other description. They also used very similar methodology by normalizing the blaCTX-M abundance with 16s. I think it would be helpful for the authors to further explain how this differs from prior work and is additive to the literature. The conclusions from that manuscript were almost identical as well. Would favor a re-write to address this existing work and contrast the submitted manuscript.

Reviewer #3:

Niehus and colleagues present an original study in which they investigate the impact of different antibiotics on ESBL resistance within-host in Enterobacteriaceae carriers. The data comes from 133 patients from Romania, Serbia and Italy and consists in longitudinal series of rectal swabs (with a median of 5 swabs per patient). Swabs were analysed in order to quantify the abundance of blaCTX-M and 16SrRNA. For both indicators, the within- and between-host variabilities, and the within-host dynamics are analysed using Bayesian state-space models. Parameters associated to the selection exerted by different classes of antibiotics on the dynamics of blaCTX-M, 16SrRNA and blaCTX-M/16SrRNA are estimated.

This is a very exciting and well written article. This is, to my knowledge, the first study attempting to analyse, using hypothesis-driven mechanistic models, the dynamics of within-host resistance genes and the impact of various classes of antibiotics. I believe that the study and results presented here represent a major contribution to the field.

– It would be good for the reader to detail and clarify how each indicator should be interpreted. My understanding is for example that a decrease in blaCTX-M/16SrRNA would more or less be interpreted as a decrease in resistance rate. However, it is not clear to me how important is the contribution of Enterobacteriaceae species to the global quantity of 16SrRNA.

– My main comment is related to the simulation study. It would be good to provide a validation of the model before running simulations, for example by doing out of fit assessment? Is that possible? Would you have enough power?

– It would be interesting to estimate whether the antibiotic classes have different delays of impact.

– I am not an expert but my understanding is that the aptitude to induce degradation of the flora may not be fully characterized by the "broad spectrum" characteristic. Other characteristics define whether the antibiotic destroys or not non-pathogenic (and anaerobic?) bacteria. It would be good to discuss more that aspect.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for submitting your article "Quantifying antibiotic impact on within-patient dynamics of extended-spectrum β-lactamase resistance" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Wendy Garrett as the Senior Editor.

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

Summary:

Antimicrobial resistance is recognised by the World Health Organization and others as one of the most pressing threats to global health. However, our understanding of the ecological processes that link antimicrobial use to the emergence and spread of resistance is limited. In particular (and surprisingly for such a fundamental question), a quantitative understanding of the within-host processes relating clinical antibiotic exposure to within-patient resistance abundance is lacking. For bacteria that colonize the human gut, such within-host processes are likely to play a major role in mediating the relationship between antimicrobial use and the prevalence of resistance in the wider population. This paper focuses on resistance in Enterobacteriaceae conferred by extended-spectrum β-lactamase (ESBL) production. Such ESBL-producing organisms are responsible for a high and increasing burden of disease globally. The authors aimed to determine the effects of typical antimicrobial exposures in hospitalized patients on the dynamics of resistance gene abundance and total bacterial load. The authors asked: if such effects exist, how do they vary between different treatment regimens, and what is the predicted impact of antibiotic exposure on persistence of patient colonization with ESBL-producing organisms?

Essential revisions:

The manuscript remains of high interest because of the conceptual innovation of this work related to the model.

However, reviewers still do not fully understand the dynamical model, which in our view was one of the key points of novelty of this work because it makes the link between mechanistic processes and the noisy observations. Accordingly, further revisions of this section should be completed before the paper can be accepted for publication.

It would be helpful to link the data to the notation so we know which symbols correspond to data. For clarity, even after doing this please note which are the hidden-state variables (in the dynamic model; it's listed in the observation/process noise model).

Technically, the notation a, b refers to a closed interval so infinity should not be listed in this form. We would suggest simply the phrase "improper uniform prior".

Equation 4 and text: it appears that you haven't "added them on a log scale" because the sum that you have computed goes into the mean of the normal distributions in Equation 4.

What does xt=1 = xt=1a0c1c2 have to do with any of the terms in (4) ? The text in remains inscrutable compared to the equations and notation. There is mention of the model "looping" through time points – how does this work? Please expand the text and be explicit about the likelihood at each stage, what the looping steps are, whether the terms are added to the mean or are multiplicative, and so on. Please define xi,j,g=ratio, xi,j,g=16S and all notation.

When you say that you "fitted the model with STAN" please give the posterior decomposition, and the likelihood. State which parameters are estimated. The notation of the model should be connected to Figure 4 and 5 where the model's inferences are shown.

Despite the response to review where you agree that it will be difficult to assess "abundance", you then include the following statement. It seems that the authors agree with the limitations around abundance in the response to review but then do not fully incorporate into the discussion. Please consider revising. Example: “Surprisingly, despite the relatively broad antibacterial spectrum of cefuroxime and ceftriaxone, there was no evidence that exposure to these antibiotics reduced 16S rRNA abundance (although we note that a "broad" spectrum is not defined in relation to the microbiota but to bacterial species of clinical importance).”

Essential revisions:

The reviewers brought up the previous publication, Meletiadis, J. et al., 2017, in which a good portion of the cases were reported, and where some of the same conclusions were reached.

The authors should clarify the novel contributions of this work in light of the previous manuscript. It seems that the key claimed results here are (1) variability, Figures 1, 2; (2) effect of treatment on relative resistance; (3) effect of treatment on abundance and (4) the model with fitting, and consequently the estimates re duration of carriage. Each has limitations.

Please set these results, with the novelty and limitations, in context.

This is an important point. As we mention in our Introduction, a subset of the data we use were also used in a publication by Meletiadis et al., in 2017. Their work used standard statistical methods (Fisher exact test, classification regression tree, stepwise logistic regression) on collapsed data (each patient’s time series was converted to a zero or one) to ask, is there a statistical association between ceftriaxone treatment and blaCTX-M normalised by total bacteria? Our study addresses a set of related but different questions and makes use of the full, uncollapsed data. First we are concerned with the differential effects of a range of antibiotics (not only one) both on blaCTX-M and on total bacterial load. Second, our focus is on developing a mechanistic understanding of the impact of antibiotics on the dynamics of both blaCTX-M and 16S rRNA, rather than considering static statistical associations. Third, through formulating and fitting a dynamic model to the data we aimed to estimate ecologically important parameters (strength of selection, speed of resistance waning) which are of interest in themselves. Finally, we sought to make predictions about the impact of different antibiotic treatment regimens on persistence of blaCTX-M. We have now changed the final paragraph of our Introduction to clarify how our goals are different from the previous work. We also added a new section in the Discussion that highlights the methodological novelties and advantages where we previously only stated the limitations.

The modelling is clearly new and I think it is important. In that direction, the clarity should be improved. I particularly noted the very long paragraph in subsection “Dynamic within-host model”, and the benefit that explaining the model and its assumptions (and limitations) for non-statisticians, non-modellers would likely bring.

We very much agree with this point. In response, we have rewritten the “Dynamic within-host model” section by breaking it into clearer sections, and by more clearly explaining the novelty of the model and more clearly listing the model assumptions. We also make sure to use language that is accessible for non-statisticians and by explaining some mathematical terminology.

Please clarify re interpretation of 16s, abundance as noted by reviewers.

As noted by reviewer 1, it would be interesting to validate 16S as an appropriate measure of cell density. Our study did not include an internal validation of the 16S rRNA measurements, however other studies have previously shown the ability of qPCR to quantify bacterial amounts in fecal samples (Rinttilä, 2004). We now mention this previous validation and the limitation that we do not include an internal validation in the “Quantitative PCR” section. We further added in our Discussion a comment on newer methods that quantify specifically viable bacterial cells (van Frankenhuyzen et al., 2013), which may guide future studies like ours.

Related to this, reviewer 2 notes that it may be difficult to draw conclusions about antibiotic effects on 16S rRNA given variability in DNA extraction efficiency. We agree that the variability of DNA extraction is expected to be large. Therefore our model was specifically designed to directly estimate variability that affects repeat measurements (through noise in the swab efficiency and DNA extraction) and to separate this variability from real effects on the underlying true 16S rRNA quantity. We realised that in our previous manuscript we only mentioned the swab procedure as a source of noise affecting the DNA concentration in the qPCR machine. We correct this now by being more explicit in various places for our new manuscript:

· “noise through the swab process and gene extraction”

· “observation noise (stochastic effects impacting the swab efficiency, DNA extraction or the qPCR measurements)”

· “variation in 16S rRNA would be caused mainly by the swab procedure and DNA extraction (and other steps in the protocol)”

· “error introduced by the swab procedure, DNA extraction etc.”

· “noise from qPCR measurement or from the DNA extraction process” “observation uncertainty – for instance through the swab procedure, DNA extraction, or qPCR process”

Reviewer 2 also states that one would anticipate antibiotics to cause shifts in relative species abundances as opposed to reductions in overall cell numbers. We fully share this expectation, in particular for antibiotics that target only a small subset of species present in the microbiome. Indeed, almost all of our inferred 16S effects support this hypothesis. We now take this up more explicitly in our Discussion by saying: The stable 16S abundance could be due to rapid overgrowth of non-susceptible strains which replace bacteria killed by these antibiotics (Hildebrand et al., 2019), leading to shifts in species composition rather than total species abundance.

Please improve Figures 1 and 2. While "hieroglyphics" may go a bit far, I agree that Figure 1 is visually striking but hard to see and interpret. Perhaps in Figure 1 you could group the time curves by similarity, or present an expanded view of some representative curves, with the full set in an appendix. Or something like that.

We have now changed Figures 1 and 2 to make them easier to interpret. In particular, we followed reviewer 1’s comment on Figure 2 and made it more concise by removing the additional bars of the lower part that indicated between and within patient time series variability. They were confusing as they showed a variability on the same scale as the sequence abundance, and between and within patient variability can easily be summarised and reported as numbers, which is what we do now.

For Figure 1, we have now produced a new version that contains a random selection of 50% of the study patients only. We additionally improved the previous figure with all patients by choosing a colourblind-friendly colour scheme and by increasing the plot area so that gene dynamics can be seen more clearly. We do agree with reviewer 2 (“this figure is interesting to look at the variability across the data set and would favor keeping”) that this version of the figure is a nice way to show the entire dataset. As reviewer 1 mentions, through eyeballing it appears impossible to spot waxing and waning in response to antibiotic treatment or to spot differences between treated and untreated patients, and this is exactly the purpose of this figure i.e. to show how noise is dominating the dynamics. Because we think that both are effective ways of presenting the data we attach both versions of Figure 1 for the editor to make the final decision on which one to include should the manuscript be accepted.

For ‘Figure 1 alternative’ the adjusted caption is: “Time series plots demonstrating the diverse range of dynamical patterns of blaCTX-M resistance gene abundance across patients. For visualisation purpose, a subset of 49 patients was selected uniformly at random from all 132 patients with more than one rectal swab. The x-axis scale is identical across panels, the length of one week is given for scale in the top-left corner. Timelines are ordered by length. The y-axis scale differs between panels, with the space between vertical grey lines representing a 10-fold change in the absolute blaCTX-M gene abundance (measured in copy numbers). The left-hand side shows patients who received antibiotic treatment (n=113), and the two right-hand side columns are patients without antibiotic treatment (n=19). For clarity, we show only the twelve most frequently used antibiotics in distinct colours and other antibiotics in light grey.”

Reviewer #1:

The motivation for this study and its goals are certainly important and timely. The collateral effects of antibiotic treatment on the distribution of bacteria species and strains in the enteric flora (microbiome) of humans and the incidence of antibiotic resistance genes and resistance encoding plasmids are subjects of considerable importance, epidemiologically as well as clinically. In their Introduction, the authors do a fine job of presenting this and particularly so for the extended-spectrum β-lactamase (ESBL)-producing Enterobacteriacae and the blaCTX family of resistance genes, that are the focus of their study.

This investigation certainly has the virtue of being extensive; fecal samples were taken sequentially from 133 hospitalized patients from three countries (Romania, Serbia and Italy) for a median of five samples from each patient. This virtue does have a downside; these patients were hospitalized for different reasons and treated orally of intravenously with 10 different antibiotics. The results presented suggest that some antibiotics are more likely than others to affect the distribution and abundance of these enteric bacteria and the blaCRX-M genes, their study is correlative rather than mechanistic. From the results presented, it's not clear how the different modes and purpose of treatment the within-patient dynamics ESBL resistance. From what they say in their Discussion, the authors are very much aware of this limitation of this study.

Indeed, we are aware of the limitations of our observational data regarding causal conclusions, including a potential for selection bias and confounding through patient differences that affect both the treatment decision and resistance dynamics, as well as time-varying confounding. We now add a statement on this in the Discussion.

Although their observation some antibiotics are more likely than others to affect the enteric microbiome and the frequency of blaCRX-M genes could have epidemiological and clinical implication, they don't really consider these potential implications of their results. It's also not clear how this information can be used.

Even though our Figure 5 shows the impacts of different alternative antibiotic choices on the duration of resistance carriage (a clinically interesting and relevant outcome) we refrain from drawing any clinical conclusions or recommendations from our results. This is because: 1) our study only looks at the effects on a single type of resistance; and 2) it is not clear what the individual patient or the community effects are of the observed changes to the blaCTX-M dynamics.

Presumably the choice of antibiotics, mode and frequency of administration in these hospitalized patients is based on the nature of the infection and information about the susceptibility of the target bacteria to the drugs employed for treatment. While in an ideal world, consideration should to the collateral effects of treatment; we are far from that ideal world. This is particularly so for hospitalized patients. It would be of some interest to do an analogous study for the treatment of community-acquired infections (about 90% of antibiotic use in humans), where information about the collateral effects of treatment with different antibiotics may have broader and implications that could be implemented.

This point is well taken and the obvious reason for studying hospital patients is the relative ease of following up patients. There are however existing studies on the way that explore similar questions in a community setting.

Reviewer #2:

This is an interesting and important piece of work demonstrating the correlation between certain antimicrobial exposure in patients and associated increase in abundance of blaCTX-M. The gene increase is then normalized to 16s abundance for overall internal QC for extraction efficiency and bacterial content. The work appears to be very well done and the methods and approach to analysis appear to support the conclusions of the authors. The figures are easy to interpret and additive to the manuscript (not sure if Figure 1 is critical but it is interesting to look at the variability across the data set and would favor keeping). Also the discussion includes a portion focused on relative biomass and bacterial presence in the gut of patients on antimicrobials but it is not clear they would be able to infer this information from the variability in swab collection and extraction efficiency. Would favor revising the discussion to address this issue.

One concern is the amount this manuscript overlaps with another previously published manuscript (1). The authors state in the Introduction section that they used a subset of the data but in reviewing the other manuscript it appears that there are 133 patients in this analysis and 122 from the other description. They also used very similar methodology by normalizing the blaCTX-M abundance with 16s. I think it would be helpful for the authors to further explain how this differs from prior work and is additive to the literature. The conclusions from that manuscript were almost identical as well. Would favor a re-write to address this existing work and contrast the submitted manuscript.

This is the first major point raised by the editor and we addressed this above.

Reviewer #3:

Niehus and colleagues present an original study in which they investigate the impact of different antibiotics on ESBL resistance within-host in Enterobacteriaceae carriers. The data comes from 133 patients from Romania, Serbia and Italy and consists in longitudinal series of rectal swabs (with a median of 5 swabs per patient). Swabs were analysed in order to quantify the abundance of blaCTX-M and 16SrRNA. For both indicators, the within- and between-host variabilities, and the within-host dynamics are analysed using Bayesian state-space models. Parameters associated to the selection exerted by different classes of antibiotics on the dynamics of blaCTX-M, 16SrRNA and blaCTX-M/16SrRNA are estimated.

This is a very exciting and well written article. This is, to my knowledge, the first study attempting to analyse, using hypothesis-driven mechanistic models, the dynamics of within-host resistance genes and the impact of various classes of antibiotics. I believe that the study and results presented here represent a major contribution to the field.

– It would be good for the reader to detail and clarify how each indicator should be interpreted. My understanding is for example that a decrease in blaCTX-M/16SrRNA would more or less be interpreted as a decrease in resistance rate. However, it is not clear to me how important is the contribution of Enterobacteriaceae species to the global quantity of 16SrRNA.

This is an interesting and important point, and our previous version of the manuscript was not clear enough about why we chose to model effects on blaCTX-M / 16S rRNA. As the reviewer suggests, we expect the blaCTX-M / 16S rRNA ratio to be approximately propotional to the fraction of bacterial cells harbouring the blaCTX-M gene. We then model antibiotic effects on the ratio blaCTX-M / 16S rRNA (representing effects on the relative fitness of blaCTX-M carriers) and effects on the total bacterial abundance as approximated by 16S. But our model also gives the resulting dynamics of absolute blaCTX-M abundance, which is important for predicting its extinction and persistence.

– My main comment is related to the simulation study. It would be good to provide a validation of the model before running simulations, for example by doing out of fit assessment? Is that possible? Would you have enough power?

We fully agree that model validation is important. In our study we therefore performed model comparison in which we compared our model with a model without antibiotic effects – we find that the data supports effects of antibiotics. For this we now performed a leave-one-out cross-validation approximation. We approximate because the number of patients treated with the same antibiotics is small. We have revised the manuscript so that we explain this more clearly now.

– It would be interesting to estimate whether the antibiotic classes have different delays of impact.

This is an interesting point that is on our (long) list of possible extensions to our model. However, while interesting, we consider this to be outside the focus of this study and leave this question to further work on this and similar datasets.

– I am not an expert but my understanding is that the aptitude to induce degradation of the flora may not be fully characterized by the "broad spectrum" characteristic. Other characteristics define whether the antibiotic destroys or not non-pathogenic (and anaerobic?) bacteria. It would be good to discuss more that aspect.

This a great point: the term “broad”-spectrum describes an antibiotic’s killing-spectrum in terms of known clinical pathogens. This may or may not be related to the extent of impact on the entire microbiota biomass. We now clarify this in the Discussion (first paragraph).

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Essential revisions:

The manuscript remains of high interest because of the conceptual innovation of this work related to the model.

However, reviewers still do not fully understand the dynamical model, which in our view was one of the key points of novelty of this work because it makes the link between mechanistic processes and the noisy observations. Accordingly, further revisions of this section should be completed before the paper can be accepted for publication.

It would be helpful to link the data to the notation so we know which symbols correspond to data. For clarity, even after doing this please note which are the hidden-state variables (in the dynamic model; it's listed in the observation/process noise model).

We have rewritten major parts of the dynamic model description. First, we have explained all symbols in Equation 4 within the paragraph below (previously symbols where explained in two separate parts). Second, we have included a new paragraph that lists all the symbols that represent data, and all the symbols representing estimated variables. We have been explicit about which symbols denote hidden-state variables.

Technically, the notation a, b refers to a closed interval so infinity should not be listed in this form. We would suggest simply the phrase "improper uniform prior".

This is a good point, and we now use the phrase “improper uniform prior” instead.

Equation 4 and text: it appears that you haven't "added them on a log scale" because the sum that you have computed goes into the mean of the normal distributions in Equation 4.

What does xt=1 = xt=1a0c1c2 have to do with any of the terms in (4) ? The text in remains inscrutable compared to the equations and notation. There is mention of the model "looping" through time points – how does this work? Please expand the text and be explicit about the likelihood at each stage, what the looping steps are, whether the terms are added to the mean or are multiplicative, and so on. Please define xi,j,g=ratio, xi,j,g=16S and all notation.

We agree that this part was not well explained. In response, we have substantially changed this part of our model description. The new description is explicit about which parameters are on a log-scale, and what the real scale reflects (gene copy numbers). To clarify the way antibiotic effects work, we give a numerical example, in which we used the same symbols as our Equation 4, instead of introducing new ones. We have also clarified how the sums of the definition of f(abx) work, and we explain each likelihood term in our model, and also give the full posterior distribution (see below). Finally, we have defined all notation, including xi,j,g=ratio, xi,j,g=16S.

When you say that you "fitted the model with STAN" please give the posterior decomposition, and the likelihood. State which parameters are estimated. The notation of the model should be connected to Figure 4 and 5 where the model's inferences are shown.

It was not entirely clear to us what was meant by ‘posterior decomposition’. Based on what Neiswanger et al. (https://arxiv.org/pdf/1510.04163.pdf) refer to as posterior decomposition, we now (in addition to Equation 4) give the full posterior distribution as the product of different terms (Equation 5). But in addition we also explain our three different likelihoods in terms of decomposing the total amount of variation in the data into variation attributable to different processes. We finally also added the marginal posterior distributions (together with prior distributions) of the antibiotic effect parameters (cz,b) as an additional Supplementary Figure. To address the second point, we now explicitly list all parameters of Equation 4 that are estimated. Finally, to link methods to figures, we added an explanation of how the posterior draws of cz,b were transformed before plotting them in Figure 4. We added a reference to Figure 4 at that point, and we added two references to Figure 5 into the description of the posterior predictions.

Despite the response to review where you agree that it will be difficult to assess "abundance", you then include the following statement. It seems that the authors agree with the limitations around abundance in the response to review but then do not fully incorporate into the discussion. Please consider revising. Example: “Surprisingly, despite the relatively broad antibacterial spectrum of cefuroxime and ceftriaxone, there was no evidence that exposure to these antibiotics reduced 16S rRNA abundance (although we note that a "broad" spectrum is not defined in relation to the microbiota but to bacterial species of clinical importance).”

Based on this comment, we have now decided to remove any interpretation of antibiotic effects on 16S gene abundance. We think that this also helps the discussion be more focused on the CTX-M resistance gene.

Supplementary materials
• elife-49206-repstand3.pdf info     save_alt
• elife-49206-repstand2.pdf info     save_alt
• elife-49206-repstand1.pdf info     save_alt
• elife-49206-transrepform.pdf info     save_alt

Citing articles via
https://www.researchpad.co/tools/openurl?pubtype=article&doi=10.7554/eLife.49206&title=Quantifying antibiotic impact on within-patient dynamics of extended-spectrum beta-lactamase resistance&author=Rene Niehus,Esther van Kleef,Yin Mo,Agata Turlej-Rogacka,Christine Lammens,Yehuda Carmeli,Herman Goossens,Evelina Tacconelli,Biljana Carevic,Liliana Preotescu,Surbhi Malhotra-Kumar,Ben S Cooper,Wendy S Garrett,Caroline Colijn,Caroline Colijn,Amy Mathers,Lulla Opatowski,&keyword=antibiotic resistance,within-host dynamics,resistance carriage,extended-spectrum beta-lactamase,state-space model,gut microbiota,Human,&subject=Research Article,Epidemiology and Global Health,Microbiology and Infectious Disease,