## Abstract

The assumption of the fast binding of transcription factors (TFs) to promoters is a typical point in studies of synthetic genetic circuits functioning in bacteria. Although the assumption is effective for simplifying the models, it becomes questionable in the light of *in vivo* measurements of the times TF spends searching for its cognate DNA sites. We investigated the dynamics of the full idealized model of the paradigmatic genetic oscillator, the repressilator, using deterministic mathematical modelling and stochastic simulations. We found (using experimentally approved parameter values) that decreases in the TF binding rate changes the type of transition between steady state and oscillation. As a result, this gives rise to the hysteresis region in the parameter space, where both the steady state and the oscillation coexist. We further show that the hysteresis is persistent over a considerable range of the parameter values, but the presence of the oscillations is limited by the low rate of TF dimer degradation. Finally, the stochastic simulation of the model confirms the hysteresis with switching between the two attractors, resulting in highly skewed period distributions. Moreover, intrinsic noise stipulates trains of large-amplitude modulations around the stable steady state outside the hysteresis region, which makes the period distributions bimodal.

## 1. Introduction

Living cells are multi-stable and multi-rhythmic dynamical systems. Coexisting multiple stable dynamic behaviours underlie biological diversity and cell differentiation. However, the knowledge of how cells sustain their diverse dynamic behaviours, especially in the presence of molecular noise, is still limited. During the 1970–1980s, there were many studies of cell variability in biochemistry and in cell surface processes. However, the focus has shifted today as the variability is mainly addressed in genetics with synthetic gene circuits being the main source of models to study (e.g. [1]).

Here, we use one of the prominent examples of genetic oscillators [2], known as the repressilator [3], which is a three-element genetic circular network, in which each of the elements unidirectionally represses the adjacent one (figure 1).

Thus, gene *G*_{1} is a source of mRNA *M*_{1} (transcription) that codes for the monomer protein *P*_{1} produced in the translation process. The proteins *P*_{1} participate in reversible homo-dimerization producing dimers *D*_{1}, which are repressing transcription factors (TFs) for the next gene *G*_{2}. Similarly, gene *G*_{2} produces (via transcription, translation and dimerization) dimer molecules *D*_{2}, repressing transcription from *G*_{3}. Finally, TF molecules *D*_{3}, produced from gene *G*_{3}, repress production from *G*_{1}, completing the cycle (figure 1).

The idea of the repressilator is old [4], but it has been realized experimentally in live bacterial cells (the three-gene network was incorporated into a small circular DNA molecule called plasmid) and modelled by means of differential equations and using stochastic simulations only in 2000 [3].

For many years, it was well accepted to use the quasi-steady-state approximation (QSSA) for the multimerization of TFs and regulatory reactions on the promoter sites as a standard assumption for the deterministic modelling in genetics [5,6]. The approximation, sometimes called ‘adiabatic’ [7], results in important model simplification, giving the reduced dimension of the model and the appearance of the Hill function describing transcription regulation. Using this approach, the repressilator was formulated as a 6-dim system of ordinary differential equations (ODE), which exhibits sustained oscillations emerging via the super-critical Andronov–Hopf bifurcation (superAHB) [3], which has been shown analytically in the limits of QSSA [8,9]. This bifurcation (a qualitative change in the dynamical behaviour as parameters of the system are varied) indicates the presence of the stable steady state, which loses stability through the bifurcation and gives rise to small-amplitude oscillation (limit cycle attractor). A repressilator composed of a larger number of genes (both odd and even) was also shown to have oscillation emerging via superAHB [10].

Discrete stochastic approaches are widely used in the dynamical studies of gene regulatory networks [11–13]. However, the stochastic approach implies an additional statistical analysis and time-consuming simulations, thus, urging many authors to use reduced models formulated with QSSA. Using the Hill and Michaelis–Menten functions for the compact description of enzymatic activities is a well-accepted technique that has been successfully applied to many genetic circuits, demonstrating results for the average steady-state characteristics in accordance with the real systems. Nevertheless, when fluctuations are allowed, the models with these nonlinear functions as transcription rates fail to reproduce the corresponding results [14].

Finally, in the standard reduction of deterministic ODEs according to the QSSA, the change of TF concentrations is considered as a slow process. However, protein concentrations depend not only on the slow reactions (e.g. translation and degradation), but also on the two fast reactions (dimerization and binding/unbinding from promoters). It has been shown in [15] that there is a correct procedure for the reduction, where the protein concentrations are considered as a mixture of both slow and fast variables, resulting in the perfect quantitative agreement between the full and reduced models. On the other hand, the standard reduction of the repressilator leads to great quantitative distortions in values of the period and amplitude of oscillations as well as to a significant shift of the dynamical regime boundaries in the parameter space [15].

