# Memory shapes microbial populations

Chaitanya S Gokhale, Stefano Giaimo, Philippe Remigi

► To cite this version:

Chaitanya S Gokhale, Stefano Giaimo, Philippe Remigi. Memory shapes microbial populations. PLoS Computational Biology, 2021, 17 (10), 16 p. 10.1371/journal.pcbi.1009431 . hal-03385871

HAL Id: hal-03385871

https://hal.science/hal-03385871

Submitted on 19 Oct 2021

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# Memory shapes microbial populations

Chaitanya S. Gokhale<sup>1*</sup>, Stefano Giaimo<sup>2</sup>, Philippe Remigi<sup>3</sup>

1 Research Group for Theoretical Models of Eco-evolutionary Dynamics, Department of Evolutionary Theory, Max-Planck Institute for Evolutionary Biology, Plön, Germany, 2 Department of Evolutionary Theory, Max-Planck Institute for Evolutionary Biology, Plön, Germany, 3 LIPME, Université de Toulouse, INRAE, CNRS, Castanet-Tolosan, France

* gokhale@evolbio.mpg.de

<figure><img src="image_2.png" /></figure>

## Abstract

Correct decision making is fundamental for all living organisms to thrive under environmental changes. The patterns of environmental variation and the quality of available information define the most favourable strategy among multiple options, from randomly adopting a phenotypic state to sensing and reacting to environmental cues. Cellular memory—the ability to track and condition the time to switch to a different phenotypic state—can help withstand environmental fluctuations. How does memory manifest itself in unicellular organisms? We describe the population-wide consequences of phenotypic memory in microbes through a combination of deterministic modelling and stochastic simulations. Moving beyond binary switching models, our work highlights the need to consider a broader range of switching behaviours when describing microbial adaptive strategies. We show that memory in individual cells generates patterns at the population level coherent with overshoots and non-exponential lag times distributions experimentally observed in phenotypically heterogeneous populations. We emphasise the implications of our work in understanding antibiotic tolerance and, in general, bacterial survival under fluctuating environments.

<figure><img src="image_3.png" /></figure>

OPEN ACCESS

**Citation:** Gokhale CS, Giaimo S, Remigi P (2021) Memory shapes microbial populations. PLoS Comput Biol 17(10): e1009431. https://doi.org/10.1371/journal.pcbi.1009431

**Editor:** James O'Dwyer, University of Illinois at Urbana-Champaign, UNITED STATES

**Received:** March 11, 2021

**Accepted:** September 8, 2021

**Published:** October 1, 2021

**Peer Review History:** PLOS recognizes the benefits of transparency in the peer review process; therefore, we enable the publication of all of the content of peer review and author responses alongside final, published articles. The editorial history of this article is available here: https://doi.org/10.1371/journal.pcbi.1009431

**Copyright:** © 2021 Gokhale et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** https://github.com/tecoevo/phenotypic_heterogeneity.

