Analysis and control of genetic toggle switches subject to periodic multi-input stimulation
AAnalysis and control of genetic toggle switchessubject to periodic multi-input stimulation
Davide Fiore , Agostino Guarino , Mario di Bernardo , Abstract —In this letter, we analyze a genetic toggle switchrecently studied in the literature where the expression of tworepressor proteins can be tuned by controlling two differentinputs, namely the concentration of two inducer molecules inthe growth medium of the cells. Specifically, we investigate thedynamics of this system when subject to pulse-width modulated(PWM) input. We provide an analytical model that captures qual-itatively the experimental observations reported in the literatureand approximates its asymptotic behavior. We also discuss theeffect that the system parameters have on the prediction accuracyof the model. Moreover, we propose a possible external controlstrategy to regulate the mean value of the fluorescence of thereporter proteins when the cells are subject to such periodicforcing.
Index Terms —Systems biology, Genetic regulatory systems,Modeling
I. I
NTRODUCTION T HE genetic toggle switch is a fundamental componentin synthetic biology as it plays a major role in celldifferentiation and decision making [1], [2]. Its importancecomes from its ability to endow host cells with memory ofsome previous stimulus reporting this information as highexpression rate of a specific repressor protein [3]–[5].The genetic toggle switch as first designed in [3] consists oftwo repressor proteins, both repressing each other’s promoter,so that only one protein is fully expressed at any time.From a modelling viewpoint, the genetic toggle switch is abistable dynamical system, possessing two stable equilibria,each associated to a fully expressed protein, and a saddleequilibrium point, whose stable manifold is the boundaryseparating the basins of attraction of the other two.Different approaches have been presented to control the re-sponse of genetic toggles switches. Examples include methodsbased on piecewise affine approximations [6], pulse shapingof the external inputs based on monotone systems theory [7],and the analysis of the stationary probability distributions ofthe outputs in different working conditions [8].Recently, in [9] the problem has been studied of dynamically“balancing” a genetic toggle switch (based on the LacI/TetRpromoters in E.coli, schematically shown in Figure 1) inan undecided state somewhere in between its two stableequilibrium points. The expression level of the two repressing Davide Fiore, Agostino Guarino and Mario di Bernardo are withthe Department of Electrical Engineering and Information Technology,University of Naples Federico II, Via Claudio 21, 80125 Naples, Italy. [email protected], [email protected] Mario di Bernardo is also with the Department of Engineering Math-ematics, University of Bristol, University Walk, BS8 1TR Bristol, U.K. [email protected]
Fig. 1. Genetic toggle switch embedded in
E.coli considered in [9] (Figurereused under Creative Commons license). proteins can be controlled by regulating the concentrationof two inducer molecules, aTc and IPTG. The former, aTc,binds to TetR, increasing the rate of production of LacI, andtherefore causing the cell to commit to the stable equilibriumpoint corresponding to high expression of LacI (high LacI/lowTetR). The latter, IPTG, binds instead to LacI, causing thecommitment of the cell to the other stable equilibrium point(high TetR/low LacI). From a dynamical systems viewpoint,varying the two input signals causes the occurrence of twosaddle-node bifurcations changing the phase portrait of thesystem from bistability to monostability (Figure 2).In their work Lugagne et al. [9] focus on both the problemof controlling a single cell and that of taming the behavior ofthe whole population. Their approach is based on consideringthe toggle switch as a multi-input control system and isaimed at using both inputs to keep the switch evolving ina neighborhood of its saddle point; a problem they proposeas a test-bed scenario in synthetic biology similar to that ofstabilizing an inverted pendulum in classical control.When implementing single cell control, the fluorescencelevel of the reporter proteins in a single cell are measured andcompared to their reference values. Two different classes ofcontrollers were used in [9], PI and bang-bang, both designedindependently for each control input (aTc and IPTG). Using PIcontrollers on both input channels, it is possible to make thesingle cell evolve (oscillate) near the saddle point. Althoughthe controlled cell follows (on average) the desired reference,the rest of the population is observed to drift away, converginginstead to some other equilibrium point.Surprisingly, it is reported in [9] that this undesired effect isabsent when the single cell is controlled by two independentbang-bang inputs with the rest of the population exhibitingan evolution similar to the target cell in this case. To furtherexplore this effect, the authors then consider an open-loop a r X i v : . [ q - b i o . M N ] S e p eriodic stimulation (two mutually exclusive pulse waves withprescribed width) to control the whole population. Again thewhole population is shown to converge to some periodicorbit surrounding the saddle point with a remarkable levelof coherence in terms of both mean and standard deviationdespite cell-to-cell variability and other phenotypic differencesbetween cells.Using an in-silico model this effect is explained in [9] asdue to the phase portrait of the forced system periodicallychanging from one presenting a unique high-LacI equilibriumpoint to another with a unique high-TetR equilibrium point.Heuristically, this results in an average phase portrait havinga unique attractor in between the former two given that,as conjectured in [9], the cell dynamics and the periodicexcitation act on different time-scales. Also, changing thecharacteristics of the periodic PWM forcing (such as period,width and amplitude of the pulses) shifts the position of theaverage attractor causing cells to evolve towards a differenttarget solution.Despite providing some qualitative explanation of the ex-perimental observations, several open questions remain. Forinstance, what causes the massive reduction in standard devi-ation between different cells in the population and what theperiod/duty cycle should be of the control inputs to achievethe desired behavior. Also, the challenge remains of designingbetter multi-input feedback strategies to control populations ofhost cells endowed with synthetic toggle switches.In this letter, we address some of these open problemsby providing an analytical investigation of the phenomenareported in [9]. We start by deriving a quasi-steady statemodel of the toggle-switch system proposed therein. Usingformal averaging techniques for nonlinear systems [10], wederive an autonomous average vector field , whose solutions,under some conditions, approximate those of the original time-varying system. To simplify the analysis, we assume that thediffusion of the inducer molecules across the cell membraneis instantaneous .We prove that if the average vector field has a uniqueattracting equilibrium point ¯ x av , whose position in state spacedepends on the duty cycle D and on the amplitude of theforcing pulse waves u aTc ( t ) and u IPTG ( t ) , then every solutionof the original time-varying system asymptotically convergesto a periodic orbit in some neighborhood of ¯ x av . We compareour model predictions with the experimental observationsmade in [9] and with the mean-value trajectories of the originalmodel proposed therein. We use the model and its analysisto provide some indications on how the parameters of thetoggle switch may be tuned to enhance its response to theclass of periodic inputs of interest, and exploit the results tosynthesize an external control strategy to regulate the mean-value of the measured fluorescence of the reporter proteinsin the cell at some desired value. We wish to emphasize thatthe analysis provided in this letter can be instrumental for thedesign of further control strategies for this particularly relevantclass of synthetic devices and to investigate the effects at thepopulation level of different types of periodic stimuli to thecells. T e t R Fig. 2. Nullclines of the toggle switch system (8). Main picture: bistability:two stable and one saddle equilibrium points. Reference values aT c =20 ng / ml , IP T G = 0 .
25 mM . Insets: a) monostability: unique highLacI/low TetR equilibrium point. aT c = 50 ng / ml , IP T G = 0 .
25 mM ;b) monostability: unique high TetR/low LacI equilibrium point. aT c =20 ng / ml , IP T G = 0 .
50 mM
II. M
ATHEMATICAL MODEL OF THE TOGGLE SWITCH
A. Transcription-translation model
The deterministic model of the toggle switch that we startfrom can be given as follows [9] d mRNA
LacI dt = κ m0L + κ mL (cid:16) TetRθ
TetR · aTc/θ aTC ) η aTc (cid:17) η TetR − g mL · mRNA LacI (1) d mRNA
TetR dt = κ m0T + κ mT (cid:16) LacIθ
LacI · IPTG/θ
IPTG ) η IPTG (cid:17) η LacI − g mT · mRNA TetR (2) d LacIdt = κ pL · mRNA LacI − g pL · LacI (3) d T etRdt = κ pT · mRNA TetR − g pT · T etR (4)
In the above equations the variables denote concentrations ofmolecules inside the cell, and the parameters κ m0L / T , κ mL / T , κ pL / T , g mL / T , g pL / T are leakage transcription, transcription,translation, mRNA degradation, and protein degradation rates,respectively. All parameter values are provided in Supplemen-tary Table 1 and are also the same used in [9].The inducer molecules diffuse in and out of the cell acrossthe membrane with non-symmetrical exchange dynamics mod-eled by d aT cdt = (cid:40) k inaTc ( u aTc − aT c ) , if u aTc > aT ck outaTc ( u aTc − aT c ) , if u aTc ≤ aT c , (5) d IP T Gdt = (cid:40) k inIPTG ( u IPTG − IP T G ) , if u IPTG > IP T Gk outIPTG ( u IPTG − IP T G ) , if u IPTG ≤ IP T G , (6) where aT c and
IP T G denote the concentrations of theinducer molecules inside the cell, while u aTc and u IPTG thosein the growth medium.
10 20 30 40 50aTc0.00.20.40.60.81.0 w w Fig. 3. Top: Static nonlinear functions w ( aT c ) and w ( IP T G ) as in (11)and (12). Bottom: Pulse wave s q ( t ) : period , duty cycle D ∈ [0 , . B. Quasi-steady state model
Assuming that the concentrations of the mRNA moleculesreach steady state more rapidly than their correspondingproteins, that LacI and TetR proteins degrade at the same rate,that is g pL = g pT = g p , and using the following dimensionlessvariables (similarly as done in [11], [12]) t (cid:48) = g p t, x = LacIθ
LacI , x = T etRθ
TetR , (7)we obtain the following nondimensional quasi-steady statemodel of the genetic toggle switch dx dt (cid:48) = k + k x · w ( t (cid:48) /g p ) − x dx dt (cid:48) = k + k x · w ( t (cid:48) /g p ) − x (8)where k = κ m0L κ pL g mL θ LacI g p , k = κ mL κ pL g mL θ LacI g p , (9)and k = κ m0T κ pT g mT θ TetR g p , k = κ mT κ pT g mT θ TetR g p , (10)are dimensionless parameters, and we have set η LacI = η TetR = 2 . The steps of the previous derivation are reportedin the Supplementary Material.The nonlinear functions w ( t ) and w ( t ) in (8) take intoaccount the static relationship between the repressor protein(TetR or LacI) and their regulator molecule (aTc or IPTG,respectively). They are shown in Figure 3 and are defined as w ( aT c ( t )) = 1 (cid:16) (cid:16) aT c ( t ) θ aTC (cid:17) η aTc (cid:17) η TetR (11) w ( IP T G ( t )) = 1 (cid:16) (cid:16) IP T G ( t ) θ IPTG (cid:17) η IPTG (cid:17) η LacI (12)System (8) with the static relations (11)-(12) and diffusiondynamics across the cell membrane (5)-(6) can be representedin block form as in Figure 4. The cell membrane acts asa linear (non-symmetrical) first order low-pass filter for thesignals u aTc ( t ) and u IPTG ( t ) with a cut-off frequency that Fig. 4. Block diagram of system (8) with diffusion dynamics across the cellmembrane (5)-(6). depends on the diffusion exchange rates k in / outaTc and k in / outIPTG .Hence, aT c ( t ) and IP T G ( t ) are filtered version of theirrespective input signals whose attenuation depends both onthe cut-off frequency and on their spectral density.In our analysis we make the following simplifying assump-tion. Assumption 1:
The diffusion dynamics of the inducermolecules, aTc and IPTG, across the cell membrane is in-stantaneous, that is aT c ( t ) = u aTc ( t ) , (13) IP T G ( t ) = u IPTG ( t ) , (14)for every t ≥ t .Later in Section IV, we will compare our results derivedfrom system (8) under the above Assumption 1 with thesolutions of the complete toggle switch model (1)-(4) withmore realistic diffusion dynamics given by (5)-(6).III. A VERAGING ANALYSIS OF THE TOGGLE SWITCHUNDER
PWM
INPUT SIGNALS
A. Forcing signals
Following [9], the concentrations of the inducers in thegrowth medium are selected as two mutually exclusive pulsewaves of period T , duty cycle D ∈ [0 , and amplitude ¯ u aTc and ¯ u IPTG , respectively, that is u aTc ( t ) = ¯ u aTc · (1 − s q ( t/T )) (15) u IPTG ( t ) = ¯ u IPTG · s q ( t/T ) (16)where s q ( t ) is the pulse wave taking values 0 and 1, withperiod and duty cycle D , reported in Figure 3. In theexperiments described in [9], the amplitude ¯ u aTc and ¯ u IPTG were allowed to take values between and
100 ng / ml , and and , respectively.Note that D = 0 corresponds to “high aTc/no IPTG” inthe growth medium which in turns results in full steady-stateexpression of LacI (high x ). Likewise, D = 1 correspondsto “no aTc/high IPTG” yielding full expression of TetR (high x ). Therefore, the duty cycle can be used to control the ratiobetween the activation time of the two monostable systemsassociated to the presence or absence of the two inducermolecules whose nullclines are shown in the insets in Figure2. Under Assumption 1 it follows that w ( t ) = w ( aT c ( t )) = w (¯ u aTc · (1 − s q ( t/T )))= ¯ w + (1 − ¯ w ) · s q ( t/T ) , (17)here ¯ w = w (¯ u aTc ) , and w ( t ) = w ( IP T G ( t )) = w (¯ u IPTG · s q ( t/T ))= 1 − (1 − ¯ w ) · s q ( t/T ) , (18)where ¯ w = w (¯ u IPTG ) . Therefore, w i ( t ) is a pulse wavetaking values between and ¯ w i . B. Average vector field
By rescaling time setting τ = t (cid:48) T g p , system (8) can be recastas dx dτ = ε (cid:20) k + k x · w ( τ T ) − x (cid:21) dx dτ = ε (cid:20) k + k x · w ( τ T ) − x (cid:21) (19)with ε = T g p . The vector field in (19) is time-varying in τ with period , and it is now in a form amenable for periodicaveraging analysis (see Supplementary Material).In particular, the average vector field, say f av ( x ) , can beobtained by integrating the vector field in (19) over a period,yielding f av , ( x ) = 11 (cid:90) (cid:18) k + k x · w ( τ T ) − x (cid:19) dτ = k + k (cid:32)(cid:90) D
11 + x · dτ + (cid:90) D
11 + x · ¯ w dτ (cid:33) − x = k + k (cid:18) D x + 1 − D x · ¯ w (cid:19) − x , where we used (17), and similarly for f av , ( x ) , f av , ( x ) = 11 (cid:90) (cid:18) k + k x · w ( τ T ) − x (cid:19) dτ = k + k (cid:32)(cid:90) D
11 + x · ¯ w dτ + (cid:90) D
11 + x · dτ (cid:33) − x = k + k (cid:18) D x · ¯ w + 1 − D x (cid:19) − x , where we used (18).Hence, the resulting average system is dx dτ = ε (cid:20) k + k (cid:18) D x + 1 − D x · ¯ w (cid:19) − x (cid:21) dx dτ = ε (cid:20) k + k (cid:18) D x · ¯ w + 1 − D x (cid:19) − x (cid:21) (20)Let x ( τ, ε ) and x av ( ετ ) denote the solutions to (19) and(20), respectively. Assume ¯ x av is an exponentially stable equi-librium point of the average system (20). Let Ω be a compactsubset of its basin of attraction, and assume x av (0) ∈ Ω , and x (0 , ε ) − x av (0) = O ( ε ) . Then, from [10, Theorem 10.4],there exists a positive parameter ε ∗ = T ∗ g p such that for all < ε < ε ∗ x ( τ, ε ) − x av ( ετ ) = O ( ε ) (21)for all τ > . That is, solutions x ( τ, ε ) to system (19)can be approximated by solutions x av ( ετ ) to (20) with anerror that is proportional to ε . As a consequence, if ¯ x av
500 1000 1500 2000 2500 3000 3500 LacI20040060080010001200TetR
Fig. 5. Equilibrium points ¯ x av of (20) as a function of duty cycle D rescaledin arbitrary fluorescence units using (7). Each dot represents the location ofthe unique stable equilibrium point of system (20) evaluated for D takingvalues in the interval [0 , with increments of . . is the unique equilibrium point of system (20), then for all < ε < ε ∗ system (19) has a unique, exponentially stable,periodic solution ¯ x ( τ, ε ) in a O ( ε ) -neighborhood of ¯ x av .The number and position in state space of the equilibriumpoints ¯ x av of the average system (20) depend on the specificchoice of the amplitudes ¯ u aTc and ¯ u IPTG of the pulse waves,and also on the value of the duty cycle D . For example, forthe reference values ¯ u aTc = 50 ng / ml and ¯ u IPTG = 0 . ,system (20) is monostable and the position of the equilibriumpoint ¯ x av varies monotonically with D as reported in Figure 5(blue dots). Hence, given certain values of ¯ u aTc and ¯ u IPTG , itis possible to move the position of ¯ x av on the correspondingcurve by varying D (Supplementary Figure S1).The phase portrait of the average system (20) together witha representative solution of the time-varying system (19) for D equal to . are depicted in Figure 6(a), while for D equal to . and . are reported in Supplementary FigureS2. The parameter ε has been set to . which corresponds toa forcing period T = ε/g p ≈ , and the system has beensimulated for t f = τ f T ≈ · . Larger valuesof ε correspond to larger values of the forcing period T . Inturn, from (21), this also implies that the solution x ( τ, ε ) of(19) will asymptotically converge to a periodic solution ¯ x ( τ, ε ) contained in a larger set (Figure 6(b)), and hence to a worseapproximation (see also Supplementary Figure S6 for theirtime evolution). IV. D IFFUSION EFFECTS
The analysis in the previous section was conducted un-der Assumption 1. As already mentioned before, the cellmembrane acts as a low-pass filter, hence, when Assumption1 is dropped, aT c ( t ) and IP T G ( t ) will not anymore beideal pulse waves but their filtered versions through the cellmembrane. Therefore, in order for the average system (20)to continue being a good approximation of the actual cellresponse, the cut-off frequency of the two low-pass filtersshould be sufficiently higher than the fundamental frequency
20 40 60 80 100 120010203040 x x (a) D = 0 . , T ≈ ( ε = 0 . ). x x (b) D = 0 . , T ≈
180 min ( ε = 3 ).Fig. 6. Background: phase portrait of the average system (20). Red line:the solution of the time-varying system (19) with ¯ u aTc = 50 ng / ml and ¯ u IPTG = 0 . from initial condition [1 , T . /T of the input pulse waves. However, due to the inevitableattenuation of high-frequency harmonics, there will always bea mismatch between the actual mean response of the cell andthe value predicted by (20).The effects of relaxing Assumption 1 on the time responseof system (19) can be observed in Supplementary FiguresS4-S5. The mean steady-state response of the complete four-dimensional system (1)-(4) with diffusion dynamics (5)-(6)is compared in Figure 7, and the corresponding equilibriumpoint ¯ x av ( D ) predicted by the autonomous two-dimensionalaverage system (20), for a representative value of the PWMamplitudes and different values of D (see SupplementaryFigure S3 for a different choice). Although as expected thereis no perfect matching between the two, the observed behavioris well captured by the average system. Note that in regulationproblems, this mismatch can be compensated by designing anadequate feedback action.When, on the other hand, the cut-off frequency of one ofthe filters is lower than the frequency /T of the input pulsewaves, the input signal will be highly attenuated, resulting in Fig. 7. Orange dots: Mean-value, evaluated at regime, of the response ofsystem (1)-(4) (with membrane dynamics (5)-(6)) to PWM inputs with T =240 min and varying D from . to . with increments of . . Bluedots: corresponding equilibrium point ¯ x av ( D ) of system (20) rescaled in a.u.using (7). Amplitude of pulse waves set to ¯ u aTc = 50 ng / ml and ¯ u IPTG =0 . . the simple regulation of the toggle switch to either one of thestable equilibrium points (a phenomenon that was reported inthe experiments described in [9, Supplementary Figure 8]). Asimilar phenomenon can also occur when the duty cycle isclose to or . Indeed, close to these values, the amplitude ofthe harmonics of the pulse wave is | a n | = (cid:12)(cid:12) un π sin( nπD ) (cid:12)(cid:12) ≈ uD , therefore low-frequency harmonics will have amplitudessimilar to those of high-frequency ones, and the pulse wavewill be highly attenuated.V. P ERSPECTIVES FOR CONTROL
We wish to emphasize that the analytical results derivedhere can be exploited for the synthesis of external controllersto regulate the mean-value of the output response of thegenetic toggle switch. Specifically, we propose the controlschematic shown in Figure 8 which is currently under de-velopment and will be presented elsewhere. Indeed, as donein Figure 5, it is possible to numerically compute ¯ x av as afunction of ¯ u aTc , ¯ u IPTG and D , and get interpolating curves Γ ¯ u aTc , ¯ u IPTG ( D ) . From these one can then obtain, for givenvalues of the amplitude ¯ u aTc and ¯ u IPTG , the duty cycle D ref corresponding to the desired average set-point ¯ x refav , thatis D ref = Γ − u aTc , ¯ u IPTG (¯ x refav ) . The mismatch e between themeasured mean-value of the plant outputs and ¯ x refav is thenprojected by π onto the curve Γ ¯ u aTc , ¯ u IPTG and compensatedby a PI controller. The control scheme should also take intoaccount the effects of the sampling time and of the slowtransients. VI. C
ONCLUSIONS
We derived and analyzed a model to capture the responseof the genetic toggle switch to mutually exclusive PWMinputs observed experimentally in [9]. The analysis was basedon the assumption that the diffusion of inducer moleculesacross the cell membrane is instantaneous. From this, usingthe periodic averaging method for nonlinear systems, we ig. 8. External controller for the regulation of the mean-response of a genetictoggle switch. derived an autonomous vector field that describes the dynamicsof the mean-value of the periodic solutions of the originalsystem. After discussing the predictions of the model underthe assumption of instantaneous diffusion, we relaxed thisassumption so that the input signals become filtered versionsof themselves worsening the predictions.However, even if it is not possible to eliminate the atten-uation due to the cell membrane, our analysis shows that tomitigate its effects the frequency /T of the input pulse wavesshould be chosen sufficiently lower than the cut-off frequencyof the low-pass membrane filter, and extreme values of theduty cycle D should be avoided. At the same time, we find thatto avoid large oscillations around ¯ x av , the parameter ε = T g p ,that is the ratio between the time-scales of the forcing inputsand system dynamics, should be taken as small as possible,e.g., for fixed T , by cooling down the temperature of thegrowth medium and thus reducing the cell growth rate andtherefore g p .Future work will be aimed at quantifying the effects of theattenuation of the input signals due to the cell membrane toimprove the predictions of our model, and at implementingand validating (in-silico and in-vivo) external controllers,also capable of modulating the ON/OFF values of the pulsewaves. Furthermore, we also plan to investigate the effect thatdifferent classes of periodic forcing could have on the varianceof the response of a population of cells with extrinsic noise.ACKNOWLEDGMENT The authors wish to acknowledge support from the research project COSY-BIO (Control Engineering of Biological Systems for Reliable SyntheticBiology Applications) funded by the European Union’s Horizon 2020 researchand innovation programme under grant agreement No 766840. R EFERENCES[1] U. Alon,
An introduction to systems biology: design principles ofbiological circuits . CRC press, 2006.[2] L. Chen, R. Wang, C. Li, and K. Aihara,
Modeling biomolecularnetworks in cells . Springer-Verlag London, 2010.[3] T. S. Gardner, C. R. Cantor, and J. J. Collins, “Construction of a genetictoggle switch in Escherichia coli,”
Nature , vol. 403, no. 6767, pp. 339–342, 2000.[4] T. Tian and K. Burrage, “Stochastic models for regulatory networksof the genetic toggle switch,”
Proceedings of the National Academy ofSciences , vol. 103, no. 22, pp. 8372–8377, 2006.[5] M. Wu, R.-Q. Su, X. Li, T. Ellis, Y.-C. Lai, and X. Wang, “Engineeringof regulated stochastic cell fate determination,”
Proceedings of theNational Academy of Sciences , vol. 110, no. 26, pp. 10 610–10 615,2013.[6] M. Chaves and J.-L. Gouz´e, “Exact control of genetic networks in aqualitative framework: the bistable switch example,”
Automatica , vol. 47,no. 6, pp. 1105–1112, 2011. [7] A. Sootla, D. Oyarz´un, D. Angeli, and G.-B. Stan, “Shaping pulses tocontrol bistable systems: Analysis, computation and counterexamples,”
Automatica , vol. 63, pp. 254–264, 2016.[8] A. Petrides and G. Vinnicombe, “Understanding the discrete genetictoggle switch phenomena using a discrete ‘nullcline’ construct inspiredby the Markov chain tree theorem,” in
Proc. of the 56th IEEE Conferenceon Decision and Control , 2017, pp. 1614–1621.[9] J.-B. Lugagne, S. S. Carrillo, M. Kirch, A. K¨ohler, G. Batt, andP. Hersen, “Balancing a genetic toggle switch by real-time feedbackcontrol and periodic forcing,”
Nature Communications , vol. 8, no. 1, p.1671, 2017.[10] H. K. Khalil,
Nonlinear systems , 3rd ed. Prentice Hall, 2002.[11] A. Kuznetsov, M. Kærn, and N. Kopell, “Synchrony in a population ofhysteresis-based genetic oscillators,”
SIAM Journal on Applied Mathe-matics , vol. 65, no. 2, pp. 392–425, 2004.[12] E. V. Nikolaev and E. D. Sontag, “Quorum-sensing synchronization ofsynthetic toggle switches: A design based on monotone dynamical sys-tems theory,”
PLoS Computational Biology , vol. 12, no. 4, p. e1004881,2016.
UPPLEMENTARY M ATERIAL
Periodic averaging
We recall here that, from [10, Theorem 10.4], the periodic averaging method says that the solutions of the system ˙ x = εf ( t, x, ε ) (i)where f ( · ) is sufficiently smooth with respect to ( x, ε ) , and T -periodic and measurable in t , can be approximated by an autonomous average system ˙ x = εf av ( x ) (ii)where f av ( x ) = T (cid:82) T f ( s, x, ds. More precisely, let x ( t, ε ) and x av ( εt ) denote the solutions of (i) and (ii), respectively. If system (ii) has an exponentiallystable equilibrium point ¯ x av , then there exist positive constants ε ∗ and k such that, for all < ε < ε ∗ , system (i) has a unique, exponentially stable, T -periodicsolution ¯ x ( t, ε ) in a O ( ε ) -neighborhood of ¯ x av , that is (cid:107) ¯ x ( t, ε ) − ¯ x av (cid:107) ≤ kε . Moreover, if the initial conditions are such that x (0 , ε ) − x av (0) = O ( ε ) ,then x ( t, ε ) − x av ( εt ) = O ( ε ) , for all t ≥ . Nondimensionalization of system (1)-(4)
Equations (1) and (2) can be rewritten as d mRNA
LacI dt = κ m0L + κ mL (cid:16) TetRθ
TetR (cid:17) η TetR · w ( t ) − g mL · mRNA LacI d mRNA
TetR dt = κ m0T + κ mT (cid:16) LacIθ
LacI (cid:17) η LacI · w ( t ) − g mT · mRNA TetR where w ( t ) and w ( t ) are defined in (11)-(12).Now, taking into account that mRNA molecules are degraded faster than other molecules, we can obtain a quasi-steady stateapproximation by setting d mRNA LacI dt = 0 and d mRNA TetR dt = 0 , yielding mRNA LacI = κ m0L g mL + κ mL g mL
11 + (cid:16)
TetRθ
TetR (cid:17) η TetR · w ( t ) mRNA TetR = κ m0T g mT + κ mT g mT
11 + (cid:16)
LacIθ
LacI (cid:17) η LacI · w ( t ) Assuming that LacI and TetR proteins degrade at the same rate, that is g pL = g pT = g p , that η LacI = η TetR = 2 , and using thedimensionless state variables and time t (cid:48) = g p t , x = LacIθ
LacI , x = T etRθ
TetR , we obtain from (3) dx dt (cid:48) = κ pL θ LacI g p mRNA LacI − x = κ pL θ LacI g p (cid:34) κ m0L g mL + κ mL g mL
11 + x · w ( t (cid:48) /g p ) (cid:35) − x = k + k x · w ( t (cid:48) /g p ) − x with k and k as in (9), and from (4) dx dt (cid:48) = κ pT θ TetR g p mRNA TetR − x = κ pT θ TetR g p (cid:34) κ m0T g mT + κ mT g mT
11 + x · w ( t (cid:48) /g p ) (cid:35) − x = k + k x · w ( t (cid:48) /g p ) − x with k , k as in (10).
00 1000 1500 2000 2500 3000 3500 LacI20040060080010001200TetR (a) Equilibrium points for ¯ u aTc = 100 ng / ml and differentvalues of ¯ u IPTG .
500 1000 1500 2000 2500 3000 3500 LacI20040060080010001200TetR (b) Equilibrium points for ¯ u IPTG = 1 mM and different valuesof ¯ u aTc .Fig. 9. Equilibrium points ¯ x av of (20) as a function of duty cycle D rescaled in arbitrary fluorescence units using (7). Each dot represents the location ofthe unique stable equilibrium point of system (20) evaluated for D taking values in the interval [0 , with increments of . . x x (a) D = 0 . , T ≈ ( ε = 0 . ). x x (b) D = 0 . , T ≈ ( ε = 0 . ).Fig. 10. Background: phase portrait of the average system (20). Red line: the solution of the time-varying system (19) with ¯ u aTc = 50 ng / ml and ¯ u IPTG = 0 . from initial condition [1 , T .
500 1000 1500 2000 2500 3000 3500 LacI020040060080010001200TetR
Fig. 11. Orange dots: Mean-value, evaluated at regime, of the response of system (1)-(4) (with membrane dynamics (5)-(6)) to PWM inputs with T = 240 min and varying D from . to . with increments of . . Blue dots: corresponding equilibrium point ¯ x av ( D ) of system (20). Amplitude of pulse wavesset to ¯ u aTc = 100 ng / ml and ¯ u IPTG = 1 mM .(a) With instantaneous diffusion (b) With diffusion dynamics (5)-(6)Fig. 12. Time evolution of the time-varying system (19) (in solid lines) and of the average system (20) (dashed lines) from initial conditions [1 , T with ¯ u aTc = 50 ng / ml , ¯ u IPTG = 0 . , T = 180 min , D = 0 . .a) With instantaneous diffusion (b) With diffusion dynamics (5)-(6)Fig. 13. Time evolution of the time-varying system (19) (in solid lines) and of the average system (20) (dashed lines) from initial conditions [1 , T with ¯ u aTc = 50 ng / ml , ¯ u IPTG = 0 . , T = 180 min , D = 0 . .(a) With T = 180 min , that is ε ≈ . (b) With T = 45 min , that is ε ≈ . Fig. 14. Time evolution of the time-varying system (19) (in solid lines) and of the average system (20) (dashed lines) from initial conditions [1 , T with ¯ u aTc = 50 ng / ml , ¯ u IPTG = 0 . , D = 0 . , with instantaneous diffusion (Assumption 1).ABLE IP ARAMETERS OF THE TOGGLE SWITCH MODEL ( AS REPORTED IN [9])Transcription rates(mRNA min − ) κ m0L . · − κ m0T . · − κ mL . κ mT . Translation rates(a.u. mRNA − min − ) κ pL . · − κ pT . Degradation rates(min − ) g mL . · − g mT . · − g pL . · − g pT . · − plac regulationby LacI θ LacI . a.u. η LacI θ IPTG . · − mM η IPTG ptet regulationby TetR θ TetR . a.u. η TetR θ aTc . ng/ml η aTc IPTG exchange rate(min − ) k inIPTG . · − k outIPTG . · − aTc exchange rate(min − ) k inaTc . · − k outaTc . · −−