In this work, we study the dynamics of the repressilator without using QSSA bearing in mind the recent experimental results about the complex way in which TF (LacI dimer) searches for its cognate site on DNA [16]. The *in vivo* searching times are estimated as several minutes, which is comparable with the mRNA and TF degradation times, which are of the order of 5 and 15 min, respectively. Recently, the specific-site searching time for all of the 180 known *Escherichia coli* TFs were analysed and times between 160 s and 1550 s were determined [17].

We explore the minimal but detailed model of the repressilator with seven elementary reactions explicitly described as a 12-dim ODE system to discover deterministic dynamical regimes. Stochastic simulations are used to examine the role of internal noise due to the discreteness of the molecular populations. We find that the non-adiabatic character of the system, which contains a more or less slow subsystem of TF–DNA interaction, opens the possibility for oscillation emergence via the subcritical AHB (subAHB), which may give rise to the coexistence of the steady state and oscillation, i.e. both dynamical manifestations are present for the same parameter values. A non-adiabatic repressilator was studied previously in [18] mainly by means of stochastic simulations, but the dynamical properties of the system were not consistently explored.

Parametric analysis reveals the subAHB in a broad range of control parameters' values, especially, if the degradation of TF dimers is taken into account. A repressilator formulated in this way becomes very sensitive to even small changes in the dimer degradation rate, which, to our best knowledge, has not been reported previously in the literature.

The subAHB reveals the hysteretic properties of the model, that is, two stable dynamic regimes, steady state and limit cycle, coexist in the parameter space opening the essential possibility for switching between the two, e.g. due to stochastic effects. The hysteresis is further confirmed with the stochastic simulation, where we have shown the existence of two distinct dynamical behaviours for a single parameter set. Although the ultimate coexistence of the two dynamical attractors has been shown, the means for switching and control of the behaviours requires further investigation.

Furthermore, it has been shown that the slow TF–DNA interaction results in slow and large-amplitude stochastic modulations of the stable steady state, which leads to bimodal character of the return time (period) distributions. However, in the hysteretic region, the lifetime of the limit cycle is significantly longer than that of the stable steady state, which is manifested by the return time distribution biased towards longer periods.

## 2. Model

We use an idealized model of the repressilator circuit [15,19] with an explicit account for the protein dimerization process and reactions on the promoter site (figure 1). For all three genes, we assume promoters, mRNA, proteins and dimers to be identical in all respects, including the same parameter values for similar reactions.

The corresponding ODE system is [15]
2.1
where the concentration notations are the same as in figure 1 and *i* = 1, 2, 3, *k* = 3, 1, 2 and *j* = 2, 3, 1.

We opted to use parameters that were used in recent publications on closely related subjects [13,15,20] (and references therein). The default parameter values are shown in table 1.

Transcription rate *r*_{m}, monomer and dimer degradation rates *d*_{p} and *δ*, respectively, undergo significant variation throughout this study, as the exact estimation of the parameters is difficult in live cells.

Finally, we vary *k*_{r} and *k*_{ur} to switch between adiabatic and non-adiabatic conditions. The actual values for these parameters are not fixed to stress that the effects observed in this study do not result from the specific choice of parameter values.

## 3. Material and methods

We study the ODE system (2.1) by means of bifurcation analysis [23]. Furthermore, we use the stochastic simulation algorithm [11] to introduce intrinsic noise to the system. We smooth the stochastic time series using the moving average algorithm. To calculate return times distributions (referred to as period distributions in this work) from the averaged time series, we define the Poincaré section level at the middle between maximum and minimum points (see electronic supplementary material, S1).

The results of the bifurcation analysis are presented in the form of the bifurcation diagram. For the one-parameter continuation, the diagram contains the measurable entity (e.g. variable level or period) as a function of the parameter. On our one-parameter diagrams, we show the level of the protein *P*_{1} (maximum level in case of oscillations) as a function of the corresponding parameter. Note that the stable/unstable solutions are shown in solid/dashed lines.

For the two-parameter continuation of the system solutions, bifurcation diagrams contain one parameter as a function of the other at the bifurcation point of the system.

## 4. Deterministic results

