Slow protein fluctuations explain the emergence of growth phenotypes and persistence in clonal bacterial populations

Abstract

One of the most challenging problems in microbiology is to understand how a small fraction of microbes that resists killing by antibiotics can emerge in a population of genetically identical cells, the phenomenon known as persistence or drug tolerance. Its characteristic signature is the biphasic kill curve, whereby microbes exposed to a bactericidal agent are initially killed very rapidly but then much more slowly. Here we relate this problem to the more general problem of understanding the emergence of distinct growth phenotypes in clonal populations. We address the problem mathematically by adopting the framework of the phenomenon of so-called weak ergodicity breaking, well known in dynamical physical systems, which we extend to the biological context. We show analytically and by direct stochastic simulations that distinct growth phenotypes can emerge as a consequence of slow-down of stochastic fluctuations in the expression of a gene controlling growth rate. In the regime of fast gene transcription, the system is ergodic, the growth rate distribution is unimodal, and accounts for one phenotype only. In contrast, at slow transcription and fast translation, weakly non-ergodic components emerge, the population distribution of growth rates becomes bimodal, and two distinct growth phenotypes are identified. When coupled to the well-established growth rate dependence of antibiotic killing, this model describes the observed fast and slow killing phases, and reproduces much of the phenomenology of bacterial persistence. The model has major implications for efforts to develop control strategies for persistent infections.

Figure 1. The epigenetic landscape.

The probability distribution (upper panel) and the corresponding epigenetic landscape (lower panel) are shown for a 1-dimensional system, characterized by the expression levels of one single gene. The bimodal structure of the probability distribution of the cell states induces a dual landscape, similar to the potential energy landscape in Hamiltonian systems, in which the modes of the probability distribution are mapped into metastable states. The profile of the landscape is a manifestation of the gene dynamics, which the cell can explore at equilibrium, driven by stochastic fluctuations. For fluctuations small enough (light blue arrows), the cell remains confined in the basin of attraction around each metastable state, while strong enough fluctuations (red arrows) will make the cell hop from one basin of attraction to an adjacent one. If these transitions are rare, and the sojourn time within a basin of attraction is comparable with the observational time, we identify all possible states within that basin of attraction as one single (noisy) phenotype.

Figure 2. Protein, division time, and growth rate distributions in the ergodic regime.

Parameter values are k₁ = 3⋅10⁻², k₂ = 0.35, γ₁ = 0.04, γ₂ = 4⋅10⁻³ (all units in sec⁻¹), T₀ = 2100 (sec), κ = 0.01 (nM) ⁻¹, V₀ = 1.7 fl. Histograms are the result of direct stochastic simulations (see text for details). The full curve in panel (a) corresponds to the Gamma distribution (8) with parameters a = 7.14 as from (9) and b = k₂ ln2/γ₁ = 6.06, and no other fitting parameters. The ln2 here comes from averaging the cellular volume over the cell cycle , . Panel (b) shows the histogram for the division time distribution, obtained by direct measurement of T_div during the simulation. The histogram for the growth rate distribution (panel (c)) is obtained from the measured T_div by computing μ = ln2/T_div. In the ergodic regime, the growth rate distribution is characterized by one mode only, corresponding to one single phenotype.

Figure 3. Protein, division time, and growth rate distributions in the slow fluctuations case.

Parameter values for the simulation (a,b,c) are k₁ = 1.6⋅10⁻⁵, k₂ = 1.0, γ₁ = 0.01, γ₂ = 4⋅10⁻⁵ (all units in sec⁻¹), T₀ = 2100 (sec), κ = 0.01 (nM) ⁻¹, V₀ = 1.7 fl. The full curve in panel (a) corresponds to the Gamma distribution (8) with parameters a = 0.045 and b = 96.78 fitted as described in the text. The full curves in panels (b) and (c) correspond respectively to the division time distribution, Eq. (12), and growth rate distribution, Eq. (13), evaluated with the same parameters. Parameter values for the simulation (d,e,f) are k₁ = 1.6⋅10⁻⁴, k₂ = 1.0, γ₁ = 0.01, γ₂ = 4⋅10⁻⁵ (all units in sec⁻¹), T₀ = 2100 (sec), κ = 0.01 (nM) ⁻¹, V₀ = 1.7 fl. The full curve in panel (d) corresponds again to the Gamma distribution (8) with parameters a = 0.5 and b = 112.68 fitted as described in the text, and the full curves in panels (e) and (f) correspond to Eqs. (12) and (13), evaluated with these same parameters. The good agreement between direct simulations and the predictions (12) and (13) supports the validity of the slow fluctuations approximation, leading to Eq. (10). The peak on the right in the growth rate distribution corresponds to the majority of cells growing at the maximal growth rate.

