Degradation rate uniformity determines success of oscillations in repressive feedback regulatory networks

Ring oscillators are biochemical circuits consisting of a ring of interactions capable of sustained oscillations. The nonlinear interactions between genes hinder the analytical insight into their function, usually requiring computational exploration. Here, we show that, despite the apparent complexity, the stability of the unique steady state in an incoherent feedback ring depends only on the degradation rates and a single parameter summarizing the feedback of the circuit. Concretely, we show that the range of regulatory parameters that yield oscillatory behaviour is maximized when the degradation rates are equal. Strikingly, this result holds independently of the regulatory functions used or number of genes. We also derive properties of the oscillations as a function of the degradation rates and number of nodes forming the ring. Finally, we explore the role of mRNA dynamics by applying the generic results to the specific case with two naturally different degradation timescales.


Unattainability of chaotic regimes
One of the difficulties in extending the results to more than three dimensions is the possibility of chaotic attractors in the system. Nevertheless, monotone systems such as the ones described in equation (6), follow a Poincaré-Bendixson-type result [Mallet-Paret and Smith, 1990]. This implies that the ω-limit set of an orbit is either a steady state or a limit cycle. According to theorem 4.1 of [Mallet-Paret and Smith, 1990], the Poincaré-Bendixson-type result applies when R n + is positively invariant for monotone cyclic feedback systems defining a negative feedback loop that contain a unique critical point x * with det(−J) > 0. In our system, at the unique fixed point det(−J) = (1 + A) N n=1 δ n > 0 (see eq. (9)). Additionally it is easy to show that R n + is positively invariant since at each hyperplane x n = 0 the dynamics of the n-th component followṡ x n = δ n f n (x n−1 ) > 0 (see eq. (6)). As a result of this theorem, the system can only converge to a steady state or a limit cycle. The possible bifurcations in the N-dimensional space have a 2-dimensional analogue, excluding bifurcations such as period doubling. Therefore, the change of stability of any N-dimensional ring with N n=1 λ n = 0 (see eq. (9)) can only correspond to a Hopf bifurcation.

Derivation of the Implicit Equation
Introducing the Hopf bifurcation condition λ = iα in Equation 9 we obtain, Taking the argument of both sides and using the identity Similarly, taking the squared moduli of both sides of Equation S.1, The left hand side is increasing in α, so to find the smallest value of A for which there is a pair of imaginary eigenvalues, we must find the smallest α which satisfies Equation S.2. This will happen when the sum of the arguments equals π (k = 0), because each of the terms in the sum in Equation S.2 runs from 0 to π/2 as α increases. This reduces the Implicit Equation

Increase in oscillation probability after introducing a new species
The change of the parameter size for which the system oscillates as a function of the N + 1th degradation rate of a ring of N + 1 species can be obtained from Eq. 10 as > 0 if and only if N n=1 sin 2 θ n (cot θ N +1 − cot θ n ) > 0. Since increasing δ N +1 , while keeping the other degradation rates fixed, decreases θ N +1 and increases the other θ n , every term in this sum is increasing in δ N +1 . This means that A decreases with δ N +1 until it reaches a certain value (which must be within the range of values of the other degradation rates). Above this value A increases. Thus the probability of oscillations increases from zero when δ N +1 is zero to a maximum value for some optimal value of δ N +1 (intermediate between the δ n , n = 1, . . . , N ) and then decreases to the same probability as for the network without species N + 1 as δ N +1 tends to infinity.

