Phenotypic-dependent variability and the emergence of tolerance in bacterial
populations

José Camacho Mateu1, ¶, Matteo Sireci2, ¶, Miguel A. Muñoz2, *,

1 Departamento de Matemáticas, Universidad Carlos III de Madrid, Leganés,
Spain

2 Departamento de Electromagnetismo y Física de la Materia and

Instituto Carlos I de Física Teórica y Computacional. Universidad de Granada.
Granada, Spain

¶ These authors contributed equally to this work.

*<EMAIL_ADDRESS>
## Abstract

Ecological and evolutionary dynamics have been historically regarded as
unfolding at broadly separated timescales. However, these two types of
processes are nowadays well-documented to intersperse much more tightly than
traditionally assumed, especially in communities of microorganisms. Advancing
the development of mathematical and computational approaches to shed novel
light onto eco-evolutionary problems is a challenge of utmost relevance. With
this motivation in mind, here we scrutinize recent experimental results
showing evidence of rapid evolution of tolerance by lag in bacterial
populations that are periodically exposed to antibiotic stress in laboratory
conditions. In particular, the distribution of single-cell lag times —i.e.,
the times that individual bacteria from the community remain in a dormant
state to cope with stress— evolves its average value to approximately fit the
antibiotic-exposure time. Moreover, the distribution develops right-skewed
heavy tails, revealing the presence of individuals with anomalously large lag
times. Here, we develop a parsimonious individual-based model mimicking the
actual demographic processes of the experimental setup. Individuals are
characterized by a single phenotypic trait: their intrinsic lag time, which is
transmitted with variation to the progeny. The model —in a version in which
the amplitude of phenotypic variations grows with the parent’s lag time— is
able to reproduce quite well the key empirical observations. Furthermore, we
develop a general mathematical framework allowing us to describe with good
accuracy the properties of the stochastic model by means of a macroscopic
equation, which generalizes the Crow-Kimura equation in population genetics.
Even if the model does not account for all the biological mechanisms (e.g.,
genetic changes) in a detailed way —i.e., it is a phenomenological one— it
sheds light onto the eco-evolutionary dynamics of the problem and can be
helpful to design strategies to hinder the emergence of tolerance in bacterial
communities. From a broader perspective, this work represents a benchmark for
the mathematical framework designed to tackle much more general eco-
evolutionary problems, thus paving the road to further research avenues.

## Author summary

Problems in which ecological and evolutionary changes occur at similar
timescales and feedback into each other are ubiquitous and of outmost
importance, especially in microbiology. A particularly relevant problem is
that of the emergence of tolerance to antibiotics by lag, that has been
recently shown to emerge very fast in bacterial (E. coli) populations under
controlled laboratory conditions. Here, we present a computational individual-
based model, allowing us to reproduce empirical observations and, also,
introduce a very general analytical framework to rationalize such results. We
believe that our combined computational and analytical approach may inform the
development of well-informed strategies to mitigate the emergence of bacterial
tolerance and resistance to antibiotics and, more generally, can help shedding
light onto more general eco-evolutionary problems.

## Introduction

The extraordinary ability of species to adapt and survive in unpredictably-
changing and unfavorable environments is certainly one of the most astonishing
features among the many wonders of the phenomenon that we call life. Such
adaptations can occur at extremely fast temporal scales thus interspersing
ecological and evolutionary processes [1, 2, 3]. A widely spread surviving
strategy is _latency_ or _dormancy_ , i.e., the possibility for organisms to
enter a period of reduced metabolic activity and non-replication adopted
during adverse environmental conditions [4, 5, 6, 7, 8]. Examples of dormancy
can be found across kingdoms, with examples ranging from microorganisms such
as viruses, bacteria or fungi [9, 10, 11, 12] to plants [13, 14] and animals
[7]. During the latency period the organism is said to be in a latent or
dormant state and the time it takes to wake up is referred to as “lag time” or
simply “lag”. Entering and exiting a dormant state are not cost-free
processes, since individuals may require of a specific metabolic machinery for
performing such transitions and/or the development of specifically-devised
“resting structures” [15, 16, 4, 5, 17]. The exit from the dormant state can
occur either as a response to environmental signals or cues [4, 5, 18, 10] or,
alternatively, in a stochastic way [19, 20, 21, 22, 23]. As a matter of fact,
the duration of the lag intervals often varies widely between conspecific
individuals and even between genetically identical organisms exposed to the
very same environmental conditions [24, 10, 25, 19]. Such a variability is
retained as an example of _phenotypic diversification_ or bet-hedging strategy
[26, 27] that confers a crucial competitive advantage in unpredictable and
rapidly changing environments, thus compensating the above-mentioned
individual costs and providing important benefits to the community as a whole
[4, 5, 10, 28, 29, 14].

Although, as already stated, latency is a widespread phenomenon, bacterial
communities constitute the most suitable playground for quantitative analysis
of latency owing to their diversity, fast life cycle, and the well-controlled
conditions in which they can grow and proliferate in the laboratory [30, 31,
32, 33]. Actually, latency was first described by Müller back in $1895$ as an
explanation for the observed irregularities in the growth rate of bacterial
cultures in his laboratory [34]. In recent years it has been realized that
bacterial latency is a more complex and rich phenomenon than previously
thought. Indeed, paraphrasing a recent review on the subject, the lag phase is
“dynamic, organized, adaptive, and evolvable” [10].

Bacterial latency is at the root of _tolerance_ to antibiotics as, rather
often, bactericidal antibiotics act during the reproduction stage and thus, by
entering a dormant state, bacteria become transiently insensitive to
antibiotics. Let us recall that bacterial _tolerance_ is not to be confused
with bacterial _resistance_ [35]. While _resistance_ refers to the ability of
organisms to grow within a medium with antibiotics, provided these are not in
high concentrations, _tolerance_ is the ability to transiently overcome
antibiotics, even at very high concentrations, provided the exposition time is
not too large [36, 24, 35]. The strengths of these two complementary surviving
strategies are quantified, respectively, in terms of quantities: (i) the
minimum inhibitory concentration (MIC) of drug that must be supplied to stop
the population growth —a quantity that is significantly increased in resistant
strains [35, 37, 25, 24]— and (ii) the minimum duration to kill $99\%$ of the
cells $MDK_{99}$, which is increased in tolerant strains [38].

While the importance of bacterial resistance has long been recognized, studies
underlining the crucial role played by tolerance are less frequent and more
recent [35, 37, 25, 24]. An important caveat is that, while resistance is
specific to one or a few antibiotics, tolerance is generically effective for a
large diversity of them, leading to survival even under intensive multidrug
treatment [35, 25, 24]. Moreover, there exists firm evidence that tolerance is
the first response to antibiotic stress [37], facilitating the later
appearance of resistance [25]. Therefore, understanding the emergence of
tolerance is crucial for the development of more effective therapies aimed at
dealing with recalcitrant infections and possibly preventing them. Aimed at
shedding light on these issues, here we present an eco-evolutionary approach
to analyze the emergence of tolerance by lag in bacterial communities under
controlled laboratory experiments. In particular, we scrutinize the conditions
under which modified lag-time distributions evolve as a response to stressful
environments and investigate the origin of the experimentally-observed broad
heavy tails in lag-time distributions (see below).

Beside this specific focus, the present work has a broader breath. The example
of rapid evolution of lag-time distributions is used as a test to prove a
theoretical framework that we are presently developing. Our framework is
similar in spirit to existing approaches such as the theory of _“adaptive
dynamics”_ and related models in population genetics [39, 40, 41, 42], but
aims at reconciling and generalizing them.

As a historical sidenote, let us recall that _adaptive dynamics_ (AD) was born
as a generalization of evolutionary game theory [43] to allow for a set of
strategies that is continuously varying and, upon which selection acts. AD led
to the satisfactory explanation of intriguing phenomena such as evolutionary
branching [39, 40, 44, 45], speciation [46, 47, 48], diversification [49, 50],
the emergence of altruism and cooperation [51, 52], and the evolution of
dispersal [53]. Importantly, its foundations are also mathematically well-
established [54]. However, in spite of its very successful history, AD in its
standard formulation has some limitations that make it not directly applicable
to complex situations such as the one we aim at describing here:

(i) First of all, in its standard formulation, populations are considered as
monomorphic, i.e. point-like in phenotypic space; thus it does not allow for
phenotypically-structured populations (see however [55, 56]).

(ii) The “macroscopic equations” of AD for the populations are not easily
connected to microscopic birth-death processes in individual-based models
[57].

(iii) Variations are assumed to be small, typically Gaussian-distributed and
independent of the parent’s phenotypic state.

(iv) Variations are considered to be rare: “after every mutational event, the
ecological dynamics has time to equilibrate and reach a new ecological
attractor” [58]. In other words, a separation is assumed between ecological
and evolutionary timescales, while in microbial communities, such processes
may occur in concomitance. Such a convergence of characteristic timescales is
the hallmark of eco-evolutionary dynamics [1, 2, 3] and is at the basis of
fascinating phenomena such as eco-evolutionary tunneling [58, 59] and other
rapid evolutionary phenomena [60, 61, 32, 62, 63] which are difficult to
account for in the standard formulation of adaptive dynamics.

In what follows, we employ a theoretical framework in the spirit of
statistical mechanics that aims to fill the gap between theory,
phenomenological models and, most importantly, experiments. Thus, our approach
—which is similar in spirit to previous work on bacterial quorum-sensing by E.
Frey and collaborators [64, 65]— implements a number of extensions with
respect to standard AD [66], as it makes explicit the connection between
individuals (microscale), and community dynamics (macroscale), introduces a
general variation kernel, allows for large and phenotypic-state-dependent
variations, etc. These extensions allow us to study phenotypic diversity
within a well-characterized eco-evolutionary framework. A full account of this
general theoretical framework will be presented elsewhere [67]. Let us finally
emphasize that many of the above questions and extensions have been already
tackled in the mathematical literature, at a formal level [54, 68]. Yet, to
the best of our knowledge, these results have have been confined to rigorous
analyses of toy models and have not fully percolated through the biological
and physical literature.

The paper is organized as follows: in the first section we discuss in detail
the experimental setup and empirical findings object of our study; then, we
introduce a stochastic individual-based model implementing phenotypic
variability and inheritability to account for experimental results. We present
an extensive set of both computational and analytical results for it,
discussing in particular the conditions under which the mathematical results
deviate from computational ones. Finally, we discuss the implications of our
work both from a biological viewpoint and how it contributes to the
understanding of the evolution of heterogeneous phenotypic distributions, as
well as from a more general eco-evolutionary perspective.

### Empirical observations: rapid evolution of lag-time distributions

For the sake of concreteness, we focus on recent experimental results on the
rapid evolution of tolerance in populations of Escherichia coli in laboratory
batch cultures in Balaban’s lab [24]. In particular, a bacterial population is
periodically exposed to antibiotics (amplicillin) in very high concentrations
(much larger than the MIC) during a fixed-duration time interval $T_{a}$
(e.g., $T_{a}=3,5,$ or $8$ hours). After antibiotic exposure the system is
washed and the surviving population is regrown in a fresh medium during a time
interval $T$ (with $T=23h-T_{a}$). The antibiotics/fresh-medium cycle is
iterated at least $8$ or $10$ times. Results are averaged over $2$
experimental realizations for each $T_{a}$ and the resulting maximal carrying
capacity is about $10^{9}$ individuals (we refer to [24] for further
biological and experimental details).

