Stability analysis of an artificial biomolecular oscillator with non-cooperative regulatory interactions

ABSTRACT Oscillators are essential to fuel autonomous behaviours in molecular systems. Artificial oscillators built with programmable biological molecules such as DNA and RNA are generally easy to build and tune, and can serve as timers for biological computation and regulation. We describe a new artificial nucleic acid biochemical reaction network, and we demonstrate its capacity to exhibit oscillatory solutions. This network can be built in vitro using nucleic acids and three bacteriophage enzymes, and has the potential to be implemented in cells. Numerical simulations suggest that oscillations occur in a realistic range of reaction rates and concentrations.


Introduction
All organisms require timing circuits to orchestrate processes related to their life cycle, such as cell growth, metabolism, and division [36]. By building molecular timers from the bottom up, we have an opportunity to understand the design requirements to programme periodic biochemical behaviours. In addition, synthetic oscillators are useful components to direct autonomous molecular operations in vivo and in vitro [11,13,30,33,35].
In vitro nucleic acid oscillators can be built with a small number of parts, and their behaviour is quantitatively predictable [15,16,19,22,35]. Nucleic acids have become molecular building blocks for a variety of logic and dynamic circuits, because their thermodynamic and kinetic interactions can be programmed by choosing their sequence content with rational optimization algorithms. Existing nucleic acid oscillators however cannot be ported to the cellular environment, because they rely on the presence of multiple singlestranded or partially single-stranded DNA species, which are incompatible with the cellular machinery [15,16,19,22]. Here we describe a new nucleic acid oscillator architecture that has the potential to overcome this limitation, as it does not require single-stranded DNA molecules. A particularly interesting aspect of our circuit is that all regulatory interactions are non-cooperative. Therefore, the corresponding model does not include Hill-type nonlinearities, present in the majority of models for molecular non-equilibrium circuits. Our oscillator comprises three polymerases, two of which mutually regulate each other (Figure 1(a )). The interactions among enzymes are defined by four synthetic genes and four RNA species (Figure 1(b )). The activity of two of the enzymes is modulated by RNA species that serve as inhibitors or activators. The third enzyme species controls the baseline production of two of the RNA species, and has a net effect of counteracting the mutual regulation of the other two enzymes. For instance, let us consider the pathway by which enzyme E 2 is inhibited by enzyme E 1 and activated by enzyme E 3 . E 1 produces RNA species R 1 by transcribing gene g 1 ; R 1 binds to and inhibits enzyme E 2 , converting it to inactive enzyme E * 2 (a reaction experimentally demonstrated, for instance, in [23,24]). RNA species R 4 (transcribed by E 3 ) counteracts this pathway and causes reactivation of E 2 (conversion of E * 2 to E 2 ), because it is designed to displace R 1 bound to E 2 , and to titrate free R 1 as well. Similar reactions generate inhibition and activation pathways for E 1 (due to E 3 and to E 2 , respectively). Overall, these interactions contribute to creating a negative feedback loop. This system can be experimentally implemented using T7, T3, and SP6 bacteriophage RNA polymerases [20,21,31], which can be purchased off-the-shelf from many vendors. RNA sequences (known as aptamers [12]) that bind to bacteriophage RNA polymerases and work as inhibitors have been experimentally characterized [23,24]. RNA activators can be designed as strands whose sequences are complementary to the sequences of the inhibitors via the mechanism of strand displacement and strand titration [18,37]. We describe this system by means of ordinary differential equations (ODEs) built using the law of mass action, starting from a list of chemical reactions reported in Section 2. We demonstrate that the system is a candidate oscillator due to the sign pattern of its Jacobian matrix [4,5]; in particular we show that the system admits transitions to instability that are exclusively oscillatory.
Our analysis relies on monotone systems theory (background is provided in Section 3) and the theory of invariant sets. In Section 4 we study the capacity of this dynamical system to structurally exhibit sustained oscillations whenever it becomes unstable, in view of its particular Jacobian structure; this approach can be applied to a variety of chemical reaction networks, as we have shown, for instance, in the context of other titrationbased regulatory networks [9]. Structural (namely, parameter-free) results can greatly help unravel the functioning of biological systems, which are affected by intrinsic uncertainties and variabilities in their parameters, but can nonetheless exhibit an extraordinary robustness and resilience [3]. We conclude with a numerical bifurcation analysis and study of period and amplitude as a function of variations in individual parameters, showing that for realistic reaction rates the system exhibits oscillatory behaviours (Section 5). We previously described a two-enzyme oscillator relying on RNA aptamers [2,10]; we claim that a three-enzyme system is more tunable, and simulation results indicate that in a certain region of parameter space its amplitude can be modulated independently of the period.

A three-enzyme oscillator regulated by first-and second-order reactions
In the following, capital letters represent chemical species and the corresponding lowercase letters represent species concentrations (e.g. species A has concentration a). Our three-node oscillator is described by the biochemical reactions below. Reactions are grouped in two sets corresponding to functional modules ( Figure 2), whose common external input is E 3 . For simplicity we assume a common degradation rate for all products R i , i = 1, . . . , 4.

Module1 :
Module2 : The differential equations describing Module 1 are: (1) Figure 2. Schematic of the interconnections between reaction Modules 1 and 2, with enzyme concentrations as inputs and outputs.
The differential equations describing Module 2 are: The total concentration of E 1 and E 2 is assumed to be constant, and equal to e tot 1 and e tot 2 , respectively; hence, mass conservation laws yield e * 1 = e tot 1 − e 1 and e * 2 = e tot 2 − e 2 . The two modules are interconnected and form a feedback loop: Module 1 (associated with variables r 1 , r 3 and e 2 ) receives input e 1 from Module 2; in turn, Module 2 (associated with variables r 2 , r 4 and e 1 ) receives input e 2 from Module 1. Both modules receive input e 3 , which we assume is constant ( Figure 2); we assume that the timescale at which e 3 binds to a gene and transcribes RNA is fast relative to the other timescales in the system, so that its dynamics can be neglected; this assumption is sensible for short transcripts (30-60 bases). In the next sections we demonstrate that transitions to instability in this system can occur exclusively due to a pair of complex conjugate eigenvalues crossing the imaginary axis, hence sustained oscillations necessarily arise whenever the system is driven to instability. From numerical simulations it is apparent that the system can actually be destabilized, for suitable parameter choices, and is therefore a good candidate oscillator.

Background
We summarize several background notions that are required to introduce our main results in Section 4. Additional information can be found in references [4,5]. Consider a system: where μ is a real-valued parameter and f (·, ·) is a sufficiently smooth function, continuous in μ, satisfying the following Assumptions for every admissible value of μ.
Assumption 1: All the solutions of (3) are globally uniformly asymptotically bounded in the compact set S ⊂ R n .
Assumption 2: ∂f i /∂x j is either always positive, always negative, or always null in the considered domain.
Due to the monotonicity of f i (·) with respect to each argument x j , the Jacobian matrix J of system (3) is sign definite. Definition 3.1: Given a system with a sign-definite Jacobian J, its structure is the sign pattern matrix = sign[J]. The structure of system (3) is assumed to be invariant with respect to μ. Assumption 1 ensures that an equilibrium exists; all the following definitions refer to this equilibrium, which is, in general, a function of μ: f (x μ , μ) = 0. We assume thatx μ depends continuously on μ. Note that a suitable change of coordinates always allows us to shift the equilibrium to the origin, without affecting our analysis.

Definition 3.2: System (3) undergoes a transition to instability (TI)
at μ = μ * iff its Jacobian matrix J(x μ ) is asymptotically stable in a left neighbourhood of μ * , and unstable in a right neighbourhood. 1 A TI is simple if at most a single real eigenvalue or a single pair of complex conjugate eigenvalues crosses the imaginary axis. (3) undergoes an oscillatory transition to instability (OTI) at μ = μ * iff its Jacobian matrix J(x μ * ) has a single pair of pure imaginary eigenvalues, while all the other eigenvalues have negative real part:
We now provide general definitions for candidate oscillatory and multistationary systems. We consider system (3), with its given structure (invariant with respect to μ), under Assumptions 1, 2 and 3.

Definition 3.4:
A system of the form (3), with structure , is structurally a candidate (1) oscillator in the weak sense iff it admits an OTI for some μ = μ * ; (2) oscillator in the strong sense iff every simple TI (if any) is an OTI; Necessary and sufficient conditions characterizing strong and weak oscillators/ multistationary systems are provided in [4] in terms of cycles in the structure graph. We associate matrix with a directed n-node graph, whose arcs are positive (+1), negative (−1), or zero depending on the sign of the corresponding matrix entries. Definition 3.5: Given a graph, a cycle is an oriented, closed sequence of distinct nodes connected by distinct directed arcs. A cycle is negative (positive) if the number of negative arcs is odd (even). The order of a cycle is the number of arcs involved in the cycle. We say a system is critical when all negative cycles (if any) are of order two.

