Understanding mechanisms of bacterial eradication is critically important for overcoming failures of antibiotic treatments. Current studies suggest that the clearance of large bacterial populations proceeds deterministically, while for smaller populations, the stochastic effects become more relevant. Here, we develop a theoretical approach to investigate the bacterial population dynamics under the effect of antibiotic drugs using a method of first-passage processes. It allows us to explicitly evaluate the most important characteristics of bacterial clearance dynamics such as extinction probabilities and extinction times. The new meaning of minimal inhibitory concentrations for stochastic clearance of bacterial populations is also discussed. In addition, we investigate the effect of fluctuations in population growth rates on the dynamics of bacterial eradication. It is found that extinction probabilities and extinction times generally do not correlate with each other when random fluctuations in the growth rates are taking place. Unexpectedly, for a significant range of parameters, the extinction times increase due to these fluctuations, indicating a slowing in the bacterial clearance dynamics. It is argued that this might be one of the initial steps in the pathway for the development of antibiotic resistance. Furthermore, it is suggested that extinction times is a convenient measure of bacterial tolerance.

1. Introduction

The rise of pathogenic bacteria that are resistant to antibiotics is one of the major global health threats of the twenty-first century. High mortality rates and increasing healthcare costs associated with fighting bacterial infections call for designing new effective therapeutic strategies [1,2]. A major challenge in overcoming treatment failures comes from ineffective eradication of antibiotic-susceptible bacteria [3–5]. Despite the introduction and wide application of a very large range of antibiotics since the 1940s, important aspects of how antibiotics clear bacterial population at all levels (molecular, cellular and population) remain unclear. A deeper understanding of the underlying dynamics of bacterial clearance requires not only extensive laboratory studies but also the development of new theoretical approaches to investigate the bacterial response to antibiotics [6].

Majority of current experimental and theoretical studies focus on the eradication of initially large quantities of bacteria [7–9], and it was shown that a deterministic picture describes well the decrease in these bacterial populations [9,10]. In this deterministic framework, the dynamics of bacterial populations exposed to an antibiotic is characterized by a minimum inhibitory concentration (MIC), the minimal drug concentration required to inhibit bacterial growth [9–11]. The MIC can be regarded as a threshold on the antibiotic concentration such that only above the MIC a bacterial population can undergo full extinction, while for concentrations below the MIC the infection will never disappear.

However, it can be argued that it is also critically important to investigate the clearance dynamics for small bacterial populations. Failure to completely eradicate a population of bacteria can have two main consequences. First, even a small number of surviving bacteria can restore infections [12–14]. There are indications that as few as 10–100 bacterial cells, such as Salmonella or Shigella, are enough to restart the infection. Second, certain strains of surviving cells may develop antibiotic resistance, which, in turn, can complicate subsequent therapies [15–17]. Therefore, the effective treatment of infections requires not only the reduction of a large population number to a small number but also the complete eradication of the bacterial population [18–20].

Despite earlier technical problems [7,8], recent experiments were able to quantitatively investigate the antibiotic-induced clearance of small bacterial populations [21]. It was demonstrated that stochastic factors play much more important roles under these conditions. For example, Coates et al. [21] showed that even in sub-MIC antibiotic concentrations, bacterial populations decline with non-zero probability. This means that under the same conditions some populations experience growth with cells continuously dividing, while other populations quickly become extinct. A Markovian probabilistic birth-and-death model was introduced to uncover the relationship between the extinction probability and the antibiotic concentration [21]. This stochastic approach predicted that antibiotics induce fluctuations in bacterial population numbers. These fluctuations, in turn, lead to the stochastic nature of the clearance of small bacterial populations.

Although the Markovian model developed by Coates et al. successfully described the experimental observations, it could not predict an extinction time, i.e. the mean time at which the given number of bacterial cells will go to zero. This is a very important property of bacterial population clearance dynamics because it gives a better measure of the efficiency of the antibiotic treatments than the extinction probability. One could use an analogy with thermodynamic and kinetic descriptions of chemical processes. Thermodynamics gives the probability for the process to happen, but if the process is actually taking place in real time, it is determined by kinetic rates. In our language, this means that the large extinction probability might not always correlate with fast removal of bacterial infection. While the extinction probability can give a qualitative measure of the bacterial population dynamics, the extinction time is much more useful in the quantitative characterization of bacterial resistance and tolerance. It seems that the development of new drugs and new therapies in fighting against bacteria should use this quantity as a measure of their success.

In this study, we developed a discrete-state stochastic model of the antibiotic-induced clearance of bacteria that employs a method of first-passage probabilities. This method has been successfully used to analyse multiple processes in chemistry, physics and biology [22–24]. It allows us to quantitatively describe the stochastic dynamics of bacterial eradication by explicitly calculating extinction probabilities and extinction times and clarifying the physical meaning of the MIC. Our method is also applied to investigate the effect of fluctuations in the growth rates on the stochastic clearance of bacterial populations. These fluctuations can be attributed to various environmental factors such as the availability of nutrients, changes in osmolarity and other factors [25]. Our results suggest that these fluctuations influence the extinction probabilities and extinction times differently. There is a large range of antibiotic concentrations when the extinction times increase due to fluctuations, and this corresponds to the slowdown of the dynamics of bacterial eradication. We speculate that this might be a first step in developing antibiotic resistance. It is also argued that extinction times is a convenient new measure of bacterial tolerance.

2. Model

2.1. Stochastic clearance with a constant growth rate

