Chaos and its Degradation-Promoting-Based Control in an Antithetic Integral Feedback Circuit

— This letter deals with a novel variant of antithetic integral feedback controller (AIFC) motifs which can feature robust perfect adaptation, a pervasive (desired) ability in natural (synthetic) biomolecular circuits, when coupled with a wide class of process networks to be regulated. Using the separation of time-scales in the proposed kind of AIFC, here we ﬁnd a reduced-order controller that captures the governing slow part of the original solutions under suitable assumptions. Inspired by R¨ossler systems, we then make use of such a simpler controller to show that the antithetic circuit can exhibit chaotic behaviors with strange attractors, where the bifurcation from a homeostatic state to chaotic orbits can happen, e.g., when considering saturated Hill-type reactions for the actuation. Addition of degradation terms to the controller species, whether naturally due to dilution or exogenously using protein tags, is showcased by simulation results to be an effective solution in suppressing deterministic chaos and aperiodic oscillations. In the same vein, we recapitulate the recently introduced antithetic rein IFC motif and conﬁrm that the promotion of degradation by a rein mechanism also can control chaos and improve the stability of closed-loop circuit.


I. INTRODUCTION
R OBUST perfect adaptation (RPA), as a pervasive feature observed in evolved natural circuits, from cellular to organismal level, provides the living network with a structural architecture that maintains its functionality (e.g., the regulation of an output variable at given targets) regardless of unexpected perturbations and environmental uncertainties [1]- [4].To be more specific, such a (perfect) homeostatic property refers to the ability of a biological system in robustly returning the state of its output species back into the (same) prestimulus baseline, even when the external stimulation is still present [2], [5], [6].
This is proven to be achievable by structuring an integral feedback control (IFC) mechanism in the circuit of interest [5], [7], [8].In particular, the antithetic IFC motif (AIFC), initially developed in [3], has gained a significant attention in recent years [6], [9]- [12] due to its universal [4] capability of structuring IFC (and, hence, ensuring RPA) in control of a generally unknown biochemical process network.As an improved version of this motif, the authors of [13] have recently added an extra inhibitory reaction from a controller species to the target species being regulated in the process side, the addition of which it is believed to enhance the stability of resulting closed-loop circuit.They called it antithetic rein IFC (ARIFC) motif and proved its functionality for the case of controlling a simplified (two-species) gene-expression model.
Observation of chaos, extreme sensitivity to slight changes, has long been a classical phenomena of interest explored both experimentally and theoretically in a wide range of fields, from electrical [14] to biochemical systems [15].Thus far, however, there has been little reports on chaos tailored to synthetic biology contexts [16].In this regard, studying the functionality of genetically embeddable IFC structures in the presence of rich dynamics and chaotic attractors is of central interest to this letter.More specifically, we provide a case study based on a modified kind of minimal AIFC motifs to address the question of whether sustained nonperiodic oscillations can coexist with such RPA-achieving regulatory circuits.
To do so, in Section II, we initially sketch the involved chemical reaction network (CRN) and the associated dynamic model of the considered variant of AIFC.Motivated by the previous work [6], we next utilize geometrical singular perturbation tools [17]- [19] to derive a (zero-order approximated) quasi-steady state approximation (QSSA) for this motif under suitable assumptions, treating the original (full) circuit as a singular singularly perturbed system (SSPS).As the next step, based on the Rössler systems [20], we propose the CRN of a two-species autocatalytic process network whose closed-loop circuit, where the previously derived reduced model acts as its controller part, can generate chaos with spiral-type attractors.As such a dissipative oscillatory response indeed represents the slow solution of its original AIFC counterpart circuit, we may expect that the full antithetic system also behaves similarly.Chaos analysis in both of the reduced and full (autonomous) systems comes then into consideration to study this matter with more details, presented in Section IV.
From the fact that the concentrations of biomolecules found in naturally-evolved living circuits are being "diluted" out as the host cell's volume increases [10], it can be immediate that in the fast-growing cells the dilution effect may cause memory leakage in the stored integral variable, thereby giving rise to an imperfect adapting circuit [4], [6].Viewed in this way, Section V discusses the effect of such a leakage on the observed chaotic oscillations and interestingly show that it acts as a natural chaos control mechanism.This holds also true if the degradation terms has been added to the controller species exogenously, e.g., by the addition of protein degradation tags.More, this section suggests another biomolecular mechanism to suppress chaos, where we modify our chaotic circuits in such a manner that they encompass a rein structure similar to [13].Here we indeed try to recapitulate an enzymatic version of the ARIFC, increasing the rein reaction rate of which can help precluding chaotic orbits.Section VI concludes the letter.

II. CONSIDERED BIOCHEMICAL SYSTEMS
A. Full (AIFC-based) closed-loop circuit Consider that the process being controlled, say F (t), is consisted of n arbitrary, distinct species, say {X 1 , ..., X n }, connected to each other through an unknown, complex reaction network R x with uncertain parameters (rate constants) and reaction rates.Now, assume that the dynamic model of F (t) can be generally described by some nonautonomous (sufficiently smooth) vector field , where the state variable vector ≥0 represents the corresponding concentration of the species evolving over time with t ∈ [t 0 , t 1 ] Note the parametric space of the system that is not constrained in this setting to be necessarily time-invariant.Here we connect a biochemically constructible exogenous controller block (the AIFC-based motif under study) to this process and prove that robust (perfect) homeostasis can be ensured at steady-state level of the target species X n (i.e., the output concentration x n (t) as t → ∞) in the resulting closed-loop structure, regardless of what reactions constitute F (t) and how exactly they are mapped.To exhibit RPA, it implicitly requires the (positive) equilibrium(s) of closed-loop to be stable.This controller block, in a minimal manner, only incorporates two distinct (controller) species, called as Z 1 and Z 2 , for the desired integral feedback strategy to be structured in the closed-loop system.This happens by taking place a specific annihilation mechanism, consider the comparison (outflow) reaction R 3 = Z 1 + Z 2 η − → ∅, in which the two species sequester [21] each other identically with a rate η ∈ R >0 .It produces some inactivated species (the biomolecules Z 1 :Z 2 ) not affecting other reaction rates inside the network of interest (indicated symbolically by "∅", where ∅ represents the neutral annihilation product leaving the reacting system).
Choose the time-varying (differentiable) reference input signal µ(t) : R → R ≥0 , or the control input in a traditional sense, to be externally encoded into the closed-loop circuit as the propensity of firing an inflow reaction in the form R 2 = ∅ µ − → Z 1 , corresponding to which the species Z 1 is constitutively being produced.Then, let the process network to be sensed in the controller block by the sensing reaction R 5 = X n θ2 −→ X n + Z 2 , in which X n (catalytically) charges Z 2 with the sensor gain θ 2 ∈ R >0 .Moreover, let the controller block to convey the (computed) control signal, which encodes the time integration of output deviation from the set-point (i.e., (x n (t) − µ(t)/θ 2 )dt) as molecular dynamics, through taking place the (enzymatic) actuation reaction Where the output species is deemed as the substrate being converted to an enzyme-substrate complex (the species C), and where the species Z 2 from the controller side acts as the (total) enzyme quantity catalyzing the degradation of X n , with a catalytic rate θ 1 ∈ R >0 .The association (dissociation) rate of such an actuation reaction is assumed to be the positive scalar k f (k r ).
These reactions, taken together, shape such a closed-loop circuit shown in Fig. 1a, in which the two species Z 2 and X n establish a local negative feedback loop nearby.They can be represented as the set In this scenario, by assuming (deterministic) mass-action rate equations and also that the substrates in the actuation reaction are in considerable excess w.r.t the enzymes (making it legitimate to apply the QSSA of R 1 [22]), the dynamical behavior of the considered closed-loop control system can then be modeled by the following evolution equations where e i ∈ Z n is the standard basis of R n and K n , given by representing its enzyme affinity.Note the effect of diluting cellular context in the controller part that is approximated in (1b) and (1c) by some constant degradation terms with rates γ c1 , γ c2 ∈ R ≥0 .It implicitly assumes that the host cell grows exponentially in size.In the case that we neglect the effect of such a controller degradation, a hidden IFC architecture lied on the molecular concentrations of Z 1 and Z 2 can be further unraveled by formulating the dynamics of the integral variable ζ = z 2 − z 1 , i.e. ζ = θ 2 x n − µ(t).This mechanism apparently forces x n (t) to precisely track the set-point µ/θ 2 at steady state and helps it to (perfectly) reject a wide class of disturbances that may  (c) Target closed-loop output, x 2 (t), of both the full and reduced systems, obtained from ( 5) and ( 6) respectively for the default set P d and initial conditions I d , in two different scenarios of controller leakage.
perturb the system (1) meanwhile (where, equivalently, RPA ensues), the internal model principle says [7].It additionally necessitates the changing input µ to vary in a reasonably slow manner.Notice also that the positiveness of the species' concentrations in this circuit naturally imposes the positive orthant of R n+2 to be a positively invariant set for this system, so long as the prescribed initial conditions are also picked from the same set (i.e., [x(t 0 ), The CRN of considered (two-node) process network can be shown by the reaction set R 1b shows the schematic of this (autocatalytic) circuit.Also, the forward and reverse rate pairs for R x 1 and R x 5 are {α f 1 , α r 1 } and {α f 2 , α r 2 }, respectively.
III. MODEL REDUCTION ANALYSIS In this section, we analytically show that the solutions of (full) system (1) for large selections of η can be approximated over a finite time interval t = [t 0 , t 1 ] by that of the following reduced-order system (whose state variables are distinguished by bar marks) provided that some assumptions are met.
Assumption 1. Suppose that F and µ(t) are sufficiently smooth functions with respect to all their variables so that a unique, well-behaved solution of the initial value problem (2) exists for some . Assume this holds also for the full system (1) with In this letter, we shall restrict our attention to a particular class of processes that can together with the selected controller parameters hold the following assumption valid for the closedloop response (2).Assumption 2. For the obtained solutions from (2) with z(t 0 ) = z0 and xi (t 0 ) = x0 i , the state variable z varies for all t ∈ [t 0 , t 1 ] in a compact set S z ⊂ R ≥0 whose minimum is in the interior of positive orthant.
To prove that such a QSSA exists for the full system, we reformulate (1) as a standard singularly perturbed system (SPS) [17] by taking the perturbing (parasitic) parameter, called ∈ R >0 , as := 1/η.In this formalism, z i are regarded as the fast (exhausting) variables while x i (t) are considered to be the slow ones with dynamics independent of the rate constant η.Under Assumption 1 and for sufficiently small , one may find a slow (integral) manifold for the dynamics of SPS (3) with the form = M (X , t, ), Z = [z 1 , z 2 ] T , if it can be shown that the boundary layer (BL) solution of system (3) admits at least one isolated root which is asymptotically stable uniformly in t and x for the prescribed z 0 1 and z 0 2 .Such a BL system can be obtained from freezing the slow modes ( = 0) and dilating time (τ := t/ ), while treating the slow variables and time as parameters.However, a simple computation reveals that the Jacobian of BL dynamics, given by dZ/dτ = [−z 1 z 2 , −z 1 z 2 ], is singular throughout the state space.This property is a direct result of special annihilation mechanism employed to establish IFC, where each molecule of Z 1 binds only to one molecule of its sequestration pair Z 2 .
Applying a proper state transformation may be of help to regularize the SSPS problem (3).Particularly, we introduce the change of variable where z(t) is now deemed as a slow variable (because has disappeared from ż) with initial condition z(t 0 ) = z 0 2 − z 0 1 .Theorem 1. Hold Assumptions 1 and 2 true.Then, there exists a sufficiently large η * ≥ η 0 such that for η ∈ [η * , ∞) the following relations between time responses of the original system (1) and its QSSA obtained from (2) remain valid, provided that Proof: The inner layer dynamics of the transformed system (4) can be readily determined as dz 1 /dτ = −z 1 (z + z 1 ), which admits two different isolated quasi-steady states for z 1 (t), that are, z 1 = M 0 = 0 and z 1 = −z.Its Jacobian also can be computed as −∂(z 2  1 + zz 1 )/∂z 1 = −2z 1 − z.Assuming z > c as a fixed parameter where c is positive (Assumption 2), it is a straightforward exercise to solve this one-dimensional differential equation and see that z 1 = 0 is an exponentially stable equilibrium for this BL subsystem, whose domain of attraction includes all non-negative z 1 (t 0 ).Then, by referring to the Tikhonov's theorem (See [17, Th. 3.1], [18], [19]), one can easily follow the desired conclusion.

IV. CHAOS ANALYSIS
In this section, we consider the proposed circuit x as the process network being controlled in (2) and describe (with the aid of QSSAs [22] for R x 1 and R x 5 ) the dynamic model of the resulting reduced closed-loop circuit as where  1 = 75 nMh −1 , K 1 = 0.03 nM, γ c1 = 0 h −1 , γ c2 = 0 h −1 } and (x 1 (t 0 ), x2 (t 0 ), z(t 0 )) = (0.035, 5.2, 2.59), respectively, and make use of these two as the default setting for providing our numerical simulation results hereafter.For this setting, the (unique) positive equilibrium point of ( 5), denoted by p * := [x * 1 , x * 2 , z * ], can be written out as p * ≈ [0.05, 6.28, 5.3], implying the stability of closed-loop as the eigenvalues of its Jacobian, 6.This point can be confirmed from Fig. 1c, which evaluates the time response of this circuit when being perfectly adapted to a step-like perturbation increasing θ 1 by 100% in magnitude.
It is seen that a (spiral-type) attractor, illustrated in Fig. 2a, absorbs the solution of this system.Indeed, decreasing K n in P d has amplified the oscillations created by ( 5b 1 in P d with other values), the two oscillating variables will jump back and forth along the two branches of a Z-shaped hysteresis-type (slow) manifold repeatedly, in which case a chaotic response yet invariant w.r.t R 3 ≥0 can be observed.Taking one step forward, we show the bifurcation diagram of ( 5) for variations on K n as the bifurcation parameter.As can be seen from Fig. 2b, when 3.2 nM ≤ K n ≤ 6.5 nM, the system works in a steady state region with one fixed point (so that RPA also takes place), while the emerging bifurcation point at K n = 3.2 nM generates a period-1 limit cycle for 0.788 nM ≤ K n ≤ 3.2 nM, whose amplitude is decreasing w.r.tK n .At K n equal to 0.788 nM, 0.384 nM, and 0.31 nM, the period-doubling cascade converges to period-2, -4, and -8 orbits, respectively, and routs to chaos when K n ≈ 0.29 nM.Further details regarding chaotic regions and rich dynamics for smaller values of K n can be determined from the magnified area of Fig. 2b (indicated by a magnifier icon), whose results are also supported in parallel by providing the pertinent finitetime largest Lyapunov exponent (LLE) of ( 5), computed using the Wolf's method [23].Here the observation time, integration step size, and initial conditions used to calculate LLEs are 0 ≤ t ≤ 5000 h, ∆t = 0.001 h, and I d , respectively.We shall also note that the obtained state variable z(t) in such a chaotic setting has a minimum greater than c = 0.4 for all 0 ≤ t ≤ 5000 h.Thus, Assumption 2 is already met if we take t 0 = 0 and t 1 = 5000 h.
As the last step, by comparing ( 5) to (2), the analogous full (AIFC-based) circuit of (5), in the form (1), can be noticed as It is easy to check that both the reduced model (5) and full system (6) satisfy Assumption 1. Now, one can choose a nonnegative initial condition for z 1 and set x 0 i = x0 i , z 0 2 = z0 +z 0 1 .Then, Theorem 1 in Section III declares that, when 0 ≤ t ≤ 5000 h, the fast dynamics of ( 6) for (sufficiently) large η will be attracted (rapidly enough) to a slow manifold (z 1 = 0) on which the evolution of (slow) dynamics can be tightly captured by (5).Here Z 1 acts as a short-lived intermediate species (governing the transient) being waned faster as η becomes stronger.By this explanation it becomes clear that if the annihilation mechanism (R 3 ) occurs faster than the switching rates k c 1 and κ, we then may expect to have (for P d , K n = 0.15, η 1, and z 1 (t 0 ) ≥ 0) the slow solution of (6) generating chaos after a transient.
To confirm this point, we set x(t 0 ) = x(t 0 ) (where x(t 0 ) was already provided in I d ), z 1 (t 0 ) = 0.01, and z 2 (t 0 ) = 2.6 to illustrate the simulation results in Fig. 3.Note the trend in Fig. 3b when η ≈ k c 1 .Here a (Hopf) bifurcation at η ≈ 14 nM −1 h −1 routs to chaos for the first time.See also Fig.  for further explorations.Although there may not be available valid, accurate in vivo values for the affinity η at hand, at least in vitro studies initially suggest that, using some engineered protein-protein interactions displaying sequestration, e.g., by homodimerizations of (C62)Gcn4 2 or heterodimerizations of FosC and JunC (with resulting binding rates vary from 26.64 to 57.6 nM −1 h −1 ), the chaotic range for η evoked by Fig. 3 can be of biological relavance [21,Supplementary Table SII].