Once the cycles are completed, Fridman et al. [24] isolated some individuals
from the surviving community and by regrowing them in a fresh medium they
found that the distribution $P(\tau)$ of lag times $\tau$ —i.e. the time
individual dormant cells take to start generating a new colony after
innoculation into a fresh medium— changes from its ancestral shape to a
modified one, shifted towards larger $\tau$ values. More specifically, the
mean value grew to a value that approximately matches the duration of the
antibiotic-exposure time interval, $T_{a}$ (see [24]). This modified lag-time
distribution entails an increase in the survival probability under exposure to
ampicillin but, also, to antibiotics of a different bactericidal class such as
norfloxacin, for the same time period. Furthermore, mutations were identified
in diverse genes, some of them known to be related with regulatory circuits
controlling the lag-time distribution, such as the toxin–antitoxin one [69].
Subsequently, after many cycles, the population was also observed to develop
resistance to ampicillin [24]. Thus, the conclusion is that non-specific
tolerance —stemming from lag— emerges in a very rapid way as a first adaptive
change/response to antibiotic stress. More in general, these results reveal
that the adaptive process is so fast that ecological and evolutionary
processes occur at comparable timescales [63, 70, 71, 72, 73].

The experimentally-determined lag-time distributions reveal another intriguing
aspect that —to the best of our knowledge— has not been extensively analyzed
so far: their variance is also significantly increased as $T_{a}$ grows and,
related to this, the resulting mean value of the distribution is always larger
than its median [74]. This is an indication that, as a matter of fact, the
empirically-obtained lag-time distributions are skewed and exhibit heavy
tails, including phenotypes with anomalously-large lag times —much larger than
$T_{a}$—, especially for large $T_{a}$’s. This observation is surprising as,
under such controlled lab conditions, one could naively expect to find lag-
time distributions sharply peaked around the optimal time value, $T_{a}$,
since, ideally, the best possible strategy would be to wake-up right after
antibiotics are removed and any further delay comes at the price of a
reduction of the overall growth rate or fitness. Fridman  et al. proposed that
the increase in the variance might suggest a past selection for a bet-hedging
strategy in natural unpredictable environments; however, anomalously-large
lag-time values were not present in the original wild-type population. The
authors also suggested that there could be constraints at the molecular level
imposing the mean and the variance of the lag-time distribution to increase
concomitantly [24, 69], a possibility that inspired us and that we will
carefully scrutinize from a theoretical and computational perspective in what
follows.

## Model building

Fig 1: Sketch of the main ingredients of the individual-based stochastic
model. Each individual bacterium ($i$) is characterized by its phenotypic
state, lag time $\tau_{i}$ and experiences demographic processes. (A) In the
presence of antibiotics, bacteria can stochastically switch between the
dormant and the growing state (at transition rates $s$ and $1/\tau_{i}$,
respectively); growing individuals can also attempt reproduction (at a “birth”
rate $b$) and be immediately killed by the action of antibiotics (as
bactericidal antibiotics usually act during duplication attempts). (B) In the
fresh medium, dormant bacteria can wake up at a rate $1/\tau_{i}$, that
depends on their intrinsic (phenotypic) lag time; on the other hand, growing
cells can reproduce asexually by duplication; the resulting offspring inherit
the characteristic time scale with some variation, as specified by a function
$\beta$. (C) Two possible types of variation functions $\beta$: in the
additive case (top), the standard deviation is constant, i.e. independent of
the initial state $\tau_{i}$, while in the multiplicative case (bottom) the
standard deviation is assumed to grow linearly with the parent’s lag time
$\tau_{i}$. (D) Sketch of the environmental variation, alternating
periodically between antibiotic exposure (time $T_{a}$) and a fresh medium
($T_{max}-T_{a}$).

Aimed at shedding light onto these empirical findings, here we propose an
individual-based stochastic model for phenotypic adaptation in which each
single individual cell can be either in an “awake” or in a “dormant” state
[75, 22, 23, 76, 77] (see Fig 1 for a sketch of the model). Mimicking the
experimental protocol of Fridman et al.— a population of such individuals is
exposed to alternating adverse and favorable conditions with durations $T_{a}$
and $23h-T_{a}$, respectively (a function $\eta(t)$ labels the environmental
state at any given time $t$: $\eta(t)=-1$ in the presence of antibiotics and
$\eta(t)=+1$ in the fresh medium).

The model assumes that each awake cell is able to sense the environment and
respond to it by regulating its state: they can sense the presence of
antibiotics and enter the dormant state at rate $s$, while such a machinery is
assumed to be turned off during dormancy. In S1 Text, Sec. S5, we also
consider a generalization of the model in which awake individuals can also
enter the dormant state as a response to other sources of stress such as
starvation [78]. Indeed, the wake-up is assumed to occur as a result of a
Markovian stochastic process; each individual bacteria $i$ is _phenotypically_
characterized by its intrinsic typical mean lag time $\tau_{i}$ meaning that,
it wakes up stochastically at a constant transition rate $1/\tau_{i}$.
Therefore, the time $t$ at which a dormant cell actually resumes growth is a
random variable distributed as $P(t|\tau_{i})=e^{-t/\tau_{i}}/\tau_{i}$, with
mean value $\tau_{i}$ [79, 80]. In the last section we discuss recent
alternatives to Markovian processes, i.e. including some form of “memory”,
which can give rise to non-exponential residence times, to describe this type
of waking-up phenomena [81, 82].

Awake individuals are exposed to stochastic demographic processes: they
attempt asexual reproduction (i.e., duplication) at a constant birth rate $b$
and die spontaneously at rate $d$ (that we fix to $0$ without loss of
generality). Reproduction attempts are successful in the fresh medium while,
in the presence of antibiotics, they just lead to the parent’s death and its
removal from the community. Following this dynamics, the population can freely
grow, until its size reaches a maximal carrying capacity $K$. Once this limit
has been reached, the population enters a saturated regime, within which each
new birth is immediately compensated by a random killing (much as in the Moran
process [83]).

Importantly, in parallel with the above demographic processes, the model
implements an evolutionary/adaptive dynamics. The phenotypic state $\tau_{i}$
of each successfully dividing individual is transmitted, with possible
variation, to its progeny. In particular, the two offspring resulting from
duplication have phenotypic states $\tau_{i}+\xi_{1}$ and $\tau_{i}+\xi_{2}$,
respectively, where $\xi_{1}$ and $\xi_{2}$ are the phenotypic stochastic
variations, sampled from some probability distribution, that we generically
call $\beta(\xi;\tau_{i})$ and that, in the more general case, can be state-
dependent, i.e. depend on $\tau_{i}$. More specifically, we implemented two
different variants of the model, depending of the standard deviation of the
probability distribution $\beta(\xi;\tau_{i})$:

  * •

The _additive_ model, with a standard deviation, $\alpha_{A}$, common to all
phenotypes.

  * •

The _multiplicative_ model, with a state-dependent standard deviation,
$\alpha_{M}\tau_{i}$, for individuals with intrinsic lag time $\tau_{i}$,
where $\alpha_{M}$ is a constant (see Methods).

Observe that in the multiplicative case, the larger the parent’s lag time the
larger the possible amplitude of variations, in a sort of rich-get-richer or
Matthew-effect mechanism, well-known in the theory of stochastic processes to
generate heavy tails [84, 85, 86, 87, 88, 89, 90, 91, 92, 93]. As a motivation
for this choice, let us mention that there is solid evidence that the genetic
circuits involved in the regulation of the lag-time distribution (such as the
toxin-antitoxin one), can indirectly produce this type of fluctuations at the
phenotypic level [69]. Furthermore, similar phenotypic-variation kernels have
been argued to arise from non-linear effects in the way genotypic changes
(mutations) are manifested into phenotypic variability (see e.g. [94, 95]).

## Analytical (mean-field) theory

Before delving into computational analyses of the model, let us present a
mathematical framework allowing us to obtain theoretical insight. Readers not
particularly interested in analytical approaches can safely skip this section,
and just be aware that it is possible to mathematically understand all the
forthcoming computational results.

The previous Markovian stochastic individual-based model is mathematically
defined as a “many-particle” Master equation ruling the time evolution of the
joint probability-distribution functions for the whole set of all “particles”
( i.e., cells). The resulting master equation can be simulated computationally
by employing the Gillespie Algorithm (see below and S1 Text, Sec. S2A, for
details) [79, 80, 96]. However, as it is often the case for such many-particle
Master equations, it is hard to handle analytically in an exact way. Thus, in
order to gain quantitative understanding beyond purely computational analyses,
here we develop an approximation —which becomes exact in the limit of
infinitely large population sizes [64, 65]— that allows us to derive a
macroscopic (or “mean-field”) description of the stochastic model in terms of
the probability density of finding an individual at any given phenotypic
state, $\tau$ (i.e. the “one-particle” probability density). The mean-field
approach that we employ in what follows is just a first example of a much more
general framework that we will expose in detail elsewhere [67].

A first step toward the derivation of a _macroscopic equation_ relies on a
marginalization of the many-particle probability-distribution function to
obtain a one-particle probability density (see S1 Text, Sec. S2B). The
resulting marginalized distribution function encapsulates the probability
density $\phi(\tau,t)$ that a randomly sampled individual at time $t$ has lag
time $\tau$. This probability —that needs to be normalized, so that
$\int^{\infty}_{0}\phi(\tau,t)d\tau=1$— can be decomposed in two contributions
$\phi(\tau,t)=\phi_{G}(\tau,t)+\phi_{D}(\tau,t)$ representing, respectively,
the relative fraction of individuals in growing (G) and dormant (D) states.
Observe that these two densities are not probability distributions and thus
they are not normalized to unity separately. In the limit of infinitely-large
population sizes, the evolution of the probability density for individuals in
the growing state, ${\phi}_{G}(\tau,t)$ is ruled by the following equation
(details of the derivation can be found in SI, S2A):

$\displaystyle\partial_{t}\phi_{G}(\tau,t)$ $\displaystyle=$
$\displaystyle\frac{1+\eta(t)}{2}\left[-b\phi_{G}(\tau,t)+2b\int^{\infty}_{0}d\tilde{\tau}\beta(\tau-\tilde{\tau};\tilde{\tau})\phi_{G}(\tilde{\tau},t)\right.$
(1) $\displaystyle-$
$\displaystyle\left.b\phi_{G}(\tau,t)\int^{\infty}_{0}d\tilde{\tau}\phi_{G}(\tilde{\tau},t)\right]$
$\displaystyle-$
$\displaystyle\frac{1-\eta(t)}{2}\left[b\Big{(}1-\int^{\infty}_{0}d\tilde{\tau}\phi_{G}(\tilde{\tau},t)\Big{)}+s\right]\phi_{G}(\tau,t)+\frac{1}{\tau}\phi_{D}(\tau,t).$

Even if this equation might look cumbersome, its different terms have a rather
intuitive interpretation:

  * •

_In the fresh medium_ (terms proportional to $1+\eta(t)$): (i) the first term
represents the negative probability flow stemming from growing individuals
with generic phenotypic trait $\tau$ that reproduce (at rate $b$) and change
to any other arbitrary phenotypic state; (ii) the second represents the
positive contribution of reproducing individuals (at rate $b$) with any
arbitrary trait $\tilde{\tau}$, for which one of the two resulting offspring
jumps to $\tau$ (controlled by the function
$\beta(\tau-\tilde{\tau};\tilde{\tau})$); (iii) the third _selection term_
stems from the normalization of the overall probability density: if the
population size grows because any individual with arbitrary trait
$\tilde{\tau}$ reproduces successfully (at rate $b$), then the relative
probability to observe phenotype $\tau$ decreases to keep the overall
probability-density conserved.

  * •

_In the presence of antibiotics_ (terms proportional to $1-\eta(t)$): (i) the
first term represents the rate at which growing individuals that attempt
reproduction (at rate $b$) are killed by antibiotics; (ii) the second term is
a selection term, fully analogous to the above-discussed one: when any
arbitrary individual dies the overall probability density at $\tau$ increases;
(iii) the third term represents the outflow of individuals entering the
dormant state at rate $s$.

  * •

_In both environments_ (no dependence on $\eta(t)$): the only term,
proportional to the rate $1/\tau$, describes the probability inflow stemming
from dormant individuals that become awake.

Similarly, the equation for the density of individuals in the dormant state is

