Multistable decision switches for flexible control of epigenetic differentiation.

It is now recognized that molecular circuits with positive feedback can induce two different gene expression states (bistability) under the very same cellular conditions. Whether, and how, cells make use of the coexistence of a larger number of stable states (multistability) is however largely unknown. Here, we first examine how autoregulation, a common attribute of genetic master regulators, facilitates multistability in two-component circuits. A systematic exploration of these modules' parameter space reveals two classes of molecular switches, involving transitions in bistable (progression switches) or multistable (decision switches) regimes. We demonstrate the potential of decision switches for multifaceted stimulus processing, including strength, duration, and flexible discrimination. These tasks enhance response specificity, help to store short-term memories of recent signaling events, stabilize transient gene expression, and enable stochastic fate commitment. The relevance of these circuits is further supported by biological data, because we find them in numerous developmental scenarios. Indeed, many of the presented information-processing features of decision switches could ultimately demonstrate a more flexible control of epigenetic differentiation.


Effective models; motivation
To describe the dynamics of many of the control circuits listed in Table 1 (main text), we consider two interacting regulators which are able to affect each other's production, while positively autoregulating their own expression. We further suppose these regulators to be transcription factors acting as homodimers, as is the case in bHLH proteins (e.g. Ac, Sc, Pax6, MyoG in Tables 1 and S1), leucine zipper factors or some types of homeodomain proteins, such as the POU factor Oct4 [1] or the caudal related protein Cdx2 [2]. Although many of the feedback interactions found in cell differentiation contexts are mediated by additional transcription factors, signals, etc., we seek to represent the basic dynamical features of these switches with a plausible mathematical model and a minimal set of parameters of biological significance. For instance, indirect regulation with the aid of an intermediate species can be simplified, under the assumption that this intermediate is in equilibrium, to the same mathematical equations discussed here. Thus, we start by writing the biochemical reactions for production and degradation of four molecular species: protein and mRNA of each of the two components (X, Y ) of the control module. These reactions can be separated in fast and slow.

Fast equilibrium reactions
We assume that DNA-binding reactions and dimerization are fast. Furthermore, they are lumped together in a single reaction term under the equilibrium assumption. Since each component is regulated by its own product as well as by the partner's, we consider in principle different promoter binding sites for X or Y molecules. These sites can be independent, interacting or completely overlapping as is the case in both prokaryotic [3] and eukaryotic transcription regulation [4], giving rise to a potential combinatorial logic of gene expression [5,6].
Here the K's represent products of equilibrium constants, being this the product of binding (protein-DNA) and dimerization constants, in Eqs. (1), or joint binding constants, K i xy (i = x, y) in Eqs. (2). In the case of two independent promoter sites, K i xy = K i x × K i y , while for overlapping promoters, only one species can be bound at a time and K i xy = 0.

Slow equations
Transcription, translation and degradation of both components are considered as slow reactions. Explicitly, and analogous equations for the Y species. Here β i (i = x, y) are the corresponding basal transcription rates, ρ i and ν i the auto and crossregulation strengths respectively, and µ i the joint regulatory strength. For independent promoter sites µ i = ρ i × ν i . With the constraint that the total number of promoters (copy number) is fixed, P T i = P i + P i X + P i Y 2 + P i X 2 Y 2 , and using the fast equilibrium reactions [Eqs. (1)-(2)] we obtain the time evolution of the four molecular species considered here as These are differential equations for the total number of proteins per cell. To convert to cellular protein concentration one needs to divide by cell volume, [X] = X/V . The time evolution for each species concentration is now where [P T i ] is the total promoter concentration for each gene, and the deterministic rate constants have been appropriately rescaled by cell volume for the series of bimolecular reactions, i.e., k i j = K i j ×V 2 and k i ij = K i ij ×V 4 , with {i, j} = {x, y}. Cell growth and division can be taken into account in an approximate manner by assuming an exponential law for growth [8], The growth rate is given by k g = ln(α)/T ccyc , where T ccyc is the mean division time and α a scale parameter giving the maximum increase in cell volume (typically α = 2). Then cell division occurs when V (t = T ccyc ) = αV 0 . Assuming that the total number of promoters scales proportionally with the cell volume, and using for the time evolution of the molecular , the resultant equations are identical to Eqs. (5) except that the degradation rates for mRNA and protein are now given by δ mx,my = δ mx,my +k g , δ x,y = δ x,y +k g , which are effectively taken as a single parameter. For stable proteins, δ x,y k g and degradation is just due to dilution by cell growth and division. Protein degradation can be faster (e.g., by proteases or by sequestration or inactivation with different molecular species), which implies that δ x,y k g and the growth term can be neglected.