**Funding:** This study received funding from the Max-Planck-Gesellschaft received by C.S.G. and S. G. Funding from the French Laboratory of Excellence (Agence Nationale de la Recherche

## Author summary

While being genetically the same, a population of cells can show phenotypic variability even under homogeneous environments. Often advantageous under heterogeneous environments, this *phenotypic heterogeneity* is highly relevant in the studies of antibiotic resistance evolution and cancer resurgence. Numerous theoretical models exist applying a simple model of phenotypic switching. Experimental measurements on phenotypic heterogeneity have increased in precision over the past decade, and the simple models are inadequate to explain the new observations. In this paper, we explore the role of cellular memory as a crucial component of phenotypic switching. We see that memory helps account for the hitherto unexplained observations and fundamentally extend our understanding of phenotypic heterogeneity.
project TULIP (ANR-10-LABX-41) was received by P.R. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared that no competing interests exist.

Introduction

In an ideal world, living organisms would be able to adapt instantly and reliably to changing environmental conditions to maximise their instantaneous performance. However, conditions may change abruptly and unpredictably, making it ineffective to mount a specific rapid response. Also, some responses require the synthesis of complex molecules (secretion systems or capsules in bacteria) or entering a physiological state (dormancy) that cannot be reverted instantaneously, if need be. Switching to a new phenotype may thus commit the cell to a response that lasts for a specific timescale different from the duration of the changing environment. Besides phenotypic switching, diversity in response rates can result in intricate patterns of phenotypic heterogeneity. We postulate that instantiating memory, i.e. some form of mechanistic regulation of the time spent in a state before switching to another, will affect population composition, size and ultimately fitness. While classical models of phenotypic heterogeneity use simple on-off switches, recent results stemming from accurate, modern experimental methods require us to delve deeper into the dynamics. Here we provide a first step towards a more realistic theoretical framework for phenotypic switching.

When environmental fluctuations show stereotypical patterns, unicellular organisms may harness this temporal information to adjust their mode of phenotypic adaptation to match, or even anticipate, environmental fluctuations [1]. Such strategies can be embedded in genetic regulatory networks [2, 3] or arise from epigenetic phenotypic switches [4, 5]. These forms of fitness optimization by associative learning can arguably be assimilated to memory-based processes [6, 7]. Far from being a phenomenon involving cognition or learning, memory in bacteria may emerge as a component of phenotypic heterogeneity [8–12]. In homogeneous environments, non-genetic individuality can arise through fluctuations in the concentration of signaling molecules during transcriptional bursts [13], unequal partitioning at cell division [14] or via other epigenetic mechanisms [15–19]. When these events operate on molecules involved in ultrasensitive responses, they open the possibility that a population of genetically identical cells cultivated in homogeneous environment might split into two (or more) sub-populations of cells exhibiting qualitatively different phenotypes, often characterized by the on-off expression switch of some genes [20]. Numerous studies have identified bimodal distributions of gene expression levels in a variety of organisms, and it was shown that such ‘in-built’ mechanisms for phenotypic variation help bacterial populations adapt to harsh environments [19, 21–23].

Nevertheless, in most cases, the conditions and the dynamics at which cells alternate between different phenotypic states are still obscure. If switching from one phenotype to the other only relies on one stochastic event (that typically follows a Poisson process), one can predict that residence time in each phenotypic state should be exponentially distributed [24]. This hypothesis underpins a large body of theoretical work dealing with the evolution of phenotypic heterogeneity in bacteria [22, 25–27]. But this constant phenotypic switch, however attractive for its modeling simplicity, is not a generic model when approaching the problem of intergenerational memory. For example, it is inadequate when explaining phenomena such as broad time-lag distributions and both over- and undershooting behaviour of specific phenotypes observed in bacteria or eukaryotic organisms [4, 28, 29]. Recent technological developments now allow measuring actual switching rates and phenotypic residence times of individual cells grown in a stable environment. By studying the motility/chaining phenotypic switch in *Bacillus subtilis*, Norman et al. [30] showed that residence time in the motile state was exponentially distributed, whereas it followed a gamma distribution in the chaining state. Similar observations of non-exponentially distributed residence times are observed for cells exiting lag phases [4, 29]. Explaining these behaviours at the molecular level requires considering the entire
molecular network controlling the switch and its cellular context. For example, excitable gene regulatory networks can display activity pulses of stereotyped duration [31]. In the case of the *B. subtilis* motility/chaining switch, the high-affinity binding of two proteins forming an inactive complex, coupled to their dilution at each cell division, explain the timed pulses of chaining-associated gene expression [32]. In this context, memory can thus be defined as the ability of unicellular organisms to track the time already spent in a given phenotypic state. This sensitivity allows cells to adjust their switching rates according to some ‘internal clock’ so that the time spent in this state follows a gamma distribution [24]. In contrast, cases where the probability of switching remains constant with time, are deemed memory-less.

Working within this conceptual framework [24, 30], we propose a mechanistic approach to model memory that does not abdicate the simplicity of earlier work. It would be possible to directly implement memory by preselecting an appropriate non-exponential waiting time distribution. However, that does not help decipher a mechanism to how the distribution materials. Hence, first, we derive deterministic dynamics from first-principles at the level of individual cells and track the behaviour of a cell population emerging from different switching genotypes. Using this model, we can consistently explain the experimental observation of over/undershoots and point to a possible way to understand wide time (non-exponential) lag distributions. Then, paying particular attention to the characterisation of transient states of the system, we discuss their potential consequences on the fitness of a lineage in the presence of fluctuating environments. Significantly, fluctuations across stressed and relaxed environments, as in specific antibiotic treatment regimes, highlight the applicability of our approach.

# Results

## Building memories

We simulate the ecological dynamics of isogenic populations of unicellular organisms that can exist in two different phenotypic states: ‘on’ and ‘off’. Switching from ‘off’ to ‘on’ is unidirectional and stochastic and occurs at rate $µ$. After turning ‘on’, cells cascade via a deterministic, multi-step process through *n* compartments, eventually returning to the ‘off’ phenotype. These compartments represent internal molecular substates (potential), that may correspond to various biological processes, such as the dilution of cytoplasmic or membrane proteins [8, 11, 13, 14], or the sequential realization of independent molecular reactions collectively required to trigger a response [33–35].

Immediately after turning ‘on’, cells have the highest potential. While retaining the ‘on’ state, cells gradually lose potential by transitioning through the successive compartments at a leaching rate of $ε$. In our model, this movement reflects a decrease in protein concentration in a simplistic manner, whereas more complicated forms can be formulated (with bumps or plateaus on the landscape, see S1 Text). Cell flow through compartments ($i = n, ..., 1$) is a representation of phenotypic dynamics, while growth (cell birth and death) dynamics occur separately at rates $b_i$ and $d_i$, respectively. Admittedly, birth and death can themselves lead to cell movement from one compartment to another (via symmetric or asymmetric division of the protein concentration in the daughters) [36]. We have explored such a model where cells differentiate into the downstream compartments separately and is available online on GitHub. However, to sort out the effects of leaching from growth and keep the model simple, in this manuscript, we have not included such scenarios.
Fig 1 visualizes this compartmental model for $n = 4$. The reactions in which an individual cell $X_i$ in compartment $i$ is involved are then,

$$ 
\begin{aligned}
X_i & \xrightarrow{b_i} X_i + X_i \\
X_i & \xrightarrow{d_i} \emptyset \\
X_i & \xrightarrow{\epsilon} X_{i-1}, \quad i = n, \dots, 1 \\
X_0 & \xrightarrow{\mu} X_n,
\end{aligned} 
\quad (1) $$

These reactions have a deterministic counterpart in the following linear system of differential equations,

$$ 
\begin{aligned}
\dot{x}_0 &= (b_0 - d_0)x_0 - \mu x_0 + \epsilon x_1 \\
\dot{x}_n &= (b_n - d_n)x_n - \epsilon x_n + \mu x_0 \\
\dot{x}_i &= (b_i - d_i)x_i - \epsilon x_i + \epsilon x_{i+1},
\end{aligned} 
\quad (2) $$

which describe the time evolution of the abundances $x_i$ of cells in the compartments with $i = 1, \dots, n-1$. We envisage that once a cell is in an ‘off’ state, it stays ‘off’ until an external or internal perturbation of sufficiently large magnitude switches it back ‘on’. The origin of perturbations may be abiotic or biotic [21, 37]. For $n = 1$, we recover the classically studied memoryless switch. In what follows, we consistently compare this more classic model with two states and two compartments, which shows no memory according to our definition, with the case of our interest $n > 1$, i.e. two states yet multiple ‘on’ compartments, which can be shown to generate a form of memory that correlates with the magnitude of $n$. In the subsequent sections, we compare the performances of the two models ($n = 1$ versus $n > 1$) with particular attention to the transient phase. We ensure that the models can be considered equivalent from the point of view of long term properties by tuning the $\epsilon$ parameter accordingly.

However, we end this section by showing how our model meets the minimal requirement of producing memory. Tracking the amount of time a cell spends in a particular compartment makes this evident. Following the trajectories of individual cells, qualitatively different distributions of the time spent in the ‘on’ state emerge as the number of compartments $n$ increases. In particular, departure from an exponential distribution, the hallmark of memory [30] (Fig 1B) in this context, is observed. For a constant leaching rate and negligible growth parameters, the number of compartments acts like a timer (a deterministic time as termed in [30] for the residence time in the motile state). The overall density function can be captured by a combination of multiple exponential waiting times which results in a gamma distribution with a shape parameter given by the memory size $n$ (Fig 1B). This shows that our model satisfies the minimal requirement of generating memory as a deviation from an exponentially distributed waiting time before switching to ‘on’. As the magnitude of this deviation relates directly to $n$, we conveniently refer to this parameter as the length or size of memory. Codes for implementing our algorithm as well as for reproducing the relevant figures are available at GitHub.

## Asymptotic properties

The asymptotic properties of the model were already well described and do not represent the primary focus of our study [30]. However, to make further progress, we recapitulate them briefly. As the more classic two-state two-compartment model, our model is linear and Markovian: the current system state entirely determines cell dynamics. The equilibrium cell distribution in the compartments can be recovered from the dominant eigenvector of the constant matrix that captures the system in Eq 2 (S1 Text). The corresponding eigenvalue corresponds
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Fig 1.** Growth and phenotypic dynamics of ‘off’ and ‘on’ cells. (A) The growth dynamics of cells is given by a birth process and a death process. A cell divides into two cells at rate $b_i$ and dies at rate $d_i$. The index $i$ refers to the phenotypic state of being ‘off’ or ‘on’. The ‘off’ cells are the cells in the resting state which is the zeroth compartment $i = 0$. Due to a trigger (an internal constant or a dependence on frequency/density/environmental state) the cells can be turned ‘on’ at rate $\mu$ and they jump to the $n^{th}$ compartment (here $n = 4$). The cells do not stay in the ‘on’ state but slowly decay back to the ‘off’/resting state via passaging at rate $\epsilon$ through a number of intermediate phenotypic stages, the different substates of the ‘on’ state. The number $n$ of compartments in the ‘on’ state modulates memory. With $n = 1$, there is zero memory, while for $n > 1$ memory increases in the number of compartments, i.e. the increase is here understood as the increase in</figcaption></figure>
the average time spent in the ‘on’ state relative to the average of the exponentially distributed time spent in this state when $n = 1$. Thus we intend to make use of this mechanistic interpretation of memory instead of an assumption of a separation between the amount of time a cell spends in ‘on’ or ‘off’ state. (B) Distribution of times a cell stays ‘on’ before switching ‘off’. Starting with a single cell in the completely ‘on’ state we ask how long it takes for the cell to reach the ‘off’ state. Assuming extremely low replication and death rates $b_i = d_i = 10^{-6}$, the process is governed by the leaching rates $\epsilon = 0.3$. For 1, 2 and 4 ‘on’ compartments we start the Gillespie simulation with one cell in the first ‘on’ state. Once the cells reaches the ‘off’ state we stop the simulation. The normalised histogram of waiting times is then plotted. The probability density function of a gamma distribution with the shape parameter given by the number of ‘on’ compartments and the rate $\epsilon = 0.3$ provides the theoretical estimates (black lines) for the simulated 1, 2 and 4 ‘on’ compartments. Increasing the memory, ($n = 6, 8$ and 10) flattens the distribution of the time spent in the ‘on’ state.

https://doi.org/10.1371/journal.pcbi.1009431.g001

to asymptotic population growth. Under complete symmetry in the growth dynamics ($b_i = b$ and $d_i = d$) we can get a simple, closed-form expression for this equilibrium for any number $n$ of compartments. We begin with a focus on this symmetric case, where the asymptotic growth rate $b-d$ of the population is independent of the number of compartments. The assumption is relaxed further in the manuscript. Our interest is in the effect of changes in the leaching rate $\epsilon$ and memory $n$. We observe that as $\epsilon$ decreases or $n$ increases, both the time a cell spends in the ‘on’ state and the equilibrium fraction of ‘on’ cells increases (see S1 Text). For a small number of compartments or faster transition through ‘on’ states, the system is dominated by ‘off’ cells.

Transient properties

Bacteria in natural environments are frequently exposed to changing conditions. These fluctuations can often keep populations from reaching evolutionary or ecological steady states, and transient dynamics become crucial. Although less amenable to analytical treatment, we shall focus on them requiring greater reliance on computational exploration. While under constant switching ($\mu$), a decrease in $\epsilon$ is qualitatively equivalent to an increase in $n$ in terms of asymptotic properties (i.e. the equilibrium cell distribution and the time spent in either state), this equivalence breaks down when considering transient dynamics. As the number of compartments increases, the frequencies of the ‘on’ and ‘off’ cells do not approach their equilibrium values directly. Instead, in the case of ‘on’ cells, the frequencies overshoot and undershoot, in the case of ‘off’ cells. The frequency dynamics oscillates when approaching equilibrium (Fig 2A).

Interestingly, overshooting and oscillations in the transients of cell frequencies as seen in our model recapitulate observations in cancer cells dynamics [28] and in bacterial growth rate recovery following antibiotic treatments [38]. This transient effect directly relates to the length of memory. While modelling the ‘off’ and ‘on’ states as two compartments would appear more parsimonious, it falls short of replicating empirical results. Adding compartments intensifies transients, an effect attributable to the spectrum of the matrix model underlying Eq 2 (see S1 Text). The presence of multiple compartments introduces oscillations due to complex subdominant eigenvalues, absent when there is a single ‘on’ compartment. Increasing the number of compartments also magnifies the influence of complex subdominant eigenvalues on cell dynamics, as the real parts of such eigenvalues get closer to the real part to the dominant eigenvalue. Importantly, this analysis is possible while preserving linear dynamics and, thus, a form of simplicity. This, we suggest, would be lost if non-exponentially distributed waiting times were coded directly into a two-state two-compartment model via the introduction of non-linear terms in a coupled system of two differential equations.

Cells switch back to the ‘off’ state after spending a certain amount of time in the ‘on’ state. The presence of multiple compartments extends memory and affects population composition
<figure><img src="image_6.png" /><figcaption>**Fig 2.** *Transient dynamics of multi-state memory. (A).* Multi-state memory and overshoots. Assuming full symmetry in growth dynamics between the compartments ($b_i = b = 1.0$ and $d_i = d = 0.98$) and a switch rate of $\mu = 0.2$ and $\epsilon = 0.5$ we show the population dynamics as well as the population composition for fixed time period. The instantaneous growth rate of the two population with different memory sizes is the same $g = 1/t_{max} \log(N_{final}/N_{initial})$. However, with a larger memory, the fraction of 'on' cells in the final population is higher. On the way to the equilibrium, multiple compartments result in overshoot dynamics as the 'on' cells first need to seed them while depleting the 'off' state. (B). *Transient magnitude and number of compartments.* To characterise the oscillations, we plot the peak-to-peak amplitude (the amplitude between the first peak, i.e. overshoot, and the second peak, i.e. undershoot) against the equilibrium value of 'on' cells as the number $n + 1$ of compartments increases ($n$ 'on' compartments and one 'off' compartment, gray numbers). The equilibrium fraction of cells in the 'on' state is given by $n\mu/(\epsilon + n\mu)$ (S1 Text). Keeping $\mu$ fixed, the peak-to-peak amplitude is obtained by varying the leaching rate $\epsilon$ and by varying the number $n$ of compartments in the multi-compartment system, which has constant $\epsilon$. As the equilibrium fraction of 'on' cells approaches unity, boundary effects prevail, i.e. transient frequency cannot exceed 1, and oscillations' amplitude gets compressed. The right panel shows the temporal dynamics of 'on' cells, all leading to the same equilibrium of 0.8 but for different leaching rate $\epsilon$ (vertical line in the left panel). The corresponding memory sizes are $n = 5, 15, 25, 35$ respectively. Increasing the number of compartments leads to more pronounced oscillations, not observed in the system with a single 'on' compartment regardless of the value of $\epsilon$ and the equilibrium value of 'on' cells. Codes for generating these panels and the general algorithm are available on GitHub.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1009431.g002

as well as population density when evaluated at a specific time-point (Fig 2A). To capture the transient dynamics when increasing memory, we estimate the peak-to-peak amplitude (the difference between the largest overshoot and undershoot in the frequency of 'on' cells) (Fig 2B). Typically we see a monotonic increase in the relationship between the amplitude and memory length. The relationship wanes as the equilibrium value of 'on' cells approaches 1 due to boundary effects Fig 2B). More generally, the magnitude of the over/undershoot depends on the initial conditions and equilibrium of the system.

