Time-Delayed Models of Gene Regulatory Networks

We discuss different mathematical models of gene regulatory networks as relevant to the onset and development of cancer. After discussion of alternative modelling approaches, we use a paradigmatic two-gene network to focus on the role played by time delays in the dynamics of gene regulatory networks. We contrast the dynamics of the reduced model arising in the limit of fast mRNA dynamics with that of the full model. The review concludes with the discussion of some open problems.


Introduction
Cancer is a complex disease, triggered by multiple mutations in various genes and exacerbated by a number of different behavioural and environmental factors. Some risk factors associated with possible onset and development of cancer are preventable, such as inappropriate diet, physical inactivity, smoking, and drinking [1], while other causes include pathogens (HPV16 and HPV18 are known to cause up to 70% of cervical cancer cases [2]), as well as genetic predisposition. Many studies have focussed on identifying efficient genetic cancer biomarkers, such as specific genes and groups of genes associated with significant number of cases of breast cancer [3] and prostate [4] and pancreatic cancer [5]. Despite this progress, due to significant complexity associated with mutations of various cancer genes, many molecular mechanisms of oncogenesis remain poorly understood.
Recent advances in microarray and high-throughput sequencing technologies have provided pathways for measuring the expression of thousands of genes and mapping most crucial genes and groups of genes controlling different types of cancer. The networks of interactions between DNA, RNA, proteins, and molecules are defined as gene regulatory networks (GRNs). GRNs play a major role in a large number of normal life processes, including cell differentiation, metabolism, the cell cycle, and signal transduction; hence, significant efforts have been made to develop mathematical techniques for their analysis [6][7][8].
GRNs are usually formalised as networks (undirected or directed) where the nodes represent individual genes, proteins, and so forth and the edges correspond to some form of regulation between the nodes. In order to make progress in understanding the onset and development of cancer, as well as to develop effective drug targets, it is essential to be able to reconstruct GRNs pertinent to particular types of cancer from available data. Yeh et al. [9] have used a -nearest-neighbours algorithm to identify GRNs correlated with cancer, tumour grade, and stage in prostate cancer. As an alternative approach, Bonnet et al. [10] have utilised LeMoNe (Learning Module Networks) algorithms to derive GRNs from gene and mRNA expression, as measured in lymphoblastoid cell lines of prostate cancer patients. A rule-based algorithm has been successfully used to determine GRNs in colon cancer [11], and similar kinds of networks have been identified from microarray data using neural fuzzy networks [12]. Madhamshettiwar et al. [13] discuss different approaches to infer GRNs in ovarian cancer, as well as the potential of using these GRNs to develop optimal drug targets. Bayesian network techniques have been employed to construct GRNs from microarray data for breast cancer [14]. In a recent paper, Emmert-Streib et al. [15] have successfully used BC3Net inference algorithm to analyse a large-scale breast cancer gene expression dataset and reconstruct the associated GRN.
In the next section we survey and compare different approaches to model the dynamics of GRNs, making an emphasis on particular biological features that can be best represented by each of the methods. Section 3 focuses on the role of transcriptional and translational time delays in GRN models, and we show how such delays can be introduced in a paradigmatic two-gene activator-inhibitor GRN. Depending on a particular biological regime in which a given GRN is operating, it is often possible to encounter a situation where there is a significant separation of time scales due to, for instance, very fast mRNA dynamics compared to other characteristic time scales. In such a case it is possible to perform dimensional reduction and concentrate on the dynamics of a smaller number of variables. In Section 4 we analyse such reduced model and show how one can derive analytical conditions that lead to a transition from a stable steady state to stable periodic oscillations that are impossible in the model system without the time delays. Section 5 extends the analysis to a full nonlinear system to illustrate differences in stability conditions. The review concludes in Section 6 with a discussion of main results and some open questions.