$\displaystyle\partial_{t}\phi_{D}(\tau,t)$ $\displaystyle=$
$\displaystyle-\eta(t)b\phi_{D}(\tau,t)\int^{\infty}_{0}d\tilde{\tau}\phi_{G}(\tilde{\tau},t)-\frac{1}{\tau}\phi_{D}(\tau,t)+\frac{1-\eta(t)}{2}s\phi_{G}(\tau,t)$
(2)

where the first (selection) term stems from the overall probability
conservation when the population either grows or shrinks (negative or positive
signs, respectively), and the remaining two terms have the opposite meaning
(and signs) of their respective counterparts in Eq.(1).

In order to make further analytical progress, in the case in which variations
are assumed to be small, it is possible to introduce a further (“diffusive” or
“Kimura”) approximation as often done in population genetics as well as in
adaptive or evolutionary mathematical approaches [97]. More specifically, one
can perform a standard (Kramers-Moyal) expansion of the master equation by
assuming that jumps in the phenotypic space are relatively small [79, 80],
i.e. expanding the function $beta$ in Taylor series around $0$. After some
simple algebra (see S1 Text, Sec. S3) one obtains a particularly simple
expression for the overall probability distribution:

$\displaystyle\partial_{t}\phi(\tau,t)$ $\displaystyle=$
$\displaystyle\eta(t)\left[f(\tau,t)-\bar{f}(t)\right]\phi(\tau,t)$ (3)
$\displaystyle-$
$\displaystyle(\eta(t)+1)\left[\partial_{\tau}\theta(\tau)f(\tau,t)\phi(\tau,t)-\frac{1}{2}\partial_{\tau}^{2}\sigma^{2}(\tau)f(\tau,t)\phi(\tau,t)\right]$

where the _“ effective fitness function”_ $f(\tau,t)\equiv
b\phi_{G}(\tau,t)/\phi(\tau,t)$ and its population average
$\bar{f}(t)=\int^{\infty}_{0}d\tau f(\tau,t)\phi(\tau,t)$ have been
introduced, and where $\theta(\tau)$ and $\sigma^{2}(\tau)$ are the first and
second cumulants of the variation function $\beta$ (in first approximation we
can assume $\theta(\tau)=0$, while $\sigma^{2}(\tau)=\alpha_{A}^{2}$ for the
additive case and $\sigma^{2}(\tau)=\alpha_{M}^{2}\tau^{2}$ for the
multiplicative case). Observe that, remarkably, this last equation is a
generalization of the celebrated continuous-time _Crow-Kimura equation_ of
population genetics [98], also called _selection-mutation equation_ [99, 100].
In particular, notice that the dynamics of the probability density is exposed
to the combined action of the process of selection (first term in Eq.(3),
which is nothing but the _replicator equation_ [43, 101]) and mutation, as
specified by the drifts in the second line. This type of equations, combining
replicator dynamics with Fokker-Planck type of terms —even if with a slightly
different interpretation— have been also studied by Sato $\&$ Kaneko and Mora
$\&$ Walzak [102, 103, 104]. The main —and crucial— differences between Eq.(3)
and the standard Crow-Kimura equation are:

  * •

The fitness function appears in the mutation terms —whereas in the standard
Crow-Kimura equation the diffusion term would read
$\partial^{2}_{\tau}\phi(\tau)$— thus correlating reproduction rates and
mutation amplitudes. Observe that here variations are always associated with
reproduction events, as typically in bacteria and viruses, in such a way that
a higher fitness rate implies a higher mutation rate.

  * •

There is a general dependence on the cumulants of the variation kernel that,
in general, can be trait-dependent and asymmetric.

These generalizations are essential ingredients to capture the essence of our
Markovian model as we will see and, to the best of our knowledge, have not
been carefully analyzed in the past. From here on, we refer to Eq.(3) as the
_generalized Crow-Kimura (GCK) equation_.

## Results

Fig 2: Tuning the only free parameter to match empirical results.(A-B) Mean of
the lag-time distribution $K_{1}$, measured right at the end of each
antibiotic cycle, as a function of the number of cycles. Computational results
are shown for antibiotic duration $T_{a}=3h$ for both the additive (A) and the
multiplicative (B) versions of the model. The different curves (color coded)
correspond to different mutation amplitudes $\alpha_{A}$ (in A) and
$\alpha_{M}$ (in B), respectively. The dashed vertical lines indicate the
$10th$ cycle, when experiments stop. Remarkably, the mean lag time strongly
depends on the mutational amplitude, both in the transient regime and in the
asymptotic state. We implement an algorithmic search to tune the only free
parameter (either $\alpha_{A}$ or $\alpha_{M}$) to best fit the experimental
mean lag times for all values of $T_{a}$ together and, in particular, their
experimentally-reported linear dependence on the antibiotic exposure time
$T_{a}$ (see Methods). (C) Mean of the lag-time distribution as a function of
$T_{a}$ for the model (squares for additive and triangles for multiplicative
versions of the model) tuned to reproduce experimental values (yellow
symbols). While empirical data are available for $T_{a}=3,5$ and $8h$, the
model can be analyzed for generic values of $T_{a}$. The solid line indicates
the linear dependence between the mean and lag-time distribution,
$K_{1}=T_{a}$, while the horizontal dashed line represents the mean lag time
of the ancestral population. Parameter values: $K=10^{5}$, $\alpha_{A}=0.16h$,
$\alpha_{M}=0.048$, $b=2.4h^{-1}$, $d=3.6\cdot 10^{-5}h^{-1}$, $s=0.12\
h^{-1}$, $T_{fresh}=23h-T_{a}$, $10$ cycles (see Methods). Fig 3: Lag-time
probability distributions: theory and simulations. (A-B) Lag-time distribution
after $10$ cycles as obtained in the simulation of the individual-based model
in both the additive (A) and the multiplicative (B) case, for different
antibiotic-exposure periods, $T_{a}=3,5$ and $8h$ (marked with different
colours). Solid and dotted vertical lines indicate, respectively, the median
and mean of the corresponding distributions (a large separation between these
two indicators reflects asymmetries in the distribution such as the emergence
of a heavy tail to the right). Dashed lines represent results from the
numerical integration of the GCK equation, Eq.3, using the same parameters and
external conditions. Observe that the multiplicative model generates much
larger tails, reproducing the experimental phenomenology better than the
additive one. (C) Initial lag-time distribution mimicking the experimentally
observed one for the ancestral population. (D) Variance, $K_{2}$, and
difference between mean and median, $K_{1}-median$, of the lag-time
distribution as a function of $T_{a}$ in the additive (blue squared symbols)
and multiplicative (red triangular symbols) versions of the model. $K_{2}$
grows with the antibiotic exposure time in the multiplicative case, while in
the additive case it remains nearly constant. The difference between the mean
and the median is very small in the additive case, while it increases with
$T_{a}$ almost monotonically in the multiplicative one. In summary, the
multiplicative model generates a distribution with a variance that grows with
the mean, as well as heavy tails, reproducing well the key experimental
findings. Parameter values are as in Fig 2.

In order to scrutinize whether the proposed adaptive stochastic model can
account for the key empirical findings of Fridman _et al._ [24], we perform
both (i) extensive computational simulations and (ii) numerical studies of the
mean-field macroscopic equation, Eq.(3).

  * •

Computational simulations rely on the Gillespie algorithm [96], which allows
us to simulate exactly the master equation defining the stochastic model. In
all cases, we consider at least $10^{3}$ independent realizations to derive
statistically-robust results. Without loss of generality and owing to
computational costs, the maximal population size or carrying capacity is fixed
to $K=10^{5}$.

  * •

On the other hand, for analytical approaches, in spite of the relatively
simple form of Eq.(3) owing to its non-linear nature and to the time-
variability of environmental conditions $\eta(t)$, it is not possible to solve
it analytically in a closed way and, thus, it becomes mandatory to resort to
numerical-integration schemes. In particular, from this equation —or, more
precisely, from integration of its two additive components: Eq.(1) and Eq.(2)—
one can derive the time-dependent as well as the asymptotic lag-time
distributions and, from them, monitor the leading moments or cumulants as a
function of time.

Further details of both computational simulations and numerical integration of
the macroscopic equation can be found in the Methods section as well as in the
S1B. In what follows we present together both types of analyses, underlining
where the mean-field approach works well and where its predictions deviate
from direct simulations.

### Transient dynamics: determining variational amplitudes

Parameter values in the model are fixed to agree as much as possible with the
empirical ones measured by Fridman et al. [24] (see Methods). In particular,
we used (i) the same set of environmental-period durations $T_{a}$ and
$T_{max}-T_{a}$ as in the antibiotics/fresh-medium cycle, (ii) the
experimentally measured reproduction rate in the fresh medium, (iii) the
empirical “falling-asleep” rate $s$, as well as (iv) the same number of
antibiotic cycles (ten) as in the experimental setup. Initially all
individuals are assumed to have small intrinsic lag-time values of $\tau$. In
particular, we consider a truncated normal distribution with mean value and
variance as in the actual ancestral population in the experiments
($\langle\tau\rangle^{exp.}=\ 1.0\ \pm 0.2\ h$).

Employing this set of experimentally-constrained parameter values and initial
conditions, we ran stochastic simulations in which the whole population
expanded and then shrank following the periodically alternating environments.
Along this dynamical cyclic process the distribution of $\tau$ values across
the population varies in time; in particular, we monitored the histogram of
$\tau$ values and obtained the corresponding probability distributions right
at the end of each antibiotic cycle, just before regrowth, as in the
experiments.

Fig 2A and 2B show the evolution of the mean (i.e. the first cumulant,
$K_{1}$) across cycles, while Fig 3A and 3B illustrate the full distribution
and higher-order cumulants after $10$ cycles. Observe that the value of
$K_{1}$ after $10$ cycles depends on the choice made for the only remaining
free parameter, i.e. the variation-amplitude parameter $\alpha_{A}$ or
$\alpha_{M}$, for the additive or multiplicative versions of the model,
respectively. In order to tune either of these parameters, we imposed that
$K_{1}(10)$ reproduces in the closest way the experimentally determined
values, as measured right before the $10th$ regrowth cycle. This tuning
procedure leads to $\alpha_{A}=0.16(1)h$ and $\alpha_{M}=0.048(1)$ for the
additive and multiplicative cases, respectively (parentheses indicate
uncertainty in the last digit); these are the values that best reproduce the
empirical findings in the sense of least-square deviation from the available
empirical data for different $T_{a}$’s (see Fig 2C).

Let us remark that both variants of the model are able to reproduce the key
experimental feature of generating mean lag times close to $T_{a}$ (observe,
however, that there is always a small deviation in the case $T_{a}=8~{}h$, for
which even experimentally, $K_{1}\approx 10~{}h>T_{a}$). Nevertheless, as
illustrated in Fig 3A and 3B there are significant differences between the two
variants. In particular, the additive model fails to reproduce the following
empirical observations:

  * •

$T_{a}$-dependent variances,

  * •

large differences between median and mean values, and

  * •

strongly skewed distributions with large tails.

For instance, in the experiments, for $T_{a}=8~{}h$, the difference between
the mean and the median is $1.1(1)~{}h$ while in the additive model is
$0.15~{}h$, i.e. about one order of magnitude smaller. Furthermore, in the
experiments, lag times of up to $30~{}h$ are observed, while in the additive
model values above $\approx 15$ are exponentially cancelled; i.e. they have an
extremely low (negligible) probability to be observed. This is also
illustrated in Fig 3D where the second and third cumulants (variance and
skewness) of the distribution after $10$ cycles are plotted as a function of
$T_{a}$. Observe that both cumulants remain almost constant, revealing the
absence of heavy tails for large values of $T_{a}$.

