Bistability: Requirements on Cell-Volume, Protein Diffusion, and Thermodynamics
11 Bistability: Requirements on Cell-Volume, Protein Diffusion,and Thermodynamics
Robert G. Endres ∗ Department of Life Sciences & Centre for Integrative Systems Biology and Bioinformatics, London,United Kingdom* [email protected]
Abstract
Bistability is considered wide-spread among bacteria and eukaryotic cells, useful e.g. for enzyme induc-tion, bet hedging, and epigenetic switching. However, this phenomenon has mostly been described withdeterministic dynamic or well-mixed stochastic models. Here, we map known biological bistable systemsonto the well-characterized biochemical Schl¨ogl model, using analytical calculations and stochastic spatio-temporal simulations. In addition to network architecture and strong thermodynamic driving away fromequilibrium, we show that bistability requires fine-tuning towards small cell volumes (or compartments)and fast protein diffusion (well mixing). Bistability is thus fragile and hence may be restricted to smallbacteria and eukaryotic nuclei, with switching triggered by volume changes during the cell cycle. For largevolumes, single cells generally loose their ability for bistable switching and instead undergo a first-orderphase transition.
Introduction
Isogenetic cell populations show remarkable heterogeneity due to unavoidable molecular noise - bacteriaare either induced or uninduced to produce enzymes for utilizing a particular sugar [1], or enter cel-lular programs such as competence during starvation in a reversible, switch-like manner [2]. In higherorganisms examples of bistability are maturation in developing oocytes in
Xenopus frog embryos [3],Hedgehog signaling in stem cells [4], and phosphorylation-dephosphorylation cycles, e.g. as occurring inmitogen-activated protein kinase (MAPK) cascades [5]. Bistable pathway designs have also been exploredin synthetic biology [6–10]. In analogy with physical bistable systems such as ferromagnets, biologicalcellular systems can indeed exhibit hysteresis, indicative of a system’s memory of past conditions [1,6–8].Functional advantages of bistability include bet-hedging strategies, decision-making, specialization, andmechanisms for epigenetic inheritance, all increasing the species’ fitness [11,12]. However, these phenom-ena have mostly been described with deterministic dynamic models or well-mixed stochastic models. Itis unclear if bistability predicted by the deterministic model always corresponds to a bimodal probabilitydistribution in the stochastic approach [13]. Furthermore, the influence of slow protein diffusion and lo-calization inside the cytoplasm (bacteria) or nucleus (eukaryotes) is often neglected. Whether bistabilityis robust to such perturbations is unclear.The question of the role of reaction volume in well-mixed bistable chemical reactions has a longhistory, e.g. [13–17]. Particularly noteworthy, Keizer’s paradox says that microscopic and macroscopicdescriptions can yield different predictions [13]. In the macroscopic description the steady-state ( t → ∞ )is considered after taking the infinite volume limit ( V → ∞ ), while in the microscopic description theopposite order of limits is taken. Since the orders are not always interchangeable [18–20], unexpectedresults can occur. For instance, in the logistic growth equation species extinction occurs in the micro-scopic description, while the macroscopic description predicts a stable finite steady-state population [21].As a consequence, in bistable systems it is generally not possible to derive Fokker-Planck or Langevinequations that produces a behavior in accordance with the master equation [13, 22]. Derived potentialsdetermining the weights of the states are incorrect. More sophisticated approximations (or modifiedFokker-Planck approaches) are required to capture rare large fluctuations [22, 23], which ultimately de- a r X i v : . [ q - b i o . S C ] M a r termine the switching between states. However, the biological implications of these issues on cell-fatedecisions have been rather unexplored, with some exceptions [24, 25].Furthermore, two recent papers address the effects of diffusion on bistability and switching of states.Zuk et al. considered a one-dimensional (1D) and a hexagonal 2D lattice model [16], while T˘anase-Nicolaand Lubensky considered an 1D M -compartment model with hopping between the M compartments torepresent diffusion [17]. Based on their results, when the system size is small such systems are effectivelywell-mixed and transitions are driven solely by stochastic fluctuations in line with the well-mixed masterequation. However, when the system is spatially extended the more stable state spreads out in spaceand overtakes the more unstable state by the mechanism of traveling waves. Interestingly, in presence ofdiffusion the stability of steady states in the extended system is determined by the deterministic (mean-field) potential, which also describes the speed of the traveling waves. However, it is unclear if theseresults also hold in 3D, for small volumes comparable to nuclei and cytoplasms in cells, and using morerealistic particle-based approaches.Living cells are open molecular systems, characterized by chemical driving forces and free-energydissipation [26, 27]. Here, we map known biological bistable systems onto the well-characterized non-equilibrium biochemical Schl¨ogl model [14] (recently reviewed in [13]), allowing us to obtain analyticalresults for the well-mixed case. For slow diffusion we use stochastic spatio-temporal simulations. Inaddition to network architecture and strong thermodynamic driving away from equilibrium, we show thatbistability requires fine-tuning towards small cell volumes (or compartments) and fast protein diffusion(well mixing). Bistability is thus fragile and hence may provide upper limits on cell or nuclear sizes. Forincreasing volume, a separation of time scales occurs and switching does not only become infinitesimally(exponentially) rare but the weights of the states shift as well. Although states do not disappear perse, weights can disappear, leading effectively to monostability. Hence, single cells loose their ability forreversible bistable switching and instead undergo a first-order phase transition similar to mesoscopicphysical systems. Strict cell and nuclear size control may provide a protective molecular environmentfor bistability. Indeed, our analysis of previously published time-lapse movies of bacteria indicates thatvolume changes during cell growth and division may function as triggers for switching. Results
Mapping of bistable systems onto Schl¨ogl model
Bistability is driven by high-energy fuel molecules such as ATP and sources of precursor molecules [28].Here, we choose the self-activating gene, whose protein product binds its own promoter region to cooper-atively activate its own transcription as a dimer (see Fig. 1A, mRNA is not explicitly modeled here). Inaddition to ATP required for charging synthetase with amino acids and tRNAs, high-energy moleculesinvolved are nucleotide triphosphates during transcription and GTP during translation [29]. Additionally,we consider the phosphorylation-dephosphorylation cycle with the phosphorylation reaction catalyzed bykinase K and the dephosphorylation reaction catalyzed by phosphatase (inhibitor) I (Fig. 1B) [28]. Alsoin this case there is positive feedback from product P p to its production. The mean-field equations, givenby ordinary differential equations (ODEs), describe the temporal dynamics of the average protein concen-trations, valid in the limit of large volume (and hence large protein copy numbers). At steady state, whenall time-derivatives are zero, the equations produce a bistable bifurcation diagram for suitable parameterregimes (see Fig. 1C for a schematic). The control parameter ( x -axis) is a parameter of the model, e.g. arate constant, and the output ( y -axis) is the target protein or the phosphorylated protein concentration,respectively.The chemical reactions for the self-activating gene and the phosphorylation-dephosphorylation cyclewith their rate constants are shown in Figs. 1D and E, respectively. In Fig. 1D the first equation Figure 1. Mapping of bistable systems onto Schl¨ogl model. (A) Self-activating gene withcooperativity. (B) Phosphorylation-dephosphorylation cycle. (C) Schematic bifurcation diagram withbistable regime indicated by vertical dashed lines. (D, E) Chemical reactions corresponding to (A) and(B), respectively. (D) S is substrate (nucleotides for mRNA and amino acids for protein etc.) and P isprotein product. (E) Quantities I , K , P ( P p ), and P i are the inhibitor, kinase, (phosphorylated)protein, and inorganic phosphate, respectively. (F) Chemical reactions of Schl¨ogl model withconcentrations A and B adjustable parameters. For mapping reactions in (D) onto reactions in (F) genespecies needs to be absorbed into rate constants, and S and P identified with A/B and X , respectively.For mapping (E) onto (F) I , K , ADP, and P need to be absorbed into rate constants, and P p identifiedwith X , P i with A , and ATP with B .indicates basal expression, while the second indicates cooperative self-activation with S the substrate(e.g. nucleotides for mRNA and amino acids for protein production). In Fig. 1E, the substrate is givenby ATP, which is converted to ADP. P i is inorganic phosphate produced during dephosphorylation, whichagain is converted into ATP by the cell. Note, while the reverse reactions from protein to substrate S (Fig. 1D) and protein phosphorylation by I or protein dephosphorylation by K are extremely unlikely,they technically are nonzero and need to be included for thermodynamic consistency. Importantly, theindividual reactions can be mapped onto the well-characterized single-species Schl¨ogl model, in whichmolecular concentrations A and B are fixed to drive the reactions out of equilibrium.The mapping is justified based on the one-to-one correspondence of the molecular reactions (see Fig.1F). For this, however, to work the biological examples would need to be implemented by mass-actionkinetics instead of more realistic enzyme-driven kinetics. For instance, the self-activating gene mightbe implemented by dp/dt = a + bp / ( K + p ) − τ − p to describe cooperative self-induction with Hillcoefficient 2, protein life time τ , and additional parameters a, b and K . While for the self-activatinggene p -dependent production is to lowest order ∼ p similar to the Schl¨ogl model (with rate constant k − ), its reverse rate is assumed to be zero as the forward rate is highly driven by several enzymaticsteps. In contrast, in the Schl¨ogl model the reverse rate is assumed to be non-zero (with rate constant k +2 ). Similarly, while degradation in the Schl¨ogl model has a reverse rate (“accidental” productionfrom constituents via rate constant k +1 ) degradation in gene regulation is either implemented by activedegradation or dilution during cell division, both of which have negligible reverse rates. As a result, themacroscopic equation provided by the Schl¨ogl model is a third-order polynomial with rather large reversereactions due to the absence of enzymatically driven reactions. Macroscopic perspective
In the limit of large volume and hence large molecule numbers, the Schl¨ogl model is described by ODE dxdt = − k +2 x (cid:124) (cid:123)(cid:122) (cid:125) w +2 + k − Bx (cid:124) (cid:123)(cid:122) (cid:125) w − − k − x (cid:124) (cid:123)(cid:122) (cid:125) w − + k +1 A (cid:124) (cid:123)(cid:122) (cid:125) w +1 (1)with x the molecular concentration. Once this limit is taken, time can be sent to infinity. The result-ing steady-state bifurcation diagram is shown in Fig. 2A for standard parameters (see Materials andMethods), with concentration B chosen the control parameter. Two saddle-node bifurcations (SNs) in-dicate the creation/destruction of steady states, with a range of bistability described by ∆ B in between.However, the macroscopic perspective makes no prediction about the relative stability of the two stablesteady states (black and blue curves with the unstable steady state shown in red). In particular, dotransitions simply become rarer with increasing volume so that the state attractors become increasinglydeep but the relative weights of the states intact, or do the weights of the states change as well, leadingeffectively to loss of bistability? Furthermore, is there a thermodynamic selection principle for the moststable steady state? According to the second law of thermodynamics entropy is maximized in a closedsystem at equilibrium. Does a similar extremal principle hold for nonequilibrium steady states? The rateof entropy production describes how much heat is dissipated per time at steady state (and hence is alower bound of how much energy is consumed to maintain the state [13, 30]): dsdt = (cid:88) i =1 ( w + i − w − i ) log (cid:18) w + i w − i (cid:19) ≥ . (2)Intuitively, the entropy production is the net flux (difference between forward and backward fluxes)times the difference in chemical potential between products and educts (log term), summed over allthe reactions. Eq. 2 thus effectively describes how quickly the maximum entropy state is reached, ifleft to equilibrate. Prigogine and co-workers argued for a minimal rate of entropy production, at leastnear equilibrium [31], while others argued for maximal rate of entropy production [32, 33]. Fig. 2Bshows the macroscopic rate of entropy production. The red arrow indicates equilibrium at which therate is zero. Notably the high state always has the higher entropy production. This can easily beunderstood [34] since overall in the Schl¨ogl model A is converted to B (and vice versa, see Fig. 1F).Hence, ds/dt = ∆ G/T · dA/dt x with ∆ G the change in free energy for the overall reaction, T thetemperature of the bath and dA/dt a linear function of x . As a result, if minimal entropy production isthe rule, then the low state should be selected, while maximal entropy production would dictate that thehigh state is more stable. As eigenvalues of the Jacobian only contain information about local stability ofa fixed point and not about global stability across multiple fixed points [17, 35], further discussion needsto be postponed until the next section. Fig. 2C summarizes the phase diagram (see S1 Text), showingmonostable (low or high state only) and bistable regions in line with a cusp catastrophe. Limit V → ∞ ,shown by red lines, is relevant for the macroscopic description. Microscopic well-mixed perspective
When first taking the long-time limit for a fixed finite volume to obtain the steady-state distributionand then increasing the volume, we obtain a very different picture of bistability. Assuming a well-mixedmicroenvironment and thus neglecting diffusion (illustrated in Fig. 3A), we can employ the one-stepchemical master equation to describe the probability distribution in time ddt p ( X ; t ) = ± (cid:88) i = ± [ W i ( X − ν i | X ) p ( X − ν i ; t ) − W − i ( X | X − ν i ) p ( X ; t )] (3)with X the molecule copy number and W +1 ( X | X +1) = k +1 AV, W − ( X | X −
1) = k − X, W +2 ( X | X −
1) = k +2 X ( X − X − /V , and W − ( X | X + 1) = k − BX ( X − /V the volume-dependent transitionrates. In Eq. 3, the sum is over both the forward ( i = +1 , +2) and backward ( i = − , −
2) reactions with ν ± = ± ν ± = − ν ± [30]. Using this description, we first simulate the master equation using theGillespie algorithm [36] and confirm stochastic switching between low and high stable states (Fig. 3B).However, at steady state setting dp/dt = 0, the probability distribution can analytically be derived usinga recursive relation leading to p ( x ) = N ( x ) exp[ − V Φ( x )] with potential [22]Φ( x ) = x (ln x −
1) + x ln (cid:18) k +2 x + k − k − Bx + k +1 A (cid:19) +2 (cid:115) k − k +2 arctan (cid:32)(cid:115) k +2 k − x (cid:33) − (cid:115) k +1 Ak − B arctan (cid:32)(cid:115) k − Bk +1 A x (cid:33) (4)and volume-independent prefactor N ( x ) = k +2 x + k − Z √ x [ x + k +1 A/ ( k − B )] , (5)where Z is a normalization constant, x = X/V and p ( x ) = V p ( X ) (see S1 Text for details). Byconstruction, the stochastic potential also has minima (a maximum) at the stable (unstable) deterministicsteady states (state) with d Φ dx = − ln (cid:18) w +1 + w − w − + w +2 (cid:19) (6)equal to zero at steady state ( w +1 + w − = w − + w +2 ) and d Φ /dx having the correct sign (see S1Text for details and S1 Fig. for a plot of Φ( x )). Fig. 3C shows indeed that the resulting distribution of x has the expected bimodality.We are now in a position to address the relative stability of the steady states, in particular of the twostable states. Fig. 3D shows the probabilities evaluated at the deterministic steady states, indicating acrossing of the stable states (exchange of stability) with the probability of the unstable (metastable) stateconsistently below the probabilities of the stable states. A more precise picture emerges when plottingthe transition rates for switching between the stable steady states in Fig. 3E [22], showing coexistenceof the two stable states at B ∼ .
7. As derived in S1 Text, the rates depend exponentially on thevolume (as expected). However, due to normalization and the volume-independence of the prefactor, themore stable of the two becomes increasingly selected for larger and larger volumes, leading effectively tomonostability. Fig. 3F shows that a Maxwell-type construction (MC) is required to establish the point ofstability exchange, well known from the classical Van der Waals gas (see S1 Text for details) [13–15]. Sincethe two states have different entropy productions (Fig. 2B, which can also be confirmed by calculatingthe microscopic entropy production defined in S1 Text), we obtain a discontinuity at this point, indicativeof a first-order phase transition. Fig. 3G shows indeed a sharpening of the molecular fluctuations at thecritical point for increasing volume. Hence, mesoscopic cells can loose their ability for bistability ( i.e. abimodal distribution) with increasing volume.The strong volume-dependent of bistability can also be seen in the phase diagram in Fig. 2C (see [37]).For small volumes (black lines) the region of bistability can significantly deviate from the correspondingregion in the macroscopic limit (red lines). For instance, a point in parameter space with strong bistabilityin the microscopic system ( B = 3 . V = 10) is borderline bistable in the macroscopic limit (cf. Fig.3C for B = 4 . Microscopic perspective with diffusion
Diffusion introduces inhomogeneous distributions of molecules, with diffusion particularly slow in thecrowded intracellular environment (Fig. 4A). For this purpose we turn to the stochastic
Smoldyn simu-lation package for implementing particle-based reaction-diffusion systems in a box (Fig. 4B; see [38] andMaterials and Methods for further details). The third-order reaction (see Fig. 1F) needs to be convertedinto two second-order reactions since no two events can exactly occur at the same time. (We call thismodel the generalized Schl¨ogl model.) This conversion requires introducing of a dimer species X withadditional rate constants k +3 and k − as illustrated in Fig. 4C. For k +3 = k − the steady-state valuesremain unchanged in the macroscopic limit (see S1 Text). For reasonable diffusion constants (see Mate-rials and Methods for parameter values), we indeed observe stochastic switching, resulting in a bimodaldistribution for species X (Fig. 4D). We then compared Smoldyn simulations in detail with Gillespiesimulations of the generalized and conventional reaction systems, including convergence for rare stateswith increasing simulation time, as well as effects of diffusion and dimerization reactions on bimodaldistribution (see S1 Text and S2-S4 Fig.). From these tests we conclude that
Smoldyn simulations of thegeneralized system accurately produce bistable behavior, allowing us to study the effects of diffusion andvolume on bistability.Decreasing the diffusion constants of both molecular species by an order of magnitude, suitable formacromolecular complexes or membrane-bound proteins [39], leads to strongly fluctuating molecularconcentrations (illustrated by the molecule cluster enclosed by red dashed line in Fig. 4B) and reducedmolecule numbers in the high state (Fig. 4E). When instead increasing the reaction volume by just afactor 2, the high state is strongly induced (Fig. 4F). This result resembles the destruction of bistabilityobserved in the macroscopic limit (Fig. 3D-F).In Fig. 4A-F the molecules are able to react anywhere in the reaction volume. However, in cells, e.g.for a self-activating gene, transcription occurs localized at the DNA molecule (Fig. 4G). To investigate theeffect of localization on bistability we use a spherical cellular compartment (representing e.g. a bacterialcell or a eukaryotic nucleus) in which we introduce a small cylinder to represent the DNA molecule. Theproduction can only occur in this cylinder (Fig. 4H). In contrast, degradation can occur anywhere inthe cellular compartment. Fig. 4I shows that bistability is destroyed with localized production, evenfor drastically increased production rates and diffusion constants, which would easily produce bistabilityunder well-mixed circumstances. The broad distribution in Fig. 4I may thus be caused by strong localfluctuations in molecule number (illustrated by molecule cluster enclosed by red dashed line in Fig.4H). Note that the appearance of DNA as a single copy is markedly different from the conventional orgeneralized Schl¨ogl model, in which the molecule numbers scale with volume. Next, we will explore thereasons for the breakdown of bistability with inhomogeneity.Fig. 5 shows a systematic exploration of bistability from diffusive
Smoldyn simulations, conductedsimilar to Fig. 4D. Fig. 5A shows little evidence of bistability with the system either in the low or highstate. Diffusion causes strong fluctuations in molecule numbers (and hence clustering) as as demonstratedby the radial pair-correlation function g ( r ) in Fig. 5B (not to be confused with the spike in fluctuationsat the critical point of the well-mixed system in Fig. 3G). For small molecule-molecule distances r , weobtain g ( r ) (cid:29)
1, which is the larger the slower diffusion. In contrast, random distributions of moleculesdo not show clustering. While the next section investigates the role of such fluctuations in the loss ofbistability, our findings are summarized in Fig. 5C, which shows the bistable range ∆ B for the well-mixedcase and the inhomogeneous case with finite diffusion constants. Here the system is considered bistable ifsimulations started from low and high states exhibit at least one reversible switch within the simulationtime (see figure caption for additional details). The narrow range, especially for finite diffusion constants,suggests that bistability is a fragile property, which needs protection. Indeed, bacterial cell volume andnuclear volume in eukaryotic cells are tightly regulated (e.g. nuclear volume does not simply scale withDNA content [40]). When converting to physical units, our predicted bistability regions fall nicely intoexperimentally observed cell volumes (shaded areas in Fig. 5C). Importantly, as volume varies duringcell growth and division, such changes in volume may function as a pacemaker or trigger for phenotypicswitching. Mechanism of bistability reduction by diffusion
Increasing the volume shifts the weights of the states leading to an effective loss of bistability (although theminima of the stochastic potential coincide with the deterministic model for sufficiently large volumes).How does slow diffusion affect bistability? There are two main potential reasons for the reduction ofbistability with diffusion: (1) Diffusion may increase the barrier of switching so that bistability is harderto achieve or observe, both in simulations and experiments, or (2) diffusion may destabilize one of thestable states. In these mechanisms local fluctuations in molecule numbers may play a role as well, e.g. byintroducing damaging heterogeneity or by nucleating traveling waves so that the more stable state canspread effectively and surpass the unstable state.To rule out (1) longer and longer simulations can be conducted to guarantee convergence. S2 Fig.shows indeed that simulations are well converged, even for weakly populated states. This shows thatdiffusion does not significantly change the barrier height. To investigate (2) we use the method from [41]to renormalize the second-order rate constants of the generalized Schl¨ogl model ( k ± , k ± ; cf. Fig. 4C)by diffusion (see Materials and Methods). This allows us to effectively include diffusion in the well-mixed model without having to conduct particle-based simulations. Fig. 6A shows that histograms fromGillespie simulations of the well-mixed Schl¨ogl model with renormalized reactions match well resultsfrom Smoldyn simulations (small Kullback-Leibler divergence). In contrast, Gillespie simulations withoutrenormalized reactions do not match well. In particular, without renormalization the switch to the highstate occurs at smaller B values. Thus, achieving bistability is easier without diffusion as it requiresless thermodynamic driving. These results are summarized in Fig. 6B by the bifurcation diagram of themacroscopic model of the conventional Schl¨ogl model (Eq. 1) with renormalized rate constants k ± (note k ± do not affect the steady-state probability distribution as they are equal, see S1 Text). Specifically,this figure shows a significant delay in achieving bistability with increasing B value. Completely remov-ing first and second terms in macroscopic Eq. 1 leads to the complete collapse of bistability and a trulymonostable state around x = k +1 A/k − ≈ .
17 in line with simulations (see S5 Fig.).What role might the fluctuations observed in Figs. 4B and 5B play? Following ideas from extendedbistable spatial systems [16, 17], fluctuations may nucleate traveling waves, which then spread by dif-fusion. Although our main interest are small systems most relevant to cell biology, we extended thesimulation box in one of the spatial directions (Fig. 7A). Kymographs from simulations with standardparameters, run for different B values, show the spreading of the more stable state when initially startedin the unstable state. Near co-existence at B ∼ .
7, traveling waves exist which do not change the statepermanently, but ripple through the box. Although wave velocities can technically be obtained from theslope in the kymographs they are highly variable and hard to determine objectively due to small moleculenumbers.Taken together, slow diffusion makes reaching bistability harder as second- (and higher-order) reac-tions are impaired - molecules have difficulties encountering each other to produce nonlinear behavior.Fluctuations may lead to traveling waves in more extended spatial systems, which provides a mechanismfor the more stable state to overtake the less stable state.
Experimental prediction on switching with cell-volume changes
Our models make strong predictions on the effect of cell volume. One obvious prediction is that when asystem is tuned towards the bistable regime (which becomes harder and harder to achieve for increasingvolume), switching between the two states becomes increasingly rare. This is the well-known transitionfrom the stochastic to the macroscopic, deterministic limit, and was recently demonstrated using time-lapse microscopy. In budding yeast (
Saccharomyces cerevisiae ) the variability in the G1 phase (i.e. thetime from division to budding) is reduced with increased ploidy (copies of chromosomes) [42]. Similarly,the switching time for turning the Pho starvation program off under reversal of phosphate limitationis reduced with ploidy [25]. As the volume also scales with ploidy, the protein concentrations stayapproximately constant, thus reducing cell intrinsic noise and hence stochastic transitions between cellularstates.Our results, however, are more specific. They suggest increased unidirectional switching and hencemonostability and decision-making in growing and dividing cells. In fact, growth towards cell divisionleads to a volume increase by a factor two, which may cause cells to select a steady state (see Fig. 4F).Hence, for cellular parameters below the critical point, the low state will be selected, while for cellularparameters above the critical point, the high, induced state will be selected. Fig. 8A,B show the expres-sion of bistable reporters as a function of time, specifically LacY in
E. coli [43] and ComK (or ComG) in
B. subtilis [2]. The latter is indicator for competence (while competence is strictly speaking an excitatorypathway, the core module of ComK is bistable with an exit mechanism based on ComS [2]).Fig. 8B shows that switching to the high state appears during growth, while switching to the lowstate occurs immediately after cell division when the volume has shrunk suddenly to half its value.Note it is unlikely that the rise (drop) in fluorescence intensity is simply due to switching induced bygene duplication (halving) as the concentrations stay roughly constant due to accompanying volumechanges. Furthermore, the ratio of chromosome-replication and cell-division times is known to be about2:1 [44, 45]. Since, by visual inspection, cell division takes about 10% of the cell-cycle time ( T c ) in thetime-lapse movies, chromosome replication takes about 20% of it. This duration is short compared tothe rise-in-intensity phase, suggesting a different mechanism. Switching is still stochastic as shown bythe two yellow daughter cells - one induces the lac operon, while the other does not. In support for ourproposed scenario, spontaneous switching is extremely rare (for the lac operon estimated to be around0 .
004 per cell cycle in presence of 40 µ M inducer TMG [46]). Hence, volume changes during cell growthand division may instead be the main drivers, like a pacemaker, for switching.
Discussion
We presented a nonequilibrium thermodynamic model of bistability, relying on molecular stochasticityand chemical energy for switching and decision-making. To cover a large class of bistable systems, includ-ing self-activating genes with cooperativity and phosphorylation-dephosphorylation cycles, we mappedminimal models for these onto the well-characterized nonequilibrium Schl¨ogl model. Bistability and itshallmark of hysteresis are generic behaviors that are the same from one system to the next regardlessof details. Indeed, this property is shared with ferro-magnets and mutually repressing genes (toggleswitch) [10, 47]. Our approach is markedly different from recent deterministic approaches to postulatemultistability in signaling cascades, which neglect the physical effect of cell volume and molecular diffu-sion [48]. Deterministic approaches often predict complex dynamics with multiple attractors. However,when the volume is sufficiently large, such behaviors can disappear. Not only does switching becomeincreasingly rare, but also the weights shift and ultimately favor one of the states. Hence, bacterial cellsand eukaryotic nuclei, and cell compartments in general, may represent protectorates of complex bi- andmultistable behavior [47]. In contrast, mesoscopic cells are “boring”, unable to display complex behavior.Slow diffusion, caused by molecular crowding and localization, is a killer of bistability and cellsneed to deal with this issue. This is because slow diffusion selectively penalizes second- and higher-order reactions and hence nonlinearity. Consistent with our study, ultrasensitivity in MAPK cascades isdestroyed for slow diffusion due to rebinding of enzymes to their substrate [49], stressing the fundamentalimportance of diffusion in theoretical predictions of bistability. Phase domains and their movement arewell known from the Ginzburg-Landau equation for phase transitions - this equation is in fact similarto the Schl¨ogl model with diffusion (albeit in absence of stochastic effects). How can cells cope withthe negative effects of diffusion? While adjustment of diffusion constants is difficult [50], cells could usesmall transcription factors to speed up diffusion. Up to about 110 kDa, the mean diffusion coefficientfalls close to the Einstein-Stokes prediction for a viscous fluid [50]. This suggests that proteins up tothis size do not encounter significant diffusion barriers due to macromolecular crowding or a meshwork ofmacromolecular structures in the cytoplasm. Indeed, the repressor LacI of the
E. coli lac system, masterregulator ComK of the
B. subtilis competence system, and transcription factor Gal80 of the gal systemin budding yeast are only 38 .
6, 22 .
4, and 48 . µm /s (based on scaling relation in [51]). Another option for the cellis to tune the viscosity of its cytoplasm below a glass-transition point where metabolism-driven activemixing produces superdiffusive environments [52]. Note, however, that cells need to protect themselvesfrom sudden shrinkage of the cytoplasm or drastic increases in concentration by osmotic shock.Due to small volumes and slow diffusion in cells, bistability can only occur in a narrow range ofparameter space and thus may require fine-tuning [53] or a pacemaker. Consistent with this notion, ouranalysis of time-lapse microscopy movies shows that volume changes during the cell cycle may triggerswitching events (Fig. 8). Such assistance might be necessary since spontaneous switching can beextremely rare, likely caused by rare bursts in gene expression [43, 46, 54]. For instance, the switchingrate of lac system was estimated to be only 0 .
004 per cell cycle (in presence of 40 µ M inducer TMG) [46],and diauxic shifts take on average 2 hours [51] (see [9, 55, 56] for additional examples of slow switchingmuch beyond the cell cycle period). Taken together, switching may hence be more likely to occur via athresholding mechanism [53] as implied by the Maxwell-type construction. To clarify the details of thetrigger mechanism, more experimental investigation will be needed using modern microfluidic designs forcontinuous imaging of cells over very long times.While some of the issues raised here have individually been discussed for the Schl¨ogl model before[13–17], this has hardly been done in the context of biology. In particular, bistability and diffusion restrictvolumes to the sizes of bacteria or eukaryotic nuclei, and volume changes during the cell cycle may beexploited to robustly change cellular states. Our particle-based simulations focus on small 3D volumes,most relevant to cellular nuclei or cytoplasms, and hence are markedly different from recent studies [16,17].The latter addressed role of volume and diffusion on bistability in extended 1 and 2D lattice models,respectively. Our loss of bistability for slow diffusion is largely determined by renormalization of second-order rate constants, which ultimately leads to a collapse of the system to the low state for very smalldiffusion constants. When extending the system in one of the spatial dimensions, we observe the onsetof traveling waves, which may ultimately determine switching rates for even larger systems. Specifically,for extended systems the wave velocity is determined by the deterministic (and not the stochastic)potential [16,17]. Such quasi-1D extended states may biological be relevant in filamentous bacterial cells.Traveling waves are indeed observed in
Smoldyn simulations of the Min-system, a small biochemicalpathway which allows cells to determine their middle for accurate cell division [57].0Bistability is fascinating due to its connection with nonequilibrium physics, first-order phase transi-tions, and decision-making in cells. Our results show that volume shifts the weights of the states relativeto each other but not the steady-state values directly. Below a critical value the low state is selected,while above a critical value the high state is selected. In contrast, diffusion shifts the steady-state values .For sufficiently slow diffusion, only the low state survives. While widely studied there are a numberof open questions. One pressing question is whether epigenetic information is inherited in the spirit ofHopfield’s content-addressable memory [58]. Since the seminal work by Novick and Weiner in 1957, suchinheritance seems indeed to apply to the Lac system [59], while Fig. 8 shows that the daughter cells donot generally inherit the competence state (this might be due to the protease MecA [2]). Furthermore, cellvolume changes and their effect on bistability also have many biomedical implications. Examples includeviruses such as bacteriophage lambda [55] and HIV [60], as well as pancreatic β cells, responsible for glu-cose sensing and insulin production. These cells undergo large size changes, e.g. during pregnancy [61].Thermodynamics may shed new light on their regulatory mechanisms. Materials and Methods
Schl¨ogl model
In 1972 Schl¨ogl proposed two chemical reaction models for nonequilibrium phase transitions [14]. Oneexample shows a phase transition of first order, while another shows a phase transition of second order.When diffusion is included in the the first-order transition, coexistence of two phases in different spatialdomains may occur. For spherical domains the coexistence indicates the onset of the transition similarlyto the vapor pressure in droplets or bubbles. The volume dependence has been discussed early, e.g.in [37]. The related Keizer’s paradox has mostly been discussed in context of Schl¨ogl’s first model, butalso for the logistic growth equation [21], showing its relevance to a wide sectrum of systems. In thesesystems, large rare fluctuations have severe consequences. Keizer’s paradox says that deterministic andstochastic approaches can lead to fundamentally different results. In particular, the deterministic modelconsiders the infinite volume limit ( V → ∞ ) before considering the steady-state limit ( t → ∞ ). Hence,no transitions between steady states are allowed. The state of the system only depends on the initialcondition, which seems unphysical. In contrast, the stochastic model, by taking the steady-state limitfirst, can always settle in its lowest state, and thus is ultimately favored over the deterministic approach.Gaspard [30] and later Qian [13] took the perspective of open chemical system, and analyzed the secondSchl¨ogl model in terms of fluxes and entropy production. The Schl¨ogl model can be considered thesimplest bistable system but has not yet been verified or implemented experimentally by suitable chemicalreactions. Implementation and parameters
The chemical reactions of the Schl¨ogl model can be found in Fig. 1F. Macroscopically (for an infinitevolume) this model can be described by ODE Eq. 1. For finite volume but infinitely fast diffusion (wellmixed case) the master equation (Eq. 3) can be used or Gillespie simulations [36]. For the solution ofthe master equation and a derivation of the transition rates see S1 Text. The macroscopic (microscopic)entropy production formula is provided in Eq. 2 (in S1 Text). The chemical reactions for the generalizedSchl¨ogl model are given in Fig. 4C with further details provided in S1 Text.Stochastic spatio-temporal simulations with diffusion were conducted with
Smoldyn software version2.28 as described in [38]. Briefly, it is a particle-based, fixed time-step, space-continuous stochasticalgorithm for reaction-diffusion systems in various geometries based on Smoluchowski reaction dynamics[62]. Both simulations and Smoluchowski theory only apply to reactions up to second order. Givenreaction rate constants, diffusion constants, and time step,
Smoldyn determines reaction radii, i.e. binding1and unbinding radii. The binding radii correspond to the encounter complex, formed by diffusion. Onceformed, the reaction occurs. Intrinsic rate constants are strictly constant, i.e. independent of diffusion,for low particle densities and activation-limited reactions (see manual for details). This is checked andconfirmed by
Smoldyn at the beginning of each simulation. Under these conditions, the effects of rateconstants and diffusion constants on bistability can independently be explored.The renormalization of the second-order rate constants to effectively include diffusion into the well-mixed generalized Schl¨ogl model is done following [41]. Using k D = 4 πσD to describe the encounter oftwo molecules by diffusion, the following rate constants are obtained k (cid:48)± i = k ± i k D k + i + k D (7)with i = 2 , k ± i comparable to Smoldyn ’s intrinsic rate constants, product of cross section σ , and average diffusion constant D set to 0 . k +1 A = 0 . , k − = 3, and k +2 = k − = 1. Only concentration B , diffusion constants,and volume V were varied as indicated in figure captions. For Smoldyn simulations we additionally used k +3 = k − = 1 (Fig. 4D-F) and V = 2 . V DNA = 1 . , D = 30 ( X ) and 10 ( X ), k = k = 50 and B = 50 (Fig. 4I). Simulation time was generally 10 ,
000 unless specified differently. To convert to unitsin Fig. 5C, we set length and time scales to respective µm and s . We further express concentration in[nM], with 1nM corresponding to 1 (1000) molecules in a typical bacterial (eukaryotic HeLa) cell, andvolume as V = η µm with η a scaling number of order 1. Rate constants then become k +1 A = 5 / (6 η )[nM/s], k − = 1 [ s − ], k +2 = B [ s − ], k − = k +3 = 6 η/
10 [(nM s) − ], and k − = 1 [ s − ]. For geneexpression, k +1 A corresponds to typical basal expression rates found in bacteria [63]. Simulation analysis
The radial pair-correlation functions in Fig. 5B were calculated from 10 simulation snap shots of monomer X positions r i in a 3D cube of volume V , using g ( r ) = VN ( N −
1) 14 πr a N (cid:88) i =1 N (cid:88) j (cid:54) = i =1 I ij ( r − a < || r j − r i || ≤ r ) (8)for r ≥ a with a the mesh spacing, assuming periodic boundary conditions. I is either one if its argumentis true or zero if false. N is the current monomer number. For plot g ( r ) was then averaged over snap-shots for different N . Note since N varies, g ( r ) can systematically deviate from 1. As a control, randompositions were produced for which g ( r ) ∼ D KL = (cid:88) i p R ( x i ) ln (cid:20) p R ( x i ) p T ( x i ) (cid:21) , (9)where p R ( x i ) and p T ( x i ) are the reference and test distributions, respectively. The sum in Eq. 9 is overbins in x space. Image analysis
Fluorescence intensities of bacterial cells were extracted from time-lapse movies of fluorescence microscopy(
E. coli from [43] and
B. subtilis from [2]). The active-contour method from [64] as an
ImageJ plugin wasused to detect the cell boundaries as ridges, and pixel intensities inside of the contours were collected2using a simple custom written
Mathematica code. The background intensities including the phase contrastintensities were subtracted and only the intensities from the fluorescence channels were plotted in Fig.8B. Intensity density plotted in S6 Fig. is calculated as total intensity of a cell divided by the cell currentfocal area.
Acknowledgments
We thankfully acknowledge Hong Qian for constructive comments on the manuscript, and Nikola Ojkicfor help with the image analysis. Financial support was provided by Leverhulme-Trust Grant N. RPG-181 and European Research Council Starting-Grant N. 280492-PPHPI. This work was initiated at theworkshop Information, probability and inference in systems biology (IPISB 2013) in Edinburgh.
Supporting Information Legends
S1 Text:
Supporting text with mathematical derivations and additional explanations.
S1 Figure:
Comparison of deterministic and stochastic potentials. For concentration x the deterministicpotentials Ψ( x ) in green is calculated with Eq. 44 in S1 Text, while the stochastic potential Φ( x ) in blueis calculated with main-text Eq. 4 (Eq. 19 of S1 Text). (A) Arrow indicates that for B = 3 . B = 3 . B = 3 . B = 4 . S2 Figure:
Bimodal distribution with rare high state (low weight) for B = 3 . Smoldyn simu-lations converge for increasing simulation time. Remaining parameters were chosen as in Fig. 4D with x the monomer concentration. (A-C) Simulation time is increased from t = 1 , , , Smoldyn simulation for t = 50 ,
000 as reference distribution (see S1 Text for details).
S3 Figure:
Comparison of
Smoldyn simulations and Gillespie simulations of generalized Schl¨ogl modelfor increasing diffusion constant D for B = 3 . t = 10 , x the monomer concentration. (A) Gillespie simulation. (B) Smoldyn simulation for D = 3 ( X ) and 1 ( X ). (C) Smoldyn simulation for faster diffusion with D = 15 ( X ) and 5 ( X ). (D)Kullback-Leibler divergence between each of the two Smoldyn simulations and Gillespie simulation asreference distribution (see S1 Text for details).
S4 Figure:
Distributions from
Smoldyn simulations become increasingly similar to results from Gillespiesimulations of conventional Schl¨ogl model for increasing dimerization rate constants k +3 = k − as indi-cated, as well as fast diffusion (parameters as in S3C Fig.) with enlarged B = 3 .
7. (A) Gillespie simulationfor B = 3 .
2. (B)
Smoldyn simulation for k +3 = k − = 1. (C) Smoldyn simulation for k +3 = k − = 3.(D) Kullback-Leibler divergence between each of the two Smoldyn simulations and Gillespie simulationas reference distribution (see S1 Text for details).
S5 Figure:
Drastic reduction of diffusion collapses bistabiliy to predicted low monostable state fordifferent B values as indicated (see main text section “Microscopic perspective with diffusion”). Scalingfactor multiplies set of diffusion constants, D = 3 ( X ) and 1 ( X ). Remaining parameters as in S3B Fig.3 S6 Figure:
Switching may be triggered by cell-volume changes. Image analysis similar to Fig. 8Bbut with total intensity normalized by cell area to provide the intensity density.
References
1. Ozbudak EM, Thattai M, Lim HN, Shraiman BI, van Oudenaarden A. Multistability in the lactoseutilization network of
Escherichia coli . Nature, 2004; 427: 737-740.2. S¨uel GM, Garcia-Ojalvo J, Liberman LM, Elowitz MB. An excitable gene regulatory circuit inducestransient cellular differentiation. Nature, 2006; 440: 545-550.3. Xiong W, Ferrell Jr JE. A positive-feedback-based bistable ’memory module’ that governs a cell fatedecision. Nature, 2003; 426: 460-465.4. Lai K, Robertson MJ, Schaffer DV. The sonic hedgehog signaling system as a bistable genetic switch.Biophys J, 2004; 86: 2748-2757.5. Markevich NI, Hoek JB, Kholodenko BN. Signaling switches and bistability arising from multisitephosphorylation in protein kinase cascades. J Cell Biol, 2004; 164: 353-359.6. Angeli D, Ferrell JE, Sontag ED. Detection of multistability, bifurcations, and hysteresis in a largeclass of biological positive-feedback systems. Proc Nat Acad Sci USA, 2004; 101: 1822-1827.7. Pomerening JR, Sontag ED, Ferrell Jr JE. Building a cell cycle oscillator: hysteresis and bistabilityin the activation of Cdc2. Nat Cell Biol, 2003; 5: 346-351.8. Kramer BP, Fussenegger M. Hysteresis in a synthetic mammalian gene network. Proc Nat Acad SciUSA, 2005; 102: 9517-9522.9. Shimizu Y, Tsuru S, Ito Y, Ying B-W, Yomob T. Stochastic switching induced adaptation in astarved
Escherichia coli population. PLoS One, 2011; 6: e23953.10. Gardner TS, Cantor CR, Collins JJ. Construction of a genetic toggle switch in
Escherichia coli.Nature, 2000; 403, 339-342.11. Veening J-W, Stewart EJ, Berngruber TW, Taddei F, Kuipers OP, Hamoen LW. Bet-hedging andepigenetic inheritance in bacterial cell development. Proc Natl Acad Sci USA, 2008; 105: 4393-4398.12. Ferrell Jr JE. Bistability, bifurcations, and Waddingtons epigenetic landscape. Curr Biol, 2012; 22:R458-R466.13. Vellela M, Qian H. Stochastic dynamics and nonequilibrium thermodynamics of a bistable chemicalsystem: the Schl¨ogl model revisited. J Royal Soc Interface, 2009; 6: 925-940.14. Schl¨ogl F. Chemical reaction models for nonequilibrium phase transitions. Z Physik, 1972; 253:147-161.15. Ge H, Qian H. Thermodynamic limit of a nonequilibrium steady state: Maxwell-type constructionfor a bistable biochemical system. Phys Rev Lett, 2009; 103: 148103.16. Zuk PJ, Kocha´nczyk M, Jaruszewicz J, Bednorz W, Lipniacki T. Dynamics of a stochastic spatiallyextended system predicted by comparing deterministic and stochastic attractors of the correspondingbirth-death process. Phys Biol, 2012; 9: 055002.417. T˘anase-Nicola S, Lubensky DK. Exchange of stability as a function of system size in a nonequilib-rium system. Phys Rev E, 2012; 86: 040103.18. Kurtz TG. Limit theorems for sequences of jump Markov processes approximating ordinary differ-ential equations. J Appl Prob, 1971; 8:344-356.19. Kurtz TG. The relationship between stochastic and deterministic models for chemical reactions. JChem Phys, 1972; 57: 2976-2978.20. Luo J-L, van der Broeck C, Nicolis G. Stability criteria and fluctuations around nonequilibriumstates. Z Phys B, 1984; 56: 165-170.21. Vellela M, Qian H. A quasistationary analysis of a stochastic chemical reaction: Keizer’s paradox.Bull Math Biol, 2007; 69: 1727.22. Hanggi P, Grabert H, Talkner P, Thomas H. Bistable systems: master equation versus Fokker-Planck modelling. Phys Rev A, 1984; 29: 371-378.23. Dykman MI, Mori E, Ross J, Hunt PM. Large fluctuations and optimal paths in chemical kinetics.J Chem Phys, 1994; 100, 5735-5749.24. Mehta P, Mukhopadhyay R, Wingreen NS. Exponential sensitivity of noise-driven switching ingenetic networks. Phys Biol, 2008; 5: 026005.25. Vardi N, Levy S, Assaf M, Carmi M, Barkai N. Budding yeast escapes commitment to the phosphatestarvation program using gene expression noise. Curr Biol, 2013; 23: 2051-2057.26. Schr¨odinger E. What is Life? Mind and matter. Cambridge: Cambridge University Press; 1944.27. von Bertalanffy L. The theory of open systems in physics and biology. Science, 1950; 111: 23-29.28. Qian H, Reluga TC. Nonequilibrium thermodynamics and nonlinear kinetics in a cellular signalingswitch. Phys Rev Lett, 2005; 94: 028101.29. Hopfield JJ. Kinetic proofreading: a new mechanism for reducing errors in biosynthetic processesrequiring high specificity. Proc Natl Acad Sci USA, 1974; 71: 4135-4139.30. Gaspard P. Fluctuation theorem for nonequilibrium reactions. J Chem Phys, 2004; 120: 8898-8905.31. Glansdorff PG, Prigogine I. Thermodynamic theory of structure, stability and fluctuations. London:Wiley-Interscience; 1971.32. Sawada Y. A thermodynamic variational principle in nonlinear non-equilibrium phenomena. ProgTheor Phys, 1981; 66: 68-76.33. Dewar RC. Maximum entropy production and the fluctuation theorem. J Phys A Math Gen, 2005;38: L371-L381.34. Andresen B, Zimmermann EC, Ross J. Objections to a proposal on the rate of entropy productionin systems far from equilibrium. J Chem Phys, 1984; 81: 4676-4678.35. Nicolis G, Lefever R. Comment on the kinetic potential and the Maxwell construction in non-equilibrium chemical phase transitions. Phys Lett A, 1977; 62: 469-471.36. Gillespie DT. Exact stochastic simulation of coupled chemical reactions. J Phys Chem, 1977; 81:2340-2361.537. Ebeling W, Schimansky-Geier L. Stochastic dynamics of a bistable reaction system. Physica, 1979.98A: 587-600.38. Andrews SS, Addy NJ, Brent R, Arkin AP. Detailed simulations of cell biology with Smoldyn 2.1.PLoS Comp Biol, 2010; 6: e1000705.39. Kalwarczyk T, Tabaka M, Holyst R. Biologistics - diffusion coefficients for complete proteome of
Escherichia coli . Bioinformatics, 2012; 28: 2971-2978.40. Marshall WF, Young KD, Swaffer M, Wood E, Nurse P, Kimura A, et al.
What determines cellsize? BMC Biol, 2012; 10: 101.41. Kaizu K, de Ronde W, Paijmans J, Takahashi K, Tostevin F, ten Wolde PR. The Berg-Purcelllimit revisited. Biophys J, 2014; 106: 976-985.42. Di Talia S, Skotheim JM, Bean JM, Siggia ED, Cross FR. The effects of molecular noise and sizecontrol on variability in the budding yeast cell cycle. Nature, 2007; 448: 947-951.43. Choi PJ, Cai L, Frieda K, Xie XS. A stochastic singlemolecule event triggers phenotype switchingof a bacterial cell. Science, 2008; 322: 442-446.44. Botello E, Nordstrom K. Effects of chromosome underreplication on cell division in
Escherichiacoli . J Bacteriol, 1998; 180: 6364-6374.45. Zaritsky A, Wang P, Vischer NOE. Instructive simulation of the bacterial cell division cycle. Mi-crobiol, 2011; 157: 1876-1885.46. Choi PJ, Xie XS, Shakhnovich EI. Stochastic switching in gene networks can occur by a single-molecule event or many molecular steps. J Mol Biol, 2010; 396: 230-244.47. Laughlin RB, Pines D, Schmalian J, Stojkovic BP, Wolynes P. The middle way. Proc Natl AcadSci USA, 2010; 97: 32-37.48. Barik D, Baumann WT, Paul MR, Novak B, Tyson JJ. A model of yeast cell-cycle regulation basedon multisite phosphorylation. Mol Syst Biol, 2010; 6: 405.49. Takahashi K, Tanase-Nicola S, ten Wolde PR. Spatiotemporal correlations can drastically changethe response of a MAPK pathway. Proc Natl Acad Sci USA, 2010; 107: 2473-2478.50. Gregor T, Bialek W, de Ruyter van Steveninck RR, Tank DW, Wieschaus EF. Diffusion and scalingduring early embryonic pattern formation. Proc Natl Acad Sci USA, 2005; 102: 18403-18407.51. Nenninger A, Mastroianni G, Mullineaux CW. Size dependence of protein diffusion in the cytoplasmof
Escherichia coli . J Bacteriol, 2010; 192: 4535-4540.52. Parry BR, Surovtsev IV, Cabeen MT, O’Hern CS, Dufresne ER, Jacobs-Wagner C. The bacterialcytoplasm has glass-like properties and is fluidized by metabolic activity. Cell, 2014; 156: 183-194.53. Hermsen R, Erickson DW, Hwa T. Speed, sensitivity, and bistability in auto-activating signalingcircuits. PLoS Comput Biol, 2011; 7: e1002265.54. Boulineau S, Tostevin F, Kiviet DJ, ten Wolde PR, Nghe P, Tans SJ. Single-cell dynamics revealssustained growth during diauxic shifts. PLoS One, 2013; 8: e61686.55. Zong C, So LH, Sepulveda LA, Skinner SO, Golding I. Lysogen stability is determined by thefrequency of activity bursts from the fate-determining gene. Mol Syst Biol, 2010; 6: 440.656. Acar M, Mettetal JT, van Oudenaarden A. Stochastic switching as a survival strategy in fluctuatingenvironments. Nat Genet, 2008; 40: 471-475.57. Hoffmann M, Schwarz US. Oscillations of Min-proteins in micropatterned environments: a three-dimensional particle based stochastic simulation approach. Soft Matter, 2014; 10, 2388-2396.58. Hopfield JJ. Neurons with graded response have collective computational properties like those oftwo state neurons. Proc Natl Acad Sci USA, 1984; 81: 3088-3092.59. Novick A, Weiner M. Enzyme induction as an all-ornone phenomenon. Proc Natl Acad Sci USA,1957; 43: 553-566.60. Singh A, Weinberger LS. Stochastic gene expression as a molecular switch for viral latency. CurrOpin Microbiol, 2009; 12: 460-466.61. Granot Z, Swisa A, Magenheim J, Stolovich-Rain M, Fujimoto W, E Manduchi E, et al.
LKB1regulates pancreatic beta cell size, polarity, and function. Cell Metabol, 2009; 10: 296-308.62. Andrews SS, Bray D. Stochastic simulation of chemical reactions with spatial resolution and singlemolecule detail. Phys Biol, 2004; 1: 137-151.63. Friedman N, Cai L, Xie XS. Linking stochastic dynamics to population distribution: an analyticalframework of gene expression. Phys Rev Lett, 2006; 97: 168302.64. Smith MB, Li H, Shen T, Huang X, Yusuf E, Vavylonis D. Segmentation and tracking of cytoskeletalfilaments using open active contours. Cytoskeleton, 2010; 67, 693-705.7
Figure 2. Properties of macroscopic bistable system. (A) Bifurcation diagram x ( B ) with thelow stable steady state in blue, the unstable steady state (saddle point) in red, and high stable steadystate in black for standard parameters defined in Materials and Methods. Black arrow indicates bistableregime. (B) Corresponding entropy production rate as defined in Eq. 2. (C) Phase diagram showingmonostable states (only low or high state) and bistable regions in β - γ plane with combinationparameters β = k − k +2 / ( k − B ) and γ = k +1 Ak / ( k − B ). Two phase diagrams correspond to v = ( k − B/k +2 ) V given by 37 (exemplified by combination V = 10 , B = 3 . ∞ (red lines). The latter corresponds to the macroscopic mean-field model. SPindicates point ( β, γ ) = (0 . , .
14) corresponding to standard parameters with B = 3 . Figure 3. Well-mixed bistable system. (A) Schematic of well-mixed system with volume V (diffusion constant D is infinitely large). (B) Exemplar time trace for x = X/V from Gillespie algorithmfor standard parameters with V = 10 and B = 4 .
0. (C) Exact probability distribution p ( x ) at steadystate from master Eq. 3 for V = 10 (dark symbols) and 30 (light symbols) with B = 4 .
0. (D) Values of p ( x ) evaluated at three steady states for different values of B . (E) Transition rates from a modifiedFokker-Planck approximation valid for large V (first-mean passage time; see S1 Text for details). Redarrows indicate exchange of stability. (F) Maxwell-like construction (MC), indicating coexistencebetween two phases (low and high states) at B ∼ .
7, defined by equal transition rates in (E). At thiscritical value of B a first-order phase transition occurs (see Text S1 for an analytical derivation based onsimpler potential). (G) Relative strength of fluctuations (standard deviation over mean) as a function of B for V = 30 (solid line), 50 (dashed line), and 100 (dotted line). (Inset) Unnormalized variances.9 Figure 4. Bistable system with diffusion. (A) Schematic of diffusing molecules in volume V . (B)Snapshot of cubic reaction volume for generalized Schl¨ogl model as simulated with Smoldyn software [38]. Shown are monomers X in red and dimers X in green. Clustering is illustrated by reddashed outline. (C) Chemical reactions of generalized Schl¨ogl model. (D) Time trace (left) andhistogram (right) of x = X/V from simulation for D = 3 (for X ) and 1 ( X ), V = 10, k +3 = k − = 1,and B = 3 .
7. (E) and (F) Effects of reduced (times 0 .
1) diffusion constants (E) and increased (times 2)volume (F). In (E) B = 3 . Smoldyn . Shown are monomers in red and dimers in greenwith illustration of clustering by red dashed outline. (I) Histogram of monomer concentration x fromsimulation for V = 2 .
14 and V DNA = 1 . D = 30 ( X ) and 10 ( X ), k +1 = k +2 = 50 and B = 50.0 Figure 5. Fragility of bistability. (A) Histograms of monomer concentration x as a function ofcontrol parameter B (from 3 . . .
1) with other parameters chosen as in Fig. 4D. (B)Radial pair-correlation function for D = 3 ( X ) and 1 ( X ) (dashed black line) and D = 0 . X ) and 0 . X ) (solid black line) compared with random distribution (red line; see Material and Methods fordetails). (C) Range of B values with visible bimodal distribution from master equation (well-mixed,black line) and (inhomogeneous) Smoldyn simulations for D = 3 ( X ) and 1 ( X ) (blue line) and D = 30( X ) and 10 ( X ) (red line) in units of µm /s , the latter being typical protein diffusion constants in thecytoplasm. System was classified bistable when 10,000-long simulations (see Materials and Methods)started in low and high states showed at least one reversible switch. Note that parameters are convertedto physical units here (see Materials and Methods for details). Hence, a 10,000-long-simulationcorresponds to a duration of 2.78h, which is a very conservative estimate of cell-division times inbacteria and yeast. Shaded areas indicate bacterial (dark) and eukaryotic nuclear (light) volumes forcomparison.1 Figure 6. Diffusion can be included by renormalization of second-order rate constants. (A)Histograms of monomer frequency for different B values from simulations with Smoldyn software [38](first row) and Gillespie simulations of the generalized Schl¨ogl model with (second row) and without(third row) renormalized rate constants of second-order reactions. Standard parameters were used withvolume V = 10. The Kullback-Leiber divergence ( D KL ) shows the closer correspondence of therenormalized reactions than the normal reactions with Smoldyn (forth row). For details onrenormalization and calculation of D KL , see Materials and Methods. (B) Corresponding macroscopicbifurcation diagram of deterministic ordinary-differential equation model using renormalized rateconstants k ± to illustrate effect of diffusion. This shows that diffusion delays entry into bistable regimefor increasing B .2 Figure 7. Onset of traveling waves in spatially extended system. (A) Snapshot of elongatedreaction volume for generalized Schl¨ogl model as simulated with
Smoldyn software [38]. Shown aremonomers X in red and dimers X in green. (B) Kymographs of monomer numbers along major axis ofsimulation box (distance) as a function of simulation time. For this purpose box was divided into 20equal sized bins. Parameter values: Standard parameters were chosen with volume of simulation box V = 10 · . · . B values as indicate in subpanels of (B), and other parameters as in Fig. 4D.Steepness of white dashed lines illustrates magnitude of wave velocity.3 Figure 8. Switching may be triggered by cell-volume changes. (A) Snapshots from time-lapsefluorescence microscopy: (left) lacY-gfp of
E. coli in yellow [43], (middle) PcomK-cfp of
B. subtilis inpurple, and (right) PcomG-cfp of
B. subtilis in red [2] with time in units of cell-cycle time T c . (B) Totalfluorescence intensities inside cell contours normalized to the maximal observed total intensity of a cell(see Materials and Methods for details) with color-coding same as in panel (A). Two yellow daughtercells are shown by solid and dashed lines. Note also the appearance of multiple red and purple daughtercells right after cell division in competence. (Inset) Normalized cell lengths over time in units ofmaximal cell length L max . S6 Fig. shows same for intensity density, i.e.i.e.