Tunability of the Dual Feedback Genetic Oscillator Modeled by the Asymmetry in Transcription and Translation

TL;DR

This paper presents a simplified delay differential equation model of a dual feedback genetic oscillator, revealing how asymmetries in gene expression parameters can tune and improve its robustness and tunability.

Contribution

The authors develop a reduced DDE model that captures oscillator tunability and predict how asymmetries in gene expression affect its dynamics, guiding experimental redesign.

Findings

Oscillator period is sensitive to DNA copy number asymmetry.

Translation rate asymmetry mimics effects of plasmid copy number differences.

Modulating mRNA degradation asymmetry enhances tunability and robustness.

Abstract

Oscillatory gene circuits are ubiquitous to biology and are involved in fundamental processes of cell cycle, circadian rhythms and developmental systems. The synthesis of small, non-natural oscillatory genetic circuits have been increasingly used to test fundamental principles of genetic network dynamics. A recently developed fast, tunable genetic oscillator by Stricker et al.[23] has demonstrated robustness and tunability of oscillatory behavior by combining positive and negative feedback loops. This oscillator combining lacI (negative) and araC (positive) feedback loops, was however modeled using multiple layers of differential equations to capture the molecular complexity of regulation, in order to explain the experimentally measured oscillations. We have developed a reduced model based on delay differential equations (DDEs) of this dual feedback loop oscillator, that reproduces the tunability of oscillator period and amplitude based on the concentration of the two inducers isopropyl b-D-1-thiogalactopyranoside (IPTG) and arabinose. Previous work had predicted a need for an asymmetry in copy numbers of activator (araC) and repressor (lacI) genes encoded on plasmids. We use our reduced model to redesign the network by comparing the effect of asymmetry in gene expression at the level of (a) DNA copy numbers and the rates of (b) mRNA translation and (c) degradation. We find the minimal period of the oscillator is sensitive to DNA copy number asymmetry, but translation rate asymmetry has an identical effect as plasmid copy numbers, while modulating the asymmetry in mRNA degradation can improve the tunability of period of the oscillator, together with increased robustness to replication 'noise' and influence of the host cell cycle. Thus, our model predicts experimentally testable principles to redesign a potentially more robust oscillatory genetic network.

Tables

Table 1. Table 1: The values of parameters used in the model are based on either previous reports [ 23 ] , or varied.

Parameter	Value	Description	Reference
b	$0.36$ min-1	Basal transcription rate	[23]
$\alpha$	$20$	Transcription activation	[23]
$n_a$	$50$ and varied	araC DNA copy number	[23], this study
$n_r$	$25$ and varied	lacI DNA copy number	” ”
$k_m^a,k_m^r$	$0.54$ min-1	Degradation rate of mRNA	” ”
${\sigma }_A$	$90$ min-1	Translation rate (RBS efficiency) of araC	[23],
${\sigma }_R$	$90$ min-1	Translation rate (RBS efficiency) of lacI	[23], this study
$k_f$	$0.9$ min-1	Rate of folding of proteins	” ”
${\gamma }_0$	$1080$ min-1	Maximal degradation	” ”

Table 2. Table 2: Rates of some equilibrium reactions are taken from Stricker et al. [ 23 ]

Parameter	Value	Description
$k_1^b$	1.8 $s^{-1}$	Backward rate of promoter binding binding by AraC/LacI
$k_d^a,k_d^r,k_t^r$	0.018 molecules${}_{}{}^{-1}\cdot s^{-1}$	Forward rates of dimerisation and tetramerisation
$k_{-d}^a,k_{-d}^r,k_{-t}$	0.00018 $s^{-1}$	Backward rates of dimerisation and tetramerisation

Equations

$${\dot{m}}_x=\frac{n_{x}b(1+\alpha k_1{A}_{{\tau }_1})}{(1+k_1{A}_{{\tau }_1})(1+k_2{R}_{{\tau }_2}{)}^2}-k_m^x{m}_x\text{,}$$

$$\text{where }x\in \{a,r\}$$

$$k_1=k_1^{min}+(k_1^{max}-k_1^{min})\frac{[ara{]}^2}{{\mu }^2+[ara{]}^2}.\frac{{\nu }^2}{{\nu }^2+[iptg{]}^2}$$

$$k_2=k_2^{min}+(k_2^{max}-k_2^{min})\frac{{\lambda }^2}{{\lambda }^2+[iptg{]}^2}$$

$${\dot{X}}_{uf}={\sigma }_X{m}_X-k_f{X}_{uf}-k_X{X}_{uf}\text{,}$$