On the other hand, the multiplicative model is able to reproduce not only the
experimental values of the mean but also —with no additional parameter nor
fine tuning— (i) the existence of large lag-time variances that increase with
$T_{a}$, (ii) the above-mentioned large differences between the mean and the
median ($1.3(1)$ in this case), as well as (iii) heavily skewed lag-time
distributions that strongly resemble the empirically measured ones (see Fig 2
in Fridman _et al._ [24]). In particular, lag times of the order of $30~{}h$
have a non-negligible probability to be observed for $T_{a}=8~{}h$ within the
multiplicative version of the model, after $10$ cycles. The resulting
probability and the corresponding cumulants (see Fig 3D depend strongly on
$T_{a}$.

Importantly, the previous results are quite robust against changes in the
model. In particular, if growing cells are allowed to switch to dormancy in
response of starvation, the mean lag time increases, as expected, but the
qualitative shape of the lag-time distribution remains unchanged (see S6 ).
Hence, just by modifying accordingly the parameter $\alpha$ allows one to
recover the same conclusions.

As a word of caution let us emphasize that the distributions in Fig 3 are not
obtained exactly in the same way as the experimental ones. The first are
distributions of characteristic times $\tau$ (inverse of intrinsic transition
rates) while the second are the actual lag times $t$ measured after regrown in
a fresh medium. Actually, the characteristic time $\tau$, in our model, is
just a proxy for the actual time that it takes for the colony formed by such
an individual to be observable or detectable in actual experimental setups.
Below we discuss this issue more extensively as well as the possible
limitations it implies and extensions of the modelling approach to circumvent
them.

Let us also underline that Fig 3 reports not only the results of direct
simulations but also the theoretical predictions (dashed lines) derived from
numerical integration of the macroscopic equations for the two different
cases. The agreement with simulation results is remarkably accurate; the
origin of the existing small discrepancies will be analyzed in detail in a
forthcoming section.

Thus, the main conclusion of these computational and theoretical analyses is
that _state-dependent (multiplicative) variability is needed in order to
account for the empirically observed key features of the lag-time
distributions emerging after a few antibiotic/fresh-medium cycles._ Once this
variant of the model is chosen, a good agreement with experimental findings if
obtained by fitting the only free parameter: the amplitude of variations.

Fig 4: Characterization of the asymptotic state in the multiplicative version
of the model . (A) Approach to the dynamic asymptotic state for the
multiplicative case, as resulting from the integration of Eq.(1) and Eq.(2)
for $T_{a}=3h$ ($t$ indicates overall time as measured in hours). The
different curves correspond to three different values of the variation
amplitude (from bottom to top: $\alpha_{M}=0.0048,0.01,0.048$). The difference
between this plot and Fig 2 is that $K_{1}$ is measured at different times
within the cycle and not just right at the end of antibiotic exposure (see
(B)). The vertical dashed line marks the $10th$ cycle at which the experiment
stopped. Observe that the steady-state mean value, the oscillations amplitude,
and the relaxation time depend on the variation amplitude $\alpha_{M}$. (B)
Mean lag time within a single cycle ($T\in[0,23]h.$) in a asymptotic state.
During the killing phase (antibiotic exposure), i.e. for $t<T_{a}$, the mean
lag time increases to maximize the number of dormant individuals; then, in the
fresh medium the mean relaxes back to the initial value. (C) Lag-time
probability distribution —as derived from theory (dashed lines) and
computationally (solid lines)— at the start of the cycle (leftmost curve) and
when antibiotics are removed (rightmost curve); in the asymptotic state the
system oscillates between these two limiting probability distributions, both
of them exhibiting heavy tails. (D) Evolution of the three first cumulants
$K_{1}$, $K_{2}$, $K_{3}$ (mean, variance, and skewness, respectively) within
a asymptotic cycle (both theoretical and computational results are shown).
Observe in (C) and (D) that the theory correctly predicts the properties of
the distribution but there are some small errors due to finite size effects.
Fig 5: Population dynamics .  Abundances $N_{D}$ and $N_{G}$ of dormant (A)
and growing (B) cells, respectively, along a full cycle in the asymptotic
regime (reached after only three cycles) for the multiplicative model (the
curves are the result of averaging over many independent realizations for
different $T_{a}$ as color coded; observe the semi-logarithmic scale). Dormant
cells abundances reach a maximum value after approximately one hour of
exposition to antibiotics, almost independently of $T_{a}$, and then start a
slow decrease, while $N_{G}$ exhibits an opposite trend: it rapidly decreases
and reaches a minimum at $T=T_{a}$ (see inset), after which it grows
exponentially fast until the carrying capacity is reached. (C) The total
number of cells $N=N_{G}+N_{D}$ is plotted along the cycle: for all values of
$T_{a}$ the absolute minimum is reached near $T_{a}$ (as clearly seen in the
inset). (D) The fraction of dormant cells relative to the total number is
maximal nearby $T_{a}$ and decreases when antibiotics are removed.

### Asymptotic state

Even if experimental results are available for a fixed and limited number
($10$) of antibiotic-exposition cycles, the already-calibrated model allows us
to scrutinize the possible emergence of asymptotic states after a much-larger
number of cycles. In other words, it is possible to go beyond the experimental
limits and analyze the fate of the population. In this sense, the experimental
results can be seen as a “transient adaptation” to the environment, while the
evolutionary cycle would be completed only when an asymptotic (evolutionary
stable) state is reached. Let us remark that the asymptotic state is
necessarily a _periodic_ one, as the phenotypic distributions vary at
different instants of the cycle, i.e. the asymptotic distributions —measured
at arbitrary times within the cycle— exhibit periodic oscillations in its
shape, tracking the perpetual environmental cyclic changes. This is
illustrated in Fig 4, showing results obtained by numerically integrating the
macroscopic equations, Eq.(1) and Eq.(2). First of all, it shows periodic
oscillations of the mean lag time $K_{1}$; as shown in panel (A) it first
increases from its initial value $K_{1}=1$ and then, eventually, reaches an
oscillatory steady state. More specifically, as clearly seen in the zoomed
plot of panel (B), within the steady state, the maximum mean value within each
cycle is reached right before antibiotics removal. This is an expected result
as in the first part of the cycle, i.e. during the “killing phase”, the
presence of antibiotics induces a selective pressure towards increasing the
mean lag-time value because delaying the exit from the lag phase provides
protection from the antibiotics. On the other hand, in the fresh medium
(growing phase) the selective pressure quickly reduces the mean lag time to
foster fast growth and increased fitness. Thus, summing up, the periodic
alternation of environmental conditions induces a stable periodic change in
the mean lag-time value.

Actually, it is not only the mean that changes periodically, but the whole
probability distribution that varies cyclically. This is illustrated in Fig 4C
and 4D which shows computational and theoretical results for the lag-time
probability distribution and its first cumulants, $K_{1},K_{2}$ and $K_{3}$,
for the multiplicative case (similar plots for the additive case are shown in
S1 Text, Sec. S6). Observe, in particular, in panel C, that the distribution
oscillates between two extreme or limiting cases corresponding to the times of
antibiotics inoculation and antibiotics removal, respectively. This effect can
be more vividly seen in the S1 and S2 Videos.

Let us also highlight that the probability distributions exhibit non-Gaussian
tails and are right-skewed. In particular, to make these observations more
quantitative, Fig 4D shows the variance, $K_{2}$, and the skewness, $K_{3}$,
along the cycle in the steady state. Notice also the very-good —though not
perfect— agreement between computational results and theoretical estimates
(dashed lines in Fig 4C and 4D).

Furthermore, let us emphasize that, importantly, the amplitude of the
variations —as controlled by the parameter $\alpha_{M}$ (or, similarly,
$\alpha_{A}$ for the additive case)— has a non-trivial effect on both the
transient and the asymptotic behavior. In particular, the value of such
amplitude not only affects the mean value of lag times after $10$ cycles —as
illustrated by the plateau of the oscillations in Fig 4A and 4B— but also (i)
its asymptotic value, i.e. the mean lag time, (ii) the amplitude of the
oscillations across a cycle in the steady state, and (iii) the relaxation time
to the asymptotic state (i.e., the speed of evolution). This is due to the
heavy tails of the distribution: increasing the amplitude of variations
directly increments the variance of lag times, but this also enlarges the
left-skewness of the distribution, feeding-back to the mean value. Therefore,
the eco-evolutionary attractor is shaped both by selection and mutation,
departing from the classical evolutionary scenario, as e.g. in adaptive-
dynamics, in which the amplitude of the variations just affects the variance
of the resulting distribution but not the overall attractor.

Finally, we complement our observations with the population structure
dynamics, i.e. proportion, minimum and maximum of dormant and awake cell
numbers, in Fig 5. In particular, panels (A) and (B) show the abundances of
dormant and awake cells as function of time along an asymptotic cycle. Observe
that the number of dormant cells reaches a maximum, $N_{D,max}$, after one
hour, independently of the antibiotic duration time $T_{a}$, while its height
is proportional to this parameter. On the other hand, the position of minimum
of growing cells number, $N_{G,min}$, scales with $T_{a}$ and its magnitude
decreases correspondently. In Fig 5 we also show and discuss the dependence
the total number of cells $N=N_{G}+N_{D}$ (panel C) as well as the relative
fraction of dormant individuals along a full cycle in the stationary state
(which is reached after onley a few (three) antibiottic cycles).

For the sake of completeness, let us also emphasize that both versions of the
model are able to generate $MDK_{99}$ values that grow as a function of the
number of antibotic cycles, converging to an asymptotic-state value; at the
end of the tenth cycle simulations compare well with empirical observations
for different values of $T_{a}$ (see S10 and S11 Figs; observe that the
largest difference appears for $T_{a}=8$, a case for which also $K_{1}$
deviates slightly from $T_{a}$ in the experiments).

### Deviations between theory and simulations: finite-size effects

Thus far, we have reported results stemming from computational analyses of the
individual based model as well as from numerical integration of the associated
macroscopic theory, i.e. the GCK equation. Small but systematic discrepancies
between theory and simulations are evident, see for example Fig 4C and 4D. Let
us here discuss the origin of such differences.

The theoretical approach relies on two different approximations: (i) on the
one hand it considers the small-variation approximation to include just the
first two moments of the variation function (i.e. a diffusion approximation);
(ii) on the other hand, in order to derive the macroscopic GCK equation, one
needs to neglect correlations between individuals, a type of mean-field
approximation that, as usual, is expected to be exact only in the _infinite
population-size limit_ [79, 80]. In S1 Text, Sec. S4, we show computational
evidence that the small mutation approximation is not a significant source of
errors; hence, the discrepancies necessarily stem from finite-size effects.
Indeed, in the present experimental set up, there is a bottleneck at the end
of each antibiotics cycle, when there is a small number of surviving
individuals, thus limiting the validity of the mean-field approximation in
such a regime. As a matter of fact, one can clearly see from Fig 4D that the
largest discrepancies appear around the end of the killing phase, when the
population is the smallest. Note also that the main features of the dynamics
in phenotypic space are reproduction and variation: i.e., offspring are
similar to their progeny. But reproduction events occur only within the awake
(growing) sub-population; the full population-size, involving also dormant
ones, is not the most relevant quantity to gauge finite-size effects.
Therefore, in order to minimize the discrepancies between theory and
simulations it does not suffice to consider larger population sizes: even for
huge values of the carrying capacity $K$, we find that the population at the
end of the killing phase is always rather small and, hence, exposed to large
demographic fluctuations, i.e. to finite-size effects.

To put these observations on more quantitative bases, we define a parameter
$\delta$ as the deviation between theoretical and computational results for
the mean lag-time value after antibiotic exposure and monitor its dependence
on the minimal size of the awake population (i.e. right at the end of the
antibiotic phase). Fig 6 illustrates that: (A) the deviation grows with the
antibiotic-exposure time $T_{a}$, whereas (B) the minimum awake subpopulation
size decreases with $T_{a}$. Combining these two pieces of information one can
see (C) that the deviation parameter $\delta$ decreases as the minimum
subpopulation of awake individuals increases. Unfortunately, the convergence
to zero of this last curve is very slow, and thus, it is computationally very
expensive to remove finite-size effects.

Finally, let us remark that we leave for future work the formulation of an
extension of the mathematical theory accounting for finite-size effects [105,
106, 107], including corrections to the GCK equation.

Fig 6: Analysis of the deviations between simulations and theory. The
parameter $\delta$ is defined as the difference between the mean lag times
(right at the end of the antibiotic cycle) in the theoretical approach and in
computer simulations. (A) Double-logarithmic plot reporting the dependence of
the error parameter $\delta$ on the antibiotic time exposure $T_{a}$ for both
the additive (blue dots) and the multiplicative (red dots) versions of the
model; in either case, the larger the exposure time the larger the error. (B)
Minimum number of awake individual during the cycle in the asymptotic regime,
$N_{G,min}$, as function of $T_{a}$ in double-logarithmic scale. As expected,
the larger the exposure time the smaller the number of surviving individuals.
(C) Combining the data from (A) and (B) it follows that $\delta$ decreases
with increasing $N_{G,min}$, meaning that deviations between theory (expected
to be exact for infinitely large population sizes) and computational results
stem from finite-population-size effects. Notice that errors are smaller in
the multiplicative version of the model.

## Conclusions, discussion and perspectives

Summary of results and conclusions. We have presented a mathematical and
computational model to quantitatively analyze the emergence and evolution of
tolerance by lag in bacteria. Our first goal was to reproduce the main results
reported in the laboratory experiments of Fridman _et al._ [24] in which the
authors found a very fast evolution of tolerance by lag in a community of
_Escherichia coli_ bacteria periodically exposed to an antibiotics/fresh-
medium cycle. In particular, after a relatively small number of such cycles,
there is a clear change in the individual-cell lag-time distribution with its
mean value evolving to match the duration of antibiotic exposure. This is
remarkable, and demonstrates that tolerance by lag is the first and generic
strategy adopted by bacteria to survive under harsh environmental conditions
such as the presence of antibiotics. A second key empirical finding is that
concomitantly with the evolution of the mean lag time, also the variance of
lag times is significantly increased for longer antibiotic-exposure periods:
i.e. the harsher the conditions the more diversified the lag times within the
population. More generally, the full lag-time distribution becomes wider and
develops a heavy tail for sufficiently large times. This means that there
exist individual phenotypes that are clearly sub-optimal under the strictly
controlled laboratory conditions and most-likely reflects a bet-hedging
strategy, preparing the community to survive under even harsher conditions
(i.e. longer stressful periods).

To shed light onto these observations we developed a stochastic individual-
based model assuming that individuals are characterized by an intrinsic lag
time, setting the “typical” time at which such individual stochastically wakes
up after dormancy. This phenotypic trait is transmitted to the progeny with
possible variation. By considering a protocol analogous to the experimental
one (i.e. alternating antibiotic exposure and fresh medium growth) the model
is able to produce a distribution of characteristic lag times across the
population that reproduces quite well the empirical results in all cases by
tuning a single parameter value. In particular, the emerging lag-time
distributions have a mean that matches the period of antibiotic exposure
$T_{a}$, an increase of the mean and variance with $T_{a}$, as well as a large
difference between the mean and the median, which result from the appearance
of heavy tails in the lag-time distributions. Nevertheless, it is important to
underline that the distributions that the model generates are just a proxy for
the empirically-determined ones, where the actual times in which individual
bacteria give rise to new and detectable (i.e. visible with the available
technology) colonies are measured.

Importantly, in order to account for all the above empirical phenomenology,
the model needs to assume _multiplicative variations_ , i.e. that the
variability between the parent’s trait and those of its offspring increases
(linearly) with the parent’s lag time: the larger the parent’s lag time the
larger the possible variation. This multiplicative process — at the roots of
the emerging heavy tails in the lag-time distribution— resembles the so-called
_rich-get-richer_ mechanism of the _Matthew effect_ [84, 86, 87, 91, 85]. This
type of variations implements an effective dependence between the parent’s
trait value and the variation amplitude, that was hypothesized as a possible
mechanism behind the experimental results and that could stem from a highly
non-linear map between genotypic changes and their phenotypic manifestations
[24, 69].

Notably, our analyses reveal that the amplitude of variations affects not just
the variance ($K_{2}$) of the resulting lag-time distribution, but also its
mean ($K_{1}$) as well as other higher-order cumulants such as the skewness
($K_{3}$). This is in blatant contrast with standard approaches to
evolutionary or adaptive dynamics, in which the “mutational amplitude” only
influences the “broadness” ($K_{2}$) of the distribution of traits in
phenotypic space, but does not alter the attractor of the dynamics (e.g.
$K_{1}$). Thus, the introduction of state-dependent (multiplicative)
variability constitutes a step forward into our understanding of how simple
adaptive/evolutionary processes can generate complex outcomes.

Let us finally mention that our model describes rapid evolution, where
ecological and mutational time scales are comparable. This interplay between
ecological and evolutionary processes is explicit in the asymptotic state: it
is not an “evolutionary stable state” but a _“non-equilibrium evolutionary
stable state.”_ By non-equilibrium we mean that the detailed-balance condition
—a requirement of equilibrium states [108]— is violated and thus, there are
net probability fluxes in phenotypic space. These correspond to adaptive
oscillations in phenotypic space. Key properties of such a state (oscillation
plateaus, amplitudes, etc.) depend on the mutational amplitude, i.e., the
amplitude of variations determine the eco-evolutionary attractor. In future
work we will scrutinize much in depth non-equilibrium characteristic
properties, such as non-vanishing entropy-production of these type of complex
eco-evolutionary processes [109, 23, 110].

Advantages and limitations of the phenotypic-modeling approach. As already
underlined, the present model assumes adaptation at a phenotypic level. Is
this a biologically realistic assumption? The answer to this question, in
principle, is affirmative but some caveats are in order.

First of all, let us recall that a large part of the theoretical work on
evolutionary dynamics and adaptation developed during the last decades focuses
on phenotypic adaptation. For instance, in the theory of adaptive dynamics,
individuals are always characterized by some phenotypic trait or set of traits
which is subject to selection and transmission to the progeny with variation
[39, 40, 41, 46, 47] (see also e.g. [71, 72]). In general, this is the most
parsimonious way of modeling adaptation as the details of the genotypic-
phenotypic mapping are usually highly non-linear or simply unknown (see e.g.
[95, 111, 112, 113, 114]).

On the one hand, adaptation beyond genetic changes —for example epigenetic
adaptation— is a well-documented phenomenon in the bacterial world [115] and
is the focus of intense research activity [116, 117, 27, 118, 119]. For
instance, recent work explores “the evolutionary advantage of heritable
phenotypic heterogeneity”, which suggests that evolutionary mechanisms at a
phenotypic level, such as the ones employed in our approach, might be
biologically favored with respect to more-standard genetic mechanisms, under
certain circumstances [120]. In particular, such phenotypic variability can
provide a faster and more flexible type of response than the one associated
with traditional genetic mutations.

Nevertheless, it is important to underline that Fridman _et al._ found
empirical evidence that —in their specific setup— genetic mutations were
always present in the evolved strains. In particular, they found mutations in
genes controlling the so-called toxin-antitoxin circuit, mediating the
response to antibiotic stress [24]. This regulatory circuit is known to lead
to “multiplicative fluctuations” in the lag-time distribution at the
phenotypic level [69]. Thus, strictly speaking, our modeling approach
constitutes an effective or phenomenological approximation to the more complex
biology of this problem.

This observation opens promising and exciting avenues for future research to
shed light on how broad probability distributions of lag times —possibly with
heavy tails— can be actually encoded in phenotypic or genetic models.
Actually, scale-free (power-law) distributions of bacterial lag times have
been recently reported in a specifically-devised experimental setup [121].
Similarly to our conclusions, this work also emphasizes that a broad
distribution of individual-cell waking-up rates is needed to generate non-
exponential decays of the overall lag-time distribution.

Similarly, another exciting possibility would be to develop computational
models akin to the phenotypic one proposed here but implementing genetic
circuitry; i.e. models where the phenotype is the (possibly stochastic)
outcome of an underlying regulated genetic process and where the object of
selection are not specific lag times but their whole distributions as
genetically encoded.

Future developments and perspectives. In future research, we would like to
further delve onto several aspects, both biological and theoretical, of the
present work. As a first step, we leave for forthcoming work the analysis of
the pertinent question of how similar systems respond to randomly fluctuating
environments as opposed to periodically changing ones; do they develop heavier
tails to cope with such unpredictable conditions in a sort of bet-hedging
strategy? How do the statistical features of the environmental variability
translate into the emerging lag-time distributions? [22, 27, 82, 122, 28,
123].

From a more theoretical perspective, we leave for an impending work the
formulation of an extension of our approach that fully accounts for finite-
size effects, thus introducing the next-to-leading order corrections to the
generalized-Crow-Kimura macroscopic equation accounting for demographic
fluctuations. Within this context, treating the variation-amplitude itself as
an evolving trait is also a potentially fruitful route for further studies.

Finally, as a long-term project we plan to develop models and analytical
approaches, similar to the ones explored here, but focusing on genetic
evolution, employing explicit genotypic-phenotypic mappings, rather than just
on phenotypic changes. In particular, by introducing this further layer of
complexity it would be possible to generate more general types of single-cell
lag-time distributions, not limited to exponential ones as the purely
Markovian approach considered here. Let us recall that a more general
stochastic non-Markovian framework —i.e., including memory effects (see e.g.
[81, 82, 124])— is a challenging goal that promises to be very pertinent and
relevant for many diverse problems in which the control of time is important.

## Methods

### Numerical values of the parameters

In order to fix parameter values we employed the experimental values and
measurements in [24] as closely as possible. The number of bacteria involved
in the experiment reaches values of the order of $\sim 10^{9}$; however this
number is prohibitively large for computer simulations and we fixed a maximum
carrying capacity of $K=10^{5}$, verifying that results do not depend strongly
on such a choice (see finite-size effects section). Initially the number of
cells in the growing state is fixed to be equal to the carrying capacity $K$;
thus no cell is initially in the dormant state). The doubling time of both the
ancestral and the evolved populations is $25\pm 0.3\min$; thus on average
every single bacteria attempts reproduction at a rate $b=1/25\min=2.4\
h^{-1}$. The death rate for (natural) causes (i.e. other than antibiotics) is
$d=3.6\cdot 10^{-5}\ h^{-1}$. The awakening rate is given by the inverse of
the characteristic time $a=1/\tau$ [75]. The initial condition (ancestral or
wild population) was randomly sampled from a truncated Gaussian peaked at
$\tau=0$. Since the empirical ancestral distribution is narrow and close to
the origin [24] (mean lag time $\langle\tau\rangle^{exp.}_{0}=\ 1\ \pm 0.2\
h$) we fix the standard deviation of the truncated Gaussian distribution to
$\sigma=1\ h\ 16\ \min$ in such a way that $\langle\tau\rangle^{sim.}_{0}\sim
1\ h$. Neither the exit rate from dormancy $s$ nor the amplitude of the
mutations, $\alpha_{A}$ and $\alpha_{M}$, can be experimentally measured, but
we can fix them indirectly (the rest of the parameters are kept fixed with the
values specified above). First, $s$ can be chosen using the experimental
information that for the ancestral population $MDK_{99}\sim 2.55h$. Hence, we
leave the (simulated) ancestral population in the antibiotic phase until the
$99\%$ becomes extinct; averaging over different initial conditions, we found
that $s=0.12\ h^{-1}$ is a good approximation. To fix the constants
$\alpha_{A}$ and $\alpha_{M}$ we performed simulations for diverse values of
such parameters and looked for those that best reproduce the experimental
tendency after $10$ exposure cycles for the different exposure times under
consideration (in particular, we performed a least-square deviation analysis
to match the straight-line $\langle\tau\rangle=T_{a}$ when performing a linear
interpolation for all $T_{a}$’s). A systematic sweep of the values of the
first two significant digits led us to $\alpha_{A}=0.16h$ and
$\alpha_{M}=0.048$.

