Bull Math Biol DOI 10.1007/s11538-015-0073-9 ORIGINAL ARTICLE

A Modelling Framework for Gene Regulatory Networks Including Transcription and Translation R. Edwards · A. Machina · G. McGregor · P. van den Driessche

Received: 2 June 2014 / Accepted: 2 March 2015 © Society for Mathematical Biology 2015

Abstract Qualitative models of gene regulatory networks have generally considered transcription factors to regulate directly the expression of other transcription factors, without any intermediate variables. In fact, gene expression always involves transcription, which produces mRNA molecules, followed by translation, which produces protein molecules, which can then act as transcription factors for other genes (in some cases after post-transcriptional modifications). Suppressing these multiple steps implicitly assumes that the qualitative behaviour does not depend on them. Here we explore a class of expanded models that explicitly includes both transcription and translation, keeping track of both mRNA and protein concentrations. We mainly deal with regulation functions that are steep sigmoids or step functions, as is often done in protein-only models. We find that flow cannot be constrained to switching domains, though there can still be asymptotic approach to singular stationary points (fixed points in the vicinity of switching thresholds). This avoids the thorny issue of singular flow,

This work was partially supported by Discovery Grants from the Natural Sciences and Engineering Research Council (NSERC) of Canada. R. Edwards (B) · A. Machina · P. van den Driessche Department of Mathematics and Statistics, University of Victoria, STN CSC, PO Box 1700, Victoria, BC V8W 2Y2, Canada e-mail: [email protected] A. Machina e-mail: [email protected] P. van den Driessche e-mail: [email protected] G. McGregor Department of Mathematics and Statistics, McGill University, 805 Sherbrooke Street West, Montreal, QC H3A 0B9, Canada e-mail: [email protected]

123

R. Edwards et al.

but leads to somewhat more complicated possibilities for flow between threshold crossings. In the infinitely fast limit of either mRNA or protein rates, we find that solutions converge uniformly to solutions of the corresponding protein-only model on arbitrary finite time intervals. This leaves open the possibility that the limit system (with one type of variable infinitely fast) may have different asymptotic behaviour, and indeed, we find an example in which stability of a fixed point in the protein-only model is lost in the expanded model. Our results thus show that including mRNA as a variable may change the behaviour of solutions. Keywords Gene regulation · Piecewise-linear · Singular perturbation · mRNA-protein model · Transcription–translation model Mathematics Subject Classification

92C42 · 34A36 · 92C40

1 Introduction Despite tremendous advances in technology in recent years, modelling of gene regulatory networks is still hampered by the difficulty of determining rate constants, so that quantitative models are unreliable unless robust to parameter variation. Another issue is that the principles governing the behaviour of networks in general are not necessarily discoverable from any one particular network, no matter how well its fine details are described. These concerns have led to the development of qualitative models of gene networks, in which an attempt is made to assess qualitative behaviour somewhat independently of details, such as precise parameter values and biochemical kinetics. Thomas (1973), Glass and Kauffman (1973) and others began to explore this approach in the 1970s, using simplified model frameworks with the aim of capturing important qualitative features of the behaviour of gene networks and revealing fundamental principles of their design. Glass and Pasternack (1978), for example, showed that a class of networks that is broader than, but includes, negative feedback loops must show oscillatory behaviour, either damped or sustained according to a condition that could be calculated. A considerable literature has been built up over the intervening years, developing and elaborating this theory. Much of this work made some fairly drastic simplifying assumptions, including: – modelling only a single product of each gene, which is described as directly acting as a transcription factor in the regulation of other genes, – a step function limit of sigmoidal interactions, based loosely on the Michaelis– Menten and Hill function kinetics that arise in modelling enzyme-mediated or catalytic reaction kinetics, such as occurs in gene expression, – a single effective threshold of each gene product concentration, even when it regulates multiple other genes, – equal degradation rates of all gene products, – no autoregulation. The class of systems in which at least the first two of these assumptions are made, without necessarily requiring the others, have been called ‘Glass networks’. These simplifying assumptions allowed some interesting mathematical work, which laid a

123

Transcription–Translation Networks

foundation for the qualitative study of gene network dynamics. Of course, it is hoped that the insights gained still hold when some or all of these simplifying assumptions are relaxed, and more realistic models are used. Developments in synthetic biology has made this work all the more compelling and important. For instance, the groundbreaking papers by Elowitz and Leibler (2000) on the ‘Repressilator’ and by Gardner et al. (2000), on the genetic toggle switch (both published in the same issue of Nature), showed that simple networks built in the laboratory did, indeed, demonstrate the fundamental behaviours predicted by qualitative models. In more complex network structures, however, the correspondence is less clear-cut and considerable effort has gone into extending the analysis of the simplest class of qualitative models to more general classes in which one or more of the assumptions above are relaxed, though the mathematics becomes more difficult. Multiple thresholds are the least problematic, and a number of studies have been able to effectively extend results to that situation (e.g. Farcot 2006). In fact, it has been shown that, mathematically, such networks can be embedded into a higher dimensional system with a single threshold per variable in such a way that the dynamics are essentially equivalent (Edwards et al. 2012). Unequal degradation rates do not cause a real problem for local computation of solution trajectories, but make identification of stable periodic solutions more difficult, though some results along those lines have now been obtained (Farcot and Gouzé 2010; Edwards and Ironi 2014). Allowing effective autoregulation means that singular situations have to be resolved when gene product concentrations become confined at or near a threshold value for an interval of time. This may occur because a single transcription factor simply inhibits its own production, but may also occur when pairs of variables have damped oscillations about a threshold intersection, such that in the step function limit, they converge to the threshold intersection in finite time (Edwards 2000). Such finite time convergence again requires autoregulation. These issues of singular behaviour seem, in one sense, to be mathematical artifacts, since we know that no real regulatory function has infinite slope. However, it has been shown that interesting behaviour of sigmoidal systems may be exposed by this approach (Machina et al. 2013b). There has been a controversy about how to handle these singular situations. Gouzé and colleagues used Filippov theory for discontinuous vector fields to produce set-valued solutions (Gouzé and Sari 2002; Casey et al. 2006), while Plahte and colleagues developed an approach based on singular perturbation theory (Plahte and Kjøglum 2005; Ironi et al. 2011). Both methods had their problems. The usual Filippov approach produces spurious solutions in some situations, while the classical singular perturbation theory can only be applied in a restricted class of problems. Recent work by Machina and colleagues (Machina and Ponosov 2011; Machina et al. 2013a, b) has clarified the relationships between these approaches and removed some of the difficulties. There exists a more appropriate definition of Filippov solution that avoids most of the spurious solutions and corresponds more closely to the singular perturbation solutions. Furthermore, there exists a generalization of the classical singular perturbation theory that avoids its restricted application. However, in the latter case, non-uniqueness can still arise in the limit of step function interactions, much like in the Filippov theory, though now it can

123

R. Edwards et al.

be seen to arise from densely interwoven basins of attraction of different attractors and sensitive dependence on initial conditions in systems with steep sigmoidal interactions. The singular perturbation approach to the handling of singular dynamics, in principle, obtains results about the behaviour of steep sigmoidal systems, using the step function limit only as a baseline from which to perturb. The work of Ironi et al. (2011), for example, allows a bound to be placed on the steepness parameter, above which the qualitative behaviour predicted must exist. Even here, though, the sigmoid interactions often have to be very steep, in biological terms, for the results to hold rigorously. There have also been numerical investigations of how behaviour changes as the steepness parameter is altered (Lewis and Glass 1991, 1992; Killough and Edwards 2005; Wilds and Glass 2009). Typically, if the step function network has complex behaviour, that behaviour is simplified by making the interactions less steep. Very little work has been done in the qualitative model literature, however, on relaxing the first assumption above, and including in the model multiple interaction steps, a notable exception being a recent paper of Gedeon et al. (2012). Of course, in more realistic biochemical kinetic models, several steps are involved in producing a transcription factor that can play a regulatory role. The expression of a gene comprises transcription of copies of an mRNA molecule, characteristic of the gene, whose code is then translated into a protein, by building a corresponding string of amino acids, which are available in the cell. This protein may itself act as a transcription factor for another gene, but it may also need to undergo further modifications, such as dimerization, phosphorylation or other post-translational reactions, before it is effective as a regulator. In qualitative models, all of this is generally suppressed, and the situation is treated as if each gene product (protein), or combination of gene products, acts directly on the gene products of other genes by modifying their expression (by either repression or activation). The singular situations, including autoregulation, that arise in these simplified models, of course, hide the intermediate steps in every regulation process. Thus, another way to approach the problem of autoregulation, or fast damped oscillation in negative feedback between two or more genes, which lead to the singularities in the step function limit of the simple qualitative models, is to expand the framework of the qualitative models to account for multiple steps in the regulatory process. For example, one can model the transcription of mRNA and then translation into protein as two separate steps, and one can also model post-translational processes as separate steps. Then, even when a gene regulates itself, the autoregulation is split into multiple steps. No variable appears in its own production term, and in that sense, no variable is directly self-regulating, and one might expect that the singularities, although present in the step function limit, should only be passed through transparently (transversally) or approached asymptotically. Although this makes the models more complicated, it may avoid the somewhat artificial mathematical difficulties associated with the singular dynamics and may still be mathematically tractable, while being truer to the biochemistry. One question that arises is whether the dynamics of these models converges to the dynamics of the simpler one-step models in the limit where all but one of the multiple steps are infinitely fast. If so, this provides at least some justification for the protein-only models. Another question that arises is how to analyse the new type of model. From observations, the timescales of the dynamics of mRNA and protein

123

