Multiple transcription auto regulatory loops can act as robust oscillators and decision-making motifs.

Response time decides how fast a gene can react against an external signal at the transcription level in a signalling cascade. The steady state protein levels of the responding genes decide the coupling between two consecutive members of a signalling cascade. A negative autoregulatory loop (NARL) present in a transcription factor network can speed up the response time of the regulated gene at the cost of reduced steady state protein level. We present here a multi NARL motif which can be tuned for both the steady state protein level as well as response time in the required direction. Remarkably, there exists an optimum Hill coefficient nop t ≅ 4 at which the response time of the NARL motif is at minimum. When the Hill coefficient is n < nopt , then under strong binding conditions, one can raise the steady state protein level by increasing the gene copy number with almost no change in the response time of the multi NARL motif. Using detailed computational analysis, we show that the coupled multi NARL and positive auto regulatory loop (PARL) motifs can act as an oscillator as well as decision making component which are robust against extrinsic fluctuations in the control parameters. We further demonstrate that the period of oscillation of the coupled multi NARL-PARL dual feedback oscillator can also be fine-tuned by the gene copy number apart from the inducer concentration. We finally demonstrate robustness of bistable dual feedback decision making motifs with multi autoregulatory loop component.

Introduction

Gene regulation at the transcription level is critical to maintain the precise intracellular concentrations of various proteins involved in gene signalling, quorum sensing, genetic networks, bistable switches, synthetic gene circuits and gene memory devices [1], [2], [3], [4], [5], [6], [7]. A large number of transcription factors regulate their own expression via binding their cis-acting regulatory elements, leading to self-regulatory loops [8], [9], [10], [11]. Transcriptional auto-regulatory loops play critical roles in fine-tuning the steady-state spatiotemporal protein concentrations across various cell types [2], [8], [9], [12], [13], [14], [15]. This is important for proper functioning, stable growth, development and differentiation of any organism. Transcriptional autoregulation can be either positive or negative. Positive autoregulatory loops (PARLs) are mainly involved in the transcription memory along cell division apart from their involvement in the robust oscillatory and decision making motifs [7].

In negative auto regulation, the protein product of the gene binds the cis-regulatory element associated with its own promoter and decreases its transcription level. The response time of a gene is defined as the average time required to generate half of its steady state protein concentration [5], [16], [17]. Response time () and the steady state protein levels are the key factors which decide the functionality of a signalling cascade. Response time decides the speed of reaction to an external signal and the output steady state protein level decides the subsequent successful passage of the signal through the cascade. For an unregulated gene expression of a stable protein, the response time can be given as where is the recycling rate of the protein product [1], [2], [3], [4], [5], [6]. When the protein is stable over several generations of the cell, then the response time equals one generation time (gt) since the volume doubles at the time of cell division and the protein concentration becomes half of the original steady state value. This means that . When we measure the time in terms of lifetimes of protein, i.e., as , then in the dimensionless time space [18]. Negative auto regulatory loops (NARLs, Fig. 1A) in the transcription factor networks are known to decrease the response times and noise levels [5], [8], [19], [20] and also act as robust genetic oscillators [21].

Fig. 1

Multi negative autoregulation loops (NARL) design. Here (x, m, p) are the concentrations of protein-bound promoter, mRNA and protein-product, respectively. A. Single NARL. Here the protein product of gene G binds and downregulates its own promoter. The concentration of mRNA is m and the protein concentration is p (mol/lit). The transcription rate is km (mol/lit/s), the translation rate is k (1/s), the recycling rate of mRNA into nucleotides (n. a) is (1/s) and the recycling rate of protein into amino acids (a. a) is (1/s). The rate of binding of n protein molecules (Hill coefficient) with the promoter is k ((mol/lit)− s−1) and the dissociation of promoter-protein complex is described by the dissociation rate constant k (1/s). Here D denotes the free promoter and D denotes the protein bound promoter. Further, d = [D] + [D] is the total promoter concentration. The gene will be turned on only when the promoter is in free form. B. Multiple NARLs motifs. Here there are z identical NARLs embedded into a cellular system. The mRNA and protein products of these NARLs will be added to a common reservoir which in turn negatively autoregulate the promoters of the NARL genes. All the z NARLs are coupled through their promoter state dynamics.

Although NARLs can decrease the response times and noise levels, they suffer from a drastically reduced steady state protein concentration, much below the corresponding unregulated or positively autoregulated gene expression level. For example, the NARL motif that was constructed by Rosenfield et. al. [8] showed almost one-tenth of the unregulated steady state protein level. This decrease also results in a decreased amplitude of the genetic oscillators built solely on NARLs modules. Tuneable dual feedback oscillators with OR type logic control on the respective promoters were constructed recently by coupling NARL and PARL motifs to overcome such flaws of NARLs [22]. In this dual feedback motif, the NARL (repressor) motif negatively regulates the PARL motif and the PARL (activator) motif positively regulates the NARL motif. Here the PARL motif is coupled to the NARL motif via OR logic. In case of OR logic control, the target promoter will have two independent binding sites for two different transcription factor protein molecules. In particular, the promoter of the NARL motif will have independent binding sites for both the activator and repressor proteins. In case of AND logic, the target promoter will have a binding site for the heterodimer formed out of the two regulating protein molecules.

Unlike the self-sustained oscillations of the standard NARL motif, the oscillations produced by the dual feedback motifs were stable only over few cell generations. The cell-to-cell variability across the daughter cells produced upon cell division could be one of the main reasons for the breakdown of the oscillations [23] over cell generations. Remarkably, this instability issue seems to be an inherent property of the chaotic dynamics of the OR type logic control over the promoters upon perturbations in the control parameters corresponding to the promoter state dynamics [21]. One can overcome this chaotic nature of the gene oscillators via designing complicated AND type coupling of autoregulatory motifs which requires the engineering of protein–protein interactions. In these context, there is an urge and demand to redesign the standard NARL motif in to a robust genetic oscillator with reduced response time without compromising on the steady state protein levels.

The steady state protein levels of NARLs can be increased in two possible ways, viz. via 1) reducing the affinity of their protein products towards their own promoters and 2) increasing the binding Hill coefficient. When the affinity of the promoter towards its own promoter is decreased, then the negative feedback strength will be decreased which results in the increased protein levels. However, this will also increase the response times in parallel. Alternatively, increasing the Hill coefficient will also increase the steady state protein levels along with increasing the response times. Dual feedback type motifs suffer from chaotic promoter state dynamics that is highly sensitive to the temporal perturbations in the control parameters. In this article, we present a multi ARLs network design to overcome the reduced protein levels of NARLs without compromising on their response times. We will show that the multi ARL motifs are resistant against extrinsic perturbations in the control parameters and can act as robust genetic oscillators as well as bistable decision making units.