### Variation Functions

We consider two different variation kernels for lag-time variations
$\delta=\tau-\tilde{\tau}$: the additive one, $\beta_{A}(\delta;\tilde{\tau})$
and multiplicative one, $\beta_{M}(\delta;\tilde{\tau})$. Both of them are
probability density functions of $\delta$, normalized in
$[-\tilde{\tau},\infty]$ and may depend on the initial phenotype
$\tilde{\tau}$. In particular, the additive case reads:
$\beta_{A}(\delta;\tilde{\tau})=e^{-\frac{\delta^{2}}{2\alpha_{A}^{2}}}/Z_{A}(\tilde{\tau})$
with
$Z_{A}(\tau)=\alpha_{A}\frac{\sqrt{\pi}}{\sqrt{2}}Erfc(-\frac{\tau}{\sqrt{2}\alpha_{A}})$
(where $Erfc$ stands for the complementary error function), while in the
multiplicative case, we consider
$\beta_{M}(\delta;\tilde{\tau})=e^{-\frac{\delta^{2}}{2\alpha_{M}\tilde{\tau}}}/Z_{M}(\tilde{\tau})$
with
$Z_{M}(\tau)=\alpha_{M}\tau\sqrt{\frac{\pi}{2}}Erfc(-\frac{1}{\sqrt{2}\alpha_{M}})$.