## Environmental variation

While cell lineages can stochastically switch phenotypes even in static conditions [39], such behaviour might be evolutionarily advantageous under fluctuating environments [40]. As the
environment changes, bacteria can hedge their bets via the well-documented phenomena of persistence [41, 42], a salient case being bacterial response to antibiotic treatment strategies.

To explore the potential effects of memory on bacterial fitness, we consider the case of per-
sisters when subjected to transient antibiotic exposure. Persistence is a common phenomenon
in bacteria where a subpopulation of bacterial cells does not grow but can survive antibiotic
treatments. Persisters can arise stochastically via epigenetic switching, but environmental con-
ditions (nutrient limitation, high cell densities or antibiotic treatments) can also induce their
formation (‘triggered’ persistence or type I persisters; [37]). We track population composition
and fitness (as proxied by population size) of bacterial lineages under fluctuating environ-
ments, consisting of transient exposure to antibiotics. We assign ‘off’ cells to the ‘normal’ phys-
iological state (i.e. growing in permissive conditions, dying under antibiotic treatment). On
the other hand, ‘on’ cells are persisters (no growth in either environment and tolerant to anti-
biotic treatment; see description of parameters in Fig 3 and example dynamics in S1 Text).
While the leaching rate is kept constant, switching from ‘off’ to ‘on’ is triggered only under
stress. We consider an initial population of ‘off’ cells, representing a stable equilibrium state
under conducive growth conditions.