Theory

Multi autoregulatory motifs in the bacterial systems

Let us consider z identical copies of a given NARL present inside a single bacterial cell. Each NARL motif consists of a constitutive gene along with the upstream cis-regulatory module (CRM) where the protein product binds and downregulates its own synthesis. We denote the concentrations of the protein bound promoter, mRNA and protein of i copy of the multi NARL motif (Fig. 1B) as (x, m, p), respectively (measured in mol/lit). The promoter concentration of each copy of the multi NARL module inside the cell is d. Therefore, the total promoter concentration of all the NARLs is zd. The mRNA and protein products of all the gene copies are one and the same and added to a common reservoir, i.e., cell volume. Equilibrium statistical mechanics treatment of such multi copy gene expression systems has been carried out earlier [24], [25], [26], [27]. In these studies, the system was assumed to be at the steady state and the corresponding steady state promoter occupancies by the regulatory TF proteins in the presence of non-specific binding sites of the chromosomal DNA were computed using partition function formalism which were then used to compute the fold-change in the gene expression levels. Here the “fold-change” corresponding to lac-repressor system is defined as the ratio between the level of expression in the presence of repressor and the level of expression in the absence of repressor. The absence of repressor scenario represents the unregulated gene expression. Particularly, earlier studies [24], [25] have shown for a single copy lac repressor system that, where R is the number of repressor molecules present inside the cell, N is the number of non-specific binding sites present on the entire chromosomal DNA and is the repressor binding affinity measured in number of kBT. Here there is no autoregulation and the lac repressor gene that is induced by aTc, negatively regulates the YFP reporter gene. Actually, when the system is far away from the steady state, then the temporal expression trajectory of each copy of the same gene will be different from each other and the extent of deviation from the steady state will vary among those gene copies. Further, equilibrium statistical mechanics treatment of the binding of TFs with their cognate sites on DNA will not capture the underlying non-equilibrium one-dimensional (1D) and three-dimensional (3D) diffusion mediated search mechanism [28]. One should note that the specific and non-specific binding are the integral components of 1D-3D diffusion mediated search processes. When the system is far from the steady state, upon independent and parallel expression of all the z number of NARL genes, the emerging protein products temporally interact with all the promoters in an indistinguishable manner. Under this condition, the differential rate equations associated with the cumulative dynamics of such multi NARL system can be written in the dimensionless form as follows:

In this set of nonlinear ordinary differential equations which are not solvable analytically, the dimensionless variables are defined as follows:

The transcription rate is k (mol/lit/s), translation rate is k (1/s), recycling rate of mRNA is (1/s) and recycling rate of protein is (1/s). The rate of binding of n protein molecules (Hill coefficient) with the promoter is k ((mol/lit)− s−1) which accounts for both specific and nonspecific binding components along with 1D-3D diffusion mediated search mechanism [28] and the dissociation of promoter-protein complex is described by the dissociation rate constant k (1/s). Here we use summing formalism to obtain the overall promoter occupancy corresponding to the multi NARLs motif rather than the partition function formalism used in Refs. [24], [25], [27] that is applicable for the equilibrium systems. In our non-equilibrium model, the overall promoter occupancy by TFs is described the dynamical variable X. The steady state mRNA and protein levels corresponding to an unregulated single gene expression are g and h (mol/lit) respectively. Various dimensionless parameters along with the rescaled time variable associated with Eqs. 1a-c are defined as follows.here the parameter v describes the binding kinetics of protein with its own promoter, w describes the relative stabilities of mRNA and protein, describes the binding-unbinding dynamics of protein with its own promoter and describes the coupling between the protein synthesis and the binding of protein with the promoter. The binding strength will be inversely connected to . Particularly, for lac repressor system, is positively correlated with the site-specific binding energy corresponding to the operator sequence. All the parameters defined in Eqs. 3a-e are the same for all the identical copies of NARLs. We add subscript S to denote the steady states of cumulative variables X, M, P, derived by setting the derivatives in Eqs. 1 to 0, and which can be written as follows:

The first observation of such multi NARLs motif is the rising of the steady state protein levels with respect to the copy number z in a square root manner which is evident from the limiting strong binding conditions as in Eqs. 4a-d. Let us assume that there are z isolated copies of the same NARL gene inside a cell. Under non-interacting conditions, the total steady state protein level will be z times the steady state protein level of the individual NARL motif i.e. the total steady state protein levels will increase linearly with the gene copy number under non-interacting conditions. However, Eq. 4c clearly suggests a square root scaling of the total steady state protein level of multi NARL motif with the copy number z under open and interacting conditions. In case of non-autoregulatory systems similar to the one considered in Ref. [24], one can conclude that the overall steady state protein level will linearly increase with z especially when the number of repressor protein molecules is higher than the gene copy number. For example, the fold change level increases over 100 times for 100 number of promoters and 256 number of lac repressors (inset of Fig. 4 in Ref. [24]. When and then one finds that and the integral solution to Eqs. 1 can be implicitly written as follows.

Fig. 4

Statistical properties of the response times of multi NARLs motif. Numerical integration of Eqs. 1 was done with = 10−8 as in Eqs. 9. Common settings are w = 4.86 × 10−2, v = 1.3 × 10−3 and σ = 101.2. Here , w 10−3 and 10−3 and statistical properties were computed over 105 trajectories. A-C. Here z = 1. Clearly, increasing n and will increase the variance of the response times. There exists an optimum n at which the variance attains minimum. D-F. Here  = 8 × 10−4. G-I. Here . Clearly, increase in z will decrease the variance of the response times. Increase in will increase the variance of the response times and there exists an optimum Hill coefficient at which the variance of the response time attains a minimum.

Here we are dividing by ln (2) to convert the timescale into the number of generation times (gt) since by definition one generation time of the cell = ln (2) in the dimensionless space for a stable unregulated protein product and is the time required by the gene expression machinery to generate a fraction of the steady state protein concentration. Here r = ½ corresponds to the standard response time. When then the integrals in Eqs. 5a-b can be approximately inverted and one can express P as a function of as follows.

In line with Ref. [24], the steady state fold change corresponding to the NARL motif can be defined as = the steady state protein levels in the presence of auto regulation / steady state protein levels in the absence of auto regulation. With this definition when then we find the fold change since the total steady state protein level of z number of gene copies in the absence of auto regulation is . This means that when n = 1 then the fold change of NARL motif decreases with the gene copy number in a square root manner as fold-change . This decrease in fold-change with respect to an increase in z is similar to the titration effects of non-autoregulatory multi copy genes observed in Ref. [24] where the fold change increases along the gene copy number contrasting from the NARL motif. Using Eqs. 5 and 6 one can show that the response time of the multi NARLs will not be affected by the copy number z under strong binding conditions as follows.

When n = 1, then one recovers the known result for z = 1 as gt. Eqs. 6 and 7 clearly suggest that when and then the steady state protein levels of multi NARLs can be increased by raising z as without much change in their response times, especially when . Using detailed numerical simulation of Eqs. 1 we will show for finite and physiologically relevant values of that an increase in z will actually decreases the response times in addition to increasing the steady state protein levels under certain conditions.

Eukaryotic multi autoregulatory motifs

Unlike bacterial systems where the transcription and translation take place in the cytoplasm, in case of eukaryotic systems, the transcription occurs inside the nucleus whereas the translation occurs in the cytoplasm. Hence produced protein product will be generally in the inactive form which needs to be converted to the active form before it can bind the promoter. In this background, to understand the robust oscillatory behaviour of multi NARL motif, we considered the following set of coupled differential rate equations:

In this equation, P represents the normalized overall concentration of the inactive form of the protein product, where u is the concentration of the active protein product of the i NARL gene copy which can bind its own promoter and p is the concentration of the corresponding inactive form of the protein product. Since the system of Eqs. 1 can only produce asymptotic spirals (Murugan, 2014), we increased the nonlinearity of the dynamical system by adding the active-inactive conversion step of the protein product along with raising the Hill coefficient. One can also achieve this via introducing delay differential equation or Michaelis-Menten type enzyme mediated mRNA and protein recycling modes which are common in eukaryotic systems. Various dimensionless parameters are defined as follows.here (1/second) are the forward and reverse rates associated with the inactive-active form conversion of the protein product and (1/second) is the decay rate associated with the active protein product. In particular, (v, w, ε) are singular perturbation parameters by definition since they multiply the derivative terms and (κ, μ, λ, χ) are all ordinary perturbation parameters. The Jacobian matrix associated with the linearization of the system of nonlinear ordinary differential Eqs. 8 corresponding to the dynamical variables (X, M, P, U) can be derived as follows (Murugan, 2014):

This Jacobian matrix needs to evaluated at the steady state values of (X, M, P, U)S which can be computed from the following steady state expressions.

Since the dynamical variables (X, M, P, U) are finite, the system of Eqs. 8 can produce self-sustained oscillations only when it deviates from the Routh-Hurwitz criteria for the stability that also warrants at least one complex eigen value of the Jacobian matrix given in Eq. 10 with positive real component evaluated at the steady state. Earlier studies have shown (Murugan, 2014) that there exists a critical Hill coefficient for any given set of parameters (v, w, ε, κ, σ, μ, λ, χ) at which the system of Eqs. 8 can produce self-sustained oscillations. This means that a critical level of nonlinearity is required to produce oscillations and by raising the Hill coefficient one can increase the nonlinearity of the system of Eqs. 8. However, one should note that such obtained critical Hill coefficient will be meaningful only when it is close to the physiological relevant range of values. For example, typical value of Hill coefficient for the lac repressor system will be n = 4.

Dual feedback motifs with multi NARL and PARL components

The dual feedback motifs involve two different genes which autoregulate themselves as well as cross regulate each other as shown in Fig. 6. We denote the configuration of a dual feedback motif by “C1C2C3C4” where C1 and C3 represents the autoregulation types (+1 or −1 corresponding to positive or negative) of gene 1 and gene 2 respectively. The regulation of gene 1 by the protein product of gene 2 is denoted by C2. Whereas, C4 represents the type of regulation of gene 2 by the protein product of gene 1. By this definition, the dual feedback oscillator motif constructed in Ref. (Stricker et al., 2008) will be denoted as “-1 + 1 + 1–1” configuration where gene 1 is the repressor and gene 2 is the activator.

Fig. 6

Dual feedback oscillator motif. Here the repressor gene (R) negative regulates its own expression apart from repressing the activator gene. Whereas, the activator gene (A) positive regulates its own expression apart from enhancing the expression of the repressor gene. Here m, p, and u represent the concentrations of mRNA, inactive and active form of protein product of gene H = (A, R), k and k are the transcription and translation rate of gene H = (A, R), k is the forward rate constant associated with the binding of the active form of protein product of gene T = (A, R) with the promoter of gene S = (A, R) and k is the respective dissociation rate constant, γ, γ and γ are the recycling rates corresponding to mRNA, active and inactive forms of protein product of gene H = (A, R), λ and λ are the forward and reverse rates associated with the conversion of inactive to active protein product of gene H = (A, R). One can denote a given dual feedback motif as, “C1C2C3C4” where C1-4 denote the signs of regulation. In this nomenclature, C1 and C2 correspond to gene R and C3 and C4 correspond to gene A. For example, the configuration “-1 + 1 + 1–1” is the dual feedback oscillator considered in Ref. (Stricker et al., 2008). Here, C1, C3 are the type of autoregulation of the genes R and A respectively. C2 is the type of cross regulation of gene R by the protein product of gene A and C4 is the type of cross regulation of gene A by the protein product of R. The motifs defined by “+1–1 + 1–1”, “+1–1 + 1 + 1” and “-1–1 + 1 + 1” are bistable decision making motifs. Depending on the parameter settings, these motifs robustly evolve into a state where one gene is turned on and the other one is tuned off.

To understand the effect of gene copy number on the oscillatory behaviour of dual feedback type motifs, we considered the synthetic one constructed in Ref. (Stricker et al., 2008). Here there are two different transcription factor genes viz. activator (A, ) and repressor (R, ) cross regulate each other. The activator gene (A) positively regulates its own transcription apart from upregulating the repressor gene. Whereas, the repressor gene (R) negatively self-regulates its own expression apart from downregulating the activator gene. The estimated copy number of activator gene seems to be z = 50 and for the repressor it is z = 25 and their binding Hill coefficients are n = 2 and n = 4 (Stricker et al., 2008). Therefore, the dual feedback motif considered in Ref. (Stricker et al., 2008) is actually a coupled one with multi NARL and PARL components. Remarkably, both the activator and repressor genes are controlled by both their protein products via an OR type logic. We denote the fraction of the promoter of the gene H occupied by the active protein product of gene G by X. Clearly, the total fraction of the promoter of gene H bound with the protein products of H and G will be X = X + X where the subscripts can be (H, G) = (A, R). We denote the normalized concentrations of mRNA, inactive and active forms of protein product of gene H respectively by M, P and U. With these settings, the differential rate equations associated with the dual feedback oscillatory motif with configuration (-1 + 1 + 1–1) can be written in the dimensionless form as follows (see Appendix B and Table 1 for details).

Table 1

Parameters used for the simulation of dual feedback oscillator.

Parameter	Definition and settings
XGH	Fraction of the promoter of gene G bound with the protein product of gene H. (G, H) = (R, A). R = repressor, A = activator.
XH	XR = XRR + XRA, XA = XAA + XAR
MH	Scaled mRNA level.
PH	Scaled level of inactive form of protein.
UH	Scaled level of active form of protein.
nH	nA = 2, nR = 4, binding Hill coefficients.
zH	zA = 50, zR = 25, gene copy numbers.
vH	vA = 0.02, vR = 0.05 (we assume that v_HQ=v_H for H, Q = A, R)
wH	wA = 0.1, wR = 15.1
ƐH	ƐA = 1.9, ƐR = 12.1
χH	χA = 0.9, χR = 2.9 (we assume that χ_HQ=χ_Q for H, Q = A, R)
λH	λA = 1.1, λR = 1
μH	μA = 0.1, μR = 0.004 (we assume that μ_HQ=μ_Q for H, Q = A, R)
κH	κA = 0, κR = 0

Dimensionless form of rate equations for the activator (A) gene regulation

Dimensionless form of rate equations for the repressor (R) gene regulation

In Eqs. 12–13, and are the input functions associated with the transcription control of the dual feedback oscillator motif and various dimensionless parameters are defined as follows.

In Eqs. 12–13, we have defined various dynamical variables as follows.here (m, p, u) represent the concentrations of mRNA, inactive and active form of the protein product. From Eqs. 15 one finds that and the steady state concentrations of various dynamical variables are defined as follows.here m, p, and u represent respectively the concentrations of mRNA, inactive and active form of protein product of gene H, k and k are the transcription and translation rate of gene H, k is the forward rate constant associated with the binding of the active form of protein product of gene T with the promoter of gene W and k is the respective dissociation rate constant, γ, γ and γ are the recycling rates corresponding to mRNA, active and inactive forms of protein product of gene H, λ and λ are the forward and reverse rate constants associated with the conversion of inactive protein product to active protein product of gene H = (A, R). The copy number of gene H is z where H = (A, R). Since the binding sequences are similar for the same protein product, one can assume that for H, Q = A, R which greatly simplifies Eqs. 12 and 13. We further simplify by assuming that and for H, Q = A, R. The input functions for the decision making motif configuration “+1–1 + 1–1” will be, and where both the genes R and A are positive auto regulated and a they negative cross regulate each other as shown in Eqs. B9 of Appendix B. Various parameters used for numerical simulation of decision making motifs are given in Table 2.

Table 2

Parameters used for the simulation of bistable decision making motifs.

Parameter	Definition and settings (Motifs “+1–1 + 1–1”, “-1–1 + 1 + 1”, “+1–1 + 1 + 1”)
nH	nA = 4, nR = 4, binding Hill coefficients.
zH	zA = 8, zR = 5, gene copy numbers.
vH	vA = 0.001, vR = 0.01 we assume that v_HQ=v_H for H, Q = A, R)
wH	wA = 5.5, wR = 5.5
ƐH	ƐA = 2, ƐR = 2
χH	χA = 50, χR = 50 (we assume that χ_HQ=χ_Q for H, Q = A, R)
λH	λA = 2, λR = 2
μH	μA = 0.001, μR = 0.001 (we assume that μ_HQ=μ_Q for H, Q = A, R)
κH	κA = 0, κR = 0

