Population Density Modulates Drug Inhibition and Gives Rise to Potential Bistability of Treatment Outcomes for Bacterial Infections

The inoculum effect (IE) is an increase in the minimum inhibitory concentration (MIC) of an antibiotic as a function of the initial size of a microbial population. The IE has been observed in a wide range of bacteria, implying that antibiotic efficacy may depend on population density. Such density dependence could have dramatic effects on bacterial population dynamics and potential treatment strategies, but explicit measures of per capita growth as a function of density are generally not available. Instead, the IE measures MIC as a function of initial population size, and population density changes by many orders of magnitude on the timescale of the experiment. Therefore, the functional relationship between population density and antibiotic inhibition is generally not known, leaving many questions about the impact of the IE on different treatment strategies unanswered. To address these questions, here we directly measured real-time per capita growth of Enterococcus faecalis populations exposed to antibiotic at fixed population densities using multiplexed computer-automated culture devices. We show that density-dependent growth inhibition is pervasive for commonly used antibiotics, with some drugs showing increased inhibition and others decreased inhibition at high densities. For several drugs, the density dependence is mediated by changes in extracellular pH, a community-level phenomenon not previously linked with the IE. Using a simple mathematical model, we demonstrate how this density dependence can modulate population dynamics in constant drug environments. Then, we illustrate how time-dependent dosing strategies can mitigate the negative effects of density-dependence. Finally, we show that these density effects lead to bistable treatment outcomes for a wide range of antibiotic concentrations in a pharmacological model of antibiotic treatment. As a result, infections exceeding a critical density often survive otherwise effective treatments.

In the turbidostat, cell density n is held constant using feedback control with peristaltic pumps (see Methods). Once cells have grown to the appropriate density, drug is added at concentration D in to both the media feed reservoir and the culture vial. The drug concentration D in the reservoir is then given by where F is the flow rate of the pumps required to maintain constant density, V is the volume of the culture chamber, and f (n, D) is an unspecified function that accounts for the effective change in drug concentration due to cell density.
Here we take f (n, D) = − Dn j , where is a rate constant and j is an exponent that describes the kinetics of drug decay (i.e. j = 1 indicates decay that is linear in n, while j = 2 indicates quadratic decay). This form (with j = 1 or j = 2) provides a good description of the turbidostat data for the drugs in this study (Figures E and F in S1 Text).
To maintain constant cell density, the flow rate must be related to per capita growth rate g by g = F/V , and Equation S1 therefore becomes which has steady state solution (S3) The growth rate g, which we measure in dimensionless units of per capita growth rate in the absence of drug, is also dependent on the steady state drug concen-tration, D, according to the Hill-like dose response function where h is a Hill coefficient and D is measured in units of IC 50 . We can combine Equation S3 and Equation S4 to arrive at a single nonlinear equation for g (or equivalently, for D). Specifically, we have (S5) Equation S5 can be solved numerically to find g for any given values of n, D in and the parameters h, j, and . In principle, Equation S5 can have multiple solutions, but we do not find evidence for multiple solutions (g ≥ 0) over the parameter range of our measurements. The parameter h can be determined from standard low-density dose response curves (for example, h ≈ 2.1 for tigecylcine; see main text), and is a free parameter that we estimate using turbidostat data (see Table B in S1 Text).

Parameter Estimation and Model Selection
To statistically compare the the two models for drug decay (j = 1 or j = 2), we use standard model-selection techniques [68]. Specifically, we assume that the experimental errors are independent and Gaussian distributed with unknown variance σ 2 . We then calculate for each model the Akaike Information Criteria, which is given by AIC = −2 log(L(ĉ|y)) + 2m (S6) where log(L(ĉ|y)) is the log likelihood function, y is the turbidostat growth data, c is maximum likelihood estimate of the free parameters of the model (in this case, the decay rate constant and, implicitly, the unknown error variance), and m is the number of free parameters (m = 2 for both models). Note that the maximum likelihood estimate is the same as the least squares estimate for the Gaussian noise case. The AIC is an estimate of the expectation value of the relative Kullback-Leibler (KL) divergence between the fitted model and the "true mechanism" generating the observed data. The model with the lowest AIC value among a set of models is considered the best model in that it minimizes the KL divergence between the model and statistical mechanism underlying the data. For independent Gaussian errors, AIC reduces (up to an additive constant) to AIC = −N log(σ 2 ) + 2m, where N is the number of observations andσ 2 is the maximum likelihood estimate of the variance. In practice, we use a small sample estimator of AIC that includes a bias correction term AIC = −2 log(L(ĉ|y)) + 2m + 2m(m + 1) The differences in AIC values between the two models can be converted to an Akaiki weight, where δ ≡ AIC 2 − AIC 1 , where AIC j is the AIC of model j. Because exp(−δ/2) is proportional to the likelihood of the model j given the data, the weight w can be interpreted as a measure of the evidence in favor of one model or the other. Parameter estimates, AIC values, and weights for each model are given in Table SB.
2 PK/PD Model for Infection Dynamics