The system with explicitly modelled dimerization and TF (un)binding shows the emergence of oscillations via the superAHB for two different large values of the TF binding rate to the promoter site *k*_{r} (figure 2*a*). This is the expected behaviour of the repressilator as was shown in many other studies (e.g. [24,25]).

Slower searching process results in the qualitative change of the dynamics (figure 2*b*): the subAHB provides for limit cycle oscillations as well as for the coexistence (hysteresis) of the limit cycle and a branch of the stable steady-state solutions. The hysteresis is observed in a broad range of the bifurcation parameter values, which is further assessed by the two-parameter continuations of the AHB and LPC bifurcation points (figure 3).

Figure 3 shows the multi-dimensional character of super- → subAHB transition in the parameter space, i.e. transition from single-attractor dynamics to the hysteresis with the two dynamical regimes present. Although some parameters can compensate for the effect of slow TF–promoter interaction on the type of the bifurcation, the coexistence region of the steady state and the limit cycle remains significantly large. We also have to note that although the number of plasmids *g*_{tot} can move the boundaries of the bifurcation curves, the qualitative results presented in figures 2 and 3 hold (see insets in figure 3).

We have not been considering the free dimer degradation up to this point. This simplification is widely used in the simulation studies of molecular networks, because the QSSA approximation assumes one-to-one correspondence between kinetic characteristics of monomer and oligomer forms. We hypothesize that dimer degradation is much slower than that of the monomeric proteins based on the indirect evidence reported in [21,22]. Here, we represent the dimer degradation as a fraction (*δ*) of the degradation rate of monomers (*d*_{p}).

To study the effect of the dimer degradation in the limits of our model, we compute the two-parameter continuation of the AHB bifurcation at the *r*_{m}-*δ* parameter plane. In the case of strong and fast binding of TFs, the system shows superAHB up to very large values of *δ* covering the whole interval of realistic rates of the dimer degradation (figure 4). As expected, for strong adiabaticity and small enough dimer degradation rate *δ*, the system demonstrates the emergence of the limit cycle via the superAHB. Additionally, the region of oscillations expands as the dimer degradation tends to zero (figure 4). Similarly, in [22], the authors arrived at the same conclusion.

Conversely, if the binding rate of TFs further decreases, the continuation over transcription rate *r*_{m} (figure 5) reveals an unexpected limitation of the limit cycle branch. When moving from low to high *r*_{m} values, the continuation shows the emergence of the limit cycle from superAHB and its disappearance through LPC bifurcation resulted from the subAHB. Accordingly, the AHB bifurcation curve on the *r*_{m}–*δ* plane dramatically changes as shown in figure 6.

Interestingly, the non-zero values for the dimer degradation restore the subcritical character of AHB boundary, which is otherwise super-critical. Note the hysteresis region for *δ* > 0 in figure 6, although the *k*_{r} values are from the superAHB parameter set (given there is no dimer degradation).

Thus, a slight increase in the dimer degradation rate effectively enlarges the hysteresis region between subAHB and LPC curves along the transcription rate *r*_{m}. However, there is an upper limit for *δ*, after which the circuit functions only at the steady-state regime, that is, the limit cycle disappears through the Hopf–Hopf bifurcation (HH). Therefore, the more non-adiabatic the oscillator, the lower the intensity of the dimer degradation to sustain the oscillations (figure 6). To the best of our knowledge, there are no special experimental reports on dimeric TF degradation. However, some studies indirectly suggest stronger thermodynamic stability of the dimers (e.g. [21,22]), which, in turn, favours the oscillatory behaviour for the system studied here.

## 5. Stochastic simulations

We studied the effect of noise buffering by increasing the number of protein/mRNA molecules. For this, we have simulated the stochastic system varying the degradation rate of the proteins *d*_{p} and the transcription rate *r*_{m} for the adiabatic conditions. The results of the study are summarized in electronic supplementary material, table S1, and typical distributions and time series are shown in electronic supplementary material, figures S1 and S2.

From this, we conclude that the detailed model (containing explicit dimerization and TF–promoter interaction reactions) of the repressilator circuit demonstrates the noise buffering property similar to many other genetic and molecular systems as a result of the increased number of constituting molecules, e.g. proteins and mRNA [11,19].

In this study, we are to maintain the superAHB type of transition from steady state to oscillation. This enables the noise-buffering property to be studied without mixing with the effect of hysteresis. However, when varying the degradation rate, this is not fulfilled and the system demonstrates the subAHB transition for some parameter values. Nevertheless, the noise buffering property is persistent even for the subAHB type of the transition (electronic supplementary material, S2). We argue that this is due to the cyclic repressive feedbacks inherent in the repressilator, making the oscillatory regime dominant in the hysteresis region (see Discussion).

