Improvement of the memory function of a mutual repression network in a stochastic environment by negative autoregulation

Background Cellular memory is a ubiquitous function of biological systems. By generating a sustained response to a transient inductive stimulus, often due to bistability, memory is central to the robust control of many important biological processes. However, our understanding of the origins of cellular memory remains incomplete. Stochastic fluctuations that are inherent to most biological systems have been shown to hamper memory function. Yet, how stochasticity changes the behavior of genetic circuits is generally not clear from a deterministic analysis of the network alone. Here, we apply deterministic rate equations, stochastic simulations, and theoretical analyses of Fokker-Planck equations to investigate how intrinsic noise affects the memory function in a mutual repression network. Results We find that the addition of negative autoregulation improves the persistence of memory in a small gene regulatory network by reducing stochastic fluctuations. Our theoretical analyses reveal that this improved memory function stems from an increased stability of the steady states of the system. Moreover, we show how the tuning of critical network parameters can further enhance memory. Conclusions Our work illuminates the power of stochastic and theoretical approaches to understanding biological circuits, and the importance of considering stochasticity when designing synthetic circuits with memory function.


Background
Memory is ubiquitous in biological systems [1][2][3][4]. Characterized by a continued response to a transient stimulus [5], cellular memory has been found to aid in the robust control of diverse biological functions such as synaptic plasticity [6], differentiation [7,8], cell cycle transition [9], or gene regulation [10]. Memory in cellular circuitry is tightly linked to bistability, i.e. the presence of two stable steady states [11,12]. To this end, memory is achieved if an input signal evokes a switch to an alternative steady state where the system remains over time even after the input signal has disappeared [11,12].
Our general understanding of design principles of biological networks has improved dramatically during recent years [13]. Feedback loops as elementary components provide common control mechanisms of cellular networks [14]. For instance, negative feedback can facilitate adaptation and oscillation [15], while positive feedback loops, which play pivotal roles in cellular signaling [16], often promote signal amplification, bistable switches [17] and memory [11]. Indeed, a common circuit architecture that is known to give rise to memory is based on interlinked positive feedback network loops [18,19]. Several experimental and theoretical analyses of such network architectures that can give rise to memory have been reported [7,12,[20][21][22][23][24][25]. These studies have revealed many general properties underlying successful memory, for instance that ultrasensitivity may be sufficient to generate two stable states [26], and that the transition periods between the bistable states compose an important characteristic of a cellular memory module [19]. In all cases, bistability and memory can convert a transient input signal to a continued cellular response [27,28]. This property is directly visible from the characteristic hysteresis behavior observed for bistability and memory systems. Experimental examples confirm these findings, such as the presence of hysteresis behavior and bistability in a synthetic genetic circuit in Escherichia coli [18]. Based on these general principles that link network architecture to systems behavior, the blossoming field of synthetic biology has contributed many useful designs of genetic circuits based on a quantitative understanding of biological networks [29][30][31] with myriad implications for biotechnology, biocomputing and gene therapy [13,26,32]. However, the rational design of a robust memory function can be implemented via many alternative mechanisms that remain incompletely understood [21,29].
Importantly, the vast majority of genes in a cell is normally only expressed at very low levels. Due to the inherent randomness of transcription and translation, this gives rise to substantial stochastic fluctuations in cellular mRNA and protein copy numbers of lowly expressed genes [33][34][35][36][37][38][39][40]. Accordingly, the analysis and design of biological circuits cannot be based solely on deterministic properties of the network topology alone. In response, much progress has been made in better understanding the role of stochasticity to network function [41][42][43][44][45]. For instance, careful work has delineated the sources of intrinsic noise in small transcriptional networks [46] and how noise may propagate through network architectures [39,41], as well as how feedback loops regulate bistability [42,47] and intrinsic noise [48]. Similarly, extensive work has established general principles underlying the design of genetic circuits that exhibit bistability and memory [11,[42][43][44][45].
Recent examples of opposing behavior further highlight the importance of considering the contribution of stochasticity to cellular circuitry. On the one hand, it has been shown that noise can induce multimodality and even stochastic memory in a system that, according to a deterministic description, lacks bistability [49,50]. On the other hand, stochastic fluctuations in gene expression levels often reduce or disrupt the memory function of biological networks [21,51,52]. Equally importantly, seminal technical and conceptual advances have made strong progress in the efficient and accurate approximation of multivariate nonlinear stochastic systems [53,54].
Here, we investigate how a negative autoregulation network architecture can improve the sustained memory function of a mutual repression network in a stochastic environment. Our work extends the previously published regulated mutual repression network (MRN) [12] to the regulated mutual repression network with negative autoregulation (MRN-NA) to investigate the effect of a negative autoregulation loop on memory function. The network architecture is characterized by a mutual repression cycle that can give rise to bistability and memory, adopted from two well-characterized mutual repression networks, the system consisting of LacI and TetR in E. coli that displays a bistable gene expression memory module [24], and the mutual repression of the two repressors cI and Cro that yield a bistable memory module in the bacteriophage λ [21].
To explore and compare the memory behavior of the MRN-NA network model, as well as how it compares to the MRN model, deterministic rate equations, stochastic simulations and theoretical analyses of Fokker-Planck equations were employed to identify principles of the robustness of the memory function. We demonstrate how negative autoregulation can reduce intrinsic noise and thus improve the memory function by increasing the stability of the steady states. Negative autoregulation can for instance be achieved by the ability of proteins to bind and sequester their own mRNA. Systematic analyses of successful memory as a function of central model parameters that describe the mutual repression cycle highlight principles and limits of the memory function in these mutual repression networks. Moreover, our work emphasizes the importance of considering stochasticity when designing synthetic circuits with memory function.

Results
To investigate the effect of negative autoregulation and stochasticity on network memory function, we constructed a model network following a previously established graphical notation [55,56] (Fig. 1). Specifically, we modelled a mutual repression network with negative autoregulation (MRN-NA) that extends our previously published mutual repression network (MRN) [12]. The network consists of proteins y(1), y(2) and y(3). An input signal S induces the synthesis of y(1), which activates and represses the syntheses of y(2) and y(3), respectively. The syntheses of y(2) and y(3) are mutually repressed with cooperativity. Moreover, negative autoregulation controls the synthesis reactions of y(2) and y(3) (Fig. 1). The addition of negative autoregulation permits to investigate how negative feedback loops may affect memory. While modulation of several parameter values should tune the presence of memory, the present analysis focuses on a direct back-to-back comparison of the two network architectures that only differ in the addition of negative autoregulation circuits. To compare the memory regions of MRN-NA and MRN networks across our deterministic, stochastic, and theoretical analyses, the values of the corresponding kinetic parameters as well as the steady state levels of y(2) and y(3) were conserved as much as possible in discrete parameter space [12] (Table 1; Additional file 1: Texts S1, S2) (see Methods).

MRN-NA model exhibits robust memory
Similar to our previously reported results on the memory function of the MRN [12], the MRN-NA network exhibits robust memory functionality that depends on strong cooperativity in the mutual repression cycle, as expressed by high Hill coefficients n. Thus expectedly, the deterministic solution to the time evolution of the MRN-NA model showed persistent memory after the end of the period during which the signal S was applied for both Hill coefficients of n = 7 and n = 8 (Fig. 2a). Strong hysteresis behavior was observed for both y(2) and y(3) levels as a function of increasing and decreasing S (Fig. 2b, c; Additional file 1: Figure S1). In the deterministic case, bistability is sufficient to generate memory. Stochastic simulations of the MRN-NA model as described by the birth and death processes of Eqs. 1-3 revealed strong and persistent memory both with Hill coefficients n = 7 (Fig. 2d) and n = 8 (Fig. 2e). Importantly, the MRN alone could not sustain memory for either y(2) or y(3) at n = 7 in a stochastic context [12], but the MRN-NA yielded strong and sustained memory due to the addition of negative feedback loops ( Fig. 2d; Additional file 1: Figure S2, S3).
To further understand how the addition of negative autoregulation increases the robustness of the memory function, we sought to quantify and compare the intrinsic noise in the protein levels of y(2) and y(3) in the MRN and MRN-NA models. It is well established that noise or a stochastic perturbation can flip gene expression levels from one to the other steady state [40,46]. To obtain a reliable and comparable estimate, we computed the coefficients of variation (CVs) of y(2) and y(3) during the active signal S period; the interval from simulation step 270 to 500 was chosen to allow the system to respond to the stimulus for 20 simulation steps. Our analysis of noise in protein levels y(2) and y(3) as a function of the dissociation constants K(2) = K(4) confirmed a systematic reduction of stochasticity in y(2) and y(3) levels upon the addition of negative autoregulation (Fig.  2f). While maintaining the same steady state levels, noise was dramatically reduced for y(3) (Fig. 2f). Of note, the lower steady state y(3) is more susceptible to noise. Similarly, intrinsic fluctuations in y(2) during the signal period were reduced in the MRN-NA compared to the MRN model in our stochastic simulations irrespective of the dissociation constants tested (Fig. 2f). Even more important to the stability of the steady states are the intrinsic fluctuations after the end of the signal period. The levels of intrinsic noise are overall slightly higher, likely because the signal does not stabilize gene expression levels any longer. However, our corresponding analysis reveals that the negative autoregulation reduces, as expected, variability in the system also after the stop of the signal ( Fig. 2g; Additional file 1: Figure S4). Taken together, our results suggest that the negative autoregulation played a vital role in  Tables 1 and 2. In this work, the MRN-NA is compared to the previously published MRN [12] that is identical with the exception of the negative autoregulation loops. In addition, a simplified wiring diagram is shown Table 1 List of kinetic parameters increasing the persistence of stochastic memory by reducing intrinsic noise.

Negative autoregulation enhances the memory region
Our analyses show that the MRN-NA model can exhibit memory, which however critically depends on persistent mutual repression with strong cooperativity. We next set out to systematically map its memory region as a function of Hill coefficient and dissociation constants K(2) = K(4) in comparison to the previously reported MRN [12] (Fig. 3). Successful memory was defined as sustained levels of y(2) and y(3) at or near their levels during the applied S after the end of the input signal S at simulation step 500 until the end of the simulations at simulation step 1000. Moreover, for direct comparison, the kinetic parameter k(7) was adjusted to conserve the steadystate levels between the two models [12] (Additional file 1: Text S1). Indeed, the high steady-state levels of both the models were set to the same expression values, while the low steady-state levels were tuned be as similar as possible. Memory was assessed in both the deterministic and stochastic context (see Methods). Accordingly, two memory regions were defined: deterministic memory and stochastic memory (Fig. 3).  Table 2). The simulated time evolution is shown for y(2) and y(3) at n = 7 (red and magenta lines), and at n = 8 (black and cyan lines). b The hysteresis curves of y(2) at different Hill coefficients n = 7 and n = 8 are consistent with the numerical integration of the rate equations. c The hysteresis curves of y(3) at different Hill coefficients n = 7 and n = 8 are consistent with the numerical integration of the rate equations. d Trajectories of the stochastic simulation of y(1), y(2) and y(3) at Hill coefficient n = 7. e Trajectories of the stochastic simulation of y(1), y(2) and y(3) at Hill coefficient n = 8. f, g Stochastic fluctuations in the MRN and MRN-NA models during the interval from simulation step 270 to 500, and F during the period from simulation steps 500 to 750. g The coefficients of variation (CVs) from the simulated stochastic trajectories during the signal period as a function of changing dissociation constants K(2) = K(4)at Hill coefficient n = 8. For more optimal comparability between the two models, the parameters associated to the negative autoregulation reactions were tuned to conserve the high steady state levels between both models. Shown are the CVs of y In both the MRN and MRN-NA models, deterministic memory could be observed for Hill coefficients of n = 4 or greater (Fig. 3). In contrast, stochastic memory not only required higher Hill coefficients, but also posed different limits on the MRN and MRN-NA models. Without negative autoregulation that suppresses noise, stochastic memory in the MRN model required a Hill coefficient of n = 8 or greater (Fig. 3). In the presence of negative autoregulation that reduces noise in protein levels, stochastic memory in the MRN-NA model could be observed from n = 7 (Fig. 3). Thus, the negative autoregulation relaxes the need for strong cooperativity to generate stochastic memory. Furthermore, the stochastic memory region for the MRN-NA model was found to be larger than for the MRN model (Fig. 3, red area), indicating that negative autoregulation of the MRN-NA enhances memory function by reducing intrinsic noise (Fig. 2f, g). While the deterministic memory regions in both networks are essentially the same (Fig. 3, green area), it was clearly more difficult to maintain a memory effect under noise in the MRN than the MRN-NA (Fig. 3, red area). The large discrepancy between the deterministic and stochastic memory areas indicate the general challenge of maintaining memory function in the context of stochasticity.

Stochastic potential of the MRN-NA
Prerequisite for memory is normally the presence of bistability in the system. Thus, a limited bistable regime could explain why such high levels of cooperativity in the mutual repression between y(2) and y(3) were required to yield memory. To map out the bistable regimes of the MRN-NA in comparison to the MRN [12], theoretical analyses based on the chemical master equations were employed to support the numerical simulation results.
First, the corresponding Fokker-Planck equations Eqs. 9-13 describing the time evolution of the probability densities for y(2) and y(3) levels were derived from the three rate equations Eqs. 1-3 (see Methods). We made use of the quasi-steady-state assumption to derive a onevariable equation for each system [12] (Additional file 1: Text S2). While there are inherent limitations to it [51], the quasi-steady-state assumption is widely used [46]. Since we focus on the steady-state at t → ∞, applying the quasi-steady-state assumption to y(2) or y(3) is mathematically reasonable. By using the Fokker-Planck equations, we next estimated the two-dimensional stochastic bistable regions, characterized through the presence of a double well potential, for the MRN and MRN-NA models as a function of the Hill coefficient and dissociation constants [12] (Fig. 4; Additional file 1: Figure S5; Texts S1-3).
In these theoretical analyses, both the MRN and the MRN-NA models permitted bistability starting from a Hill coefficient of n = 3 for both y(2) and y(3) (Fig. 4). However, the bistable regions in the MRN-NA were larger than those in the MRN model (Fig. 4). To display bistability in a stochastic context, the MRN model required dissociation constants K(2) = K(4) within a more constrained parameter space (Fig. 4). These analyses established that the MRN-NA more readily exhibits stochastic bistability than the MRN. Conversely, a limited bistability regime in the MRN in a stochastic context (Fig. 4) contributes to explaining why this network had a smaller stochastic memory region (Fig. 3). end, the stability of a steady state of a stochastic system can be estimated by the mean first-passage time (MFPT). The MFPT quantifies the average number of simulation steps it takes for a system to leave a favorable steady state due to random events. As such, the MFPT provides a useful description of the time-scale on which a phase transition is likely to happen [57][58][59]. Thus, because the presence of a stochastic bistable region only indicates the bimodality of gene expression as prerequisite for memory but not the presence of memory itself, a MFPT analysis can be helpful for identifying the precise conditions under which memory can be attained in a stochastic context. Here, MFPT analyses were performed between the two stable steady states under the quasi-steady-state assumption [12] (Additional file 1: Figure S5 Due to the asymmetry in the both MRN and MRN-NA models (Table 2), the MFPTs for y(2) and y(3) showed opposing behavior, as expected (Fig. 5). While the MFPT for y(2) to leave the lower steady state, T L ðyð2Þ st l Þ, is consistently longer than the corresponding time to leave the upper steady state, T U ðyð2Þ st u Þ, the system is more likely to leave the lower steady state of y(3) as indicated by T L ðy ð3Þ st l Þbeing much smaller than T U ðyð3Þ st u Þ (Fig. 5). An increase in the Hill coefficient as measure of cooperativity exacerbated this trend observed for the MFPTs of y(2) and y(3) (Fig. 5a). The prolonged MFPTs as a result of increasing Hill coefficients explained why the robustness of sustaining stochastic memory improved with increasing cooperativity. For the parameter space explored, the lower steady state of y(2) was more persistent in both models for any Hill coefficients, as evident by T L ðyð2Þ st l Þ > T U ðyð2Þ st u Þ (Fig. 5a). Moreover, all MFPTs for the MRN-NA, T L ð yð2Þ st l Þ , T U ðyð2Þ st u Þ , T L ðyð3Þ st l Þ , and T U ðyð3Þ st u Þ , were observed to be longer than their equivalents for the MRN model, for any Hill coefficient (Fig. 5a). These results confirm an increased persistence of the steady states upon introduction of negative autoregulation in MRN-NA compared to MRN, which explains the improved memory function in the MRN-NA system. Specifically, the residence times at the upper and lower steady states of y(2) and y(3) were extended to achieve the persistent memory after the stop of the signal.
Next, we investigated the effect of the dissociation constants in the mutual repression cycle on the MFPTs of y(2) and y(3) (Fig. 5b). Consistently, the MFPTs for all steady states, T L ðyð2Þ st l Þ, T L ðyð3Þ st l Þ, T U ðyð2Þ st u Þ and T U ð yð3Þ st u Þ ,were found to be longer in the MRN-NA than the MRN for any dissociation constants (Fig. 5b). For both models, all MFPTs, T L ðyð2Þ st l Þ , T U ðyð2Þ st u Þ , T L ðy ð3Þ st l Þ and T U ðyð3Þ st u Þ gradually decreased with an increase in the dissociation constant, suggesting that strong binding in the mutual repression cycle is necessary for persistent memory. We also studied the MFPTs for the MRN-NA as a function of the negative autoregulation constants k(9) = k(10) (Fig. 5c). The MFPTs for the MRN-NA at all steady states, T L ðyð2Þ st l Þ, T L ðyð3Þ st l Þ, T U ðyð2Þ st u Þ and T U ðyð3Þ st u Þ, were found to increase with increasing strength of the negative autoregulation. In summary, these MFPT analyses illuminated in detail

Probability densities of the steady-state levels in the MRN-NA
The previous analyses revealed a strong dependency of sustained memory on the robustness of the steadystates. To formally establish the probability densities associated with populating the upper and lower steady states of y(2) and y(3), we calculated the probability densities from the Fokker-Planck Equations (Eq. (15)) as a function of the Hill coefficient and dissociation constants (Fig. 6). The probability density of the upper steady state of y(2) increased while that of the corresponding lower steady state decreased with an increase in the Hill coefficient for both the MRN and the MRN-NA models (Fig. 6a). The addition of negative autoregulation in the MRN-NA system visibly made the upper steady state level of y(2) more dominant (Fig. 6a). The probability density of y(3) clearly illustrated the evolution of bistability and subsequent increased probability densities of the lower steady states with increasing Hill coefficient (Fig. 6b). At low Hill coefficient (n = 3) the MRN lacked a pronounced probability density of a lower steady state (Fig. 6b, solid blue line) while in the MRN-NA the probability density did not separate well the peaks at lower and upper levels of y(3) (Fig. 6b, dashed blue line). This result of weak bistability was consistent with our previous analysis of the bistability regions (Fig. 4). In turn, higher Hill coefficients of n = 8 and n = 12 promoted first the emergence and then separation of two clearly populated steady states (Fig. 6b). While the upper steady state of y(3) was dominant for n = 12 the system switched to a dominant lower steady state for y(3) (Fig. 6b). Moreover, the addition of negative autoregulation strongly promoted the population of the lower steady state at all Hill coefficients analyzed (Fig. 6b, dashed lines). By buffering intrinsic noise, this analysis thus confirmed that the negative feedback loops stabilized particularly the low expression states of y(3), i.e. those that are most susceptible to fluctuations.
Similarly, the probability densities of y(2) and y(3) were found to quantitatively and qualitatively also depend on the dissociation constants K(2) and K(4) in the mutual repression cycle (Fig. 6c, d). At low K(2) = K(4) = 15, the upper steady state of y(2) was dominant in the MRN and MRN-NA models with slightly higher probability density for the MRN-NA model (Fig. 6c). In general, with decreasing values of K(2) and K(4) the upper steady states became more pronounced in both models, even more so in the MRN-NA (Fig. 6c). In turn, the probability density of the lower steady state of y(3) became dominant with decreasing the dissociation constants in both the models (Fig. 6d). Moreover, lower values for K(2) and K(4) as well as negative autoregulation favored more pronounced bistability (Fig. 6d). The MRN-NA model increased the probability density of the upper steady-state of y(2) and the lower steady-state of y(3) more than the MRN model. This result can support the persistent memory of the upper steady-state of y(2) and lower steady-state of y(3) after the stop of the signal (Fig. 6d). Taken together, these results suggest a

Consistency between stochastic simulation and the Fokker-Planck equations
To verify the quasi-steady-state assumption for our case, we validated that the solutions to the Fokker-Planck equations were consistent with the probability densities obtained by stochastic simulation of the full system ( Figs. 7 and 8). The Fokker-Planck equations provided almost the same probability density as the Gillespie stochastic simulation for initial molecule numbers ofy(1) = y(2) = y(3) = 1 (Figs. 7 and 8). The main limitation of the quasi-steady-state assumption is an over-estimation of the separation of the steady states ofy (2). However, in the MRN-NA y(3) is more susceptible to fluctuations and the driver for switching between steady states. Consequently, the results obtained under the quasi-steadystate assumption are meaningful for investigating the memory function in these mutual repression networks.

Discussion
By applying numerical integration of the rate equations, stochastic simulations and theoretical analysis of the Fokker-Planck equations, we investigated if negative autoregulation can improve the memory function in a mutual repression network. Previous work had established the regulated mutual activation network (MAN) and mutual repression network (MRN) as good model systems to study fundamental properties of cellular memory [12]. Here, we have systematically decoupled contributions of additional negative autoregulation (MRN-NA) to the persistence of the memory functionality. In general and as expected, stochasticity decreased the memory function in all models. The addition of negative autoregulation in the MRN-NA however extended robust and persistent memory in both the deterministic and stochastic approaches to lower levels of cooperativity in the mutual repression cycle. Our results thus suggest that, in addition to the added negative autoregulation, the robustness of the stochastic memory of the MRN-NA could be further improved by increasing the binding strength of the repressor proteins.
In the present work, select sets of kinetic parameters that can give rise to memory were used instead of exhaustive searches for alternative parameter combinations. Despite accurately simulating the time evolution of protein concentrations in stochastic systems, the Gillespie algorithm is known for its limited power in characterizing sustained memory in more complex networks due to its computational cost that can quickly become prohibitive. To generalize some of our findings, phase diagrams of the memory region were derived as a function of the two most critical and previously identified parameters. However and as evident from the small stochastic memory regions, finely-tuned parameter combinations were necessary to achieve robust memory. This observation suggested both a biological challenge of implementing a robust memory circuit, as well as a theoretical challenge in rationalizing its mechanistic basis through exhaustive parameter perturbation.
Theoretical approaches that can complement these shortcomings and make accessible the exploration of larger parameter ranges often rely on additional assumptions. Notably, the quasi-steady state assumption of separately analyzing steady states is widely used, performs in general well and is reasonably justified in the presented work. Seminal advances have now provided advanced approximation methods that start to render even multivariate and nonlinear chemical masters equations and related Fokker-Planck equations amenable to theoretical analyses. For instance, the recently introduced linear-mapping approximation converts a nonlinear system to a linear problem via a mean-field approximation [54]. This and related techniques will vastly improve accuracy and computational efficiency for future analyses of stochasticity in more complex networks. Biology has certainly evolved more complex architectures, but in most cases the trade-offs between network complexity and function remain poorly understood. For instance, while the MAN network was found to yield robust memory based on mutual activation [12], comparable results for the mutual repression network could only be obtained after addition of negative autoregulation, i.e. by a more complex network architecture. More complex network architectures are often more expensive to maintain, but may yield advantages in fidelity of their robustness and dynamic control. To this end, several related cellular networks of increased complexity have been discovered and await further characterization. For instance, a Notch-Delta mutual repression network serves to communicate between neighboring cells [60] where an increase in Notch activity within a cell decreases that in its neighboring cell. The Notch-Delta mutual repression provides inhomogeneous or opposite protein synthesis in homogeneous cell populations that however depends on spatial changes in gene expression.
Finally, recent papers have outlined many fundamental principles of how to achieve bistability in small networks [61][62][63]. The present work here extends these findings by a comparative network analysis that delineates the effect of additional negative autoregulation on cellular memory.

Conclusions
We have shown that the addition of negative autoregulation can reduce intrinsic noise and generate persistent stochastic memory in a mutual repression network. Our mathematical comparison and theoretical analyses contribute to an improved understanding of how genetic circuits can encode biological function and may aid the rational engineering of memory networks for applications in synthetic biology, medicine, and biotechnology.

MRN-NA model
We constructed a simple model of a gene regulatory network that consisted of two genes encoding a transcription factor that we visualized according to a previously established graphical notation [55,56] (Fig. 1). Our mutual repression network with negative  (1), which in turn activated the synthesis of y(2) and repressed production of y(3). Furthermore, the synthesis reactions of both of y(2) and y(3) were mutually repressed with cooperativity. Negative autoregulation governs synthesis reactions of y(2) and y(3). This system can be bistable and exhibits memory. In formulating our model, we made deliberate use of Michaelis-Menten approximations under the assumption that substrate concentrations were in excess and association equilibria quickly attained. The model is described by the following eqs. 1-3: All parameters are described in Table 1.

Systems modeling and determination of successful memory
The time evolution of the gene expression levels of y(1), y(2) and y(3) was simulated both deterministically and stochastically. For a deterministic systems description, the ordinary differential equations (Eqs. 1-3) were solved in MATLAB (Mathworks) with standard solvers. Stochastic trajectories were simulated with the Gillespie algorithm [64]. In all model simulations, the deterministic and stochastic time courses of y(1), y(2) and y(3) were simulated for 1000 simulation time steps. The input signal S was applied from simulation step 250 to 500. Deterministic memory was defined as sustained protein levels in the numerical integration of the rate equations during the subsequent period from simulation steps 500 to 1000 after stop of the signal S. Similarly, stochastic memory was assessed as sustained protein levels in the stochastic simulations in the period from simulation steps 500 to 1000 after stop of the signal S. While arbitrary, a threshold of 1000 simulation time steps yielded a robust and readily accessible criterium to assess memory. Due to the probabilistic nature of the stochastic simulations, robust stochastic memory for a given set of parameters was defined as successful memory if 18 out of 20 stochastic simulations, i.e. 90%, yielded persistent memory [12].

Theoretical model comparisons
All parameters were set as to render the MRN and MRN-NA models as comparable as possible (Tables 1, 2) [12,65]. With the exception of the additional negative autoregulation loops, this meant using the same parameters throughout. The parameters of the negative autoregulation were tuned to conserve the steady state levels of y(2) and y(3) between the MRN and MRN-NA models (Tables 1, 2). Indeed, the high steady-state levels of both the models were set to the same values, while the low steady-state levels to be as similar as possible. Note that it is not possible to also tune the low steadystate levels to the exact same values. Given the asymmetry of the models (activation of y(2) and suppression of y(3) by S), the deterministic steady-state levels of y(2) and y(3) always show opposing behavior: when the steady-state level of y(2) is high, that of y(3) is low and vice versa. Our parameters choices conserved both the low and high steady states of y(2) and y(3) between the MRN and MRN-NA models [12] (Additional file 1: Texts S1, S2).
To systematically identify parameter choices that can give rise to successful memory, phase diagrams of the memory region as a function of the two most critical parameters, the Hill coefficient and dissociation constants in the mutual repression cycle, were computed by successive simulation at different parameter combinations.
Intrinsic noise in gene expression was quantified by computing the coefficient of variation (CV) of the levels of y(2) and y(3) during the input signal S from the stochastic simulations. To eliminate transition effects, we considered the period from simulation step 270 to 500, i.e. omitting the first 20 simulation steps of the signal period (Fig. 2f). In the same manner, we estimated the CVs of the levels of y(2) and y(3) after stop of the signal, for the period from simulation step 500 to 750 (Fig. 2g).
Finally, we consider the stochastic potential analysis. The limit of P(y, t) at t → ∞ yields P st (y), the stationary probability density function of y [12,21,59,66], which is given by: where N c is the normalization constant. Eq. (14) can be recast in the form: is called the stochastic potential of f(y) [57,58,66].

Mean first-passage time analysis
In gene expression, the stability of a steady state has to be estimated in the presence of noise. The stability of a steady state of a stochastic system can be estimated by the mean first-passage time (MFPT), which describes the expected time within which the system leaves a stable steady state due to random fluctuations. An equilibrium point can exit from its minimum potential due to the effect of noise. The exit time depends on the specific realization of the random process; this is known as the first passage time. The MFPT is the average of the first passage times over many realizations. In the context of anticipating phase shifts, the MFPT provides a useful tool to characterize the time-scale on which a phase transition is likely to occur. Let us consider y st l and y st u ( y st l < y st u ) as two steady states corresponding to a low and a high protein concentration, respectively, separated by the unstable steady state defining the potential barrier y un b (i.e., the unstable equilibrium point). The basin of attraction of the state y st u extends from y un b to +∞, as it is to the right of y st l . Let T(y) be the MFPT to state y un b starting at y > y un b . T(y) satisfies the following ordinary differential equation [21,59,67,68]: with boundary conditions: