A novel optogenetically tunable frequency modulating oscillator

Synthetic biology has enabled the creation of biological reconfigurable circuits, which perform multiple functions monopolizing a single biological machine; Such a system can switch between different behaviours in response to environmental cues. Previous work has demonstrated switchable dynamical behaviour employing reconfigurable logic gate genetic networks. Here we describe a computational framework for reconfigurable circuits in E.coli using combinations of logic gates, and also propose the biological implementation. The proposed system is an oscillator that can exhibit tunability of frequency and amplitude of oscillations. Further, the frequency of operation can be changed optogenetically. Insilico analysis revealed that two-component light systems, in response to light within a frequency range, can be used for modulating the frequency of the oscillator or stopping the oscillations altogether. Computational modelling reveals that mixing two colonies of E.coli oscillating at different frequencies generates spatial beat patterns. Further, we show that these oscillations more robustly respond to input perturbations compared to the base oscillator, to which the proposed oscillator is a modification. Compared to the base oscillator, the proposed system shows faster synchronization in a colony of cells for a larger region of the parameter space. Additionally, the proposed oscillator also exhibits lesser synchronization error in the transient period after input perturbations. This provides a strong basis for the construction of synthetic reconfigurable circuits in bacteria and other organisms, which can be scaled up to perform functions in the field of time dependent drug delivery with tunable dosages, and sets the stage for further development of circuits with synchronized population level behaviour.