$$\text{where }{X}_{uf}\in \{A_{uf},R_{uf}\}$$

$$k_R=\frac{{\gamma }_0}{c_0+\Sigma P}$$

$$k_A=g\cdot k_r$$

$$\dot{A}=\frac{1}{p_a}{k}_f{A}_{uf}-k_{a}A$$

$$\dot{R}=\frac{1}{p_r}{k}_f{R}_{uf}-k_{r}R$$

$$\frac{d{m}_X}{dt}$$

$$m_X=e^{-k_{m}t}\int n_{X}g(A_{{\tau }_1},R_{{\tau }_2})e^{k_{m}t}dt$$

$$\frac{d{X}_{uf}}{dt}=$$

$$y(t_0)-y(t_0-s)>0\text{and}y(t_0)-y(t_0+s)>0$$

$$P_{0,j}^{a/r}+A\stackrel{k_1}{\rightleftharpoons }{P}_{1,j}^{a/r}$$

$$P_{i,0}^{a/r}+R\stackrel{2k_2}{\rightleftharpoons }{P}_{i,1}^{a/r}$$

$$P_{i,1}^{a/r}+R\stackrel{\frac{1}{2}{k}_2}{\rightleftharpoons }{P}_{i,2}^{a/r}$$

$$P_{0,0}^{a/r}\stackrel{b_a}{\to }{P}_{0,0}^{a/r}+m_{a/r}$$

$$P_{1,0}^{a/r}\stackrel{\alpha b_a}{\to }{P}_{0,1}^{a/r}+m_{a/r}$$

$$m_a\stackrel{t_a}{\to }{m}_a+A_{uf}$$

$$m_r\stackrel{t_r}{\to }{m}_r+R_{uf}$$

$$A_{uf}\stackrel{k_{fa}}{\to }A$$

$$R_{uf}\stackrel{k_{fa}}{\to }R$$

$$m_a\stackrel{{\gamma }_{ma}}{\to }\varphi$$

$$m_r\stackrel{{\gamma }_{mr}}{\to }\varphi$$

$$A_{uf}\stackrel{k_a}{\to }\varphi$$

$$R_{uf}\stackrel{k_r}{\to }\varphi$$

$$A\stackrel{k_a}{\to }\varphi$$

$$R\stackrel{k_r}{\to }\varphi$$

$$\displaystyle\begin{split}P_{1,0}&=k_{1}A_{2}P_{0,0}\\ P_{0,1}&=2k_{2}R_{4}P_{0,0}\\ P_{0,2}&=\frac{k_{2}}{2}R_{4}P_{0,1}=\frac{k_{2}}{2}R_{4}\Big{(}2k_{2}R_{4}P_{0,0}\Big{)}=k^{2}_{2}R^{2}_{4}P_{0,0}\\ P_{1,1}&=\Big{(}k_{1}A_{2}\Big{)}P_{0,1}=2k_{1}A_{2}k_{2}R_{4}P_{0,0}\\ P_{1,2}&=\Big{(}k_{1}A_{2}\Big{)}P_{0,2}=k_{1}A_{2}k^{2}_{2}R^{2}_{4}P_{0,0}\end{split}$$

Full text

Tunability of the Dual Feedback Genetic Oscillator Modeled by the Asymmetry in Transcription and Translation

Yash Joshi1, Yash K. Jawale1 and Chaitanya A. Athale1,†,

1: Div. of Biology, IISER Pune, Dr. Homi Bhabha Road,

Pashan, Pune 411008, India.

$\dagger$ All correspondence to be addressed to: cathale@iiserpune.ac.in

Abstract

Oscillatory gene circuits are ubiquitous to biology and are involved in fundamental processes of cell cycle, circadian rhythms and developmental systems. The synthesis of small, non-natural oscillatory genetic circuits have been increasingly used to test fundamental principles of genetic network dynamics. A recently developed fast, tunable genetic oscillator by Stricker et al. [23] has demonstrated robustness and tunability of oscillatory behavior by combining positive and negative feedback loops. This oscillator combining lacI (negative) and araC (positive) feedback loops, was however modeled using multiple layers of differential equations to capture the molecular complexity of regulation, in order to explain the experimentally measured oscillations. We have developed a reduced model based on delay differential equations (DDEs) of this dual feedback loop oscillator, that reproduces the tunability of oscillator period and amplitude based on the concentration of the two inducers isopropyl $\beta$-D-1-thiogalactopyranoside (IPTG) and arabinose. Previous work had predicted a need for an asymmetry in copy numbers of activator (araC) and repressor (lacI) genes encoded on plasmids. We use our reduced model to redesign the network by comparing the effect of asymmetry in gene expression at the level of (a) DNA copy numbers and the rates of (b) mRNA translation and (c) degradation. We find the minimal period of the oscillator is sensitive to DNA copy number asymmetry, but translation rate asymmetry has an identical effect as plasmid copy numbers, while modulating the asymmetry in mRNA degradation can improve the tunability of period of the oscillator, together with increased robustness to replication ‘noise’ and influence of the host cell cycle. Thus, our model predicts experimentally testable principles to redesign a potentially more robust oscillatory genetic network.