Transcription–Translation Networks

concentrations are typically different, but the magnitude of this difference depends on the gene and organism. For example, Bernstein et al. (2002) find mRNA half-lives of 3–8 min for most genes in E. coli, and Mosteller et al. (1980) found protein half-lives of 2–23 h or more, also in E. coli. However, in humans, mRNA half-lives can be on the order of 2–10 h and protein half-lives from 45 min to 28 h (see BioNumbers website). Thus, the ratio of decay rates of mRNA to protein is not at all clear and can range from at least 100 to much smaller values. Gedeon et al. (2012) use a value of 10 as biologically reasonable. Paetkau et al. (2006) cite particular examples of mRNAs and proteins in particular organisms where the ratio may be close to unity, even though such cases are recognized as exceptional. Note, furthermore, that the production rates of mRNA and protein are necessarily more closely linked than the decay rates, since the protein is produced from the mRNA by translation. Here, we take the approach of accounting for two steps in each regulation, the transcription and translation, thus keeping track of both mRNA and protein concentrations. The transcription is still governed by sigmoid functions, while the translation is governed by a linear term in the differential equations. This is the model expansion considered by Gedeon and co-authors, but the conditions imposed on the general results in (Gedeon et al. 2012, Section 3.2) are restrictive. We discuss the relation between our results and theirs in Sect. 8. Post-transcriptional modifications can take many forms. Dimerization, for example, requires a quadratic term in the model equations, while an enzyme-mediated or catalytic reaction requires a Michaelis–Menton or Hill Function (sigmoid) interaction term. Since these multiple possibilities make things more complicated, we do not consider them in detail here. However, in an Appendix we consider another mathematical possibility for avoiding the singular flow: a two-step system in which both steps are mediated by sigmoid functions, thus giving a larger sigmoidal network, but one in which no variable is directly self-regulating. This might be considered to model the expression of a protein as one step and the catalytic modification of that protein as a second step, though this is only one type of reaction that can occur in gene networks. We can also consider it as a mathematical method for dealing with the singular flow in a traditional gene network model by considering an infinitely fast limit of one of the variables in each pair. In Sect. 2, we summarize the traditional protein-only model using the description of Plahte and Kjøglum, and we briefly summarize the singular perturbation analysis (Plahte and Kjøglum 2005). Then in Sect. 3 we introduce the two-step extension to make it an mRNA-protein model and give a result on the singular stationary points (SSPs) in this expanded model. In Sect. 4, we show that solutions of the original system are recovered if we take an infinitely fast limit in the rate of the protein equations. In Sect. 5, we do the same for infinitely fast mRNA rates. In the more sensitive case of SSPs in domains where two variables are switching simultaneously, we show in Sect. 6 that such points are generally unstable. Then in Sect. 7, we determine the flow from one switching point to the next, to allow integration of trajectories in regular domains. In the Discussion, we deal with the relationship between our results and those of Gedeon et al. (2012), as well as the relationship between our model and the differential delay systems of Ponosov (2005) and Shlykova et al. (2008). We also

123

R. Edwards et al.

summarize the results obtained and questions remaining unanswered. In an Appendix, we consider an alternate two-step model where both steps involve sigmoidal regulation. 2 Preliminaries We are concerned here with a class of qualitative protein-only gene network models as discussed in the introduction, which has been formalized in a number of ways. Using the notation of Plahte and Kjøglum (2005), the network equations are y˙i = Fi (Z ) − γi yi ,

i = 1, . . . , n,

(1)

with yi the concentration of the protein product of gene i, γi > 0 and Z = (Z 1 , . . . , Z n ) an n-dimensional vector of sigmoid functions Z i = S(yi , θi , q). In general, S(yi , θi , q), i = 1, . . . , n, must satisfy a number of conditions listed in (Plahte and Kjøglum 2005). For example, these conditions are satisfied by the Hill function given by 1/q

Z i = H (yi , θi , q) =

yi 1/q

yi

1/q

+ θi

,

i = 1, . . . , n,

0 < q ≤ 1.

(2)

Here, q is the steepness parameter (smaller values of q correspond to greater steepness of the sigmoid function) and θi is called the threshold, because H (yi , θi , q) ∈ (0, 1) and H (θi , θi , q) = 21 . We take S = H in the main body of the paper, but a different sigmoid function arises in the Appendix. The production rate function Fi ≥ 0 is a multilinear polynomial, i.e., an affine function with respect to each Z i . Production rates are inherently bounded, so that 0 ≤ Fi (Z ) ≤ F i , i = 1, . . . , n, where each F i , is a constant. Plahte and Kjøglum allow the degradation rate to also depend on concentrations of transcription factors so that the model equations become y˙i = Fi (Z ) − G i (Z )yi ,

i = 1, . . . , n.

Although this does not make mathematical analysis more difficult, in practice degradation rates can normally be considered constant, and G i (Z ) = γi makes our exposition more straightforward. The regulatory effect of the transcription factors occurs only through the production rates. An alternative to the model given by (1) and (2) is to approximate steep sigmoidal functions Z i by step functions, also denoted by Z i , with the jump at θi . System (1) then becomes a switched system with the phase space partitioned into regular domains by the threshold hyperplanes (walls) yi = θi . Within each of the regular domains, (1) is linear, since Z is a binary vector that has a constant value in the domain. System (1) has a unique ‘focal point’ for the value of Z appropriate to the domain, the point at which the right-hand side functions are all 0. If this focal point lies in the domain itself, it is an asymptotically stable node to which solutions converge, but if not then solution trajectories hit a threshold hyperplane or intersection of hyperplanes before

123

Transcription–Translation Networks

reaching the focal point. On the other side of this threshold or thresholds, the value of Z is different and the focal point is therefore different in general. Whenever we deal with the functions Z i as step functions (whether as limits of sigmoids or as an approximation), we will make the following assumption to avoid technical difficulties that only arise for very special sets of parameters (on a set of measure 0 in parameter space). Assumption 1 No focal point, , where i = Fi (Z )/γi for any binary vector Z , lies on a threshold hyperplane. Thus, i = θi for any i. Instead of simply replacing the sigmoid with a step function, we can attempt to deal with steep sigmoids by means of a singular perturbation from the step function limit. We will refer to a singularly perturbed system several times in the rest of this paper, not only in relation to system (1). In general, such a system involves coupled slow and fast motion, say d X i = f (X, Y ), i dt i = 1, . . . , n, (3) d Yi ε dt = gi (X, Y ), where X = (X 1 , . . . , X n ), Y = (Y1 , . . . , Yn ), and ε > 0 is a small parameter. The initial value problem is determined by initial conditions X (0) = X 0 , Y (0) = Y0 . The variables X and Y are referred to as the slow and fast variables, respectively. The solution to the initial value problem (3) is denoted by (X ε (·), Y ε (·)). The standard analysis based on Tikhonov’s theorem (see e.g. Artstein 2002) guarantees under certain conditions that the solution (X ε (·), Y ε (·)) converges as ε → 0 to the solution of the differential-algebraic system d X i = f (X, Y ), i dt 0 = gi (X, Y ),

i = 1, . . . , n,

obtained from (3) when ε = 0. In particular, a crucial condition is that for each fixed X, solutions of the differential equation dYi = gi (X, Y ), dτ

i = 1, . . . , n,

representing fast flow in a new time variable, τ = t/ε, should converge as τ → ∞ to a solution of the algebraic equation 0 = gi (X, Y ). If the fast system does not converge to a fixed point, then the Tikhonov theorem does not apply and we have to resort to a more general theory. It is useful to have a relatively simple example of our class of gene networks to which we can refer in later sections. Example 1 (Plahte and Kjøglum 2005) y˙1 = Z 1 + Z 2 − 2Z 1 Z 2 − γ1 y1 , y˙2 = 1 − Z 1 Z 2 − γ2 y2 ,

(4)

123

R. Edwards et al.

with Z 1 , Z 2 given by (2). Example parameter values chosen by Plahte and Kjøglum (2005, Fig. 2) are γ1 = 0.6, γ2 = 0.9, θ1 = θ2 = 1. Solutions with small q > 0 (i.e., steep Hill functions) are studied by means of singular perturbation analysis with reference to the limiting system in which q → 0. However, the limit system in the vicinity of threshold values is studied only in a transformed system of differential equations, where the variables near threshold (switching variables) are described by their Z i values and time is rescaled by τ = t/q. In this transformed system, there is no discontinuity. In the limit q → 0, the switching domain {y1 = θ1 , y2 < θ2 } is a ‘white wall’, i.e. the solutions depart from it. The switching domain {y1 < θ1 , y2 = θ2 } is a ‘transparent wall’, i.e. the solutions pass through it. The switching domains {y1 > θ1 , y2 = θ2 } and {y1 = θ1 , y2 > θ2 } are ‘black walls’, i.e. the solutions are attracted to these domains. In the first of these black walls, there exists a unique singular stationary point (SSP), which is at (1.5, 1) with the parameter values chosen by Plahte and Kjøglum and is asymptotically stable. For some parameter values, γ1 , γ2 , θ1 , θ2 , including those chosen above, and small enough q, there is also a stable SSP near the threshold intersection, (θ1 , θ2 ). This behaviour is depicted by Plahte and Kjøglum (2005, Fig. 2).  3 The mRNA-protein Model We now expand the n-dimensional model (1), which accounts only for the protein transcription factors that regulate other genes, to obtain a model of dimension 2n that includes both mRNA and proteins as separate variables: x˙i = Fi (Z ) − βi xi , y˙i = κi xi − αi yi ,

i = 1, . . . , n,

(5)