Introduction
The aim of synthetic biology is to design and synthesize biological networks that perform desired functions in a predictable manner. These systems then, can be used in conjunction, and be joint to one another to solve a particular real world problem. This is mostly achieved using a combination of genetic logic gates [1][2][3][4][5][6] which use the property of biological parts being analogous to digital logic within certain defined threshold values, which are similar to PLOS  the logic gates designed and used in electrical engineering. Several circuits in synthetic biology have been designed to work as expected using these combinations of digital logic to create circuits with applications in the field of pharmaceuticals [7], biofuels [8], to encode memory in cells [9][10][11][12] or to create biological analogues of electrical circuits such as transistors [13,14]. A significant chunk of work in the field synthetic biology has targetted the construction of two important types of networks, namely switches [15,16] and oscillators [17,18]. The oscillator is crucial to and lies at the center of this paper. Some of the common topologies which are known to show oscillatory behavior are [17,19]: • Goodwin Oscillator [20] • Repressilator [21] • Amplified negative feedback oscillator 1. Goodwin Oscillator-This was the first synthetic genetic oscillator to be demonstrated. It consists of a single gene with negative auto regulation, and had three equations in a simple model, for the mRNA, protein and inhibitor concentrations. This laid the framework for oscillations, as oscillations occur in the Goodwin oscillator if hill kinetics is followed with a high cooperativity coefficient [22]. The hill function approximation was shown to be valid if the quasi steady state assumption for the multisite phosphorylation of the transcriptional repressor was made.
2. Repressilator-The repressilator was the first complete ring oscillator that was demonstrated, and was a three component system, with each component repressing the next (A represses B represses C represses A). This circuit was shown to have oscillatory behaviour in silico and in vivo in E.coli, triggering a large amount of work in the field of synthetic biology, and oscillators in particular. The repressilator however, was not a synchronised oscillator, and since then, further two component oscillators such as the Danino oscillator [23] have shown synchronisation in the oscillations of a colony of 200 cells, while the Prindle oscillator [24] has shown gas-phase redox signaling between colonies to create an LCD like display among colonies that are physically separated from each other.
3. Amplified negative feedback oscillator-Leaving aside the biological implementation, this abstract topology is the one used in the Danino Oscillator. This architecture consists of a two component system, with an activator gene (A) activating its own repressor gene (B). In the Danino oscillator there is a latent phase when both the activator and repressor accumulate, thus allowing large amplitudes to be achieved for this configuration. The qualitative behavior of this oscillator is discussed in the modelling section using the mathematical model suggested in [23]. The advantage of the Danino oscillator topology over this amplified negative feedback oscillator is that it uses the small molecule AHL coupled to a constitutively expressed LuxR complex as the activator. This AHL molecule can diffuse across the cell membrane of various cells, where it can complex with the intracellular LuxR, and activate transcription inside the cell. Thus, this provides a means for cell-cell communication, which allows for synchronization of oscillations in small populations of cells over time.
Thus, we often cannot have two circuits with completely different functionalities present in the same cell, separated from one another and acting independently. Thus, for combining two circuits in one cell with switchability between the two when we desire, we need the circuits to be independent of each other in the sense that the repressors/activators for one system should not interact with the circuitry of the other. This is what is called orthogonality in synthetic biology systems [25].
In particular, consider the case of the repressilator [21] and the genetic toggle switch [15]. The repressilator used the only three classes of repressors that have been properly characterised (LacI, TetR and cI), while the genetic bi-stable toggle switch also used two of the repressors from the same list. Thus, it is not possible to add the two systems in a single cell in tandem and expect them to work independently of each other. This is due to the lack of availability of orthogonal parts that are well defined, characterized, and largely similar in functional form to the one desired or predicted from the model [26].
While there are significant advances being made and libraries of orthogonal parts for use in synthetic genetic circuits are being made [1], the use of these parts remains limited, and there is a need for combining the two circuits in a manner that while working with one aspect of the system in the circuit, the other one remains unaffected. Also there is a need to combine these circuits in such a manner that we can channel between different circuits as required. From here arises the concept of reconfigurability in biological circuits, ie circuits with switchable behaviour. These systems would be the ideal multi-taskers, switching between functionalities as the need arises. Synthetic Biology has seen the creation of certain circuits in this domain; from reconfigurable logic gate systems [23,27] [28] to systems with switchable dynamical behavior [27,29,30]. Works like [29,30] that propose a novel synthetic circuit, which can act both as an oscillator and a toggle switch have been a major influence on our work. Inspired from reconfigurable computer architecture that has been done by processing with very flexible high speed computing fabrics like field programmable gate arrays (FPGAs), genetic reconfigurable circuits could provide the solution to the limitation of using multiple circuits consisting of non-orthogonal parts, and also reduce the size of circuits, enabling multiple functionality from a single reconfigurable circuit.
Synthetic biology is at the intersection of biology, engineering and computational mathematics. Consequently, there are two distinct aspects to the design, analysis and implementation of a synthetic circuit: in-silico and in-vivo validation. Under the in-silico paradigm, we forego the specific implementation details, and focus on the abstract mathematical model and the dynamical behavior such a model allows. While for in-vivo analysis, we work on identifying the specific biological components that can be used to implement the identified topology for the synthetic network. Here, in order to create reconfigurable circuits in bacteria (E.coli), we ran simulations and provided a framework and proposed a circuit for the creation of reconfigurable networks that allow for switching of behavior dependent on the input given to the cells. We use optogenetics to implement the required behavior. Thus, the system implements an oscillator with a light dependent switch, which can work in one of two ways. The native system is a synchronized two component oscillator, which, dependent on the frequency of light shone on it, can be switched to a system which shuts off the oscillations completely, sending the system to a state of constitutive expression of A and no expression of B (A and B being the two components of the oscillator), or can create a system in which the frequency of the oscillations changes. This module is titled the 'Highly Optogenetically Tunable Frequency Modulator (HOTFM)'. Oscillators can be used as clock references, and thus, controlling the frequency of oscillation would offer certain advantages. For instance, increasing the frequency might allow faster execution of a process and vice versa. In this regard, one possible medical application could be time-controlled drug delivery. If a drug is known to function optimally if delivered at distinct regular intervals, an oscillator handling the delivery of the drug would be an attractive option. Further, if under circumstances quite stark from the ordinary, the rate of delivery needs to be increased (decreased), our approach could be a viable option. We require two additions to an oscillator to realize HOTFM. A mechanism to modulate the frequency of the oscillations and an optogenetic module to act as a switch for this mechanism.
We lay the computational framework here for the creation of such circuits, and its scale up, which would allow switching to an even wider range of responses in complex systems.

