Oscillations in well-mixed, deterministic feedback systems: Beyond ring oscillators

Highlights • I present a way of breaking down regulatory networks to find Hopf bifurcations. This helps find optimal conditions for oscillations in dynamical systems models of these networks.• In a model of negative auto-regulation of a gene by its dimeric protein, it is optimal for the monomer to degrade faster than the mRNA and the mRNA to degrade faster than the dimer.• Adding a weak positive feedback loop to a repressilator increases the probability of oscillations.• The optimal degradation rate of species in the sub-loop is higher than that of species outside it.• The opposite is true for a negative feedback sub-loop or a very strong positive feedback sub-loop.


Introduction
Oscillations in gene expression are important in a variety of biological processes, from circaidian rhythms ( Reddy and Rey, 2014 ) to somitogenesis ( Dequéant et al., 2006;Hirata et al., 2002 ) to apoptosis in response to DNA damage ( Bar-Or et al., 20 0 0;Michael and Oren, 2003 ). Determining the conditions necessary for oscillations to occur is thus important in the understanding of these processes.
A ring oscillator is typically a gene regulatory network, containing a single negative feedback loop between the biochemical species. Gene/ biochemical species 1 regulates the expression of species 2 which regulates species 3 and so on. The final species (species n ) regulates species 1, see Fig. 1 a. A species can up-regulate or down-regulate the next. Up-regulations are indicated by sharp-headed arrows and down-regulations by flatheaded arrows, see Alon (2007) for a review of network mo-E-mail address: karen.page@ucl.ac.uk tifs and De Jong (2002) for a review of mathematical models of gene regulatory networks. If the number of down-regulations in the cycle is odd, so that the feedback loop for species 1 (or any other species) on itself is negative, oscillations are expected ( Smith, 1987;Müller et al., 2006;Bu ¸s e et al., 2009;. If the number of down-regulations is even, so that the feedback loop is positive, multi-stability is expected ( Smith, 1987;Plahte et al., 1995;Cinquin and Demongeot, 2002;Soulé, 2003;Angeli et al., 2004;Müller et al., 2006;Kaufman and Soulé, 2019 ), although long-lived oscillations are possible, for example, in systems with stochasticity ( Strelkowa and Barahona, 2010 ).
Previous work has shown that when the systems giving rise to oscillations are ring oscillators, that is they contain a single negative feedback loop between the chemical species, oscillations are most likely if the degradation rates of all species are equal . Frequently, oscillations arise in regulatory systems which are more complicated and contain more than one feedback loop, see for example ( Panovska-Griffiths et al., 2013;Perez Carrasco et al., 2018;Verd et al., 2017;2018;Goodfellow et al., 2014 ). Here I investigate the conditions necessary for oscillations in these multiple feedback loop systems. The study is limited to ordinary differential equation models.
The general theory for predicting the onset of oscillations is presented in Section 2 . Section 3 concerns the case of networks consisting of a loop and a sub-loop and explores optimal conditions for oscillations. Section 3.1 concerns the example of a single auto-repressive gene, and the species modelled are its mRNA concentration and the concentrations of its protein in the form of a monomer and of a dimer, see Fig. 1 b. Section 3.2 concerns threespecies systems with a loop and sub-loop, such as the ACDC motif, see Fig. 1 c. Section 3.3 reverts to the n -species case and the optimal conditions for oscillations. Section 4 concerns networks with two interconnected rings, exploring a four-species example, which cannot be described as a loop with a sub-loop. Finally the results are discussed in Section 5 .

The feedback system
In the most general case of a feedback system, there are n species and any species can regulate any other. This can be represented by a directed graph, with two types of edges. Sharp arrows indicate that the biochemical species at the base of the arrow positively influences the other, while flat-headed arrows indicate that the biochemical species at the base negatively influences the other. In gene regulatory networks, transcription of a gene is typically regulated by the protein of the other gene. Proteins that do this are referred to as transcription factors. Therefore, to a good approximation, it is the mRNA of the second gene that is directly regulated by the protein of the first. In some models, this is dealt with directly (see for example Elowitz and Leibler, 20 0 0;Page and Perez Carrasco, 2018 ). It is common, however, to assume that the dynamics of mRNA are fast compared to those of protein so that their concentrations are at quasi-steady state. Here we not only consider transcriptional regulation, but also consider reactions between biochemical species, for example dimerisation of a protein.
The time evolution of the concentrations of the biochemical species of interest can be described by differential equations. If it is assumed that the biochemicals are well-mixed within the cell and sufficiently abundant, then ordinary differential equations can be employed (see Discussion for a description of more refined models in which this assumption is not made). Here we assume, in addition, that all species, be they mRNAs or proteins, degrade at a rate which is linear in their concentration. Within the general theory, it would be fairly straightforward to replace this assumption with degradation at a rate which is monotonically increasing in concentration (see for example Alrikaby, 2018 ).
As in Page and Perez Carrasco (2018) , we model the concentrations of the species with a system of ordinary differential equations: where d i is the degradation rate of species i and d i f i describes its regulation by the other species. f i is some function of the concentration of the n -species. In the case of a ring oscillator, f i depends only on x (i −1) mod n (see Fig. 1 a). In addition, the number of downregulations is odd.
The steady states of the system 1 are given by solutions to the system of equations: A straightforward argument shows that, for the ring oscillator, there is a unique steady state (see for example Page and Perez Carrasco, 2018 ). Rings in which there are an even number of down-regulations typically show multi-stability and are more complicated to analyse (see e.g. Smith, 1987;Müller et al., 2006;Strelkowa and Barahona, 2010 ). We assume that since we model chemical reactions, transcription and translation only, the dynamics will be bounded. We look for the onset of oscillations by locating Hopf bifurcations. According to linear analysis, these occur when the eigenvalues of the Jacobian matrix at the steady state are imaginary.
We consider a Hopf bifurcation in such a system and derive the characteristic equation, setting the eigenvalue to be i α: where all the partial derivatives are taken at the relevant steady state and δ is the Kronecker delta. For a given permutation σ , let us divide the elements of { 1 , 2 , . . . , n } into those for which j = σ j , which we call fixed points, and those for which this is not true, which we call moving points. The contribution from the term corresponding to σ in the determinant is only nonzero if the moving points of σ correspond to a set of mutually exclusive feedback loops. It is straightforward to see this: σ can be decomposed into a product of cycles and the term in the determinant corresponding to σ is a product of the products over these cycles. For the moving points, the product for each cycle is simply the products of the derivative of f j with respect to x σ j within the cycle. These are only nonzero if species σ j regulates species j for each j in the cycle and therefore the cycle forms a feedback loop. Therefore, further assuming that there is not direct self-regulation of species ( for all i ), the characteristic equation is given by: 0 = S ∈ set of sets of mutually exclusive feedback loops where l S is the number of loops in S and H S is the product of the gradients of the functions representing each link in each loop in S at the steady state. Hence positive feedback loops give rise to positive parameters, H S , and negative feedback loops give rise to negative parameters (the overall sign of H S is the product of the signs of the individual loops). The parameter is multiplied by a factor (1 + iα/d j ) for each value of j representing a species not in the set of feedback loops. For example, the original ring contains all species and is therefore not multiplied by any of these factors. We must also include the empty set in the sum, with H ∅ = 1 . Eq.
(2) provides a systematic way of organising the characteristic equation of an ODE model of a regulatory network. It will facilitate the identification of Hopf bifurcations.

Analysis of system with a loop and a sub-loop
Let us start by considering the case where there is a ring feedback involving all the species which is a negative feedback and in which there is one feedback loop other than the ring and the empty set, we refer to H for the whole ring as F and H for this other loop as G . We refer to the set of species involved in the other feedback loop as W . We define θ j = tan −1 (α/d j ) , for j = 1 , ., n, with θ j ∈ [0, π /2). Then the characteristic equation becomes When G was zero, this had a solution for minimal | F | when n j=1 θ j = π . We can show that if G is greater than zero, the value of n j=1 θ j , and hence α and hence | F | at the bifurcation is increased, whereas if G is negative it is decreased.
Since in the case with just a ring, oscillations were favoured by equal degradation rates and since there is no means to distinguish within the characteristic equation between species involved in the same feedback loops, we assume that the degradation rates within W are the same and the degradations of the other species are the same as each other, but potentially different from those of the species in W . We call angles within W, θ , and those outside W, φ. Then Consider the sine rule applied to the triangle in the complex plane whose vertices are the origin and the points −G and −G + This means that Suppose, G is fixed and we want to minimise the value of | F | at the bifurcation. We first pause and study a couple of simple examples, with n = 3 and | W | = 2 .

Simplest case incorporating a reversible reaction
Let us consider a loop with three species with one repression and two activations. One of the activations represents a dimerisation reaction, which is reversible, so there is also a reverse positive link from the dimer (A) to the monomer (C). Crucially this gives us a unique steady state. If the third species (B) is supposed to be the corresponding mRNA which is transcribed to form C, then this can be considered a model of negative auto-regulation of a gene.
The dynamical system is given by where f is a decreasing function and g is an increasing function. Typically, models of gene regulation employ linear functions for g , assuming that protein is produced at a rate proportional to the level of mRNA. k denotes the forward reaction rate for the dimerisation reaction and l the reverse rate. d a , d b and d c are the degradation rates of dimer, mRNA (or whatever this species is) and monomer respectively. Let us denote the steady state of the system of Eq. (10) by ( a * , b * , c * ) and for the sake of brevity, denote f ( a * , b * , c * ) by f and g ( a * , b * , c * ) by g . The characteristic equation is given by At the Hopf bifurcations, λ = iα. Taking imaginary parts: Taking real parts gives: Varying d b , the modulus of f g is minimised when Using this optimal value of d b , we get This implies We can see that this is minimised when l / d c and d a / d c are as small as possible. In the limit as they tend to zero, we get Therefore to optimise the chance of oscillations we want the dimer to be as stable and long lived as possible, whilst the rate at which monomers produce the dimer should be much larger than the dimeric degradation rate, but fourfold smaller than the monomeric degradation rate. At steady state, kc 2 * = (l + d a ) a * . In order that l + d a << 4 kc * , we require that c * < < a * , so at steady state, the majority of the molecules should exist in the dimerised state.
Using d c = 4 kc * in the optimal value for d b , we get d b = (2 d a + l) d c . This means that the fastest degradation rate is of the monomer, followed by the other species (which might be the corresponding mRNA), followed by the dimer. The frequency of the oscillations at the bifurcation will be 2 d b d c .
We illustrate the results with numerical simulations for an example system. For the sake of simplicity, we consider that g is linear and f is piecewise linear. The function f takes value f max for a < 0.9, it takes value 0 for a > 1.1 and between these values it is linear, decreasing from f max to 0. This means that f ( a ) only takes two values: 0 or −5 * f max . g(b) = mb, so that the production rate of the monomer is proportional to the level of species B , which might, for example represent the corresponding mRNA, whereby this is a model of auto-repressive gene regulation. We If d c = 4 kc * = 10 , this gives ideal conditions for oscillations, which will occur approximately if | f g | > 8. If we set k = 25 , this is approximately true. As d c diverges from 4 kc * , oscillations become progressively less likely. We therefore fix parameters as mentioned, set f max = 1 . 0 , m = 5 . 0 and vary d c .
In Fig. 2 , we show numerically that indeed for intermediate values of d c there are oscillations. There are no oscillations when d c is too small or too large. It appears that there are two Hopf bifurcations at d c ≈ 1.25 and d c ≈ 82. Whether these are sub-or supercritical warrants further investigation (see Discussion).
Including a reversible reaction substantially changes the conditions for oscillations to occur. As before, they occur when the regulation functions are sufficiently steep at the steady state, but this time they are most likely to occur when the degradation rates differ between species.
Couching this in terms of the general system, we can no longer assume that the angles in the loop are the same, since one of the gradients is 2 kc * (see Discussion). Thus we must take the slightly more general form with three angles θ a , θ b and θ c , where tan θ a = α l+ d a The equation for tan θ b gives This implies: which is the same result as above ( Eq. (12) ).
If we have two species involved in a negative feedback loop, e.g. when the protein corresponding to a gene represses the gene's transcription, then in the simplest well-mixed model, there can be no oscillations, since there are only two species. If, however, the protein dimerises before acting as a repressor, there can be oscillations, since the number of species is now three. Two-species systems can show oscillations in models in which there are delays, see Monk (2003) , Jensen et al. (2003) , Lewis (2003) , and Momiji and Monk (2008) . It is possible to think of the two-species delay system as a higher dimensional system without delays. Delays may be due to the time taken for transcription, for example. In that case, it would alternatively be possible to model the levels of transcripts of different lengths and not explicitly use delays. Similarly, oscillations are possible in models in which there are two chemical species and diffusion ( Glass and Kauffman, 1972;Shymko and Glass, 1974;Sturrock et al., 2011;. However, in the papers by Sturrock and co-workers, the number of mathematical species is at least four, since the protein and mRNA have a nuclear and a cytoplasmic pool. In Glass and Kauffman (1972) and Shymko and Glass (1974) a two-species system is given which exhibits oscillations when there is diffusion, but the authors demonstrate in the appendix of Glass and Kauffman (1972) that a simpler ordinary differential equation model with only two-compartments shows the same behaviour. We also note that if we explicitly model the state of the enhancer or promoter which controls the gene's expression, then this introduces a third species into the model of the single auto-repressive gene. It is common to assume that binding and unbinding to the DNA is fast compared to the dynamics of mRNA and protein. Assuming the binding is at quasi-equilibrium reduces the dimension of the resulting system, so that in the well-mixed system without delays, oscillations are impossible. However if DNA binding and unbinding occur on a timescale similar to transcription, translation and decay of protein and mRNA, then the system is three-dimensional and oscillations are possible ( François and Hakim, 2005;Morant et al., 2009 ). Stochastic oscillations can even occur in such a threedimensional systems when the regulation of transcription is linear, which would preclude oscillations in the deterministic system ( Wang et al., 2014 ).  1 . 01 , 1 . 02 , . . . , 2 . 00 . In (c), we zoom in on the second bifurcation, using values d c = 70 . 1 , 70 . 2 , . . . , 85 . 0 . In (d), we show the period of the oscillations against d c near the first bifurcation. In (e), we show the period of the oscillations against d c near the second bifurcation. We compute the period by identifying timepoints, in [50 0,10 0 0], at which the concentration of monomer is less than its mean concentration (over [50 0,10 0 0]), but at the next timepoint that is no longer true. We then divide the time difference between the first and last such timepoints by the number of such timepoints. In (d)-(e) we mark with a red circle the value of d c at the bifurcation and the period of oscillations there, as predicted by linear analysis. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The ACDC network and other three-species systems with a loop and a sub-loop
The ACDC network is a simple gene circuit consisting of a negative feedback loop between three species, with one reverse link, making a positive feedback loop between two of the species. It is hence a superposition of a repressilator and a toggle switch ( Panovska-Griffiths et al., 2013;Perez Carrasco et al., 2018 ). It can display either switches between expression of the three genes as an input parameter (e.g. a morphogen signal) is changed or oscillations. For some values of the parameters, a stable steady state and a stable limit cycle can coexist. When there are stochastic fluctuations, this can lead to spontaneous switching on and off of oscillations. It has been proposed that the gap gene system, which patterns the anterior-posterior axis in the Drosophila melanogaster embryo, consists of three linked ACDC circuits, two of which operate in the DC (switch-like) regime and one in the AC (oscillatory) regime ( Verd, 2016 ). More complex networks showing multiple negative and positive feedback loops have been proposed to control circaidian rhythms in Drosophila ( Leloup and Goldbeter, 1998 ).
This simplest case of a system with a loop and a sub-loop represents the ACDC network if G is positive.
Let t = tan θ and λ = tan φ where d in is the degradation rate of species in the sub-loop and d out is the degradation rate of species not in the sub-loop. In the notation of Section 3 . and Now we have G − 1 + t 2 = 2 λ , so t 2 = 1 − G + 2 λ . Therefore This is minimised at 4[1 + √ G < 1. This means that when G > 0, i.e. the second loop represents a positive feedback loop, λ > 1, so that the degradation of species in the second loop is faster than that of those outside the second loop. In the case, G > 1, we get λ < 2 / (G − 1) and hence | F | > G − 1 . | F | is minimised at G − 1 when λ = 2 G −1 . Therefore the optimal value of λ is > 1 provided G < 3. Therefore with a positive feedback sub-loop, oscillations are favoured when the degradation of the species in the sub-loop is faster that of those outside the loop, provided the feedback loop is not too strong. If it is very strong then oscillations are favoured when the degradation rate of the species in the sub-loop is slower than those outside. The oscillations that will form when the ratio of the degradation rates is close to 2 G −1 , have very low frequency. When 0 < G < 1, the oscillations have an angular frequency which decreases as G increases from √ 3 d in (when G = 0 , optimally d in is the same as the degradation rate of the species outside the sub-loop).
By contrast, when G < 0, i.e. the sub-loop represents a negative feedback loop, λ < 1, so that the degradation of species in the subloop is slower than that of those outside the sub-loop. In terms of angular contribution to the oscillations, those species involved in coherent sub-loops (which have the same sign as the outer loop) contribute a greater angle than those outside the sub-loop, whereas those involved in incoherent loops have a lesser contribution than those outside the sub-loop. An exception is when there is a very strong incoherent loop. We note that for G < 0, the angular frequency of oscillations decreases when compared to the degradation rate of the species outside the sub-loop, until it tends to a value equal to that degradation rate. The degradation rate of the species within the sub-loop becomes much smaller, so the oscillations increase in frequency compared to it. We note that in the case G < 0, the steady state remains unique. In the case G > 0, there may be more than one steady state.
The value of | F | at the bifurcation is minimised when G is slightly greater than 1. Thus, counterintuitively, the probability of oscillations is maximised when there is a positive feedback loop embedded in the negative feedback loop. The ACDC network gives rise to oscillations more readily than the repressilator, provided that G < 9. By contrast, when there is an embedded negative feedback loop, oscillations are less likely than in the repressilator. For 1 < G < 3, the probability of oscillations is greater than for 0 < G < 1.
We illustrate this in Fig. 3 by showing the results of numerical simulations. We plot the maximal amplitude of oscillations (maximum concentration minus minimum concentration in the time window [50 0,10 0 0] for a gene regulatory network with a negative feedback ring (with repression functions equal to 2 for x < 1/2, 3 − 2 x for 1/2 < x < 3/2 and 0 for x > 3/2) and a second feedback loop (with a reverse repression or activation given by Max (1 − k + kx, 0) , which multiplies the other input to that species), see Fig. 3 a. This is the ACDC motif, if the reverse link is a repression (see Fig. 1 c). Thus | F | = 8 and G = −2 k takes various values. We use the same initial condition, x 1 = 2 , x 2 = 1 . 1 and x 3 = 1 . 1 . We see that sustained oscillations are possible for G > 0. We also plot the period of oscillations against G ( Fig. 3 b) and show for comparison the analytically derived values of G at the bifurcation and the analytically computed period 2 π / α. Note: We find that for all values of G < 0 used in Fig. 3 , the oscillations in simulations are actually decaying, but when G is just below zero, the amplitude is still significant at t = 500 .

n -Species case
Returning to the general case. We again want to minimise the value of | F | at the bifurcation keeping G fixed.
When G is large, in order for | F | to be reasonably small, we require from Eq. (4) , G ≈ (1 + i tan θ ) | W | , so that sin | W | θ ≈ 0 and therefore | W | θ ≈ r π and sec | W | θ (−1) r ≈ G . This implies that | sec θ | is very large, so that θ ≈ π /2, which means | W | = 2 r. This only works if (−1) r = sgnG . If this is true, tiny changes in θ can lead to large changes in the argument of −G + (1 + i tan θ ) | W | . Therefore the argument of (1 + i tan φ) n −| W | is relatively unconstrained and, in order that | F | is as small as possible, sec φ should be minimised. Therefore φ should be small. Thus if G is large, and | W | is 2 r , where (−1) r = sgnG, the degradation rates of species within the sub-loop should be small and those of species outside the sub-loop large. The system is thus approximately equivalent to the sub-loop only, with the species outside it in quasiequilibrium with those inside. We have said that for the species in the loop, θ should be close to π /2, which means that their degradation rates should be close to zero. This, however, is at the Hopf bifurcation. What will happen is that, when the sub-loop has a very large negative coefficient, the system will oscillate for most val-ues of the degradation rates. When the degradation rates are very small, oscillations will cease. For information on what happens when there is a strong positive feedback loop, see Strelkowa and Barahona (2010) .

Two interconnected rings
Many cases of two interconnected rings can be viewed as a large ring with a single sub-ring, which we have already studied. A simple example, which cannot, is obtained by extending the ACDC network, such that the reverse link added to the repressilator forms another repressilator with a fourth gene (see Fig. 1

d).
Although there is an outer ring ( ACDB in the figure) and a subring ( BC in the figure), there are more sub-loops to consider, corresponding to the three-species repressilators. The equations can be written as: where all functions f, g, h and j are decreasing in their arguments.
We note that this system can have multiple steady states. In the example version that we use for numerical simulation, below, there are five steady states.
The characteristic equation at any of these steady states is given by where all derivatives are evaluated at the relevant steady state. Therefore 1 We see that this is symmetric with respect to d b and d c , so we assume the optimal condition for oscillations has d b = d c .
Without loss of generality, we set d c = d b = 1 , and substituting in we get: If we further assume that the two negative feedback loops have equal coefficients ( f h a g c = j g d h b ) , then species a and d are involved in the characteristic equation in symmetric ways, so we assume the probability of oscillations will be optimised if d a = d d = δ. Setting f h a /h b = −L = j g d /g c , we therefore get: Taking the imaginary part gives and observing that the left hand side of the equation is positive gives: L + Lδ − δ. If δ is very small, then the angular frequency of oscillations at the Hopf bifurcation will be approximately √ Lδ ( Ld a d b in dimensional terms), which is intermediate between the degradation rates of the fast and the slow species.
In Fig. 4 , we show numerically that in a simple example sys- there are oscillations when δ is small, but these cease when δ gets big enough.
This system exhibits multi-stability because of the existence of positive feedback loops as well as negative feedback loops. Therefore a more careful investigation of its bifurcations is warranted, see Discussion. In addition, the oscillations appear to occur around steady states which break the symmetry of the equations, i.e. the steady state values of A and D are different, as are those of B and C . It is therefore difficult to impose the symmetry condition that the repressilator coefficients are equal. A piecewise linear model could be employed in order to create a direct numerical comparison with the analytic results. This is beyond the scope of the current study. Nevertheless, simulations do show that for high δ oscillations are switched off. This is to be expected, since the system is then effectively a bistable switch between species B and C . Other parameters are as described in the main text. We solve the ODEs using the Matlab ode solver ode23s. We use initial conditions, a = 1 .
The structure of this four-species model is similar to a model of the p53-Mdm2 network response to DNA damage presented in Abou-Jaoudé et al. (2009) andOuattara et al. (2010) . The difference there is that the links from C to D and from D to B ( Fig. 1 d) are positive (see Fig. 1 of Abou-Jaoudé et al., 2009 ) and their direction is reversed relative to the model presented here.

Conclusions and discussion
In this paper, I have presented a systematic way ( Eq. (2) ) of organising the characteristic equation for the eigenvalues of the Jacobian matrix which determine the stability of steady states in regulatory networks. This facilitates the identification of Hopf bifurcations. I have used this to determine optimal conditions for oscillations in regulatory networks.
The chance of oscillations in a ring oscillator is maximised if the degradation rates of the species are equal . This is independent of all details of the negative feedback loop. When the system has more than one feedback loop the situation is more complicated. For example, when there are two loops, one of which is negative and contains all the species, then if the other loop is also negative, the optimal degradation rates of the species within it are smaller than those of the species outside it. If the second loop is a positive feedback loop, then, provided it is not too strong, the opposite is true and species outside the loop should degrade more slowly.
Counterintuitively, adding a positive feedback loop to a repressilator increases the probability of oscillations, provided that the loop is not too strong, whereas adding another negative feedback loop decreases the probability of oscillations. The fact that positive feedback loops could enhance oscillations in a repressilator was already known for specific kinetic functions, see for example ( Ananthasubramaniam and Herzel, 2014 ).
In Section 3.1 , it is shown that the optimal conditions for oscillations for the autoregulatory gene model have the majority of the protein in its dimerised form at the steady state. The dimeric degradation rate and the rate at which the dimer dissociates to form monomers should both be small. The degradation rate of the mRNA should optimally be intermediate and the monomeric protein degradation rate large. This is unlikely to be typical of autoregulatory genes, since mRNA degradation rates are typically larger than protein degradation rates (e.g. Hoffmann et al., 2002;Nelson et al., 2004 ) and may suggest that oscillations are most likely when some process explicitly removes monomeric protein.
Degradation rates of Hes1 protein and mRNA are however similar (e.g. Jensen et al., 2003;Hirata et al., 2002 ). Optimally d c = 4 kc * , which means that, at the steady state, only one in every three monomers dimerises whilst the other two degrade (see Eq. (10 )).
The optimal mRNA degradation rate is then (2 d a + l) d c . A simple closed form expression in terms of the parameters of the model is given for the angular frequency of the oscillations close to the bifurcation ( Eq. (11) ). This simplifies under the optimal assumptions listed above to α = 2 d b d c . Thus, under these conditions, the period of oscillations is intermediate between the timescale of degradation of the mRNA and the protein monomer. If, instead of assuming optimal conditions for oscillations, we simply assume that the dimeric protein does not degrade directly, then the angular frequency at the bifurcation is given by Again assuming the optimal condition for one in every three monomers to form a dimer, gives α = d c (2 d b + l) and assuming that the dimeric unbinding rate is small compared to the mRNA degradation rate yields the same angular frequency as before of α = 2 d b d c , without assuming that the mRNA is longer lived than the monomeric protein. If the dimer less stable, this should speed up oscillations. A dimeric model of autoregulation of Hes1 with protein and mRNA degradation rates of approximately 1/(25 mins.) would predict an oscillation period at the bifurcation of 111 mins., assuming one in every three monomeric proteins dimerises and the unbinding rate of dimers is very low. This is close to the actual period (see e.g. Yoshiura et al., 2007 ). Much more rapid oscillations would be possible if the dimer were less stable. See Momiji and Monk (2008) and Sturrock et al. (2014) for more detailed models incorporating nuclear and cytoplasmic components and transcriptional and translational time delays. Whilst the detailed features of sub-cellular location are no doubt important, an angular frequency of √ 2 times the geometric mean of the mRNA and protein degradation rates may serve as a useful estimate.
In Section 4 , we predict that for two linked repressilators (or equivalent negative feedback loops), see Fig. 1 d, of equal strength (i.e. feedback coefficient), conditions for oscillations are optimised if the species forming the two-species positive feedback loop ( B and C in the figure) degrade faster than the other species. A and D should have equal degradation rates, as should B and C . If B and C degrade very fast compared to the other two species, then the angular frequency of oscillations at the bifurcation will be α = Ld a d b , where −L is the coefficient of the negative feedback loop divided by the coefficient of the two-species positive feedback loop. This warrants further investigation, see below.
The dynamics of systems with more than one feedback loop are more complicated than those with a single feedback loop. There can be multiple stable steady states (as is true for a ring with a positive feedback sub-loop). There can be the co-existence of a stable steady state and a stable limit cycle, for example for the simple ACDC circuit, Perez  . Finally there can be chaos, Zhang et al. (2012) , which is not possible for the negative feedback ring systems, Page and Perez Carrasco (2018) .
In this paper, we only consider ordinary differential equation models of regulatory networks. Much interesting work has shown how oscillations can arise when there are only two biological species, but when there are delays ( Lewis, 2003;Monk, 2003;Jensen et al., 2003;Momiji and Monk, 2008 ), multiple compartments ( Glass and Kauffman, 1972 ) and/ or diffusion ( Busenberg and Mahaffy, 1985;Shymko and Glass, 1974;Cangiani and Natalini, 2010;Sturrock et al., 2011;. As mentioned in Section 3 , these cases are equivalent or similar in nature to ordinary differential equation models with a higher number of mathematical species. Nevertheless, sub-cellular location of biochemical species is clearly an important factor in the regulation of gene expression (e.g. Hoffmann et al., 2002;Nelson et al., 2004 ). The environment of the cells is also noisy and so stochastic effects are important. For examples of numerical algorithms to simulate stochastic gene regulation and the implications of intrinsic noise in negative feedback systems, see Hegland et al. (2007) and Zeron and Santillán (2010) .
Further work could address the inclusion of nonlinear degradation or other nonlinear negative self-regulatory terms for the species. This would allow the auto-regulatory gene model to be couched in terms of the general theory of Section 2 and would also facilitate application to a wider range of real gene regulatory models. Characterising the Hopf bifurcations as supercritical or subcritical would also be interesting. This is somewhat laborious, however, even for the case of a ring oscillator (see Page and Perez Carrasco, 2018 ). A better numerical bifurcation analysis of the fourspecies model and similar models, such as that of Abou-Jaoudé et al. (2009) and Ouattara et al. (2010) , would also be desirable. The multiple steady states render finding Hopf bifurcations more difficult in this model. Finally, future work could extend the results on optimal conditions for oscillations to systems with delays, diffusion and stochasticity.