We start our analysis by considering a simple stochastic model for the clearance of bacteria as shown in figure 1a. Our goal is to obtain a minimal theoretical description of bacterial clearance dynamics. For this reason, the model is characterized by only two parameters: the rate of cell growth λ and the rate of cell death ϕ (figure 1a). The bacterial growth rate is generally controlled by environmental factors such as the availability of nutrients, temperature, osmotic pressure and other factors [25]. When exposed to antibiotics, the cell growth rate can also depend on the antibiotic concentration [26]. Also, as a possible mechanism of antibiotic resistance, bacteria can sequester, degrade or modify antibiotics [27,28]. This, in turn, might significantly complicate the relationships between the cell death and growth rates. For the sake of simplicity, we assume that the cell growth rate is independent of antibiotic concentration and remains constant over different generations, while the cell death rate, ϕ, is controlled by the antibiotic concentration. It is also assumed here that if the bacterial population reaches size N, the organism hosting the bacteria will change its behaviour and it goes into another metabolic state. This different metabolic state might correspond, for example, to the situation when bacteria produce so many toxins that they damage normal cell membranes and/or inhibit normal protein synthesis in the host organisms. It might also describe a state in which the host organism mounts an active immune response, strongly modifying all metabolic processes in the system [29]. Also, the bacteria might achieve antibiotic resistance or the organism might even die from the infection. This is known as a fixation. Figure 1. (a) Schematic of the single growth-rate model for the clearance of bacteria. Each state n (n = 0, 1, …, N) represents a bacterial population with n cells. The states 0 and N correspond to the bacterial eradication (no cells in the system), and the fixation, respectively. From each state n, the bacterial population can change to the state n + 1 (growth) with a total rate nλ, or it can jump to the state n − 1 (shrinking) with a total rate nϕ. We define the normalized death rate, x as the ratio of death rate and growth rate, x = ϕ/λ. Analytical calculations of extinction probabilities (b) over n − x parameter space for N = 50; (c) for a specific mid-size inoculum (n = N/2) over N − x parameter space. (Online version in colour.)

To describe dynamical transitions in the system, we define F n (t) as a probability density function to clear the system from infection at time t if the initial population number (so-called inoculum size) is equal to n (1 ≤ n ≤ N − 1). The temporal evolution of this probability function is governed by the following backward master equation [23,24]:

d F n ( t ) d t = n ϕ F n − 1 ( t ) + n λ F n + 1 ( t ) − n ( λ + ϕ ) F n ( t ) . 2.1

F n ( s ) ~ = ∫ 0 ∞ F n ( t ) e − s t

( s n + λ + ϕ ) F n ( s ) ~ = ϕ F n − 1 ( s ) ~ + λ F n + 1 ( s ) ~ . 2.2

F n ~ ( s ) ≃ f n − s b n . 2.3

F n ~ ( s = 0 ) = f n

f n = x N − x n x N − 1 , 2.4

Introducing the Laplace transform of this function,, we transform the backward master equation intoBecause we are mostly interested in the stationary dynamic behaviour at long times (→ 0), the following expansion can be written:Thenyields the first-passage probability of bacterial clearance or simply the extinction probability for the bacterial population with inoculum size. It can be shown that the extinction probability is given by (see appendix A.1 for details)where a parametercan be viewed as an effective death rate for the bacterial population normalized over the growth rate. Since in our model it is assumed that the growth rate does not depend on the death rate, the extinction probability is determined only by the ratio ofand λ.

Our analytical results for the extinction probability are presented in figure 1. The dependence of the bacterial clearance probability (from equation (2.4)) on the initial size of the bacterial population and on the values of x is given in figure 1b. For x = 1, which corresponds to ϕ = λ, the cell growth and death rates are equal to each other. In this state, which in the deterministic picture of bacterial clearance is described as the MIC, the extinction probability linearly decreases with the inoculum size, f n = (N − n)/N. In this case, the growth and the death rates are the same, and the probability of bacterial clearance is proportional to the relative distance from the initial state n to the fixation state N. The smaller the inoculum size, the larger the probability to eradicate the infection. But even for n = 1, the extinction probability is not equal to one [f 1 (x = 1) = (N − 1)/N < 1]. For x < 1 (sub-MIC conditions), the extinction probability is a decaying function of the inoculum size n. In this case, the growth rate is faster than the death rate, and the larger the inoculum size, the harder for the system to reach the total eradication of the infection (n = 0 state). One could also see this more clearly in the limit of x → 0 and N → ∞ when we have f n ≃ xn. This implies that even for sub-MIC conditions (low antibiotic concentrations) the extinction probability is never equal to zero, which is a clear signature of the stochastic effects in the bacterial clearance dynamics. The situation is different for x > 1 (large antibiotic concentrations), when the extinction probability is always close to one except in the region near the fixation state N. This can be also seen from the case of x ≫ 1 and N → ∞ when we obtain f n ≃ 1 − xn−N. This result suggests that even for concentrations above MIC the extinction probability is never equal to one, which is again due to the stochastic fluctuations in the system. Our analytical calculations are verified with Monte Carlo computer simulations, in which we used typical doubling times associated with bacteria E. coli, in range from 20 to 300 min [25].

The stochastic effects of the bacterial clearance can be understood better if we consider the extinction probability of a specific inoculum size (n = N/2), equally distant from the state n = 0 (eradication) and n = N (fixation), which is plotted in figure 1c. One can see that the dependence of the extinction probability on x follows a logistic sigmoid curve. The steepness of the curve at the midpoint (x = 1) is controlled by the values of n and N. In other words, for x < 1 the extinction probability is still non-zero, while for x > 1, it is still less than one. Therefore, x = 1 does not satisfy the classical definition of the MIC as the MIC required for clearance. Thus, we need to calculate the effective saturation value for which the extinction probability becomes very high and realistically not much different from one. This might be viewed as an effective MIC for stochastic bacterial clearance. This saturation point is given by (see appendix A.2)

x sat = 1 + N n ( N − n ) . 2.5

sat

For the special case/2, this equation yields= 1 + 4/. Therefore, asincreases the steepness of the curve becomes sharper, such that the extinction probability becomes insensitive to population number while it is ultrasensitive with respect to. In this case, the large population alleviates the stochastic effects in the bacterial clearance, and= 1 yields the MIC, as expected.

