Electronic Implementation of a Repressilator with Quorum Sensing Feedback

We investigate the dynamics of a synthetic genetic repressilator with quorum sensing feedback. In a basic genetic ring oscillator network in which three genes inhibit each other in unidirectional manner, an additional quorum sensing feedback loop stimulates the activity of a chosen gene providing competition between inhibitory and stimulatory activities localized in that gene. Numerical simulations show several interesting dynamics, multi-stability of limit cycle with stable steady-state, multi-stability of different stable steady-states, limit cycle with period-doubling and reverse period-doubling, and infinite period bifurcation transitions for both increasing and decreasing strength of quorum sensing feedback. We design an electronic analog of the repressilator with quorum sensing feedback and reproduce, in experiment, the numerically predicted dynamical features of the system. Noise amplification near infinite period bifurcation is also observed. An important feature of the electronic design is the accessibility and control of the important system parameters.


Introduction
Biological networks such as transcriptional genetic networks and metabolic networks are complex in structure and dynamical behaviors [1]. A systematic understanding of the functional behaviors of such networks is challenging both from theoretical and experimental viewpoints. Alternatively, an engineering approach [2][3][4] has been undertaken, in recent times, to simulate desired biological functions by designing small genetic units. The main target is to integrate large numbers of such units to derive complex biological functions [5][6][7]. This is like engineering small integrated chips to build a computer to perform a desired function. One of the important discoveries, in this direction, is the design of an artificial genetic network known as a repressilator [4] which consists of a ring of three genes inhibiting each other in cyclic order that shows oscillatory behaviors. Theoretical studies of the single deterministic, stochastic, and even electronic repressilator have attracted attention of many investigators [8][9][10][11][12][13][14][15][16][17].
On the other hand, quorum sensing (QS) [18,19] is a typical process of communication in a bacterial colony, and it has been used as a mechanism for coupling synthetic genetic oscillators [20][21][22][23][24][25]. The QS is accomplished by diffusive exchange of small autoinducer (AI) molecules which participate in the intercellular coupling as well as in self-feedback. Recently, the effective diffusion of AI has been demonstrated experimentally to induce coherence of oscillation in a small colony of synthetic gene units [25]. It is important to understand the effect of the AI feedback on a single repressilator before considering the design of a network of multiple units. We focus on the effect of the self-feedback that QS may provide to a single isolated repressilator and try to understand its dynamics. Details of the model are reported recently [26].
Here, we report a variety of interesting dynamics in the isolated repressilator with a QS type feedback such as multi-stability of limit cycle with stable steady-state (SS), multi-stability of different stable steady-states, limit cycle (LC) with period-doubling (PD) and reverse period-doubling (RPD) and infinite period bifurcation (IPB) transitions for both increasing and decreasing strength of quorum sensing feedback. A noise amplification near the IPB is also observed. Our results are based on both numerical simulation and experimental measurement. Experiments are accomplished by designing an electronic circuit analog of the repressilator with QS feedback. The electronic design is carefully made so that the key system parameters are accessible and controllable, which is not easy in real biological experiments. The electronic analog of the synthetic genetic network is thereby able to reproduce the important dynamical features predicted numerically. The electronic bench-work provides an actual physical system, including the presence of intrinsic system noise, external noise, and device mismatch and, also produces results free of numerical artifact. As such, we use both numerical simulations and circuit measurements as complementary approaches for investigating the dynamics of the particular network topology. This knowledge about the dynamics of the network may be used later for planned biological experiments.
We present the model and the numerical results first, followed by the electronic circuit and experimental measurements. Details of the circuit analysis and derivation of the relation between circuit parameters and model parameters are given in Appendix S1.

Model and Numerical Results
We first investigate the dynamics of a single repressilator with QS feedback using numerical tools. The starting point is the repressilator with QS feedback [15] as shown in Fig. 1. The three genes in the loop produce mRNA (a, b, c) and proteins (A, B, C), and they inhibit each other in cyclic order by the preceding gene. The QS feedback is maintained by the AI produced (rate k S1 ) by the protein B while the AI communicates with the external environment and activates (rate k) production of protein C which, in turn, reduces the concentration of protein A resulting in activation of protein B production. The protein B plays a dual role of direct inhibition of protein C synthesis and AI-dependent activation of protein C synthesis, resulting in complex dynamics of the repressilator.
The original model of a repressilator [4,15] used re-scaled dimensionless quantities for rate constants and concentrations. We reduce the model for the case of fast mRNA kinetics (a, b and c are assumed in steady state with their respective inhibitors C, A, B, so that da/dt = db/dt = dc/dt < 0). The resulting equations for the protein concentrations and AI concentration S are, where the b i is the ratio of protein decay rate to mRNA decay rate, a accounts for the maximum transcription rate in the absence of an inhibitor, n is the Hill coefficient for inhibition, and k S0 is the ratio of AI decay rate to mRNA decay rate. The diffusion coefficient g depends on the permeability of the membrane to the AI molecule. S ext is the concentration of AI in the external medium which is taken to be zero in our case of a single isolated repressilator. We choose parameter values similar to previously used values shown to be experimentally reasonable taking into account realistic biochemical rates and binding affinities [4,15,23].
We explore a wide range of values for rates a and k.
The strength of QS feedback is determined mainly by the activation rate parameter k. Recently, it has been shown [17] that the dynamics of coupled repressilators depend strongly on the time scale ratio b. Therefore we investigate how the dynamics of a single repressilator with QS feedback depends on k for a small value b = 0.1, which is close to biologically motivated value, and for the larger value b = 1 which assumes the same time scales for both the proteins and the mRNA. A b-value near 1 was realized in the biological repressilator [4].

Numerical Results
We first present the numerical results of available dynamical regimes as stable steady state (SS) and limit cycle (LC) as calculated from Eq.(1) using the XPPAUTO [27] for b 1 = b 2 = b 3 = b = 0.1 (slow protein kinetics) and n = 3.3. Figure 2a shows the period of LC states under different transcription rates a, as well as the protein B concentration of SS for a = 30, both as functions of k. For the a = 30 LC, increasing k from 0 causes the activation of protein C production to become so large that a stable SS with high concentration of protein B emerges (upper red line) near k = 14. This high-B-SS coexists with the stable LC over a broad interval of larger k, from about 22 to 36. For k in the range 15 to 22 there is a break in LC continuity where the LC is unstable. Near k = 36 the LC transitions via Andronov-Hopf bifurcation (HB) to a stable SS with low concentration of B (lower red line in Fig. 2a). This second SS is stable over a small range of k (36 to 43) in which it coexists with the high-B-SS, thus behaving as a switch. The operation of the 3-gene ring as a bistable genetic switch over this range of k may be of practical interest. At higher k (.43), only the high-B-SS is stable. Figure 2b shows protein B dynamics, LC amplitude (thick colored lines, left and lower portion of graph) and SS level (thin colored lines, upper and right), for different transcription rates a as a function of k. For each a the strong dependence of the LC amplitude on k suppresses oscillations resulting in the HB to the low-B-SS near k = 36. The low-B-SS level is about 2.5 for all a, whereas the high-B-SS level is approximately a as predicted by Eq. (1b) since protein A is small. For small a (#25) the k-dependence of LC is continuous to the HB in Figs. 2a and 2b because of the low amplitude of LC, but for larger a ($30) the LC becomes unstable and transitions via IPB to the high-B-SS because the larger a has larger amplitude oscillations which allow the phase point to reach the SS. Hysteresis is predicted in Figs. 2a and 2b by the overlapping k-range of the LC and SS for a = 30. The steeply increasing periods in Fig. 2a indicate the formation of the infinite period regimes between two values of k for a $30. IPBs are realized by increasing k from small values as well as by decreasing it from large values. The mechanism of this phenomenon in 4D space deserves separate investigation.
Next we consider b = 1 (protein kinetics same as mRNA and AI). The behavior of SS is not affected by this choice but the evolution of the LC is changed dramatically. Figure 3 shows predictions of LC dynamics for various a as k increases. When a is less than 47 (results not shown here), the LC continuation is continuous period-1 to the HB. However, for larger a the dynamics become more complex. For a = 64 Fig. 3a shows that the LC continuation (green curve) remains stable over the k interval but is more complex due to the appearance of PD followed by RPD bifurcations. Increasing k beyond the HB causes a transition to SS (lower red curve) with all proteins low which coexists with the high-B-SS creating a genetic bistable switch as occurred for b = 0.1. Further increase of k causes a jump to the high-B-SS (upper red curve). For larger a ( Fig. 3b) a break in continuity of LC appears, in addition to more complex PD/RPD cascades. For a = 200, n = 3.43 (Fig. 3c) a second break in LC continuation occurs. Figure 4 shows details of how the LC amplitude changes for all three proteins as k increases for the case b = 1, a = 117, n = 3. The amplitudes of period-1 LCs decrease sharply as k goes from 0 to 15. Bifurcations to period-2 occur at k = 16 and back to period-1 at k = 25, followed by decreasing LC amplitude to SS via HB at k = 40.
Here we summarize the numerical findings. Increasing the strength of QS feedback has two main effects on the repressilator dynamics. Firstly, when the feedback is increased beyond a critical value, a stable steady state (SS) appears with protein B (which controls the AI production) concentration becoming high and the other two protein concentrations becoming low. Multistability exists here between SS and LC. Secondly, as the feedback is further increased, the LC evolves (possibly through unstable regions) with shrinking amplitude until it converges via HB to a different SS characterized by low concentration of all three proteins. Interestingly, the existence of this second SS was previously unknown and allows multistability between the two different SS. The particular dynamics of the evolution of the LC depends on the parameter values. Large transcription rates cause the LC to become unstable for a specific range of QS feedback due to the appearance of IPBs from a LC to a SS. These IPBs occur for both increasing feedback strength from low values and for decreasing feedback strength from high values. A hysteresis between the LC and the SS is seen. A faster kinetics (b = 0.1) of mRNA and QS signaling molecules compared to repressilator protein results simply in period-1 LC dynamics. However, when the QS autoinducer and repressilator proteins have similar time scales (b = 1), then the LC may have complex dynamical behavior which is characterized by the appearance of PD and RPD.

Electronic Circuit and Measurements
Next we report the circuit realization of the genetic network with quorum sensing feedback. Previously we reported an electronic version of the repressilator (e-Rep) [28] which consists of three genes coupled in a ring, inhibiting each other in a cyclic order. The strength of this electronic analog is its accuracy in reproducing the Hill function nonlinearity that controls the rate of transcription. The circuit allows an effective control of the system parameters. As a result, it reproduced the basic properties of an isolated repressilator. Here we make further improvement upon the Hill function nonlinearity in the circuit design and include additional circuitry to imitate the QS feedback. In this way we study how the addition of QS feedback to an isolated e-Rep changes its dynamical behavior and we verify the numerical results. The circuit uses the reduced 4-dimensional model (Eq. 1) of the repressilator in which only three equations for protein dynamics and one for AI are simulated since the mRNA kinetics are assumed to be fast enough to keep mRNA in quasi-equilibrium with its inhibitor.
One can couple many of the e-Rep to simulate various desired functional behaviors thus enabling prediction or design of collective behaviors quantitatively by controlling the circuit parameters analogous to genetic network's parameters. The predicted behaviors can be used for decisive biological experiments. This is a major benefit in studies of e-Rep circuits over the biological experiments where a priori knowledge of the system parameters is not available. The control of the system parameters is a very difficult task in a biological experiment which has additional complications due to its fluctuating medium. We can include additive or multiplicative noise in e-Rep experiments to imitate some types of the fluctuations encountered in real biological experiments. Another target of designing e-Rep is to design new electronic devices (e.g., logic gate, switch) [29] based on the knowledge of synthetic genetic networks. We implement the model (Fig. 1) in an electronic circuit based on the e-Rep circuit used in Ref. 28. Figure 5 shows the electronic  analog of a single gene used here, which is a slight modification of the circuit used in Ref. 28 with an improved Hill function nonlinearity. The relation between the circuit components and the model parameters n, a, b, k, k S0 , k S1 , g are given below. We use k S0 = 1, k S1 = 0.025, and g = 1, in this paper.
The equation for the voltage V i across the capacitor C i in Fig. 5 is obtained by using the Kirchhoff's laws, where I t is the transistor current which depends on the inhibitory input V i21 . The ratio V i /V th is the normalized state variable for protein concentration (i = A,B,C). The relation between the model parameters (a, b) and the circuit parameters are already derived [28] as reproduced here, a = I max R C /V th and b i = C 0 /C i where C 0 defines the dimensionless time t/(R C C 0 ). The Hill coefficient n and the parameters V th and V cth are redefined here as affected by the modification of the circuit in Ref. 28. In Appendix S1, we show that n comes from na = 2.35G 1 G 22 (compared to 1.5G 1 G 22 in Ref. 28) where G 1 and G 22 are unsaturated gains of op-amps U 1 and U 2 , maximum transistor current I max < 3 mA, maximum V i is I max R C = 2.7 V, and V cth < V th +1/(G 1 G 22 ). Numerically simulated and experimental oscillations of the 3-gene e-Rep without QS feedback for b i = (0.1, 1) are shown in Fig. 6 for comparison. Figure 7 shows the circuit for an isolated e-Rep with QS. The triangular blocks with (a,n) and subsequent R and C i are the individual genes of Fig. 5. Details of the QS feedback loop show the connections between protein concentrations B and C corresponding to the feedback loop between B and c in Fig. 1. The equation for the voltage V S across C S corresponding to the AI concentration S is

Circuit with Quorum Sensing
where V ext = 0 since R d is connected to ground in Fig. 7 indicating an isolated e-Rep. Choosing R S2 C S equal to the characteristic time R C C 0 (which sets k S0 = 1) and normalizing by V sth , and then using S = V S /V sth , Eq. (3) is converted, where t is a dimensionless time. We choose the component values such that [1-RS2Vth/(RS1Vsth)] < 1 since [RS2Vth/(RS1Vsth) ,,1 )]. The model parameters (kS1, g) are now expressed in terms of the circuit parameters,  Protein B stimulates production of auto-inducer S that activates production of C. When the voltage at S increases, the inverting opamp output drops, which turns the pnp transistor on and the npn off, thus sending current to protein C voltage. Each triangle with parameters a, n and subsequent R and C i represents the single-gene circuit in Fig. 5 We use R S1 = 100 kV and R S2 = 6.8 kV in the circuit. For a single e-Rep with QS feedback, S ext = V ext /V sth = 0.
In Appendix S1 we show that, k~I S max RC V th~0 :35RC