Methods

To obtain the physiologically relevant values of , we considered the digitized dataset from Ref. (Rosenfeld et al., 2002). Upon close observation over this dataset, one finds that gt and gt which were earlier (Murugan and Kreiman, 2011) used as the exit conditions for the iterative numerical integration algorithm to obtain the parameters of Eqs. 1. The earlier obtained (Murugan and Kreiman, 2011) theoretical estimates of the parameters satisfying these criteria were w = 0.1, v = 5 × 10−4, µ = 7 × 10−3 and σ = 4. Here, the corresponding parameters were obtained via nonlinear least square (NLS) fitting of this dataset with Eqs. 1 using Marquart-Levenberg algorithm (Marquardt, 1963) with the fixed settings n = 1 and z = 1. Hence obtained parameter values are w = 4.86 × 10−2; v = 1.3 × 10−3, µ = 8 × 10−4 and σ = 101.2. The NLS fit along with the digitized data with and without fixing n and z values are shown in Fig. 2.

Fig. 2

A. Nonlinear least square fitting of the experimental data from Ref. (Rosenfeld et al., 2002) to Eqs. 1 with fixed n = 1 and z = 1 using Marquardt-Levenberg algorithm. The obtained parameters were w = 4.86 × 10−2, v = 1.3 × 10−3, µ = 8 × 10−4 and σ = 101.2 where the standard errors of these fit values are < 10−5 with root mean square error (RMSE) = 8.1 × 10−4 and R2 = 0.9989. Since the fitting results are dependent on the sample size, after digitizing the experimental data we filled the gaps between the consecutive data points using interpolation method so that the final sample interval was Δτ = 10−4 where the data points are within and is measured in terms of number of generation times (gt). Since the experimental trajectory started at an offset of P/PS 0.08 at τ = 0, we subtracted this value from the digitized data before entering it into the fitting routine. B. Nonlinear least square fitting with fixed z = 1 and varying n. The fit parameters here are n = 4, w = 1.7 × 10−3, v = 8 × 10−4, µ = 5.1 × 10−3 and σ = 142.4. In this case, there is a significant deviation of the prediction from the observed dataset especially at the initial times. C. both n and z were varied. n = 4, z = 2, w = 1.8 × 10−4, v = 5 × 10−4, µ = 4.7 × 10−4 and σ = 186.6. D. n was fixed at 1 and z was varied, z = 4, w = 9.8 × 10−2, v = 4 × 10−2, µ = 8.7 × 10−3 and σ = 9.5.