Proposition 3.6: A non-critical system is a candidate oscillator in the weak sense if and only
if its structure has at least one negative cycle (necessarily of order greater than two).

Proposition 3.7:
A non-critical system is a candidate oscillator in the strong sense if and only if its structure has only negative cycles.
Proofs for Propositions 3.6 and 3.7 can be found in [4].

Remark 1:
The results above are verified as well if we drop Assumption 1 and we restrict our analysis to solutions that belong to a compact positively invariant set S, with a nonempty interior and with no equilibrium points on the boundary.
The graph-based results in [4] have been generalized in [5] to the case of systems that are the sign definite interconnection of subsystems that are either monotone or anti-monotone (as in the case of our system). We provide below the definitions of monotone and antimonotone system.

Definition 3.8: A systemẋ
where u(·) ∈ R is a scalar, time varying input, is input-to-state monotone if, denoting as x 1 (t) and x 2 (t) the solutions of the system corresponding to inputs u 1 (t) and u 2 (t), the fact that where inequalities are intended to hold componentwise. The system is input-to-state antimonotone if the input has the opposite effect on the state, that is, if A simple characterization of input-to-state monotonicity and anti-monotonicity [1,28] can be provided by exploiting the concept of Metzler matrix: a matrix is Metzler if its elements satisfy a ij ≥ 0, ∀(i, j) such that i = j. Theorem 3.9: System (4) is input-to-state monotone if its Jacobian matrix J = ∂f /∂x is a Metzler matrix and ∂f /∂u ≥ 0 componentwise. Conversely, system (4) is input-tostate anti-monotone if its Jacobian matrix J = ∂f /∂x is a Metzler matrix and ∂f /∂u ≤ 0 componentwise.
A more general concept, which we will use in the following, is given by monotonicity (or anti-monotonicity) with respect to a given signature tuple (s 1 , . . . , s n ), where s i = 1 or −1 for all i [14]: this amounts to requiring that, after changing the sign of the state variables aŝ x i = s i x i for all i, the system becomes monotone (or anti-monotone). Hence, Theorem 3.9 applies to the system in the new coordinates.

Proof:
The existence of the compact invariant set S ρ where the solutions of the system are globally uniformly asymptotically bounded (Proposition 4.1) implies the existence of a steady-state [25,26,29].
We later demonstrate that this steady state is unique.