Mathematical Models of Gene Regulatory Networks
In the analysis of gene regulatory networks and their dynamics, the first step is the identification of key modules or components and possible relations between them, which is often done by interrogating available expression data. Once the topology of the GRN has been fixed, the next step in modelling the dynamics is making realistic assumptions about specific rules that govern the expression of particular genes. Depending on the level of understanding of underlying processes, the complexity of the GRN under investigation, and the specific questions to be addressed, there are several methodologically different approaches that can be employed. Endy and Brent [16] and Hasty et al. [17] discuss biological underpinnings for studying and modelling GRNs, while excellent reviews by de Jong [7], Bernot et al. [8], Tušek and Kurtanjek [18], and Hecker et al. [19] give an overview of mathematical and statistical techniques that have been successfully used to model GRNs, and some of these methods are discussed below.

Boolean Networks.
Some of the first models developed for modelling GRNs were the so-called Boolean networks [20][21][22], where the states of all genes participating in the interactions are represented by binary variables having the values of ON and OFF, or 1 and 0, with the possibility of either synchronous or asynchronous update rules for the nodes. Boolean logic rules are then used to approximate regulatory control of gene expression [23], with updates of binary states of all genes taking place simultaneously [24]. Boolean networks approach has been extended in several directions to provide a better approximation of real GRNs. Shmulevich et al. [25] have proposed a probabilistic analogue of Boolean networks to account for stochastic nature of many processes involved in gene expression. Silvescu and Honavar [26] have proposed temporal Boolean networks, where the next state of genes in the networks is determined not only by their current state but also by a fixed number of their previous states, which effectively allows one to take into account some history of transitions in GRN. Recently, Boolean network models of GRNs have been compared to models based on ordinary differential equations (ODEs), and, in fact, it has been shown that some Boolean models can be rigorously derived as coarse-grained analogues of some ODE models [27]. Significant advantage of using Boolean networks to model GRNs lies in the fact that they allow one to consider networks with a very large number of nodes. At the same time, there are several deficiencies in this approach. The first one concerns the fact that since the gene states only admit the values of ON or OFF, this formalism does not take into account intermediate stages of gene expression [28]. Another issue is that GRNs modelled by Boolean networks can exhibit behaviour not observed in real life; hence, special care has to be taken when choosing the class of admissible Boolean functions [29].

Fuzzy
Methods. Due to intrinsic imprecision and uncertainty associated with gene expression data, it may be appropriate to move away from precise rules of Boolean logic in favour of machine learning techniques based on fuzzy logic. The basic idea is that, rather than trying to reconstruct some assumed fixed gene network topology, one considers the whole family of possible networks with all possible distributions of links between nodes. The problem lies in using actual data to assign appropriate probabilities to each of these configurations, so that for a given input the fuzzy network would provide an output that most resembles actual data. A significant advantage of fuzzy logic for inferring the structure of GRNs lies in their ability to rely on already available knowledge of biological relations between different nodes in the network and, at the same time, being able to recover important previously unknown connections. On the other hand, fuzzy methods for GRN inference are characterised by a high level of computational complexity.
To give a few examples, fuzzy approach has been used to analyse microarray data from the yeast cell cycle and to recover a set of GRNs, with -nearest-neighbour algorithm being used to replace missing data [30]. Woolf and Wang [31] have used -means clustering algorithm to reconstruct and evaluate GRNs for Saccharomyces cerevisiae. In this approach, groups of coregulated genes are considered clusters, and clustering algorithm is then used to detect cluster centres. Volkert and Mahlis [32] have used a smooth response surface algorithm to recover GRNs from gene expression data for Saccharomyces cerevisiae. Approaches based on an artificial bee colony search algorithm have allowed the reconstruction of a GRN in Escherichia coli [33]. A very recent review by Al Qazlan et al. [34] gives an overview of different fuzzy methods, as well as their combinations with other approaches, such as ordinary differential equations, with the purpose of optimising data mining of gene expression and microarray datasets to recover GRNs.