Theoretical calculations also predict that the extinction probability strongly depends on inoculum size and on its relative distance to the fixation state N, as illustrated in figure 1b. For n = 1, the dependence on x is linear for small antibiotic concentrations (x < 1), while n = N − 1 is almost zero for x < 1 and it is slowly approaching one for larger antibiotic concentrations. These different behaviours are again a consequence of the stochastic nature of bacterial population clearance. A critically important property of bacterial eradication is how long does it take to clear the infection from the host, which is known as the extinction time. This timescale is crucial for the development of new therapies and it can be also useful in quantifying bacterial tolerance, which is the ability of a bacterial population to survive at longer periods of time exposed to antibiotics [30]. Our first-passage probabilities method is a powerful tool to evaluate this quantity. We define T n as a mean first-passage time to reach the extinction state (n = 0) from the inoculum of size n, and this is exactly the extinction time. Using the probability density function F n (t), it can be written as

T n = ∫ t = 0 ∞ t F n ( t ) d t ∫ t = 0 ∞ F n ( t ) d t . 2.6

T n = − ∂ F n ~ / ∂ s | s = 0 F n ~ ( s = 0 ) = b n f n . 2.7

T n = 1 λ ( x N − x n ) ( x − 1 ) [ 1 − x n 1 − x N ∑ k = 1 N − 1 ( x N − x k ) ( x N − k − 1 ) k − ∑ k = 1 n − 1 ( x N − x k ) ( x n − k − 1 ) k ] . 2.8

T n = 1 λ ( N − n ) [ n N ∑ k = 1 N − 1 ( N − k ) 2 k − ∑ k = 1 n − 1 ( N − k ) ( n − k ) k ] . 2.9

T n = 1 λ [ x n − 1 x − 1 ln ( x x − 1 ) − ∑ k = 1 n − 1 ( 1 k ∑ j = 0 n − k − 1 x j ) ] , 2.10

T n ≃ 1 λ [ 1 n + x n + 1 + ⋯ ] . 2.11

Using the Laplace transform and equation (2.3), we obtainAs explained in appendix A.1, the extinction time is explicitly given byIt can be shown that for= 1 that the expression for the extinction time takes the formFor> 1 and→ ∞, the extinction times are given by (see appendix A.1)while for→ 0 we have

The results of our calculations for the extinction times are presented in figure 2. As expected, it takes longer to clear the infection for larger inoculum sizes (figure 2). For large antibiotic concentrations (x > 1), the extinction time is shorter and it depends weaker on the inoculum size n. For small antibiotic concentrations (x < 1), the time to eradicate the infection is larger and it is more sensitive to the inoculum size. More interesting behaviour is observed when we analyse the extinction time for different antibiotic concentrations (figure 2b). A non-monotonic behaviour as a function of x is predicted, and the largest extinction time is observed for MIC conditions (x = 1). Increasing the antibiotic concentrations (x > 1) shortens the time for bacterial clearance because the drive to infection eradication becomes stronger. However, the surprising observation is that lowering the antibiotic concentrations below MIC (x < 1) can also accelerate the bacterial clearance despite the fact that the probability of clearance decreases. This can be explained by the following arguments. At these conditions, only those bacterial populations lead to the full eradication that rapidly shrink. If it is not fast, the shrinking of the bacterial population will be reversed and the infection will spread again. This is another non-trivial signature of the stochastic effects in the bacterial clearance dynamics. However, we should also emphasize here that increasing parameter N lowers the stochastic effects. Figure 2. Analytical calculations for the extinction times (in minutes): (a) as a function of the inoculum size for three different values of x; and (b) as a function of the parameter x for different inoculum sizes (n = 10, 25 and 40). In all calculations N = 50 and λ = 1/60 min−1 were used. (Online version in colour.)

Our analysis of extinction times allows us to reinterpret the meaning of MIC. For N → ∞ from (2.10) we conclude that the extinction time diverges logarithmically for x → 1, and it becomes infinite for x < 1. This suggests a new more practical definition of the MIC (x = 1). It is the antibiotic concentration at which the extinction time is maximal (for finite bacterial populations), or it is the antibiotic concentration below which the extinction times diverge (for N → ∞). This analysis also suggests that, from the practical point of view, to eliminate the infection it is important to apply the antibiotic concentrations that significantly differ from the MIC to avoid the slowdown in the dynamics.

It is interesting to compare our theoretical predictions with experimental measurements of stochastic bacterial clearance [21]. In these experiments, the stochastic population dynamics of bacteria exposed to bactericidal drugs have been monitored starting from single E. coli bacteria for sub-MIC conditions (x = 0.8) and for concentrations above the MIC (x = 1.2). It was also estimated that the growth rate is λ ≃ 1/100 min−1. Then using (2.10) and (2.11), we predict that for both cases, x = 0.8 and x = 1.2, the extinction times are close to 200 min, which agrees well with these experimental observations.

2.2. Stochastic clearance with fluctuations in the growth rate

Although the mechanisms of the development of antibiotic resistance remain not fully understood, recent studies suggest that random fluctuations of various parameters can stimulate the bacterial tolerance to antibiotic drugs [6,31]. Bacterial population dynamics are subject to intrinsic stochasticity because of variations in gene expression and subject to extrinsic stochasticity due to environmental variations. For example, single-cell experiments have shown that the duration of the cell cycle is subject to random fluctuations [25,32]. In this case, cell cycle duration follows a distribution with certain variance. We can investigate the effect of growth-rate fluctuations on the bacterial clearance dynamics using our theoretical first-passage probabilities method. To do so, we introduce a simple model as shown in figure 3. It is assumed that the infection can spread with two growth rates, λ 1 and λ 2 , while the death rate ϕ is assumed to be the same in both populations. The system can stochastically transition between two different growth regimes with rates δ and γ (figure 3). For the sake of simplicity, in calculations, we assume that δ = γ. Similar deterministic models for population dynamics in fluctuating environments have been already discussed [33–35]. Figure 3. Schematic of the model for the clearance of bacteria with fluctuating growth rates. The model comprises two coupled lattices. At each state n on lattice 1 (lattice 2), population can jump to state n + 1 with growth rate nλ 1 (nλ 2 ). Death rates are equal along the lattices. Also, δ and γ are rates to transition between lattices.

In this model, we define F n ( 1 ) ( t ) and F n ( 2 ) ( t ) as the probability density functions to clear the system from infection if the bacterial population starts with n cells while growing with the rate λ 1 or λ 2 , respectively. The temporal evolution of these probability functions is governed by the following backward master equations:

d F n ( 1 ) ( t ) d t = n ϕ F n − 1 ( 1 ) ( t ) + n λ 1 F n + 1 ( 1 ) ( t ) + n γ F n ( 2 ) ( t ) − ( n δ + n ϕ + n λ 1 ) F n ( 1 ) ( t ) 2.12

d F n ( 2 ) ( t ) d t = n ϕ F n − 1 ( 2 ) ( t ) + n λ 2 F n + 1 ( 2 ) ( t ) + n δ F n ( 1 ) ( t ) − ( n γ + n ϕ + n λ 2 ) F n ( 2 ) ( t ) . 2.13

andIn general, it is difficult to obtain a full analytical solution for this problem for arbitrary. However, exact solutions for simple cases with= 2 and= 3 can be derived (see appendix A.3 for details).

To better understand the effects of fluctuation on the dynamics of clearance, it is convenient to compare the fluctuating growth model (rates λ 1 and λ 2 ) presented in figure 3 with a single growth-rate model with λ = (λ 1 + λ 2 )/2 presented in figure 1a. Since the average growth rates in both cases are the same, the possible differences in the dynamics properties for bacterial clearance are coming from the fluctuations. To quantify this effect, we define a function r n ( f ) as the ratio of the extinction probabilities predicted by the fluctuating-growth model and by the single growth-rate model:

r n ( f ) = f n ( avg ) f n = f n ( 1 ) + f n ( 2 ) 2 f n . 2.14

r n ( T )

r n ( T ) = T n ( avg ) T n = T n ( 1 ) + T n ( 2 ) 2 T n . 2.15

r n ( f ) > 1

r n ( T ) > 1

Similarly, one can define a functionfor the ratio of extinction timesIf, then it means that fluctuations increase the extinction probability, whileindicates that fluctuations increase the extinction times.

As shown in appendix A.3, the fluctuating growth-rate model has been solved exactly to evaluate extinction probabilities and extinction times for N = 3, and the results are presented in figure 4. It is found that r n ( f ) is always larger than one (figure 4a), which indicates that in the bacterial population with fluctuations in the growth rate the probability of eradication of infection is always larger than in the single growth population. The effect is stronger for not very large antibiotic concentrations and for slow transitions between two growth regimes. It can be argued that switching transitions opens new pathways for the eradication of the bacteria, and this should increase the extinction probability. At the same time, increasing the amplitude of the switching transition rates leads to an effective equilibrium single growth-rate regime with the growth rate given by the average between two dynamic regimes, and this clearly does not increase the extinction probability. Figure 4b presents the ratio of extinction times, and our theory predicts that r n ( T ) > 1 , i.e. fluctuations in the growth rates unexpectedly slow down the bacterial clearance dynamics, in contrast to expectations from the extinction probabilities. The effect is stronger for not very large antibiotic concentrations and it disappears for x → ∞. It is also strong for weak fluctuation rates between two growth regimes. This surprising result can be explained by noting that due to weak transition rates the system can be effectively trapped in the regime with smaller death rates, and this should slow down the bacterial clearance dynamics. Figure 4. Analytical calculations of dynamic properties for N = 3 model. (a) The ratio of the extinction probabilities as a function of x for three different values of the transition rates γ and δ. (b) The ratio of the extinction times as a function of x for three different values of transition rates γ and δ. In all calculations, n = 2, λ 1 = 3 60 min − 1 , λ 2 = 0.3 60 min − 1 and λ = (λ 1 + λ 2 )/2 were used. (Online version in colour.)

More realistic situations of bacterial population dynamics require consideration of systems with large N. Because analytical calculations cannot be done for these cases, we explored Monte Carlo computer simulations to evaluate the dynamic properties of stochastic bacterial clearance. The results are presented in figure 5. One can see that for relatively small antibiotic concentrations (x < 1), the fluctuations in the growth rate increase the extinction probability (figure 5a). In this case, which is generally unfavourable for the eradication of infection, opening new pathways should help to clear the infection. This is because the system can spend half of the time in the dynamic regimes with smaller death rates, which helps to fight the infection better. However, the situation changes for large antibiotic concentrations (x > 1), when the fluctuations decrease the extinction probability. In this case, due to switching transitions, the system spends half of the time in the dynamic regime where it is more difficult to eradicate the infection. Figure 5. Predictions of Monte Carlo computer simulations for the fluctuating growth-rate and single growth-rate models. (a) The ratio of extinction probabilities as a function of x for three different values of n; and (b) the ratio of extinction times as a function of x for three different values of n. In simulations, the following parameters were used: N = 20, λ 1 = 3 60 min − 1 , λ 2 = 0.3 60 min − 1 and δ = γ = 0.165 60 min − 1 . (Online version in colour.)

A more complex picture is observed when we analyse the ratio of extinction time (figure 5b). It is found that for small antibiotic concentrations and for very large antibiotic concentrations the fluctuations in the growth rates lead to slower bacterial clearance dynamics. Only for intermediate antibiotic concentrations around MIC (x ∼ 1), fluctuations might accelerate the removal of infection. Apparently, opening new pathways for x < 1 and x ≫ 1 regions lowers the drive to eradicate the infection because the system spends more time switching between different dynamic regimes and not shrinking the bacterial populations.