We subject a growing population of cells to a set of 84 environmental sequences where bac-
teria are exposed to drug treatment sequences (each horizontal line in Fig 3A corresponds to a
different environmental sequence). The total duration of the procedure is kept constant, and
the duration of drug treatment varies, with exposure to the permissive environment before
and after drug exposure (also see S1 Text). Using these environmental sequences, we compute
the fitness of a set of lineages with different memory sizes for each condition. We observe a
non-linear relationship between memory size and fitness (Fig 3B). While there is a general
trend that more memory is beneficial when drug treatment increases, local maxima emerge at
intermediate drug treatment lengths (Fig 4). This effect arises from the trade-off between

<figure><img src="image_7.png" /></figure>

**Fig 3.** Growth performance of lineages with different memory lengths exposed to treatment of varying length. Cell lineages with different memory sizes (n = (1, ..., 31) + 1) are temporarily exposed to antibiotic treatment. Each row in the grid is a sequence of seasons from left to right. Seasons with no-treatment (white) and treatment (black) where each lasts for one unit of time. All sequences begin at t = 0 with 1000 'off' cells and last until tmax = 165. Under the no-treatment season the growth rate of 'off' cells is boff - doff = 1 - 0.98 = 0.02 and the switch is inactive μno treatment = 0. As cells encounter the treatment season (black squares), the switch is triggered, μno treatment = 0.2 and the death rate of 'off' cells increases, boff - doff = 1 - 1.02 = -0.02. The 'on' cells are produced but they do not grow, bon - don = 0. The 'on' cells leach through the memory compartments at rate ϵ = 0.25. At the end of the sequence (t = 165), the absolute growth rate, gm = log (Nmemory(tmax)/Nmemory(0)), was computed. A memory-less process (two compartments with n = 1 and an appropriate ϵ) can generate the same stable fraction of 'on' cells under sustained treatment. We subtract the absolute growth rate of such a memory-less lineage gmless = log (Nmless(tmax)/Nmless(0)) from gm to estimate the relative fitness of having memory r = (gm - gmless)/tmax.

https://doi.org/10.1371/journal.pcbi.1009431.g003
<figure><img src="image_8.png" /><figcaption>**Fig 4.** Lineage adaptation under treatment. Top panels show the effect of memory size on absolute and relative lineage fitness. Each of the lines in the top panels corresponds to each of the 84 sequences from Fig 3. The colorscheme of the left and right top panels is the same as from Fig 3B and 3C with the 42^nd^ sequence highlighted in bold. We report the frequency of 'on' cells and the population size in the bottom panels. The results are shown for a chosen number of memory sizes (n) + 1 ('off' compartment). The non-linear relationship between the final 'on' cell frequency and the final population size is shown in the last panel as n + 1 ranges from 2 to 32.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1009431.g004

producing a high proportion of 'on' cells (occurring with higher memories) and exiting rapidly from the 'on' state when the environment switches from drug to permissive (occurring with smaller memories). For extended treatment lengths, longer memory in such cases staves off an inevitable population collapse.

Memory is coupled to the transients and the final frequency of 'on' cells observed at the end of each sequence. To disentangle the respective effect of these two factors, we compared the fitness of a lineage with memory (n > 1) to a lineage without memory (n = 1). To make lineages comparable, the memory lineage has a different leaching rate which results in the same equilibrium 'on' frequency under sustained treatment Fig 3C. The relative growth rate is the difference between the growth rate of a lineage with memory and without memory following http://myxo.css.msu.edu/ecoli/srvsrf.html. When the treatment length is short, longer memory is disadvantageous as compared to larger leaching rates. The disadvantage arises because the cells stay in the 'on' state even after the short treatment has elapsed. For intermediate treatment
lengths, more memory allows (i) to produce more ‘on’ cells and (ii) longer residual time in the ‘on’ state, then becoming advantageous. However, long residual times corresponding to longer lag times could conflict with the total sequence length where the fitness is measured. Hence, very long memory is also not helpful as the cells then take much longer to exit the ‘on’ state than the sequence length. Thus, we see the presence of local maxima in memories (driven partly by the absolute fitness of the memoryless lineage, Fig 4). Under lengthy sustained treatment, all lineages are at the ‘on’ equilibrium, and memory length is inconsequential. This analysis reveals that memory outperforms ‘classical’ memoryless switching against an equal equilibrium value of ‘on’ cells for a range of conditions.

