Modeling the bacterial dynamics in the gut microbiota following an antibiotic-induced perturbation.

Recent studies have highlighted the importance of ecological interactions in dysbiosis of gut microbiota, but few focused on their role in antibiotic-induced perturbations. We used the data from the CEREMI trial in which 22 healthy volunteers received a 3-day course of ceftriaxone or cefotaxime antibiotics. Fecal samples were analyzed by 16S rRNA gene profiling, and the total bacterial counts were determined in each sample by flux cytometry. As the gut exposure to antibiotics could not be experimentally measured despite a marked impact on the gut microbiota, it was reconstructed using the counts of susceptible Escherichia coli. The dynamics of absolute counts of bacterial families were analyzed using a generalized Lotka-Volterra equations and nonlinear mixed effect modeling. Bacterial interactions were studied using a stepwise approach. Two negative and three positive interactions were identified. Introducing bacterial interactions in the modeling approach better fitted the data, and provided different estimates of antibiotic effects on each bacterial family than a simple model without interaction. The time to return to 95% of the baseline counts was significantly longer in ceftriaxone-treated individuals than in cefotaxime-treated subjects for two bacterial families: Akkermansiaceae (median [range]: 11.3 days [0; 180.0] vs. 4.2 days [0; 25.6], p = 0.027) and Tannerellaceae (13.7 days [6.1; 180.0] vs. 6.2 days [5.4; 17.3], p = 0.003). Taking bacterial interaction as well as individual antibiotic exposure profile into account improves the analysis of antibiotic-induced dysbiosis.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

Whereas the composition of the bacterial gut microbiota is highly variable among individuals, ecological interactions contribute to the altered bacterial dynamics following antibiotics administration.

WHAT QUESTION DID THIS STUDY ADDRESS?

We investigated whether combining the Lotka‐Volterra equations with nonlinear mixed effect approach enables to better quantify antibiotics‐induced perturbation of bacterial taxa.

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

We developed an analytical framework to study the bacterial dynamics following an external perturbation by antibiotics. Bacterial interactions need to be accounted for when studying the effect of antibiotics on individual bacterial taxa.

HOW MIGHT THIS CHANGE DRUG DISCOVERY, DEVELOPMENT, AND/OR THERAPEUTICS?

When the association between bacterial composition in gut microbiota and antibiotic exposure is investigated, especially regarding the emergence of antibiotic resistance, our approach may be used to quantify both an antibiotic effect and identify relevant bacterial interactions.

INTRODUCTION

Molecular approaches to genome analysis increased our knowledge of the composition of microbiota and its role in host health. It is now recognized that the intestinal microbiota and its host function in mutualism, each partner contributing to the balance and health of the other. However, what precisely constitutes a healthy microbiota is still poorly defined. Even if there is a conserved balance between healthy subjects in the taxonomic groups of bacteria constituting the microbiota, its composition in terms of species, subspecies, or clades might vary considerably with the origin of the subjects, their genetic background, or their diet and life habits. To study interactions between the different members of the bacterial ecosystem, the variability in the composition of gut microbiota in healthy subjects needs to be considered.

A prolific literature relates the composition of the gut microbiota to a large number of acute and chronic disease states, such as inflammatory bowel disease, diabetes mellitus, and or mental disorders. , , , The role of many disruptive factors has been highlighted, antibiotics being one of the most important. , , , For decades, they have been recognized to profoundly disturb the ecological equilibrium of the gut microbiota, a phenomenon called dysbiosis. Generally, such a shift is associated with modifications in the bacterial composition of the ecosystem and a reduction in alpha biodiversity, which are likely to allow the establishment of a new balance within the bacterial ecosystem and the colonization of the ecosystem by resistant bacterial. , , , , Accordingly, understanding the altered bacterial dynamics following antibiotics administration is critical to characterize their harmful consequences. The most natural mathematical framework for studying bacterial dynamics is the generalized Lotka‐Volterra model, which can account for both interactions and the effect of exogeneous factors. , , Nevertheless, the use of this model requires to have absolute counts of bacterial species, and not only relative abundance that are typically provided by sequencing methods. Recently, an approach for the flow cytometric enumeration of microbial cells in a fecal sample has been proposed. Combined with the metagenomic analysis of bacterial communities, this approach allows to determine the absolute counts of each bacterial taxa for the study of the interplay between bacterial taxonomic groups and of the impact of exogeneous perturbations of the ecosystem.

