Dynamic phenotypic heterogeneity generated by delayed genetic oscillations

Eukaryotes and prokaryotes exploit the ability of genetically identical cells to exhibit different phenotypes in order to enhance their survival. However, the mechanisms by which cells transition from one phenotype to another remain unclear. Canonical models of this dynamic posit that molecular fluctuations provide the noise that drives the cell out of one stable state and into another. Stochastic processes generated by canonical models should, therefore, be good descriptors of phenotype dynamics and between-state transitions should become more likely at greater noise amplitude, for instance at higher extracellular temperatures. To test these predictions, we observed temporal expression dynamics of the promoter of a flagellum gene, fliC, in a microfluidic device using Salmonella enterica serovar Typhimurium and green fluorescent protein (GFP). Our observations show that while cells can exhibit multistable phenotypes, including stable fliC-OFF and fliC-ON states characterised by low and high GFP levels, respectively, between-state transitions can exhibit oscillatory dynamics whose return statistics do not conform to canonical theories. For example, here the fliC-ON state was more frequent following a temperature increase. To better understand our data we developed different dynamical frameworks to predict fliC expression data. We conclude that a stochastic dynamical system tailored to the genetic network of fliC is better suited to our data than prior theories where dynamical features, like oscillations and pulsing, are driven by inevitable delays in the post-translational regulation of fliC. Thus, while transcriptional noise promotes phenotypic heterogeneity, as we show here, regular features like oscillatory heterogeneity can result from delays that fundamental molecular processes impose upon a cell’s gene regulatory architecture.


13
Phenotypic heterogeneity -the expression of different phenotypes by clonal individuals living in a homogeneous environment - Figure 1. A) Network diagram representing the interaction between the three-tier operon regulating the expression of the flagellum. B) Bifurcation phase diagram of a gene regulatory model showing that the system presents bistability for certain parameter values (in this case, α V , the rate of interaction between Y diV and FlhDC). Solid lines represent stable steady states, the dotted is an unstable branch and the shaded area the parameter region where the system presents bistability. C) Theoretical frequency distribution of the expression of fliC simulated using the stochastic model with 1000 cells. As expected in a bistable toggle-switch, the resulting histogram of gene expression presents bimodality. D) Stochastic simulation of the bistable toggle-switch model (black solid line represents the expression level of FliAZ, Ω = 0.05). Note how at t = 0 the state of the system is in the fliC-OFF state, until transcriptional noise drives the system towards the fliC-ON steady state. Genes are denoted with rectangles and proteins with circles. heterogeneity in flagellar gene expression. However, the observed phenotype distributions can also be described by different 102 underlying dynamics at the single-cell level. For instance, it has been reported that high-levels of noise in a graded switch 103 with a sigmoidal response can also produce bimodal distributions in gene expression data, 47 so too can complex non-linear 104 dynamics, notably oscillatory behaviour. 48,49 To avoid making inferences about the underlying gene regulatory network that 105 produces the observed population-level distribution of gene expression, we use time-resolved gene expression data taken from 106 individual cells grown in microfluidic devices observed under a microscope. To exploit our time-resolved data, we proceeded to analyse the temporal dynamics of gene expression. Figure 2 shows the 118 expression profiles of a mother cell and her daughters before they are washed out of the channel. (The cell with the oldest cell 119 pole resides at the bottom of each channel, and we refer to this cell as the mother cell. The daughters are young-pole cells that 120 emerge when the mother cell divides. The daughter cells are located above the mother cell in the channel and are eventually 121 removed from the microfluidic device). Notably, the expression profile of the daughter and mother cells are correlated at the 122 moment of division ( Figure 2C), suggesting that transitions between different states are not driven by stochastic partitioning 123 during cell division. Moreover, the GFP expression profile of the daughter cells seems to follow the dynamics of the mother 124 cell until both trajectories eventually diverge ( Figure 2B).

125
As predicted by the toggle-switch model, a subset of cells are not expressing fluorescence for the entire duration of the 126 experiment and therefore the state of the network is always in the fliC-OFF steady state in those cells see ( Figure 3A-B).

127
Conversely, cells expressing fliC exhibit complex temporal patterns that fluctuate between fliC-ON and fliC-OFF states ( Figure   128 3A-B). We used this dataset to estimate the expression profile of fliC (as measured by GFP intensity) and estimated a population-129 level histogram of GFP expression. As Figure 3C shows, the latter is a bimodal distribution of GFP intensity as predicted by the   164 Figure S4 shows numerical solutions of the delayed system of ordinary differential equations, with fixed delays but of 165 different durations. For small delays, τ, the system approaches the fliC-ON stable stationary point in an oscillatory manner, but 166 for larger delays, the dynamical system undergoes a bifurcation that creates dynamical oscillations in the expression of FliAZ.

167
The fact that solutions of time-delay differential equations exhibit oscillatory behaviour is unsurprising from a theoretical 168 perspective; it is well-known in dynamical systems theory that delays induce instabilities that may produce oscillations in implementing different values of Ω. 195 We then performed stochastic simulations with initial conditions in the fliC-OFF state and found, as expected, that the 196 mean escape time from the fliC-OFF basin of attraction correlates negatively with Ω ( Figure S5C). Indeed, noise-driven escape  Figure 4C).

201
To test this, we obtained single-cell mean fluorescent intensities in experiments performed at a range of temperatures.
202 Figure 5A shows the deviation from population-level mean fluorescence intensities exhibited by all cells as a stacked dataset.
suitably augmented with initial conditions and where k i denotes translation rate, γ i mRNA degradation rate, δ i protein Numerical simulations of the model presented in this paper were performed using standard numerical solvers of Matlab and 276 NumPy with parameter values described in Table S2.

278
We used a Gillespie algorithm to obtain numerical realisations of the stochastic model. Table 1  optimisation algorithm based on a minimal information criterion with n as the optimisation parameter showed that n = 7 was 296 the number of states that best approximated our data (see Figure S3A for an example). We then used the Causal State Splitting

297
Reconstruction algorithm (CSSR) 74 to estimate a n × n matrix of transition probabilities between the quantised states and their 298 quantised between-state trajectories. Finally, by initialising the system into any one of these discrete states and iterating the 299 resulting Markovian model, we could produce synthetic timeseries data with the same between-state transition probabilities as 300 the single-cell data obtained experimentally.

14/21
Transition probabilities (log10)  Figure S6. These are histograms of observed states that result when increasing noise intensity in a stochastic model based on a mathematical model of phenotypic 'particles' moving in a bistable, 2-well potential: note, because the time states spend close to each phenotype (at x = +1 and x = −1) reduces as noise increases (ordered from left to right) the system spends more time in intermediate states (near zero x = 0). Finally, as noise increases further (the rightmost plot) the histogram approaches a near-uniform distribution wherein the two phenotypes can no longer be distinguished in the data. These data should be contrasted against the experimentally observed outcomes from single-cell observations in Figure 5C that do not have this form. Simulations for this supplementary figure use the stochastic oscillator d 2 x/dt 2 +V (x) = δ dx/dt + σ dW where x(t) is said to be the phenotype or state, additive noise intensity σ is the temperature proxy, V (·) is a 2-well potential, dissipation is δ > 0 and W is Brownian motion (Matlab codes are supplied). Also see Figure S7 below for more details.  Figure S7. Exemplar stochastic trajectories (x(t)) realised from numerical solutions of d 2 x/dt 2 +V (x) = δ dx/dt + σ dW that form the basis of the histograms in Figure S6. They do not resemble the empirical fliC trajectories in Figure 3 and the main text explains that the statistics of return times between states is one way of demonstrating this lack of resemblance formally.