where yi is the concentration of the protein product of gene i, xi is the concentration of the ith mRNA, Z = (Z 1 , . . . , Z n ) with Z i = H (yi , θi , q), and 0 < q 0, αi > 0, κi ≥ 0. Note that since y˙i is independent of Z i , all switching hyperplanes yi = θi are ‘transparent’ (solutions just pass through switching hyperplanes). For smooth sigmoids, the point ( ακ11 y1∗ , . . . , ακnn yn∗ , y1∗ , . . . , yn∗ ) is a stationary solution of (5) if and only if (y1∗ , . . . , yn∗ ) is a stationary solution of y˙i =

κi αi

Fi (Z 1 , . . . , Z n ) − βi yi ,

i = 1, . . . , n.

(6)

For switched systems, stationary points are traditionally defined as the limit as q → 0 of the corresponding stationary points of the continuous systems with smooth sigmoids. Thus, the same relation holds between the stationary solutions of the switched systems (5) and (6) with the step regulatory functions (q → 0). The problem of locating the stationary solutions of switched systems (5) is then reduced to locating the

123

Transcription–Translation Networks

stationary solutions of (6), which is well studied in the literature (Plahte and Kjøglum 2005; Casey et al. 2006). Stationary points in regular domains are straightforward and are always asymptotically stable. Below we provide a stability result for the system (5) in case of only one variable yi assuming its threshold value. For notational simplicity, we assume that y1 = θ1 but the result is valid for any fixed yi = θi . Theorem 1 Suppose that y ∗ = (θ1 , y2∗ , . . . , yn∗ ), yi∗ = θi , i = 2, . . . , n, is a SSP of (6) in the limit as q → 0. Then (x ∗ , y ∗ ) with x ∗ = ( ακ11 θ1 , ακ22 y2∗ , . . . , ακnn yn∗ ) is an asymptotically stable (respectively, unstable) SSP of (5) in the limit as q → 0 if and only if y1∗ = θ1 is an asymptotically stable (respectively, unstable) solution of (6) with i = 1. Proof We separate the singular variables x1 , y1 and rewrite (5) as the following 2n dimensional system x˙1 y˙1 x˙r y˙r

= F1 (Z 1 , Z 2 , . . . , Z n ) − β1 x1 = f 1 (x1 , y), = κ1 x1 − α1 y1 = f 2 (x1 , y1 ), = Fr (Z 1 , Z 2 , . . . , Z n ) − βr xr = fr +1 (xr , y), = κr xr − αr yr = fr +n (xr , yr ), r = 2, . . . , n.

(7)

We look at the Jacobian matrix J (q) of the system (7) at (x ∗ , y ∗ ) for a small q > 0, written by columns ( ∂∂xf1 , ∂∂yf1 , ∂∂xf2 , . . . , ∂∂xfn , ∂∂yf2 , . . . , ∂∂yfn ) ⎛

∂ f1 ∂ x1 ∂ f2 ∂ x1

⎜ ⎜ ⎜ ⎜ 0 ⎜ ⎜ ··· J (q) = ⎜ ⎜ ⎜ 0 ⎜ ⎜ 0 ⎜ ⎝ ··· 0

∂ f1 ∂ y1 ∂ f2 ∂ y1 ∂ f3 ∂ y1

···

∂ f n+1 ∂ y1

0 ··· 0

0 0 −β2 ··· 0 κ2 ··· 0

··· 0 ··· 0 ··· 0 ··· ··· · · · −βn ··· 0 ··· ··· · · · κn

∂ f1 ∂ y2

0 ∂ f3 ∂ y2

···

∂ f n+1 ∂ y2

−α2 ··· 0

· · · ∂∂ yf1n ··· 0 ∂ f3 · · · ∂ yn ··· ··· · · · ∂∂fn+1 yn ··· 0 ··· ··· · · · −αn

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟. ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

(8)

Note that since yr is away from its threshold, ddZy r , r = 2, . . . , n, vanishes when r q → 0, and, in fact, does so exponentially fast, so the entries in the last n − 1 columns of the first n + 1 rows of the matrix are all negligible. The characteristic equation of (8) is thus    ∂ f1 − λ ∂ f1    ∂ x1 ∂ y1  (−β2 −λ) · · · (−βn −λ)(−α2 −λ) · · · (−αn −λ)+ζ (q) = 0 (9)  ∂ f2 ∂ f2   ∂x − λ ∂y 1

1

where limq→0 ζ (q)/q = 0, since ddZy r → 0 as q → 0 exponentially fast. The terms r of J (q) that approach ∞ as q → 0 are those in the second column. These take the  1 form ∂∂yf1i = ∂∂ ZF1i ddZy 1 , with ddZy 1  y ∗ =θ = 4qθ , so these grow only proportionally to q1 . 1 1

1

1

1

123

R. Edwards et al.

Thus, the terms of the characteristic equation involving ζ (q) become insignificant as q → 0. Also, βi > 0 and αi > 0 for all i = 1, . . . , n. Thus, the stability analysis at (x ∗ , y ∗ ) reduces to studying the 2 × 2 Jacobian

J (q) =

∂ f1 ∂ x1 ∂ f2 ∂ x1

∂ f1 ∂ y1 ∂ f2 ∂ y1



=

−β1 ∂∂ yf11 κ1 −α1

 (10)

with trace( J (q)) = −β1 − α1 < 0 and det( J (q)) = β1 α1 − κ1 ∂∂ yf11 , where  ∂ f1 ∂ F1 d Z 1 d Z1  1 ∗ ∗ ∂ y1 = ∂ Z 1 dy1 at (x , y ). Since dy1 y1∗ =θ1 = 4qθ1 → ∞ as q → 0, the sign of det( J (q)) is determined by the sign of ∂ F1 . Thus, (x ∗ , y ∗ ) is asymptotically stable ∂ Z1

(respectively, unstable) if and only if ∂∂ ZF11 < 0 (respectively, > 0), which is equivalent to (θ1 , y2∗ , . . . , yn∗ ) being an asymptotically stable (respectively, unstable) stationary point of (6) in the limit q → 0.  Remark 1 There exists a sufficiently small q¯ > 0 such that for all 0 < q < q, ¯ the result of Theorem 1 holds qualitatively with an asymptotically stable stationary point, (x ∗ , y ∗ ) with y1∗ near θ1 . Remark 2 The result of Theorem 1 holds for system (5) even when the protein reactions (translation) are much faster than the mRNA reactions (transcription) and vice versa. The situation with infinitely fast reactions is considered below in Sects. 4 and 5. In the next two sections, we assume different timescales for the protein and mRNA reactions. Our results in the extreme case of infinitely fast rates of one of the variable types (protein or mRNA) provide a partial justification for using the protein-only equations (1), but it should be noted that asymptotic behaviour may still be different. In practice, the actual ratio of mRNA to protein rates may not be extremely small or extremely large. Estimates of this ratio in different contexts vary (see, for example, Garcia-Ojalvon et al. 2004, p. 10,958, 2nd paragraph; Gedeon et al. 2012). 4 Infinitely Fast Protein Rates In this section, we consider a situation where the protein translation rates are much faster than the mRNA transcription rates. Such infinitely fast protein dynamics are modelled by introducing a small parameter ε → 0. The singular perturbation theory is then applied with ε as the small parameter. We show that the limit dynamics of the expanded system after a short transient is represented by (1) with an appropriate choice of Fi . Consider the system of dimension 2n as in (5) but with 1ε scaling the rate of the protein dynamics. x˙i = Fi (Z ) − βi xi , i = 1, . . . , n, (11) y˙i = 1ε (κi xi − αi yi ), where Z = (Z 1 , . . . , Z n ), Z i = H (yi , θi , q). We denote by (x ε , y ε ) a solution to (11) for ε > 0. The initial value problem is determined by the initial conditions xi (0) = xi,0 , yi (0) = yi,0 .

123

Transcription–Translation Networks

Let q > 0 be fixed. The system (11) represents a singular perturbation problem as ε → 0 with the fast flow yi and the slow flow xi . The pseudo-stationary solution yi∗ = κi xi /αi of the fast flow is obtained from κi xi − αi yi = 0. It is easy to check that the stationary solution yi∗ is asymptotically stable for each fixed xi . Thus, all assumptions of the Tikhonov theorem are fulfilled, and the slow flow is therefore governed by (12) x˙i = Fi (Z 1∗ , . . . , Z n∗ ) − βi xi , where Z i∗ = H (yi∗ , θi , q) = H (κi xi /αi , θi , q) is the Hill function. Note that due to the properties of the Hill function, H (κi xi /αi , θi , q) = H (xi , θi αi /κi , q). The system (12) can be rewritten as x˙i = Fi ( Z1, . . . , Z n ) − βi xi ,

i = 1, . . . , n,

(13)

where Z i = H (xi , θi , q), with θi = θi αi /κi , or as y˙i∗ =

κi αi

Fi (Z 1∗ , . . . , Z n∗ ) − βi yi∗ ,

i = 1, . . . , n.

(14)

By the Tikhonov theorem, we have the following result (see, e.g., Artstein 2002, Theorem 2.1). Theorem 2 Let q > 0 and an initial condition (x0 , y0 ) be fixed. The solutions xiε (t), i = 1, . . . , n, of (11) converge as ε → 0 to the solution xi (t) of the system (13) uniformly on intervals of the form [0, T ]. The solutions yiε (t), i = 1, . . . , n, converge as ε → 0 to the solution yi∗ (t) of the system (14) uniformly on intervals of the form [δ, T ] for δ > 0. On intervals τ ∈ [0, S] with S > 0 fixed and τ = t/ε, the trajectories yiε (τ ) converge uniformly as ε → 0 to

 κi xi,0 κi xi,0 −αi τ yi (τ ) = + yi,0 − , e αi αi the solution of