Here, we applied the generalized Lotka‐Volterra equations on absolute counts of bacterial families to study the dynamics of the bacterial population after an antibiotic perturbation, considering both the gut exposure to antibiotics and the ecological interactions within the bacterial ecosystem.

METHODS AND MATERIALS

Study population and sample collection

We used the data of the CEREMI trial (ClinicalTrials.gov identifier NCT02659033), a prospective, randomized open‐label clinical trial conducted at the Clinical Investigation Center of the Bichat‐Claude Bernard Hospital (Paris, France) from March 2016 to August 2017. The trial was sponsored by Assistance Publique—Hȏpitaux de Paris and approved by French Health Authorities and by the Independent Ethics Committee Ȋle‐de‐France 1. All procedures were conducted in compliance with good clinical practice and the Declaration of Helsinki. Full details of the trial have been reported elsewhere. Healthy volunteers aged between 18 and 65 years old without exposure to antibiotics in the preceding 3 months nor history of hospitalization in the last 6 months were prospectively included after obtention of their informed consent. A total of 22 healthy volunteers were randomly assigned (1:1 ratio) to receive for 3 days either 1 g of ceftriaxone once a day (n = 11) or 1 g of cefotaxime three times a day. Antibiotic treatments were administered as 30‐min intravenous infusions. A total of 13 fecal samples were collected: before the beginning of treatment at days −15, −7, and −1; every day during treatment at days 1, 2, and 3, and after the end of treatment at days 4, 7, 10, 15, 30, 90, and 180.

Determination of bacterial counts in feces