1 Introduction

The ubiquity of oscillatory genetic networks suggests a central role in biology, resulting in extensive experimentally and theoretical studies as seen in case of the cell cycle clocks [15, 14, 17, 18, 5], circadian rhythms [9, 22] and developmental clocks in embryogenesis [8, 11, 25]. These oscillatory networks appear to have been selected for tunability and robustness as a part of the general ‘homeostatic’ mechanisms of such physiological processes critical to living systems. However, understanding the design principles of such networks raises challenges with due to their complexity. Increasingly, synthetic biology of small genetic networks has become an important alternative approach to gaining a fundamental understanding of fundamental principles of gene regulation driving such oscillators using small and relatively tractable genetic networks such as single gene negative-feedback systems [2], the three component ‘repressilator’ [4] and cell-free two- and three-stage gene cascades [16]. Based on theoretical studies of naturally occurring genetic oscillators, they have been broadly classified into either negative feedback loops or coupled positive and negative loops [5]. Indeed comparative modeling has demonstrated that while a minimal negative feedback loop network can produce oscillations, robustness and tunability are improved by the addition of other loops, which could explain the evolutionary selection of such complex networks [12].

A canonical example of a synthetic dual feedback loop - positive and negative- oscillator that has rapidly become a standard mode is the araC and lacI genetic oscillator with expression determined by a dual input $p_{lac/ara-1}$ promoter [23]. The genes are regulated by their own protein products (feedback) - activation by AraC protein in the presence of arabinose (positive loop) and repression by LacI protein in the absence of IPTG (negative loop). A model with twenty-seven coupled ordinary differential equations (ODEs) was developed to match the experimental findings [23], since a simple minimal model of dual feedback loops [6], failed to reproduce the oscillatory response of the system to parameter changes. The many intermediate reactions in the model such as the relatively slower rate of mRNA production, protein folding, protein multimerization and promoter binding play an important role in the experimental validation of this model. While explicit models of detailed molecular mechanisms are physically more realistic than minimal models, they also lead to an ‘explosion’ of the number of parameters and variables. One solution to develop a simple model that captures the essential nature of the process, while reproducing measurable dynamics of the system is through the use of delay differential equations [27]. However, it is not clear whether a delay differential equation model could be used to describe the dual-input oscillator.

A search of general design principles of genetic oscillators points to delays and noise as important features, alongside with topology as playing an important role [21]. The explicit use of delay differential equations for modeling genetic networks is seen in oscillator models of the cell cycle [5] and somitogenesis clock [11] and lac operon dynamics [27]. Indeed a cell-free extract based study of the lac-ara dual feedback oscillator have demonstrated that protein translation can serve as a bottleneck in the dynamics of the oscillator [16]. Hence, such a separation of time-scales seen in experiment and the utility of delays in oscillatory network models, suggests such an approach could help reduce the model complexity of the lac-ara oscillator.

Here, we describe a novel model of the lactose-arabinose dual feedback loop oscillator, that was first developed by Stricker et al. [23]. Our model consists of six delay differential equations, that take into account the canonical components of gene expression: DNA promoter states, RNA expression and protein translation and stability. Using two delay terms to represent intermediate states, we can quantitatively reproduce the experimentally reported inducer-dependent oscillatory behaviour, based on a DNA copy number asymmetry of the positive and negative feedback loops. Based on this model, we redesign of the network to examine whether we can improve the tunability of the oscillator, based on modifications at the level of (a) DNA copy numbers, (b) RNA degradation rates and (c) protein production rates.

2 Model

We have modeled the oscillator based on the coupled dynamics of (a) DNA transcription to RNA based on promoter state dynamics, (b) translation of RNA to protein and (c) protein folding (Fig. 1). The derivation and assumptions of this reduced model are described in detail in Appendix 5.2.

Transcription is modeled in terms of the dynamics of two mRNA species, encoding the AraC activator ($m_{a}$) and LacI repressor ($m_{r}$). The general equation of mRNA copies ($m_{x}$) for activator ($x=a$) and repressor ($x=r$) is:

