1. Epidemiology and Global Health
  2. Microbiology and Infectious Disease
Download icon

Quantifying antibiotic impact on within-patient dynamics of extended-spectrum beta-lactamase resistance

  1. Rene Niehus  Is a corresponding author
  2. Esther van Kleef
  3. Yin Mo
  4. Agata Turlej-Rogacka
  5. Christine Lammens
  6. Yehuda Carmeli
  7. Herman Goossens
  8. Evelina Tacconelli
  9. Biljana Carevic
  10. Liliana Preotescu
  11. Surbhi Malhotra-Kumar
  12. Ben S Cooper
  1. University of Oxford, United Kingdom
  2. National Institute for Public Health and theEnvironment, Netherlands
  3. University of Antwerp, Belgium
  4. Tel-Aviv University, Israel
  5. University of Tuebingen, Germany
  6. Infectious Diseases, University of Verona, Italy
  7. Clinical Centre of Serbia, Serbia
  8. Matei Balş National Institute for Infectious Diseases, Romania
Research Article
  • Cited 0
  • Views 792
  • Annotations
Cite this article as: eLife 2020;9:e49206 doi: 10.7554/eLife.49206

Abstract

Antibiotic-induced perturbation of the human gut flora is expected to play an important role in mediating the relationship between antibiotic use and the population prevalence of antibiotic resistance in bacteria, but little is known about how antibiotics affect within-host resistance dynamics. Here we develop a data-driven model of the within-host dynamics of extended-spectrum beta-lactamase (ESBL) producing Enterobacteriaceae. We use blaCTX-M (the most widespread ESBL gene family) and 16S rRNA (a proxy for bacterial load) abundance data from 833 rectal swabs from 133 ESBL-positive patients followed up in a prospective cohort study in three European hospitals. We find that cefuroxime and ceftriaxone are associated with increased blaCTX-M abundance during treatment (21% and 10% daily increase, respectively), while treatment with meropenem, piperacillin-tazobactam, and oral ciprofloxacin is associated with decreased blaCTX-M (8% daily decrease for all). The model predicts that typical antibiotic exposures can have substantial long-term effects on blaCTX-M carriage duration.

eLife digest

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.

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 with 1 supplement see all
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.

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.

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
  1. 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
Overview of antibiotic treatments showing the ten most used antibiotics in this patient cohort.

Intravenous (iv), oral (or), and intramuscular (im) route administration is given in percent of treatment days.

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 σ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 with 1 supplement see all
Variability of 16S abundance and blaCTX-M abundance within and between patient time series using a Bayesian hierarchical model.

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).

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).

Association of antibiotic use with change in relative resistance (abundance of blaCTX-M divided by abundance of 16S rRNA).

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)).

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 with 3 supplements see all
Estimated effects of different antibiotics on within-host dynamics from a multivariable model.

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.

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).

Simulated predictions of blaCTX-M carriage duration under different antibiotic treatments.

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.

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.

Though credible intervals were wide, meropenem, piperacillin-tazobactam, and oral ciprofloxacin – all common agents for treating hospital-acquired infections (Lautenbach et al., 2001; Masterton et al., 2003; Paterson, 2006) – were associated with reductions in blaCTX-M abundance. All three are broad-spectrum antibiotics with activity against ESBL producers in the absence of specific co-resistance, suggesting that this association may at least in part be a causal effect. These antibiotics were also associated with a negative effect on relative resistance (blaCTX-M divided by 16S rRNA gene abundance). This observation can be explained by a general reduction of bacterial biomass that leads the blaCTX-M abundance to drop below detection levels. In line with this, our simulations suggested that a typical course of meropenem or of piperacillin-tazobactam would reduce blaCTX-M carriage duration relative to no treatment by about 40%, and each course reduces blaCTX-M carriage duration by about 70% relative to a combined course of ceftriaxone plus amikacin. Pharmacokinetic models suggest that bacteriocidal serum concentrations of ceftriaxone persist for relatively short time periods after treatment (typically between 1 hr and 4d, Garot et al., 2011). In contrast, our model predicts that treatment effects on the gut flora can be much longer lasting (on the order of weeks). Also, a 14 day course of meropenem is predicted to reduce ESBL resistance carriage duration relative to a shortened course (5 days) by approximately 70%. While these findings suggest suppression of ESBL producing bacteria by use of carbapenem as a measure to reduce the risk of infection in high-risk patients, such a measure clearly needs to be balanced against the considerable risk of selecting for carbapenem resistance. Finally, we also found that adding oral ciprofloxacin to amoxicillin-clavulanic acid or ceftriaxone was associated with reductions in predicted median carriage duration of ESBL-producing bacteria by approximately 50%. This is in line with recent work indicating that hospital monotherapy with cephalosporins is more strongly association with later ESBL carriage relative to combination therapy (Tacconelli et al., 2020). Administration of oral fluoroquinolone to reduce faecal load of ESBL-producers in asymptomatic carriers has been used in outbreak settings with ESBL E. coli and K. pneumoniae (Paterson et al., 2001), where risk for fluoroquinolones resistance was low. However, the wide variation in rates of ciprofloxacin resistance amongst ESBL-producing Enterobacteriaceae across settings (Winokur et al., 2001) is likely to limit the generalisability of this finding. Although oral ciprofloxacin showed an association with reduced blaCTX-M abundance, intravenous ciprofloxacin showed near zero effect. Antibiotic selection for resistance due to antibiotics with different routes of administration has been previously explored in a mouse model, suggesting that, indeed, oral drug administration has stronger selective effect on resistance than intravenous administration (Zhang et al., 2013), but similar studies for humans are lacking. Delineating the relationship between the various routes of antibiotic administration and resistance selection will be important for a better understanding of advantages and disadvantages of different routes of administration.

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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:

(1) qi,j,k,gN(si,j,g,σqpcr),si,j,gN(xi,j,g,σswabg),xi,j+1,gN(si,j,g,σbiog),

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 qi,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: si,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 xi,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 σqpcr, the qPCR machine error (assumed to be the same for blaCTX-M and 16S rRNA), σswabg, the swab variation of the genetic target g, and σbiog, 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 σqpcr, σswabg, and σbiog. 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

Request a detailed protocol

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 xi,j (actual gene abundance in time point j for patient i) from the previous model, and it took the form

(2) xi,j,gN(μi,g,σwithin,g),μi,gN(μg,σbetween,g),

where mui,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 σ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 σ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:

(3) cv=esln2-1,

where sln is the estimated standard deviation of the log-scaled data.

Association of antibiotic treatment and changes in resistance

Request a detailed protocol

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

Request a detailed protocol

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 σ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:

(4) qi,j,g,kN(si,j,g,σqpcri,j,g),si,j,g=16SN(xi,j,b=16S,σswab),xi,j,b=ratio=si,j,g=CTXMsi,j,g=16S,xi,j+1,bN(xi,j,b+f(abx)i,j,b,tj+1tjσbio,b),f(abx)i,j,b=tjtj+11(ab+z=1nzcz,byz,t),

where g denotes either blaCTX-M or 16S rRNA, and b denotes either relative resistance (blaCTX-M/16S rRNA) or 16S rRNA. Then, qi,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 σqpcri,j,g, and it is estimated as described above. si,j,g is the mean of the qPCR measurements in sample j of patient i, and genetic target g, and xi,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 σ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 xi,j,b=ratio, which on a log-scale is computed as the difference (si,j,g=CTX-M-si,j,g=16S). Further, tj denotes the calendar day of sample j in patient i, and tj+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 σbio,b=ratio and σbio,b=16S, respectively, with tj+1-tj 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(abx)i,j,b, which is the change in the expected value of xi,j,b between sample j and j+1. In the definition of f(abx)i,j,b in line 5 of Equation 4, ab denotes the neutral growth or loss of g, cz,b denotes the effect of antibiotic z on b, and yz,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 (ab) and adds to this the summed effect of all antibiotics given on day t. This is computed for each calendar day from tj until a day before the subsequent sample (tj+1-1). The effects of all of these days are then summed up. Note, that xi,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 xi,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 10xi=1,j=1,b=16S=100 copy numbers and a neutral trend of ab=16S=-0.5. Suppose on the day of this sample (ti=1,j=1) two antibiotics z=1 and z=2 are given with effects cz=1,b=16S=+0.5 and cz=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, qi,j,g,k and yz,t correspond to measured data, σqpcri,j,g is computed from this data (see above). All other parameters are estimated: the hidden-state variables are si,j,g, xi,j,g=16S, and xi,j,g=ratio, the noise variables are σswab, σbio,g=16S, and σbio,g=ratio, and finally the variables describing the ecology are ag, cz,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 ab and cz,b, and σbio. All parameters (including all hidden-states) are estimated simultaneously. We can express the posterior distribution over all estimated parameters (Θ), which is conditional on the set of all data (D), and the prior over the parameter space p(Θ). For readability, here we only keep necessary subscripts:

(5) p(Θ|D)=p(s,x,a,c,σswab,σbio|q,σqpcr,y),p(Θ|D)p(Θ)i,j,g,kp(qi,j,g,k|si,j,g,σqpcri,j,g)i,jp(si,j,g=16S|xi,j,b=16S,σswab)i,j,bp(xi,j+1,b|cb,ab,σbio,b,xi,j,b,y).

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 cz,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 cz,g parameters (Figure 4—figure supplement 2). Diagnostic plots for the MCMC sampling of the cz,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.

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
    Increasing prevalence of ESBL-producing Enterobacteriaceae in Europe
    1. TM Coque
    2. F Baquero
    3. R Canton
    (2008)
    Euro Surveillance : Bulletin Europeen Sur Les Maladies Transmissibles = European Communicable Disease Bulletin 13:19044.
  10. 10
  11. 11
    Enterobacteriaceae
    1. MS Donnenberg
    (1979)
    In: M. D Bennett, editors. Principles and Practice of Infectious Diseases. Elsevier. pp. 2503–2517.
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
    An Introduction to Stochastic Processes in Physics
    1. DS Lemons
    2. P Langevin
    (2002)
    JHU Press.
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
    Mandell, Douglas, and Bennett’s Principles and Practice of Infectious Diseases
    1. JR Masci
    (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. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
    R: A Language and Environment for Statistical Computing
    1. R Development Core Team
    (2016)
    R Foundation for Statistical Computing, Vienna, Austria.
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
    Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory
    1. S Watanabe
    (2010)
    Journal of Machine Learning Research 11:3571–3594.
  52. 52
  53. 53
  54. 54
  55. 55

Decision letter

  1. Wendy S Garrett
    Senior Editor; Harvard T.H. Chan School of Public Health, United States
  2. Caroline Colijn
    Reviewing Editor; Simon Fraser University, Canada
  3. Amy Mathers
    Reviewer
  4. Lulla Opatowski
    Reviewer; Institut Pasteur, France

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).”

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

Author response

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.

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

Article and author information

Author details

  1. Rene Niehus

    University of Oxford, Oxford, United Kingdom
    Present address
    Harvard T.H. Chan School of Public Health, Harvard University, Cambridge, United States
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology
    For correspondence
    rniehus@hsph.harvard.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6751-4124
  2. Esther van Kleef

    National Institute for Public Health and theEnvironment, Bilthoven, Netherlands
    Contribution
    Conceptualization, Resources, Data curation, Formal analysis, Validation, Investigation, Visualization, Methodology, Project administration
    Competing interests
    No competing interests declared
  3. Yin Mo

    University of Oxford, Oxford, United Kingdom
    Contribution
    Conceptualization, Validation, Methodology
    Competing interests
    No competing interests declared
  4. Agata Turlej-Rogacka

    University of Antwerp, Antwerp, Belgium
    Contribution
    Resources, Validation, Investigation, Methodology
    Competing interests
    No competing interests declared
  5. Christine Lammens

    University of Antwerp, Antwerp, Belgium
    Contribution
    Resources, Data curation, Validation, Methodology
    Competing interests
    No competing interests declared
  6. Yehuda Carmeli

    Tel-Aviv University, Tel-Aviv, Israel
    Contribution
    Conceptualization, Resources, Supervision, Funding acquisition, Validation, Investigation, Methodology, Project administration
    Competing interests
    No competing interests declared
  7. Herman Goossens

    University of Antwerp, Antwerp, Belgium
    Contribution
    Conceptualization, Supervision, Funding acquisition, Validation, Investigation, Methodology, Project administration
    Competing interests
    No competing interests declared
  8. Evelina Tacconelli

    1. University of Tuebingen, Tuebingen, Germany
    2. Infectious Diseases, University of Verona, Verona, Italy
    Contribution
    Conceptualization, Supervision, Validation, Investigation, Methodology, Project administration
    Competing interests
    No competing interests declared
  9. Biljana Carevic

    Clinical Centre of Serbia, Belgrade, Serbia
    Contribution
    Resources, Investigation, Methodology
    Competing interests
    No competing interests declared
  10. Liliana Preotescu

    Matei Balş National Institute for Infectious Diseases, Bucharest, Romania
    Contribution
    Resources, Investigation
    Competing interests
    No competing interests declared
  11. Surbhi Malhotra-Kumar

    University of Antwerp, Antwerp, Belgium
    Contribution
    Conceptualization, Resources, Supervision, Funding acquisition, Validation, Investigation, Methodology, Project administration
    Competing interests
    No competing interests declared
  12. Ben S Cooper

    University of Oxford, Oxford, United Kingdom
    Contribution
    Conceptualization, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Methodology, Project administration
    Competing interests
    Reviewing Editor, eLife
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9445-7217

Funding

Medical Research Council

  • Rene Niehus
  • Ben S Cooper

Wellcome Trust (Major Overseas Programme, 106698/Z/14/Z)

  • Ben S Cooper