Ordinary and Delay Differential Equation Models.
A very powerful and mathematically insightful methodology for analysis of GRNs is based on nonlinear ordinary or delay differential equations (ODEs or DDEs). In this approach, a gene regulatory network is represented by concentrations of different mRNAs and proteins, and the dynamics can be written as a system of ODEs or DDEs using the law of mass action for individual reactions [6]. Some of the earliest results on ODE models of gene regulation go back to Goodwin [35,36], who introduced and studied a negative feedback loop involving the concentrations of mRNA, an enzyme, and a metabolite. It has been later shown that a negative feedback loop is absolutely essential to ensure the existence of stable periodic solutions, while positive feedback is required for multistationarity [37,38]. This approach was subsequently generalised and expanded [39][40][41][42]; reviews by Smolen et al. [24], de Jong [7], and Hecker et al. [19] discuss some of these models based on systems of nonlinear ODEs. A very important aspect of all these models is a regulation function that controls the rates of gene expression. In light of experimental evidence suggesting monotonic sigmoidal shape of regulation functions [43], a conventional choice for this function is given by the Hill function [44][45][46]. Weiss [47] has discussed various chemical mechanisms associated with the Hill function, including different kinds of ligand binding, and a more recent review of the uses of the Hill function in GRN models can be found in [48].
In order to more accurately represent a switch-like behaviour of the gene expression, several authors have developed models of GRNs using piecewise-linear differential equations, in which the continuous Hill function is replaced by a discontinuous step function [49][50][51][52][53][54]. Besides regular steady states, the piecewise-linear models also allow for singular steady states, which, although important for representing homeostasis in GRNs, are complex to analyse due to discontinuities at the thresholds [55,56]. Polynikis et al. [44] discuss various features of piecewise-linear ODE models and different dynamical regimes that can be exhibited in these models, including possible periodic solutions, sharp-threshold dynamics, and the comparison with models based on continuous regulation function.
In terms of applications to cancer, ODE models have explained aberrant dynamics of the NF-B transcription factor linked to oncogenesis, tumour progression, and resistance to therapy, as well as the dynamics of I B-NF-B [57,58]. Another example is the analysis of the feedback loop between the tumour suppressor p53 and the oncogene Mdm2 [59] and the single-cell response of p53 to radiation-induced DNA damage [60]. There is a clinical evidence suggesting that different components of the PI3K/AKT pathway can lead to aberrant cell growth, metastatic competence, and therapy resistance, and some progress has been made in modelling this pathway and identifying inhibitors responsible for the regulation of PI3K/AKT signalling [61]. Cheng et al. [62] and Edelman et al. [63] give a number of examples of the uses of differential equation based models for the analysis of GRNs in cancer.
Another aspect that has to be properly accounted for in dynamical models is the fact that transcription and translation during gene expression often take place over nonnegligible time periods. Monk [64] has shown how time delays can cause oscillatory gene expression and provide insights into the dynamics of interactions between p53 and Mdm2 proteins associated with cancer suppression. Subsequent research has focused on the role of time delays in GRN dynamics [65][66][67][68][69]. Xiao and Cao [70] have analysed a Hopf bifurcation in a gene network with two transcriptional delays, which occurs when the sum of the delays passes through a critical value, and shown how the amplitude and period of oscillations of gene expression change with the time delays. Due to the fact that it may not be practically possible to identify discrete transcription/translation time delays, a better alternative would be to use models with distributed delay [71]. Models with time delays have been used to understand the regulation of feedback loops involving transcription factors E2F and Myc, known oncogenes, and possible tumour suppressors [72,73]. Ribeiro et al. [74] have developed a delayed stochastic simulation algorithm for analysis of the p53-Mdm2 feedback loop whose malfunction is associated with 50% of cancers. Sequences of multiple reactions with unknown intermediate kinetics can also be successfully analysed using time-delayed models [75,76].