Model Definition
To model infection dynamics, we assume a simple PK/PD model similar to those in [6,41]. Specifically, we assume cell density n and effective antibiotic concentration D are governed by where g(D) is a growth rate function that accounts for the effects of antibiotic, C is the carrying capacity of the environment, and k d is the natural decay rate of the antibiotic in the clinical setting of interest (e.g. in the bloodstream of the patient). To account for the measured density dependence of drug efficacy, the equation for D involves a density dependent decay term similar to Equation S2. The growth rate function g(D) describes the impact of drug on division rate and is given by [6,41] where g max is the maximum growth rate of the population (which we can set equal to 1 without loss of generality by measuring time in units of g −1 max ), g min < 0 is the maximum rate of kill of the drug, h is a Hill-like steepness coefficient, and K 0 is the minimal inhibitory concentration (MIC) of the drug. The growth function, Equation S11, is a monotonically decreasing sigmoidal function that equals zero at D = K 0 and asymptotically approaches g min for high drug concentrations.

Periodic Drug Dosing
To mimic clinical dosing protocols [6,41], we assume that drug is added at a concentration D 0 at periodic intervals with length T . In what follows, we outline an adiabatic approximation to this model valid when the dynamics of D(t) occur on a faster timescale than those for n(t). This approximation provides a qualitative description of the system dynamics and is quantitatively accurate when |g min | << 1.

Adiabatic Approximation
To make analytical progress, we assume that n(t) changes slowly on the timescale of D(t). As a result, the drug decays exponentially over each period according to where k(n) ≡ k d − n j , and D(t) is reset to D 0 at the start of the next dosing period. Because we are interested in the long-time behavior of the system spanning many periods T , we approximate the dynamics of Equation S10 by inserting Equation S12 into Equation S10 for ∂n ∂t and then averaging over one period to arrive at where brackets represent an average over one period T . This amounts to an adiabatic approximation that averages over dynamics on the faster timescale governing D(t). We expect this approximation to be valid when the system requires multiple periods T to reach its steady state (n = 0 or n = C), and we find that it is increasingly accurate when |g min | << 1, which achieves the desired separation of timescales. Specifically, we have where β(n) ≡ k(n)T . Equation S14 follows from straightforward integration of the growth function over one period followed by simple algebraic manipulation (note that this form is valid only for g min < 0, but the general case can be written in terms of logarithms rather than the tanh −1 function).

Fixed Points and Linear Stability
Using Equations S13, S14, one can show that the system has three fixed points (n 0 , n 1 , n 2 ) given by n 0 = 0, where γ(n) is given by The fixed points occur when the population is extinct (n 0 ) or has reached carrying capacity (n 2 ). However, the density-dependence of the MIC gives rise to a third non-trivial fixed point (n 1 ) given by the (physically relevant) solutions to D 0 = K 0 γ(n), which ensure that g(n) = 0. By linearizing Equation S13 around the fixed points in Equation 15, it is straightforward to show that n 0 is stable for D 0 > K 0 γ(0) and unstable for D 0 < K 0 γ(0). Similarly, fixed point n 2 is stable for D 0 < K 0 γ(C) and unstable for D 0 > K 0 γ(C). Finally, the physically meaningful values of the nontrivial fixed point n 1 are unstable for > 0. On the other hand, if < 0, the non-trivial fixed point becomes stable over the physically relevant range.