Danino oscillator architecture
The Danino oscillator [23] is a two component oscillator consisting of an activator and a repressor, LuxI and AiiA respectively, Fig 1. LuxI produces small molecules called AHL which bind to the LuxR protein to form a transcriptional activator complex. This complex activates LuxI and AiiA, both of which are placed downstream of the Lux promoter. Whereas, AiiA The topology remains largely the same, except that B undergoes negative autoregulation via self repression. It is suggested that here that this self-repression be mediated by optogenetics, using two component light response regulation systems such as [31,32].
https://doi.org/10.1371/journal.pone.0183242.g001 encodes a lactonase which degrades the AHL molecules, thereby leading to an indirect repression of LuxI. The novelty of the Danino oscillator lies in its ability to generate synchronized oscillations in a colony of cells. This happens via AHL's capability to diffuse across the cell membrane and be exchanged between cells akin to a communication signal for synchronizing the oscillators in each cell.
The oscillations in the Danino oscillator can be divided into three phases: AHL buildup, AHL degradation and the resting phases. These phases can be seen in Fig 2. During the buildup phase AHL molecules accumulate in the system, and once the concentration reaches sufficient levels, there is a burst in the concentrations of both LuxI and AiiA. Concurrently, the increased AiiA values begin to degrade the accumulated AHL molecules. After the AHL levels have been brought down to a resting state, both LuxI and AiiA concentrations start to fall, owing to enzymatic degradation; following which we again have the buildup phase. In this manner the oscillations continue.

Danino oscillator: Mathematical model
The Danino oscillator is modelled using a set of delay differential equations consisiting of four equations, one for each AiiA, LuxI, internal AHL and external AHL concentrations [23]. We recapitulate the model for the Danino oscillator here from Eqs 1-4. This is a deterministic model for intracellular concentrations of LuxI(I), AiiA(A), internal AHL(H i ) and external In Eqs 1 and 2 the term accounts for the time-delayed effect of AHL on AiiA and LuxI through the concentration of the internal AHL H τ (t) = H i (t-τ). Description for the various parameters in Eqs 1-4 can be found in [23]. A summary of the parameters is given as follows:

HOTFM: Architecture
The Danino oscillator has a linear relationship between amplitude of oscillations and the time period [23], owing to the fact that the time period is dictated by the time it takes for the AiiA and LuxI concentrations to degrade. One possible way to bring down the time period of oscillations would be to reduce the length of this degradation period. Thus, we could, in principle, modulate the amplitude to modulate the frequency of oscillation. We propose to do this by adding negative autoregulation on AiiA Fig 1. This transforms the amplified negative feedback oscillator topology of the Danino oscillator into a smolen oscillator [19]. The Smolen oscillator is known to have scope for tunability of the frequency response [19,33].
To control the negative regulation of AiiA we need a reliable way to switch on and off this pathway. For this we use optogentics as a flexible and easy mechanism. Optogenetics is a method that is quickly gaining traction for regulation of gene expression. This is due to the fact that such systems offer dynamic switching between the active and inactive state (light systems typically switch from ON to OFF state in the order of seconds to minutes) [34], and the input can be varied instantaneously by varying the intensity of light as opposed to chemical systems where a difference in chemical inducer concentration requires time to be diffused into the system and achieve a uniform concentration gradient. Thus, optogenetic systems have found applications in several fields such as in silico feedback [35,36] and control [31,37], understanding cell signaling pathways [38] and neuroscience [39] The self-repression pathway in our proposed system thus, is put under the control of a repressor, such as TetR; and this repressor is controlled by a two component light response system [31], which is an optogenetic module. So, whenever the oscillator is required to change its frequency of operation, we can, in principle, shine light on the cells and the self-repression pathway (mediated by the repressor) would be activated. Further, in order to ensure that the repression is a self-repression, i.e. in proportion to the amount of AiiA present, we have the two genes of the two component system placed under the same controls as that of AiiA, with the promoter and degradation tag used for both AiiA and the two component system being the same.

HOTFM: Mathematical modelling
The model for HOTFM is based on that of the Danino oscillator, Eqs 1-4 [23]. To incorporate optogenetically triggered negative autoregulation on AiiA, we add the following two equations to the set of differential equations for the Danino oscillator: and modify Eq 1, where,  4 ] accounts for the limitation in the protein production owing to nutrient limitation and waste accumulation [23]. Both the two component system and AiiA are repressed by the repressor Eqs 6 and 7 via hill type repression.