[TABLE]

The two parts of the RHS of the equation represent the production and degradation terms. The feedback of activator dimers (A) and repressor tetramers (R) determine mRNA levels with the subscripts $\tau_{1}$ and $\tau_{2}$ indicating respective delays due to the formation of intermediate states. The terms $A_{\tau_{1}}\equiv A(t-\tau_{1})$ and $R_{\tau_{2}}\equiv R(t-\tau_{2})$, represent the equivalence of the number of protein dimers where $t$ is current time. The copy numbers of the genes encoding araC and lacI, $n_{a}$ and $n_{r}$ respectively, act as multiplication factors. We assume an asymmetry in plasmid copy numbers, i.e. $n_{a}\neq n_{r}$, based on the reported requirement of copy number differences of the genes in the original experiment by Stricker et al. [23]. Transcription rate is determined by $b$, the basal transcription rate when there is no activator bound to the promoter, $\alpha$ the multiplicative increase in transcription rate when activator (A) is bound and the degradation rates $k^{a}_{m}$ and $k^{a}_{m}$ of $m_{a}$ and $m_{r}$ respectively. The terms $k_{1}$ and $k_{2}$ are equilibrium binding rates of the A dimer and R tetramer binding to the promoter. These are modeled as variables based on a previously described expression [23]. The equilibrium binding rate of AraC binding to the promoter, $k_{1}$, depends on the concentrations of the inducers arabinose ([ara]) and IPTG ([iptg]) as:

[TABLE]

and the range of activator binding to the promoter is set by $k^{min}_{1}=0\text{ molecules}^{-1}$, $k^{max}_{1}=1\text{ molecules}^{-1}$ and $\nu=1.8$ mM is a scaling parameter. The LacI-promoter binding rate $k_{2}$ is determined by:

[TABLE]

where while the repressor binding range is determined by $k^{min}_{2}=0.01\text{ molecules}^{-1}$, $k^{max}_{2}=0.2\text{ molecules}^{-1}$ and $\mu=2.5\%$ (in % w/v) and $\lambda=0.035$ mM is a scaling parameter. The cooperativity of inducer driven binding of the proteins to the promoter is modeled by a second order dependence of both $k_{1}$ and $k_{2}$ on [ara], [iptg], $\nu$ and $\lambda$.

The mRNAs are translated into unfolded proteins $A_{uf}$ and $R_{uf}$, the monomeric activator and repressor respectively. The dynamics of translation of the unfolded protein species ($X_{uf}$) is given by the general expression-

[TABLE]

where $\rho_{X}$ is the translation rate of X that corresponds to $A_{uf}$ and $R_{uf}$. This is also referred to as the ribosome binding site (RBS) efficiency in the synthetic biology literature [10]. The folding rate $k_{f}$ is constant and common to both proteins, while the degradation rate of the unfolded protein $k_{X}$ is variable due the ssrA tag that results in rapid degradation [23] by proteolysis mediated by ClpXp [3]. To account for the fact that the number of ClpXp molecules available to degrade proteins is limiting, the degradation rate for repressor is given as-

[TABLE]

while that for the activator is given as-

[TABLE]

Here, $\gamma_{0}$ is the maximal degradation rate, $c_{0}$ = 0.1 is the concentration of proteins at which the rate of ClpXP is half-maximal and $g$= 2.5 represents the differential degradation of two proteins, i.e. the activator is degraded faster than the repressor. The total copy number of all proteins in this system, $\Sigma P=A_{uf}+R_{uf}+A+R$, thus inversely determines the degradation rate.

Dimers of AraC ($A$) and tetramers of LacI ($R$) are the active form of the proteins that regulate gene expression. We assume the monomers to be in rapid equilibrium with their respective dimers and tetramers, as a result of which the dynamics of AraC dimers are described by-

[TABLE]

and the dynamics of LacI tetramers are described by-

[TABLE]

The protein concentration is thus determined by folding rate $k_{f}$ of unfolded proteins ($A_{uf}$ and $R_{uf}$). The factors $p_{a}=2$ and $p_{r}=4$ account for the dimerization and tetramerization of the activator and repressor respectively, required to account for the rapid equilibrium of folded monomers with their multimers. As a consequence, at any time point, the dimers are half ($1/p_{a}$) the number of monomers and tetramers are 1/4th ($1/p_{r}$) the number of monomers at the same time.

The simplifications of the model include assuming equilibrium of protein dimerization and tetramerization and promoter binding of the activator and inhibitor. Since each of these processes is expected to multiple intermediate states, we account for these in our equations by an explicit delay. Such an approach of using delay to account for intermediate processes (typically transcription) without explicitly modeling them has been used successfully used to model negative feedback oscillators [13] and the lac operon [28].