### Measuring lag-time distributions

In order to determine lag-time distributions, we computed histograms in
phenotypic space, as discretized in bins of size $\Delta\tau=10^{-2}$ and
averaged over many realizations of the process. In the asymptotic steady
state, similar histograms were computed at different times along the
antibiotic/fresh-medium cycle (e.g. right after antibiotic inoculation or
after antibiotic removal). To obtain results for the transient state we
determined the histogram after running for $10$ cycles. On the other hand, to
determine the steady state, we started measuring after $300$ cycles (to make
sure that a steady state has been reached) and then collect statistics up to
cycle $1500$, at intervals of $10$ cycles to avoid correlations. We repeated
the process for $30$ realizations and calculated the histogram as well as the
associated cumulants.

### Numerical Integration of the macroscopic equation

The parameter set and initial condition for numerical integration of the mean-
field equations are the same as specified above. Numerical integration was
carried out using the finite differences method. In particular $\partial_{t}$
was approximated using first order forward differences, $\partial^{2}_{\tau}$
using second-order centered differences, and integrals were approximated as
Riemann sums [125]. The numerical integration steps used in the figures are:
$h_{t}=10^{-6}$ and $h_{\tau}=10^{-2}$ . Note that, when we calculate the
probability distributions during the simulation, we must use the same bin size
to be able to correctly compare with the theoretical distributions later. We
used absorbing boundary conditions $\phi(\tau=0,t)=\phi(\tau=\tau_{max},t)=0$,
where $\tau_{max}$ is the limit of the phenotype space considered for the
numerical integration, in particular: $\tau_{max}=60h$.

Acknowledgments: We are grateful to N. Balaban for sharing data with us as
well as for her very valuable insights and comments. We also thank R. Rubio de
Las Casas, S. Manrubia, G. B. Morales, J.A. Bonachela, J. Grilli, and V.
Buendía for stimulating and fruitful discussions.

