Biological Signal Processing with a Genetic Toggle Switch

Complex gene regulation requires responses that depend not only on the current levels of input signals but also on signals received in the past. In digital electronics, logic circuits with this property are referred to as sequential logic, in contrast to the simpler combinatorial logic without such internal memory. In molecular biology, memory is implemented in various forms such as biochemical modification of proteins or multistable gene circuits, but the design of the regulatory interface, which processes the input signals and the memory content, is often not well understood. Here, we explore design constraints for such regulatory interfaces using coarse-grained nonlinear models and stochastic simulations of detailed biochemical reaction networks. We test different designs for biological analogs of the most versatile memory element in digital electronics, the JK-latch. Our analysis shows that simple protein-protein interactions and protein-DNA binding are sufficient, in principle, to implement genetic circuits with the capabilities of a JK-latch. However, it also exposes fundamental limitations to its reliability, due to the fact that biological signal processing is asynchronous, in contrast to most digital electronics systems that feature a central clock to orchestrate the timing of all operations. We describe a seemingly natural way to improve the reliability by invoking the master-slave concept from digital electronics design. This concept could be useful to interpret the design of natural regulatory circuits, and for the design of synthetic biological systems.

1 Derivation of a reduced model for the JK-latch with independent heterodimer binding sites The reactions describing the genetic JK-latch, presented in Table S1 and S2 can be classified into three different categories: protein-protein reactions (including protein degradation), protein-DNA reactions and gene expression (including transcription and degradation of mRNA and translation into proteins).
A full deterministic model can be derived by setting up a rate equation for the mean concentration of each biochemical species in the reactions. A first step towards a reduced model is to assume the occupation states of the promoters (that is the protein-DNA reactions) to be in equilibrium with the current transcription factor concentrations at all times. Since all TF's act as repressors, we assume that the only promoter state contributing to mRNA production is the unoccupied one. Therefore the rate equation for mRNA of, for instance, gene A is, where ν m A is the maximal transcription rate, λ m the mRNA degradation rate and P A (t) denotes the promoter activity function, that is, the probability of the promoter to be unoccupied and free to bind RNA polymerase at a certain point in time. Hence, assuming all protein-DNA reactions to be in equilibrium, P A (t) describes the equilibrium probability to find the promoter unoccupied as a function of current transcription factor concentration. This probability can also be calculated by thermodynamic models, corresponding to [1,2]. Their specific form for the genetic JK-latch is defined in the boxes of Figs. 2 and 3.
As a further simplification, we assume all protein-protein reactions to rapidly reach equilibrium with the current total concentration of reacting proteins. For instance, for the concentrations of all protein species containing A and K this leads to the following set of equations 1 where we can additionally set K = K tot − KA. After solving this system for the concentrations of homoand heterodimers all promoter activity functions in Figs. 2 and 3 of the main text can then be expressed in terms of total concentrations, i.e., P A (B 2 , KA) = P A (A tot , B tot , K tot , J tot ).
Altogether, we are left with two equations left to describe the total concentration of each gene product.
For instance, the equations for gene product A are As a last step, we set Eq. (6) to zero thereby assuming mRNA translation to be fast with respect to transcription (we will discuss this assumption in the next section). Solving the resulting equation for m A and substituting it in Eq. (5) leads to with an effective maximal expression rate ν A = ν m A ν p A /λ m . Thus, each gene product in our models can be described by a single effective equation.
2 Derivation of a reduced model for the genetic J-K latch with mutually exclusive heterodimer binding sites Making the same quasi-equilibrium assumptions for protein-protein and protein-DNA as in Section I, the reduced model of the genetic J-K latch with mutually exclusive binding sites for the heterodimers KA and JB reads where only the operator occupancy functions of the heterodimers now also depend on the respective other heterodimer: (9c) However, in order to introduce a delay required for a successful toggle operation, the unbinding kinetics of the overlapping operator sites needs to be slow (see main text), such that the quasi-equilibrium assumption that lead to Eqs. (9a) and (9b) is not strictly valid. To describe the slow dynamics of the overlapping operator complex we set up the master equation for the occupancy states of the overlapping operators. There are three states to be accounted for: (i) the binding site for KA is occupied, (ii) the binding site for JB is occupied and (iii) both binding sites are unoccupied. We denote the probabilities of these states as q KA , q JB and q 0 , respectively. The probabilities equivalent to Eqs. (9b) and (9c) to find the operator in an unoccupied state are then given by O KA = 1 − q KA and O JB = 1 − q JB . Taking k on to be the on rate of both operators and k off the respective off-rate, the master equations read: These equations are controlled by the external variables KA and JB . For given values of these variables, the system has a unique fixed point, which we denote by (q * KA ,q * JB ,q * 0 ). In the fixed point, the right hand side of Eq. (10) becomes zero, thus we can rewrite Eqs. (10a) and (10b) aṡ Furthermore, we want to investigate Eqs. (10) under a constant toggle signal. Therefore, it is reasonable to assume that both heterodimers KA or JB are abundant in the system at all times. Then, the effective on-rates k on KA and k on JB are much faster than the off-rate k off and thus the probability q 0 to find the overlapping operator sites unoccupied is approximately zero. Additionally, substituting q KA = 1 − O KA and q JB = 1 − O JB , the approximated equations can be written aṡ We can identify the fixed points O * KA and O * JB with the quasi-equilibrium operator occupancy functions defined in Eqs. (9b) and (9c), which are functions of the time dependent heterodimer concentration.
These equations can then be solved by variation of the constant: with "memory kernel" g k (τ ) = k off e −k off τ . Instead of accounting for the initial conditions in an own term, we formally assume that for all negative times the system is held fixed at the initial point and integrate infinitely long into the past. This auxiliary assumption also simplifies the initial condition problem for the delayed dynamical system: for a distributed delay, instead of a single initial value, strictly the entire history prior to t = 0 has to be defined. Since we are interested in the qualitative long term behavior of the solution, we assume this history to be a single point from which the system is released at t = 0. Putting this result back into Eq. (8) leads to the delayed differential equations with exponentially distributed delay introduced in Fig. 3 Model 2.
Additional delays in the reduced model. For the following stability analysis, we include the additional delay caused by the translation process. While this is done for the sake of completeness, this additional delay is about ten times smaller than the delay caused by the overlapping heterodimer operators and is therefore not considered in the discussion of the main text. The additional delay arises from the reaction chain between the beginning of mRNA transcription and completion of protein translation. In principle, these kinds of delay alone can already lead to oscillations [4,5]. Consider the rate equation for an arbitrary protein concentration A, with regulated transcription, linear degradation and an explicit equation for its mRNA denoted by m A : Here, P A (t) denotes an arbitrary promoter activity function of gene A. Solving Eq. (15) and treating the initial conditions as above leads to Putting this back into Eq. 14 leads to with an effective expression rate ν A = ν p ν m /λ m and a delay kernel g λ (τ ) = λ m e −λmτ . Thus, taking the dynamics of mRNA into account, leads to an additional distributed delay acting on the entire promoter activity function. In a reduced model the protein concentration therefore responds to a change in promoter activity on a timescale given by the mRNA degradation rate. This is generally valid for any gene controlled by regulated recruitment. There are, however, mechanisms of gene regulation, involving active degradation of mRNA [6,7]. In that case, the dynamics of mRNA have to be accounted for by its full dynamic description.