3 Results

3.1 Minimal model can reproduce oscillatory dynamics from previous experiments

Our model is based on parameters taken from previous reports (Table 1), leaving only two free parameters, the delay terms $\tau_{1}$ and $\tau_{2}$ for the activator and repressor respectively. Since the activator $A$ must only transition through one step, dimerization, to be active, while the active repressor $R$ undergoes two intermediate reactions of dimerization and tetramerization [23], we assume $\tau_{2}=2\tau_{1}$. Having reduced the free parameter to only one free parameter $\tau_{1}$, we estimate the delay in monomer to dimer transition to be in the order of minutes, based on previous reports [13]. We then numerically optimized $\tau_{1}$ by minimizing the deviation between predicted and experimentally reported period of the oscillations by Stricker et al. [23], for [IPTG] = 2mM and 0.7% arabinose. A scan of varying it from $\tau_{1}$ from 0 to 5 minutes with 0.25 step-size resulted in an optimal value of $\tau=1.25$, at which the sustained oscillations in the concentrations of AraC and LacI mRNA and protein complexes are observed (Fig. 2(a)), with 0.7% arabinose and 2 mM IPTG as inducer concentrations. The qualitative nature of the oscillations is comparable to that of the detailed model of Stricker et al. in previous work [23]. The oscillations show the same period and amplitude, independent of initial conditions.

The so called ‘degrade-and-fire’ (DF) behaviour of the oscillator observed [13] is consistent with the mechanism of a basal mRNA transcription rate ($b$) producing low concentrations of protein, that are constantly degraded by proteolysis, determined by $\gamma_{0}$ and $c_{0}$. Since the production of AraC dimers is faster than the LacI tetramers, the activator protein induces network expression, initially at a slow rate, and then at increasing rates via the positive feedback loop (‘fire’). The repressor protein concentration is however also gradually increasing, and once it reaches a threshold required for promoter repression, the negative feedback loop begins to inhibit transcription. This negative feedback reduces mRNA levels, and combined with the mRNA degradation rates $k_{m}^{a}$ and $k_{m}^{r}$, turns the network off (‘degrade’). The basal transcription rate then once more induces the cycle, resulting in oscillations.

We find our minimal delay differential equation model predicts that the time period of the oscillator responds to IPTG and arabinose induction differently (Fig. 2 (b)). The period increases for 0 to 1 mM [IPTG] and then decrease for higher values when arabinose is constant (0.7%). On the other hand, increasing arabinose concentrations from 0 to 3%, with constant IPTG (2 mM) results in an increase and saturation of the period. Surprisingly, these model predictions closely match previously reported experimentally determined values [23]. Indeed both the period and amplitude of oscillations show similar behavior over wider ranges of both inducer concentrations (Fig. 2 (c)).

Since the ability of the ‘fire, inhibit and degrade’ model to produce oscillations depends on an asymmetry in activator and repressor molecules, we tested the sensitivity of the model to both the extent of asymmetry and whether DNA asymmetry can be mimicked by translational asymmetry.

3.2 Effect of gene copy number asymmetry on oscillations

In previous experimental and theoretical work where the lac-ara oscillator was first developed [23], the genes encoding the activator and repressor were expressed from two different plasmids that were maintained in cells at different copy-numbers. The gene copy numbers of DNA molecules encoding activator was 50 and repressor was 25, a ratio of $n_{a}:n_{r}=2:1$. Our model uses the same values in order to reproduce experimental data (Fig. 2). However, both cell to cell variability due to molecular states, as well the differences in the stage of division, could result in variations in plasmid copy numbers. Additionally, the precision of absolute plasmid copy number control cannot be ruled out. As a result, we have examined the effect of varying the asymmetry of the DNA copies of activator ($n_{a}$) and repressor ($n_{r}$) on the oscillatory dynamics. We find that when the activator gene copies are two-fold in excess of or equal to those of the repressor, the oscillations are rapid and amplitudes comparable. However, a two-fold excess of repressor copies results a small increase in the period of oscillation of both activator protein $A$ (Fig. 3(a)) and repressor $R$ (Fig. 3(b)). Indeed the asymmetry, $\rho_{C}=(n_{a}-n_{r})/(n_{a}+n_{r})$, appears to drive the increased amplitude of that gene, i.e. if the asymmetry involves an excess of lacI DNA, then the repressor protein $R$ has a higher amplitude than the activator $A$, and vice-a-versa.