Stochastic Models.
Experimental evidence suggests that significant stochastic fluctuations are observed during gene expression and regulation; hence, in many cases it is paramount to use stochastic models for studying GRN dynamics [77,78]. Even in the absence of extrinsic noise associated with variability in different environmental factors, there are several fundamental processes responsible for intrinsic stochasticity of gene expression [79,80]. One of these is the process of initiation of transcription, which starts by first forming an elongation complex by binding RNA polymerase (RNAp) to the promoter region of the gene, and there is a significant variation in the duration of elongation processes between different transcription events [81][82][83][84]. Binding of RNAp to the promoter regions of different genes results in switching of these genes on and off, thus either blocking or facilitating further transcription, which gives another major source of noise in GRNs. Stochasticity in expression of individual gene results in stochastic behaviour of larger genetic circuits and GRNs [77,85]. Some of the early work on stochastic gene expression emerged from experiments in synthetic biology [86,87] that demonstrated how stochasticity can result in sustained oscillations, and significant amount of research has been subsequently done both theoretically and experimentally on the analysis of stochastic (and delayed) oscillations in gene regulatory networks [80,[88][89][90]. Zavala and Marquez-Lago have recently considered delay-induced oscillations in deterministic and stochastic models of single-cell gene expression, highlighting important differences between these two types of models and associated behaviours [91].
Besides being an intrinsic feature of biological dynamics, stochasticity has proved to be important in the context of engineered genetic switches [86,92]. de Jong [7] and El Samad et al. [93] discuss various methods for modelling stochastic GRN models, including stochastic master equation and various stochastic simulation algorithms. Bratsun et al. [90] have developed an algorithm for analysis of non-Markovian dynamics in GRNs with time delays and showed that these delays are able to induce oscillatory dynamics in the case where deterministic models do not exhibit oscillations. This methodology was later improved, and several exact stochastic simulation algorithms have been developed for simulations of time-delayed models [94,95]. A review by Ribeiro [84] discusses various techniques for simulating stochastic time-delayed dynamics of gene expression, and very recently Jansen and Pfaffelhuber [96] have reviewed the role of delay distribution in the stochastic dynamics during gene expression.
Another way to approach stochasticity in the analysis and reconstruction of GRNs is by using so-called Bayesian networks [97], where gene expression values are represented as random variables and relations between them are probabilistic. Learning techniques for Bayesian networks [98,99] allow one to combine expression data with a priori knowledge to deduce the structure of GRN that best matches the available expression data. Friedman et al. [97] have developed an algorithm for deriving Bayesian networks that circumvents a dimensionality problem, and this method has been used to analyse the cell cycle data for S. cerevisiae containing numerous measurements of mRNA expression levels [100]. Out of 800 genes it was possible to identify a few genes controlling the regulation of cell cycle processes.
The rest of this paper is devoted to consideration of the effects of transcriptional and translational time delays on the dynamics of GRNs. In the next section we introduce the timedelayed model of a two-gene activation-inhibition network together with its quasi-steady state simplification and establish the well-posedness of both models. Section 4 contains the derivation of analytical conditions for stability and Hopf bifurcation in the case of very fast mRNA dynamics, while in Section 5 the analysis is extended to the full time-delayed system. The paper concludes in Section 6 with discussion of results and future research directions.

Time-Delayed Models: Derivation and Positivity
To motivate the analysis of time-delayed effects on gene regulatory dynamics, following Polynikis et al. [44], we consider an activation-inhibition two-gene GRN consisting of two genes and , which are assumed to have no effect on their own expression; at the same time, protein is assumed to activate the expression of gene , while protein inhibits the expression of gene . This is one of the fundamental motifs, which has been shown to be functionally relevant in GRNs [62,101]. Denoting the concentrations of proteins and as and and concentrations of transcribed mRNAs as and , the following system of equations can be derived for the dynamics of this GRN [44]: where are the maximum transcription rates, are the translation rates, are the mRNA degradation rates, and are the protein degradation rates for = , . Equations (1) are called the complete nonlinear model (CNM). To make further analytical progress, the activation and inhibition where and are known as activation and inhibition coefficients and the integer parameters and , known as Hill coefficients, determine the steepness of Hill curves [6]. The parameters and give the values of protein concentrations and , at which the corresponding Hill function achieves half of its maximum value. Depending on the values of transcription rates, this would then lead to a significant increase in the respective mRNAs regulated by these proteins [8,44].
Due to the fact that the dynamics of mRNA is normally much faster than that of related proteins, one can use a quasisteady state assumption to simplify CNM (1) by reducing the number of equations. Effectively, this means assuming that mRNAs have already reached their steady state concentrations, that is, takinġ≈ 0, = , in CNM (1), and then focusing on the dynamics of proteins only, as given by the following simplified nonlinear model (SNM): where = , = . (4) Polynikis et al. [44] have shown that while the CNM exhibits Hopf bifurcation of a positive equilibrium, leading to persistent oscillations, in the case of the SNM model this behaviour can disappear. They have also demonstrated an important role played by the Hill coefficients, as well as the separation of timescales between mRNA and proteins, with a larger scale separation favouring a stable equilibrium rather than oscillatory behaviour. While the transcription and translation may be faster than characteristic times associated with significant changes in protein concentrations (of the order of 5 minutes for transcription + translation and 1 hour for a 50% change in the concentration of translated protein for E. coli [6]), these are, in fact, multistep processes consisting of thousands of consecutive chemical reactions. Hence, the duration of transcription and translation is nonnegligible when considered in the context of GRN dynamics [84,96] and has to be correctly accounted for in mathematical models. To analyse the effects where and are the delays during transcription of mRNAs and and and are the delays during translation of proteins and , respectively. This model will be referred to as the delayed complete nonlinear model (DCNM). Similar to the case of instantaneous transcription and translation, the quasi-steady state assumption simplifies system (5) with parameters and defined in (4). Before proceeding with the analysis, one has to augment models (5) and (6) with the appropriate initial conditions and establish that these models are well posed; that is, their solutions remain nonnegative for all time to ensure their biological feasibility. The initial conditions for DCNM model (5) are given by where max = max ( , , , ) and ( ) ∈ ([− max , 0], R) with ( ) ≥ 0 (− max ≤ ≤ 0, = 1, . . . , 4) and similarly for DSNM model (6). Here, It is further assumed that (0) > 0 and (0) > 0 to ensure that at least some amount of proteins will be produced. We now prove that solution ( ( ), ( ), ( ), ( )) of DCNM model (5) with the initial condition (7) is positive for all > 0. This result can be proven by contradiction, following the methodology used in [102]. As a first step, let us show that ( ) ≥ 0 for all > 0. Let 1 > 0 be the first time when ( 1 ) ( 1 ) = 0; assuming that ( 1 ) = 0 implies ( ) ≥ 0 for all ∈ [0; 1 ] and since 1 is the first time when ( 1 ) = 0, this also means ( 1 )/ ≤ 0; that is, the function ( ) is decreasing at = 1 . On the other hand, evaluating the second equation of system (5) at = 1 yields which gives a contradiction. Since (0) > 0, this implies ( ) > 0 for all > 0. Now that the positivity of ( ) has been established, let 2 > 0 be the first time when ( 2 ) = 0. In order for this to happen, one must have ( 2 )/ ≤ 0; that is, the function ( ) should be decreasing at = 2 . At the same time, evaluating the last equation of system (5) at = 2 yields which gives a contradiction and, therefore, ( ) > 0 for all > 0. In a similar manner, the positivity of ( ) implies the positivity of ( ), which in turn implies the positivity of ( ).
where satisfies the polynomial equation: and we used the notation = , = .
Even for realistically small values of Hill coefficients, such as = 2, 3 [103] or = 4-8 [104], (12) is too complicated to allow one to analytically find closed form expressions for and other state variables. Despite not having explicit formulae for possible steady states ( , , , ), one can still perform the analysis of stability in terms of system parameters, and such results would be valid for the values of steady state variables that can be accurately and efficiently determined through numerical solution of the polynomial equation (12).

Analysis of the Delayed Simplified Nonlinear Model (DSNM)
In order to gain some first insights into the role of transcriptional and translational delays on the dynamics of GRN, we focus on the behaviour of the delayed simplified nonlinear model (DSNM) (6). To reduce the number of free parameters in the model, we introduce the new variables: which transform the first equation of system (6) intȯ The second equation of system (6) evaluated at − − has the forṁ is the new combined time delay. The equation for characteristic eigenvalues of the linearisation near a steady state ( , ) of system (18) has the form where In limit = 0, this equation reduces to the quadratic equation [44]: whose roots always have negative real parts, since > 0, > 0, and DSNM > 0. This implies that, for = 0, the steady state ( , ) is stable for any values of parameters. To investigate whether this steady state can lose stability for > 0, one can note that = 0 is not a solution of the characteristic equation (20). Hence, the only possible way that the steady state ( , ) can lose its stability is when a pair of complex conjugate eigenvalues crosses the imaginary axis. In the light of this observation, one can look for eigenvalues of (20) in form = for some real > 0. Substituting this into (20) and separating into real and imaginary parts gives Squaring and adding these two equations yields the following equation for = 2 : which can be solved to give the critical frequency as One should note that 2 0 will only admit real values, provided < DSNM , which implies that, for ≥ DSNM , the steady state ( , ) is stable for all values of the time delay . Note that ℎ ( ) = 2 + 2 + 2 > 0 for any ≥ 0.
The critical value of the time delay can be found from (23), which gives where 0 is determined by (25) and arctan corresponds to the principal value of arctan. When = 0, , the characteristic equation ( Substituting the expressions for cos( 0 0, ) and sin( 0 0, ) from system (23) gives Hence, the eigenvalues of the characteristic equation cross the imaginary axis at = 0 (here, 0 = 0,0 ) and never cross back for higher values of . Thus, we have proved the following result.  In Figure 2 we show how the stability boundary varies depending on the parameters and and the time delay . One observes that, for sufficiently high values of , the range of possible values of for which the steady state is unstable is significantly reduced, thus making the system more prone to support a stable positive equilibrium rather than exhibit oscillations. At the Hopf bifurcation, the associated critical value of the time delay monotonically increases with the parameter . At the same time, there is a minimum value of the time delay , such that for smaller than this value the steady state ( , ) is stable for any value of .
In a similar way, the effects of the transcription rates and are illustrated in Figure 3, which shows that the critical transcription rate of the inhibitor increases with decreasing , and, similar to Figure 2, below certain value of , the steady state ( , ) is stable for any value of . Qualitatively similar dependence is observed between the critical value of and the transcription rate , though this dependence is not completely monotonic. Figure 4 demonstrates how increasing the overall time delay results in a Hopf bifurcation of the steady state ( , ) and the emergence of a stable periodic orbit. The shift between individual time series for and can be interpreted in the same way as in predator-prey or activatorinhibitor systems [105]. In accordance with Theorem 1, once the stability of the steady state ( , ) is lost, it can never be regained for higher values of , so the system will be exhibiting oscillatory behaviour. This result highlights the significance of correct mathematical representation of the transcription and translation processes, since inclusion of transcriptional and translational delays can lead to sustained periodic oscillations even in the simplified model, where such oscillations were impossible when the time delays were neglected.
The characteristic equation (31) can be recast in the form where = + + + , At = 0, (33) reduces to a quartic By the Routh-Hurwitz criterion [105], the necessary and sufficient conditions for all roots of (35) to have negative real parts are given by From the fact that all system parameters are positive and using the definitions of , , and in (34), it follows that Δ 1 > 0 and Δ 2 > 0 for any values of the parameters. Since + DCNM > 0, it is sufficient to require Δ 4 > 0 to ensure that condition Δ 3 > 0 is also satisfied. This leads to the following result.
Lemma 2. Let = 0. The steady state ( , , , ) of system (5) is stable whenever the condition From now on, we will assume that the condition in Lemma 2 holds and analyse whether stability can be lost as increases. Since both and DCNM are positive, this means that = 0 is not a root of the characteristic equation (33), so once again the stability can only be lost through a possible Hopf bifurcation. To investigate this possibility, we look for solutions of (33) in the form = for some real > 0. Substituting this into (33) and separating into the real and imaginary parts gives Squaring and adding these equations yields a quartic equation as follows: where = 2 and = 2 − 2 , Without loss of generality, suppose that (38) has four positive real roots, denoted by 1 , 2 , 3 , 4 , respectively, which would give four possible values of : Define and then 0 is the first value of > 0 such that the characteristic equation (33) has a pair of purely imaginary roots. We have the following result. (38). Then the steady state ( , , , ) of system (5) is stable for 0 ≤ < 0 and unstable for > 0 and undergoes a Hopf bifurcation at = 0 .

Theorem 3. Suppose the conditions of Lemma 2 hold and
Proof. The conclusion of Lemma 2 ensures that the steady state ( , , , ) of system (5) is stable at = 0, and the fact that the roots of the characteristic equation (33) depend continuously on implies that the steady state ( , , , ) is also stable for sufficiently small positive values of . Since 0 is the first positive , for which the characteristic eigenvalues lie on the imaginary axis, in order to verify whether or not the steady state actually loses stability at = 0 , one has to compute the sign of Re( )/ | = 0 . Let ( ) = ( ) + ( ) be the root of the characteristic equation (33) Using the expressions for cos( 0 0 ) and sin( 0 0 ) from (37) gives (45) which means that at = 0 a pair of complex conjugate eigenvalues of the characteristic equation (33) crosses the imaginary axis with a positive speed. This implies that the steady state ( , , , ) of system (5) does lose its stability at = 0 . Figure 5 shows the stability boundary of the steady state ( , , , ) of system (5) depending on the transcription rates and and the total time delay . In a manner similar to that for the simplified model, the critical value of the transcription rate at the Hopf bifurcation reduces with increasing . However, a major difference from the DSNM model, as shown in Figure 3, is that now the Hopf bifurcation can take place even at = 0, as the DCNM system is able to support sustained oscillations [44]. In Figure 6 we illustrate the transition from a stable steady state ( , , , ) to a stable periodic solution around this steady state as the time delay passes through the critical value of = 0 .

Conclusions
In this review we have discussed various mathematical models for the analysis of GRNs and focussed on the role played by the transcriptional and translational time delays in the dynamics of a two-gene activator-inhibitor GRN. By reducing the model to the one with a single time delay, we have considered possible behaviour in the quasi-steady state approximation of very fast mRNA dynamics, which has resulted in a lower-dimensional system of DDEs. Due to the presence of time delays, even this simplified model is able to exhibit loss of stability of the positive equilibrium through a Hopf bifurcation and a subsequent emergence of sustained periodic oscillations, which was not possible in the absence of the time delays, as discussed in Polynikis et al. [44]. We have found analytically the boundary of the Hopf bifurcation depending on the total time delay and other system parameters and illustrated different types of behaviour by direct numerical simulations. Our results suggest that once the positive steady state loses its stability, it can never regain it for higher values of the time delay.
We have also studied the stability of a positive steady state in the full system and showed that this steady state can also undergo a Hopf bifurcation depending on the time delay and system parameters. Our analysis extends an earlier result of Polynikis et al. [44] by showing how the critical values of the parameters at the Hopf boundary change when the time delay increases from zero. Numerical simulations have illustrated the transition from a stable positive steady state to a stable periodic solution as the time delay exceeds its critical value.
The work presented in the paper can be extended in several interesting and important research directions. One possibility would be to account for the fact that in most experiments the transcriptional and translational time delay are not fixed but rather obey some form of a delay distribution. Recent work on the effects on delay distribution on system dynamics [106][107][108] has shown that, even for the same mean delay, details of the distribution can also play an important role. He and Cao [71] have used Lyapunov functional approach to derive conditions for global stability of equilibria in some types of GRNs with distributed delays, and it would be insightful to investigate the possibility of extending this methodology to other types of GRNs and various types of delay kernels. Alternatively, one could use the framework of a master stability function for systems with distributed delays [109] to study possible synchronization dynamics in GRNs with a large number of proteins involved.
As it has already been mentioned, in some cases gene expression behaviour is characterised by a switch-like behaviour that can be better modelled using piecewise-linear rather than continuous transcription functions [50,53]. Whilst some preliminary work has been done recently on the analysis of piecewise-linear systems with discrete time delays, primarily in engineering applications [110][111][112], the dynamics of GRNs with piecewise-linear transcription functions and transcriptional/translational delays have remained completely unexplored. Further inclusion of distributed delays would make such models mathematically very challenging, but it could provide a new level of understanding of GRN dynamics.
Besides providing insights into the dynamics of GRNs, there are several practical ways in which models similar to the one described in this review are helpful in monitoring and treatment of cancer. GRN models based on differential equations coupled with other techniques, such as machine learning and Bayesian networks, have proved effective in identifying specific oncogenes that can be used as biomarkers or drug targets [14,62,[113][114][115][116]. Similar kinds of models are useful for modelling cancer cell growth and understanding interactions between tumour growth and immune response and for analysis of the effects of chemotherapy (or immunotherapy) and  drug resistance [62,63,117,118]. The methodology described in this review can be directly used to improve the performance of these models by elucidating the role of transcriptional and translational time delays in GRN dynamics and its impact on various aspects of cancer onset and development.