We use the following Euler iterative scheme to numerically integrate Eqs. 1 with a given set of parameters along with the initial conditions at . We set appropriate time step to capture the promoter state dynamics along with the dynamics of mRNA and protein synthesis and degradation.

The steady state protein concentration for an arbitrary value of the Hill coefficient n and copy number z was computed by numerically solving for PS. The response time of multi NARLs motif was computed from the integral trajectory of Eqs. 17 via finding the point of time at which the cumulative expression machinery generates rP amount of protein product starting from zero protein level. To understand the stochastic nature of the response times defined by Eqs. 1, one can consider the following set of chemical Langevin equations 1 ([18], [30]).

In these equations, we have defined the parameters , and and various noise terms corresponding to the dynamical variables (X, M, P) in terms of individual genes and reaction steps of multi NARL module can be approximated (see Appendix A for details) as follows.

The mean and the variance associated with the various delta correlated gaussian white noise terms are defined as follows.

To stochastically simulate Eqs. 18–20, we use the following Euler iterative scheme.

In this equation, q1…q6 are random numbers drawn from the standard normal distribution and the initial conditions are set as (X, M, P) = (0, 0, 0). To obtain the critical Hill coefficient associated with the oscillatory behaviour of Eqs. 9, the eigenvalues of the corresponding Jacobian matrix given in Eq. 10 was evaluated at the steady state given by Eqs. 11. Presence of complex roots with positive real part suggests the possibility of self-sustained oscillations which will be further checked by numerical integration of Eqs. 9. The numerical workflow corresponding to the fixed z = 1 and z = 8 is as follows.

start from the physiologically relevant values of the parameters (v, w, ε, κ, μ, λ, χ, n) = (5 × 10−5, 1.1, 1, 0, 10−5, 1,10, 8).

iterate w inside (0.1, 2), v inside (10−6, 5 × 10−4) and μ inside (10−6, 10−5) by fixing other parameters at constant values.

at each iteration, the Jacobian matrix will be evaluated at the steady state that is computed using Eqs. 11 and the corresponding characteristic polynomial and eigenvalues will be computed.

upon finding complex eigenvalues with positive real part, iteration will be stopped.

the obtained set of parameters will be checked for the presence of self-sustained oscillations by numerical integration of the system of Eqs. 10 using Euler scheme as in Eqs. 19 with these parameters and the corresponding amplitude and period will be evaluated.

similar methodology can be used to obtain the critical Hill coefficient that is required to generate self-sustained oscillations. Here the Hill coefficient will be iterated starting from n = 1 by fixing all the other parameters at constant values.

Perturbation of parameters (v, w) does not change the steady state levels of various dynamical variables, which is evident from Eqs. 10. However, changes in μ as well as z will change the steady state values of the dynamical variables. Therefore, upon computing the trajectories inside the perturbation time window, one needs to consider the change in the steady state values of the dynamical variables (X, M, P, U) before the normalization step. Similar computational workflow was used to numerically integrate the dynamical rate equations corresponding to the dual feedback oscillatory motif described by Eqs. 12–13. Numerical simulation of dual feedback decision making motifs will be similar to Eqs. 12–13 with modified input functions as given in Eqs. B9 of Appendix B.

Results

Physiologically relevant parameter values

For a given gene of multi NARLs motif, there is d = 1 promoter copy, h 103 molecules and g 102 molecules. In the τ time space, we find the mRNA recycling rate γ 1/w and the protein recycling rate γ 1. We assume that under in vivo conditions, the protein interacts with its own promoter via three-dimensional diffusion with a rate k 106 (mol/lit)−1 s−1. The concentration of a single specific binding site or a protein molecule inside the E. Coli cell with volume 10−15 lit will be 2 × 10−9 mol/lit and the number of collisions that can happen between a single transcription factor protein product with its own promoter will be on the order of k 10−3 molecules−1 s−1 (Murugan, 2007). Here we have used the scaling, 1 molecule 2 nM inside the cellular volume. For a protein lifetime of 60 min, we find γ 3 × 10−4 s−1. This means that in the τ space k (10−3/γp) molecules−1. Using these values, one finally obtains , and . Upon nonlinear least square fitting of Eqs. 1 with fixed (n, z) = 1 over the digitized and normalized experimental data on the negatively auto-regulated E.coli TetR system (Rosenfeld et al., 2002) using Marquardt-Levenberg algorithm as in Fig. 2A, one finds that (Methods). We used these parameter fit values as the standards of individual NARL motif to evaluate various dynamical properties of the multi NARLs motif in the context of speeding up the response times. Remarkably, when n and z are allowed to vary along the nonlinear regression fit process of the same dataset, then one obtains n = 4 and z = 2 as the best fit values along with slight changes over the earlier best fit parameters (w, v, μ, σ) as demonstrated in Fig. 2B-D.

Variation of the response time and the steady state protein levels with respect to changes in z, n and are demonstrated in Fig. 3. Though an increase in increases the steady state protein level in a monotonic manner (Fig. 3B, 3D and 3F), variation of these parameters influences the response times in a complicated manner (Fig. 3A, 3C, 3E). Remarkably, increase in z increases the steady state protein level without much change in the response times (Fig. 3A-B) under strong binding conditions especially when as shown in Fig. 3A and B. Under such conditions, there exists an optimum Hill coefficient at which the response time attains a minimum that is demonstrated in Fig. 3C. Increase in with fixed μ monotonically increases the steady state protein levels as shown in Fig. 3D. From Fig. 3C one can conclude that when z = 1, then . For higher values of z and under weak binding conditions as , the optimum Hill coefficient asymptotically shifts towards . Interestingly, increasing z can decrease the response time when n < n. However, when n > n, then increasing the copy number increases the response time as shown in Fig. 3C. Fig. 3E suggests that for z = 1, the optimum n occurs only when .

Fig. 3