V. CHAOS SUPPRESSION MECHANISMS
This section introduces two different mechanisms which can control the observed chaos.The first approach is to add degradation terms to the controller species and the second one considers a rein feedback mechanism as effective solution.
In Section IV, given that γ c1 and γ c2 are set zero in P d , we have indeed supposed no degradation terms for the circuits of controller blocks.Controller degradation can be a direct result of cellular growth or might be determined more by, e.g., synthetically addition of some protein degradation tags.In both cases, it leads to a leaky integration [10] producing steady-state tracking errors (Follow this from Fig. 1c).Although having controller degradation may deteriorate the performance of circuit's adaptation mechanism, it is shown to be effective in enhancing the closed-loop stability [6], [9], a crucial prerequisite for RPA to be achieved.To have an intuition, form the gradient flow of the simpler system (5) and readily see that it holds ∇f ≤ ϑ − γ c2 , introducing γ c2 as a dissipative element individually increasing which steers the flow of biological system to a null-volume attractor.
Let us now simply add an extra actuation reaction, say R θ3 : Z 1 θ3 − → Z 1 + X 1 , from the controller species Z 1 to X 1 , taking into account that a new (positive) feedback loop comprising the reactions R θ3 , R x 5 , R 5 , and R 3 will now be a part of the closed-loop system.Assuming massaction rate kinetics for R θ3 , the addition of such a unimolecular activation will only affect the dynamics in (1a) as ẋ = F(x, t) − e n θ 1 z 2 x n /(K n + x n ) + e 1 θ 3 z 1 .It is now fairy easy to check that the resulting motif (depicted in Fig. 6a) follows a structure similar to the ARIFC motif in [13], with the only distinction that the rein mechanism here (R 1 ) is modeled by Hill reaction rates.One can follow this from (1a).Can this modification, even in a chaotic setting, still bring the stability improvements for the closed-loop circuit in the regime of high rein gains, i.e. θ 1 → ∞? Fig. 6b tries to address this question where the process being controlled is kept intact (R x in Section II-B).As can be seen (for η = 10 2 and θ 3 = 1), increasing the rein constant rate can lead the ARIFC-based circuitry to achieve asymptotic stability and avoid strange behaviors.Fig. 6c further highlights this fact for a wide range of variations on θ 3 .This point puts forth the idea that the promotion of degradation even by utilizing a rein structure can also play the role of a natural chaos control scheme.

VI. CONCLUSION
In this letter, we considered a modified version of the universal (standard) AIFC motif.Using singular perturbation methodologies and under suitable assumptions, we obtained a QSSA for the solutions of this motif when the sequestration binding η is sufficiently strong, considered as the associated reduced model.To be more specific, we indeed proved through a nonlinear approach that, for a particular class of solutions, the abundance of species X i and Z 2 from (1) can be made (arbitrarily) close to the values of xi and z obtained from (2), respectively, where the sequestration binding η is (sufficiently) strong and the initial conditions satisfy x 0 i = x0 i , z 0 1 ≥ 0, and z 0 2 = z0 + z 0 1 .This particular class of solutions shall be characterized by imposing Assumptions 1 and 2 on the dynamic model (2).
Benefiting from such a finding, we next suggested a chaotic three-species (autocatalytic) circuit based on the previously obtained reduced-order controller to show that its corresponding AIFC-based circuit also can generate chaos for η 1. Being limited to the considered biomolecular control systems, we lastly studied the underlying role of degradation promotion in damping chaotic oscillations and precluding aperiodic orbits.
One may incidentally look at the regulatory effects of employed rein IFC on a time-average level basis and notice that it still exhibits the ability to maintain homeostasis in face of undesired disturbances, even in a chaotic fashion.This is illustrated by Fig. 6d, where Avg(S(t)) = 1 t−t0 t t0 S(τ )dτ .It is seen that taking average over the output x 2 can filter even nonperiodic fluctuations about the set-point, whereby showing stringent resilience w.r.t.perturbations.Similar topics were previously taken into account by [8], [16] to address such phenomena, however, putting further effort to analytically study the performance of biological IFC schemes in a more general, chaotic setup might be of interest for future work.