Analysing the dynamic properties of the fluctuating growth-rate model, two important observations can be made. First, the extinction probability and extinction time generally do not correlate with each other when the system experiences fluctuations between different growth regimes. Second, turning on the fluctuations in the growth rates of bacteria can significantly increase the tolerance to antibiotic drugs for a large range of parameters. It seems reasonable to speculate that bacteria might explore this option in fighting against antibiotics.

We theoretically investigated the clearance of bacterial populations under the effect of antibiotic drugs by concentrating on stochastic aspects of this process. To understand better the mechanisms of eradication of infection, a method of first-passage probabilities is introduced. This allows us to obtain a comprehensive description of bacterial clearance dynamics. Two important dynamic features, extinction probabilities and extinction times, are explicitly calculated. We also clarified the physical meaning of MIC in the systems where the stochasticity is more relevant. Furthermore, using our method, we investigated the effect of fluctuations in the growth rates on the bacterial population dynamics, and we find that these random fluctuations affect differently extinction probabilities and extinction times. For the single growth-rate model, our analysis shows that extinction probabilities depend strongly on the antibiotic concentration, the inoculum size and the distance to the fixation state N. But the stochastic effects show up in observations that, even for concentrations above the MIC, the extinction probabilities are not equal to one, while for concentrations below the MIC the extinction probabilities are not equal to zero. More complex behaviour is observed for extinction times. For finite-size bacterial populations, the extinction times show non-monotonic dependence on the antibiotic concentrations with the maximum at the MIC. The unexpected acceleration in the eradication of infection for concentrations below the MIC is explained by the fact that the successful events, which are rare at these conditions, must proceed very fast. For infinitely large bacterial populations, our calculations show that the extinction times increase with lowering of antibiotic concentration and diverge for concentrations at the MIC and sub-MIC. These properties of extinction times provide an additional way of defining the conditions corresponding to MIC.

By introducing a stochastic model in which bacteria can randomly switch between two growth rates, we investigated the effect of environment fluctuations in the bacterial clearance dynamics. Our analytical and computer simulations results predict that these switchings increase the extinction probabilities for low antibiotic concentrations and decrease them for high antibiotic concentrations. However, the effect of fluctuations in the growth rates on extinction times is more complex. With the exception of the intermediate concentrations around MIC, random switchings slow down the bacterial clearance dynamics.

Our calculations lead to several important conclusions. Extinction probabilities and extinction times generally do not correlate with each other, so it is dangerous to make predictions on bacterial population dynamics by considering only the extinction probabilities as typically done in the field. There is a significant range of parameters when the fluctuations in the growth rates lead to the overall slowing down in the eradication of the infection. Bacterial response to antibiotics is a complex process, which depends on genetic and environmental factors [36]. Some bacterial strains are difficult to eradicate because their clearance needs a higher level of antibiotics that are toxic to hosts. Such bacteria are commonly known as antibiotic-resistant. It is a very challenging task to uncover the mechanisms of the development of bacterial resistance. Our results suggest that one of the first steps in the resistance pathway might be due to fluctuations in growth rates, which would give bacteria additional time to find another means to avoid the effect of antibiotic drugs. Although at this moment, this is just pure speculation, it will be interesting to investigate this possibility with experimental methods and more advanced theoretical approaches.

Even at concentrations above the MIC, some bacteria survive short-term exposure to antibiotics before being affected by it. This ability of a bacterial population is known as tolerance [37]. In contrast to resistance, which is quantified by the MIC, tolerance is poorly characterized. The most commonly used approach for quantifying tolerance is the measurement of time–kill curves [38]. Recently, a new metric for bacterial tolerance has been introduced [30]. This new metric, known as the minimum duration for killing 99% of the population, MDK 99 , can be evaluated by statistical analysis of measurements. Our theoretical method provides the extinction time as a new measure of bacterial tolerance. The advantage of this approach is that it takes into account the stochastic features of the population dynamics and it gives the average dynamic property of the bacterial clearance, which might be much more useful for practical applications.

Data accessibility

This article has no additional data.

Authors' contributions

H.T. and A.B.K. designed the research; H.T. and A.B.K. performed the research; H.T. and A.B.K. wrote the manuscript.

Competing interests

We declare we have no competing interests.

Funding

The work was supported by the Welch Foundation (C-1559), by the NSF (CHE-1664218), and by the Center for Theoretical Biological Physics sponsored by the NSF (PHY-1427654).