Dynamics of multi NARL loops. Numerical integration of Eqs. 1 was done with = 10−8 as in Eqs. 8. Common simulation settings are w = 4.86 × 10−2, v = 1.3 × 10−3 and σ = 101.2. A-B. Here n = 1. Clearly, increase in will increase both the steady state protein level P as well as the response time measured in number of generation times (1 gt = ln(2) in the τ space). Increase in z will increase the steady state protein level. The response time is slightly dependent on z only under strong binding condition especially when µ < 10−2. C-D. Here µ = 8 × 10−4. There exists an optimum n at which the response time attains a minimum. When z = 1, then the optimum n = 2. When z > 10, then the optimum Hill coefficient approach n = 4. E-F. Here z = 1 and the optimum Hill coefficient seems to be dependent on µ under strong binding conditions.

The computed statistical properties of the response time of the multi NARLs module are shown in Fig. 4. Increase in with fixed z = 1 increases the variance of the response times in a turn over manner with an optimum Hill coefficient at which the variance attains a minimum (Fig. 4A). Remarkably, when , then there exists an optimum Hill coefficient which shifts to under weak binding conditions as at which the variance and the Fano factor associated with the response times attains a minimum (Fig. 4A, 4B, 4C). This optimum point of n disappears as z increases (Fig. 4D) and at high z values, the variance decreases monotonically with increase in n. Increasing z monotonically decreases the variance, coefficient of variation and the Fano factor of the response times especially when as shown in Fig. 4D, 4E and 4F. Furthermore, when n = 1, then there exists an optimum at which the variance, CV and Fano factor of the distribution of the response times attain local minima (Fig. 4G, 4H, 4I).

Motifs with multi NARL components can be robust oscillators

The oscillatory dynamics of the single and multi NARLs motif are demonstrated in Fig. 5A-B. Here the set of parameters which can produce self-sustained oscillations were obtained by repeated iterative workflow as given in the simulation methods section. To understand the robustness of oscillations, perturbations in the control parameters (w, v, µ) were introduced as temporal rectangular pulses. Change in these parameters, will eventually change or increase the required critical Hill coefficient. As a result, the oscillatory orbit of the perturbed system breaks down in to asymptotic spirals as demonstrated in Fig. 5C–H. When the oscillatory motif is robust against perturbations, then it should return back to either the original unperturbed orbit or some other self-sustained orbit upon removal of the perturbations. Clearly, the single NARL motif is not robust against fluctuations, especially in the singular perturbation parameters (v, w). The multi NARLs motif is robust against perturbations in all the parameters (w, v, µ). Results suggest that one can increase the robustness of the genetic oscillator by raising the copy number z as demonstrated in Fig. 5D, F, H where z = 8. Remarkably, the single NARL motif is still robust enough against the perturbations in µ as shown in Fig. 5G. Here perturbations in µ can be well correlated with the fluctuations in the promoter state binding-unbinding dynamics.

Fig. 5

Multi NARL motif as a robust genetic oscillator. Common settings in A, C, E, G are z = 1, w = 1.1, v = 5 × 10−5, κ = 0, χ = 10, ε = 1, λ = 1, µ = 10−5 which corresponds to the critical n = 8. Here period is around 10.19 and the steady state protein level is around 0.27. Common settings in B, D, F, H are z = 8, w = 1.6, v = 4 × 10−5, κ = 0, σ = 1, χ = 10, ε = 1, λ = 1, µ = 10−6 which corresponds to the critical n = 7. Here period is around 10.38 and the steady state protein level is around 0.23. Perturbations in various parameters start at τ = 50 and ends at τ = 100. A, B. unperturbed trajectories. C, D. perturbation is introduced in w. Upon perturbation, the critical Hill coefficient corresponding to z = 1 increases from 8 to 12. For z = 8, from 7 it increases to 10. E, F. perturbation is introduced in v. Upon perturbation, the critical Hill coefficient corresponding to z = 1 increases from 8 to 785. For z = 8, from 7 it increases to 11. G, H. perturbation is introduced in µ. Upon perturbation, the critical Hill coefficient corresponding to z = 1 increases from 8 to 14. For z = 8, from 7 it increases to 9. Single NARL motif is more sensitive to the perturbations in (w, v) than the multi NARL motif since the oscillatory behaviour is lost upon removal of perturbation.

Dual feedback oscillator and decision-making motifs with multi NARLs and PARLs

The effects of gene copy number on the oscillatory behaviour of the dual feedback motif constructed in Ref. (Stricker et al., 2008) (Fig. 6) is demonstrated in Fig. 7, Fig. 8 with the set of parameter values given in Table 1. Here the activator ara C gene positive self-regulates its own transcription apart from upregulating the repressor lac I gene. Whereas, the repressor gene negative self-regulates its own transcription apart from downregulating the ara C gene. The promoters of both the genes are regulated by the protein products of both the genes via an OR type logic. Experimental estimation of the gene copy numbers are z = 50 and z = 25 with the corresponding Hill coefficients n = 2 and n = 4 (Stricker et al., 2008). Fig. 7B-D and 8B-D clearly suggest that 1) increase in the copy numbers z and z of the dual feedback motif described in Fig. 6 increases the robustness of oscillation against temporal perturbations in various critical parameters associated with repressor as well as activator genes and, 2) the period and amplitude of oscillations can also be tuned by varying the gene copy numbers apart from the inducer concentration which is evident from Fig. 7A-D and 8A-D. 3) The motif with z = 1 and z = 1 more sensitive to the temporal perturbations in the binding parameters µ and µ than the coupled multi NARL and PARL one which is evident from Fig. 3, Fig. 3. Robustness behaviour of the decision-making dual feedback motif in the presence of multiple auto regulatory components is demonstrated in Figs. A-D.

Fig. 7

Variation of the oscillatory behaviour of the dual feedback motif with respect to changes in the parameters of repressor gene R. Here the unperturbed parameter settings are n = 2, n = 4, v = 0.04, v = 0.05, w = 0.1, w = 15.1, Ɛ = 1.9, Ɛ = 12.1, χ = 0.9, χ = 2.9, λ = 0.4, λ = 1, μ = 0.1, μ = 0.004, κ = 0, κ = 0. The subscript ‘A’ denotes the activator and ‘R’ represents the repressor gene. Here we assumed for H, Q = A, R to simplify Eqs. 12 and 13. Index 1 represents the temporal perturbation in the respective parameter, 2 represents the oscillatory patter corresponding to zA = 1 and zR = 1 and 3 represents z = 50 and z = 25. A1-3. time dependent perturbations introduced in μR. B1-3. time dependent perturbation introduced in the parameter w of the repressor gene. C1-3. time dependent perturbation introduced in the parameter vR of the repressor gene. D1-3. time dependent perturbation introduced in the parameter εR of the repressor gene.