RK V th :
The circuit approximates the S-function S/(1+S) by a piecewise linear function, min(S/(1+S cr ),1) where S cr is the value of S where the piecewise linear function crosses S/(1+S). We typically choose S cr = 0.25. We show that S cr determines [see Appendix S1] the feedback resistor R Sf by

Circuit Measurements
Now we present measurements from the e-Rep with QS feedback. Numerical predictions in Fig. 2 for n = 3.3 and b = 0.1 show the break in k-continuation of the LC (IPB regime) that occurs when a is increased. Simulations also show that reducing n causes an increase in the critical value of a for the break in kcontinuation of LC. Decreasing n from 3.3 to 3.0 moves the critical value to between 45 and 50 compared to between 25 and 30 respectively in Fig. 2. We reproduce some of the interesting numerical results for b = 0.1, in experiment, as shown in Fig. 8.
The circuit measurements of LC period and protein B LC amplitude as a function of k were made for a = 46 and 60. The resulting continuous curves (green) and curves-with-break (blue) shown are in good qualitative agreement with the simulated behaviors in Figs. 2a and 2b. The IPB region is clear in Fig. 8 for a = 60. Figure 9 shows the LC for a = 46 for the two values k = 14 and 19, demonstrating the decrease in both period and amplitude of LC. We were not able to increase k high enough, in experiment, to observe the HB due to a practical limit on k as discussed in Appendix S1.
The break in k-continuation of LC in Fig. 8 (blue) for b = 0.1, a = 60, n = 3.0 was investigated by slowly changing k (increasing from low value as well as decreasing from high value) until transition of LC to SS was observed. These transitions as shown in Fig. 10 are the IPBs, occurring for increasing k at 9.5 (Fig. 10a) and for decreasing k at 15.8 (Fig. 10b). The hysteresis as predicted in Fig. 2 is also observed in Fig. 10 where the transition from SS to LC occurs for decreasing k at 7.9 (Fig. 10c) compared to 9.5 for LC to SS. The range of k for the coexistence of LC and SS due to hysteresis is clearly from 7.9 to 9.5. For b = 0.1, the only LC observed is period-1, in agreement with the predictions in Fig. 2. For b = 1, a forward and backward period-doubling of LC as shown by numerical simulations in Fig. 4 is realized in circuit measurements for the case b = 1, a = 117, n = 3. Figure 11 shows the evolution of the LC including period-1 and period-2. Plots of measured data are shown instead of oscilloscope images so that   Fig. 4. A steep decrease in amplitude of LC is apparent as k increases from 0 to 12, as is the period-2 LC for k = 19. The small amplitude LC for k = 43 is close to the HB to SS. The lower red line in Fig. 4 for k = 6 to 11 corresponds to a small peak that appears in the plateau shown in the k = 8 case in Fig. 11.
For b = 1 with larger a, a break in k-continuation of the LC is observed as predicted in Fig. 3b similar to what is seen with b = 0.1. Figure 12a shows the IPB for decreasing k from large values for b = 1, a = 138, n = 3. More complex dynamics predicted for larger n (Fig. 3b and 3c) including period-4 LC are seen (Fig. 12b) for b = 1, a = 77, n = 3.4.