## References

  * 1. Hendry AP.  Eco-evolutionary dynamics.  Princeton university press; 2016.
  * 2. Schoener TW.  The Newest Synthesis: Understanding the Interplay of Evolutionary and Ecological Dynamics.  Science. 2011;331(6016):426–429.
  * 3. Pelletier F, Garant D, Hendry AP.  Eco-evolutionary dynamics.  Philosophical Transactions of the Royal Society B: Biological Sciences. 2009;364(1523):1483–1489.
  * 4. Lennon JT, Jones SE.  Microbial seed banks: the ecological and evolutionary implications of dormancy.  Nature reviews microbiology. 2011;9(2):119–130.
  * 5. Shoemaker WR, Lennon JT.  Evolution with a seed bank: the population genetic consequences of microbial dormancy.  Evolutionary applications. 2018;11(1):60–75.
  * 6. Fraser D, Kaern M.  A chance at survival: gene expression noise and phenotypic diversification strategies.  Molecular microbiology. 2009;71(6):1333–1340.
  * 7. Guppy M, Withers P.  Metabolic depression in animals: physiological perspectives and biochemical generalizations.  Biological Reviews. 1999;74(1):1–40.
  * 8. Sturm A, Dworkin J.  Phenotypic diversity as a mechanism to exit cellular dormancy.  Current Biology. 2015;25(17):2272–2277.
  * 9. Barton ES, White DW, Cathelyn JS, Brett-McClellan KA, Engle M, Diamond MS, et al.  Herpesvirus latency confers symbiotic protection from bacterial infection.  Nature. 2007;447(7142):326–329.
  * 10. Bertrand RL.  Lag phase is a dynamic, organized, adaptive, and evolvable period that prepares bacteria for cell division.  Journal of bacteriology. 2019;201(7):e00697–18.
  * 11. van Vliet S.  Bacterial dormancy: how to decide when to wake up.  Current Biology. 2015;25(17):R753–R755.
  * 12. Rittershaus ES, Baek SH, Sassetti CM.  The normalcy of dormancy: common themes in microbial quiescence.  Cell host & microbe. 2013;13(6):643–651.
  * 13. Childs DZ, Metcalf C, Rees M.  Evolutionary bet-hedging in the real world: empirical evidence and challenges revealed by plants.  Proceedings of the Royal Society B: Biological Sciences. 2010;277(1697):3055–3064.
  * 14. Gremer JR, Venable DL.  Bet hedging in desert winter annual plants: optimal germination strategies in a variable environment.  Ecology Letters. 2014;17(3):380–387.
  * 15. Rechinger KB, Siegumfeldt H, Svendsen I, Jakobsen M.  Early protein synthesis of Lactobacillus delbrueckii ssp. bulgaricus in milk revealed by methionine labeling and two-dimensional gel electrophoresis.  Electrophoresis: An International Journal. 2000;21(13):2660–2669.
  * 16. Larsen N, Boye M, Siegumfeldt H, Jakobsen M.  Differential expression of proteins and genes in the lag phase of Lactococcus lactis subsp. lactis grown in synthetic medium and reconstituted skim milk.  Appl Environ Microbiol. 2006;72(2):1173–1179.
  * 17. Van Bodegom P.  Microbial maintenance: a critical review on its quantification.  Microbial ecology. 2007;53(4):513–523.
  * 18. Paredes-Sabja D, Setlow P, Sarker MR.  Germination of spores of Bacillales and Clostridiales species: mechanisms and proteins involved.  Trends in microbiology. 2011;19(2):85–94.
  * 19. Wright ES, Vetsigian KH.  Stochastic exits from dormancy give rise to heavy-tailed distributions of descendants in bacterial populations.  Molecular Ecology. 2019;28(17):3915–3928.
  * 20. Epstein SS.  Microbial awakenings.  Nature. 2009;457(7233):1083–1083.
  * 21. Buerger S, Spoering A, Gavrish E, Leslin C, Ling L, Epstein S.  Microbial scout hypothesis, stochastic exit from dormancy, and the nature of slow growers.  Appl Environ Microbiol. 2012;78(9):3221–3228.
  * 22. Kussell E, Leibler S.  Phenotypic Diversity, Population Growth, and Information in Fluctuating Environments.  Science. 2005;309(5743):2075–2078.
  * 23. Leibler S, Kussell E.  Individual histories and selection in heterogeneous populations.  Proceedings of the National Academy of Sciences. 2010;107(29):13183–13188.
  * 24. Fridman O, Goldberg A, Ronin I, Shoresh N, Balaban NQ.  Optimization of lag time underlies antibiotic tolerance in evolved bacterial populations.  Nature. 2014;513:418–421.
  * 25. Levin-Reisman I, Ronin I, Gefen O, Braniss I, Shoresh N, Balaban NQ.  Antibiotic tolerance facilitates the evolution of resistance.  Science. 2017;355(6327):826–830.
  * 26. Philippi T, Seger J.  Hedging one’s evolutionary bets, revisited.  Trends in Ecology and Evolution. 1989;4(2):41 – 44.
  * 27. Veening JW, Smits WK, Kuipers OP.  Bistability, Epigenetics, and Bet-Hedging in Bacteria.  Annual Review of Microbiology. 2008;62(1):193–210.
  * 28. Villa Martín P, Muñoz MA, Pigolotti S.  Bet-hedging strategies in expanding populations.  PLoS computational biology. 2019;15(4):e1006529.
  * 29. Xu Y, Vetsigian K.  Phenotypic variability and community interactions of germinating Streptomyces spores.  Scientific reports. 2017;7(1):1–13.
  * 30. Kussell E.  Evolution in Microbes.  Annual Review of Biophysics. 2013;42(1):493–514.
  * 31. Lenski RE.  Experimental evolution and the dynamics of adaptation and genome evolution in microbial populations.  The ISME Journal: Multidisciplinary Journal of Microbial Ecology. 2017;11(10):2181–2194.
  * 32. Good B, McDonald M, Barrick J, Lenski R, Desai M.  The dynamics of molecular evolution over 60,000 generations.  Nature. 2017;551(7678):45–50.
  * 33. Elena SF, Lenski RE.  Evolution experiments with microorganisms: the dynamics and genetic bases of adaptation.  Nature Reviews Genetics. 2003;4(6):457–469.
  * 34. Müller M.  Ueber den Einfluss von Fiebertemperaturen auf die Wachsthumsgeschwindigkeit und die Virulenz des Typhus.  Zeitschrift für Hygiene und Infektionskrankheiten, medizinische Mikrobiologie, Immunologie und Virologie. 1895;20:245.
  * 35. Brauner A, Fridman O, Gefen O, Balaban NQ.  Distinguishing between resistance, tolerance and persistence to antibiotic treatment.  Nature Reviews Microbiology. 2016;14(5):320.
  * 36. Kohanski MA, Dwyer DJ, Collins JJ.  How antibiotics kill bacteria: from targets to networks.  Nature Reviews Microbiology. 2010;8(6):423–435.
  * 37. Li B, Qiu Y, Shi H, Yin H.  The importance of lag time extension in determining bacterial resistance to antibiotics.  Analyst. 2016;141(10):3059–3067.
  * 38. Brauner A, Shoresh N, Fridman O, Balaban NQ.  An experimental framework for quantifying bacterial tolerance.  Biophysical journal. 2017;112(12):2664–2671.
  * 39. Geritz SAH, Metz JAJ, Kisdi E, Meszéna G.  Dynamics of Adaptation and Evolutionary Branching.  Phys Rev Lett. 1997;78:2024–2027.
  * 40. Geritz SAH, Metz JAJ, Kisdi E, Meszéna G.  Evolutionarily singular strategies and the adaptive growth and branching of the evolutionary tree.  Evolutionary Ecology. 1998;12:35–37.
  * 41. O D.  Diekmann O. 2004. A beginner’s guide to adaptive dynamics school.  Banach Center Publications; 2004.
  * 42. Dieckmann U, Metz JA.  Surprising evolutionary predictions from enhanced ecological realism.  Theoretical Population Biology. 2006;69(3):263–281.
  * 43. Hofbauer J, Sigmund K.  Evolutionary game dynamics.  Bulletin of the American Mathematical Society. 2003;40(4):479–519.
  * 44. Dieckmann U, Law R.  The dynamical theory of coevolution: a derivation from stochastic ecological processes.  Journal of mathematical biology. 1996;34(5-6):579–612.
  * 45. Doebeli M, Dieckmann U.  Evolutionary branching and sympatric speciation caused by different types of ecological interactions.  The american naturalist. 2000;156(S4):S77–S101.
  * 46. Dieckmann DM Ulf.  On the origin of species by sympatric speciation.  Nature. 1999;400:354–357.
  * 47. Doebeli M, Dieckmann U.  Speciation along environmental gradients.  Nature. 2003;421(6920):259–264.
  * 48. Dieckmann U, Doebeli M, Metz JA, Tautz D.  Adaptive Speciation.  Cambridge Studies in Adaptive Dynamics. Cambridge University Press; 2004\.
  * 49. Spencer CC, Tyerman J, Bertrand M, Doebeli M.  Adaptation increases the likelihood of diversification in an experimental bacterial lineage.  Proceedings of the National Academy of Sciences. 2008;105(5):1585–1589.
  * 50. Friesen ML, Saxer G, Travisano M, Doebeli M.  Experimental evidence for sympartric ecological diversification due to frequency-dependent competition in escherirchia coli.  Evolution. 2004;58(2):245 – 260.
  * 51. Fletcher JA, Doebeli M.  A simple and general explanation for the evolution of altruism.  Proceedings of the Royal Society B: Biological Sciences. 2009;276(1654):13–19.
  * 52. Doebeli M, Hauert C, Killingback T.  The Evolutionary Origin of Cooperators and Defectors. 2004;306(5697):859–862.
  * 53. Doebeli M, Ruxton GD.  Evolution of disperal rates in metapopulations models: branching and cyclic dynamics in phenotype space.  Evolution;51(6):1730–1741.
  * 54. Champagnat N, Ferrière R, Méléard S.  Unifying evolutionary dynamics: From individual stochastic processes to macroscopic models.  Theoretical Population Biology. 2006;69(3):297 – 321.
  * 55. Henriques GJB, Ito K, Hauert C, Doebeli M.  On the importance of evolving phenotype distributions on evolutionary diversification.  PLOS Computational Biology. 2021;17(2):1–21.
  * 56. Rubin IN, Doebeli M.  Rethinking the evolution of specialization: A model for the evolution of phenotypic heterogeneity.  Journal of Theoretical Biology. 2017;435:248–264.
  * 57. Doebeli M, Yaroslav.  Towards a mechanistic foundation of evolutionary theory.  Elife. 2017;6(:e23804).
  * 58. Kotil K Seyfullah Enes Vetsigian.  Emergence of evolutionarily stable communities through eco-evolutionary tunnelling.  Nature Ecology and Evolution. 2018;.
  * 59. Iwasa Y, Michor F, Nowak MA.  Stochastic Tunnels in Evolutionary Dynamics.  Genetics. 2004;166(3):1571–1579.
  * 60. Hiltunen T, Virta M, Laine AL.  Antibiotic resistance in the wild: an eco-evolutionary perspective.  Philosophical Transactions of the Royal Society B: Biological Sciences. 2017;372(1712):20160039.
  * 61. Cordero OX, Polz MF.  Explaining microbial genomic diversity in light of evolutionary ecology.  Nature reviews Microbiology. 2014;12(4):263—273.
  * 62. Sanchez A, Gore J.  Feedback between Population and Evolutionary Dynamics Determines the Fate of Social Microbial Populations.  PLOS Biology. 2013;11(4):1–9.
  * 63. Callahan BJ, Fukami T, Fisher DS.  Rapid evolution of adaptive niche construction in experimental microbial populations.  Evolution. 2014;68(11):3307–3316.
  * 64. Bauer M, Knebel J, Lechner M, Pickl P, Frey E.  Ecological feedback in quorum-sensing microbial populations can induce heterogeneous production of autoinducers.  Elife. 2017;6:e25773.
  * 65. Frey E, Knebel J, Pickl P.  Mean-field equation for a stochastic many-particle model of quorum-sensing microbial populations.  arXiv preprint arXiv:180205307. 2018;.
  * 66. Geritz SA, Metz JA, Kisdi É, Meszéna G.  Dynamics of adaptation and evolutionary branching.  Physical Review Letters. 1997;78(10):2024.
  * 67. Sireci M, Muñoz MA.  Statistical physics of phenotypic evolution: Adaptive dynamics and beyond;.
  * 68. Buerger R.  The Mathematical Theory of Selection, Recombination, and Mutation.  John Wiley & Sons; 2000.
  * 69. Rotem E, Loinger A, Ronin I, Levin-Reisman I, Gabay C, Shoresh N, et al.  Regulation of phenotypic variability by a threshold-based mechanism underlies bacterial persistence.  Proceedings of the National Academy of Sciences. 2010;107(28):12541–12546.
  * 70. Messer PW, Ellner SP, Hairston Jr NG.  Can population genetics adapt to rapid evolution?  Trends in Genetics. 2016;32(7):408–418.
  * 71. Bonachela JA, Wortel MT, Stenseth NC.  Eco-evolutionary Red Queen dynamics regulate biodiversity in a metabolite-driven microbial system.  Scientific reports. 2017;7(1):1–9.
  * 72. Vallina SM, Martinez-Garcia R, Smith SL, Bonachela JA.  Models in microbial ecology.  In: Encyclopedia of Microbiology. Elsevier; 2019. p. 211–246.
  * 73. Wortel MT, Peters H, Bonachela JA, Stenseth NC.  Continual evolution through coupled fast and slow feedbacks.  Proceedings National Academy of Sciences. 2020;117(8):4234–4242.
  * 74. Levin-Reisman I, Gefen O, Fridman O, Ronin I, Shwa D, Sheftel H, et al.  Automated imaging with ScanLag reveals previously undetectable bacterial growth phenotypes.  Nature methods. 2010;7(9):737—739.
  * 75. Balaban NQ, Merrin J, Chait R, Kowalik L, Leibler S.  Bacterial Persistence as a Phenotypic Switch.  Science. 2004;305(5690):1622–1625.
  * 76. Baranyi J.  Stochastic modelling of bacterial lag phase.  International journal of food microbiology. 2002;73(2-3):203–206.
  * 77. Métris A, Le Marc Y, Elfwing A, Ballagi A, Baranyi J.  Modelling the variability of lag times and the first generation times of single cells of E. coli.  International journal of food microbiology. 2005;100(1-3):13–19.
  * 78. Moreno-Gámez S, Kiviet DJ, Vulin C, Schlegel S, Schlegel K, van Doorn GS, et al.  Wide lag time distributions break a trade-off between reproduction and survival in bacteria.  Proceedings of the National Academy of Sciences. 2020;117(31):18729–18736.  doi:10.1073/pnas.2003331117.
  * 79. Van Kampen NG.  Stochastic processes in physics and chemistry. vol. 1.  Elsevier; 1992.
  * 80. Gardiner C.  Stochastic Methods: A Handbook for the Natural and Social Sciences.  Springer Series in Synergetics. Springer; 2009.
  * 81. Norman TM, Lord ND, Paulsson J, Losick R.  Stochastic switching of cell fate in microbes.  Annual review of microbiology. 2015;69:381–403.
  * 82. Himeoka Y, Mitarai N.  When to wake up? The optimal waking-up strategies for starvation-induced persistence.  PLoS computational biology. 2021;17(2):e1008655.
  * 83. Moran PAP.  Random processes in genetics.  Mathematical Proceedings of the Cambridge Philosophical Society. 1958;54(1):60–71.  doi:10.1017/S0305004100033193.
  * 84. Sornette D.  Multiplicative processes and power laws.  Phys Rev E. 1998;57(4):4811.
  * 85. Barabási AL, Albert R.  Emergence of scaling in random networks.  Science. 1999;286(5439):509–512.
  * 86. Mitzenmacher M.  A Brief History of Generative Models for Power Law and Lognormal Distributions.  Internet Mathematics. 2002;1:226–251.
  * 87. Newman MEJ.  Power laws, Pareto distributions and Zipf’s law.  Contemp Phys. 2005;46(5):323–351.
  * 88. Manrubia SC, Zanette DH.  Stochastic multiplicative processes with reset events.  Physical Review E. 1999;59(5):4945.
  * 89. Grinstein G, Muñoz MA, Tu Y.  Phase structure of systems with multiplicative noise.  Physical Review Letters. 1996;76(23):4376.
  * 90. Tu Y, Grinstein G, Muñoz MA.  Systems with multiplicative noise: critical behavior from KPZ equation and numerics.  Physical Review Letters. 1997;78(2):274.
  * 91. Muñoz M. Multiplicative noise in non-equilibrium phase transitions: A tutorial; in Advances in Condensed Matter and Statistical Physics, Ed. E. Koroutcheva and R. Cuerno; 2004.
  * 92. Muñoz MA, Colaiori F, Castellano C.  Mean-field limit of systems with multiplicative noise.  Physical Review E. 2005;72(5):056102.
  * 93. Genovese W, Muñoz MA.  Recent results on multiplicative noise.  Physical Review E. 1999;60(1):69.
  * 94. van Nimwegen E.  Influenza Escapes Immunity Along Neutral Networks.  Science. 2006;314(5807):1884–1886.
  * 95. Koelle K, Cobey S, Grenfell B, Pascual M.  Epochal Evolution Shapes the Phylodynamics of Interpandemic Influenza A (H3N2) in Humans.  Science. 2006;314(5807):1898–1903.
  * 96. Gillespie DT.  A general method for numerically simulating the stochastic time evolution of coupled chemical reactions.  Journal of computational physics. 1976;22(4):403–434.
  * 97. Kimura M.  Diffusion models in population genetics.  Journal of Applied Probability. 1964;1(2):177–232.
  * 98. Crow JF, Kimura M, et al.  An introduction to population genetics theory.  An introduction to population genetics theory. 1970;.
  * 99. Hofbauer J.  The selection mutation equation.  Journal of mathematical biology. 1985;23(1):41–53.
  * 100. Page KM, Nowak MA.  Unifying evolutionary dynamics.  Journal of theoretical biology. 2002;219(1):93–98.
  * 101. Nowak MA.  Evolutionary dynamics: exploring the equations of life.  Harvard University Press; 2006.
  * 102. Sato K, Kaneko K.  On the distribution of state values of reproducing cells.  Physical biology. 2006;3(1):74.
  * 103. Sato K, Kaneko K.  Evolution equation of phenotype distribution: General formulation and application to error catastrophe.  Physical Review E. 2007;75(6):061909.
  * 104. Mora T, Walczak AM.  Effect of phenotypic selection on stochastic gene expression.  The Journal of Physical chemistry B. 2013;117(42):13194–13205.
  * 105. Cremer J, Melbinger A, Frey E.  Evolutionary and population dynamics: a coupled approach.  Physical Review E. 2011;84(5):051921.
  * 106. Cremer J, Melbinger A, Frey E.  Growth dynamics and the evolution of cooperation in microbial populations.  Scientific reports. 2012;2(1):1–6.
  * 107. Wakano JY, Iwasa Y.  Evolutionary Branching in a Finite Population: Deterministic Branching vs. Stochastic Branching.  Genetics. 2013;193(1):229–241.
  * 108. Gnesotto FS, Mura F, Gladrow J, Broedersz CP.  Broken detailed balance and non-equilibrium dynamics in living systems: a review.  Reports on Progress in Physics. 2018;81(6):066601.
  * 109. Mustonen V, Lässig M.  Fitness flux and ubiquity of adaptive evolution.  Proceedings of the National Academy of Sciences. 2010;107(9):4248–4253.
  * 110. Kussell E, Vucelja M.  Non-equilibrium physics and evolution—adaptation, extinction, and ecology: a Key Issues review.  Reports on Progress in Physics. 2014;77(10):102602.
  * 111. Dalziel AC, Rogers SM, Schulte PM.  Linking genotypes to phenotypes and fitness: how mechanistic biology can inform molecular ecology.  Molecular ecology. 2009;18(24):4997–5017.
  * 112. Arias CF, Catalán P, Manrubia S, Cuesta JA.  toyLIFE: a computational framework to study the multi-level organisation of the genotype-phenotype map.  Scientific reports. 2014;4(1):1–10.
  * 113. Manrubia S, Cuesta JA.  Distribution of genotype network sizes in sequence-to-structure genotype–phenotype maps.  Journal of The Royal Society Interface. 2017;14(129):20160976.
  * 114. Manrubia S, Cuesta JA, Aguirre J, Ahnert SE, Altenberg L, Cano AV, et al.  From genotypes to organisms: State-of-the-art and perspectives of a cornerstone in evolutionary dynamics.  arXiv preprint arXiv:200200363. 2020;.
  * 115. Novick A, Weiner M.  Enzyme induction as an all-or-none phenomenon.  Proceedings of the National Academy of Sciences (USA). 1957;43(7):553.
  * 116. Casadesús J, Low D.  Epigenetic gene regulation in the bacterial world.  Microbiology and molecular biology reviews. 2006;70(3):830–856.
  * 117. Casadesús J, Low DA.  Programmed heterogeneity: epigenetic mechanisms in bacteria.  Journal of Biological Chemistry. 2013;288(20):13929–13935.
  * 118. Ronin I, Katsowich N, Rosenshine I, Balaban NQ.  A long-term epigenetic memory switch controls bacterial virulence bimodality.  Elife. 2017;6:e19599.
  * 119. Mutlu A, Trauth S, Ziesack M, Nagler K, Bergeest JP, Rohr K, et al.  Phenotypic memory in Bacillus subtilis links dormancy entry and exit by a spore quantity-quality tradeoff.  Nature Comm. 2018;9(1):1–12.
  * 120. Carja O, Plotkin JB.  The evolutionary advantage of heritable phenotypic heterogeneity.  Scientific reports. 2017;7(1):1–12.
  * 121. Şimşek E, Kim M.  Power-law tail in lag time distribution underlies bacterial persistence.  Proceedings of the National Academy of Sciences. 2019;116(36):17635–17640.
  * 122. Hidalgo J, Pigolotti S, Muñoz MA.  Stochasticity enhances the gaining of bet-hedging strategies in contact-process-like dynamics.  Phys Rev E. 2015;91:032114.
  * 123. Villa Martin P, Hidalgo J, Rubio de Casas R, Muñoz MA.  Eco-evolutionary model of rapid phenotypic diversification in species-rich communities.  PLoS computational biology. 2016;12(10):e1005139.
  * 124. Zhang J, Zhou T.  Computation of stationary distributions in stochastic models of cellular processes with molecular memory.  bioRxiv. 2019; p. 521575.
  * 125. Pearson CE.  Numerical methods in engineering & science.  CRC Press; 1986.