Discussion and conclusion

Studies on phenotypic heterogeneity and persistence in bacteria have typically relied on models assuming two distinct cellular states, ‘off’ and ‘on’ with constant switching rates. We have argued that this assumption hinders us from understanding the rich dynamics observed in empirical studies. Keeping the two states model while decomposing the ‘on’ state into multiple, interconnected ‘on’ substates performs better in reconciling observations. The key to this reconciliation is that having more than one ‘on’ state generates memory. In this context, memory is the ability to revert to ‘off’ while accounting for the time spent in ‘on’ [24].

Multiple ‘on’ states are also implemented in a recent theoretical study studying the optimal switch time dependence on the likelihood and duration of a hazardous environment [43]. This study’s overall approach was to look at evolutionary stable states but ignores transient population dynamics and stochasticity, an integral component of experimental findings. Our contribution fills this gap. We show that memory, originating from a process going through multiple substates, appears to explain otherwise anomalous observations.

Take, for example, the time it takes to resume division for a single cell taken out of a tolerant population that has just exited a prolonged antibiotic treatment [37]. This duration is the so-called lag time. Under artificial selection experiments, cells adapt their lag time to match treatment duration [4, 29] while displaying within-population variability. This variability increases with the mean lag time (see S1 Text). Initially, this observation seems counter-intuitive, as the best strategy would be for all cells to resume growth as soon as the treatment is over. Moreno-Gámez et al. [29] showed that heterogeneous lag times could promote survival to transient antibiotic treatments while having a negligible cost on population growth. In these studies, the lag time was shown to have a gamma distribution, which is not compatible with a memory-less process of switching from ‘on’ to ‘off’. Our model finds that a cell in any ‘on’ state will exhibit a gamma-distributed lag time, coherent with observations.

It is possible to study how the lag time evolves in response to treatment likelihood and length [43]. Our framework, although not tackling this evolutionary problem directly, is suggestive of the fact that the variability in lag times upon which selection can act may potentially arise through multiple sources: changes in the number of ‘on’ states/compartments, in the leaching rate, in both or competition between lineages with variable properties of the gamma distribution. However, for a fuller understanding of the evolutionary forces acting on these traits, explicit considerations of their potential trade-offs will be required.

In experiments, the details of the growth dynamics often play an essential role. For example, cells may be first let grow exponentially, then starved into a stationary phase and finally sampled to regrow on abundant, fresh media. In our model, growth dynamics were minimized on purpose to focus on phenotypic dynamics. Moreover, while sampling from a single ‘on’ compartment would produce gamma-distributed lag times, sampling from multiple ‘on’ compartments would result in a lag time that is not gamma distributed, although emerging from a
combination of different gamma distributions (see S1 Text). Since our model abstracts bacte-
rial populations’ natural growth dynamics, we do not produce a qualitative match between the
experimental data on lag times. To approximate experiments, it may be possible to design a
two-compartment model with different waiting time distributions that recapitulate some of
the properties that we present in this study. However, this would beg to reason the mecha-
nisms resulting in the different distributions. Instead, we aim to postulate a mechanistic pro-
cess that encompasses the internal “to-do” list before switching phenotype as suggested in
[30]. Also, the compartments allow us to mathematically connect the number of ‘on’ states to
the observed transient dynamics quantitatively.

Expressly, the transients of the dynamics can be understood in terms of eigenvalue analysis.
Complex eigenvalues introduce oscillations in the dynamics (S1 Text). Larger memory sizes
correspond to more solutions in the imaginary space, which are reflected in the dynamics with
more oscillations (Fig 2). The magnitude of the effect hinges on how fast cells experience the
memory (leaching rate) and the initial switch rate. For any leaching rate, however, as the mem-
ory increases, the oscillations (captured by the peak-to-peak amplitude) increase but only up
to a limit (Fig 2B). The difference is solely in the fraction of the ‘on’ cells and not the popula-
tion size. The result rests on our assumption of ‘on’ and ‘off’ cells having the same birth and
death rates. Forgoing this assumption would lead to further divergence between the composi-
tion and size of lineages with different memories. Whether fitness depends on the composition
or size of a lineage, memory will bring unique dynamical properties that might impact survival.
The ecological context will set the timing of fitness evaluation, realising Darwinian selection
on lineages (Fig 2).

Phenotypic heterogeneity while advantageous for a lineage [27, 44] can be a nuisance when population expansion is harmful. Since Hobby and Bigger [41, 42] persisters have been a fly in the ointment for antibiotic treatment only exacerbated now by the antibiotic crisis. We have presented the structure of the use of multi-state memory and its application to antibiotic tolerance. The persisters, as defined in our case, are a subpopulation of tolerant cells appropriately defined as “heterotolerant” [37]. While tolerance does not affect minimum inhibitory concentration, the duration of treatment will be crucial for eradicating bacteria. Another application of our approach could be in understanding the persistence and dormancy mechanisms in cancer populations when underlined by phenotypic switching [45]. Quiescent subclones are a persistent problem leading to cancer relapse [46–48]. The subclone population dynamics also show overshooting [28, 49]. A tunable and evolvable memory could then underlie the mechanisms of dormancy categorised as adaptive response [50]. Thus considering memory brings a different time-scale that can be exploited for steering evolution to control pathogenic populations.

Numerous extensions of our approach are possible. For example, we have focused on het-
erogeneity resulting from only environmentally triggered switches, a property observed in sev-
eral experimental models [10, 12, 37, 51]. Alternative switching mechanisms may depend on
the intrinsic properties of the population, such as density or composition. Similarly, the leach-
ing mechanism is predetermined and constant across the compartments. Such processes may
be under the complex, joint control of the organism as well as the environment. Although
detailed experimental characterisation of bidirectional switching behaviours remains rare (due
to technical challenges), we expect that memory-based switching is not an exception.

Oscillations involving excitation and decay can be identified across life. For example, excit-
able genetic switches are found in bacteria [52], belonging to circadian clocks in cyanobacteria
[53], in other non-photosynthetic bacteria [54], in Drosophila [55] and in plants [56, 57]. The
molecular mechanisms identified in these systems point towards a malleable duration of the
oscillations by mutations [52, 58], modifications of the regulatory network [32] or degradation
and destabilisation of protein complexes [55]. Similar to the external cycling antibiotic stimu-
lus explored in this study, the explicit dependence on external signals, flies' chronobiology has
been proposed as a tool to understand adaptation to variable environmental conditions. Thus
the adaptive evolution of molecular processes underpinning the oscillations of the excitable
systems such as our model is conceivable.