HOTFM: Proposed biological implementation
The core of the HOTFM consists of an architecture similar to the Danino oscillator [23], with the two components being LuxI (A) and AiiA (B), that function in the same manner as described in Danino Oscillator Architecture. In this system, however, the LuxI gene is under the Lux promoter, while the AiiA gene is under a hybrid Lux-Tet promoter, which is activated by AHL but repressed by a repressor such as TetR. This modification has been added in order to facilitate an auto self-repression step on the AiiA component. To achieve this, we propose to express two two component light response systems, namely the CcaS-CcaR and the UirS-UirR systems [31,32], both under the control of hybrid Lux-Tet promoters, in order to ensure that at any point in time, the level of expression of these light response systems is identical to the expression of AiiA, which is crucial for self-repression. Additionally, it may be noted, that these two systems are orthogonal to each other as the wavelength of light for activation for the two systems are different, with the CcaS-CcaR system being activated under green light and deactivated under red light, while the UirS-UirR is activated under UV-violet light and also turned off under green light. In this cascade, we also express two copies of a repressor, such as TetR, one under the CPCG2 promoter driven by phosphorylated CcaR and the other under the csiR1 promoter, driven by phosphorylated UirR. Oscillations are monitored by a sfGFP reporter, connected downstream of a hybrid Lux-Tet promoter, which measures AiiA levels.
The circuit diagram is shown in Fig 3. In the figure, TC1 and TC2 represent the genes of the two component light response systems (TC), of which TC2 ultimately acts as the activator. The other symbols have their usual meaning.
The proposed system then, would work as follows-In the native dark state, when no light is being shone on the system, the cells would have a behavior similar to that described by [23], with a burst of transcription leading to a buildup of luxI, AHL and AiiA, which would quench the AHL and bring the system to its native state and repeat, driving the oscillations. In our system, the same behavior occurs, but along with this, the CcaS, CcaR, UirS and UirR proteins are being produced. These proteins however, are in an inactivated state, due to a lack of activating light. From here, our system can be switched into one of two states. If green light of wavelength 520nm is shone upon the system, it drives a conformational change in the membrane bound sensory kinase CcaS, which then phosphorylates the response regulator CcaR, which can now activate transcription of genes placed under the CPCG2 promoter. This CcaR therefore, activates the production of the repressor, which represses the production of aiiA transcriptionally. This regulation, coupled with AHL degradation that is catalyzed by AiiA bring down the levels of AHL faster than the native oscillator (without the light response system), thereby not allowing the oscillations to go up as high as the former case, bringing down the mean amplitude and consequentially increasing the frequency of the oscillations. This behavior is in accordance with the simulations by [23], which showed that the oscillator had a linear relationship between amplitude and time period. The other state of the HOTFM can be driven by shining UV-Violet light of wavelength 405 nm, which would activate the UirS-UirR system in a similar manner as the CcaS-CcaR system, leading to activation of the repressor under the csiR1 promoter. However, since the dynamics of the UirS-UirR system are different from that of the CcaS-CcaR system, the same repressor would have a different effect on the oscillations Results and Discussions. In this case, due to a shift in the dynamics of the light response system, the repressor has a higher activation under the same intensity of light, and thus represses all the AiiA, leading to a stop in the oscillations. This sends AiiA to a constitutive OFF state, while the LuxI now is free of repression and can self-amplify itself and go to a constitutive ON state, thus killing the oscillations.

Synchronization
Synchronization is one of the properties central to the Danino oscillator [23], as it presented the first biological oscillatory circuit capable of synchronizing over a small colony of around 200 cells. Previously developed oscillators [21,33] had single cells oscillating in either mircrofluidic chambers or tailored environments without synchronization. There are many concepts of synchronization, such as amplitude and frequency synchronization. The Danino oscillator presented a case of phase synchronization, wherein the oscillations among different cells, are present initially in different phases (due to stochasticity). Over time, the oscillators synchronize, due to coupling of oscillators, when they align to a common phase in their oscillations spontaneously [40]. This general definition of phase synchronization holds true across all domains, and has applications over vast areas of physics, biology and other engineering fields. The system discussed in this paper in particular, the Danino oscillator, from which the HOTFM architecture has been inspired, this coupling is achieved via a quorum sensing molecule called AHL. The AHL, produced in the system due to the presence of the LuxI gene in the circuit, diffuses across the cell membrane, entering neighboring cells, leading to coupling. This diffusion of AHL appears in the model as a partial second order derivative in the equation for external AHL, governed by Fick's second law of diffusion. This is a unique feature bestowed upon the cells only due to the oscillator circuit being added to them, as the E.coli cells naturally do not produce AHL or any molecule which could act as a means to couple the cells. When the oscillator circuit is added to them, AHL is produced due to the transcription and translation of the LuxI gene which is critical for oscillations, and produces AHL leading to coupling. This circuit, thus, when transformed into any cell would cause AHL generation and coupling, which is the key factor for synchronization.