dyi (τ ) = κi xi,0 − αi yi (τ ) with yi (0) = yi,0 . Finally, dτ lim lim yiε (εS) = κi xi,0 /αi .

S→∞ ε→0

Remark 3 Note that Theorem 2 only discusses the situation  → 0 and q > 0, where the limit solution is given by (13) and (14). However, the systems (13) and (14), where Z i are Hill functions dependent on q, are well studied by several authors. For the Z i∗ , analysis and calculations of the limit solutions of (13) and (14) as q → 0, we refer the reader to Plahte and Kjøglum 2005; Gouzé and Sari 2002, and Machina et al. 2013a. Remark 4 Theorem 2 shows that y approaches its pseudo-equilibrium quickly and then for an arbitrarily long finite time interval stays close to the pseudo-equilibrium, which evolves slowly with x. Solutions, xi , to system (11) converge as  → 0 to solutions of system (13), which is of the same form as the original system, (1), uniformly on finite time intervals of arbitrary length, but not necessarily on t ∈ (0, ∞). Thus, asymptotics

123

R. Edwards et al.

may differ when ε > 0. If system (13) as q → 0 has an asymptotically stable fixed point, then the point may not be stable in system (11) for any finite ε (see Remark 6 and Example 4 in Sect. 6, where the system does not explicitly contain  but admits arbitrarily small parameters αi > 0, κi > 0, i = 1, 2). Example 2 Consider the following system, which is an expansion of the system in Example 1. x˙1 = Z 1 + Z 2 − 2Z 1 Z 2 − β1 x1 , x˙2 = 1 − Z 1 Z 2 − β2 x2 , (15) y˙1 = 1ε (κ1 x1 − α1 y1 ), y˙2 = 1ε (κ2 x2 − α2 y2 ).

where Z i = H (yi , θi , q), i = 1, 2, and 0 < ε 0. The system (16) can be studied following the analysis of Plahte and Kjøglum (2005) and in fact can be rescaled to put it into the form of our Example 1 but with different values of parameters. Thus, it has the same white, black and transparent walls, and an SSP in the wall as in Example 1, for appropriate parameter values (Fig. 1). While the original system (4) in Example 1 has an asymptotically stable SSP at the threshold intersection, Theorem 2 does not allow us to conclude that system (15) 4 4

y2∗

3

3

2

y2 2

1 1 0 0

1

2

y1∗

(a)

3

4 0

0

1

2

y1

3

4

(b)

Fig. 1 a The phase space of (16) with a few trajectories plotted. The small parameter is q = 0.01, and the other parameter values are as follows: α1 = 1, α2 = 2, β1 = 0.6, β2 = 0.9, κ1 = 3, κ2 = 4, θ1 = θ2 = 1. b Evolution in the y1 , y2 plane of two trajectories of (15) with different initial values: x0 = (1, 2), y0 = (1, 0.5) for the solid line, and x0 = (1.3, 0.75), y0 = (0.5, 1.5) for the dashed line. The small parameter values are q = 0.01 and ε = 0.02, and other parameter values are as in (a). Initially, the trajectory in the y1 , y2 plane quickly reaches the vicinity of the point (κ1 x1,0 /α1 , κ2 x2.0 /α2 ), which is (3, 4) for the solid line and (3.9, 1.5) for the dashed line. Then it follows the dynamics of (16) while approaching the SSP near (2.25, 1)

123

Transcription–Translation Networks

has a corresponding asymptotically stable stationary point that converges to that of system (4) as ε → 0. We investigate this question further in Sect. 6.  5 Infinitely Fast mRNA Rates Now we consider the situation where the mRNA transcription rates are much faster than the protein translation rates. We show that the limit dynamics of the expanded system with a constant G i is given by (1), with an appropriate choice of Fi and G i for that system. Consider the system of dimension 2n as in (5) but with 1ε scaling the rate of the mRNA dynamics. x˙i = 1ε (Fi (Z ) − βi xi ), y˙i = κi xi − αi yi ,

i = 1, . . . , n,

(17)

where Z = (Z 1 , . . . , Z n ), Z i = H (yi , θi , q). Let q > 0 be fixed. The system (17) represents a singular perturbation problem as ε → 0 with the fast flow xi and the slow flow yi . The pseudo-stationary solution xi∗ = Fi (Z )/βi of the fast flow is obtained from Fi (Z ) − βi xi = 0. It is easy to check that the stationary solution xi∗ is asymptotically stable for each fixed yi . Thus, all assumptions of the Tikhonov theorem are fulfilled, and the slow flow is therefore governed by Fi (Z ) − αi yi . (18) y˙i = κi βi The system (18) is of the form (1) with Fi (Z ) replaced by By the Tikhonov theorem, the following result holds.

κi βi

Fi (Z ).

Theorem 3 Let q > 0 and an initial condition (x0 , y0 ) be fixed. The solutions yiε (t), i = 1, . . . , n, of (17) converge uniformly on any finite time interval as ε → 0 to yi , which is a solution to the system (18). Remark 5 The same remarks apply here as Remarks 3 and 4 following Theorem 2. The limit solutions of system (18) as q → 0 are well studied by Plahte and Kjøglum 2005; Gouzé and Sari 2002; Machina et al. 2013a. Solutions to system (17) converge to solutions of system (18) on finite time intervals of arbitrary length, but the result does not tell us about the asymptotic behaviour for any ε > 0. Example 3 Consider the following system, which is an expansion of the system in Example 1. x˙1 = 1ε (Z 1 + Z 2 − 2Z 1 Z 2 − β1 x1 ), x˙2 = 1ε (1 − Z 1 Z 2 − β2 x2 ), y˙1 = κ1 x1 − α1 y1 , y˙2 = κ2 x2 − α2 y2 .

(19)

where Z i = H (yi , θi , q), i = 1, 2, 0 < ε 0, αi > 0, i = 1, 2, are small, the system (21) represents the situation with different timescales for the dynamics of mRNA and protein concentrations. In particular, we are interested in stability of fixed points of (21) where y1 and y2 are simultaneously switching. Thus, we seek fixed points with the coordinates (x1∗ , x2∗ , Z 1∗ , Z 2∗ ), with Z i∗ ∈ (0, 1). We therefore rewrite (21) as

123

Transcription–Translation Networks

x˙1 = F1 (Z ) − β1 x1 x˙2 = F2 (Z ) − β2x2  Z˙ 1 = Z 1(1−Z 1)q κ1 x1 − α1 θ1 Z˙ 2 =

Z qθ1 1−Z1 1

Z 2(1−Z 2)

qθ2

q Z2 1−Z 2

Z1 1−Z 1

q  (22)

  q  Z2 κ2 x2 − αi θ2 1−Z , 2 ∗q

αi θi Z i ∗ with q > 0. Solving Z˙ i = 0, we obtain xi∗ = κ (1−Z ∗ q , where Z i are roots of the i i) system F1 (Z ) − β1 x1∗ = 0 (23) F2 (Z ) − β2 x2∗ = 0.

In the limit q → 0, notice that xi∗ = ακi iθi , which greatly simplifies the root finding problem in (23). Assuming that such a fixed point (x1∗ , x2∗ , Z 1∗ , Z 2∗ ) exists, we seek conditions on the parameters of the expanded system, κi and αi , that guarantee stability or instability at the fixed point. We analyse the Jacobian at P = (x1∗ , x2∗ , Z 1∗ , Z 2∗ ), ⎡

⎤ −β1 0 ∂∂ Zx˙11 ∂∂ Zx˙12 ⎢ 0 −β ∂ x˙2 ∂ x˙2 ⎥ 2 ∂ Z1 ∂ Z2 ⎥ ⎢ J (P) = ⎢ A1 ⎥ 0 −α1 0 ⎦ ⎣ q 0 Aq2 0 −α2 , where Ai =

κi Z i∗ (1−Z i∗ )

∗ q . Z θi 1−Zi ∗

We are interested in stability when q is very small. If λ

i

denotes an eigenvalue of J (P), we write the characteristic polynomial in terms of √ μ = λ q, multiply through by q 2 and take the limit as q → 0, to obtain

μ4 −



 ∂ x˙2 ∂ x˙2 ∂ x˙1 ∂ x˙1 ∂ x˙2 ∂ x˙1 A2 + A1 μ2 − A1 A2 − = 0. ∂ Z2 ∂ Z1 ∂ Z1 ∂ Z2 ∂ Z2 ∂ Z1

(24)

Since the roots of a polynomial are continuous in their coefficients, we can make the roots of the characteristic polynomial in μ as close to the roots of (24) as we wish by taking q small enough. This implies that whenever there is a root of (24) with a nonzero real part, we can find a q small enough that the characteristic polynomial in μ will also have a root with a nonzero real part. More importantly, since if λ is a root of (24), then −λ is also a root of (24), we can find a q small enough such that the characteristic polynomial in μ will have a root with a positive real part. Finally, since the roots of the characteristic polynomial in λ are positive scalar multiples of roots of the characteristic polynomial in μ, whenever the characteristic polynomial in μ has a root with a positive real part, there is an eigenvalue of J (P) with a positive real part. Therefore, for a given fixed point P of (22), with Z i ∈ (0, 1), we know that P is unstable when q < δ for some δ > 0, if there exists a root of (24) with a nonzero real part. Thus, the only situation in which we do not immediately ascertain instability for small enough q is when all solutions of (24) have pure imaginary roots. In that case, we do not know the stability of the equilibria of (22).

123

R. Edwards et al.

Remark 6 We have shown that any possible SSPs of system (21) located in codimension 2 switching domains are generically unstable, except for possibly very specific cases when stability of the equilibria is unknown. Note that this is true for arbitrary positive values of αi , κi , including cases when the timescales of mRNA and protein are very different. However, the corresponding SSPs of the protein-only model reduced from (21) may in general be stable, as in the Plahte–Kjøglum system (2005) and the example below. Our results agree with the analysis of Polynikis et al. (2009) for protein-only systems without self-regulation. Under the assumption that mRNA xi is not regulated by Z i , they show that for some values of αi > 0, κi > 0 the system (21) exhibits oscillations around an unstable singular stationary point. The amplitude of these oscillations decreases as αi and κi decrease. In the limit αi → 0, κi → 0, the equilibria predicted by their reduced protein-only model are always stable; see Polynikis et al. (2009, Sections 5.2, 5.3). Similar results were found earlier by Tyson and Othmer (1978), also for systems without self-regulation. Example 4 We now apply this result to the mRNA-protein expansion of the Plahte– Kjøglum system, x˙1 = Z 1 + Z 2 − 2Z 1 Z 2 − β1 x1 , x˙2 = 1 − Z 1 Z 2 − β2 x2 , (25) y˙1 = κ1 x1 − α1 y1 , y˙2 = κ2 x2 − α2 y2 , with the same parameters as in the original system: β1 = 0.6, β2 = 0.9 and θi = 1. We ask whether there exist values of κi and αi such that the dynamics of the expanded system differ from those in the original Plahte–Kjøglum system (4). First, we look for conditions on κi and αi ensuring the existence of fixed points of the type (x1∗ , x2∗ , Z 1∗ , Z 2∗ ), with Z i∗ ∈ (0, 1). Following the singular perturbation approach, making the fast time substitution τ = qt with x = ∂∂τx and taking q → 0, the system becomes x1 = 0, x2 = 0, Z 1 = Z 2 =

Z 1 (1−Z 1 ) (κ1 x1 y1 Z 2 (1−Z 2 ) (κ2 x2 y2

− α1 y1 ), − α2 y2 ).

Fixed points with Z i∗ ∈ (0, 1) occur when yi = θi , which yields xi∗ = ακi iθi . Returning to slow time, putting these values for xi∗ into the first two equations above gives x˙1 = Z 1 + Z 2 − 2Z 1 Z 2 − β1 x1∗ , x˙2 = 1 − Z 1 Z 2 − β2 x2∗ , Z˙ 1 = 0, Z˙ 2 = 0.

123

Transcription–Translation Networks

At a stationary point, all these derivatives must be 0. From the second equation above, Z 2∗ =

1 − β2 x2∗ , Z 1∗

(26)

and then we obtain a quadratic in Z 1∗ , F(Z 1∗ ) = Z 1∗ − (c1 + 2(1 − c2 ))Z 1∗ + (1 − c2 ) = 0, 2

(27)

where ci = βi xi∗ . The roots of F, when combined with Z 2∗ and xi∗ , are fixed points of (25). Since we are only interested in roots in the open interval (0, 1), it is clear from (26) that in order to have a stationary point c2 < 1 and thus F(0) = 1 − c2 > 0. In that case, if F(1) = c2 − c1 < 0, then F has exactly one root in (0, 1) by the intermediate value theorem. If F(1) > 0, then the number of roots in (0, 1) depends c2

on Fmin = (c2 − c1 )(1 − c2 ) − 41 : If Fmin > 0 then there are no roots in (0, 1), while if Fmin < 0 there are 2 roots in (0, 1). Thus, the conditions for existence of exactly one stationary point are c2 < 1 and c1 > c2 . The conditions for the existence of two stationary points are c1 < c2 < 1 and c2

(c2 − c1 )(1 − c2 ) − 41 < 0. Figure 3 labels the regions in the c1 , c2 parameter plane where system (25) has 0, 1 or 2 fixed points, (Z 1∗ , Z 2∗ ), in the open square (0, 1)2 . In order to compare the xi variables of the expanded system (25) to those of the original Plahte–Kjøglum system (4), note that if κi = αi , then ci = βi θi . In this case c1 = 0.6 and c2 = 0.9, which satisfy the requirements for the existence of two fixed points and give the same Z i∗ values as in the original Plahte–Kjøglum system (4). With these values for Z i∗ , the roots of the reduced characteristic equation (24) are (0.28321i, −0.28321i, 0.31456, −0.31456) for the first fixed point, and (0.17325 − Fig. 3 (Color Figure Online) Values of c1 and c2 where system (25) will have 0, 1 or 2 fixed points, (Z 1∗ , Z 2∗ ), in the open square (0, 1)2

123

R. Edwards et al.

(a)

(b)

Fig. 4 (Color Figure Online) Numerical simulations of system (25) with κ1 = α1 = 3.33, β1 = 0.6, κ2 = α2 = 5.0, β2 = 0.9, and q = 0.001. a A projection of the 4-dimensional phase space onto the x1 , x2 -plane, showing four solutions, none of which converge to the threshold intersection. b A blow-up of the region around the fixed point near (1.5, 1). In both panes, the asterisk indicates the fixed point in the wall where y1 > 1, x1 > 1, and y2 and x2 are near 1

0.24305i, 0.17325 + 0.24305i, −0.17325 − 0.24305i, −0.17325 + 0.24305i) for the second fixed point (after rounding). In both cases, there is a root with positive real part, so we can find a q > 0 small enough such that both fixed points are unstable. Numerical solutions of (25) agree with this analysis, as shown in Fig. 4, with q = 0.001. This is not consistent with the behaviour of the original Plahte–Kjøglum system (4) in the vicinity of the threshold intersection (i.e., in the fast variables, Z 1 , Z 2 ), where one fixed point is stable. The implication is that there are parameter values κi , αi (even if βi = γi from the original system) for which the SSP in (25) at y1 = θ1 , y2 = θ2 is no longer stable (see Remark 4). 

7 Analysis of the Flow in Regular Domains If timescales are of the same order (not orders of magnitude different), the protein dynamics of (5) cannot be analysed just by studying the simpler system (1). A natural question is whether the methods of analysis developed to determine the flow of system (1) through regular domains can be extended to an analysis for system (5). In this section, we take q → 0 and describe a mapping from one wall to another that allows an analytical construction of the solution of the switched system (5) for any finite number of transitions. (k) (k) Given a Boolean vector [B1 , . . . , Bn ], we start in the regular domain B (k) where (k) (k) (k) (k) Z i = Bi , i = 1, . . . , n, and denote by the constants πi = Fi (B1 , . . . , Bn ), so that initially (5) is given by x˙i = πi(k) − βi xi , y˙i = κi xi − αi yi , (k)

with yi (0) = yi

123

i = 1, . . . , n,

(28) (k)

as the initial value in the box B (k) and xi (0) = xi .

Transcription–Translation Networks

The solution to the first equation in (28) is (k)

π xi (t) = i − βi



 (k) πi (k) − xi e−βi t , βi

i = 1, . . . , n.

(29)

Inserting this into the second equation in (28) gives (k)

π y˙i = κi i + κi βi

xi(k)

(k)

π − i βi

 e−βi t − αi yi ,

which has solution

 (k) (k) κi πi πi κi (k) yi (t) = + xi − e−βi t βi αi αi − βi βi



 κi πi(k) πi(k) κi (k) (k) + yi − − xi − e−αi t , βi αi αi − βi βi

(30)

under the assumption (which we make here and in the remainder of this section) that βi = αi . Thus in a particular box B (k) given by a Boolean vector [B1(k) , . . . , Bn(k) ], the solution to (28) is given by (29) and (30). Since all walls in the model (5) are transparent ( y˙i does not depend on Z i ), the solution can be continuously extended from one box to a neighbouring box using these formulas. For example, suppose that at some time T (k) = Tsk , the solution passes through the (k) (k) wall ysk = θsk at the point ysk (T (k) ) = ysk = θsk , xi (T (k) ) = xi , and yi (T (k) ) = (k) yi , i = 1, . . . , n, and proceeds into the box B (k) . The index (k) is a counter for the intersections of the solution and walls. Consider all j such that y j = θ j are n walls adjacent to the box B (k) . For each j, we have, from (30), (k)

θj =

κjπj

βjαj

+

κj + αj − βj

(k) yj

(k)

−

κjπj

βjαj

(k) xj

(k) 

−

πj

βj

κj − αj − βj

e−β j T j

(k) xj

(k) 

−

πj

βj

e−α j T j .

(31)

 (k) (k) κjπ πj κj (k) x , and With the substitutions u = e−β j T j , a = β j αj j − θ j , b = α j −β − j βj j

 (k) (k) κjπ πj κj (k) (k) c = y j − β j αj j − α j −β x , the equality (31) can be rewritten as − j βj j 0 = a + bu + cu

αj βj

,

(32)

123

R. Edwards et al. Fig. 5 A few trajectories of (33) are plotted in its phase space. All the trajectories spiral into the singular stationary point near (1/3, 1), which is asymptotically stable by Theorem 1. The small parameter value is q = 0.01

4 3

y

2 1 0 0

1

2

x

3

4

α

which is potentially a transcendental equation (if β jj is an irrational number). If for a particular j, there  exist N ≥ 1 positive solutions of (32), say u¯  for  = 1, . . . , N , then T j = min − lnβu¯j  (the first moment of hitting a wall). 

If there exists a positive solution of (32) for some j, then the index sk+1 , specified by Tsk+1 = min{T j }, indicates the wall ysk+1 = θsk+1 that the solution crosses at step j