While theoretical developments are essential, the applied aspect can be further exploited.
Beyond antibiotics and cancer treatment, bioengineering and understanding microbial con-
sortia formation could be informed using our approach. Knowing the memory limitations
(e.g. coming from decay rates of protein complexes) involved in critical oscillatory processes
such as in chronobiology can inform us about the limits of adaptation under extreme environ-
mental events or anthropogenic changes. Especially when harsh environments and time-lags
are of importance, such as in niche construction and the evolution of multicellularity [30, 59].

Understanding gene regulatory networks on a developmental landscape, à la Waddington [45, 60], poses exponentially complex computational challenges (e.g. the explosion of possible attractors when considering multiple switches [61] and multiple phenotypes) (as also discussed in [62]). We show multiple local maxima for memory, which may change depending on the definition of fitness (population size or composition). Furthermore, this forces us to rethink our concepts about possible phenotypic states and how they are determined by a plethora of molecular constructs (see S1 Text for an extended discussion). How epigenetic memory functions over generations would focus on understanding how the phenotypic clock (molecules, appropriate histone modifications) are inherited. In our model, time, as kept by degrading molecules, is arbitrary. Inputs from chronobiology could help see how the molecular concentrations and circuits are entrained, resulting in biologically driven growth and leaching dynamics.

Phenotypic heterogeneity forces us to reassess the genotype-phenotype map fundamentally. Choosing the appropriate phenotypic response to complex and varied environments is possible via numerous processes such as environmental sensing, epigenetic triggers and controlled molecular concentrations. Such processes that interpret genetics to a large but finite phenotypic space to survive in a possibly infinite environmental space are incredibly relevant for natural selection. Theories, as described herein, coupled with experiments exploring diverse environments, will help us elucidate the variety of possible interpreting mechanisms bridging the genotype-phenotype divide.

**Supporting information**

S1 Text. Supplementary text for memory shapes microbial populations.
(PDF)

Acknowledgments

We thank Silvia De Monte and Orso Romano for helpful discussions.

**Author Contributions**

**Conceptualization:** Chaitanya S. Gokhale, Philippe Remigi.

**Formal analysis:** Chaitanya S. Gokhale, Stefano Giaimo, Philippe Remigi.

**Funding acquisition:** Chaitanya S. Gokhale, Philippe Remigi.

**Investigation:** Chaitanya S. Gokhale, Stefano Giaimo, Philippe Remigi.

**Methodology:** Chaitanya S. Gokhale, Stefano Giaimo, Philippe Remigi.
**Project administration:** Chaitanya S. Gokhale.

**Resources:** Chaitanya S. Gokhale, Philippe Remigi.

**Software:** Chaitanya S. Gokhale, Stefano Giaimo.

**Validation:** Chaitanya S. Gokhale, Stefano Giaimo.

**Visualization:** Chaitanya S. Gokhale, Stefano Giaimo, Philippe Remigi.

**Writing – original draft:** Chaitanya S. Gokhale, Stefano Giaimo, Philippe Remigi.

**Writing – review & editing:** Chaitanya S. Gokhale, Stefano Giaimo, Philippe Remigi.

## References

1. Bernhardt JR, O'Connor MI, Sunday JM, Gonzalez A. Life in fluctuating environments. Philosophical Transactions of the Royal Society B. 2020; 375(1814):20190454. https://doi.org/10.1098/rstb.2019.0454 PMID: 33131443

2. Tagkopoulos I, Liu YC, Tavazoie S. Predictive Behavior Within Microbial Genetic Networks. Science. 2008; 320(5881):1313-1317. https://doi.org/10.1126/science.1154456 PMID: 18467556

3. Mitchell A, Romano GH, Groisman B, Yona A, Dekel E, Kupiec M, et al. Adaptive prediction of environmental changes by microorganisms. Nature. 2009; 460(7252):220-224. https://doi.org/10.1038/nature08112 PMID: 19536156

4. Fridman O, Goldberg A, Ronin I, Shoresh N, Balaban NQ. Optimization of lag time underlies antibiotic tolerance in evolved bacterial populations. Nature. 2014; 513(7518):418-421. https://doi.org/10.1038/nature13469 PMID: 25043002

5. Bergh BVd, Michiels JE, Wenseleers T, Windels EM, Boer PV, Kestemont D, et al. Frequency of antibiotic application drives rapid evolutionary adaptation of Escherichia coli persistence. Nature Microbiology. 2016; 1(5):16020. https://doi.org/10.1038/nmicrobiol.2016.20 PMID: 27572640

6. Sorek M, Balaban NQ, Loewenstein Y. Stochasticity, Bistability and the Wisdom of Crowds: A Model for Associative Learning in Genetic Regulatory Networks. PLoS Computational Biology. 2013; 9(8): e1003179. https://doi.org/10.1371/journal.pcbi.1003179 PMID: 23990765

7. Lyon P. The cognitive cell: bacterial behavior reconsidered. Frontiers in Microbiology. 2015; 6:264. https://doi.org/10.3389/fmicb.2015.00264 PMID: 25926819

8. Wolf DM, Fontaine-Bodin L, Bischofs I, Price G, Keasling J, Arkin AP. Memory in Microbes: Quantifying History-Dependent Behavior in a Bacterium. PLoS ONE. 2008; 3(2):e1700. https://doi.org/10.1371/journal.pone.0001700 PMID: 18324309

9. Lambert G, Kussell E, Kussel E. Memory and Fitness Optimization of Bacteria under Fluctuating Environments. PLoS Genetics. 2014; 10(9):e1004556. https://doi.org/10.1371/journal.pgen.1004556 PMID: 25255314

10. Ronin I, Katsowich N, Rosenshine I, Balaban NQ. A long-term epigenetic memory switch controls bacterial virulence bimodality. eLife. 2017; 6:e19599. https://doi.org/10.7554/eLife.19599 PMID: 28178445

11. 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 Communications. 2018; 9(1):69. https://doi.org/10.1038/s41467-017-02477-1 PMID: 29302032

12. Kordes A, Grahl N, Koska M, Preusse M, Arce-Rodriguez A, Abraham WR, et al. Establishment of an induced memory response in Pseudomonas aeruginosa during infection of a eukaryotic host. The ISME Journal. 2019; 13(8):2018-2030. https://doi.org/10.1038/s41396-019-0412-1 PMID: 30952997

