Non-monotonic auto-regulation in single gene circuits

We theoretically study the effects of non-monotonic response curves in genetic auto-regulation by exploring the possible dynamical behaviors for such systems. Our motivation is twofold: we aim at conceiving the simplest genetic circuits for synthetic biology and at understanding the natural auto-regulation of the LrpB protein of the Sulfolobus solfataricus archaeon which exhibits non-monotonicity. We analyzed three toy models, based on mass-action kinetics, with increasing complexity and sought for oscillations and (fast) bistable switching. We performed large parameter scans and sensitivity analyses, and quantified the quality of the oscillators and switches by computing relative volumes in parameter space reproducing the sought dynamical behavior. All single gene systems need finely tuned parameters in order to oscillate, but bistable switches are more robust against parameter changes. We expected non-monotonic switches to be faster than monotonic ones, however solutions combining both auto-activation and repression in the physiological range to obtain fast switches are scarce. Our analysis shows that the Ss-LrpB system can not provide a bistable switch and that robust oscillations are unlikely. Gillespie simulations suggest that the function of the natural Ss-LrpB system is sensing via a spiking behavior, which is in line with the fact that this protein has a metabolic regulatory function and binds to a ligand.

A The models: Ordinary differential equations and quasi steady states The systems were modeled through mass action kinetics. The models take into account protein-DNA binding and unbinding, dimerization, transcription, translation and dilution/degradation of mRNA and proteins.

A.1 MDS
The monomer-dimer system contains five ordinary differential equations (ODEs):   or dimer to site i (i ∈ 1, 2, 3) n.a. cooperativity factor between sites i and j (and k)

A.2 2DS
Analogously to the system of ODEs for the MDS, we can write down the system for the 2 dimer system which contains six ODEs A summary of all variables and parameters is given in Table A. The quasi steady states of the DNA complexes and dimer are: DNA 1, qss = DNA 0, qss K 1 d, DNA 2, qss = DNA 0, qss K 2 d and

B Physiological ranges of the parameters
Parameter Boundaries Parameter Boundaries C Oscillations

C.1 Bifurcation analysis
In order to check if a system oscillates, we first look if we can find a single steady state which is unstable (at least one of the eigenvalues of the Jacobian evaluated in this steady state has a positive real part). Unstable steady states can lead to solutions that oscillate, chaotic behavior or solutions that go to infinity. To be certain the system is oscillating, we do a time series and look for two characteristics of stable oscillations: 1. Stable amplitude: the difference between consecutive maxima in the time series must be less than 1%, the same holds for consecutive minima.
2. Stable period: the difference between to consecutive distances between local maxima must be small.
We want to find the oscillatory regions around a found oscillating solution for every parameter. This means all parameters are kept constant except for one, the parameter of which we determine the bifurcation points. Therefore we use method based on the bisection method for finding roots. To find the left boundary of the oscillating region, we first check whether the system with the considered parameter equal to the left boundary of the physiological range is oscillating. If this is the case, we have found the left boundary of the oscillating region, if not, we know that the bifurcation happens in the interval [l 1 , l 2 ] where we initialize l 1 with the parameter value at the left boundary of the physiological range and l 2 with the parameter value of the oscillating solution. Next we will update the boundaries of the interval around the bifurcation until the interval is considerably small and approximate the bifurcation point by the point that is logarithmically halfway this interval. To update the interval we take the point logarithmically halfway the interval l m = 10 log 10 (l 1 )+log 10 (l 2 ) 2 (1) and check whether it is oscillating. In the case that it is oscillating we update right bound of the interval, l 2 = l m , otherwise the left bound l 1 = l m . The right boundary of the oscillating region is found in a similar way. Finding these two outer solutions does not guarantee that all solutions in the range are oscillating, since other bifurcations might happen in between. Therefore we computed time series for five logarithmically equally spaced samples to see if they were oscillating (large computation time limits the number of solutions that can be tested).