Monotonicity properties and uniqueness of equilibrium point
Now we show that the overall system is the feedback interconnection of two subsystem, corresponding to the modules defined earlier, that are, respectively, anti-monotone and monotone. This property further implies that the system admits a unique equilibrium.
We individually linearize subsystems (1) and (2) around an equilibrium point (which is guaranteed to exist), and we begin by studying each subsystem in isolation.
where the linearized state variables of each subsystems are z = [δr 1 δr 3 δe 2 ] and w = [δr 2 δr 4 δe 1 ] . We denote equilibrium values of each variable with a¯symbol (e.g.ē 1 is the equilibrium of e 1 ). The linearized dynamics are defined by matrices: The two linearized subsystems are stable, and the matrices defining their dynamics (Jacobian matrices of the nonlinear systems) are Metzler up to changes in the sign of some variables. This can be easily shown by changing sign to the first component of z and to the second component of w: z 1 := −z 1 and w 2 := −w 2 . This is equivalent to changing sign to δr 1 and δr 4 , where r 1 and r 4 are variables of the original system, and provides matrices: and  Proof: Consider systems (5) and (6), which after the sign change have matrices (7) and (8).
Since all of their off-diagonal entries are non-negative,Â 1 andÂ 2 are Metzler matrices. They are also irreducible. 2 Hurwitz stability (all the eigenvalues of the Jacobian J = ∂f /∂x(x) have a negative real part) immediately follows from the fact thatÂ 1 andÂ 2 are Metzler and diagonally dominant, with negative diagonal entries (this is a consequence of Gershgorin's circle theorem). Finally, any stable and irreducible Metzler matrix has an element-wise negative inverse (see [6] for details).
We are now ready to demonstrate the monotonicity properties of the two nonlinear modules. (1) and (2) are, respectively, input-to-state anti-monotone and monotone after the sign change in the relative variables: δr 1 = −δr 1 and δr 4 = −δr 4 . ( 9 ) Proof: This follows from Theorem 3.9, since the state matricesÂ 1 andÂ 2 are Metzler, while the input matricesB 1 andB 2 are, respectively, non-positive and non-negative.

Proposition 4.4: Systems
Monotonicity and stability have important consequences on the static input-state and input-output characteristics (input-output equilibrium conditions) and on uniqueness of the equilibrium point. Indeed, the feedback of two systems that are either monotone or anti-monotone always admits a single equilibrium point (if any).
We have shown in Proposition 4.2 that an equilibrium always exists; we prove below, for completeness, that the static input-output characteristics of the two modules are monotonic, hence such an equilibrium point is unique.

Proof:
We recall that, for a generic systemẋ = f (x, u), the steady-state characteristicx(u) is implicitly defined by We can apply the implicit function theorem to find its derivative: Consider Module 1, after the sign change in the variables in Equation (9): The inequality holds componentwise (Proposition 4.3), hence after the sign change equilibriar 1 (e 1 ),r 3 (e 1 ) andē 2 (e 1 ) are monotonically decreasing functions of e 1 .

Proof:
The system always admits a steady state, as shown in Proposition 4.2. Due to Proposition 4.5,ē 2 (e 1 ) is a decreasing function andē 1 (e 2 ) is an increasing function. Thus, the system of equations: has a unique solution.
It is possible to demonstrate that this unique equilibrium is strictly positive, and there cannot be equilibria with zero components. This claim can be proved by showing that the   Table 1. Right: Trajectories in the plane e 1 -e 2 (black) and equilibrium conditions (red and blue).
two equilibrium equations intersect for positive values of e 1 and e 2 . Then, we can show that all other variables have a positive steady-state from their equilibrium conditions, which are all derived analytically in the appendix.

The interconnected system admits exclusively oscillatory transitions to instability
Based on the properties demonstrated in the previous sections, we establish that our threeenzyme network has the appropriate structure to exhibit sustained oscillations, whenever it is driven to instability. More precisely, the network admits exclusively oscillatory transitions to instability. (1) and (2) is a strong candidate oscillator.

Proposition 4.7: The interconnection of systems
Proof: The Jacobian of the overall system, with variables ordered as (δr 1 , δr 3 , δe 2 , δr 2 , δr 4 , δe 1 ) and with the variable sign change δr 1 = −δr 1 and δr 4 = −δr 4 , highlights that the  Table 1). Each black dot in this plot indicates that the (randomly) chosen parameter vector results in oscillations. Axes are in log scale. Orange diamonds represent the nominal value of each parameter (Table 1). Each parameter was varied between one tenth to 10 times its nominal value (black diamond; nominal values listed in Table 1). Orange regions indicate oscillatory behaviour; blue regions indicate a single stable equilibrium.
system is the negative feedback interconnection of two monotone subsystems: Due to Proposition 4.1, the system satisfies Assumption 1. By inspecting the Jacobian matrix, it is apparent that Assumptions 2 and 3 are also satisfied. Therefore, the system is a strong candidate oscillator [4,5]. This means that the system can transition to instability exclusively due to a pair of complex conjugate eigenvalues crossing the imaginary axis (OTI) and yielding oscillatory dynamics.

Numerical analysis
Model (1)-(2) was integrated using the MATLAB routine ode23. Bifurcation analysis, period and amplitude computation was also done writing MATLAB scripts ad hoc.
In the numerical analysis that follows, we choose nominal parameters ( Table 1) that are compatible with reaction rates measured in nucleic acid strand displacement reactions and in vitro transcription. An example solution trajectory for Model (1)-(2), integrated with the nominal parameters, is shown in Figure 3.

Randomized parameter sampling
First, we selected random values for the parameters sampling from a uniform distribution in the interval 10 −2 to 10 2 times the nominal parameter value (Table 1). We locate peaks and wells of the oscillations and compute period and amplitude as averaged over all the measured peaks and wells. A trajectory is classified as oscillatory if at least three oscillations are measured, if the period of the trajectory is between 0.5 and 40 h, and its amplitude is Figure 6. Period (h) as a function of each parameter (x axis in log scale). The period was computed numerically for damped and sustained oscillations. We classify a solution as oscillatory (damped or sustained) as long as the period is between 0.5 and 40 h, and the amplitude is larger than 1 nM. Blue circles indicate when the Jacobian has at least one pair of complex eigenvalues with negative real part (damped oscillations). Red circles indicate when the Jacobian has at least one pair of complex eigenvalues with positive real part (sustained oscillations). The parameters were changed in the range of one tenth to 10 times their nominal values. larger than 1 nM. This plot highlights that high degradation rates and low concentrations of e 1 and e 2 are associated with loss of oscillations ( Figure 4).