### 5.1. Hysteresis

The stochastic time series shown in figure 7*a*,*b* demonstrates a transition from stationary to a two-stable hysteresis region over transcriptional rate *r*_{m} along with the corresponding deterministic trajectories. Additionally, the figure contains the noise-induced modulations around the stable steady state before superAHB (figure 7*c*). The latter is served for the sake of comparison with the modulations around the steady state for the subAHB conditions (figure 7*b*).

It is seen from figure 7 that there is alternation between small-amplitude fluctuations around the steady state and larger amplitude oscillations before the hysteresis region as determined by the ODE analysis (figure 7*b*). We call the latter the *noise-induced limit cycle*, which emerges due to the stochasticity of the system and not having the deterministic counterpart. The time series of both small and large fluctuations have slow frequency modulated amplitudes, which reflect the slow process of dimer binding. Oppositely, the *stochastic limit cycle* emerges due to the very dynamical properties of the system and corresponds to the deterministic limit cycle (figure 7*a*).

The period distributions corresponding to the time series from figure 7 are shown in figure 8. The distributions are highly skewed, i.e. have prolonged right tails, due to the longer times the system spends in the vicinity of the steady state, thus contributing to the longer periods (details on the skewness study of the distributions can be found in electronic supplementary material, S3). However, before subAHB (deterministic steady state), the distribution is bimodal (figure 8*a*, dark), whereas before superAHB the distribution shape is exponential (figure 8*b*). This reflects the fact of the slow frequency modulation of the amplitude in the time series before the hysteresis region (subAHB conditions), while superAHB conditions are devoid of this feature resulting in the exponential period distribution.

Furthermore, the hysteresis time series (figure 7*a*) mainly contains stochastic limit cycle regions rarely interspersed with jumps to the steady state, thus, the modulation and, hence, bimodality disappear (see figure 8*a* in grey). However, the distribution for this case is still significantly skewed, thus, we can find periods several times larger than the period of the stochastic limit cycle oscillations (which closely corresponds to the deterministic period shown with arrow in figure 8).

Before AHBs, the shapes of the period distributions (figure 8) are clearly distinct for the two conditions leading to oscillations. Namely, the superAHB scenario reveals an exponential-like distribution (figure 8*b*) and the subAHB scenario shows bimodality (figure 8*a*, dark).

Finally, within the oscillatory region one can expect that the superAHB conditions will have symmetric period distribution close to the normal one (see, for example, distributions in electronic supplementary material, figure S2), due to the equal probability of getting larger and smaller values than the mean period and the absence of the distortions from other stable attractors. However, the subAHB distribution is highly skewed whenever there are two stable coexisting attractors (figure 8*a*, grey).

## 6. Discussion

The repressilator is a prototypical system that has been studied as an example of the genetic ring oscillator. It is included in the list of typical genetic and biochemical motifs [26] and is even considered as a core element for a circadian oscillator [27].

The standard mathematical model of the repressilator, following the QSS approximation and based on the assumption that repressor–promoter binding is fast [3], demonstrates some conditions for the oscillations to emerge; however, it is unknown, if other conditions, leading to oscillations, can be possible. Moreover, this assumption is not definitely valid in the light of the recent experimental measurements of searching times [16].

We have studied the dynamics of the full (unreduced) model of the repressilator and found that slow repressor dimer binding (non-adiabatic conditions) changes the bifurcation type from superAHB to subAHB, which indicates the coexistence of stable attractors (limit cycle and steady state) and the possibility for alternation between the two in response to various stimuli. This novel genetic switch occupies a large region in the parameter space and the higher the degree of the repressilator's non-adiabaticity, the wider the hysteretic region.

Slow binding requires more repressors to support sustained oscillation, hence, the role of the free dimer degradation becomes more important for the dynamics than it was thought before. There are two connected effects that should be considered: (i) even a small rate of dimer degradation (as compared to monomer degradation) significantly extends the parametric region of the coexistence and (ii) the dimer degradation rate should be limited to avoid full suppression of oscillations. Perhaps, the trade-off between the degree of non-adiabaticity and the dimer degradation is the new condition for the oscillatory function of the repressilator.