Robustness
Robustness is an important concern in synthetic genetic circuits. As the field of synthetic biology grows, the number of circuits being proposed for a variety of applications is rising exponentially. However, a major difficulty associated with the development of new circuits is that they do not function as expected due to uncertain initial conditions and disturbances in the environment, which can cause the host cell to not produce the desired response [41]. Analysis of robustness, thus, is essential for the design of synthetic circuits [42][43][44]. There are several definitions of robustness, presented with different aspects in mind [45][46][47]. The basic idea behind robustness is the property of a system to be resistant and maintain its state against perturbations [48]. Here, we stick to a narrow definition of robustness, ie transient robustness. In order to achieve this, we study the regions in the parameter space of degradation rate of LuxI and AiiA, and also LuxI and AHL, in which oscillations are obtained for both the Danino oscillator and the HOTFM. In this region we observe the rate of synchronization in the two oscillations, and note down which oscillator synchronizes faster and is therefore more transiently robust over the entire parameter space. Additionaly, we look at the integrated absolute synchronization error (IASE) for quantifying the extent of deviation from synchronization; The oscillator with a lower IASE would be more robust Synchronization: HOTFM is more Robust compared to the Danino Oscillator.

Assessment of the proposed synthetic circuits
The two main hypotheses of this work are concerned with tunability of frequency of oscillations and robustness of synchronization. The former is a property at the cell level whereas the latter concerns the collective behavior of a colony of cells. Thus, we analyze the system in both these scenarios.
Single cell study. To conduct insilico assessment of the proposed synthetic circuit 'bulk' simulations are performed by dropping the diffusion term D 1 in Eq 4 as suggested in [23]. The resulting dynamical model for HOTFM Eqs 2-4, 5-7 is simulated with respect to the parameter of the repressor hill function K introduced in Eq 7 and the degradation coefficient for the repressor γ introduced in Eq 5. Numerical Integration of a system of differential equations is sensitive to the initial conditions assumed for the involved variables; and only converge to the stable solutions. Thus, it is not conducive for eliciting the qualitative dependence of the solution on the parameters of the system. In this work, numerical bifurcation analysis has been employed for this purpose; the change in behavior of the HOTFM oscillator is tracked in the two parameter space of K and γ using numerical bifurcation performed with the Matlab package DDE-BIFTOOL [49]. Further, numerical intergation is performed with the Matlab solver dde23 for obtaining representative solutions in the different regions identified via numerical bifurcaton.
Colony level study. Synchronization is a population level behavior, thus for studying robustness of synchronization a linear colony of 6 cells is assumed for both HOTFM and the Danino oscillators. For each cell in a colony a dynamical model is formulated using Eqs 1-4 for the Danino oscillator and Eqs 2-4, 5-7 for HOTFM. The diffusion term D 1 in Eq 4 is retained for this experiment since this parameter accounts for the coupling between the cells in a colony. Two two-parameter numerical bifurcations are conducted for both HOTFM and the Danino oscillators. The parameter spaces for this purpose are γ I -γ A and γ H -γ A ; where, γ A , γ I and γ H represent the enzymatic degradation for AiiA, LuxI and internal AHL respectively. For each parameter space, transient synchronization behavior is studied in the region of parameters where oscillations are feasible; further, for a comparative analyis, the region of intersection for HOTFM and the Danino oscillators is considered.
A narrow view of robustness, limited to transient behavior of synchronization, has been adopted in this work Robustness. To quantify this behavior two metrics to judge the level of synchronization in a colony of cells is used. These metrics are defined as follows: We also apply periodic boundary conditions. The n systems of differential equations are numerically integrated for 500 minutes. We define ASE as Where A i represents the concentration of AiiA at the ith cell at time t. ASE captures the absolute deviation of the colony of cells from synchronization at time t. To calculate an estimate of overall deviation we integrate ASE over the time span of numerical intergration t final , which is 500 minutes in this case, to give Integrated Absolute Sychronization Error (IASE).
The IASE values are compared between HOTFM and the Danino oscillators for ascertaining the robustness behavior for synchronization. The difference IASE Danino − IASE HOTFM is used as the comparative metric; a positive (negative) value indicating greater (lesser) robustness for HOTFM compared to the Danino oscillator.
2. Rate of Synchronizaton (r)-In addition to total deviation from synchronization, we also measure the rate at which the effect of noise introduced at time t = 0 on synchronization dies down. Equivalently, this rate gives the rate of synchronization (r). We fit an exponential Envelope to the ASE (EASE) and use the exponent of this fitted exponential envelope as an estimate of r.
The r values are compared between HOTFM and the Danino oscillators for ascertaining the robustness behavior for synchronization as well. The difference r HOTFM − r Danino is used as the comparative metric; a positive (negative) value indicating greater (lesser) robustness for HOTFM compared to the Danino oscillator. Fig 4 shows a representative example of ASE and r for HOTFM and the Danino oscillators. In these experiments the values for the parameters are adapted from [23], except for the free parameters. For the equations unique to HOTFM, Eqs 5 and 6, the hill function parameters for the two component system are adapted from [31] and [32] for CcaS-CcaR and UirS-UirR, respectively. For the study of robustness of synchronization the diffusion constant D 1 is set to a value of 100 μm 2 /s. The Matlab scripts for the bifurcation and time domain analyses have been provided as supporting information in S1 Code. Fig 5 shows the two-parameter bifurcation for HOTFM. We see two qualitatively different stable behaviors here-stable steady state and stable oscillations, separated by a branch of hopf solutions. The region below the hopf branch corresponds to the parameter values for which there exists a stable steady state; and the region above the hopf branch contains region for oscillatory behavior. It is evident that the branch of hopf solutions follows a hyperoblic trend: as the strength of the repression coefficient increases for the repressor, the degradation rate for which the system transitions from a stable steady state solution to a stable oscillatory solution decreases. Low value for the repression coefficient suggests a strong repressor, which has the potential to kill the oscillations at the AiiA node. Intuitively, this behavior can be counteracted by employing strong degradation tags that afford high degradation rates [50], thus allowing the oscillations to continue. Conversely, for weak repressors the repression coefficients take on large values, thus even weak degradation tags would suffice to ensure sustenance of oscillations; for weak repressors the system becomes equivalent to the Danino oscillator since the hill function component due to the repressor tends to unity for large values of the repression coefficient. The bifurcation plots of Fig 5 can be used to select the feasible values for the degradation tag for which a given repressor, such as TetR, would kill the oscillations.