Bifurcation analysis
Using analytical equilibrium conditions (expressions (A.1) and (A.2) reported in the appendix), we find equilibria numerically and compute the eigenvalues of Jacobian (10) at the equilibria. If at least one pair of complex conjugate eigenvalues with non-negative real part is found, the equilibrium is classified as oscillatory. We vary two parameters simultaneously, while all others are kept constant as in Table 1. Oscillatory regions are shown in orange in Figure 5, while stable regions are shown in blue.

Period and amplitude
We focus on the influence of reaction rates and total concentrations of e i on the period and amplitude. Parameters α 1 , α 2 , α 3 , α 4 , e 3 , e tot 1 and e tot 2 are particularly relevant because they are experimentally easy to change (Figure1(b)): α i , i = 1, . . . , 4, are transcription rates, Figure 7. Amplitude (nM) as a function of each parameter (x axis in log scale). We computed numerically the amplitude of the solutions, as long as they classify as damped or sustained oscillations (period between 0.5 and 40 h and amplitude larger than 1 nM). Blue circles indicate when the Jacobian has at least one pair of complex eigenvalues with negative real part (damped oscillations). Red circles indicate when the Jacobian has at least one pair of complex eigenvalue with positive real part (sustained oscillations). The parameters were changed in the range of one tenth to 10 times their nominal values. which can be tuned by mutating the promoter region; e tot 1 , e tot 2 and e 3 can be chosen by the experimenter.
We compute the period and amplitude from integrated solutions to the ODEs. As explained in Section 5.1, we locate peaks and wells of the oscillations and compute period and amplitude as averaged over all the measured peaks and wells. A trajectory is classified as oscillatory if its period is between 0.5 and 40 h, and its amplitude is larger than 1 nM. The results are shown in Figures 6 and 7, where each individual parameter is varied in the range of one tenth to 10 times its nominal value, while other parameters are held fixed at their nominal value (Table 1). Correlation between period and amplitude is shown in Figure 8. To discriminate between damped and sustained oscillations, we check the sign of the real part of complex conjugate eigenvalues of the Jacobian matrix (evaluated at the considered combination of parameters). In Figures 6-8, damped oscillations are marked by blue circles, and sustained oscillations are marked by red circles.
From Figures 6 and 7 we observe that the period can be tuned from 0 to 5 hours. Also, the parameters related to the kinetics rate can change the period up to 3 hours in the range of one tenth to 10 times their nominal value. These plots show that when varying e 3 in a range between 0.1 and 10 times its nominal value, the period remains flat. In that same range, amplitude varies significantly. We also observe that varying δ 1 between 0.1 and 10 times its nominal value, amplitudes stays flat while the period varies between 0 and 3 h. It is worth noting that the titration rates δ 1 and δ 2 do not affect drastically neither amplitude nor period, which indicates that the system performance is robust relative to variations in the titration rates.
We observe that there is a range in which parameters α 2 and α 4 could be varied to tune exclusively the period, while the amplitude remains nearly constant. Alternatively, there is a range in which parameters e tot 1 and e 3 could be varied to modulate exclusively the amplitude, keeping the period nearly unchanged (and slow).

Conclusion
We have described an artificial three-enzyme biochemical network that has the capacity to oscillate. The network is designed for in vitro implementation with nucleic acid components and bacteriophage RNA polymerases, but has the potential to be implemented in vivo as well. The polymerases transcribe synthetic genes whose RNA transcripts in turn regulate enzyme activity, generating a negative feedback loop that is necessary for oscillations (the famous Thomas' conjecture [27,32]). We analytically demonstrate that this architecture can exclusively undergo oscillatory transitions to instability, due to the structure of its Jacobian matrix. Numerical analysis shows that in a range of realistic parameters the system oscillates; simulations are useful to direct the experimental implementation of this circuit, which is currently being pursued.

Notes
1. The definition holds as well for systems transitioning to instability from the right to the left neighbourhood of μ * : just takeμ = μ * − μ as the bifurcation parameter. 2. A matrix is irreducible if there does not exist a permutation of its rows or columns that transforms it into a block triangular matrix.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work was supported by the US National Science Foundation under grant [CMMI 1266402].