Linear stability analysis on a system with distributed delays
Here we discuss linear stability analysis for nonlinear systems with multiple cascaded exponentially distributed delays, which closely follows the book of [8]. The results presented here are generally applicable to any system of that kind. The model system considered here is the simple J-K latch with overlapping heterodimer operators in Eqs. (8) and (9). As discussed in the last sections, it contains delays on the promoter activity functions P A and P B and, within those, on the operator occupation functions O KA and O JB . The memory kernels g k (τ ) and g λ (τ ) of these delays are exponential distributions with rate parameters k off and λ m , respectively.
To perform a linear stability analysis on our model system, it first is important to note that while delays change the stability of a fixed point, they do not alter their position with respect to the undelayed system. This is easy to see: consider a fixed point FP = (A * , B * ) of the undelayed system and the where D g k abbreviates a distributed delay with memory kernel g. We assume that the delayed system has stayed in this fixed point for a very long time, such that the delayed function can be written as for any normalized kernel g k (τ ). Therefore, a fixed point of the undelayed system is a fixed point of the delayed system as well.
To perform a linear stability analysis on the model system Eq. (8), we must take into account that the dynamic variables A and B occur within different delay integrals and hence need to be treated separately in the linearization around a fixed point FP . Therefore, we categorize A and B by the delays acting on them. The model system iṡ which can more formally be written aṡ Here, variables with subscript g λ occur in functions, which are delayed only by mRNA degradation, whereas variables with subscript g λ,k are delayed by mRNA degradation and the off-rate of the overlapping operators. We now can linearize the system around a fixed point FP = (A * , B * ) and obtain linear where D g λ denotes a delay with kernel λ m e −λmτ and D g k a delay with kernel k off e −k off τ . Since the position of the fixed point does not change for delayed variables, all partial derivatives are evaluated in the same point FP . Taking into account that the functionals D are linear, the system can be written as Herein, the delayed variables take the form We make the usual ansatz for a linear system:Ã = c A e zt andB = c B e zt with z ∈ C to be determined.
Upon this ansatz, the delayed variables become where L(g λ (τ ); z) denotes the Laplace transform of the delay kernel. This holds for any normalized delay kernel. Since g λ (τ ) is an exponential distribition, its Laplace transform takes a particularly simple form: Variables with two delays can be evaluated in a similar way: Putting these results back into Eqs. 25 and 26 yields a linear equation for the coefficients c A and c B : Note that all partial derivatives have to be evaluated at the fixed point under consideration. Subtracting the left hand side from the diagonal elements of the right hand side leads to the familiar form a an eigenvalue problem. Here, in contrast to a common linear stability analysis, additional polynomial terms are incorporated, representing the delays. In order for the system to have a solution, its determinant has to be zero. Thereby we obtain a characteristic equation for the exponents z.
Evaluating the determinant, we obtain a polynomial of 6th order in z. In order for the considered fixed point to be stable, each solution of Eq. 33 must have a negative real part. The timescales of the system considered here are given by τ λ = 1/λ m and τ k = 1/k off , which are the mean values of the memory kernel. These characteristic times are usually referred to as the mean delay [8]. Written in terms of the mean delay τ , the Laplace transform of an exponential distribution becomes 1/(1 + τ z). Thus, if the delay is small, its contributing terms in the characteristic equation are close to one and therefore do not change the stability of a fixed point. In the following we will employ the derived method of linear stability analysis to investigate the genetic JK-latch with overlapping heterodimer operators. (iv) If necessary, characterize the system's qualitative behavior (e.g. oscillatory) by solving its equations numerically. Employing this scheme for different parameter values we identify two important bifurcations: the system undergoes a Hopf bifurcation [9,10] if the concentration of input proteins J tot and K tot is increased (see Fig. S2), then is oscillatory for a certain concentration range until it becomes stable again. Additionally, the system exhibits another Hopf bifurcation if the mean delay time is increased (i.e. the parameter k off is decreased), while a constant toggle signal is applied. This is illustrated in Fig. S3A: at a critical operator dwell time τ crit ≈ 32 min the system becomes oscillatory with a period that increases with the dwell time.
In the form we set up the reduced model for the genetic JK-latch with overlapping heterodimer operators, we couple the dynamics of protein concentrations with the intrinsically fluctuating switching of operator states. Therefore, although the deterministic model provides valuable insights into the working principles of the circuit, the quantitative period of the stochastic system is not correctly reproduced. This is shown in Fig. 3B, where it is apparent that the average of the stochastic period is much shorter than the period of the deterministic model. Also the stochastic system displays oscillations in a regime of short heterodimer dwell times, in which the deterministic system does not oscillate at all. In fact, it is known that a small number of reactant molecules together with negative feedback and time delay in gene expression can lead to delay-induced instabilities, such that a system turns oscillatory even when its deterministic counterpart is not [4].
If the delay (dwell time) is long, the strength of the toggle signal is chosen appropriately and the signal duration is well timed, then high probabilities of a correct toggle response can be achieved -as can be seen in Fig. S4A and B.