Fig. 1 .
Fig. 1.(a) Considered AIFC structure (left) in control of an unknown process network (right).(b) The two-species network specified in Section II-B.(c) Target closed-loop output, x 2 (t), of both the full and reduced systems, obtained from (5) and (6) respectively for the default set P d and initial conditions I d , in two different scenarios of controller leakage.

Fig. 2 .
Fig. 2. Reduced system (5) can show chaos for small Kn.(a) Chaotic attractor of this system in phase space, where Kn = 0.15.A chaotic temporal evolution of the states is also shown.(b) Bifurcation diagram demonstrating the local extrema of x2 (t) for variations of Kn from 0 to 6.5.Also provided is the LLE for Kn ∈ [0, 0.625].P d and I d are the typical settings used.
) and (5c) in the variables x2 (t) and z(t), thereby giving rise to a 2D attractive limit cycle near the surface x2 -z, where a (dissipative) motion spirals out around the positive (unstable) focus p * .At this time, the equation (5a) plays the role of a threshold switch triggering which guides such a spiral flow and lifts it up along x1 .The nonlinear term −k c 1 x1 /(K 1 + x1 ) wherein reinjects such a flow, settling it back again near the spiraling-out motion.If we let the switch to work considerably faster than remaining reactions (Compare κ and k c

Fig. 3 .
Fig. 3. Full (AIFC-based) circuit also exhibits chaotic oscillations for large values of η.(a) Attractors of the system (6) for four different scales of η, which start by a limit cycle for η = 10 0 and become sparsely-and denselyfilled chaotic for stronger annihilations.(b) Bifurcation diagram showing the spectrum of oscillation amplitudes in x 2 (t), where η varies logarithmically from 10 −2 to 10 5 .The complementary LLE is drawn for η ∈ [0.5, 2000].

Fig. 4 .
Fig. 4. Comparing the bifurcation diagram (andLLE) of the reduced system (gray colors), which was previously provided in Fig.2b, with that of its corresponding full AIFC system obtained from (6) for = 10 3 (pink colors).