SUPPORTING INFORMATION:

S1 Text: Sections: S1: The Microscopic process, S1A: Master equation, S1B:
Gillespie algorithm. S2: From Microscopic to Macroscopic Process: S2A.
Marginalization, S2B: Mean Field Approximation. S3: Small Variation
Approximation. S4: Deviation between theory and simulations, S4A: Validity of
the small variation approximation and border effects, S4B: FInite-size
effects. S5: Spontaneous shifting to the dormant state. S6: Additional
Figures. S7: Movies.

S1 Fig: Validity of the small variation approximation in the additive case.
Evolution of the first three cumulants ($K_{1},K_{2},K_{3}$) in a asymptotic
cycle calculated by the numerical integration of both the general mean-field
equation (S95)(solid line) and the generalized Crow-Kimura equation
(S114)(points) for two different values of $\alpha_{A}$, in (A)
$\alpha_{A}=0.16h$ and in (B) $\alpha_{A}=0.035h$. The first variation value
is the best fit the experimental results (main text), while the second causes
negligible border effects. In both cases the small-variation approximation is
a good approximation. Parameter values: $T_{a}=3h$, the rest of the
parameters, as well as the initial conditions, are kept fixed as specified in
the main text).

S2 Fig: Range of validity of the small variation approximation in the additive
scenario. Systematic comparison of Eq.(S114) and Eq.(S95) via parameter
Eq.(S124) for different $\alpha_{A}$ values, both axis in $log-$scale. Note
that the deviation monotonically increases with $\alpha_{A}$. In particular
for $\alpha_{A}=0.16h$, the one used in the main text results, the deviation
$\delta_{st.}=(3.5\pm 0.2)\cdot 10^{-2}h$, Eq.(S124) is small enough to use
Eq.(S114). Parameter values: $T_{a}=3h$, the other as in the main text.

S3 Fig: Validity of the small variation approximation in the multiplicative
scenario. Evolution of the first three cumulants ($K_{1},K_{2},K_{3}$) in a
asymptotic cycle calculated by the numerical integration of both the general
mean-field Eq.(S95)(solid line) and the generalized Crow-Kimura
Eq.(S114)(points). The generalized Crow-Kimura equation is a good
approximation of the general mean-field equation for the $\alpha_{M}$-value
that best fits the experimental results, e.g.$\delta_{st.}=(1.5\pm 0.3)\cdot
10^{-3}h$. Parameter values: $T_{a}=3h,\alpha_{M}=0.048$, the remaining are
specified in the main text.

S4 Fig: Border effects in the variation functions. First (A-B) and second
(C-D) moment of the variation functions as function of the trait $\tau$ in
additive (left) and multiplicative (right) variation cases. In (A) the
additive case, $\theta_{A}(\tau)$ is positive in $\tau=0$ but decays rapidly
to zero as $\tau$ increases, such that it is sufficient to restrict $\tau$
axis between 0 and 1. In (C) the same is shown for the second moment
$\sigma_{A}^{2}(\tau)$. Note that the magnitude of the dependence decreases
with $\alpha_{A}$). Nevertheless, the value used in the main text,
$\alpha_{A}=0.16h$, is too large to neglect Eq.(S120) this effect in the
generalized Crow-Kimura eq. On the other hand, in the multiplicative case the
main text value, $\alpha_{M}=0.048$, is small enough to avoid the border
effects in the GCK eq. In (B) one can observe that $\theta_{A}(\tau)$ almost
vanishes, and in (C) that the exact moment $\sigma^{2}_{M}(\tau)$ coincides
quite well with the approximation $\alpha_{M}\tau^{2}$, Eq.(S123).

S5 Fig: Schematic definition of the Parameter $\delta$ (illustrated in the
multiplicative amplitude scenario). $\alpha_{M}=0.048,T_{a}=6h$

S6 Fig: A) Functional dependence of the rate to enter the dormant state in
fresh medium $s_{f}$ and the number of particles in growing state $N_{G}$. B)
and C) Lag-time probability distribution function, $P(\tau)$, at the end of
the tenth cycle (for the multiplicative case) for a constant $s_{f}$ and
$s_{f}=s_{k}\left(\text{tanh}[-cN_{G}/K]\right)$, respectively. The results in
the main text, for $s_{f}=0$, are represented by a dashed line. Observe that
in both cases the value of $K_{1},K_{2},K_{3}$ increase, but the qualitative
form of the distributions remains unchanged. Parameter values:
$h_{\tau}=0.01$, $\alpha_{M}=0.048$, $T_{a}=3h$, $c=3$, $T=23h$, the rest of
parameters are fixed as in the main text.

S7 Fig: Dynamics of the averaged population structure within one cycle.
Abundances $N_{G}$ and $N_{D}$ along the $10$th cycle for $s_{f}$ constant —A)
and C)— and $s_{f}=\tilde{s}\left(\text{tanh}[-cN_{G}/K]\right)$ —B) and D)—.
The two scenarios are very similar to each other. At the beginning of the
cycle $N_{D}$ is non-vanishing, since bacteria entered this state during the
fresh phase of the previous cycle. During the antibiotic exposure phase both,
$N_{D}$ and $N_{G}$, decrease to a minimum as the bacteria die. When the
antibiotic is removed and a fresh medium is added, $N_{G}$ grows towards the
system’s carrying capacity. Observe that the bacteria still enter the dormant
state, but the reproduction rate is much higher, in such a way that an overall
reduction of $N_{G}$ is only observed as the system approaches the carrying
capacity.

S8 Fig: Characterization of the asymptotic state in the additive version of
the model. (A) Relaxation of the mean lag-time to its asymptotic state (curves
obtained from the integration of Equation Eq.(S114) with $T_{a}=3h$ and
additive model. The different curves correspond to three different values of
the variation amplitude, $\alpha_{A}$; from the lowest to the highest:
$\alpha_{A}=0.035h$, $\alpha_{A}=0.1h$ and $\alpha_{A}=0.16h$. (B) Zoom of the
curve $\alpha_{A}=0.16h$ for one single cycle. In particular, the mean lag
time, $K_{1}$, is shown along a cycle in the asymptotic regime. At $T=0$ the
antibiotic is added and the system enters in the “killing phase”
($T\in[0,T_{a}]$). When the antibiotic is present, the system experiences a
selection pressure towards longer lag times, in consequence, $K_{1}$
increases. At $T=T_{a}$ the antibiotic is washed and the fresh medium is added
($T\in[T_{a},T_{max}]$). In this regime, the selection pressure is towards
shorter lag times and $K_{1}$ relaxes back to the initial value. (C) Lag-time
probability distribution at $T=0$ (leftmost curve) and $T=T_{a}=3h$ (rightmost
curve) as derived theoretically (Eq.(S114), dashed lines) and computationally
(dots). In the asymptotic state the system oscillates between these two
limiting probability distributions, both of them exhibiting weak tails. (D)
Evolution of the first three cumulants, $K_{1},K_{2}$ and $K_{3}$, (mean,
variance , and skewness respectively) along a cycle in the asymptotic state
(both theoretical and computational results are displayed). Observe that in
C/D the theory correctly predicts the properties of the distribution, however
there are deviations due to finite size effects.

S9 Fig: Simulated number of cells via the Gillespie algorithm —additive-
amplitude scenario, linear scale—. a) Number of cell in dormant state during a
whole cycle in the asymptotic state for different exposition times $T_{a}$. b)
Same as a) but for growing state bacteria. Obviously, the minimum number of
bacteria in growing state $N_{G,min}$ is reached at $t=T_{a}$, since, during
the antibiotic exposure phase, that number can only decrease. Parameters:
$\alpha_{A}=0.16h$, the rest are fixed as in the main text.

S10 Fig: Evolution of the $MDK_{99}$ over 10 exposure cycles. a) Additive
case. b) Multiplicative case. In our simulations, the $MDK_{99}$ is calculated
analogously to the experimental procedure. After a certain number cycles (i.e.
$\\#$ in the figure) of antibiotics-fresh environment, the evolved population
is posed back in the antibiotic phase for a long time. The maximum number of
cycles is $10$ as in the experiments. The $MDK_{99}$ is estimated by the time
necessary to kill to the $99$% of the population. Both in the additive and
multiplicative case, the $MDK_{99}$ increases with the $\\#$ and $T_{a}$.
Interestingly, in the multiplicative case by increasing $T_{a}$ the change in
$MDK_{99}$ is bigger than the additive one.

S11 Fig: Comparison of experimental and simulated $MDK_{99}$ of the evolved
population after $10$ cycles of exposure. For both the additive and
multiplicative cases, we observe that the simulated $MDK_{99}$ falls within
—or it is very close— the experimental values (i.e. the mean plus error) for
$T_{a}=3h$ and $5h$, but outside for the case of $T_{a}=8h$. This result is to
be expected given the higher noise of the experimental measurements. In
particular the experimental mean is
$\langle\tau\rangle^{exp.}_{T_{a}=8h}=10\pm 1h$ higher than the theoretical
prediction of $8h$.

S1 Video: This supporting information file contains a video showing the
evolution of the asymptotic lag-time probability distribution for the additive
version of the model.

S2 Video: This supporting information file contains a video showing the
evolution of the asymptotic lag-time probability distribution for the
multiplicative version of the model.