Figure 4. Protein, division time, and growth rate distributions in the slow fluctuations case and fast translation.

Parameter values for the simulation are k₁ = 1.0⋅10⁻⁴, k₂ = 5.0, γ₁ = 0.01, γ₂ = 4⋅10⁻⁵ (all units in sec⁻¹), T₀ = 2100 (sec), κ = 0.01 (nM) ⁻¹, V₀ = 1.7 fl. The full curve in panel (a) corresponds to the Gamma distribution (8) with parameters a = 0.69 and b = 579.8 fitted as described in the text. The full curves in panels (b) and (c) correspond to Eqs. (12) (division times distribution) and (13) (growth rate distribution) respectively, evaluated with the same parameters. In this parameter regime, weakly non-ergodic components dominate the dynamics. The second peak on the left of the growth rate distribution represents a minority of cells growing at slow growth rate, while the peak on the right corresponds instead to the majority of cells growing at the maximal growth rate.

Figure 5. Killing curves showing the phenomenon of persistence.

Killing curves result from direct Gillespie simulations and are here compared with the static disorder approximation, Eq. (16). (a) Ergodic regime. Parameters are the same as in Fig. 1, with k₀ = 5. In this case, no biphasic behaviour is apparent. The dashed red line corresponds to the slope of a single exponential killing with and 43.2 nM. Panel (b) shows different parameter sets characterizing the weakly non-ergodic regime. Parameter sets for the simulation were: (Red) k₁ = 2.0⋅10⁻⁷, k₂ = 1.0, γ₁ = 0.01, γ₂ = 4⋅10⁻⁵ (all units in sec⁻¹), T₀ = 2100 (sec), κ = 1.0 (nM) ⁻¹, V₀ = 1.7 fl, k₀ = 5; (Blue) k₁ = 1.0⋅10⁻⁶, k₂ = 1.0, γ₁ = 0.01, γ₂ = 4⋅10⁻⁵ (all units in sec⁻¹), T₀ = 2100 (sec), κ = 1.0 (nM) ⁻¹, V₀ = 1.7 fl, k₀ = 5; (Black) k₁ = 1.0⋅10⁻⁶, k₂ = 10.0, γ₁ = 0.01, γ₂ = 4⋅10⁻⁵ (all units in sec⁻¹), T₀ = 2100 (sec), κ = 1.0 (nM) ⁻¹, V₀ = 1.7 fl, k₀ = 5. The values for the mean number of bursts a, and for the mean burst size b were fitted from the corresponding protein distribution and used to evaluate the static disorder approximation (16), indicated with dashed lines. The corresponding values of a and b are reported in the legend box for ease of reading. Notice the regularly spaced jolts, more apparent during the fast killing phase, corresponding to the majority of cells dividing at regular intervals T₀. The biphasic behaviour of the killing curve depends qualitatively on both a and b. The lower the mean number of bursts a, the longer the initial killing phase, and the smaller the persister population, while the larger the mean burst size b, the flatter the persister tail. In general, within the present model persistence requires small a’s and large b’s.

Figure 6. Contour plot showing the discriminant given by Eq. (19).

In this plot κ = 0.01 (nM)⁻¹ has been assumed. Weakly non-ergodic behaviour, characterized by the emergence of two modes in the growth rate distribution, is predicted for a and b values such that is positive.

Figure 7. Behaviour of the volume growth law.

Plot of Eq. (26) for and (blue curve). Parameters were litres,  = 0.01 (nM) ⁻¹, (sec⁻¹), and the protein copy number was rescaled to with (Avogadro number). The dashed red line represents the slope associated with the exponential asymptotic growth law , shown for comparison. The shift between the two curves is due to the arbitrary prefactor in front of the exponential. The Lambert function introduces a deviation from exponential volume growth at short times.