13. Sanchez A, Choubey S, Kondev J. Regulation of Noise in Gene Expression. Annual Review of Biophysics. 2013; 42(1):469-491. https://doi.org/10.1146/annurev-biophys-083012-130401 PMID: 23527780

14. Bergmiller T, Andersson AMC, Tomasek K, Balleza E, Kiviet DJ, Hauschild R, et al. Biased partitioning of the multidrug efflux pump AcrAB-TolC underlies long-lived phenotypic heterogeneity. Science. 2017; 356(6335):311-315. https://doi.org/10.1126/science.aaf4762 PMID: 28428424

15. Spudich JL, Koshland DE. Non-genetic individuality: chance in the single cell. Nature. 1976; 262 (5568):467-471. https://doi.org/10.1038/262467a0 PMID: 958399

16. Thattai M, Oudenaarden Av. Intrinsic noise in gene regulatory networks. Proceedings of the National Academy of Sciences. 2001; 98(15):8614-8619. https://doi.org/10.1073/pnas.151588598 PMID: 11438714
17. Maamar H, Raj A, Dubnau D. Noise in Gene Expression Determines Cell Fate in Bacillus subtilis. Science. 2007; 317(5837):526-529. https://doi.org/10.1126/science.1140818 PMID: 17569828

18. Casadesús J, Low DA. Programmed heterogeneity: epigenetic mechanisms in bacteria. Journal of Biological Chemistry. 2013; 288(20):13929-13935. https://doi.org/10.1074/jbc.R113.472274

19. Ackermann M. A functional perspective on phenotypic heterogeneity in microorganisms. Nature Reviews Microbiology. 2015; 13(8):497-508. https://doi.org/10.1038/nrmicro3491 PMID: 26145732

20. Ferrell JE, Ha SH. Ultrasensitivity part I: Michaelian responses and zero-order ultrasensitivity. Trends in Biochemical Sciences. 2014; 39(10):496-503. https://doi.org/10.1016/j.tibs.2014.08.003 PMID: 25240485

21. Kussell EE, Kishony RR, Balaban NQN, Leibler SS. Bacterial persistence: a model of survival in changing environments. Genetics. 2005; 169(4):1807-1814. https://doi.org/10.1534/genetics.104.035352 PMID: 15687275

22. Kussel E, Leibler S. Phenotypic Diversity, Population Growth, and Information in Fluctuating Environments. Science. 2005; 309:2075-2078. https://doi.org/10.1126/science.1114383

23. Boxtel Cv, Heerden JHv, Nordholt N, Schmidt P, Bruggeman FJ. Taking chances and making mistakes: non-genetic phenotypic heterogeneity and its consequences for surviving in dynamic environments. Journal of The Royal Society Interface. 2017; 14(132):20170141. https://doi.org/10.1098/rsif.2017.0141 PMID: 28701503

24. Norman TM, Lord ND, Paulsson J, Losick R. Stochastic Switching of Cell Fate in Microbes. Annual Review of Microbiology. 2015; 69:381-403. https://doi.org/10.1146/annurev-micro-091213-112852 PMID: 26332088

25. Lohmar I, Meerson B. Switching between phenotypes and population extinction. Physical Review E. 2011; 84:051901. https://doi.org/10.1103/PhysRevE.84.051901 PMID: 22181438

26. Carja O, Libermann U, Feldman MW. The Evolution of Phenotypic Switching in Subdivided Populations. Genetics. 2014; 196(4):1185-1197. https://doi.org/10.1534/genetics.114.161364 PMID: 24496012

27. Carja O, Plotkin JB. The evolutionary advantage of heritable phenotypic heterogeneity. Scientific Reports. 2017; 7(1):5090. https://doi.org/10.1038/s41598-017-05214-2 PMID: 28698577

28. Sellerio AL, Ciusani E, Ben-Moshe NB, Coco S, Piccinini A, Myers CR, et al. Overshoot during phenotypic switching of cancer cell populations. Scientific reports. 2015; 5(1):15464. https://doi.org/10.1038/srep15464 PMID: 26494317

29. Moreno-Gámez S, Kiviet DJ, Vulin C, Schlegel S, Schlegel K, Doorn GSv, et al. Wide lag time distributions break a trade-off between reproduction and survival in bacteria. Proceedings of the National Academy of Sciences. 2020; p. 202003331.

30. Norman TM, Lord ND, Paulsson J, Losick R. Memory and modularity in cell-fate decision making. Nature. 2013; 503(7):481-486. https://doi.org/10.1038/nature12804 PMID: 24256735

31. Süel GM, Garcia-Ojalvo J, Liberman LM, Elowitz MB. An excitable gene regulatory circuit induces transient cellular differentiation. Nature. 2006; 440(7083):545-550. https://doi.org/10.1038/nature04588 PMID: 16554821

32. Lord ND, Norman TM, Yuan R, Bakshi S, Losick R, Paulsson J. Stochastic antagonism between two proteins governs a bacterial cell fate switch. Science. 2019; 366(6461):116-120. https://doi.org/10.1126/science.aaw4506 PMID: 31604312

33. Salazar C, Brümmer A, Alberghina L, Höfer T. Timing control in regulatory networks by multisite protein modifications. Trends in Cell Biology. 2010; 20(11):634-641. https://doi.org/10.1016/j.tcb.2010.08.012 PMID: 20869247

34. Doan T, Mendez A, Detwiler PB, Chen J, Rieke F. Multiple Phosphorylation Sites Confer Reproducibility of the Rod9s Single-Photon Responses. Science. 2006; 313(5786):530-533. https://doi.org/10.1126/science.1126612 PMID: 16873665

35. Upadhyay A, Marzoll D, Diernfellner A, Brunner M, Herzel H. Multiple random phosphorylations in clock proteins provide long delays and switches. Scientific Reports. 2020; 10(1):22224. https://doi.org/10.1038/s41598-020-79277-z PMID: 33335302

36. Werner B, Dingli D, Lenaerts T, Pacheco JM, Traulsen A. Dynamics of mutant cells in hierarchical organized tissues. PLoS Computational Biology. 2011; 7:e1002290. https://doi.org/10.1371/journal.pcbi.1002290 PMID: 22144884

37. Balaban NQ, Helaine S, Lewis K, Ackermann M, Aldridge B, Andersson DI, et al. Definitions and guidelines for research on antibiotic persistence. Nature Reviews Microbiology. 2019; 17(7):441-448. https://doi.org/10.1038/s41579-019-0196-3 PMID: 30980069