(k + 1). Otherwise, nothing switches and the solution stays in the same box B (k) as t → ∞ with xi (t) →

(k)

πi βi

(k)

(k)

κi πi βi αi > 0, monotonically. (k) (k) ysk+1 , πsk+1 , βsk+1 , κsk+1 , αsk+1 , θsk+1 )

, and yi (t) →

of (32) is When a solution u¯ = u(x ¯ sk+1 , inserted into (29) and (30) for e−βi t , the remaining coordinates xi and yi , i = 1, . . . , n, of the intersection point (x (k+1) , y (k+1) ) with the wall ysk = θsk are obtained. 2n We thus have described a mapping M : R2n + → R+ from one wall to another, i.e., (x (k+1) , y (k+1) ) = M(x (k) , y (k) ). Example 5 Consider the system x˙ = 1 − Z − βx, y˙ = κ x − αy,

(33)

where Z = H (y, θ, q), α = 1, β = 0.6, κ = 3, θ = 1. These equations would be obtained by considering the original Plahte–Kjøglum example in the vicinity of either of its black walls and expanding to the 4-dimensional mRNA-protein system. Some trajectories, calculated according to the above analysis, are plotted in Fig. 5. In this example, there is only one threshold, dividing the phase space into two boxes, and we calculated repeated transitions across the same threshold at y = 1. The stationary solution of (33) is ( ακ y ∗ , y ∗ ) (see Sect. 3), where y ∗ is the stationary solution of κ (34) y˙ = (1 − Z ) − βy. α The system (34) has a unique asymptotically stable singular stationary solution at y ∗ = θ = 1. Note that in this case, the stability of the SSP in the original protein-only model is preserved in the expanded mRNA-protein model. 

123

Transcription–Translation Networks

Figure 4b shows similar behaviour in a four-dimensional system with different parameters, where x2 and y2 spiral in to a fixed point as above.

8 Discussion We have obtained results regarding the more realistic model framework for gene regulatory networks that keeps track of both mRNA and protein concentrations and thus models both transcription and translation as separate steps. We have also drawn some comparisons with the simpler modelling framework, studied at length in the past, in which the only variables are protein concentrations and they are taken to regulate each other’s production directly. We have shown that in the limit of infinitely fast rates for either the mRNA or the protein, solutions of the expanded system converge to solutions of the proteinonly model, but for less-than-infinite rates this conclusion only holds for a finite time interval. Thus, asymptotic behaviour can differ in the expanded system and stability properties of fixed points can be changed. Our Theorem 1, however, guarantees that stability of a SSP in a co-dimension 1 switching domain of the protein-only model is preserved in the expanded model, at least for steep enough sigmoids. Thus, it is only in higher co-dimension switching domains that stability can be changed in the expanded model. In fact, our analysis of SSPs in co-dimension 2 switching domains shows that they are generally unstable, although there are special situations where the linearized system in the limit q → 0 is neutrally stable (Sect. 6, in the discussion following (24)). In these cases, we have not determined the effect on stability of the SSP of the nonlinear terms nor of the perturbation q > 0. This question remains for future work. It is clear, though, that the stability of a fixed point in a co-dimension 2 domain of the protein-only model can be lost in the expanded model, depending on the values of the parameters. This is true for similar as well as different timescales between the protein and mRNA dynamics. The overall picture is that the switching domains are more benign in the expanded model than in the protein-only model. In the infinitely steep limit (q → 0) case of the protein-only model, solutions can reach the switching threshold of a single variable (a wall) in finite time and remain stuck there for an interval of time (infinite if approaching an SSP within the switching domain). In the infinitely steep limit of the mRNA-protein model, there may still be fixed points in walls, but these are approached only asymptotically after spiralling in the corresponding mRNA and protein variables for an infinite period of time (Sect. 7, Example 5). All walls are transparent, so there is no sliding motion in switching domains, even in the limit q → 0. Co-dimension 2 switching domains are generally unstable, so solutions can normally not be trapped there, though there are special cases where we are not yet sure what happens, as mentioned above. We conjecture, though, that any stable fixed points in higher codimension switching domains are also approached only asymptotically, after an infinite period of spiralling. Because of this transparency of switching domains, and the lack of singular flow, the most mathematically difficult and sensitive situations that arise in the protein-only models are avoided in the mRNA-protein models. This promises to

123

R. Edwards et al.

make analysis easier, though the flow in regular domains (the mapping from entry wall to exit wall, which we derived in Sect. 7) is now somewhat more complicated. For example, the flow may enter and leave such a domain across the same wall. We now discuss briefly the relationship between our work and that of two other groups. One group reduced a delay system to an ordinary differential equation system with additional variables such that if only one additional variable is introduced, our mRNA-protein system is obtained. The other tackled a similar model expansion to that of the current paper, but focussed on different questions, such as the behaviour of feedback loops and an oscillating example. There is some overlap, however, and that is discussed next. 8.1 Relation to the Analysis of Gedeon, Cummins and Heys In general, our conclusions regarding stability of SSPs are consistent with those Gedeon et al. (2012), who deal with a similar expansion to include both mRNA and protein dynamics, but only with smooth sigmoids (never the infinitely steep limit). Much of their paper deals with specific example networks, but their Section 3.2 gives a general result about the stability of fixed points in the protein-only model compared to the mRNA-protein model. We agree, in particular, with their statement that the inclusion of the mRNA variables generally has a destabilizing effect on equilibria of the protein-only system or else preserves their stability. However, they also find a small parameter region in which it is possible for an unstable equilibrium to be stabilized by the inclusion of the mRNA variables (their Theorem 3.3, part 3), which we have not found. The conflict is only apparent, however, as their result does not apply to our model, except possibly in some rare special cases. Their Theorem 3.3 can be satisfied, for example, by expanding the equation x˙ = κ Z (x) − dx,

with

Z (x) =

x θ+x

(35)

(having the Hill function with q = 1) and d, κ, θ > 0. The equilibrium at x = 0 is unstable if κθ > d. The expanded system is then x˙ = cy − dx, y˙ = κ Z (x) − by,

(36)

where x and y here are the concentrations of protein and mRNA, respectively (as in Gedeon et al. 2012). The Jacobian matrix at the origin is  −d c . J (0, 0) = κ θ −b 

κ bd The origin is stable if det(J (0, 0)) > 0, i.e., if κθ < bd c . So if d < θ < c , the origin is stabilized by the expansion. However, this cannot occur with Hill

123

Transcription–Translation Networks

functions with exponent > 1 (i.e., q < 1) in this example since then the equilibrium at x = 0 of (35) is stable as is the equilibrium of the expanded system (36). The condition in their Theorem 3.3, resulting from the assumption that the equilibrium point is in the identical location in the protein-only and expanded systems, means that it does not apply in almost all cases with which we are concerned. Gedeon and co-authors also consider the infinitely fast limits, in a slightly different way, but do not deal with the steep sigmoid limit and focus on a particular example network with periodic behaviour. 8.2 Relation to Delay Systems Our model (5) fits the delay differential model considered by Ponosov (2005) and Shlykova et al. (2008). Indeed, after the variable substitution, αi X i = κi xi , and allowing a variable degradation term for the mRNA, system (5) becomes i (Z ) − G i (Z )X i , X˙ i = F y˙i = αi X i − αi yi ,

i = 1, . . . , n,

(37)

i (Z ) = κi Fi (Z ), Z = (Z 1 , . . . , Z n ) and Z i = H (yi , θi , q). with F αi Ponosov and co-authors considered the delay system i (Z ) − G i (Z )X i , X˙ i = F y˙i = (RX i )(t)

i = 1, . . . , n

(38)

with the delay operator t (RX i )(t) =

G(t − s)X i (s)ds, −∞

where the function G is a linear combination of gamma distributions (Ponosov 2005, Eqs. (4.4)–(4.6)). Our system (37) is a particular case of (38), where G(t − s) = αi e−αi (t−s) , also known in the literature as the weak generic delay kernel (see, e.g., Ponosov 2005). The delay differential models (38) were studied by Ponosov (2005) and Shlykova et al. (2008), who considered (solved) the problem of singular stationary points (see Ponosov (2005), Theorem 8.2 and Conclusion). Theorem 6.6 (localization principle) in Shlykova et al. 2008, Theorem 7.2) reduced the stability problem to singular coordinates only (their analysis is done under the assumption that sigmoids are logoid functions). Acknowledgments The authors thank two anonymous reviewers for helpful comments and for drawing their attention to the paper by Polynikis et al. (2009).

123

R. Edwards et al.

