Hidden hysteresis – population dynamics can obscure gene network dynamics

Background Positive feedback is a common motif in gene regulatory networks. It can be used in synthetic networks as an amplifier to increase the level of gene expression, as well as a nonlinear module to create bistable gene networks that display hysteresis in response to a given stimulus. Using a synthetic positive feedback-based tetracycline sensor in E. coli, we show that the population dynamics of a cell culture has a profound effect on the observed hysteretic response of a population of cells with this synthetic gene circuit. Results The amount of observable hysteresis in a cell culture harboring the gene circuit depended on the initial concentration of cells within the culture. The magnitude of the hysteresis observed was inversely related to the dilution procedure used to inoculate the subcultures; the higher the dilution of the cell culture, lower was the observed hysteresis of that culture at steady state. Although the behavior of the gene circuit in individual cells did not change significantly in the different subcultures, the proportion of cells exhibiting high levels of steady-state gene expression did change. Although the interrelated kinetics of gene expression and cell growth are unpredictable at first sight, we were able to resolve the surprising dilution-dependent hysteresis as a result of two interrelated phenomena - the stochastic switching between the ON and OFF phenotypes that led to the cumulative failure of the gene circuit over time, and the nonlinear, logistic growth of the cell in the batch culture. Conclusions These findings reinforce the fact that population dynamics cannot be ignored in analyzing the dynamics of gene networks. Indeed population dynamics may play a significant role in the manifestation of bistability and hysteresis, and is an important consideration when designing synthetic gene circuits intended for long-term application.

The Lewis model for gene expression and its reduced, dimensionless form are 8 reproduced from the main text here: 9 dg dt = k 1 s 0 − k 2 g + k 2 3 g 2 k 2 4 + g 2 (S1) Dividing the non-linear term involving g on the RHS by k 2 4 : 10 dg dt = k 1 s 0 − k 2 g + (k 2 3 /k 2 4 )g 2 1 + (g 2 /k 2 4 ) Setting x = g/k 4 , we get: Dividing throughout by k 2 3 we get: , we reduce equation S1 to the 13 non-dimensionless form as shown below: S1.2 Analyzing the dynamical system 15 The stable points of Equation S5 occur when dx dτ = 0, or when x 2 1+x 2 = −(s − rx). Figure S1 Depending on the specific values of r and s, we can get between zero and 3 solutions.

18
When three solutions exist, as shown in Figure S1(A), the two outside solutions are 19 stable, while the middle solution is unstable, making the system bistable.

20
The boundary condition for bistability arises when the unstable solution moves towards 21 one of the stable points. When the two solutions coalesce, the slope of the two curves, as 22 well as the ordinates are equal. Mathematically: Solving these above equations we obtain the following expressions for r and s 24 parameterized in x: The above equations can be plotted as a parametric plot to visualize the dependence of 26 the system behavior on r and s as seen in Figure S2. The key conclusion from Figure S3 is that the system is capable of hysteretic behavior.  improved. Firstly, we can assume a bipartite degradation term that separately accounts 44 for decay and dilution. Moreover, since we re modeling our culture as a batch process 45 we can assume that the dilution rate decreases with time as the cell reproduction rate 46 decreases. The decay rate can be assumed linear with respect to the gene product, while the dilution rate can be approximated as an exponential decay.

48
Secondly, we can reinterpret the term x equation S5 to be the concentration of LuxR, 49 rather than GFP. The expression of GFP can have its own induction rate proportional to 50 x and degradation and dilution rates proportional to its own concentration. Thus the 51 equation in S5 instead becomes a system of coupled differential equations as follows: However this system is difficult to solve numerically; we weren't able to explore its 53 behavior due to the absence of solutions of this very stiff system. Dropping the 54 exponential dilution terms we get a system that can be solved, but one which can be 55 approximated satisfactorily with the much simpler form in S5. We therefore do not The coefficients α and β determine respectively how much the presence of the OFF  Neutral competition for resources only α < 1 β < 1 Mutualism / synergy between species α > 1 β > 1 Mutual inhibition α > 1 β < 1 Parasitism by the OFF phenotype at the cost of the ON phenotype α < 1 β > 1 Parasitism by the ON phenotype at the cost of the OFF phenotype The p ON and p OF F are state change probabilities for cells turning ON and OFF 74 respectively. The system of differential equations in Equation S13 is autonomous -the 75 slopes do not depend upon t, so we can visualize the behaviors of the solutions using a 76 vector field plot as shown in Figure S4.  Figure S4: A vector field showing the direction of the slopes for the system of differential equations shown in S13. Units are normalized. As long as K OF F > K ON . Two special situations are shown: The green arrows show the trajectory of equal but small initial concentrations of cells growing at the same level until N ON + N OF F = K ON . Subsequently the vector field shifts such that the N ON population goes to zero and the N OF F population approaches K OF F . The red arrow denotes the situation when an OFF phenotype appears in a steady state culture of ON phenotypes.
The vector plot above is insensitive to the growth rate ρ, which can be readily seen since S2.2 Sensitivity analysis for the competitive growth model 84 We solved the numerical differential equations shown in Equation S13 using parameter 85 estimates for ρ and K for the ON and OFF phenotypes obtained from Figure 5 in the 86 main text. The model sensitivity was visualized for changes in dilution rate, changes in 87 the stochastic switching probability as well as sensitivity to the initial ratio of ON and 88 OFF phenotypes. The dependence of the dilution rate is shown in the main text ( Figure   89 8). The following Figure S5 shows the sensitivity to switching probabilities and initial 90 distribution of the phenotypes. Increasing p OFF Figure S5: Numerical estimates of dependence of the ON fraction solely on A) the initial ON -OFF distribution, and B) variation in the probability of switching between the ON and OFF fractions. The model behaves as expected, i.e. the rate at which the ON population fraction decreases is directly proportional to the initial distribution of the ON population, and is indirectly proportional to the probability with which an ON phenotype will switch to OFF.

S3 Solving for the constants of the dynamical system 92
The experimental data for induction of the feedback circuit presented in Figure 5 of the 93 main text can be fitted to Equation S1. Figure S6 shows the optical density inducer level is measured in ng/ml and time is measured in hours.  Table S2: Units for the terms in Equation S1 .
Corresponding to the fitted values of the constants k 1 through k 4 for given above, we can

122
The value of r < 0.5 indicates that the system should produce irreversible hysteresis, and 123 for a induction level of s 0 < 2.4 ng/ml, the system will be bistable.

Normalized Cell Count
Group A (1:1), 22% ON Group B (1:10), 28% ON Figure S8: Flow cytometry histograms for cells from groups A and B at 0.1 ng/mL aTc. Cells from each dilution group exhibit hysteresis, but group B contains a higher proportion of ON cells than group A. However, the steady-state fluorescence level (mode) of the group A ON subpopulation is greater than the mode for group B. Thus, group A has a higher overall magnitude of expression as measured by relative fluorescence, but a lower proportion of ON cells, indicating that the dilution ratio used during resuspension may be affecting the steady-state expression level.