38. Greulich P, Doležal J, Scott M, Evans MR, Allen RJ. Predicting the dynamics of bacterial growth inhibition by ribosome-targeting antibiotics. Physical Biology. 2017; 14(6):065005. https://doi.org/10.1088/1478-3975/aa8001 PMID: 28714461
39. Raj A, Oudenaarden Av. Nature, nurture, or chance: stochastic gene expression and its consequences. Cell. 2008; 135(2):216-26. https://doi.org/10.1016/j.cell.2008.09.050 PMID: 18957198

40. Acar M, Mettetal JT, van Oudenaarden A. Stochastic switching as a survival strategy in fluctuating environments. Nature Genetics. 2008; 40(4):471-475. https://doi.org/10.1038/ng.110 PMID: 18362885

41. Hobby GL, Meyer K, Chaffee E. Observations on the Mechanism of Action of Penicillin. Experimental Biology and Medicine. 1942; 50(2):281-285. https://doi.org/10.3181/00379727-50-13773

42. Bigger JW. Treatment of staphylococcal infections with penicillin—By intermittent sterilisation. Lancet. 1944; 2(6320):497-500. https://doi.org/10.1016/S0140-6736(00)74210-3

43. Himeoka Y, Mitarai N. When to wake up? The optimal waking-up strategies for starvation-induced persistence. PLOS Computational Biology. 2021; 17(2):e1008655. https://doi.org/10.1371/journal.pcbi.1008655 PMID: 33571191

44. van Gestel J, Nowak MA. Phenotypic Heterogeneity and the Evolution of Bacterial Life Cycles. PLoS Comput Biol. 2016; 12(2):1-23. https://doi.org/10.1371/journal.pcbi.1004764

45. Sahoo S, Subbalakshmi AR, Jolly MK. The fundamentals of phenotypic plasticity. In: Levine Herbert K P N V Jolly Mohit Kumar, editor. Phenotypic Switching—Implications in Biology and Medicine. Elsevier; 2020. p. 1-21.

46. Aguirre-Ghiso JA. Models, mechanisms and clinical evidence for cancer dormancy. Nature Reviews Cancer. 2007; 7(11):834-846. https://doi.org/10.1038/nrc2256 PMID: 17957189

47. Burrell RA, McGranahan N, Bartek J, Swanton C. The causes and consequences of genetic heterogeneity in cancer evolution. Nature. 2013; 501(7467):338-345. https://doi.org/10.1038/nature12625 PMID: 24048066

48. Cleary AS, Leonard TL, Gestl SA, Gunther EJ. Tumour cell heterogeneity maintained by cooperating subclones in Wnt-driven mammary cancers. Nature. 2014; 508(7494):113-117. https://doi.org/10.1038/nature13187 PMID: 24695311

49. Sisan DR, Halter M, Hubbard JB, Plant AL. Predicting rates of cell state change caused by stochastic fluctuations using a data-driven landscape model. Proceedings of the National Academy of Sciences of the United States of America. 2012; 109(47):19262-7. https://doi.org/10.1073/pnas.1207544109 PMID: 23115330

50. Miller AK, Brown JS, Enderling H, Basanta D, Whelan CJ. The Evolutionary Ecology of Dormancy in Nature and in Cancer. Frontiers in Ecology and Evolution. 2021; 9:676802. https://doi.org/10.3389/fevo.2021.676802

51. Carey JN, Mettetal EL, Roggiani M, Myers KS, Kiley PJ, Goulian M. Regulated Stochasticity in a Bacterial Signaling Network Permits Tolerance to a Rapid Environmental Change. Cell. 2018; 173(1):196-207.e14. https://doi.org/10.1016/j.cell.2018.02.005 PMID: 29502970

52. Süel GM, Kulkarni RP, Dworkin J, Garcia-Ojalvo J, Elowitz MB. Tunability and Noise Dependence in Differentiation Dynamics. Science. 2007; 315(5819):1716-1719. https://doi.org/10.1126/science.1137455 PMID: 17379809

53. Johnson CH, Zhao C, Xu Y, Mori T. Timing the day: what makes bacterial clocks tick? Nature Reviews Microbiology. 2017; 15(4):232-242. https://doi.org/10.1038/nrmicro.2016.196 PMID: 28216658

54. Eelderink-Chen Z, Bosman J, Sartor F, Dodd AN, Kovács ÁT, Merrow M. A circadian clock in a nonphosphosynthetic prokaryote. Science Advances. 2021; 7(2):eabe2086. https://doi.org/10.1126/sciadv.abe2086 PMID: 33523996

55. Helfrich-Förster C, Bertolini E, Menegazzi P. Flies as models for circadian clock adaptation to environmental challenges. European Journal of Neuroscience. 2020; 51(1):166-181. https://doi.org/10.1111/ejn.14180

56. Simpson GG, Dean C. Arabidopsis, the Rosetta Stone of Flowering Time? Science. 2002; 296(5566):285-289. https://doi.org/10.1126/science.296.5566.285

57. Song YH, Shim JS, Kinmonth-Schultz HA, Imaizumi T. Photoperiodic Flowering: Time Measurement Mechanisms in Leaves. Annual Review of Plant Biology. 2014; 66(1):1-24. https://doi.org/10.1146/annurev-arplant-043014-115555 PMID: 25534513

58. Ito-Miwa K, Furuike Y, Akiyama S, Kondo T. Tuning the circadian period of cyanobacteria up to 6.6 days by the single amino acid substitutions in KaiC. Proceedings of the National Academy of Sciences. 2020; 117(34):20926-20931. https://doi.org/10.1073/pnas.2005496117 PMID: 32747571

59. Doulcier G, Lambert A, Monte SD, Rainey PB. Eco-evolutionary dynamics of nested Darwinian populations and the emergence of community-level heredity. eLife. 2020; 9:e53433. https://doi.org/10.7554/eLife.53433 PMID: 32633717

60. Wang J, Zhang K, Xu L, Wang E. Quantifying the Waddington landscape and biological paths for development and differentiation. Proceedings of the National Academy of Sciences of the United States of America. 2011; 108(20):8257-62. https://doi.org/10.1073/pnas.1017017108 PMID: 21536909
61. Sasai M, Wolynes PG. Stochastic gene expression as a many-body problem. Proceedings of the National Academy of Sciences. 2003; 100(5):2374-2379. https://doi.org/10.1073/pnas.2627987100
PMID: 12606710

62. Ozbudak EM, Thattai M, Lim HN, Shraiman BI, Oudenaarden Av. Multistability in the lactose utilization network of Escherichia coli. Nature. 2004; 427(6976):737-740. https://doi.org/10.1038/nature02298
PMID: 14973486