Appendix A system of the form (1) with autoregulation, i.e. when x˙i depends on Z i , may contain so-called switching domains, threshold hyperplanes or intersections of hyperplanes where solutions may remain for a finite or infinite period of time. Several methods have been developed for dealing with the dynamics of the switched system (1) in switching domains (Filippov 1998; Plahte and Kjøglum 2005). Here we propose another method consisting in introducing in a special way a bigger system without autoregulation, which is a singular perturbation to the original system. In the absence of autoregulation, the problem of dealing with switching domains is avoided, and we have a system of the same class as before (i.e., a Glass network), though with twice as many variables. We show that this specifically chosen bigger system has the same dynamics in the singular perturbation limit as the system (1). To this end, we introduce the following 2n-dimensional system with n extra artificial variables yi : i = 1, . . . , n, x˙i = Fi (Z ) − G i (Z )xi , (39) y˙i = 1ε (i − αi yi ), where Z = (Z 1 , . . . , Z n ), Z i = H (yi , θ yi , q), i = H (xi , θxi , q). Both steps now involve sigmoid functions, so we call it the sigmoid-sigmoid system (as opposed to the main body of the text, in which the expanded system was a sigmoid-linear system). This might also be considered as a two-step model for gene regulation in which the first step is expression of a protein (suppressing the transcription and translation steps) and the second is a post-translational modification of the protein by a catalytic (and therefore sigmoidal) reaction. We show that when the rates y˙i are infinitely fast (ε → 0), the original dynamics of xi of the type (1) are recovered for a specific choice of sigmoid function S. Note that since x˙i ( y˙i ) is independent of i (Z i ), all switching hyperplanes xi = θxi (yi = θ yi ) are ‘transparent’ (the solution just passes through them). Let q > 0 be fixed. Assume that ε → 0. The system (39) represents a singular perturbation problem as ε → 0 with the fast flow yi and the slow flow xi . The stationary solution yi∗ = i /αi of the fast flow is obtained from i − αi yi = 0. It is easy to check that the stationary solution yi∗ is asymptotically stable for each fixed xi . Thus, all assumptions of the Tikhonov theorem are fulfilled, and the slow flow is therefore governed by (40) x˙i = Fi (Z 1∗ , . . . , Z n∗ ) − G i (Z 1∗ , . . . , Z n∗ )xi , where Z i∗ = H (yi∗ , θ yi , q) = H (i /αi , θ yi , q) = H (i , θ yi αi , q). If θ yi αi ≥ 1, then i < θ yi αi and H (i , θ yi αi , q) → 0 as q → 0 for any finite xi . Assume now that the parameters are such that θ yi αi < 1. If we think of H (i , θ yi αi , q) as a function of xi , then from the properties of the Hill function (xi , H −1 (θ yi αi , θxi , q), q) = H (xi , H (i , θ yi αi , q) = H θi , q),

123

(41)

Transcription–Translation Networks Fig. 6 The functions H (x1 , θ y1 , q) (solid line) and (x1 , θ1 , q) (dashed line) with H θ y1 = 1, θx1 = 1, q = 0.1, θ1 = α1 = 0.6 and H −1 (θ y1 α1 , θx1 , q) ≈ 1.04

1 0.8

 H H,

0.6 0.4 0.2 0

θˆ1 0

0.5

1

1.5

x1

2

(xi , where H θi , q) is a sigmoid function (not necessarily q a Hill function) of xi with the θ yi αi −1 and the steepness parameter q threshold at θi = H (θ yi αi , θxi , q) = θxi 1−θ yi αi

(see Fig. 6). Note that since θ yi αi < 1, the new threshold θi → θxi as q → 0. In a functional form, (xi , θi , q) = H

1/q 2

1/q

i 1/q

i

=

+ (θ yi αi )1/q

where i = H (xi , θxi , q) =

1/q xi 1/q 1/q xi +θxi

xi 1/q 2

xi

1/q 1/q

+ θ yi αi

1/q

(xi

1/q

+ θxi )1/q

,

(42)

.

The system (40) can be rewritten as Z1, . . . , Z n ) − Gi ( Z1, . . . , Z n )xi , x˙i = Fi (

i = 1, . . . , n,

(43)

(xi , θi , q) with θi = H −1 (θ yi αi , θxi , q) → θxi as q → 0. where Zi = H We thus have proven the following result. Theorem 4 Let q > 0 be fixed. The solutions xi , i = 1, . . . , n, of (39) converge given by (41) as ε → 0 to the solution of the system (43) with the new sigmoid H uniformly on any finite time interval. We now compare the solutions of the system (43) with solutions of a corresponding system (44) below in which the sigmoids are Hill functions. Note that the Z i in (44) are functions of xi , whereas the Z i in (39) are functions of yi . (xi , θi , q) Theorem 5 Both the solution of the system (43) with a new sigmoid Zi = H and the solution of the system x˙i = Fi (Z 1 , . . . , Z n ) − G i (Z 1 , . . . , Z n )xi ,

i = 1, . . . , n,

(44)

with the Hill function Z i = H (xi , θxi , q) converge as q → 0 to the same limit solution uniformly on any finite sequence of finite time intervals required for passing regular

123

R. Edwards et al.

and switching domains, provided that the sequence of domains for (43) and (44) are the same for sufficiently small q. Proof Starting at an initial point, we consider one domain at a time, regular or singular, with the corresponding passage time. For a regular domain, the convergence on the corresponding time interval is obvious. For a singular domain, the convergence is proved below using the Tikhonov theorem and the a priori assumption that the limit solution is in a singular domain. Then convergence on each of these finite time periods can be naturally extended to any finite sequence of these finite periods. Given an ordered set of switching variables S ⊂ {1, . . . , n}, the ordered set of regular variables R = {1, . . . , n} \ S, and a corresponding Boolean vector B R of θs , dimension |R|, consider a motion of (43) confined to the switching domain xs = s ∈ S, Z r = Br , r ∈ R. Separating singular and regular parts of the system (43) gives Z ) − G s ( Z )xs , s ∈ S x˙s = Fs ( Z ) − Gr ( Z )xr , r ∈ R. x˙r = Fr (

(45)

For i ∈ S, since Z i = H (i , θ yi αi , q) with i = H (xi , θxi , q), by the chain rule, d d d H (i , θ yi αi , q) · H (xi , θxi , q). Zi = dxi di dxi

(46)

Differentiating the inverse of the Hill Function, yi = H

−1

(Z i , θ yi , q) = θ yi

Zi 1 − Zi

q ,

0 < q ≤ 1,

(47)

gives

d 1 Z i (1 − Z i ) H (yi , θ yi , q) = · , dyi q yi where yi is given by (47). Using (48), Eq. (46) yields

(48)

d Z i (1 − Z i ) 1 i (1 − i ) Z i )(1 − i ) 1 1 Z i (1 − · · = 2· . Zi = · dxi q i q xi q xi

(49)

Using the derivative of an inverse function in the left-hand side of the first equation in (45) and inserting (49) gives the equivalent system   Z˙ s = Z s (1− Zxss)(1−s ) FS ( Z ) − G s ( Z )xs , s ∈ S, q2 (50) x˙r = Fr ( Z ) − Gr ( Z )xr , r ∈ R. Following the singular perturbation theory, the stretching transformation τ = the first equation in (50) into   Z s = Z s (1− Zxss)(1−s ) Fs ( Z ) − G s ( Z )xs , s ∈ S,

123

t q2

takes

Transcription–Translation Networks

where the prime denotes differentiation with respect to the new ‘fast’ time τ . Letting θs → q → 0 and assuming a priori that the limit solution belongs to the domain xs = Z r = Br , r ∈ R, i.e. that s = Hs (xs , θxs , q) → θxs (see lines before (42)), s ∈ S, Hs (θxs , θxs , q) = 1/2 gives the boundary layer equation  Zs )  Z s (1 − Fs ( Z ) − G s ( Zs = Z )θxs 2θxs

(51)

with Z r = Br . Similarly, see e.g. Plahte and Kjøglum (2005) for details, for the system x˙s = Fs (Z ) − G s (Z )xs , s ∈ S, x˙r = Fr (Z ) − G r (Z )xr , r ∈ R,

(52)

with Z i = H (xi , θxi , q), the boundary layer equation (now in time qt ) is

Zs =

 Z s (1 − Z s )  Fs (Z ) − G s (Z )θxs θxs

(53)

with Z r = Br . Clearly, the ω-limit sets of (51) and (53) are the same. The slow dynamics of (45) and (52), which are governed by the asymptotic behaviour of the fast dynamics, are therefore the same. We have proved that the limit dynamics in both ) are the same, at least on finite time cases (the Hill function H and the new sigmoid H intervals (because of the application of Tikhonov’s Theorem in both cases).  Note that the sequence of domains for (43) and (44) is guaranteed to be the same as long as no two variables become switching variables at the same time in the limit q → 0, which is a set of initial conditions of measure zero. 8.3 Singular Stationary Points in the Sigmoid–Sigmoid System In this section, we describe some general characteristics of fixed points in switching domains for the sigmoid–sigmoid system. We do not specify the value of ε here, as we want to keep this analysis as general as possible. We will observe that the dynamics of the expanded system are consistent with the original system for some values of the additional parameters. We also show that the resulting structure of the expanded system makes stability of fixed points in switching domains much harder to achieve. Consider the expanded system with ε > 0 but not necessarily small. Recall that a trajectory of system (39) is on a wall corresponding to variable x j if x j = θx j and all other variables, xi with i = j and yi ∀i, are off their thresholds. If there is only 1 threshold for each gene, then each variable of an n-dimensional system has 2(n − 1) walls, where each wall is defined by the other variables being above or below their respective thresholds. In systems without autoregulation, all walls are transparent, meaning that as a single variable, say x j , reaches its threshold, x j = θx j , it will pass through without its focal point,

F j (Z ) G j (Z ) ,

changing.

123

R. Edwards et al.

Remark 7 By Assumption 1, the focal point of a trajectory is always in a regular domain when the trajectory itself is in a regular domain. If a trajectory is in a switching domain, it is possible for its focal point to also be in a switching domain. We analyse a 2M-dimensional switching domain of (39). We have M x-variables switching and M y-variables switching. Using singular perturbation theory, we arrive at the system: 





1 (1−1 ) θx1



m (1−m ) m (Z )θxm θxm  Fm (Z ) − G  Z 1 (1−Z 1 ) 1 − α1 θ y1 ) εθ y1

1 = : m =

Z1 = :

Zm =

F1 (Z ) − G 1 (Z )θx1 