Fig. 8

Variation of the oscillatory behaviour of the dual feedback motif with respect to changes in the parameters of the activator gene A. Here the unperturbed parameter settings are n = 2, n = 4, v = 0.04, v = 0.05, w = 0.1, w = 15.1, Ɛ = 1.9, Ɛ = 12.1, χ = 0.9, χ = 2.9, λ = 0.4, λ = 1, μ = 0.1, μ = 0.004, κ = 0, κ = 0. The subscript ‘A’ denotes the activator and ‘R’ represents the repressor gene. Here we assumed for H, Q = A, R to simplify Eqs. 12 and 13. Index 1 represents the temporal perturbation in the respective parameter, 2 represents the oscillatory patter corresponding to zA = 1 and zR = 1 and 3 represents zA = 50 and zR = 25. A1-3. time dependent perturbations introduced in μA. B1-3. time dependent perturbation introduced in the parameter wA of the repressor gene. C1-3. time dependent perturbation introduced in the parameter vA of the repressor gene. D1-3. time dependent perturbation introduced in the parameter εA of the repressor gene.

The dynamical behaviour of the decision-making dual feedback motifs ([32], [33]) are demonstrated in Fig. 9, Fig. 10, Fig. 11. The genes R and A embedded in these motifs are such that depending on the parameter settings, one gene will be turned on and the other will be turned off. However, the stability of such decision state can be easily perturbed by fluctuations in the critical parameters as shown in Fig. 9, Fig. 10, Fig. 11A2, B2, C2 and D2. Clearly, robustness in the stability of the decision state against perturbations in the parameters (ε, w, µ) can be increased by rising the copy number of the positive as well as negative autoregulatory loops as demonstrated in Fig. 9, Fig. 10, Fig. 11A3, B3, and C3.

Fig. 9

Variation of the decision making behaviour of the dual feedback motif configuration “+1–1 + 1–1” (Appendix B) with respect to changes in the parameters of the gene A which is also applicable to gene R since there is a degeneracy in the motif architecture (R and A and interchangeable). Here the unperturbed parameter settings are n = 4, n = 4, v = 0.001, v = 0.01, w = 5.5, w = 5.5, Ɛ = 2, Ɛ = 2, χ = 50, χ = 50, λ = 2, λ = 2, μ = 0.001, μ = 0.001, κ = 0, κ = 0. The subscript ‘A’ denotes the activator and ‘R’ represents the repressor gene. Here we assumed for H, Q = A, R to simplify Eqs. 12 and 13. Index 1 represents the temporal perturbation in the respective parameter, 2 represents the oscillatory patter corresponding to zA = 1 and zR = 1 and 3 represents zA = 8 and zR = 5. A1-3. time dependent perturbations introduced in wA. B1-3. time dependent perturbation introduced in the parameter μA of the A gene. C1-3. time dependent perturbation introduced in the parameter vA of the A gene. D1-3. time dependent perturbation introduced in the parameter εA of the A gene.

Fig. 10

Variation of the decision making behaviour of the dual feedback motif configuration “+1–1 + 1 + 1” (Appendix B) with respect to changes in the parameters of the gene A which is also applicable to gene R since there is a degeneracy in the motif architecture (R and A and interchangeable). Here the unperturbed parameter settings are n = 4, n = 4, v = 0.001, v = 0.01, w = 5.5, w = 5.5, Ɛ = 2, Ɛ = 2, χ = 50, χ = 50, λ = 2, λ = 2, μ = 0.001, μ = 0.001, κ = 0, κ = 0. The subscript ‘A’ denotes the activator and ‘R’ represents the repressor gene. Here we assumed for H, Q = A, R to simplify Eqs. 12 and 13. Index 1 represents the temporal perturbation in the respective parameter, 2 represents the oscillatory patter corresponding to zA = 1 and zR = 1 and 3 represents zA = 8 and zR = 5. A1-3. time dependent perturbations introduced in wA. B1-3. time dependent perturbation introduced in the parameter μA of A gene. C1-3. time dependent perturbation introduced in the parameter vA of A gene. D1-3. time dependent perturbation introduced in the parameter εA of A gene.

Fig. 11

Variation of the decision making behaviour of the dual feedback motif configuration “-1–1 + 1 + 1” (Appendix B) with respect to changes in the parameters of the gene A which is also applicable to gene R since there is a degeneracy in the motif architecture (R and A and interchangeable). Here the unperturbed parameter settings are n = 4, n = 4, v = 0.001, v = 0.01, w = 5.5, w = 5.5, Ɛ = 2, Ɛ = 2, χ = 50, χ = 50, λ = 2, λ = 2, μ = 0.001, μ = 0.001, κ = 0, κ = 0. The subscript ‘A’ denotes the activator and ‘R’ represents the repressor gene. Here we assumed for H, Q = A, R to simplify Eqs. 12 and 13. Index 1 represents the temporal perturbation in the respective parameter, 2 represents the oscillatory patter corresponding to zA = 1 and zR = 1 and 3 represents zA = 8 and zR = 5. A1-3. time dependent perturbations introduced in w. B1-3. time dependent perturbation introduced in the parameter μR of R gene. C1-3. time dependent perturbation introduced in the parameter vR of R gene. D1-3. time dependent perturbation introduced in the parameter εR of R gene.

Discussion

Response time decides how fast a gene can respond to an external adverse signal at the transcription level in a signalling cascade. The external signal can be physical, such as light, or chemical, such as small ligand molecules. Apart from the response time, the steady state protein level of the responding gene decides the coupling between two consecutive members of a signalling cascade. Though the negative autoregulatory loop can speed up the response time of a gene, it drastically reduces the steady state protein level, which in turn weakens the subsequent signalling processes. Therefore, it is important for a cell to have a network motif of signalling cascades which can speed up the response time without compromising on the steady state protein levels. The steady state protein level of an individual NARL motif can be increased by reducing the affinity of its protein product towards its promoter. This can be achieved by tuning the specific binding sequence of the protein product such that it increases the binding / unbinding dynamics of the protein with its promoter (parameter ). Instead, one can increase the Hill coefficient n to raise the steady state protein levels, which is however not straightforward since one needs to design the protein–protein interactions apart from tuning the specific binding DNA sequence. Here increase in the Hill coefficient will in turn increases the response time since the time required to form the oligomer introduces additional time delay in the negative feedback.