HOTFM modulates the frequency and amplitude of oscillations
As hypothesized in HOTFM: Architecture, simulations confirmed the modulatory effect of adding self repression on AiiA. In regions farther away from the hopf branch within the oscillatory region, both the period and amplitude of oscillations increases Fig 6, finally tending to behaviors equivalent to the Danino oscillator. However, the period overshoots above the Danino oscillator for a while, and finally settles into the Danino oscillator; this behavior is not present in the amplitude heatmap. For illustrative purposes, we select the TetR repressor for which the repression coefficient K is 100, after suitable scaling [51].
A strategy for tuning the response of the Danino oscillator using self-repression can be suggested from Fig 6a-6c. For a given repressor, this can be achieved by varying the strength of the degradation tag. A low value for the degradation tag would temporarily kill the oscillations when appropriate light is shone on the system; such an arrangement can be employed to design an optogenetically controlled kill switch for oscillations. Once the hopf branch is crossed, oscillations survive for all values of the degrdation tag henceforth. Within this region amplitude can be traded-off for higher frequeny of oscillation (lower period of oscillation). Staying closer to the hopf branch would afford highest frequencies of operation at the expense Self-repression on AiiA as a means to tune the frequency response was conjectured from the linear relationship between the period and amplitude of oscillatons for the Danino oscillator [23]. The heatmaps in Fig 5 confirm positive correlation between the period and amplitude of oscillation for HOTFM. Fig 6c shows the variaton of the period and amplitude of oscillation along the branch of periodic solutions emanating from the hopf branch for the TetR repressor as the strength of the degradation tag is varied. Fig 6a and 6b again depict the increasing behavior of the period and amplitude of oscillations with increase in strength of the degradation tag. The linear relationship between the period and amplitude of oscillation is confirmed from Fig 6c for the TetR repressor.
A greater reconfigurability can be achieved by combining both the kill switch and frequency modulating behavior into a single circuit. As discussed in HOTFM: Proposed Biological Implementation, we can use both the CcaS-CcaR and the UirS-UirR systems in a single