Z m (1−Z m ) εθ ym



 (54)

 m − αm θ ym ) .

  ∗ with  ∗ , Z ∗ ∈ (0, 1) Theorem 6 Suppose P = 1∗ , 2∗ , . . . , m∗ , Z 1∗ , Z 2∗ , . . . , Z m j j is a fixed point of (54) and let J (P) denote the Jacobian matrix evaluated at P. If any eigenvalue of J (P) has nonzero real parts, then P is an unstable fixed point of the system, otherwise P is neutrally stable in the linearized system.   0 A ∂ Proof The Jacobian matrix evaluated at P is of the form J (P) = since ∂ ij = B 0 ∂ Z

0 ∀ i, j ∈ 1 . . . m, and ∂ Z ij = 0, ∀ i, j ∈ 1 . . . m. Block matrices of this type have eigenvalues coming by showing that if J (P)V =   in ± pairs.Thiscan be demonstrated   u −u −u λV with V = , then J (P) = −λ . Therefore, for each eigenvalue λ of v v v J (P), −λ is also an eigenvalue. This implies that the fixed point P is unstable if J (P) has at least one eigenvalue with nonzero real part, or J (P) has only pure imaginary eigenvalues, making P neutrally stable in the linearized system. Example 6 In the original Plahte–Kjøglum system (4), with γ1 = 0.6, γ2 = 0.9 and θ1 = θ2 = 1, trajectories either end up at the stable fixed point in the threshold intersection, x1 = θ1 and x2 = θ2 , or at the SSP with x1 > θ1 and x2 = θ2 . The expanded Plahte–Kjøglum system then takes the form x˙1 x˙2 y˙1 y˙2

= Z 1 + Z 2 − 2Z 1 Z 2 − β1 x1 = 1 − Z 1 Z 2 − β2 x2 = 1ε (1 − α1 y1 ) = 1ε (2 − α2 y2 ) ,

(55)

with β1 = 0.6, β2 = 0.9 and θx1 = θx2 = θ y1 = θ y2 = 1. Applying Theorem 6 and an eigenvalue analysis similar to that of Sect. 6, it can be shown that there are parameter values for the expanded system for which the threshold intersection no longer contains stable fixed points, which implies that the dynamics in

123

Transcription–Translation Networks

the expanded system do not match those seen in the original system. This conclusion is independent of the chosen value for ε > 0. The other stable fixed point in the original Plahte–Kjøglum system, (4), is on the wall {x1 > θx1 , x2 = θx2 }. We show that this point is stable in the expanded system (55) as well. Let x1 > θx1 and y1 > θ y1 . Then Z 1 = 1 = 1 and the x2 , y2 equations are independent of x1 and y1 : x˙2 = 1 − Z 2 − β2 x2 , (56) y˙2 = 1ε (2 − α2 y2 ). Applying the singular perturbation theory to (56) in a -vicinity of (θx1 , θ y1 ) gives 2 (1−2 ) (1 − Z 2 − β2 x2 ), x2 Z 2 (1−Z 2 ) ( 2 − α2 y2 ). εy2

˙2 = q q Z˙ 2 =

(57)

In the limit q → 0, x2 → θx2 and y2 → θ y2 and in fast time τ = t/q 2 = Z 2 =

2 (1−2 ) (1 − Z 2 − β2 θx2 ), θx2 Z 2 (1−Z 2 ) (2 − α2 θ y2 ). εθ y2

(58)

The stationary solution 2∗ = α2 θ y2 , Z 2∗ = 1 − β2 θx2 of (58) is neutrally stable. Each of the nested closed curves C k and the corresponding periodic solution ϕ k (·) = (2k , Z 2k ) can be associated with the invariant measure μk on [0, 1]2 by μk (A) =

1 τC k

λ({τ ∈ [0, τC k ] : ϕ k (τ ) ∈ A})

for any measurable subset A ⊂ [0, 1]2 , where λ is the Lebesgue measure on the real line. Thus, there exist infinitely many invariant measures {μk } of the equation (58). By Artstein’s theory Artstein (2002), Machina et al. (2013b) the variable x1 is governed by the differential inclusion  x˙1 ∈ {1 −

Z 2 dμk − β1 x1 | ∀μk }. Ck

The latter equation can be rewritten as

x˙1 ∈ {1 −

1 τC k

τC k Z 2k (t)dτ − β1 x1 | ∀k}.

(59)

0

Solving (58) for Z 2 gives Z 2 = 1 − β2 θx2 −

2 θx2 2 (1−2 ) .

123

R. Edwards et al.

For any closed orbit C k

Z 2av =

1 τC k

τC k Z 2k (t)dτ = 0

1 τC k

τC k (1 − β2 θx2 − 0

= 1 − β2 θx2 −

1 τC k

τC k 0

2 θx2 )dτ 2 (1 − 2 ) θx2 d2 2 (1 − 2 )

= 1 − β2 θx2 . Thus from (59) x˙1 = 1 − Z 2av − β1 x1 , where Z 2av = Z 2∗ = 1 − β2 θx2 , which is the same slow equation as for the original Plahte–Kjøglum system (Plahte and Kjøglum 2005). References Artstein Z (2002) On singularly perturbed ordinary differential equations with measure-valued limits. Math Bohem 127:139–152 Bernstein JA, Khodursky AB, Lin PH, Lin-Chao S, Cohen SN (2002) Global analysis of mRNA decay and abundance in Escherichia coli at single-gene resolution using two-color fluorescent DNA microarrays. PNAS 99:9697–9702 Bionumbers: http://bionumbers.hms.harvard.edu (Accessed 9 Nov. 2014) Casey R, de Jong H, Gouzé J-L (2006) Piecewise-linear models of genetic regulatory networks: equilibria and their stability. J Math Biol 52:27–56 Edwards R (2000) Analysis of continuous-time switching networks. Phys D 146:165–199 Edwards R, Farcot E, Foxall E (2012) Explicit construction of chaotic attractors in Glass networks. Chaos Solitons Fract 45:666–680 Edwards R, Ironi L (2014) Periodic solutions of gene networks with steep sigmoidal regulatory functions. Phys D 282:1–15 Elowitz MB, Leibler S (2000) A synthetic oscillatory network of transcriptional regulators. Nature 403:335– 338 Farcot E (2006) Geometric properties of a class of piecewise affine biological network models. J Math Biol 52:373–418 Farcot E, Gouzé J-L (2010) Limit cycles in piecewise-affine gene network models with multiple interaction loops. Int J Syst Sci 41:119–130 Filippov AF(1998) Differential equations with discontinuous right-hand sides. Nauka, Moscow,(1985) [Russian]; English trans. Kluwer, Dordrecht Garcia-Ojalvon J, Elowitz MB, Strogatz SH (2004) Modeling a synthetic multicellular clock: repressilators coupled by quorum sensing. PNAS 101:10955–10960 Gardner TS, Cantor CR, Collins JJ (2000) Construction of a genetic toggle switch in Escherichia coli. Nature 403:339–342 Gedeon T, Cummins G, Heys JJ (2012) Effect of model selection on prediction of periodic behavior in gene regulatory networks. Bull Math Biol 74:1706–1726 Glass L, Kauffman S (1973) The logical analysis of continuous non-linear biochemical control networks. J Theor Biol 39:103–129 Glass L, Pasternack JS (1978) Stable oscillations in mathematical models of biological control systems. J Math Biol 6:207–223

123

Transcription–Translation Networks Gouzé J-L, Sari T (2002) A class of piecewise linear differential equations arising in biological models. Dyn Syst 17:299–316 Ironi L, Panzeri L, Plahte E, Simoncini V (2011) Dynamics of actively regulated gene networks. Phys D 240:779–794 Killough DB, Edwards R (2005) Bifurcations in Glass networks. Int J Bifurc Chaos 15:395–423 Lewis JE, Glass L (1991) Steady states, limit cycles, and chaos in models of complex biological networks. Int J Bifurc Chaos 1:477–483 Lewis JE, Glass L (1992) Nonlinear dynamics and symbolic dynamics of neural networks. Neural Comput 4:621–642 Machina A, Edwards R, van den Driessche P (2013a) Singular dynamics in gene network models. SIAM J Appl Dyn Syst 12:95–125 Machina A, Edwards R, van den Driessche P (2013b) Sensitive dependence on initial conditions in gene networks. Chaos 23:025101 Machina A, Ponosov A (2011) Filippov solutions in the analysis of piecewise linear models describing gene regulatory networks. Nonlinear Anal 74:882–900 Mosteller RD, Goldstein RV, Nishimoto KR (1980) Metabolism of individual proteins in exponentially growing Escherichia coli. J Biol Chem 255:2524–2532 Paetkau V, Edwards R, Illner R (2006) A model for generating circadian rhythm by coupling ultradian oscillators. Theor Biol Med Model 3:12 Plahte E, Kjøglum S (2005) Analysis and generic properties of gene regulatory networks with graded response functions. Phys D 201:150–176 Polynikis A, Hogan SJ, di Bernardo M (2009) Comparing different ODE modelling approaches for gene regulatory networks. J Theor Biol 261(4):511–530 Ponosov A (2005) Gene regulatory networks and delay differential equations. Electron J Differ Equ 12:117– 141 Shlykova I, Ponosov A, Nepomnyashchikh Y, Shindiapin A (2008) A general framework for stability analysis of gene regulatory networks with delay. Electron J Differ Equ 104:1–36 Thomas R (1973) Boolean formulation of genetic control circuits. J Theor Biol 42:563–585 Tyson JJ, Othmer HG (1978) The dynamics of feedback control circuits in biochemical pathways. Prog Theor Biol 5:1–62 Wilds R, Glass L (2009) An atlas of robust, stable, high-dimensional limit cycles. Int J Bifurc Chaos 19:4055–4096

123