In this context, the multi NARL motif presented here can be an efficient design since we are able to tune both steady state protein level as well as response time in the required direction using a single parameter z under weak or strong binding conditions. Remarkably, under strong binding conditions, i.e., as , it is possible to raise the steady state protein levels with almost no change in the response time compared to the single NARL motif. The overall steady state protein level increases mainly because of the contribution of each member of the multi NARL motif. Under strong binding conditions, as , the total protein concentration will be shared among z promoters towards the saturation level in an equal manner and therefore the availability of protein molecules per promoter does not change much. As a result, the response time of the multi NARL motif will be the same as that of the single NARL and it does not change much with respect to z under strong binding conditions. In contrast, under weak binding conditions, as , there will be plenty of free protein molecules available for binding, which increases with z. As a result, increasing z drives the promoters-protein concentration towards its saturation. This will ultimately decrease the response time.

Detailed numerical analysis suggests that one can increase the robustness of the single NARL type genetic oscillator by increasing the copy number z. One of the main drawbacks of the synthetic gene oscillators such as the one that was constructed by Striker et. al. (Stricker et al., 2008) is the disappearance of oscillations after few generations ([21], [22]) of the cell compared to the negative feedback only motif. This observation could be well attributed to the extrinsic type fluctuations in the control parameters or the chaotic dynamics of the OR type logic gates involved in the coupling of dual feedback motifs (Murugan, 2014). Interestingly, the estimated copy numbers of the their dual feedback motif were approximately 25 and 50 for the repressor and activator plasmids respectively (Stricker et al., 2008). This means that the negative feedback only motif constructed in Ref. (Stricker et al., 2008) is actually a multi NARL motif with z = 25 which could be the reason why their negative feedback only motif exhibited stable oscillations over several generations compared to the dual feedback motif which additionally suffered from the chaotic dynamics arising out of the OR type logic control.

To avoid such chaos issues, one needs to construct a complicated dual feedback motif with AND type logic (Murugan, 2014) which in turn involves tedious design engineering of protein–protein interactions. In this context, the multi NARL motifs seems to be stable against the fluctuations in the singular perturbation parameters (v, w). Here perturbation in v reflects the promoter state fluctuations and perturbation in w represents the fluctuations at the transcription and translation levels. Perturbation in µ can be correlated with the fluctuations in the inducer concentration which in turn affects the binding-unbinding dynamics of the active protein-product at the promoter. Multi NARLs motif achieve this robustness via the averaging and collective entrainment effects. When there are several identical oscillatory motifs functioning simultaneously in an asynchronous mode, those motifs with broken orbits due to perturbations in the control parameters will be pulled back to the original orbit by the unaffected motifs via the averaging and entrainment effect.

The titration effects of multi copy gene expression systems have been studied under steady state conditions earlier in detail ([24], [25], [26], [27]). In Ref. (Brewster et al., 2014), the fold-change enhancement effects of multiple copies of the lac repressor controlled YFP reporter gene have been demonstrated. Particularly, the increase in the fold-change with respect to the gene copy number will be prominent when the number of lac repressor molecules inside the cell is lesser than the number of gene copy promoters. However, when the number of repressor molecules are higher enough to saturate the promoters, then the fold-change will be similar to that of a single copy gene system (Brewster et al., 2014). Contrasting from these observations, in case of multi NARL motif, we have shown that the fold change deceases by increasing the gene copy number z as fold-change . When , then the fold-change scales with z as . When , then the fold change scales with z as . Here n can be correlated with the number of protein product molecules binding each promoter and mimics the binding saturation effect that was observed in Ref. (Brewster et al., 2014).

In the derivation of Eqs. 1 for NARL motif, we have assumed that n number of protein products, simultaneously binds their own promoters and downregulate their own production. Here n is the Hill coefficient (Hill, 1910). In case of lac repressor, a tetramer binds the operator sequence and blocks the transcription dynamics of RNA polymerase. Here one can consider at least three different possible mechanisms of assembly of four lac repressor molecules at the operator sequence viz. 1) simultaneous arrival of four lac repressor monomers at the operator location via a combination of 1D and 3D diffusion route. 2) lac repressor dimerizes first and then these dimers assemble at the operator sequence via 1D3D diffusion and 3) dimers form tetramers which will then arrive at the operator sequence via 1D3D diffusion. Since bulky molecules diffuse slowly over 1D as well as 3D routes, one can conclude that simultaneous arrival of the monomers at the location of the operator via 1D3D diffusion is the most probable pathway (Murugan, 2014) of site-specific binding of n TFs at their cognate sites. The oscillatory behaviour of Eqs. 8 depends on the nonlinearity of the system that is proportional to the Hill coefficient n. When two dimers or a tetramer bind the promoter, the effective n will be<4 which will affect the oscillatory behaviour of the system.

Conclusions

The response time is defined as the time required by a gene expression machinery to synthesize half of its steady state protein level. Response time decides the fastness of a gene to respond to an external signal at the transcription level in a signalling cascade. The steady state protein level of the responding gene decides the coupling between two consecutive members of a signalling cascade. In transcription autoregulation of genes, the protein product binds the promoter and up and down regulates its own synthesis.

The negative autoregulatory loop (NARL) can speed up the response time of a gene at the cost of reduced steady state protein level. Here we have designed a novel multi NARL motif which can be tuned for both steady state protein level as well as response time. We have shown that under strong binding condition the steady state protein levels of the multi NARLs motif can be increased with almost no change in the response time.

Detailed analysis revealed the existence of an optimum Hill coefficient at which the response time, its variance and its coefficient of variation attain minima. When the Hill coefficient is set below this optimum level, the steady state protein level of the multi NARL motif can be increased by increasing the gene copy number along with decrease in the response time. Instead, when the Hill coefficient was set above the optimum value, then an increase in the copy number of multi NARL motif would increase both the steady state protein levels and response time. These aspects clearly reveal the tunability of the multi NARL motifs towards both the response time and steady state protein levels. Using detailed computational analysis, we have demonstrated that the multi NARL motif can act as an oscillator which is robust against extrinsic fluctuations in the parameters those control the promoter state dynamics. We further demonstrated that the period of oscillation of the dual feedback oscillator can be fine-tuned by the gene copy number and the observed remarkable robustness of such oscillator mainly originate from the coupled multi NARL and PARL motifs. Finally we also demonstrate the robustness of decision making dual feedback motif embedded with multiple auto regulated components.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.