Noise Amplification Near IPB
A pre-bifurcation noise amplification [30] is predicted for the repressilator LC as k approaches the IPB. We added noise to the circuit by connecting the voltage V B (analog of protein B) to a noise source instead of to ground through the 1 kV resistor as shown in Fig. 13. The equation for the voltage V B then becomes, where I t is current from the transistor. The source of noise is the breakdown of the reverse biased base-emitter junction in the npn transistor in Fig. 13. The noise circuit is supplied by a well regulated 612 V in order to avoid introducing AC noise. The normalized noise is x noise = V noise /V th giving the equation Ideally, the noise x noise should have zero mean and an amplitude given by its standard deviation. Any dc component of the noise from the reverse biased npn is blocked by the coupling capacitors. Therefore the mean value of the noise is due solely to the nominal 1 mV offset error of the second op-amp. For example, with a = 60, then V th = I max R C /a = (2.7 V)/60 = 45 mV, and the mean value of x noise is 1 mV/45 mV < 0.02. A similarly small value occurs for the minimal amplitude of x noise due to a nominal 1 mV rms noise at the op-amp output. Fourier transforms of the noise show the frequency content to be uniform up to a cut-off around 300 kHz. Using the time-scale RC = 0.1 ms gives a dimensionless period of 1/(300 kHz60.1 ms) = 0.033, which is much smaller than any of the LC periods. We conclude that the noise circuit provides a good approximation to zero mean random noise.
Fluctuations in LC period were measured for different noise amplitudes for three values of k. The highest value, k = 7.9, is close to the value for IPB when other parameters are n = 3, a = 60, b i = 0.1, k S1 = 0.025, S cr = 0.25. Results in Table 1 show the LC's increased sensitivity to noise near the IPB. Far from the IPB (k = 0) there is almost no effect on the LC period as noise is added. Close to the IPB (k = 7.9) fluctuations of the LC's period increase from 0.1 to 1.6 (standard deviation) when noise is added. The increase in period as k approaches the IPB is also apparent.