Fig 5. Two parameter bifurcation for HOTFM in the repression coefficient (K)-degradation cefficient for repressor (γ) space. Two parameter bifurcation for HOTFM against the parameters for the repressor downstream of the two component light system (CcaS-CcaR)
; the x-axis corresponds to the repression coefficient K and y-axis corresponds to the degradation coefficient γ. The boundary of the shaded region represents the branch of hopf solutions; the unshaded region below the hopf branch corresponds to parameter values for which there exits a stable steady state and oscillations are absent. Whereas, the shaded region above the hopf branch represents parameter values for which stable oscillations exist. a) The two parameter bifurcation plot with the shaded region representing a heatmap for the period of oscillation. At each point in this region, the color represents the period of oscillation. Points farther away from the hopf branch in the shaded region show higher period (lower frequency), while the points closer to the hopf branch exhibit lower period (higher frequency). b) The two parameter bifurcation plot with the shaded region representing a heatmap for the period of oscillation. At each point in this region, the color represents the period of oscillation. Points farther away from the hopf branch in the shaded region show higher amplitude, while the points closer to the hopf branch exhibit lower amplitude.
https://doi.org/10.1371/journal.pone.0183242.g005 configuration by placing two different copies of the TetR repressor downstream to the CPCG2 and csiR1 promoters respectively. The repression coefficient of TetR repressor is around 100 [50]. The dynamics of the UirS-UirR systems is qualitatively similar to the Ccas-Ccar system [31,32]. This suggests that similar to the Ccas-CcaR system we can modulate the frequeency response using a repressor with varying strength of the degradation tags. So, if we use a degradation rate of 1.5 with the TetR repressor associated to the CcaS-CcaR system and that with an appropriate value associated to the UirS-UirR, we can get a system that can exhibit three different types of response. In the absence of any light source, the system will oscillate like a normal Danino oscillator; in the presence of green light the Ccar-Ccas will get activated and as can be seen from Fig 6, the system will oscillate at a frequency higher than the Danino oscillator. If now the Ccas-Ccar system is switched off by shining red light on the system followed by Ultraviolet light, we see from (refer supplementary) that the oscillations will shut down.

