Flexible Dynamics of Two Quorum Sensing Coupled Repressilators

Genetic oscillators play important roles in cell life regulation. The regulatory efficiency usually depends strongly on the emergence of stable collective dynamic modes, which requires designing the interactions between genetic networks. We investigate the dynamics of two identical synthetic genetic repressilators coupled by an additional plasmid which implements quorum sensing (QS) in each network thereby supporting global coupling. In a basic genetic ring oscillator network in which three genes inhibit each other in unidirectional manner, QS stimulates the transcriptional activity of chosen genes providing for competition between inhibitory and stimulatory activities localized in those genes. The"promoter strength", the Hill cooperativity coefficient of transcription repression, and the coupling strength, i.e., parameters controlling the basic rates of genetic reactions, were chosen for extensive bifurcation analysis. The results are presented as a map of dynamic regimes. We found that the remarkable multistability of the anti-phase limit cycle and stable homogeneous and inhomogeneous steady states exists over broad ranges of control parameters. We studied the anti-phase limit cycle stability and the evolution of irregular oscillatory regimes in the parameter areas where the anti-phase cycle loses stability. In these regions we observed developing complex oscillations, collective chaos, and multistability between regular limit cycles and complex oscillations over uncommonly large intervals of coupling strength. QS-coupling stimulates the appearance of intrachaotic periodic windows with spatially symmetric and asymmetric partial limit cycles which, in turn, change the type of chaos from a simple anti-phase character into chaos composed of pieces of the trajectories having alternating polarity. Abstract truncated...