Scaled variables and quasi-steady state approximation
We want to reduce both the number of variables to make the system amenable for mathematical and phase plane analysis, and the number of parameters to the minimum set with biological significance. First, it is convenient to define dimensionless mRNA and protein concentrations by Eqs. (5) now become with σ i = k i j /k i i (ratio of cross over self-binding rates) and σ ij = k i ij /(k i i · k i j ) (ratio of joint over independent rates for simultaneous binding of both species).
We note that for completely overlapping promoter sites (XOR transcriptional gate) σ ij = 0, whereas for independent promoter sites σ ij = σ i × σ j , µ i = ρ i × ν i and the non-linear production rates for mRNA are the product of two Hill functions, one for each variable. Applying the standard quasi-steady state approximation (QSSA, fast mRNA dynamics compared to protein dynamics), dmx dt = dmy dt ∼ 0, one readily obtains Eqs. (1) in the main text with where b i = γ i /δ mi is the burst parameter or translational efficiency.
To gain insight on the validity of this approximation, let as assume similar production and degradation rates for both components, i.e., γ x = γ y ≡ γ, β x = β y ≡ β, δ x = δ y ≡ δ and δ mx = δ my ≡ δ m . Rescaling time by the protein degradation rate, τ = t · δ, Eqs. (7) can be written as where ∆ = δ m /δ is the relative messenger-protein degradation rate. The variation of mRNA at short time scales can be obtained by considering the protein concentration fixed. Then, from Eqs. (9) one sees that for constant protein values mRNA decays as m x,y (τ ) ∝ exp (−∆τ ) while for constant messenger concentration the decay in protein concentration is simply ∝ exp (−τ ). Therefore it is reasonable to expect that mRNA dynamics adapts much faster to protein changes when ∆ >> 1. This is illustrated in Fig. S1. For ∆ = 10, the QSSA is a fairly good approximation to the four-dimensional dynamics (compare black and red lines), while for ∆ = 1 the two-and four-dimensional trajectories (black and blue lines) separate, although they eventually converge to the same equilibrium state in the long run. Indeed, we checked that the multistability domains in the phenotypic maps shown in Figures 2 and S3 are the same when calculated with the four-variable model [11]. The effect of changing b is important when considering the stochastic dynamics [9]. In fact, experimental studies of stochastic gene expression in single cells demonstrate that translational bursting is an important source of stochasticity in prokarytic [10] and eukaryotic gene expression when coupled to transcriptional noise [12].

Stochastic simulations and parameter ranges
Provided that equilibrium and quasi-steady state approximations hold for the deterministic system, stochastic simulations using Gillespie's method [13] can be simplified in a similar vein [14,15]. Therefore we used the four-dimensional model Eqs. (5) as a starting point for the stochastic simulations. To have a full correspondence with the simplified two-dimensional model, Eq. (1) in main text, we need to specify the following parameters: (1) Regulatory strengths (ρ i , ν i , µ i ), or fold-change in promoter activity once the proper transcription factor is bound. These parameters span from 0 (total inhibition) to 100 for the case of strong activation. This range reflects orders of magnitude typically found in both prokayotic and eukaryotic gene regulation.
(2) Relative binding affinities both for individual(σ i ) and joint (σ i xy ) species. σ i is varied between (0.01,2). A value of σ i < 1 denotes higher binding affinity of a given promoter for its own protein (and then the threshold for autoregulation is smaller than for cross-interaction) while σ i > 1 corresponds to the opposite situation. σ i xy is set to 0 for most of the paper calculations (completely overlapping promoter sites) but it can range from this value to σ i × σ j (independent promoter sites).
(3) We need also to specify the scaling factors k x x and k y y (promoter binding constants for each species) in the deterministic case. Recall that these are the product of dimerization and protein-DNA association constants. Here we fix k x x = k y y = 10 −3 nM −2 (a range 10 −1 − 10 −3 is typically found in bacteria [16]). The stochastic rate constants are obtained from the deterministic ones as K x x y · σ xy and K y xy = K y y · K y x · σ yx , where V represents cell volume. In the stochastic simulations cell growth and division can be taken into account by considering V as an additional random variable obeying the equation dV /dt = k g V , until a maximum value is attained [17,18]. Here we do not take into account the contribution of cell growth to gene expression noise, and fix V = V 0 · Ω in stochastic simulations, where V 0 is a reference cell volume including Avogadro's number and Ω a scale factor giving the contribution of protein number or finite size effects to molecular noise [19]. For a given value of molecular concentrations, the greater the scaling factor Ω the larger the number of molecules, i.e., noise is reduced. In our simulations, we take V 0 = 10 9 l for simplicity (note that concentrations are in nM) which implies a cell volume of 1.7 µm 3 . We also take Ω = 10 in most of the simulations. The effect of finite size noise in the protein steady-state distributions is illustrated in Figure S2 (compare black and red lines).
(4) Relative molecular stabilities. From the adimensional analysis of the fourvariable model, Eqs. (9), we also fix the ratio of protein/mRNA degradation to a value ∆ = 10, which implies a ten times faster degradation rate for mRNAs than for protein. Then we suppose that proteins are stable and mainly decay by dilution due to cell growth. Assuming a protein half-life of one hour (δ x = δ y = 1 h −1 ) then messengers are degraded on average every six minutes, which are typical values of prokaryotes.
(5) Burst parameter or translational efficiency. We use b x = b y = 1 for most simulations (for a glimpse on the effect induced by translational efficiency on molecular noise, compare black and blue lines in Figure S2).
(6) Finally, we specify the average protein production rate a i = α i /δ i as for the two-dimensional model Eq. (1). Setting the copy-number P T x = P T y = 1, then the basal transcription rates follow from the given parameters as Note that a ten-fold increase in translational efficiency is coupled here to a ten-fold decrease in transcription and thus the average number of proteins remains the same, but the number of messenger molecules is reduced (mRNA noise). On the other hand, a ten-fold increase in volume produces also a ten-fold increase in transcription and abundance of proteins and mRNAs are increased in the same way (finite size noise). The different effects of these two factors in the width of protein distributions are shown in Figure S2.