A systematic scan across $\rho_{C}$ ranging from -0.4 to 0.4 (repressor DNA in excess to activator DNA in excess) demonstrates a continuous decrease in period, i.e. increase in oscillation frequency (Fig. 3(c)), while repressor amplitude reduces and activator amplitude increases (Fig. 3(d)). This would suggest that at some higher factor of asymmetry (e.g. $\rho_{C}=-1$), the oscillator period would become very large and not allow for ‘rapid’ oscillations, within a generation of the bacterium. Indeed it suggests that a DNA copy number asymmetry ensures fast oscillations. If the oscillator were to be re-designed with an equal copy of both genes, we proceeded to ask if the asymmetry of mRNA translation could replace the asymmetry in DNA copies.

3.3 DNA copy number and RBS efficiency have an equivalent effect on protein oscillations

We aimed to address the question whether a symmetric DNA copy number ratio ($\rho_{C}$) can be replaced with an equivalent mRNA translation asymmetry, based on ribosome binding site (RBS) affinity, where the RBS-based asymmetry is $\rho_{R}=(t_{a}-t_{r})/(t_{a}+t_{r})$. To answer this question, we attempt to derive an expression for protein translation in terms of the transcription and translation terms. We start from the transcription rate equation for the mRNA production of activator ($m_{A}$) and repressor ($m_{R}$) expressed in terms of the general equation for $m_{X}$ as-

[TABLE]

We now use the integrating factor method to write the functional form of the solution-

[TABLE]

where the notation of $k_{m}$ is assumed to correspond to $k_{m}^{a}$ and $k_{m}^{r}$. Substituting this in the protein translation equation that results in unfolded protein ($X_{uf}$) from mRNA, we obtain-

[TABLE]

Since, $n_{X}$ is a constant, it is taken out of the integration resulting in the expression of the rate of unfolded protein formation-

[TABLE]

This expression has a multiplicative factor consisting of two constants $n_{X}$ and $\rho_{X}$, demonstrating the exact equivalent of DNA copy number and RBS-efficiency. Since the production of RNA functions as a simple multiplicative factor of the DNA copy number, we further tested whether the degradation rate would potentially affect the dynamics any differently.

3.4 Robust design of oscillator based on mRNA degradation rate asymmetry

Based on these results, excessive expression of the activator appears to be important to ensure the rapid (small period) oscillations of the proteins with significant amplitude. However, since the DNA copy number asymmetry ($\rho_{C}$) based on plasmid copies can be subject to variability, an alternative means of generating this asymmetry needs to be considered. Since the effect of RBS efficiency on protein translation can be factorized in terms of DNA copy numbers (i.e. they have the same effect on oscillations), we proceeded to explore the role of mRNA degradation rate asymmetry $\rho_{D}=(k^{a}_{m}-k^{r}_{m})/(k^{a}_{m}+k^{r}_{m})$, as seen in Fig. 4(a), on oscillations.

We find the effect of increasing IPTG concentration result in an almost linear decrease in period of oscillations (higher frequency), while increasing arabinose concentrations increases the amplitude, for a fixed DNA copy numbers asymmetry of $\rho_{C}$=0.33, i.e. $n_{a}$ is 50 and $n_{r}$ is 25 (Fig. 4(b)). If on the other hand, the $\rho_{C}$ is itself varied, then we find a dependence of oscillation period varies on the asymmetry, for a given IPTG or arabinose concentration (Fig. 4(c,d)). If we now vary the mRNA degradation rate asymmetry $\rho_{D}$, we find the oscillator explores a wider period (20-100 min) as compared to the copy-number asymmetry (Fig. 4(e),(f)). This would appear to suggest that while the experimental system based on two plasmids can produce fast, tunable oscillations in proteins, their robustness to cell to cell variability can be improved by expressing the proteins at the same level with an asymmetry in activator and repressor introduced in terms of mRNA degradation rates, that result in lifetime differences.

4 Discussion

We have successfully developed a minimal model of the dual feedback loop tunable genetic oscillator, which still maintains sufficient detail to explore the regulation of the system at the levels of the canonical central dogma- DNA, RNA and protein. We have optimized a single free parameter, the delay term to obtain oscillatory dynamics that match the whole range of experimental values. Using this minimal model, we explore the role of asymmetry of activator and repressor, in terms of DNA copy number and the production and degradation rates of mRNA. We find that while rapid oscillations can be produced when the a two-fold excess of activator genes is present, this asymmetry is mathematically identical to mRNA production rates. On the other hand, mRNA degradation rate asymmetry results in a wider tunability of the oscillator. Thus, we believe with this simplified model of a canonical synthetic genetic oscillator, we predict a more robust design of the oscillator, capable of exploring a wider range of frequency and amplitudes, than previously described.