I. INTRODUCTION
Synthetic genetic networks provide researchers with the opportunity of designing biologically based circuitry for accomplishing specific functions. In principle, such circuits can be incorporated into natural cellular machinery or used in an entirely synthetic environment.
Oscillators are essential motifs for circuit design [1]. The repressilator is a synthetic genetic oscillator (GO) in the form of a ring of three genes sequentially inhibiting one another's transcription. The GO has been inserted experimentally into E. Coli [2] and extensively studied theoretically via deterministic and stochastic approaches [3][4][5].
Coupling of individual GOs is important in coordinated activity of GO populations.
Bacterial Quorum Sensing (QS) [6], which provides for cell-cell communications in bacterial populations by fast diffusion of small specific molecules (autoinducers), is a natural candidate for the role of being a manager for synthetic genetic network dynamics. This idea was effective, for instance, in constructing several synthetic multicellular systems like the ecological predator-prey [7], population control [8], and other models (see Refs. in recent review [9]).
QS-coupling has been used for GO synchronization in different mathematical models [10][11][12], as well as in experimental demonstration of a multicellular clock [13]. Recently, the QS system has been used to construct a synthetic microbial consortium with population-level oscillations [14]. Although autoinducers freely diffuse between cells with GOs inside, this type of coupling as a whole cannot be reduced to simple linear diffusive coupling which still dominates in studies of coupled oscillators. An autoinducer is not a required element of a GO, its production may be controlled by one gene of the GO networks but its target is transcription regulation of a different gene in the same or neighboring cells. Being physically diffusive, this coupling mechanism is, however, difficult to present as just one classical diffusive term in the ODE system describing the genetic circuits. To write down the coupling term mathematically it is necessary to consistently consider the particular reaction scheme.
This way to draw the coupling term seems quite natural for biochemical and genetic networks in which communication is a cumulative multistep process which includes regulative production, diffusion, binding to targets, and other important metabolic steps.
The first attempt to demonstrate theoretically the possibility of in-phase synchronization of repressilators was successfully realized [11] and checked in a simple electronic model [15], but not all collective regimes were detected within the limits of the coupling design considered. Later publication demonstrated that the use of limited parameter space masked other dynamical regimes. The anti-phase limit cycle (APLC) and complex periodic regimes were revealed when stronger repression in transcriptional regulation and greater differences in time scales of mRNA and transcription factor kinetics were taken into account [16].
In the next version [17,18] of QS-dependent cell-cell interaction, identical repressilators were coupled using a modification of the additional plasmid for the QS mechanism. The modification provided phase-repulsive interaction between oscillators which leads to a rich set of stable attractors: periodic regular anti-phase limit cycle (APLC time series are shifted and the degree of cooperativity (n) of transcription factor (x) binding to promoter. Previous publications [17,18] concentrated on the dynamics with small Hill repression n and limited values for transcription rates, time scales ratios for mRNA and proteins kinetics, and QS signaling molecule (autoinducer) activity as transcription activator. The goals of this paper are to significantly extend the main parameter areas within the limits of one model of QScoupled identical or near identical repressilators [18] to detect new dynamic behavior(s), to present the coarse-grained structure and content of the phase diagram (the map of regimes), and to investigate the robustness of multistability with respect to parameter values. We will use a reduced version of the mathematical model, as well as its electronic circuit model [19,20] adapted for two coupled 4-dim repressilators. The effectiveness of this analog version has been checked recently in the study of a single repressilator with QS feedback [21]. Use of the inherently different numerical and electronic models provides a test of the robustness of the dynamics. We consider a wide range of control parameters, even going beyond 4 the characteristic limits known for them today. The results obtained using the extended ranges may be viewed as predictions of the model for the future, when these extended limits become accessible with the fast development of synthetic genetics. The dynamical results produced by the QS-coupled oscillators used here may inspire application to other systems beyond synthetic gene circuits, such as electronic ones (see below) or chemical systems of interacting water droplets bearing oscillating BZ reactions coupled via diffusion in selective oil environments [22].
We find that for a small degree of repression of transcription factor binding (n < 3) antiphase LC is a single homogeneous periodic attractor, and that this attractor loses stability via torus bifurcation if n increases (n = 3) and gets replaced by complex/chaotic regimes over a large interval of coupling strengths (Q). This Q-continuation branch of the AP limit cycle and chaos coexists with stable homogeneous and inhomogeneous SS and IHLC providing a rich landscape of attractors in multi-parameter space.
The sizes of parameter's intervals occupied by the complex regimes are unusually large and include many periodic windows which contain both spatially symmetrical and asymmetrical stable limit cycles. The presence of asymmetrical limit cycles results in the appearance of chaotic trajectories built from pieces with randomly alternating "polarity". This type of chaotic time series differs from that of simple anti-phase chaos observed in parameter regions devoid of asymmetrical limit cycles.
The extension of the ranges for parameter analysis also reveals a new type of multistability if changes in the maximum transcription rate (α) are taken into account. For a value of the Hill coefficient near 3 a range of α gives rise to the coexistence of the complex nonperiodic oscillations and the anti-phase limit cycle with 5 return times in one period. Again, the range of this hysteresis covers large intervals of the coupling strengths and α.
A further increase in transcription cooperativity up to n = 4 leads to the recovery of stability of the in-phase limit cycle in addition to stable APLC. In-phase limit cycle is not a dominant regime but its appearance seems important to complete the entire dynamic picture demonstrated by the two coupled identical repressilators.

II. NUMERICAL AND ELECTRICAL MODELS
We investigate the dynamics of two repressilators interacting via repressive QS coupling as used previously [17,18]. Figure 1 shows a single repressilator coupled via QS to the external medium. The three genes in the loop produce mRNA (a, b, c) and proteins (A, B, C), and they impose Hill function inhibition on 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 κ in combination with Michaelis function) production of mRNA for protein C, which, in turn, reduces the concentration of protein A resulting in activation of protein B production. In this way 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, even for just a single repressilator [21].
The original models of a single repressilator [2,17] used re-scaled dimensionless quantities for rate constants and concentrations. We reduce the model for the case of fast mRNA kinetics ((a, b, 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 i = 1, 2 for the two repressilators, β j (j = 1, 2, 3) are the ratios of protein decay rate to mRNA decay rate, α 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 η depends on the permeability of the membrane to the AI molecule. The concentration of AI in the external medium is S ext and is determined according to quasi-steady-state approximation by AI produced by both repressilators (S 1 and S 2 ), and a dilution factor Q.
We also implement the model Eq. (1) in electronic circuits as described previously [19,21], with additional performance improvements [20]. By its very nature the circuit has inherent parameter mismatch and noise, and is free from any numerical artifacts. the circuit by the linear-piecewise-continuous behavior min(0.8S, 1) as described previously [21]. Therefore Eq. (1c) in simulations is replaced by  Figure 2 shows the resulting comparison of measured circuit performance and mathematical model for inhibition and activation.

III. RESULTS
A. Low Repression: n = 2.8 We start with the demonstration of the general structure of the phase diagram (map of regimes) for two coupled 4-dim repressilators (Eq. (1)) with parameters similar to those used previously [18] in a study with non-reduced coupled 7-dim repressilators. We choose protein B as the dynamical variable of interest.  and its associated HB2, which appears where the IHSS becomes unstable between LP2 and LP3. The IHLC extends from HB2 to LP5. Second, the coupling strength Q-range for multistability is shifted and increased. For κ = 12, AP is the only stable behavior for Q < 0.5, whereas for κ = 17 the regime of multistability extends down to Q = 0.1, and the Q-range for 3-state multistability has increased about 4-fold. We note that the LP1 and BP2 points in Fig. 3b are too close together to be resolved on the figure, as are LP3 and LP5 for the upper piece of the IHLC. Also, the LP5 of the IHLC extends to a slightly lower Q-value than does the LP3 of the unstable portion of the IH branch. The sensitivity of the dynamics in Fig. 3 was investigated by varying parameters α and β 1 . Figures SM3 and SM4 show that these dynamics are not restricted to a narrow parameter range.
Circuit results for n = 2.8, κ = 17.8 generally confirm the predictions in Fig. 3b. Figure 4a shows a bifurcation diagram constructed from measured voltage amplitudes for the various dynamics. The APLC-branch is stable to Q = 1.  Figure 5 screenshots show all four transitions corresponding to switching between stable branches in the bifurcation diagram as Q crosses limit points: Figure 5a shows IHLC to APLC transition for Q decreasing at 0.25 (LP5), Fig. 5b shows IHSS to HSS for Q increasing at 0.6 (LP2),  Fig. 5c indicates that the HSS LP1 is sensitive to unavoidable differences between the two repressilator circuits. We also point out that after the transitions in Figs.
5a and 5c the phase-shifts are changing to the 180 • characteristic of APLC.
We note that the Q-ranges for particular dynamics in Fig. 4a show some difference from the simulation in Fig. 3b. However, overall agreement of the structure of dynamics is demonstrated. Also, the circuit only finds stable dynamics (and transitions between stable states shown in Fig. 5), and therefore Fig. 4a does not find the unstable branches (black lines) shown in the Fig. 3b simulations.   chaotic behavior when Q reaches 0.6 whereas for the reduced 8-variable model presented here the AP branch remains stable. However, as shown below, chaotic behavior is found at higher Hill coefficient of repression n. Interestingly, the n-range of stable APLC (green) is restricted between 2.5 and 3.5. Decreasing n below 2.8 will only reduce the amplitude of the stable APLC. Therefore we increase n, expecting the introduction of unstable regimes of APLC (blue in Fig. 7) at the larger coupling strengths. Simulations also show that for n > 3.5 large amplitude APLC is predicted for smaller κ at small Q (essentially uncoupled repressilators), and small amplitude APLC is predicted for larger κ at large Q (κ suppresses and stabilizes the APLC, shown below).
We begin by increasing repression to n = 3.0, with α = 190 and κ = 15. Figure 8a shows the simulated bifurcation diagram. Most of the features are similar to those for the smaller repression case in Fig. 3b. However increasing repression to n = 3.0 has introduced the unstable region (blue) in the APLC branch between TR1 and TR2 as predicted in     shows a screen-shot of the opposite transition, from HSS to IPLC, which occurs in the circuit when Q is slowly decreased through the LP1 of the high-B-HSS. Note that the amplitude of the IPLC is larger than the HSS as predicted in Fig. 11.
It was shown above that increasing the rate of activation κ causes the LP1 of the HSS to move to lower Q-values. This decrease must then also move the end-point of the IPLC branch to lower Q-values, which can result in loss of the stable IPLC. For example, raising κ to 5 with Q=0.8 causes the n-range for stable IPLC to move to n = 3.8 and shrink to only 0.008 width (compared to width 0.25 at n = 4.25 for κ = 4 in Fig. 10 In the region of unstable APLC between the TRs (see Fig. 8a) a rich variety of behaviors are predicted numerically and observed in the circuit for n ≈ 3.0 − 3.4. Figure 13 shows  The main dynamical surprise found between TRs over a large range of parameter α (180 to 300) is the coexistence of the CO regime and periodic limit cycle with five subperiods ("return times") labeled below as 5:5LC. Figure 14 presents the phase portraits of 5:5LC and CO calculated using identical parameters. It is clearly seen that 5:5LC and this type of CO are anti-phase spatially symmetric attractors similar in structure and slightly different in amplitudes.  Figure   16 presents the distributions of return times of the coexisting attractors, calculated using Poincaré sections (B 1 = B 2 = 5) of very long trajectories in the presence of uncorrelated low-level white noise added to variables B i . Integrations were started from initial conditions corresponding to 5:5LC (Fig. 16a) or to CO (Fig. 16b). Noise amplitudes were varied to find the most appropriate value to escape the endless and fast mixing of dynamic regimes. There is an interval of noise amplitudes which stimulates the appearance of visible dispersion of all five peaks (corresponding to the return times in Fig. 16a) in dN/dT distributions, but it is too small to stimulate switching between the periodic 5:5LC and the CO. This means that despite similarities in the extent of phase relations (Fig. 14) and in the boundaries of return times (Fig. 16), the system Eq. (1) has two types of coexisting stable solutions-5:5LC and CO-over large intervals of parameters inside the region of unstable APLC. Figure 15 shows the effect increasing α has on the coexisting 5:5LC and CO attractors. A more in depth examination of the dynamics for n = 3.0 is in SM. Figure SM5 shows additional detail of the evolution of the 5:5LC as α increases. Figure SM6 shows a map of regimes in (Q, α)-space clearly showing coexistence of CO and 5:5LC, and the separation of CO and chaos emerged from PD-cascades. Figure SM7 shows an As5:5LC measured from a circuit.
In summary, raising the Hill coefficient to n = 3 opens new complex dynamic behaviors which are different from classic chaos found earlier for two coupled repressilators [17,18]:  Similar to the case of n = 3, the first regime after the TR1 bifurcation is a quasiperiodic attractor characterized by "beating" of the APLC time series producing two peaks in return times distributions. With further increase of Q, a third peak emerges suggesting the approach of chaos (all within the broken violet line in Fig. 18 from TR1 to 0.735). The more effective indicator of chaos is the emergence of a strongly folded sequential period map presented in Fig. SM8 to demonstrate the process of chaos maturation as Q grows.
Within the region between TR1 and Q = 0.735 there are many very narrow periodic windows which may be occupied by amplitude-symmetrical LCs (e.g. 13:13LC, 35:35LC and other n:n-LCs) or by slightly amplitude-asymmetrical regimes (e.g. the As4:4LC). Many of these regimes in narrow windows demonstrate period doubling bifurcations returning system to chaos. For example, the 5:5LC starts at 0.7116 as a stable symmetric 5:5LC which branches to slightly asymmetric 5:5LC at 0.7141 then starts a PD cascade to chaos at 0.7152. This regime coexists with asymmetrical 4:4LC (see their phase portraits and zoom-in Q-continuations in Figs. SM10 and SM11). In contrast to the 5:5LC, the As4:4LC from Q = 0.710 to 0.719 is "isola" because its Q-interval is limited by LP bifurcations. The chaos emerged through the region of TR1 to 0.735 is characterized as "simple" anti-phase attractor similar to that presented above in Fig. 14b for n = 3, which means that elements of the highly asymmetric behavior arise only for Q > 0.74. at 0.7 < Q < 1.1. If activation is significantly reduced (κ = 4) then α must be increased to get interesting regimes. However, these regimes are then located at Q > 1 if 100 < α < 180.
We conclude that limited but significant variations of basic parameters in 3-dim space do not discover other qualitatively new dynamic regimes compared to those described above.

IV. DISCUSSION
Using simple 4-dim mathematical and electronic models of identical ring oscillators we demonstrated that the realistic design of quorum sensing coupling proposed earlier [18] leads to remarkable multistability even for a pair of oscillators. A QS-coupling system can be incorporated into genetic circuits in different ways, acting solely as a coupling agent (as in the case of repressilators [11] and Relaxators [10]), or it may also be an integral part of an individual oscillator [12,13,24]. In any case, the additional differential equation for autoinducer contains a term describing its linear diffusion but the resulting coupling is nonlinear for the oscillators because of a delay in the autoinducer production and because the production of one of the oscillator's variables is activated nonlinearly. Therefore it is natural to expect unusual collective modes for QS-coupled oscillators, as was demonstrated in [18]. However, that study of multistability considered mainly the role of coupling strength, with the other parameters of the single oscillator like repression cooperativity and the transcription rate being fixed. Therefore, the generality of QS-induced multistability was still obscure because the origin of many new regimes can be traced to modulations in oscillation amplitudes and the "steepness" of repression which are strongly dependent on the dynamic properties of the isolated oscillators. Starting from intermediate α (α > 250 for our parameter set) the emergence of chaos at low Q and its extinction at higher Q are the result of sequential branch point and period doubling bifurcations of this 5:5LC. Each local bifurcation of our system on the (α, Q) plane for n = 3 is well known but the structure of the map of regimes as a whole is unusual to the best of our knowledge.
We found the interval of n-values where chaos in two QS-coupled repressilators is not an exotic regime but instead exists in a large 3-dim parametric space. Moreover, a further increase in n up to n = 3.15 introduces new behaviors inside the Q-interval with the unstable regimes, including regular and irregular asymmetrical ones (Figs. 17 and 19). All collective modes presented in our study are the result of competition between repression (Fig. 2a) and activation (Fig. 2b) of transcription of one specifically selected gene. These antagonistic impacts have different dependencies on the values of the variables: repression is cooperative while activation is linear with saturation. This means that variation of one parameter, e.g.
Hill repression n, may require the other parameters be tuned to escape the transition of the system to a simple attractor like HSS.
We investigated several sets of parameters and found that the evolutions of the complex regimes as a function of coupling strength are qualitatively very similar. Then we looked over the typical dynamics of chaos and discovered its basic skeleton although many details are still unclear. The most impressive result is the existence of several types of chaos: (1) simple anti-phase chaos riddled with very narrow periodic windows of symmetric regimes of the n:n-LC type; (2) asymmetric chaos resulting from period doubling of regular asymmetric regimes, e.g. As1:2LC and As2:3LC; and (3) symmetric chaos consisting of symmetric and asymmetric parts with alternating polarities. The last type of chaos is the most flexible in its specific manifestations and is closely linked to intermittency which is also observed but only near the boundaries between the regimes.
The possibility for chaos due to linear diffusion between identical regular oscillators is known since the middle of the 1970s [25][26][27][28], but the parameter regions of chaos are typically very narrow. Chaos via torus bifurcation of anti-phase limit cycle has been observed in two Rössler oscillators coupled by cross-diffusion [29]. To the best of our knowledge there is only one example of chaotization due to nonlinear coupling of identical biochemical oscillators [30] where "bichaoticity" was demonstrated in narrow parameter intervals without analysis of other regimes inside chaos.
A large cooperativity n (around 4) opens the gate in parameter space for the in-phase limit cycle which becomes stable for large Q before its transition to HSS via infinite period bifurcation. It is still not clear why phase-repulsive coupling permits the formation of APLC/IPLC switching. Intuitively it may be due to the increase in repressilator "stiffness" for large values of n. Similar frequency trigger covering a wide range of parameters was observed in the system of two FitzHugh-Nagumo oscillators coupled via the recovery variable if their stiffness was large [31].
All the results discussed above have been obtained for sets with different β i . The oscillators are identical but move along the partial limit cycles nonuniformly because the kinetics of the specific variable in one element of the ring is faster than in the other elements. Variation in β i is biologically motivated given that there are parameters' differences between the three repressilator genes, their promoters, and degradation rates of the corresponding transcription factors. This internal variety in β i may enhance the coupling-induced diversity of the attractors. To check this possibility we simulated the dynamic behavior in the 3-dim parameter cube (n, α, κ) using several sets of β i varying in the degree of equalization (e.g. 0.25, 0.15, 0.15; 0.3, 0.25, 0.15) (data not shown). The only general effect of smoothing the β i is the necessity to tune reasonably the other parameters but the qualitative picture of the phase diagrams is not changed.
No doubt, the model in Eq. (1) is a strongly reduced dynamic toy which has limited experimental support for formulation of exact forms of functions for the basic repressilators' ring and especially for the subsystem describing coupling. However, the general design formalized in Eq. (1) is realistic and realizable by the methods of synthetic genetics. Therefore, we propose that a pair of coupled repressilators may be viewed as a prototype of a flexible generator with remarkable dynamic behavior which may be useful beyond the construction of synthetic genetic networks.