Appendix A A.1. Exact solution for the single growth-rate model In this appendix, we present the details for calculations of the extinction probability and the extinction time. As given in the main text, the temporal evolution of the first-passage probability function is governed by the following backward master equation [23]: d F n ( t ) d t = n ϕ F n − 1 ( t ) + n λ F n + 1 ( t ) − n ( λ + ϕ ) F n ( t ) . A 1 Introducing the Laplace transform of this probability density function, F n ( s ) ~ = ∫ 0 ∞ F n ( t ) e − s t , we obtain ( s n + λ + ϕ ) F n ( s ) ~ = ϕ F n − 1 ( s ) ~ + λ F n + 1 ( s ) ~ . A 2 To solve this recurrence relation, it is convenient to write the following expansion: F n ~ ( s ) ≃ f n − s b n . A 3 Then F n ~ ( s = 0 ) = f n yields the extinction probability. To proceed further, we substitute (A 3) into (A 2): ( s n + λ + ϕ ) ( f n − s b n ) = ϕ ( f n − 1 − s b n − 1 ) + λ ( f n + 1 − s b n + 1 ) . A 4 Rearranging terms yields − s 2 b n n + s ( f n n − b n ( λ + ϕ ) ) + ( λ + ϕ ) f n = ϕ f n − 1 + λ f n + 1 − s ( ϕ b n − 1 + λ b n + 1 ) . A 5 Equating coefficients of s on both sides yields two recurrence relations ( λ + ϕ ) f n = ϕ f n − 1 + λ f n + 1 A 6 and f n n − ( λ + ϕ ) b n = − ϕ b n − 1 − λ b n + 1 . A 7 Equation (A 6) can be simplified as ϕ g n − 1 = λ g n , A 8 where g n = f n − f n − 1 . A 9 Solution of (A 9) is given by g n = ( ϕ λ ) n g 0 = x n g 0 , A 10 where x = ϕ/λ. To find constant g 0 , we perform summation over equation (A 9): ∑ k = 0 N − 1 g k = ∑ k = 0 N − 1 f k − ∑ k = 0 N − 1 f k + 1 = f 0 − f 1 + f 1 − f 2 + ⋯ + f N − 1 − f N = f 0 − f N = 1. A 11 Combining (A 10) and (A 11) yields ∑ k = 1 N − 1 g k = g 0 ∑ k = 1 N − 1 x k = 1. A 12 Then, g 0 is given by g 0 = 1 ( ∑ k = 1 N − 1 x k ) = x − 1 x N − 1 . A 13 Therefore, g n = x n g 0 = x n ( x − 1 ) x N − 1 . A 14 Now using (A 9), we obtain the extinction probability f n = 1 − ∑ k = 1 n − 1 g k = 1 − ( x − 1 x N − 1 ) ∑ k = 1 n − 1 x k = x N − x n x N − 1 . A 15 Introducing the Laplace transform of this probability density function,, we obtainTo solve this recurrence relation, it is convenient to write the following expansion:Thenyields the extinction probability. To proceed further, we substitute (A 3) into (A 2):Rearranging terms yieldsEquating coefficients ofon both sides yields two recurrence relationsandEquation (A 6) can be simplified aswhereSolution of (A 9) is given bywhere/λ. To find constant, we perform summation over equation (A 9):Combining (A 10) and (A 11) yieldsThen,is given byTherefore,Now using (A 9), we obtain the extinction probability To calculate b n , we use equation (A 7) f n n − ( λ + ϕ ) b n = − ϕ b n − 1 − λ b n + 1 . A 16 This recurrence relations can be simplified as 1 λ f n n = x K n − 1 − K n , A 17 where K n = b n + 1 − b n . A 18 It can be shown that the solution of equation (A 17) is given by K n = x n K 0 − 1 λ ∑ l = 0 n − 1 x l f n − l ( n − l ) . A 19 It is convenient to rewrite the summation in the following form ∑ l = 0 n − 1 x l f n − l ( n − l ) = ∑ l = 1 n x n − l f l l . A 20 Solution of the recurrence relation K n = b n+1 − b n takes the form b n = ∑ j = 0 n − 1 K j . A 21 Using boundary condition, we obtain b N = ∑ j = 0 N − 1 K j = 0 . To calculate constant K 0 , we perform summation over (A 19) ∑ j = 0 N − 1 K j = K 0 ∑ j = 0 N − 1 x j − 1 λ ∑ j = 0 N − 1 ∑ l = 1 j x j − l f l l . A 22 Thus, K 0 is given by K 0 = ∑ j = 0 N − 1 ∑ l = 1 j x j − l ( f l / l ) λ [ ( 1 − x N ) / ( 1 − x ) ] . A 23 Finally, combining (A 19) and (A 21) yields b n b n = K 0 [ 1 − x n 1 − x ] − 1 λ ∑ j = 0 n − 1 ∑ l = 1 j x j − l f l l . A 24 This recurrence relations can be simplified aswhereIt can be shown that the solution of equation (A 17) is given byIt is convenient to rewrite the summation in the following formSolution of the recurrence relationtakes the formUsing boundary condition, we obtain. To calculate constant, we perform summation over (A 19)Thus,is given byFinally, combining (A 19) and (A 21) yields Having determined f n and b n , we can now obtain the expression for the extinction time, T n = − ∂ F n ~ / ∂ s | s = 0 F n ~ ( s = 0 ) = b n f n . A 25 Using (A 15) and (A 24), we have T n = [ ( 1 − x n ) λ ( x N − x n ) ( 1 − x N ) ] ∑ j = 0 N − 1 ∑ l = 1 j x N + j − l − x j l − [ 1 λ ( x N − x n ) ] ∑ j = 0 n − 1 ∑ l = 1 j x N + j − l − x j l , A 26 which can be further simplified into T n = 1 λ ( x N − x n ) ( x − 1 ) [ 1 − x n 1 − x N ∑ k = 1 N − 1 ( x N − x k ) ( x N − k − 1 ) k − ∑ k = 1 n − 1 ( x N − x k ) ( x n − k − 1 ) k ] . A 27 When x = 1, this expressions yields T n = 1 λ ( N − n ) [ n N ∑ k = 1 N − 1 ( N − k ) 2 k − ∑ k = 1 n − 1 ( N − k ) ( n − k ) k ] . A 28 In the case of x > 1 and N → ∞, it can be shown that T n = 1 λ [ x n − 1 x − 1 ln ( x x − 1 ) − ∑ k = 1 n − 1 ( 1 k ∑ j = 0 n − k − 1 x j ) ] , A 29 while for x → 0 we have T n ≃ 1 λ [ 1 n + x n + 1 + ⋯ ] . A 30 A.2. Calculation of the saturation point for extinction probability Using (A 15) and (A 24), we havewhich can be further simplified intoWhen= 1, this expressions yieldsIn the case of> 1 and→ ∞, it can be shown thatwhile for→ 0 we have Since the extinction probability versus x follows a logistic sigmoid curve, we can define a saturation value of x for which the extinction probability saturates to higher values. There is not a unique way to define this saturation point. Here, we use a simple definition presented in [39,40]. In the simplest approximation, the saturation point is the value of x at which the straight line passing through the midpoint (x = 1), and with a slope equal to the first derivative of the extinction probability at this point, intersects with f n = 1. We start by taking the derivative of the extinction probability f n = (xN − xn)/(xN − 1) with respect to x. d f n d x | x = 1 = ( N x N − 1 − n x n − 1 ) ( x N − 1 ) − N x N − 1 ( x N − x n ) ( x N − 1 ) 2 = n ( N − n ) 2 N . A 31 Using this derivate value and coordinate of the midpoint (x = 1 and f n = 1/2), we can obtain the equation of the straight line passing from the midpoint. The equation of line is y = ax + b, where a = n(N − n)/2N. After some algebra, we obtain y = ( n ( N − n ) 2 N ) x + 1 2 − n ( N − n ) 2 N . A 32 Solution of this equation at y = 1 yields the saturation point x sat = 1 + N n ( N − n ) . A 33 This method only provides a first-order approximation for the saturation point. This approximation can be improved by evaluating higher order (second, third, or fourth) derivatives of f n . In this case, the straight line passes through the point at which higher derivatives are zero. A.3. Exact solution for the coupled-parallel lattice model Using this derivate value and coordinate of the midpoint (= 1 and= 1/2), we can obtain the equation of the straight line passing from the midpoint. The equation of line is, where)/2. After some algebra, we obtainSolution of this equation at= 1 yields the saturation pointThis method only provides a first-order approximation for the saturation point. This approximation can be improved by evaluating higher order (second, third, or fourth) derivatives of. In this case, the straight line passes through the point at which higher derivatives are zero. It is difficult to obtain a general analytical solution for equations (2.12) and (2.13). However, for the small population numbers, the exact solution can be derived. In the following, we present the details of our calculations for N = 3 model. Schematic of the coupled-parallel mode is shown in figure 6. Dynamics of this model is governed by following backward master equations: d F 1 ( 1 ) d t = ϕ F 0 + λ 1 F 2 ( 1 ) + γ F 1 ( 2 ) − ( δ + ϕ + λ 1 ) F 1 ( 1 ) , A 34 d F 1 ( 2 ) d t = ϕ F 0 + λ 2 F 2 ( 2 ) + δ F 1 ( 1 ) − ( γ + ϕ + λ 2 ) F 1 ( 2 ) , A 35 d F 2 ( 2 ) d t = 2 ϕ F 1 ( 2 ) + 2 δ F 2 ( 1 ) − ( 2 γ + 2 ϕ + 2 λ 2 ) F 2 ( 2 ) A 36 and d F 2 ( 1 ) d t = 2 ϕ F 1 ( 1 ) + 2 γ F 2 ( 2 ) − ( 2 δ + 2 ϕ + 2 λ 1 ) F 1 ( 1 ) . A 37 Performing the Laplace transform, we obtain ( s + λ 1 + ϕ + δ ) F 1 ( 1 ) ~ = ϕ + γ F 1 ( 2 ) ~ + λ 1 F 2 ( 1 ) ~ , A 38 ( s + λ 2 + ϕ + γ ) F 1 ( 2 ) ~ = ϕ + δ F 1 ( 1 ) ~ + λ 2 F 2 ( 2 ) ~ , A 39 ( s + 2 λ 1 + 2 ϕ + 2 δ ) F 2 ( 1 ) ~ = 2 ϕ F 1 ( 1 ) ~ + 2 δ F 2 ( 2 ) ~ A 40 and ( s + 2 λ 2 + 2 ϕ + 2 γ ) F 2 ( 2 ) ~ = 2 ϕ F 1 ( 2 ) ~ + 2 γ F 2 ( 1 ) ~ . A 41 Solving this system of four equations and four unknowns yields F 1 ( 1 ) ~ , F 1 ( 2 ) ~ , F 2 ( 1 ) ~ and F 2 ( 2 ) ~ . Expanding these functions in terms of s yields the extinction probabilities, f 1 ( 1 ) = ϕ ( 4 λ 2 ( γ δ + 2 γ ϕ + δ 2 + 2 δ ϕ + ϕ 2 ) + 4 ϕ ( γ + δ + ϕ ) 2 + 4 δ λ 2 2 + 4 λ 2 2 ϕ + 4 λ 1 Λ ) 4 λ 2 ϕ ( γ ( δ + 2 ϕ ) + ( δ + ϕ ) 2 ) + 4 ϕ 2 ( γ + δ + ϕ ) 2 + 4 λ 2 2 ( δ + ϕ ) 2 + Δ λ 1 + 4 λ 1 2 Ψ , A 42 f 1 ( 2 ) = ϕ ( 4 λ 2 ( γ δ + 2 γ ϕ + δ 2 + 2 δ ϕ + ϕ 2 ) + 4 ϕ ( γ + δ + ϕ ) 2 + 4 λ 1 2 ( γ + λ 2 + ϕ ) + 4 λ 1 Ω ) 4 λ 2 ϕ ( γ ( δ + 2 ϕ ) + ( δ + ϕ ) 2 ) + 4 ϕ 2 ( γ + δ + ϕ ) 2 + 4 λ 2 2 ( δ + ϕ ) 2 + Δ λ 1 + 4 λ 1 2 Ψ , A 43 f 2 ( 1 ) = 4 ϕ 2 ( γ 2 + 2 γ δ + 2 γ λ 2 + 2 γ ϕ + δ 2 + δ λ 1 + δ λ 2 + 2 δ ϕ + λ 2 ϕ + λ 2 2 + ϕ 2 ) 4 λ 2 ϕ ( γ ( δ + 2 ϕ ) + ( δ + ϕ ) 2 ) + 4 ϕ 2 ( γ + δ + ϕ ) 2 + 4 λ 2 2 ( δ + ϕ ) 2 + Δ λ 1 + 4 λ 1 2 Ψ A 44 and f 2 ( 2 ) = 4 ϕ 2 ( γ 2 + 2 γ δ + γ λ 1 + γ λ 2 + 2 γ ϕ + δ 2 + 2 δ λ 1 + 2 δ ϕ + λ 1 ϕ + λ 1 2 + ϕ 2 ) 4 λ 2 ϕ ( γ ( δ + 2 ϕ ) + ( δ + ϕ ) 2 ) + 4 ϕ 2 ( γ + δ + ϕ ) 2 + 4 λ 2 2 ( δ + ϕ ) 2 + Δ λ 1 + 4 λ 1 2 Ψ A 45 and, the extinction times T 1 ( 1 ) = − 4 λ 2 ( γ + 2 δ + 2 ϕ ) − 2 λ 1 ( 3 γ + δ + 3 λ 2 + 3 ϕ ) − 2 ( γ + δ + ϕ ) 2 − 6 ϕ ( γ + δ + ϕ ) − 2 λ 2 2 Ξ + 6 λ 2 ( δ + 2 ϕ ) ( γ + δ + ϕ ) + 6 ϕ 2 ( γ + δ + ϕ ) + 6 ϕ ( γ + δ + ϕ ) 2 + 6 λ 1 2 ( γ + λ 2 + ϕ ) + 6 λ 2 2 ( δ + ϕ ) + λ 1 Υ Θ , A 46 T 1 ( 2 ) = − 2 λ 2 ( γ + 3 δ + 3 ϕ ) − λ 1 ( 8 γ + 4 δ + 6 λ 2 + 8 ϕ ) − 2 ( γ + δ + ϕ ) 2 − 6 ϕ ( γ + δ + ϕ ) − 2 λ 1 2 A + 6 λ 2 ( δ + 2 ϕ ) ( γ + δ + ϕ ) + 6 ϕ 2 ( γ + δ + ϕ ) + 6 ϕ ( γ + δ + ϕ ) 2 + 6 λ 1 2 ( γ + λ 2 + ϕ ) + 6 λ 2 2 ( δ + ϕ ) + λ 1 Υ Θ , A 47 T 2 ( 1 ) = − 3 ( γ + δ + λ 2 + ϕ ) Δ + 2 B ( 6 λ 2 ( δ + 2 ϕ ) ( γ + δ + ϕ ) + 6 ϕ 2 ( γ + δ + ϕ ) + 6 ϕ ( γ + δ + ϕ ) 2 + 6 λ 1 2 ( γ + λ 2 + ϕ ) + 6 λ 2 2 ( δ + ϕ ) + λ 1 Υ ) Θ A 48 and T 2 ( 2 ) = − 3 ( γ + δ + λ 1 + ϕ ) 2 C + 6 λ 2 ( δ + 2 ϕ ) ( γ + δ + ϕ ) + 6 ϕ 2 ( γ + δ + ϕ ) + 6 ϕ ( γ + δ + ϕ ) 2 + 6 λ 1 2 ( γ + λ 2 + ϕ ) + 6 λ 2 2 ( δ + ϕ ) + λ 1 Υ Θ , A 49 where parameters Ψ, Δ, Λ, Θ, Υ , Ξ , Ω, A, B and C are given by: Λ = γ 2 + γ δ + 2 γ λ 2 + 2 γ ϕ + δ λ 2 + 2 δ ϕ + λ 2 ϕ + λ 2 2 + ϕ 2 , Δ = 4 ϕ ( γ 2 + γ δ + 2 γ ϕ + 2 δ ϕ + ϕ 2 ) + 4 λ 2 ( 2 γ δ + 2 γ ϕ + 2 δ ϕ + ϕ 2 ) + 8 δ λ 2 2 + 4 λ 2 2 ϕ , Ψ = γ 2 + 2 γ λ 2 + 2 γ ϕ + λ 2 ϕ + λ 2 2 + ϕ 2 , Ω = γ 2 + γ δ + γ λ 2 + 2 γ ϕ + 2 δ λ 2 + 2 δ ϕ + λ 2 ϕ + ϕ 2 , Θ = 4 λ 2 ϕ ( γ ( δ + 2 ϕ ) + ( δ + ϕ ) 2 ) + 4 ϕ 2 ( γ + δ + ϕ ) 2 + 4 λ 2 2 ( δ + ϕ ) 2 + Δ λ 1 + 4 λ 1 2 Ψ , Ξ = 4 λ 2 ( γ δ + 2 γ ϕ + δ 2 + 2 δ ϕ + ϕ 2 ) + 4 ϕ ( γ + δ + ϕ ) 2 + 4 Γ λ 1 + 4 δ λ 2 2 + 4 λ 2 2 ϕ , Υ = 12 λ 2 ( γ + δ + ϕ ) + 6 ( γ + 2 ϕ ) ( γ + δ + ϕ ) + 6 λ 2 2 , A = 4 λ 2 ( γ δ + 2 γ ϕ + δ 2 + 2 δ ϕ + ϕ 2 ) + 4 ϕ ( γ + δ + ϕ ) 2 + 4 λ 1 2 ( γ + λ 2 + ϕ ) + 4 λ 1 Ω , B = γ 2 + 2 γ δ + 2 γ λ 2 + 2 γ ϕ + δ 2 + δ λ 1 + δ λ 2 + 2 δ ϕ + λ 2 ϕ + λ 2 2 + ϕ 2 and C = γ 2 + 2 γ δ + γ λ 1 + γ λ 2 + 2 γ ϕ + δ 2 + 2 δ λ 1 + 2 δ ϕ + λ 1 ϕ + λ 1 2 + ϕ 2 . A 50 Figure 6. Schematic of the fluctuating growth-rate model for N = 3. Performing the Laplace transform, we obtainSolving this system of four equations and four unknowns yieldsand. Expanding these functions in terms ofyields the extinction probabilities,and, the extinction timeswhere parameters, Δ, Λ,andare given by:

Footnotes