Conclusions
Synthetic genetic oscillators are usually considered as prototypes of some natural genetic elements, or at least as motifs of genetic networks. Quorum sensing is a natural process of interaction or coupling between the oscillators. A 7-dimensional repressilator with quorum sensing feedback is the simplest genetic ring oscillator which demonstrates limit cycle oscillation over a very broad interval of control parameters. One of the possible motifs made from 3-genes [31] is shown in Fig. 1 and is our present candidate of interest. We made a reduction to a 4-dimensional model, in this work, that exhibits limit cycle behavior very similar to that of the original 7-dimensional repressilator. Additionally, quorum sensing is realized in the model as an additive regulatory function through an additional path in which protein B activates the production of protein C through an autoinducer. It is important that one protein (B) inhibits and activates, at the same time, the subsequent neighboring protein. This competition between inhibition and activation of the one element (C) in the ring of genes results in the coexistence of stationary and complex periodic regimes and opens new possibilities in the regulation of the repressilator.
The first manifestation of autoinducer production is the appearance of a stable steady state with high level of protein B. This steady state may be a single fixed point of the repressilator or may represent the high level of a genetic switch if the activator activity of autoinducer is large enough to allow the existence of a second stable steady state in which all proteins are at low level. The position of steady state does not depend on the time scales of Eq. (1) but oscillatory dynamics of the repressilator depends strongly on them.
For both fast and slow protein kinetics the limit cycle coexists with steady state, the cyclic behavior is suppressed by autoinducer via supercritical Andronov-Hopf bifurcation, and two infiniteperiod bifurcations exist when protein transcription (a) is increased enough to produce a break in the limit cycle's k-continuity. For  slow protein kinetics (b's = 0.1) the autoinducer feedback does not change the form of the limit cycle from period-1. This is not surprising because for small b the equation for autoinducer may be reduced resulting in a 3-dimensional system. For identical time scales (b's = 1) the equation for autoinducer may be considered as a delay mechanism in the appearance of activation of protein C production. It is expected that delay can stimulate complex oscillatory regimes. We show the cascades of period doubling bifurcations, the number of which along k-axis is increased with a. However, these cascades do not result in chaos but end in period-halving and several infinite-period bifurcations.
Further we observed noise amplification in the form of increased fluctuations of limit cycle period (protein B concentration) near the infinite-period bifurcation in the presence of additive noise.
We present an electronic circuit model for a single repressilator with quorum sensing feedback that provides a useful tool which is used in combination with numerical simulations to demonstrate rich dynamics: steady states, switch, limit cycles, infinite-period bifurcations, period doublings. These results for an isolated repressilator under the influence of QS self-feedback should be beneficial for understanding communication between coupled repressilators via QS. In addition, we demonstrate that the electronic analog may be used as an effective test bed for the studies of regulations in genetic networks, and more generally for studies of particular network topologies potentially leading to reverse engineering such as the design of useful electronic devices. Recently we have used the electronic gene in Fig. 5 for investigation of logical stochastic resonance in a 2-gene network [29]. Figure S1 Circuit simulation of inhibitory Hill function for n = 3. Output is the measured normalized current (blue) from the transistor in the single gene circuit (Fig. 5). V is the inhibition input voltage and V th accounts for the inhibitor binding affinity. Analytic Hill function 1/(1+ x n ) (red). (TIF) Figure S2 Circuit simulation of quorum sensing Sfunction. In the QS circuit (Fig. 7) the pnp transistor produces S-function activated protein C current which goes to the protein C capacitor voltage. Measured S-function circuit current (blue), piece-wise continuous linear model min(S/(1+S cr ),1) (green), and S/(1+S) (red). V corresponds to the concentration S of AI and V sth accounts for the binding affinity of AI activator. S cr = 0.25.