Synchronization: HOTFM is more robust compared to the Danino oscillator
As discussed in Assessment of the proposed Synthetic Circuits, two two-parameter bifurcations have been looked at for the robustness analysis. Figs 8,9 and Figs 10,11 show the bifurcation plots for HOTFM and the Danino oscillators in the γ I -γ A space and γ H -γ A spaces respectively. The regions above the hopf branches correspond to the parameter ranges for which oscillations are possible. For a comparative analysis, the oscillatory region common to both HOTFM and the Danino oscillators have been considered. Tables 1, 2 and Tables 3, 4 show the aggregate robustness metrics for the bifurcation analysis for HOTFM and the Danino oscillators in the γ I -γ A space and γ H -γ A spaces respectively.
Robustness analysis for IASE. It is evident from Figs 8, 10 and Tables 1, 3 that on an aggregate basis HOTFM is more robust compared to the Danino oscillator. Among the 520 (1510) points explored in the γ I -γ A (γ H -γ A ) parameter space, for 434 (1070) points (IASE Danino − IASE HOTFM ) > 0 and HOTFM is more robust, while for 86 (440) points (IASE Danino − IASE HOTFM ) < 0 and the Danino oscillator is more robust. To gauge the extent to which these estimates differ from randomly assigning robustness attribute (whether HOTFM is more robust or Danino) to the 520 (1510) points, boostrap resampling was conducted to generate the distributions as shown in Fig 12. The corresponding confidence intervals and effect sizes (cohen's d) [52] are displayed in Tables 1 and 3. The 17.25 (16.91) value of cohen's d suggests an extremely large deviation from the case of random assignment [52].
Boostrap estimates can also be generated for individual points in the parameter space

Robustness analysis for Rate of Synchronization (r). It is evident from Figs 9, 11
and Tables 2, 4 that on an aggregate basis HOTFM is more robust compared to the Danino oscillator for synchronization rate as well. Among the 520 (1510) points explored in the γ I -γ A (γ H -γ A ) parameter space, for 367 (1008) points (r Danino − r HOTFM ) < 0 and HOTFM is more robust, while for 153 (502) points (r Danino − r HOTFM ) > 0 and the Danino oscillator is more robust. To gauge the extent to which these estimates differ from randomly assigning robustness attribute (whether HOTFM is more robust or Danino) to the 520 (1510) points, boostrap resampling was conducted to generate the distributions as shown in Fig 12. The corresponding confidence intervals and effect sizes (cohen's d) [52] are displayed in Tables 2 and 4. The 9.72 (13.39) value of cohen's d suggests an extremely large deviation from the case of random assignment [52].
Similar to IASE, boostrap estimates can again be generated for individual points in the parameter space Fig 13; for a given point in the parameter space, r values for the cells in the   which the first 100 cells have the Danino circuitry, while the other 100 have the HOTFM circuit. In this experiment, the values for the parameters are adapted from [23]; for the equations unique to HOTFM, Eqs 5 and 6, the hill function parameters for the two component system are adapted from [31] for CcaS-CcaR; TetR is taken as the repressor downstream of Ccas-CcaR and the degradation rate for TetR is assumed to be 1.5. Beats were generated when the concentration of the AHL at the interface of cell numbers 100 and 101 was plotted with time. Also, had there been an AiiA molecule at the interface, the time series for that has been plotted as well. This also shows beat like patterns. This could be achieved biologically by placing membranes between cells number 100 and 101, such that the cells would not mix, and therefore the 100 cells of Danino would not diffuse into the HOTFM area and vice versa, while the AHL molecules could diffuse through it. In this section bound by two membranes, we would observe the pattern of beats in the AiiA and AHL levels S5 Video.

Conclusion
Here we have computationally demonstrated a framework for reconfigurable circuits, by providing the simulations for a system which can generate responses allowing us to switch between three qualitative behaviors. The first behavior is of the native Danino oscillator, while the second is of a frequency modulation of the oscillator (HOTFM), and the third is one that drives out the oscillations and forces the system to a steady state constant response. We have chosen light directed self-repression due to the noninvasive, fast acting and orthogonality of optogenetic systems [34,53,54]. We have also shown that, by modulating the strength of the degradation tag used [50], we can transform the output response obtained from the system for the same repressor, which also makes sense intuitively, as increasing the degradation rate of the repressor would require a higher strength of repressor to achieve the response that would ideally be done by a weaker repressor. This work is novel, as frequency modulation being shown for the first time in a biological oscillator. Another quality is that the system has cellcell communication and is therefore synchronised, which is a highly desirable quality when we begin to talk of real world applications of these circuits in fields such as that of time dependent targeted drug delivery, as we want to observe phenomenon at a bulk level. The requirement of synchronisation has seen recent developments in previous oscillators such as the repressilator achieving synchronisation as well [55]. Further, HOTFM exhibits more robust behavior in response to input perturbations with respect to synchronization in the transient phase.
Thus, here we have laid the framework for future work to be done in the field, which could have long range applications in the area of time dependent drug delivery, as the field advances and grows. In this case, the cells start to synchronise, due to which we see a peak in the frequency distribution, which oscillates with time. (MP4) S5 Video. Generation of beats. Representation of beats in time. The x-axis corresponds to the cell number and the color shows the concentration for AiiA. Blue shows low concentration and yellow high concentrations. The 100 cells on the left belong to HOTFM while the remaining 100 on the right to Danino. Initially, only middle cells for each colony have non-zero concentration and as time progresses, diffusion of AHL spreads the oscillation throughout the mix. We get beats at the interface of the two colonies. (PNG) S1 Code. Matlab scripts. The Matlab scripts used for conducting the numerical bifurcation and time intergation analyses. A Readme file with instructions on using these scripts has been included as well. (GZ) S1 File. Supporting information details. This file contains further details for the various supporting information figures and videos. Specifically, details regarding the synchronized behavior of the Danino oscillator, S1-S4 Videos, have been given. Further, the oscillation kill switch depicted in S1 Fig has been explained in some detail. (PDF)