While attempts to simplifying the lac-ara dual-feedback oscillator have been made, they have ignored the role of mRNA modeling only the DNA and protein components [26, 24]. As a result, the ability to explore novel network design is limited. In their work, multiple parameters are fit to experimental data, which could result in artifacts. In our work we have maintained identical parameters to the published model and experimental work of Stricker et al. [23], as well as used experimentally tested ranges of inducers (IPTG and arabinose), while leaving only one free parameter, the delay term. We believe thus to have improved upon previous attempts at simplification of this model oscillatory network.

While we have explored the role of asymmetry in oscillatory dynamics, it is clear from our simulations that copy number and translation symmetry can also result in oscillations. This is consistent with the results of a model of where dual-feedback oscillator with symmetric components can still oscillate as described by Maeda et al. [12]. Indeed a comparison of genetic oscillator designs suggests that symmetric gene networks while used in models are unlikely to be natural [21]. Indeed the assumption of many standard genetic oscillator models such as the ‘repressilator’ model is that mRNA degradation rates are comparable for the component genes [4]. Indeed recent work has demonstrated the possibility to achieve fine-tuned control over protein levels through mRNA degradation rate modulation, simply by varying the length of the polyA tail [1]. Here by exploring the possible differences in copy numbers, RBS efficiency and mRNA stability of the two-component network, we have explored more biologically realistic properties of networks, and find the effects of asymmetry depend on the level in the information flow at which they are implemented- DNA, RNA or protein.

The role of the host cell machinery in determining aspects of the genetic oscillator dynamics has been explored previously showing negative-feedback loop oscillators such as circadian oscillators can phase lock to the cell cycle [19], which natural oscillators avoid by using protein-phosphorylation, asynchronous replication and ‘noise [20]. This would suggest the robustness of the lac-ara oscillator seen in experiments could be the combined result of the positive feedback loop as well as the noise coming from replication asynchrony and copy number variability [23]. Experiments with the oscillator genes are localized on the same plasmid, to separate the effects due to copy number variability and other sources of noise, could potentially answer some of these questions.

Acknowledgments

We wish to acknowledge discussions with Dr. Pranay Goel regarding delay differential equations. YKJ is supported by a project assistantship from a grant BT/PR16591/BID/7/673/2016 from the Dept. of Biotechnology, Govt. of India to CAA.

5 Appendix

5.1 Numerical Simulations & Analysis

In summary, we have reduced the mathematical model given by Ref. [23], with more than 30 equations, to a simplified model having only six delay differential equations. We use two delay terms, $\tau_{1}$ (activator) and $\tau_{2}$ (repressor) and relate them as $\tau_{2}=2\tau_{1}$. This is based on the fact that activator proteins bind to DNA after dimerization, a 1-step process, but the repressor binds after tetramerization, involving two intermediate steps. Since the rate constants for all these reactions have the same value according to previous reports [23], and promoter binding is much faster than multimerization, we argue that the delay due to two consecutive reactions is twice that of a single reaction. This model was simulated using Python 3.6 and for the purpose of solving delay differential equations, PyDDE package was used. To solve a delay differential equation numerically, the solver needs to be given initial data of all the variables between time $0<t<2\tau$. To get around this problem, we solve the equations as simple ODEs between $0<t<2\tau$ and as DDEs when $t>2\tau$

A custom written peak-finding algorithm was used for estimation of time period and amplitude. This algorithm involved comparing the concentration value at each time point with that of some time before and some time after the particular point. If $y(t)$ denotes the time series of protein concentrations, for any time $t$, we check if-

[TABLE]

If the above condition is satisfied for all $s\ \in\ \ \{\delta t,2\delta t,30\delta t\}$ then we say that $t_{0}$ is a ‘peak’. Our calculations were run with $\delta t$ values of 0.2 min. The period and amplitude can be found by averaging the time difference between consecutive peaks and the values of $y(t_{0})$, respectively.

5.2 Model derivation

We have derived the model equations by simplifying the detailed reaction kinetics of (a) promoter dynamics determined by protein binding, (b) RNA transcription, (c) translation of mRNA and (d) protein folding and (e) protein multimerization. The promoter reactions in the AraC-LacI oscillatory system are as follows-

[TABLE]

[TABLE]

[TABLE]