Genetic master-slave latch
We extended the reaction system of the JK-latch to a master-slave latch by adding additional regulation of the signal genes J and K -see Table S4. This was done by firstly including homodimerization of signal proteins. J 2 and K 2 bind two operators on the promoter of the respective other gene to repress it. Additionally, gene J is repressed by heterodimers KA by binding to a third operator in its promoter region, while gene K is repressed by JB in the same way. Notably, the restrictive promoter layout of the overlapping heterodimer operators on genes A and B with very slow unbinding rates is no longer necessary in the master-slave latch and is changed to independent operators with normal unbinding rates.
With this additional regulation, signal genes J and K become bistable when their maximal transcription rate is increased externally (by the toggle signal) and thereby form the master toggle switch.
Then, the master toggle switch settles into a state that is determined by the current concentrations of heterodimers and therefore by the current state of the slave latch. Dependent on the operator strengths for homodimers and for heterodimers two kinds of erroneous responses are possible upon a toggle signal: If homodimers bind to strongly with respect to heterodimers, then the master switch is prone to stochastic fluctuations at the onset of the toggle signal and biased towards that signal protein which is first abundant enough to form homodimers. If, on the other hand, heterodimer binding is too strong, random single bursts of the repressed signal gene can be sufficient to switch the state of the master latch. Fig. S4C shows the toggle probability as a function of those key parameters for a toggle signal duration of 60 min. It can be seen that the master-slave latch is robust in that sense, that it responds correctly for a broad range of those parameters. Table S1. All reactions involved in the genetic JK-latch with overlapping heterodimer binding sites. A graphical representation of the reaction network is shown in Fig. S1. Proteins and their dimers are denoted by capital letters; transcripts of a gene X are denoted by m X . A gene X is represented by its promoter P X , which can be occupied by transcriptions factors. Each occupation state of a promoter is represented by an own species of reactants for which an empty operator is indicated by · and an occupied operator by the name of the respective transcription factor. In this notation the different binding sites are separated by the symbol |. To reduce the number of occupation state combinations the operator complex for heterodimers is separated from the other promoter states and denoted by O. To make transcription nevertheless conditional on heterodimers, we include the respective species for an empty binding site as reactant and product into the corresponding transcription reaction. This ensures that transcription requires an empty heterodimer operator to proceed but does not change its concentration. All corresponding parameter values are listed in Table S3.