Persistence of instability for A >Ã
The characteristic polynomial (Equation 8) is analytic in λ and A. Therefore by the implicit function theorem, the roots of the equation are analytic in A, except at the turning points of the polynomial. At these turning points it is straightforward to check that two negative real eigenvalues collide to form a pair of complex conjugate eigenvalues in such a way that the real and imaginary parts on each branch are continuous in A. The other eigenvalues are analytic in A at these points. We have already shown in the previous section that there are no zero eigenvalues for any value of A, therefore the real parts of an eigenvalue can only change sign as a pair of eigenvalues crossing the imaginary axis. From the multiplicity of the Implicit Equation (Equation S.2), we already saw that there will be other values of α k > α for which a pair of eigenvalues crosses the imaginary axis. In our analysis in the previous section we chose α to give the minimum value of A = A. The other points where eigenvalues cross the imaginary axis will occur at A k > A. This is clear from the LHS of Equations S.2 and 10, which are both increasing functions of α.
Nevertheless, each of the terms of the sum can contribute at most π/2 to the result. Therefore, the sum has the upper bound N n=1 tan −1 (α/δ n ) ≤ N π/2, and for finite α (necessary since A is finite), the inequality is strict. Therefore, when N = 3, 4, 5, 6 there can only be one possible solution α of the Implicit Equation; when N = 7, 8, 9, 10 there can be two values of α, etc. Concretely, the maximum number of times there are eigenvalues crossing the imaginary axis is (N + 1)/4 , where · is the floor operator.
Each time a crossing occurs, it must only consist of a single pair of eigenvalues. This can be seen by supposing the opposite case in which iα is a repeated root of the characteristic equation (Equation In this case it will also be a zero of the derivative of the characteristic function P δn−iα δn 2 +α 2 , which is impossible since each term in the final expression has negative real part. On the other hand, as A → ∞, the characteristic equation (Equation 8) requires the moduli of the roots to be arbitrarily large, so that the value of λ/δ n will eventually greatly outweigh 1 and λ will approach a value satisfying This has solutions λ k = N A N n=1 δ n ω k with k = 1, ..., N , where ω k = e (2k−1)πi/N are each of the N solutions to the N th root of −1. Therefore, the number of eigenvalues with a positive real part will be the number of ω k with arg(ω k ) ∈ (−π/2, π/2), which for N > 2 is 2 (N + 1)/4 . Thus, there are twice as many roots with positive real part as there are pairs of eigenvalues crossing the imaginary axis as A goes from 0 to ∞. This shows that every pair of eigenvalues crossing the imaginary axis must do so by changing the real part from negative to positive as A increases. Therefore, the steady state is unstable for all A >Ã.

Simulations of networks
The networks were simulated using Python custom code integrated using the PyDSTools dynamical systems environment [Clewley, 2012]. Trajectories were solved by precompiling in C a multistep Radau method. The source code has been made available as part of the Supplementary Material. Bifurcation diagrams were obtained using continuation techniques and integrations under the same Python environment.

First Lyapunov coefficient of the Hopf bifurcation
In the system with the same repressive function, f , for each species, the steady state F (x * ) = x * occurs when f (x * ) = x * . A Hopf bifurcation occurs atx, where f (x) =x and f (x) = − sec(π/N ).
In order to study the nature of the Hopf bifurcation, it is necessary to compute the First Lyapunov coefficient 1 evaluated at the critical pointx, that can be obtained as [Kuznetsov, 2013], which involves the computation of inner products of vectors in C N : a, b =ā T b. The matrix J in Equation S.7 is the Jacobian of the dynamical system at the critical point. For the system with same repressive functions and degradation rates for each species this is given by Note that for the sake of simplicity we have scaled the time with the degradation rate δ. Additionally, α = tan(π/N ) in Equation S.7 is the modulus of the pair of eigenvalues at the Hopf bifurcation and q the corresponding eigenvector where ω = − cos(π/N ) + i sin(π/N ) is an N th root of unity. On the other hand, p is the eigenvalue of J T , corresponding to eigenvalue −i tan(π/N ), that in this case is the same, p = q. The only missing ingredients to start computing 1 (x) are the definition of the bilinear and trilinear forms Now we can proceed by calculating sequentially the different terms of Equation S.7. The evaluation of the bilinear forms B(q, q) and B(q,q) is immediate, Thus, the product J −1 B(q,q) requires the knowledge of the row sums of J −1 . Since the row sums of J are the same and equal to (−1 − sec(π/N )), J is a stochastic matrix multiplied by −1 − sec(π/N ). Thus J −1 will also be a stochastic matrix multiplied by (−1 − sec(π/N )) −1 , and its product with the bilinear form will be, [cos(π/N ) + i sin(π/N )]. (S.8) Next, in order to compute the third term of Equation S.7 we need to work with the inverse (2 tan(π/N )i − J) −1 , which has rows which are permutations of each other as, This means that B(q, q) is an eigenvector of (2 tan(π/N )i − J) −1 with eigenvalue [b 1 + b 2 ω 2 + . . . + b N ω 2(N −1) ]. Thus it must also be an eigenvector of (2 tan(π/N )i − J) with eigenvalue 1 b1+b2ω 2 +...+b N ω 2(N −1) . On the other hand we know that, (2i tan(π/N ) + 1)ω 2N −2 + sec(π/N )ω 2(N −2) . . . , which means that the eigenvalue is 1 + 2i tan(π/N ) + sec(π/N )ω 2 , therefore B(q, And the third inner product of Equation S.7 can be written as, f (x) 2ω3 N 1 1 + 2i tan(π/N ) + sec(π/N )ω 2 , which has a real part, where c = cos(π/N ). Finally, we need to use the trilinear form to compute the first term of Equation S.7, Putting the three terms together, the first Lyapunov coefficient is given by (with c ≡ cos(π/N )): This quantity does not have a determined sign and therefore the Hopf bifurcation can be supercritical or subcritical depending on the details of f (x). For instance, for the repressive function f (x) = a/(1 + (1 + x − x 2 + x 3 ) h ) the sign of 1 will depend on the coefficient h (see Fig.S.1).

The Hopf bifurcation is supercritical if and only if
The coefficient multiplying f (x) 2 is negative for all N and decreases monotonically from −2/3 for N = 3 to −3/2 as N → ∞. This means that the criterion on f at the fixed point necessary for supercriticality becomes less stringent as N increases, so the Hopf bifurcation (for a given f ) is more likely to be supercritical as N increases. If, as in [Buşe et al., 2009], we consider f to be a repressive Hill function, i.e. f (x) = c/(1 + x r ), then −f (x) f (x) 2 = r 2 − 15r + 38 2(r − 5) 2 .
For r > 2, this is less than 2/3 and therefore the Hopf bifurcation is supercritical for any N ≥ 3. Note that we have shown that the Hopf bifurcation occurs when f (x) = − sec(π/N ), but for the Hill function system f (x) = −rx r /(1 +x r ), so |f (x)| < r, so for the Hopf bifurcation to occur requires r > sec(π/N ) ≥ 2. Therefore the Hopf bifurcation only occurs for r > sec(π/N ) ≥ 2 and is always supercritical. This extends the result of [Buşe et al., 2009] to the case of N species.