Signals
We model signals as molecular species external to the circuits considered. As such, they follow an independent dynamics that we approximate as a birth-death process: To create a signal pulse, we fix the degradation rate at δ s = 1h −1 . Initially, the signal production rate is α s = 0. At a given time, the signal is produced at a steady rate α s = s eq · δ s so that it quickly reaches their maximum value s eq (see insets in Fig. 5A). The signal is terminated at a later time when the source disappears and α s is set to zero again. For the deterministic equations, s is the signal concentration and α S has units of nM.h −1 . In the stochastic simulations, the noise in the signal can be varied independently from the noise in the circuit components defining the stochastic production rate α s s = α s · V s where V s is a volume factor changing the number of signal molecules. Signals are broadly considered affecting the circuit components in two possible ways: fast, inducing post-translational modifications of the X or Y proteins. In this case, the equations for (X, Y ) are only modified by an extra linear degradation/production term originated from the reactions Slow: signals may be transcription factors external to the circuit, that may act repressing or activating the (X, Y ) components. Here we assume that independent promoter sites are available for the signal molecules (although they could also operate cooperatively with any of both species). Equations are now modified according to the reactions.
where τ x,y are the signal regulation strengths (τ x,y < 1 for inhibitory signals, τ x,y > 1 for excitatory ones).
In the main text, we showed examples of signal processing with fast or posttranslational events. Similar computational properties can be obtained with slow or transcriptional signals in the appropriate regimes. However, slow signals may change the response dynamics and thus the discrimination performance of some stimulus features. An example is the discrimination of biased stimulus strength in stochastic decision switches (see Figure 4.C in main text). In a stochastic decision switch, a symmetric expression state (a stable node) becomes unstable due to a subcritical pitchfork bifurcation (and then becomes a saddle point, with the unstable direction being perpendicular to the diagonal in the phase plane). Although this bifurcation can be induced by both fast and slow signalling events, fast signals usually promote faster dynamics along the unstable direction, and this results in a faster discrimination of biased signal amplitudes. This is illustrated in Figure S6, where a decision switch operating close to the pitchfork bifurcation discriminates fast ( Figure S6.A) and slow ( Figure S6.B) inhibitory signals. One sees that not only discrimination performance of fast signals is steeper, but also signal noise affects this performance in a greater extent.

Role of autoregulation in symmetric progression switches
From the biological scenarios documented in the literature (see Table 1) it seems that genetic architectures of mutual activation of key transcription factors with autoregulation operate as standard (symmetric) progression switches. One possible reason inferred from our analysis is that, for moderate to strong crossinteractions, autoregulation should be exceedingly high to generate the coexistence of asymmetric expression states. What could be the role played by autoregulation in these type of circuits? Here we have investigated two possible reasons, both conferring flexibility to switch performance. One possibility is that, for a fixed crossinteraction strength where the circuit is not able to operate as a switch (being in a monostable regime), the modulation of autoactivatory loops may enable the system to work as a switch. This is demonstrated in Figure S9. For a standard switch without autoregulation ( Figure S9.A) if we keep the crossinteraction-strength at ν = 10, a (degradation) signal is not able to change the initial (high,high) expression state. For the same cross-interaction strength, and adding autoactivation to both circuit components, in a proper range of autoregulation strength the circuit is able to work as a switch upon the same signal ( Figure S9.B). The second possibility is that autoregulation may change the combinatory of signals to which a mutual-activation switch responds. In Figure S10, we show the response of a mutual-activation architecture placed in a (high,high) expression state to independent degradation signals in the X and Y components. In the case that there is no autoregulation ( Figure S10.A) the switch responds as a fuzzy OR gate (a signal affecting only one component above a given threshold induces a change of expression state). However, by adding autoactivation we need the simultaneous action of signals on both components (AND gate) to switch to the (low,low) expression state.