It must be stressed here that we consider the protein degradation as a result of the enzymatic activity only. However, some estimates of the monomer and dimer degradation rates are based on ‘dilution’ due to the cell growth and division. Our analysis indicates (data not shown) that the further increased values of *δ*, accounting for dilution due to cell growth, preserve the oscillatory and hysteretic dynamics of the circuit. However, considerations of cell growth and division must be more systematic: along with the dimer dilution, other factors should be considered as well (the replication-related increase in plasmid copies, variation of the free RNA polymerase concentration, molecular crowding of TFs, etc.). Thus, we leave the question for further investigation.

The stochastic dynamics of the new genetic switch may be very different from the noise-induced behaviour of the classical toggle switch [28]. The latter has two stable fixed points in the phase space in contrast to the former, whose one of two attractors is the closed cycle, which changes the conditions for switching. The coexistence of the cycle with the stable steady state raises expectation of fast switching, which, however, is not so frequent. The phase variables describing different genes at the multi-dimensional cycle are separated in the phase space due to the three alternating repressive interactions that are inherent in the repressilator. Therefore, intrinsic noise is not very efficient in throwing the phase point across the separatrix, which is a 12-dimensional unstable limit cycle emerging in the subAHB conditions. This hypothesis is further supported by the stochastic simulations with the reporter plasmid that increases the asymmetry of the circuit (see Discussion and electronic supplementary material, S4).

Additionally, the shapes of the period distributions are very different in the case of the two scenarios for oscillation emergence. The differences are very distinct and can be understood by investigating the topology of the phase space.

In the hysteresis case, the distribution is highly skewed due to the stochastic oscillations interchanging with rare leaps to the steady-state fluctuations (see figure 7*a* and figure 8*a*, grey). Additionally, in the vicinity of the subAHB, the system exhibits low-frequency modulations around the fixed point attractor, which lead to the bimodality in the period distributions (figure 7*b* and figure 8*a*, dark). This is stipulated by the existence of the full-amplitude limit cycle immediately available after the subAHB bifurcation point, whereas one has to move farther from the superAHB to get the full-amplitude oscillations. The latter condition, in turn, makes the period distribution around superAHB exponential by shape, which becomes more symmetric (close to normal distribution) as the system moves away from the bifurcation point deeper into the oscillatory region.

As stressed in the original paper [3], only 40% of the cells with the repressilator plasmids inside demonstrate oscillations. We speculate that in the limits of the model (2.1), it may be a result of the long time the cells spend in the stable steady state while in the hysteresis region. Furthermore, this is supported by the fact that the slower TF binding rate values are more realistic, indicating that the coexistence state is more probable in the real system.

On the other hand, the stochastic modulations have large enough amplitude around the steady state in the hysteresis region (as compared to the amplitude in the pure oscillatory area), hence, these modulations must be discerned in the experiment as well. However, the small number of cycles observed in the experiment (three to six cycles [3,29]) may not capture the transition regions completely, i.e. the regions where the dynamics changes from modulations around the steady state to the full-amplitude oscillations. Interestingly, the bimodality of the period distributions has been observed in other experimental studies on the repressilator, though the explanation is different [29].

Additionally, to emulate the real experimental set-up as reported in [3], we have conducted stochastic simulations of the repressilator augmented with the reactions on the reporter plasmid, containing the *gfp* gene (see the details of the system in electronic supplementary material, S4). The main result of this study is the increase in the switching frequency between the full-amplitude oscillations and fluctuations around the steady state. This takes place due to the additional outflow of one of the dimers of the repressilator to the reporter gene, which increases the asymmetry of the repressilator. Thus, we conclude that introducing the reporter plasmid favours the observation of the hysteresis (electronic supplementary material, S4).

Finally, we cannot conduct thorough quantitative comparison between the experiment in [3] and model (2.1) due to the difficulties arising from the assessment of the parameters (for example, exact number of plasmids *g*_{tot} or protein degradation rate *d*_{p}), which may significantly move the bifurcation boundaries. We also speculate that this task might be accomplished using the inference algorithms, given different shapes of the period distributions for the two distinct cases studied here. Nonetheless, the unreduced model (2.1) opens new perspectives for understanding the non-adiabatic repressilator, which has been used as a platform for addressing many problems arising in studies of genetic oscillators.

## Funding statement

This work was partially supported by grant nos. RFFI 12-02-00529 and 14-01-00196.

## Acknowledgement

We are thankful to two anonymous reviewers of this manuscript, whose comments helped to improve the quality of the work.

- Received November 28, 2014.
- Accepted January 2, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.