where $A$ represents AraC protein dimers, $R$ represents LacI protein tetramers, $P^{a/r}_{i,j}$ represent the states of promoters on the (a)ctivator/(r)epressor plasmids with $i\in(0,1)$ number of AraC dimers bound and $j\in(0,1,2)$ the number of LacI tetramers bound. In contrast to the previous work [23], we consider$k_{1}$ and $k_{2}$ to be the equilibrium constants, and equal to the ratio of rate of forward reaction to the rate of backward reaction. The mRNA transcription reactions are-

[TABLE]

[TABLE]

where $m_{a/r}$ represents the number of mRNA molecules of araC/lacI genes. When an AraC dimer is bound to the promoter, i.e. when promoter is in state $P^{a/r}_{1,0}$, transcription rate is increased by a factor of $\alpha$. When any number of LacI tetramers are bound, i.e $P^{a/r}_{i,1}$ or $P^{a/r}_{i,2}$, there is no transcription.

The reactions representing mRNA translation and protein folding are-

[TABLE]

[TABLE]

[TABLE]

[TABLE]

$A_{uf}$ and $R_{uf}$ are the unfolded activator and repressor proteins, $A$ is the dimeric AraC protein and $R$ is the tetrameric repressor protein. Here, unlike in previous work, we have assumed rapid equilibrium between folded monomer proteins, dimers and (in the case of LacI) tetramers.

Degradation reactions for molecular components-

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We have ignored protein looping, which is considered by Stricker et al. (2008) in their model.

Based on previous work where equilibrium assumptions were made while modeling genetic networks [7, 6]. Every plasmid has one promoter and that promoter can be in total six possible states. The different promoter states bound to the transcription factor dimers and tetramers is given in the following equilibrium equations-

[TABLE]

Here $A_{2}$ refers to AraC dimers (we refer to as $A$) and $R_{4}$ refers to LacI tetramers (we refer to as $R$). We have suppressed the promoter superscript $a$ and $r$ to simplify notation, but the following equations are independently valid for the promoter upstream of activator and repressor. Since the total number of all promoter states added together is constant and determined by the plasmid copy number $n$, so we impose the condition-

[TABLE]

The production of mRNA is governed by the differential equations -

[TABLE]

Rewriting $P_{1,0}$ in terms of $P_{0,0}$ using 14, and then writing $P_{0,0}$ in terms of $n$ using Eq. 15, into the above expression (Eq. 17) we get the following differential equation-

[TABLE]

Even though we have almost derived the equation stated in the main equation, it remains to determine $A_{2}$ and $R_{4}$, which can only come from the protein production equation. Hence, the feedback present in the system is also mimicked in our model. To proceed towards the dynamics of protein dimers and tetramers required, we must first describe the intermediate dynamics of newly translated but as yet unfolded proteins-

[TABLE]

Here, the first term represents the production of unfolded proteins by translation, the second and third term represent depletion of proteins due to folding and degradation respectively. We can now write a reaction for protein monomer production as-

[TABLE]

where $A^{\prime}$ and $R^{\prime}$ are the numbers of monomeric activator and repressor proteins produced by translation, respectively. Our aim is to describe the dynamics of protein dimers and tetramers. We again use the equilibrium assumption for this purpose. Let the number of monomers of AraC and LacI proteins remaining in the system after dimerisation and tetramerisation be $A^{\prime}$ and $R^{\prime}$ for activator and repressor respectively and the number of dimers of AraC and tetramers of LacI thus formed be $A_{2}$ and $R_{4}$ respectively. Since the total number of protein molecules ($A_{tot}$, $R_{tot}$) is conserved, we have-

[TABLE]

where all the components are in ‘monomer equivalents’. A differential equation for active forms of the proteins can be written for activator dimers as-

[TABLE]

and for repressor tetramers as-

[TABLE]

To simplify the notation, we again redefine $A=A_{2}$ and $R=R_{4}$. We also introduce delay into the first equation heuristically, using the argument that the intermediate equilibrium reactions of dimerisation and tetramerization take some time. The kinetics are given by the forward rates of dimerization of activator ($k_{d}^{a}$) and repressor ($k_{d}^{r}$) and tetramerisation of repressor ($k_{t}^{r}$) The promoter binding reactions may also contribute to the delay. These rates are taken from the previous study by Stricker et al. [23] are-

Thus, given the slow rate of dimerisation and tetramerisation compared to promoter binding rates (Table 2), we are justified in the assumption that all the delays in our equations are due to the protein multimerization. Based on the identity of values of $k_{d}^{x}$ and $k_{t}^{r}$, the delay in the feedback of the tetramer is twice that of the delay in the dimer. This allows us to arrive at six equations that fully describe the system, as described in the main text.

6 Supplementary material