European Union 7th Framework Programme (R-GNOSIS Integrated project, grant agreement number 241796)

  • Rene Niehus
  • Esther van Kleef
  • Agata Turlej-Rogacka
  • Christine Lammens
  • Yehuda Carmeli
  • Herman Goossens
  • Evelina Tacconelli
  • Biljana Carevic
  • Surbhi Malhotra-Kumar
  • Ben S Cooper

Department for International Development (MR/K006924/1)

  • Rene Niehus
  • Ben S Cooper

National Medical Research Council

  • Yin Mo

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

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.

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.

Senior Editor

  1. Wendy S Garrett, Harvard T.H. Chan School of Public Health, United States

Reviewing Editor

  1. Caroline Colijn, Simon Fraser University, Canada

Reviewers

  1. Amy Mathers
  2. Lulla Opatowski, Institut Pasteur, France

Publication history

  1. Received: June 12, 2019
  2. Accepted: March 22, 2020
  3. Version of Record published: May 7, 2020 (version 1)

Copyright

© 2020, Niehus et al.

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

Metrics

  • 792
    Page views
  • 132
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Epidemiology and Global Health
    2. Medicine
    Janina Isabel Steinert et al.
    Research Article Updated

    Background:

    Since 2015, the World Health Organisation (WHO) recommends immediate initiation of antiretroviral therapy (ART) for all HIV-positive patients. Epidemiological evidence points to important health benefits of immediate ART initiation; however, the policy’s impact on the economic aspects of patients' lives remains unknown.

    Methods:

    We conducted a stepped-wedge cluster-randomised controlled trial in Eswatini to determine the causal impact of immediate ART initiation on patients’ individual- and household-level economic outcomes. Fourteen healthcare facilities were non-randomly matched into pairs and then randomly allocated to transition from the standard of care (ART eligibility at CD4 counts of <350 cells/mm3 until September 2016 and <500 cells/mm3 thereafter) to the ‘Early Initiation of ART for All’ (EAAA) intervention at one of seven timepoints. Patients, healthcare personnel, and outcome assessors remained unblinded. Data were collected via standardised paper-based surveys with HIV-positive adults who were neither pregnant nor breastfeeding. Outcomes were patients’ time use, employment status, household expenditures, and household living standards.

    Results:

    A total sample of 3019 participants were interviewed over the duration of the study. The mean number of participants approached at each facility per time step varied from 4 to 112 participants. Using mixed-effects negative binomial regressions accounting for time trends and clustering at the level of the healthcare facility, we found no significant difference between study arms for any economic outcome. Specifically, the EAAA intervention had no significant effect on non-resting time use (RR = 1.00 [CI: 0.96, 1.05, p=0.93]) or income-generating time use (RR = 0.94, [CI: 0.73,1.20, p=0.61]). Employment and household expenditures decreased slightly but not significantly in the EAAA group, with risk ratios of 0.93 [CI: 0.82, 1.04, p=0.21] and 0.92 [CI: 0.79, 1.06, p=0.26], respectively. We also found no significant treatment effect on households’ asset ownership and living standards (RR = 0.96, [CI 0.92, 1.00, p=0.253]). Lastly, there was no evidence of heterogeneity in effect estimates by patients’ sex, age, education, timing of HIV diagnosis and ART initiation.

    Conclusions:

    Our findings do not provide evidence that should discourage further investments into scaling up immediate ART for all HIV patients.

    Funding:

    Funded by the Dutch Postcode Lottery in the Netherlands, Alexander von Humboldt-Stiftung (Humboldt-Stiftung), the Embassy of the Kingdom of the Netherlands in South Africa/Mozambique, British Columbia Centre of Excellence in Canada, Doctors Without Borders (MSF USA), National Center for Advancing Translational Sciences of the National Institutes of Health and Joachim Herz Foundation.

    Clinical trial number:

    NCT02909218 and NCT03789448.

    1. Epidemiology and Global Health
    2. Microbiology and Infectious Disease
    Jon C Emery et al.
    Research Article Updated

    A key unknown for SARS-CoV-2 is how asymptomatic infections contribute to transmission. We used a transmission model with asymptomatic and presymptomatic states, calibrated to data on disease onset and test frequency from the Diamond Princess cruise ship outbreak, to quantify the contribution of asymptomatic infections to transmission. The model estimated that 74% (70–78%, 95% posterior interval) of infections proceeded asymptomatically. Despite intense testing, 53% (51–56%) of infections remained undetected, most of them asymptomatic. Asymptomatic individuals were the source for 69% (20–85%) of all infections. The data did not allow identification of the infectiousness of asymptomatic infections, however low ranges (0–25%) required a net reproduction number for individuals progressing through presymptomatic and symptomatic stages of at least 15. Asymptomatic SARS-CoV-2 infections may contribute substantially to transmission. Control measures, and models projecting their potential impact, need to look beyond the symptomatic cases if they are to understand and address ongoing transmission.