Quantifying antibiotic impact on within-patient dynamics of extended-spectrum beta-lactamase resistance
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.

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.
Summary of the study.
Number of participating hospitals | 3 (Serbia, Italy, Romania) |
---|---|
Study duration | 2 years (Jan 2011 - Dec 2012) |
Inclusion criteria | Inpatients of medical and surgical wards, adults, non-pregnant, ESBL producting Enterobacteriaceae carriers (at admission) |
Number of patients followed up | 133 (including one with a single swab taken) |
Intervals between rectal swabs | two to three days |
qPCR targets | blaCTX-M (ESBL resistance gene), 16S rRNA (total bacterial load) |
Number of different antibiotics used | 35 |
-
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.
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.
Antibiotic | Number of treated patients (total n=133) | Route |
---|---|---|
Ceftriaxone | 64 | ( 98% iv , 2% im ) |
Ciprofloxacin | 34 | ( 67% iv , 33% or ) |
Metronidazole | 20 | ( 50% iv , 50% or ) |
Cefuroxime | 13 | ( 100% iv , 0% or ) |
Vancomycin | 13 | ( 86% iv , 14% or ) |
Meropenem | 10 | ( 100% iv , 0% or ) |
Amikacin | 9 | ( 100% iv , 0% or ) |
Amoxicillin-clavulanic acid | 9 | ( 57% iv , 43% or ) |
Piperacillin-tazobactam | 7 | ( 100% iv , 0% or ) |
Imipenem | 5 | ( 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 , 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.

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

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 protocolParticipants 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 protocolSamples 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 protocolDNA 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 protocolWe 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 protocolTo 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:
where denotes a given patient, denotes a swab (one per time point), denotes a qPCR measurement (multiple repeats per swab), and denotes the genetic target, either blaCTX-M or 16S rRNA. The term , represents the measured quantity of genetic target , of the kth qPCR replicate (on a log-scale) from patient , at time point . In addition, there are two hidden-state parameter vectors: is the underlying, true sequence abundance of genetic target that a qPCR assay with 100% efficiency could (in theory) measure at time point for patient , and is the actual gene abundance of genetic target , in the patient at time point for patient , before the added noise through the swab process and gene extraction. The unobserved variables of interest are , the qPCR machine error (assumed to be the same for blaCTX-M and 16S rRNA), , the swab variation of the genetic target , and , the variation of genetic target 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 , , and . 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 protocolFor 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 (actual gene abundance in time point for patient ) from the previous model, and it took the form
where is the mean abundance of genetic target (blaCTX-M or 16S gene) and patient , around which the log-scaled measurements were assumed to be normally distributed with standard deviation , the within time series variation. The mean abundances were assumed to follow a normal distribution with a population mean µ and between-patient standard deviation . 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:
where is the estimated standard deviation of the log-scaled data.
Association of antibiotic treatment and changes in resistance
Request a detailed protocolTo 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 protocolWe 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 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:
where denotes either blaCTX-M or 16S rRNA, and denotes either relative resistance (blaCTX-M/16S rRNA) or 16S rRNA. Then, is the -th qPCR result (log-scaled) of patient , measured in the sample with index , and genetic target (blaCTX-M or 16S rRNA). The qPCR standard deviation for sample and patient is given through , and it is estimated as described above. is the mean of the qPCR measurements in sample of patient , and genetic target , and is the 16S sequence abundance that is actually present at time point in the gut of patient . The error introduced by the swab procedure, DNA extraction etc. is given through , 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 , which on a log-scale is computed as the difference (). Further, denotes the calendar day of sample in patient , and denotes the calendar day of the following sample in the same patient . The daily biological variability of the blaCTX-M/16S ratio and of 16S are given through and , respectively, with 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 , which is the change in the expected value of between sample and . In the definition of in line 5 of Equation 4, denotes the neutral growth or loss of , denotes the effect of antibiotic on , and is a boolean variable indicating whether or not antibiotic was given on day . The term inside the bracket takes, for a single calendar day , the neutral growth/loss term () and adds to this the summed effect of all antibiotics given on day . This is computed for each calendar day from until a day before the subsequent sample (). The effects of all of these days are then summed up. Note, that 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 is equivalent to multiplying the exponentiated effects on the scale of copy numbers. For example, consider that genetic target in patient has at the time of sample an abundance of copy numbers and a neutral trend of . Suppose on the day of this sample () two antibiotics and are given with effects and , then one day after this sample the genetic target abundance has an expected abundance of .
In this model, and correspond to measured data, is computed from this data (see above). All other parameters are estimated: the hidden-state variables are , , and , the noise variables are , , and , and finally the variables describing the ecology are , . 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 and , and . 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 (), and the prior over the parameter space . For readability, here we only keep necessary subscripts:
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 () 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 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 parameters (Figure 4—figure supplement 2). Diagnostic plots for the MCMC sampling of the 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 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.
Data availability
The Bayesian model code in Stan, the R code to fit the model, and the data in Rdata format as well as csv format are deposited in Dryad under https://doi.org/10.5061/dryad.8vf034t.
-
Dryad Digital RepositoryQuantifying antibiotic impact on within-host dynamics of extended-spectrum beta-lactamase resistance in hospitalized patients.https://doi.org/10.5061/dryad.8vf034t
References
-
Definitions and guidelines for research on antibiotic persistenceNature Reviews Microbiology 17:441–448.https://doi.org/10.1038/s41579-019-0196-3
-
Evolutionary epidemiology models to predict the dynamics of antibiotic resistanceEvolutionary Applications 12:365–383.https://doi.org/10.1111/eva.12753
-
Growing group of Extended-Spectrum -Lactamases: the CTX-M EnzymesAntimicrobial Agents and Chemotherapy 48:1–14.https://doi.org/10.1128/AAC.48.1.1-14.2004
-
Effects of selective digestive decontamination (SDD) on the gut resistomeJournal of Antimicrobial Chemotherapy 69:2215–2223.https://doi.org/10.1093/jac/dku092
-
Stan : a probabilistic programming languageJournal of Statistical Software 76:1.https://doi.org/10.18637/jss.v076.i01
-
Increasing prevalence of ESBL-producing Enterobacteriaceae in EuropeEuro Surveillance : Bulletin Europeen Sur Les Maladies Transmissibles = European Communicable Disease Bulletin 13:19044.
-
Within-host dynamics shape antibiotic resistance in commensal bacteriaNature Ecology & Evolution 3:440–449.https://doi.org/10.1038/s41559-018-0786-x
-
BookEnterobacteriaceaeIn: Bennett M. D, editors. Principles and Practice of Infectious Diseases. Elsevier. pp. 2503–2517.
-
The role of the intestinal tract as a reservoir and source for transmission of nosocomial pathogensClinical Infectious Diseases 39:219–226.https://doi.org/10.1086/422002
-
Microbial interactions: from networks to modelsNature Reviews Microbiology 10:538–550.https://doi.org/10.1038/nrmicro2832
-
Population pharmacokinetics of ceftriaxone in critically ill septic patients: a reappraisalBritish Journal of Clinical Pharmacology 72:758–767.https://doi.org/10.1111/j.1365-2125.2011.04005.x
-
Inference from iterative simulation using multiple sequencesStatistical Science 7:457–472.https://doi.org/10.1214/ss/1177011136
-
Application of a mathematical model to prevent in vivo amplification of antibiotic-resistant bacterial populations during therapyJournal of Clinical Investigation 112:275–285.https://doi.org/10.1172/JCI200316814
-
Epidemiological investigation of fluoroquinolone resistance in infections due to extended-spectrum beta-lactamase-producing Escherichia coli and Klebsiella pneumoniaeClinical Infectious Diseases : An Official Publication of the Infectious Diseases Society of America 33:1288–1294.https://doi.org/10.1086/322667
-
Rectal swabs are suitable for quantifying the carriage load of KPC-producing carbapenem-resistant EnterobacteriaceaeAntimicrobial Agents and Chemotherapy 57:1474–1479.https://doi.org/10.1128/AAC.01275-12
-
Antimicrobial use and antimicrobial resistance: a population perspectiveEmerging Infectious Diseases 8:347–354.https://doi.org/10.3201/eid0804.010312
-
Detection of beta-lactamase-mediated resistanceJournal of Antimicrobial Chemotherapy 48 Suppl 1:59–64.https://doi.org/10.1093/jac/48.suppl_1.59
-
The population genetics of antibiotic resistance: integrating molecular mechanisms and treatment contextsNature Reviews Genetics 11:405–414.https://doi.org/10.1038/nrg2778
-
Mandell, Douglas, and Bennett’s Principles and Practice of Infectious DiseasesClinical Infectious Diseases , Mandell, Douglas, and Bennett’s Principles and Practice of Infectious Diseases, An Official Publication of the Infectious Diseases Society of America.
-
Appropriate antimicrobial treatment in nosocomial infections-the clinical challengesJournal of Hospital Infection 55 Suppl 1:1–12.https://doi.org/10.1016/S0195-6701(03)00294-9
-
Amplification of antimicrobial resistance in gut flora of patients treated with ceftriaxoneAntimicrobial Agents and Chemotherapy 61:e00473-17.https://doi.org/10.1128/AAC.00473-17
-
Ceftriaxone: a third-generation cephalosporinDrug Intelligence & Clinical Pharmacy 19:900–906.https://doi.org/10.1177/106002808501901203
-
Cefuroxime, a Beta-Lactamase-Resistant cephalosporin with a broad spectrum of Gram-Positive and -Negative activityAntimicrobial Agents and Chemotherapy 13:657–664.https://doi.org/10.1128/AAC.13.4.657
-
Recommendation for treatment of severe infections caused by Enterobacteriaceae producing extended-spectrum beta-lactamases (ESBLs)Clinical Microbiology and Infection 6:460–463.https://doi.org/10.1046/j.1469-0691.2000.00107.x
-
Control of an outbreak of infection due to extended-spectrum beta-lactamase--producing Escherichia coli in a liver transplantation unitClinical Infectious Diseases : An Official Publication of the Infectious Diseases Society of America 33:126–128.https://doi.org/10.1086/320882
-
Resistance in gram-negative Bacteria: enterobacteriaceaeThe American Journal of Medicine 119:S20–S28.https://doi.org/10.1016/j.amjmed.2006.03.013
-
SoftwareR: A Language and Environment for Statistical ComputingR Foundation for Statistical Computing, Vienna, Austria.
-
High levels of antimicrobial coresistance among Extended-Spectrum- -Lactamase-Producing EnterobacteriaceaeAntimicrobial Agents and Chemotherapy 49:2137–2139.https://doi.org/10.1128/AAC.49.5.2137-2139.2005
-
Functional relationship between bacterial cell density and the efficacy of antibioticsJournal of Antimicrobial Chemotherapy 63:745–757.https://doi.org/10.1093/jac/dkn554
-
Practical bayesian model evaluation using leave-one-out cross-validation and WAICStatistics and Computing 27:1413–1432.https://doi.org/10.1007/s11222-016-9696-4
-
Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theoryJournal of Machine Learning Research 11:3571–3594.
-
Antibiotic administration routes significantly influence the levels of antibiotic resistance in gut MicrobiotaAntimicrobial Agents and Chemotherapy 57:3659–3666.https://doi.org/10.1128/AAC.00670-13
Decision letter
-
Wendy S GarrettSenior Editor; Harvard T.H. Chan School of Public Health, United States
-
Caroline ColijnReviewing Editor; Simon Fraser University, Canada
-
Amy MathersReviewer
-
Lulla OpatowskiReviewer; 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.sa1Author 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.sa2Article and author information
Author details
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
- Wendy S Garrett, Harvard T.H. Chan School of Public Health, United States
Reviewing Editor
- Caroline Colijn, Simon Fraser University, Canada
Reviewers
- Amy Mathers
- Lulla Opatowski, Institut Pasteur, France
Publication history
- Received: June 12, 2019
- Accepted: March 22, 2020
- 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
-
- 1,483
- Page views
-
- 227
- Downloads
-
- 11
- Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Epidemiology and Global Health
Background:
Short-term forecasts of infectious disease burden can contribute to situational awareness and aid capacity planning. Based on best practice in other fields and recent insights in infectious disease epidemiology, one can maximise the predictive performance of such forecasts if multiple models are combined into an ensemble. Here, we report on the performance of ensembles in predicting COVID-19 cases and deaths across Europe between 08 March 2021 and 07 March 2022.
Methods:
We used open-source tools to develop a public European COVID-19 Forecast Hub. We invited groups globally to contribute weekly forecasts for COVID-19 cases and deaths reported by a standardised source for 32 countries over the next 1–4 weeks. Teams submitted forecasts from March 2021 using standardised quantiles of the predictive distribution. Each week we created an ensemble forecast, where each predictive quantile was calculated as the equally-weighted average (initially the mean and then from 26th July the median) of all individual models’ predictive quantiles. We measured the performance of each model using the relative Weighted Interval Score (WIS), comparing models’ forecast accuracy relative to all other models. We retrospectively explored alternative methods for ensemble forecasts, including weighted averages based on models’ past predictive performance.
Results:
Over 52 weeks, we collected forecasts from 48 unique models. We evaluated 29 models’ forecast scores in comparison to the ensemble model. We found a weekly ensemble had a consistently strong performance across countries over time. Across all horizons and locations, the ensemble performed better on relative WIS than 83% of participating models’ forecasts of incident cases (with a total N=886 predictions from 23 unique models), and 91% of participating models’ forecasts of deaths (N=763 predictions from 20 models). Across a 1–4 week time horizon, ensemble performance declined with longer forecast periods when forecasting cases, but remained stable over 4 weeks for incident death forecasts. In every forecast across 32 countries, the ensemble outperformed most contributing models when forecasting either cases or deaths, frequently outperforming all of its individual component models. Among several choices of ensemble methods we found that the most influential and best choice was to use a median average of models instead of using the mean, regardless of methods of weighting component forecast models.
Conclusions:
Our results support the use of combining forecasts from individual models into an ensemble in order to improve predictive performance across epidemiological targets and populations during infectious disease epidemics. Our findings further suggest that median ensemble methods yield better predictive performance more than ones based on means. Our findings also highlight that forecast consumers should place more weight on incident death forecasts than incident case forecasts at forecast horizons greater than 2 weeks.
Funding:
AA, BH, BL, LWa, MMa, PP, SV funded by National Institutes of Health (NIH) Grant 1R01GM109718, NSF BIG DATA Grant IIS-1633028, NSF Grant No.: OAC-1916805, NSF Expeditions in Computing Grant CCF-1918656, CCF-1917819, NSF RAPID CNS-2028004, NSF RAPID OAC-2027541, US Centers for Disease Control and Prevention 75D30119C05935, a grant from Google, University of Virginia Strategic Investment Fund award number SIF160, Defense Threat Reduction Agency (DTRA) under Contract No. HDTRA1-19-D-0007, and respectively Virginia Dept of Health Grant VDH-21-501-0141, VDH-21-501-0143, VDH-21-501-0147, VDH-21-501-0145, VDH-21-501-0146, VDH-21-501-0142, VDH-21-501-0148. AF, AMa, GL funded by SMIGE - Modelli statistici inferenziali per governare l'epidemia, FISR 2020-Covid-19 I Fase, FISR2020IP-00156, Codice Progetto: PRJ-0695. AM, BK, FD, FR, JK, JN, JZ, KN, MG, MR, MS, RB funded by Ministry of Science and Higher Education of Poland with grant 28/WFSN/2021 to the University of Warsaw. BRe, CPe, JLAz funded by Ministerio de Sanidad/ISCIII. BT, PG funded by PERISCOPE European H2020 project, contract number 101016233. CP, DL, EA, MC, SA funded by European Commission - Directorate-General for Communications Networks, Content and Technology through the contract LC-01485746, and Ministerio de Ciencia, Innovacion y Universidades and FEDER, with the project PGC2018-095456-B-I00. DE., MGu funded by Spanish Ministry of Health / REACT-UE (FEDER). DO, GF, IMi, LC funded by Laboratory Directed Research and Development program of Los Alamos National Laboratory (LANL) under project number 20200700ER. DS, ELR, GG, NGR, NW, YW funded by National Institutes of General Medical Sciences (R35GM119582; the content is solely the responsibility of the authors and does not necessarily represent the official views of NIGMS or the National Institutes of Health). FB, FP funded by InPresa, Lombardy Region, Italy. HG, KS funded by European Centre for Disease Prevention and Control. IV funded by Agencia de Qualitat i Avaluacio Sanitaries de Catalunya (AQuAS) through contract 2021-021OE. JDe, SMo, VP funded by Netzwerk Universitatsmedizin (NUM) project egePan (01KX2021). JPB, SH, TH funded by Federal Ministry of Education and Research (BMBF; grant 05M18SIA). KH, MSc, YKh funded by Project SaxoCOV, funded by the German Free State of Saxony. Presentation of data, model results and simulations also funded by the NFDI4Health Task Force COVID-19 (https://www.nfdi4health.de/task-force-covid-19-2) within the framework of a DFG-project (LO-342/17-1). LP, VE funded by Mathematical and Statistical modelling project (MUNI/A/1615/2020), Online platform for real-time monitoring, analysis and management of epidemic situations (MUNI/11/02202001/2020); VE also supported by RECETOX research infrastructure (Ministry of Education, Youth and Sports of the Czech Republic: LM2018121), the CETOCOEN EXCELLENCE (CZ.02.1.01/0.0/0.0/17-043/0009632), RECETOX RI project (CZ.02.1.01/0.0/0.0/16-013/0001761). NIB funded by Health Protection Research Unit (grant code NIHR200908). SAb, SF funded by Wellcome Trust (210758/Z/18/Z).
-
- Epidemiology and Global Health
Background:
Affectionate touch, which is vital for mental and physical health, was restricted during the Covid-19 pandemic. This study investigated the association between momentary affectionate touch and subjective well-being, as well as salivary oxytocin and cortisol in everyday life during the pandemic.
Methods:
In the first step, we measured anxiety and depression symptoms, loneliness and attitudes toward social touch in a large cross-sectional online survey (N = 1050). From this sample, N = 247 participants completed ecological momentary assessments over 2 days with six daily assessments by answering smartphone-based questions on affectionate touch and momentary mental state, and providing concomitant saliva samples for cortisol and oxytocin assessment.
Results:
Multilevel models showed that on a within-person level, affectionate touch was associated with decreased self-reported anxiety, general burden, stress, and increased oxytocin levels. On a between-person level, affectionate touch was associated with decreased cortisol levels and higher happiness. Moreover, individuals with a positive attitude toward social touch experiencing loneliness reported more mental health problems.
Conclusions:
Our results suggest that affectionate touch is linked to higher endogenous oxytocin in times of pandemic and lockdown and might buffer stress on a subjective and hormonal level. These findings might have implications for preventing mental burden during social contact restrictions.
Funding:
The study was funded by the German Research Foundation, the German Psychological Society, and German Academic Exchange Service.