C.2 Oscillatory ranges distributions
In the main paper, the mean values and mean ranges are given for the different parameters. Here different bifurcation ranges are given for different solutions (Figs A, B, C and D). We can see that some parameters can have different values in the physiological range but still be finely tuned, for example β or φ 0 for the MDS. For other parameters such as f i of the 2DS the distribution is bimodal. These characteristics are not represented in the main paper, but are unimportant for a global overview of the oscillating regions. Because there is much variation in the ranges, there is also much variation in the volumes (Fig E).

D Bistability
In order to check if a system is bistable, we look if we can find three steady states of which one is unstable (at least one of the eigenvalues of the Jacobian evaluated in this steady state has a positive real part) and the other two stable (all eigenvalues of the Jacobian evaluated in this steady state have a negative real part).

D.1 Quasi steady state approximation
In the assumption of fast binding and unbinding of the dimers to the DNA and fast dimer association and dissociation, we can equal the time derivatives of the DNA complexes and the dimer concentration d to zero and use their quasi steady state approximation in the equations for mRNA and monomer concentration m (given in Section 1): with f (m, d) = i f i DNA qss, i /DNA tot the transcription function. This function depends on the system, using the steady state expression of Section A we obtain with

D.2 Induction time of bistable systems
The induction time is a representation of the time it takes a bistable system to attain the stable high steady state when starting from the unstable intermediate steady state. Deterministically, systems  than 90% of the high steady state, The induction time ∆t is thus defined by the time it takes the system to attain a dimer number higher than 0.9d H starting from initial condition s i .
Simulations show that different choices for K i , f i , co bi and co ui -when coefficients A, B, C, D, E and F are fixed-do not influence the induction time. Other free parameters get the values dictated by Table C. Parameter Value We compare the induction time calculated by an explicit simulation as explained in this section with the approximated induction time∆t defined in the main text of the paper (Fig F). As mentioned earlier in this section, the initial condition and final state need to be defined. These will affect the simulated induction time. Depending on the exact shape of the response curve and the parameters, the approximated time is thus an over-or underestimation of the simulated induction time. This explains why there is no one-on-one match of the approximated time and the simulated time.

D.3 Bistability in parameter space
In Fig G in the regions within parameter space providing bistable switches for the different systems are represented. For 2DS and 3DS, a large proportion of parameter space leads to bistability while the MDS needs to be finely tuned to provide this dynamical property. In this figure, we used parameters A − F . In order to study the parameters listed in Table A, we selected 100 random solutions in the bistable region and performed a bifurcation analysis. This bifurcation analysis is similar to the one for oscillations, but the boundaries are checked for bistability instead of oscillations. The mean parameter value and mean width of the bistable range is given in

E SsLrpB compatibility
In order to be compatible with the natural SsLrpB system, a 3DS needs to have the same values as measured experimentally (Table D) and moreover the response curve needs to be of non-monotonic type 1 with a maximum higher than 2 and a minimum lower than 0.5 [5]. The radius of Sulfolobus solfataricus is around 1 µm, thus the volume can be estimated at 4 fL. A concentration of 1 µm in one cell corresponds therefore to 2400 molecules. Since the non-monotonicity is experimentally shown between 0 and 100 nm [5], we only considered systems up to 1000 molecules (equivalent to 300 nm), i.e. the response curve needed to meet the above constraints on non-monotonicity in the interval from 0 to 1000 dimers, only considering the integer numbers.
The measured parameters of Table D fix parameters, A, D, E and F , which are defined in the main paper, such that only 2 free parameters remain for the bistable switch search : B and C. The former is determined by f 12 , f 23 and f 13 and the latter by f 123 . For each f 13 − f 123 combination, we scanned over the only free parameters f 12 and f 23 , for which we get different response curves. These are shown in grey in Fig I. Bistable curves meeting the conditions as explained in the main paper (Eq. 1 of main paper) have a dashed red line and curves meeting the Ss-LrpB criteria have a dashed blue line. For every bistable (red) response curve, different degradation rates are tested (green lines) and the one with the smallest approximated induction time is chosen. A simulation is done with this solution to calculate the real induction time, which is shown by the yellow-green coloring in Fig 9 of the main paper.