Bistability
Interestingly, this stability analysis points to a region of bistability for K 0 γ(0) < D 0 < K 0 γ(C). The size of the bistable region is given by ∆D = K 0 (γ(C)−γ(0)). Because γ(n) ∼ exp(β(n)) for large n, the range of bistability can be extremely large, often many orders of magnitude larger than K 0 (see Figure 5). Examples of the bifurcation diagram for different values of g min and k 1/2 are shown in Figure 5 and Figure G in S1 Text. Finally, we note that bistability cannot occur when < 0, as for ampicillin. Instead, there is an intermediate range of drug concentrations for which all initial densities approach an intermediate, stable population size (Figure G in S1 Text). Table   Table A in S1 Text: Drugs Used in this Study.     Figure D in S1 Text: Buffered media does not affect drug-free growth but alters density dependence of ceftriaxone. A. Left: Growth rate of low density cultures in the absence of drug in regular media (blue) and highly buffered media (red). Growth in the presence of tigecycline at low densities (middle) and high densities (right). Growth in buffered and regular media did not differ significantly in the absence of drug, but in the presence of tigecycline, buffered media leads to slightly (low density) and then dramatically (high density) decreased growth rate. B. Relative growth rate of cells exposed to ceftriaxone (50 ug/mL) in regular media (red), strongly buffered media (black, dashed), and in low-density (OD=0.2) cultures with externally modulated pH (blue, dotted). In the latter case, buffer was supplemented with HCl to achieve pH = 7.5, 6.8, and 6.0, which correspond to pH of steady state cultures held at OD=0.2, 0.5, and 0.8, respectively. See also  Figure 2, main text). For tigecycline, the parameter ε was fit from turbidostat data, while K 0 (IC 50 ) and h were determined independently from low density dose response curves (see Figure F in S1 Text). For the remaining drugs, ε, K 0 , and h were fit simultaneously from the turbidostat data because low density dose response curves were not measured, though similar results were obtained if K 0 and h were determined from the OD=0.2 data alone.  Figure F in S1 Text: Model fits to experimental data. A. Relative growth rate of V583 populations as a function of tigecycline in regular media (BHI). Circles, growth rates measured in low-density (OD<0.2) exponentially growing populations. Solid line, fit to dose response curve ! = (1 + (!/!) ! ) !! . For tigecycline in regular media, K =19.3±1 and Hill coefficient h=2.1±0.1. We also found a similar dose response curve for tigecycline in highly buffered media, but K is slightly lower (K =11.1±0.4 and Hill coefficient h=2.1±0.2). For each concentration, we estimated growth rate by fitting the time series of optical density in the low density region to an exponential function using nonlinear least squares fitting. Inset: example low-density growth curves for [Tigecycline]=10, 40, and 60 ng/mL. Red points: OD measurements. Green lines: fit to exponential function. B. and C. Comparison of linear (dashed lines) and quadratic (solid lines) decay decay models to experimental data (circles) for tigecycline (B) and doxycycline (C). The former is best fit by the quadratic model, while the latter is best fit by the linear model according to AIC model selection (see Table B of S1 Text).  Figure G in S1 Text: PK/PD Phase Diagrams. A. Top panels show phase diagrams (i.e. bifurcation plots) for linear decay model (red) and quadratic decay model (blue) based on tigecycline growth data for low (left) and high (right) maximum kill rate. Left panel: g min =-0.05; right panel, g min =-1. Solid lines, stable fixed points of population density (theory). Dashed lines, unstable fixed points (theory). The curved dashed line is the phase boundary (separatrix) indicating the critical density above which a population will survive. The decay rate is k d = ½ in both plots. Lower panels show phase diagrams (quadratic model) for low (k d = ¼, left) and high (k d = 1, right) values of the native drug decay rate. In both plots, g min =-0.05. ε=0.9 (quadratic model) or ε=0.5 (linear model) and h=2 in all plots. B. Main panel: Theoretical (solid and dashed lines) and numerical (shaded region) phase diagrams indicate treatment outcomes in PK/PD model as a function of initial cell density (ranging from 0 to the carrying capacity, C=1.3) and initial antibiotic concentration D 0 . Solid red lines, stable fixed points of population density (theory). Dashed red lines, unstable fixed points (theory). Shaded regions indicate that initial population grows, while unshaded regions indicate that initial population shrinks in numerical simulations of the PK/PD model. Upper right inset: numerical solution of PK/PD equations for five different initial densities (indicated by black squares on the phase diagram). Lower inset: temporal dynamics of antibiotic concentration. For numerical phase diagram and simulations, g min =-0.1 (units of hrs -1 ) and ε=-0.2.