Total counts of Gram‐negative enteric bacilli (GNB) were determined on all samples by plating serial dilutions of fresh stools on Drigalski agar (bioMérieux, Marcy l'Etoile, France). Third generation resistant (3GC) GNB were detected and counted on ChromID ESBL agar (bioMérieux) and biplate ESBL agar (AES Chemunex; 37°C for 48 h). Full details are available in the previous results. For each sample, 3GC‐resistant GNB counts were designed as the smallest value obtained from these two agar media.

Targeted metagenomic analysis of the intestinal microbiota

Samples obtained at days −15, −1, 4, 7, 10, 30, and 180 were analyzed by 16S rRNA gene profiling. Microbial DNA was extracted using QIAamp DNA stool Kit (Qiagen). The V4 region of the 16S rRNA gene was then amplified using primers selected from previous research. The polymerase chain reaction (PCR) amplicons (292 bp) were sequenced on an Illumina MiSeq platform according to the manufacturer’s specifications (Illumina). Demultiplexing and quality filtering were performed using mothur. No mismatch was allowed in the primers’ sequences in both forward and reverse primers, and quality‐filtering was performed by truncating bases at the 3′ end with Phred quality score less than 25. Paired‐end read assembly was then performed with 100% overlap identity.

Following these pre‐processing steps, chimera sequences were detected and eliminated using Vsearch. Then, clustering of similar sequences (97% identity threshold for an affiliation at the genus level) was performed through an open‐reference Operational Taxonomic Unit (OTU) picking process and complete‐linkage method, finally creating groups of OTUs. A final OTU cleaning step corresponding to the elimination of singletons and doubletons was performed. OTUs were then aggregated at the family level, and the relative abundance of each bacterial family was determined.

Determination of total bacterial counts by flux cytometry

Aliquots containing 200 mg of feces were diluted 20,000 times in a physiological solution (8.5 g/L NaCl). Samples were filtered for removing debris from the fecal solutions using a sterile syringe filter (pore size 5 μm; Sartorius Stedim Biotech GmbH). Then, 1 ml of the microbial cell suspension obtained was stained with 1 μl SYBR Green I (1:100 dilution in dimethylsulfoxide; shaded 15 min incubation at 37°C; 10,000 concentrate; Thermo Fisher Scientific). The flow cytometry analysis of the microbial cells present in the suspension was performed using a BD LSRFortessa flow cytometer (BD Biosciences). Fluorescence events were monitored using the FITC filter 505LP 530/30 nm and perCP filter 635LP 695/40 nm optical detectors. Forward and sideways‐scattered light was also collected with a threshold value of 200. The BD FACSDiva software was used to gate and separate the microbial fluorescence events on the FL1–FL3 density plot from the fecal sample background. The gated fluorescence events were evaluated on the forward–sideways density plot, to exclude remaining background events and to obtain an accurate microbial cell count with FlowJo software for a secondary analysis.

Absolute counts of each bacterial family were determined in each sample by multiplying the relative abundance of the family obtained by 16 rRNA gene profiling by the total number of bacterial obtained by flux cytometry. Absolute counts were log10‐transformed in all analyses. If a family was not detected in a sample, its absolute count was imputed to the half of the minimal absolute count observed in the sample.

Selection of families considered in the ecological model

The model of bacterial family’s dynamics was constructed including bacterial interactions. To limit the number of bacterial families modeled, we arbitrarily restricted this analysis to bacterial families whose relative abundance was greater than 1% in one sample (whatever the sampling time) for at least eight healthy volunteers, and which were significantly impacted by the antibiotic treatment. To that end, we computed for each healthy volunteer the area under the curve of the change from baseline of the absolute counts for each bacterial family between baseline and day 30 by the trapezoid rule, and compared it to zero using a bilateral Wilcoxon test, with a bilateral type‐I error of 0.05. We defined baseline as day 0, and the baseline sample was defined as the sample obtained at day −1. If this sample at day −1 was not available, the sample obtained at day −7 was considered as the baseline, or the one obtained at day −15 if this latter was also missing. Bacterial families that were not significantly impacted by the antibiotic treatment were used as covariate in the modeling analysis.

Reconstruction of the gut antibiotic exposure using Gram‐negative bacilli kinetics

As the gut exposure to antibiotics could not be experimentally measured due to the putative action of fecal betalactamases, we developed a simple model using the counts of susceptible Gram‐negative bacilli to reconstruct the antibiotic gut exposure. As the total counts of Gram‐positive bacteria were not available, we assumed that the interaction between Gram‐negative and Gram‐positive bacteria were negligible. We assumed that the count of susceptible GNB, , was at steady‐state prior to treatment initiation, and assumed that antibiotics act on the rate of bacterial death in a linear fashion, with a rate noted : where is the net growth rate of GNB, is the carrying capacity of the gut, and is the gut antibiotic exposure. As , cannot be measured in feces, it was here modeled as a dimensionless quantity ranging between 0 (no exposure) and 1 (maximal exposure). Following intravenous administration, increased with a first‐order absorption rate, , and maximal antibiotic gut concentration value was achieved at the end of the last antibiotic infusion, (Equation 2). Afterward, gut exposure declined with a first‐order term elimination, (Equation 3) such that:

Lotka‐Volterra model structure for bacterial families

The effect of antibiotics on reducing the counts of each bacterial family was modeled using a Generalized Lotka‐Volterra equations, which is based on a logistic bacterial growth model and takes into account the ecological interactions between the bacterial families in the feces. In this model, the evolution of the absolute count of the k bacterial family over time can be written as Equation 4 below, where is the intrinsic growth rate of the k bacterial family, is the interaction term of bacterial family k′ on bacterial family k, is a dimensionless quantity denoting the antibiotic gut exposure, and is the killing rate of the antibiotic on the k bacterial family. We assumed that for each bacterial family, both antibiotics had the same killing rate . , Assuming that steady‐state is achieved before the beginning of treatment, we parameterized the model such that the intrinsic growth rate of the k bacterial family is equal to the sum of the interactions on k family as in Equation 5 below, where is the absolute count of the k′ bacterial family at baseline.

Statistical model

Nonlinear mixed effects models were used to analyze each bacterial count over time, which was defined as follows. Let denote the vector of observations for all responses and the vector of observations for the family (e.g., corresponds to Akkermansiaceae, to Bifidobacteriaceae, etc.) for the individual . Let denote the global structural model characterizing all families, based on a system of ordinary differential equations, similar for all individuals. Then one can define , where is the component of the global model describing the family, is the vector of individual parameters, is the vector of sampling times, and is the vector of residual errors for the family in individual . Each individual parameter can be decomposed as a fixed effect , which represents the mean value of the parameter in the population, and a random effect where accounts for the interindividual variability. Assuming an exponential random effect model, the individual parameters are modeled as: . Regarding the residual error model, we assumed that constant error model where is a ‐diagonal matrix with elements equal to

Parameter estimation and model building

Estimation of population parameters was performed using the stochastic approximation expectation maximization algorithm (SAEM), implemented in MONOLIX 2018R2 (Lixoft, Orsay, France, www.lixoft.eu), a software devoted to parameter estimation by maximum likelihood in nonlinear mixed effect models. This algorithm handles data below the limit of quantification.

Reconstruction of the antibiotic gut exposure

We first reconstructed the antibiotic gut exposure by fitting the counts of susceptible Gram‐negative bacilli (see Equations 1 and 2). In this model, interindividual variability was introduced on the gut elimination rate of antibiotics , on the killing rate , and on the counts of susceptible GNB at baseline . In order to take account of the different half‐life between ceftriaxone and cefotaxime, we tested a group effect on the elimination rate, , using the likelihood ratio test. Empirical Bayes estimates of parameters from the final model of antibiotic gut exposure were used to predict the individual profiles of gut antibiotic exposure during and after treatment, which were used in the dynamic model of bacterial families.

Ecological model building

We used the empirical Bayes estimates of individual parameters for the antibiotic gut exposure model to reconstruct antibiotic gut exposure profile for each subject, and its impact on the evolution of bacterial ecology. Each bacterial family dynamics was first analyzed separately, without interaction terms. We then added interaction terms between bacterial families using a stepwise approach with a p value for the likelihood‐ratio test of 0.05 for forward selection and a p value of 0.001 for backward deletion. The likelihood was estimated by linearization.

Model evaluation

Model evaluation was conducted by several goodness‐of‐fit plots: individual fits, plots of predictions versus observations, distribution of the individual weighted residuals (IWRESs), and normalized prediction distribution errors (NPDEs) versus time and versus model predictions, as well as the visual predictive check (VPC). NPDE and VPC were generated using 500 Monte Carlo simulations.

Measures of antibiotics impact on the bacterial families

We used the empirical Bayes estimates of the individual parameters to predict the evolution of each modeled bacterial family over time, for each healthy volunteer, and computed the following parameters: maximal loss of each bacterial family after the beginning of the antibiotic treatment, and time at which each bacterial family returned to 95% of its baseline counts. Differences between antibiotic treatment groups were compared using the bilateral Wilcoxon test.

RESULTS

Gut antibiotic exposure

As fecal concentrations of antibiotics were largely undetectable among participants, we reconstructed the gut antibiotic exposure using the fecal counts of susceptible GNB. A total of 286 measurements were available. Estimated population parameters are presented in the Table S1. This model fitted the counts of GNB well, and the goodness‐of‐fit plots did not reveal a significant model misspecification (Figure S1), thereby providing good confidence in the prediction of the individual gut exposure. The elimination rate of ceftriaxone was significantly lower than that of cefotaxime (0.21/day, relative standard error [RSE], 23.9% vs. 1.01/day, RSE 41.5%, p = 0.0005), resulting in empirical Bayes estimates showing that the time to eliminate 90% of exposure was 12.7 days (range, 5.3–26.1) in ceftriaxone and 4.9 days (range, 4.0–7.3) in cefotaxime, respectively (p = 0.00012). In GNB counts, the median time to return to 95% of baseline was 3.3 days longer in ceftriaxone than in cefotaxime (7.9 days, range, 5.4–14.4 vs. 4.6 days, range, 3.5; 8.0, respectively, p = 0.0006).

Selection of families for modeling

Of the 116 bacterial families identified in the fecal samples, 15 had a prevalence above 1% in at least one sample in eight subjects and were accordingly further studied. The evolution of their relative and absolute counts over time are presented in the Figure 1 and Figure S2. Overall, Lachnospiraceae, Bacteroidaceae, Ruminococcaceae, and Prevotellaceae were the most abundant families, accounting for 36.8% of bacterial taxa in the samples collected before antibiotic administration. Of these 15 families, seven (Peptostreptococcaceae, Bifidobacteriaceae, Akkermansiaceae, Tannerellaceae, Eggerthellaceae, Veillonellaceae, and Lachnospiraceae) were significantly disrupted by the antibiotic treatment, and were considered as dependent variables in the modeling analysis (Table S2 and Figure S3).

FIGURE 1

Evolution of the absolute counts of the most prevalent taxa over time, in the 22 healthy volunteers of the CEREMI trial. The evolution of the absolute counts of the seven bacterial families which were significantly impacted by the antibiotic treatment

Ecological model of bacterial interactions and antibiotics effect

Using a generalized Lotka‐Volterra model, we built an ecological model to study the interplay between bacterial families and the impact of antibiotics on their dynamics. Of the 21 pairwise interactions between bacterial families studied, five were retained in the final model (Table 1 and Figure 2). Full results of the model selection are summarized in Table S3. Goodness‐of‐fit plots for the final model were satisfactory, with no significant trends in the residuals (Figure 3 and Figures S4 and S5).

TABLE 1

Estimated population parameters and their RSEs for the final ecological model of the seven bacterial families significantly impacted by the antibiotic treatment

	Lac.	Vei.	Bif.	Tan.	Akk.	Pep.	Egg.
Population parameters
N_0(log10 CFU/g)	10.50 (0.43)	9.86 (0.64)	9.54 (0.68)	9.16 (1.34)	8.97 (1.42)	8.97 (0.83)	9.08 (0.99)
μ_kk(×10^−9 log10 CFU/g/day)	3.55 (37.5)	0.27 (21.1)	0.18 (21.1)	0.06 (66.6)	5.31 (32.4)	0.60 (24.0)	18.4 (38.2)
γ(day−1)	26.2 (26.8)	1.42 (26.8)	1.07 (16.8)	2.63 (11.2)	1.88 (34.9)	1.05 (14.8)	1.32 (263)
Standard deviations of the exponential random effects
ω_N_0	0.19 (18.0)	0.26 (17.8)	0.15 (21.4)	0.43 (21.2)	0.57 (16.1)	0.27 (18.6)	0.39 (16.8)
ω_μ_kk	0.71 (29.5)	0.31 (84.4)	0.22 (35.0)	0.78 (122)	0.46 (44.8)	0.55 (32.9)	0.56 (37.8)
ω_γ	0.85 (23.6)	0 fix	0 fix	0 fix	0 fix	0.35 (36.0)	0 fix
Residual errors
σ_inter,N (log10 CFU/g)	0.23 (6.45)	0.28 (7.01)	0.55 (6.20)	0.69 (6.95)	0.34 (7.53)	0.44 (6.66)	0.35 (6.85)

Abbreviations: , interaction term within the bacterial families (intrinsic death rate); , additive residual error; Akk., Akkermansiaceae; Bif., Bifidobacteriaceae; Chr., Christensenellaceae; Egg., Eggerthellaceae; Lac., Lachnospiraceae; , absolute count of the bacterial family at baseline; Pep., Peptostreptococcaceae; Rik., Rikenellaceae; RSE, relative standard error; Tan., Tannerellaceae; Vei., Veillonellaceae; γ, antibiotic killing rate; , standard deviation of the exponential random effects of each fixed effect parameters.

FIGURE 2

Interactions network between the seven bacterial families in the final model. Red lines indicate inhibition of bacterial growth, whereas green arrows indicate stimulation of bacterial growth

FIGURE 3

Predicted individual profiles of absolute counts of the seven modeled bacterial families, in the 22 healthy volunteers of the CEREMI trial. Colored solid lines correspond to absolute counts of each bacterial family. Dots are observed absolute counts for each bacterial family

The analysis of the interaction terms of bacterial families showed that Peptostreptococcaeae and Lachnospiraceae exerted a stimulating activity on Eggerthellaceae growth, whereas Bifidobacteriaceae inhibited it. The growth of Akkermansiaceae, was in turn stimulated by Eggerthellaceae. Bifidobacteriaceae also inhibited the growth of Tannerellaceae. The estimated effect of antibiotics on these bacterial families is presented in Table 1. The highest effect was observed on Lachnospiraceae, estimated to 26.2/day (RSE, 26.8%), and minimal for Bifidobacteriaceae (1.1/day, RSE, 16.8%) and Peptostreptococcaeae (1.1/day, RSE, 14.8%).

In order to assess the impact of these bacterial interaction, we removed all interaction terms from the model. The Bayesian Information Criteria was deteriorated by 65.6 points, showing that these interactions provided a better fit to the data. Ignoring these interactions also led to substantial changes in the estimation of the antibiotic effects, toward either a smaller (Lachnospiraceae, Veillonellaceae, and Tannerellaceae) or a larger effect (Bifidobacteriaceae, Akkermansiaceae, Peptostreptococcaeae, and Eggerthellaceae; Figure 4 and Table S4).

FIGURE 4

Ratios of individual estimates of antibiotic killing rates effects on the seven modeled bacterial families in the model without interactions between bacterial families and in the final model with interactions, in the 22 healthy volunteers of the CEREMI trial. Whiskers represent 10th and 90th percentiles

Impact of antibiotic exposure on bacterial kinetics

We used the individual parameters estimates of the model to predict the kinetics of the seven bacterial families. The highest median maximal loss from baseline was observed for Tannerellaceae, with a value of 1.7 CFU/g (range, 1.2–2.6), whereas the lowest median maximal loss was observed for the Lachnospiraceae with a value of 0.1 CFU/g (range, 0.04–3.1). Following the assumption used of a similar effect for both antibiotics, no significant difference was observed between the two antibiotics in the magnitude of bacterial loss (Table 2).

TABLE 2

Measures of the impact of ceftriaxone and cefotaxime on the bacterial counts in the gut microbiota, in the 22 healthy volunteers of the CEREMI trial

Family	All groups (n = 22)	Ceftriaxone (n = 11)	Cefotaxime (n = 11)	p Value
Maximal reduction from baseline (log10 CFU/g)
Lachnospiraceae	0.1 [0.04; 3.1]	0.1 [0.04; 3.1]	0.1 [0.06; 0.9]	0.67
Veillonellaceae	0.5 [0.1; 3.0]	0.7 [0.1; 2.9]	0.5 [0.1; 3.0]	0.72
Bifidobacteriaceae	0.7 [0.5; 1.2]	0.7 [0.5; 1.2]	0.7 [0.6; 0.9]	0.97
Tannerellaceae	1.7 [1.2; 2.6]	1.9 [1.2; 2.6]	1.6 [1.4; 2.4]	0.34
Akkermansiaceae	1.3 [0.3; 6.2]	1.3 [0.3; 5.0]	1.4 [0.4; 6.2]	0.97
Peptostreptococcaceae	0.8 [0.4; 1.9]	1.0 [0.4; 1.9]	0.7 [0.4; 1.9]	0.41
Eggerthellaceae	0.4 [0.02; 4.4]	0.3 [0.02; 4.4]	0.9 [0.04; 4.0]	0.82
Time to return to 95% of baseline counts (days)
Lachnospiraceae	0 [0; 4.7]	0 [0; 4.7]	0 [0; 2.9]	0.48
Veillonellaceae	1.9 [0;15.5]	4.9 [0; 15.5]	0 [0; 7.9]	0.31
Bifidobacteriaceae	6.2 [0; 17.7]	7.1 [0; 17.7]	5.6 [4.2; 9.5]	0.071
Tannerellaceae	8.5 [5.4; 180.0]	13.7 [6.1; 180.0]	6.2 [5.4; 17.3]	0.0030
Akkermansiaceae	5.4 [0; 180.0]	11.3 [0; 180.0]	4.2 [0; 25.6]	0.027
Peptostreptococcaceae	9.3 [0; 28.6]	10.3 [0; 28.6]	5.9 [0; 20.8]	0.14
Eggerthellaceae	0 [0; 28.7]	0 [0; 28.7]	3.5 [0; 24.2]	0.77

Note: Results are presented as median [min; max]. p‐Value refer to the comparison of the two antibiotic treatment groups, using the Wilcoxon test.

The time needed to return to 95% of baseline counts was variable among bacterial families, with a delay that could be as large as 9.3 days for Peptostreptococcaceae; there was also an impact of the type of treatment, with ceftriaxone leading to a more sustained perturbation for Akkermansiaceae and Tannerellaceae than cefotaxime (Table 2).

DISCUSSION

In recent years, several attempts have been made to quantify the antibiotic‐induced dysbiosis of the gut microbiota. , , The significant role of the ecological interactions in the gut microbiota following an external stimulus has been highlighted. , Yet, only few studies evaluated the impact of antibiotics on the bacterial ecosystem of the gut microbiota considering both gut antibiotic exposure and ecological interactions. In addition, these studies focused on the evolution of bacterial counts using relative abundances, assuming therefore that the total number of bacteria in the gut remained constant after treatment initiation.

We used a generalized Lotka‐Volterra model to decipher the effect of antibiotics on the evolution of bacterial groups at the family taxonomic level and their interactions. We found that ignoring these interactions could lead to spurious estimates of antibiotic effects of bacterial families, that could be either under‐ or overestimated. Although interactions at the family level only provide a crude estimate of the genuine biological processes among bacterial species, this suggests that bacterial interactions need to be taken into account when studying the impact of antibiotics on the gut microbiota.

We previously used next‐generation sequencing methods to study the relationship between antibiotic gut exposure and evolution of synthetic markers of bacterial diversity within the gut microbiome. Given the very high number of taxa that can be identified in a single sample and as no general rule is used regarding bacterial taxa that should be included in the analyses, we arbitrarily restricted the modeling analysis to the most abundant and prevalent taxa, and to those which were significantly impacted by antibiotics among them. Here, we show that the impact exerted by antibiotics does not have the same magnitude on the different bacterial families. Future studies including more individuals and more different antibiotics could be used to support the choice of antibiotics that minimize the impact on some specific bacterial species that are associated with pathological outcome, such as Akkermansia muciniphila in the metabolic syndrome or in the response to anticancer treatment. ,

Although both antibiotics share similar clinical indications, only few studies compared the impact of ceftriaxone and cefotaxime on the gut microbiota in humans. , They suggested a higher impact for ceftriaxone than for cefotaxime, but ceftriaxone was administered at a higher daily dose and cefotaxime at a lower daily dose than in the usual clinical practice. The assumption of a maximal concentration reached at the end of the treatment was supported by a previous modeling analysis of moxifloxacin, a quinolone antibiotic that is not hydrolyzed by betalactamases. In this model, a more mechanistic model was developed to model drug pharmacokinetics in plasma and gut. Consistent with the basic model used here, fecal concentrations increased steadily during the 5‐day treatment period, and were then eliminated with a first‐order process after treatment cessation. Unfortunately, the absence of fecal concentrations precluded the analysis of a delay between the end of the antibiotic treatment and the decrease of the antibiotic gut exposure. Although such delay was observed in the model of moxifloxacin pharmacokinetics in feces, it is not expected in the present case with betalactam antibiotics. Indeed, the quick degradation into inactive metabolites by endogenous fecal betalactamases led to undetectable concentrations in feces at all sampling times.

We also used a simple mathematical form to relate the antibiotic gut exposure to bacterial killing rate, that was assumed similar for ceftriaxone and cefotaxime. This assumption was based on the absence of difference both on specific bacteria studied using phenotypic methods, and on the bacterial diversity of the gut microbiota. Consistently, we observed that the maximal reduction of the absolute counts of the bacterial families studied was not different between treatment groups. However, consistent with literature, we found that the gut elimination of ceftriaxone was five times slower than for cefotaxime. , Consequently, the time to return to 95% of the baseline counts was significantly longer in ceftriaxone‐treated subjects than for cefotaxime‐treated subjects.

In a previous analysis of the data using global measures of antibiotics impact on the microbiota, we did not find any significant differences between the two antibiotics. Using a more precise modeling approach, we evidenced here a higher impact of ceftriaxone on the bacterial composition of the gut microbiota for some bacterial families. This might translate into a higher impact in terms of selection of antibiotic resistance. Indeed, previous studies reported that the time window during which the microbiota is exposed to sub‐inhibitory concentrations may favor the selection of resistance and induce genetic modifications in response to antibiotic selective pressure. The very low level of resistance observed in the CEREMI trial, probably related to the low selective pressure on the healthy volunteers included in the analysis, did unfortunately not allow to investigate this specific question.

Finally, we provided here a modeling framework to study the impact of antibiotics on the intestinal microbiome, and delineate the interactions within taxonomic bacterial groups within the gut microbiota. Nonlinear mixed‐effect modeling coupled to the generalized Lotka‐Volterra equations appears as a powerful tool to evaluate interactions among bacterial groups and to quantify their susceptibility to antibiotic. Such framework could particularly be useful to assess the impact of typical antibiotic treatments on the microbiome, including for various treatment duration or dosing schedule.

CONFLICT OF INTEREST

The authors declared no competing interests for this work. As Editor‐in‐Chief of CPT: Pharmacometrics & Systems Pharmacology, France Mentré was not involved in the review or decision process for this paper.

AUTHOR CONTRIBUTIONS

C.B., J. Gue, and J. Guk wrote the manuscript. C.B., N.G., F.M., and X.D. designed the research. C.B., X.D., N.G., A.B.‐N., M.M., O.T., E.D., O.C., E.R., and C.d'H. performed the research. C.B., J. Gue, J. Guk, and F.M. analyzed the data.

Supporting Information S1

Click here for additional data file.