JK-latch with overlapping heterodimer operators Promoter and operator states
Additional reactions in J-K latch without overlapping heterodimer operators Table S2. With independent heterodimer operators, heterodimers can bind simultaneously to their binding sites. All other reactions are the same as in the JK-latch with overlapping heterodimer operators, listed in Table S1.

Parameter Value Description and References
Transcription ν m J , ν m K (0.05 − 5) min −1 strong inducible promoter [11]; induction can be achieved e.g. by upstream binding activators [12] or via small non-coding RNAs [13] Dimerization s on J , s on K 0.2 nM −1 min −1 assumed to be diffusion limited [17] s off J , s off  Table S5. Additional parameters used in the genetic master-slave latch. Instead of the symbol k for rates in the slave circuit, the symbol s is used for rates of all addtional processes (dimerization and protein-DNA binding of signal proteins) in the master circuit. Since a delay in heterodimer binding is no longer needed, the assumption of very slow unbinding kinetics of heterodimers from promoters of genes A and B has been released.  Fig. (A). The root is critical in that sense that it is the only solution of the characteristic equation, which becomes positive for certain values of J tot /K tot . Starting from a low concentration of input proteins the system initially has two stable (solid line) and one unstable (dotted line) fixed point. The upper stable fixed point represents the ON state, whereas the lower one represents the OFF state. As the concentration of input proteins is increased, the two stable fixed points loose their stability at a critical concentration c 1 , after which delay induced oscillations commence (orange lines). At input protein concentration c 2 the three fixed points collapse to one. That, however, does not alter the system's oscillatory behavior. At concentration c 3 the system's single fixed point becomes stable again. This is due to depletion of homodimers by increasingly forming heterodimers, which eventually abolishes the switch-like behavior of genes A and B. . At a critical dwell time τ crit ≈ 32 min the system undergoes a Hopf bifurcation and thereafter oscillates with a period that is approximately linear in τ . (B) Average period of the stochastic race condition as a function of the mean dwell time τ . In comparison to (A) the stochastic system oscillates even for small dwell times (although with very low amplitude) and has a shorter period than the deterministic model. The latter is due to the fact that the deterministic model couples the dynamics of protein concentrations with probabilities of operator switching, which has a distorting effect on the oscillatory dynamics.

J-K latch (exclusive heterodimer binding)
Master-Slave latch Figure S4. Toggle probability p toggle as a function of key parameters. Contours indicate equal probability to toggle successfully into the complementary state. Each data point (a grid of 30 by 30 per plot) is estimated by testing the final state of 5000 simulation runs of the respective full stochastic model. (A) Dependence of p toggle on the duration T and the strength of the toggle signal in the JK-latch. Here the strength of the toggle signal is tuned by a concerted variation of the transcription rates of genes J and K, i.e., ν J = ν K = ν. (B) Dependence of p toggle on the duration T of the toggle signal and the off-rate k off for unbinding the overlapping operator sites. As expected, p toggle increases as the delay (dwell time of heterodimers on the binding site) is increased. (C) Dependence of p toggle in the master-slave latch on the additional rates s off J2 , s off K2 for unbinding the homodimer operators and s off JB , s off KA for unbinding the heterodimer operators in the additional toggle switch (master latch). The master-slave latch is robust in that sense that p toggle is high for a broad range of parameters.

Figure S5
Dynamics of the deterministic model of the master-slave latch. (A) Phase diagram of the toggle switch as a function of the effective maximal transcription ratesν A andν B . The parameters are chosen such that the system is in the bistable regime in the absence of input signals (point O) and the circuit is set to the ON state initially. The curves indicate dynamic changes ofν A andν B , incurred by applying the toggle signal (simultaneous expression of both input genes J and K) for 200 min (green solid curve) and then releasing it (red curve). In contrast to the JK-latch, even if the toggle signal is applied a long time, the system enters the correct monostable regime (A low, B high), switches to the OFF state and returns to the bistable region without approaching the other monostable regime. In particular, the master-slave latch does not oscillate, even under a continuous toggle signal.