Computational and Mathematical Methods in Medicine
Volume 2015 (2015), Article ID 347273, 16 pages
http://dx.doi.org/10.1155/2015/347273

Review Article

Time-Delayed Models of Gene Regulatory Networks

K. Parmar,¹ K. B. Blyuss,¹ Y. N. Kyrychko,¹ and S. J. Hogan²

¹Department of Mathematics, University of Sussex, Falmer, Brighton BN1 9QH, UK
²Department of Engineering Mathematics, University of Bristol, Bristol BS8 1TR, UK

Received 14 May 2015; Revised 31 August 2015; Accepted 14 September 2015

Academic Editor: Moisés Santillán

Copyright © 2015 K. Parmar et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

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.

1. 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–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.

2. 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.

2.1. Boolean Networks

Some of the first models developed for modelling GRNs were the so-called Boolean networks [20–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].

2.2. 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.

2.3. 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–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–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–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 IB-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–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].

2.4. 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–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–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 time-delayed 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.

3. 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 functions in system (1) can be written as the following Hill functions:

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 quasi-steady 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, taking , in CNM (1), and then focusing on the dynamics of proteins only, as given by the following simplified nonlinear model (SNM):

where

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 of transcriptional and translational time delays we introduce the following model:

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) to the following delayed simplified nonlinear model (DSNM):

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 and with , and similarly for DSNM model (6). Here, is the Banach space of continuous mappings of interval onto . It is further assumed that and 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 . This result can be proven by contradiction, following the methodology used in [102]. As a first step, let us show that for all . Let be the first time when ; assuming that implies for all and since is the first time when , this also means ; that is, the function is decreasing at . On the other hand, evaluating the second equation of system (5) at yields

which gives a contradiction. Since , this implies for all . Now that the positivity of has been established, let be the first time when . In order for this to happen, one must have ; that is, the function should be decreasing at . At the same time, evaluating the last equation of system (5) at yields

which gives a contradiction and, therefore, for all . In a similar manner, the positivity of implies the positivity of , which in turn implies the positivity of . Hence, all solutions , , , and of DCNM model (5) are positive for all . The same approach can be employed to show positivity of solutions of DSNM model (6).

Steady states of the DCNM model can be found as roots of the following system of algebraic equations:

This gives

where satisfies the polynomial equation:

and we used the notation

Even for realistically small values of Hill coefficients, such as [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).

4. 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) into

The second equation of system (6) evaluated at has the form

and in terms of new variables (14) this can be rewritten as

Thus, system (6) takes form

where

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 , this equation reduces to the quadratic equation [44]:

whose roots always have negative real parts, since , , and . This implies that, for , the steady state is stable for any values of parameters. To investigate whether this steady state can lose stability for , one can note that 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 . Substituting this into (20) and separating into real and imaginary parts gives

Squaring and adding these two equations yields the following equation for :

which can be solved to give the critical frequency as

One should note that will only admit real values, provided , which implies that, for , the steady state is stable for all values of the time delay . Note that

The critical value of the time delay can be found from (23), which gives

where is determined by (25) and arctan corresponds to the principal value of arctan. When , the characteristic equation (20) has a pair of purely imaginary roots. To determine whether or not these roots do indeed cross the imaginary axis, we consider as a root of (20) near , satisfying , , and . Substituting into (20) and differentiating with respect to yields

From this equation, one can find

Substituting the expressions for and from system (23) gives

Hence, the eigenvalues of the characteristic equation cross the imaginary axis at (here, ) and never cross back for higher values of . Thus, we have proved the following result.

Theorem 1. If , the steady state of DSNM system (18) is stable for all values of the time delay . If , this steady state is stable for and unstable for and undergoes Hopf bifurcation at .

Figure 1 illustrates the stability boundary of the steady state of the DSNM system (18) depending on the time delay and the protein degradation rates and , with the parameter values taken from Polynikis et al. [44]. This Figure suggests that, for any fixed value of one of such rates, there is only a limited range of positive values of the other degradation rate, for which, at a given time delay , the positive equilibrium is unstable. For sufficiently high values of and , this steady state is stable regardless of the value of the time delay , confirming the result proved in Theorem 1.

(a)

(b)

Figure 1: Stability boundary of the steady state of DSNM system (18). The steady state is stable below the surface in (a) and to the left of the boundary curves shown in (b). Parameter values are , , , and .

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 .

(a)

(b)

(c)

(d)

Figure 2: Stability boundary of the steady state of DSNM system (18). The steady state is stable below the surface in (a) and (c) and to the left of the boundary curves shown in (b) and (d). Parameter values are , , and .

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.

(a)

(b)

Figure 3: Stability boundary of the steady state of DSNM system (18). The steady state is stable below the surface in (a), and to the left of the boundary curves shown in (b). Parameter values are , , and .

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 activator-inhibitor 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.

(a)

(b)

Figure 4: Numerical solution of DSNM system (18): (a) ; (b) . Parameter values are , , , and . The critical time delay is .

5. Analysis of the Delayed Complete Nonlinear Model (DCNM)

Linearisation of the full nonlinear DCNM model (5) near the steady state results in the following characteristic equation:

where

It immediately follows from the form of the characteristic equation (31) that stability of the steady state is determined not by individual transcriptional and translational delays but rather by their overall combined duration. In the case , the characteristic equation of the DCNM model reduces to the one analysed in Polynikis et al. [44].

The characteristic equation (31) can be recast in the form

where

At , (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 and for any values of the parameters. Since , it is sufficient to require to ensure that condition is also satisfied. This leads to the following result.

Lemma 2. Let . The steady state of system (5) is stable whenever the condition holds.

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 are positive, this means that 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 . 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 and

Without loss of generality, suppose that (38) has four positive real roots, denoted by , respectively, which would give four possible values of :

Dividing the two equations in (37) gives

Define

and then is the first value of such that the characteristic equation (33) has a pair of purely imaginary roots. We have the following result.

Theorem 3. Suppose the conditions of Lemma 2 hold and , where is defined in (38). Then the steady state of system (5) is stable for and unstable for and undergoes a Hopf bifurcation at .

Proof. The conclusion of Lemma 2 ensures that the steady state of system (5) is stable at , 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 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 , one has to compute the sign of . Let be the root of the characteristic equation (33) near , satisfying and . Substituting into (33) and differentiating both sides with respect to gives

This implies, with ,

Using the expressions for and from (37) gives

which means that at 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 .

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 , 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 .

(a)

(b)

(c)

(d)

Figure 5: Stability boundary of the steady state of DCNM system (5). The steady state is stable below the surface in (a), (c) and below the boundary curves shown in (b), (d). Parameter values: , , and .

(a)

(b)

Figure 6: Numerical solution of DCNM system (5): (a) ; (b) . Parameter values: , , , , and . The critical time delay is .

6. 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–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–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–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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

G. Danaei, S. Vander Hoorn, A. D. Lopez, C. J. L. Murray, and M. Ezzati, “Causes of cancer in the world: comparative risk assessment of nine behavioural and environmental risk factors,” The Lancet, vol. 366, no. 9499, pp. 1784–1793, 2005.

M. Boshart, L. Gissmann, H. Ikenberg, A. Kleinheinz, W. Scheurlen, and H. zur Hausen, “A new type of papillomavirus DNA, its presence in genital cancer biopsies and in cell lines derived from cervical cancer,” The EMBO Journal, vol. 3, no. 5, pp. 1151–1157, 1984.

D. F. Easton, “How many more breast cancer predisposition genes are there?” Breast Cancer Research, vol. 1, no. 1, pp. 14–17, 1999.

G. Thomas, K. B. Jacobs, M. Yeager et al., “Multiple loci identified in a genome-wide association study of prostate cancer,” Nature Genetics, vol. 40, no. 3, pp. 310–315, 2008.

B. M. Wolpin, C. Rizzato, P. Kraft et al., “Genome-wide association study identifies multiple susceptibility loci for pancreatic cancer,” Nature Genetics, vol. 46, no. 9, pp. 994–1000, 2014.

U. Alon, An Introduction to Systems Biology, Champan & Hall, Boca Raton, Fla, USA, 2007.

H. de Jong, “Modeling and simulation of genetic regulatory systems: a literature review,” Journal of Computational Biology, vol. 9, no. 1, pp. 67–103, 2002.

G. Bernot, J.-P. Comet, A. Richard, M. Chaves, J.-L. Gouzé, and R. Dayan, “Modeling and analysis of gene regulatory networks,” in Modeling in Computational Biology and Biomedicine, F. Cazals and P. Kornprobst, Eds., pp. 47–80, Springer, Berlin, Germany, 2013.

H.-Y. Yeh, S.-W. Cheng, Y.-C. Lin, C.-Y. Yeh, S.-F. Lin, and V.-W. Soo, “Identifying significant genetic regulatory networks in the prostate cancer from microarray data based on transcription factor analysis and conditional independency,” BMC Medical Genomics, vol. 2, article 70, 2009.

E. Bonnet, T. Michoel, and Y. van de Peer, “Prediction of a gene regulatory network linked to prostate cancer from gene expression, microRNA and clinical data,” Bioinformatics, vol. 26, no. 18, Article ID btq395, pp. i638–i644, 2010.

X. Wang and O. Gotoh, “Inference of cancer-specific gene regulatory networks using soft computing rules,” Gene Regulation and Systems Biology, vol. 2010, no. 4, pp. 19–34, 2010.

S. Vineetha, C. Chandra Shekara Bhat, and S. M. Idicula, “Gene regulatory network from microarray data of colon cancer patients using TSK-type recurrent neural fuzzy network,” Gene, vol. 506, no. 2, pp. 408–416, 2012.

P. B. Madhamshettiwar, S. R. Maetschke, M. J. Davis, A. Reverter, and M. A. Ragan, “Gene regulatory network inference: evaluation and application to ovarian cancer allows the prioritization of drug targets,” Genome Medicine, vol. 4, no. 5, article 41, 2012.

F. K. Ahmad, S. Deris, and N. H. Othman, “The inference of breast cancer metastasis through gene regulatory networks,” Journal of Biomedical Informatics, vol. 45, no. 2, pp. 350–362, 2012.

F. Emmert-Streib, R. D. M. Simoes, P. Mullan, B. Haibe-Kains, and M. Dehmer, “The gene regulatory network for breast cancer: integrated regulatory landscape of cancer hallmarks,” Frontiers in Genetics, vol. 5, article 15, 2014.

D. Endy and R. Brent, “Modelling cellular behaviour,” Nature, vol. 409, no. 6818, pp. 391–395, 2001.

J. Hasty, D. McMillen, F. Isaacs, and J. J. Collins, “Computational studies of gene regulatory networks: in numero molecular biology,” Nature Reviews Genetics, vol. 2, no. 4, pp. 268–279, 2001.

A. Tušek and Ž. Kurtanjek, “Mathematical modelling of gene regulatory networks,” in Applied Biological Engineering—Principles and Practice, chapter 5, InTech, Rijeka, Croatia, 2012.

M. Hecker, S. Lambeck, S. Toepfer, E. van Someren, and R. Guthke, “Gene regulatory network inference: data integration in dynamic models—a review,” BioSystems, vol. 96, no. 1, pp. 86–103, 2009.

S. A. Kauffman, “Metabolic stability and epigenesis in randomly constructed genetic nets,” Journal of Theoretical Biology, vol. 22, no. 3, pp. 437–467, 1969.

R. Thomas, “Boolean formalization of genetic control circuits,” Journal of Theoretical Biology, vol. 42, no. 3, pp. 563–585, 1973.

R. Somogyi and C. A. Sniegoski, “Modeling the complexity of genetic networks: understanding multigenic and pleiotropic regulation,” Complexity, vol. 1, no. 6, pp. 45–63, 1996.

R. Albert, “Boolean modelling of genetic regulatory networks,” in Complex Networks, vol. 650 of Lecture Notes in Physics, pp. 459–481, Springer, Berlin, Germany, 2004.

P. Smolen, D. A. Baxter, and J. H. Byrne, “Modeling transcriptional control in gene networks—methods, recent results, and future directions,” Bulletin of Mathematical Biology, vol. 62, no. 2, pp. 247–292, 2000.

I. Shmulevich, E. R. Dougherty, and W. Zhang, “From Boolean to probabilistic Boolean networks as models of genetic regulatory networks,” Proceedings of the IEEE, vol. 90, no. 11, pp. 1778–1792, 2002.

A. Silvescu and V. Honavar, “Temporal Boolean network models of genetic networks and their inference from gene expression time series,” Complex Systems, vol. 13, no. 1, pp. 61–78, 2001.

M. Davidich and S. Bornholdt, “The transition from differential equations to Boolean networks: a case study in simplifying a regulatory network model,” Journal of Theoretical Biology, vol. 255, no. 3, pp. 269–277, 2008.

R. Xu, G. K. Venayagamoorthy, and D. C. Wunsch II, “Modeling of gene regulatory networks with hybrid differential evolution and particle swarm optimization,” Neural Networks, vol. 20, no. 8, pp. 917–927, 2007.

F. Karlsson and M. Hörnquist, “Order or chaos in Boolean gene networks depends on the mean fraction of canalizing functions,” Physica A, vol. 384, no. 2, pp. 747–757, 2007.

G. N. Brock, V. Pihur, and L. Kubatko, “Detecting gene regulatory networks from microarray data using fuzzy logic,” in Fuzzy Systems in Bioinformatics and Computational Biology, Y. Jin and L. Wang, Eds., vol. 242 of Studies in Fuzziness and Soft Computing, pp. 141–163, Springer, Berlin, Germany, 2009.

P. J. Woolf and Y. Wang, “A fuzzy logic approach to analyzing gene expression data,” Physiological Genomics, vol. 3, no. 1, pp. 9–15, 2000.

L. G. Volkert and N. Malhis, “An efficient method for fuzzy identification of regulatory events in gene expression time series data,” in Proceedings of the IEEE Symposium on Computational Intelligence in Bioinformatics and Computational Biology (CIBCB '04), pp. 17–24, IEEE, La Jolla, Claif, USA, October 2004.

P. Das, P. Rakshit, A. Konar, M. Nasipuri, and A. K. Nagar, “Fuzzy relational system for identification of gene regulatory network,” in Proceedings of the World Congress in Computer Science, Computer Engineering, and Applied Computing (WORLDCOMP '09), Las Vegas, Nev, USA, July 2009.

T. A. Al Qazlan, A. Hamdi-Cherif, and C. Kara-Mohamed, “State of the art of fuzzy methods for gene regulatory networks inference,” The Scientific World Journal, vol. 2015, Article ID 148010, 11 pages, 2015.

B. C. Goodwin, Temporal Organization in Cells, Academic Press, New York, NY, USA, 1963.

B. C. Goodwin, “Oscillatory behavior in enzymatic control processes,” in Advances in Enzyme Regulation, G. Weber, Ed., pp. 425–438, Pergamon Press, Oxford, UK, 1965.

J.-L. Gouzé, “Positive and negative circuits in dynamical systems,” Journal of Biological Systems, vol. 6, no. 1, pp. 11–15, 1998.

E. H. Snoussi, “Necessary conditions for multistationarity and stable periodicity,” Journal of Biological Systems, vol. 6, no. 1, pp. 3–9, 1998.

J. J. Tyson and H. G. Othmer, “The dynamics of feedback control circuits in biochemical pathways,” Progress in Theoretical Biology, vol. 5, pp. 1–62, 1978.

A. D. Keller, “Specifying epigenetic states with autoregulatory transcription factors,” Journal of Theoretical Biology, vol. 170, no. 2, pp. 175–181, 1994.

R. Heinrich and S. Schuster, The Regulation of Cellular Systems, Chapman & Hall, New York, NY, USA, 1996.

J. L. Cherry and F. R. Adler, “How to make a biological switch,” Journal of Theoretical Biology, vol. 203, no. 2, pp. 117–133, 2000.

G. Yagil, “Quantitative aspects of protein induction,” Current Topics in Cellular Regulation, vol. 9, pp. 183–235, 1975.

A. Polynikis, S. J. Hogan, and M. di Bernardo, “Comparing different ODE modelling approaches for gene regulatory networks,” Journal of Theoretical Biology, vol. 261, no. 4, pp. 511–530, 2009.

S. Widder, J. Schicho, and P. Schuster, “Dynamic patterns of gene regulation. I. Simple two-gene systems,” The Journal of Theoretical Biology, vol. 246, no. 3, pp. 395–419, 2007.

S. Goutelle, M. Maurin, F. Rougier et al., “The Hill equation: a review of its capabilities in pharmacological modelling,” Fundamental and Clinical Pharmacology, vol. 22, no. 6, pp. 633–648, 2008.

J. N. Weiss, “The Hill equation revisited: uses and misuses,” The FASEB Journal, vol. 11, no. 11, pp. 835–841, 1997.

M. Santillán, “On the use of the Hill functions in mathematical models of gene regulatory networks,” Mathematical Modelling of Natural Phenomena, vol. 3, no. 2, pp. 85–97, 2008.

L. Glass and S. A. Kauffman, “The logical analysis of continuous, non-linear biochemical control networks,” Journal of Theoretical Biology, vol. 39, no. 1, pp. 103–129, 1973.

L. Glass, “Classification of biological networks by their qualitative dynamics,” Journal of Theoretical Biology, vol. 54, no. 1, pp. 85–107, 1975.

E. Plahte, T. Mestl, and S. W. Omholt, “A methodological basis for description and analysis of systems with complex switch-like interactions,” Journal of Mathematical Biology, vol. 36, no. 4, pp. 321–348, 1998.

E. Plahte and S. Kjøglum, “Analysis and generic properties of gene regulatory networks with graded response functions,” Physica D, vol. 201, no. 1-2, pp. 150–176, 2005.

R. Casey, H. de Jong, and J.-L. Gouzé, “Piecewise-linear models of genetic regulatory networks: equilibria and their stability,” Journal of Mathematical Biology, vol. 52, no. 1, pp. 27–56, 2006.

J. Gebert, N. Radde, and G.-W. Weber, “Modeling gene regulatory networks with piecewise linear differential equations,” European Journal of Operational Research, vol. 181, no. 3, pp. 1148–1165, 2007.

E. Plahte, T. Mestl, and S. W. Omholt, “Global analysis of steady points for systems of differential equations with sigmoid interactions,” Dynamics and Stability of Systems, vol. 9, no. 4, pp. 275–291, 1994.

H. El Snoussi and R. Thomas, “Logical identification of all steady states: the concept of feedback loop characteristic states,” Bulletin of Mathematical Biology, vol. 55, no. 5, pp. 973–991, 1993.

A. Hoffmann, A. Levchenko, M. L. Scott, and D. Baltimore, “The IκB-NF-κB signaling module: temporal control and selective gene activation,” Science, vol. 298, no. 5596, pp. 1241–1245, 2002.

S. L. Werner, D. Barken, and A. Hoffmann, “Stimulus specificity of gene expression programs determined by temporal control of IKK activity,” Science, vol. 309, no. 5742, pp. 1857–1861, 2005.

N. Geva-Zatorsky, N. Rosenfeld, S. Itzkovitz et al., “Oscillations and variability in the p53 system,” Molecular Systems Biology, vol. 2, no. 1, Article ID 2006.0033, 2006.

L. Ma, J. Wagner, J. J. Rice, W. Hu, A. J. Levine, and G. A. Stolovitzky, “A plausible model for the digital response of p53 to DNA damage,” Proceedings of the National Academy of Sciences of the United States of America, vol. 102, no. 40, pp. 14266–14271, 2005.

B. T. Hennessy, D. L. Smith, P. T. Ram, Y. Lu, and G. B. Mills, “Exploiting the PI3K/AKT pathway for cancer drug discovery,” Nature Reviews Drug Discovery, vol. 4, no. 12, pp. 988–1004, 2005.

T. M. K. Cheng, S. Gulati, R. Agius, and P. A. Bates, “Understanding cancer mechanisms through network dynamics,” Briefings in Functional Genomics, vol. 11, no. 6, Article ID els025, pp. 543–560, 2012.

L. B. Edelman, J. A. Eddy, and N. D. Price, “In silico models of cancer,” Wiley Interdisciplinary Reviews: Systems Biology and Medicine, vol. 2, no. 4, pp. 438–459, 2010.

N. A. M. Monk, “Oscillatory expression of Hes1, p53, and NF-κB driven by transcriptional time delays,” Current Biology, vol. 13, no. 16, pp. 1409–1413, 2003.

A. Wan and X. Zou, “Hopf bifurcation analysis for a model of genetic regulatory system with delay,” Journal of Mathematical Analysis and Applications, vol. 356, no. 2, pp. 464–476, 2009.

K. Wang, L. Wang, Z. Teng, and H. Jiang, “Stability and bifurcation of genetic regulatory networks with delays,” Neurocomputing, vol. 73, no. 16–18, pp. 2882–2892, 2010.

M. Bodnar, U. Foryś, and J. Poleszczuk, “Analysis of biochemical reactions models with delays,” Journal of Mathematical Analysis and Applications, vol. 376, no. 1, pp. 74–83, 2011.

Z. Wang, Z. Liu, and R. Yuan, “Stability and bifurcation in a gene regulatory network model with delay,” Zeitschrift für Angewandte Mathematik und Mechanik, vol. 92, no. 4, pp. 290–303, 2012.

M. Xiao, W. X. Zheng, and J. Cao, “Stability and bifurcation of genetic regulatory networks with small RNAs and multiple delays,” International Journal of Computer Mathematics, vol. 91, no. 5, pp. 907–927, 2014.

M. Xiao and J. Cao, “Genetic oscillation deduced from Hopf bifurcation in a genetic regulatory network with delays,” Mathematical Biosciences, vol. 215, no. 1, pp. 55–63, 2008.

W. He and J. Cao, “Robust stability of genetic regulatory networks with distributed delay,” Cognitive Neurodynamics, vol. 2, no. 4, pp. 355–361, 2008.

B. D. Aguda, Y. Kim, M. G. Piper-Hunter, A. Friedman, and C. B. Marsh, “MicroRNA regulation of a cancer network: consequences of the feedback loops involving miR-17-92, E2F, and Myc,” Proceedings of the National Academy of Sciences of the United States of America, vol. 105, no. 50, pp. 19678–19683, 2008.

B. D. Aguda, “Modeling microRNA-transcription factor networks in cancer,” in MicroRNA Cancer Regulation: Advanced Concepts, Bioinformatics and Systems Biology Tools, U. Schmitz, O. Wolkenhauer, and J. Vera, Eds., vol. 774 of Advances in Experimental Medicine and Biology, pp. 149–167, Springer, Dordrecht, The Netherlands, 2013.

A. S. Ribeiro, J. J. Grefenstette, and S. A. Kauffman, “Modeling gene regulatory networks with delayed stochastic dynamics,” in Handbook of Research on Computational Methodologies in Gene Regulatory Networks, S. Das, D. Caragea, S. M. Welch, and W. H. Hsu, Eds., pp. 198–218, Medical Information Science Reference, New York, NY, USA, 2010.

N. MacDonald, “Time lag in a model of a biochemical reaction sequence with end product inhibition,” Journal of Theoretical Biology, vol. 67, no. 3, pp. 549–556, 1977.

I. R. Epstein and Y. Luo, “Differential delay equations in chemical kinetics. Nonlinear models: the cross-shaped phase diagram and the Oregonator,” The Journal of Chemical Physics, vol. 95, no. 1, pp. 244–254, 1991.

J. Paulsson, “Models of stochastic gene expression,” Physics of Life Reviews, vol. 2, no. 2, pp. 157–175, 2005.

S. Rosenfeld, “Stochastic cooperativity in non-linear dynamics of genetic regulatory networks,” Mathematical Biosciences, vol. 210, no. 1, pp. 121–142, 2007.

A. Raj and A. van Oudenaarden, “Nature, nurture, or chance: stochastic gene expression and its consequences,” Cell, vol. 135, no. 2, pp. 216–226, 2008.

P. S. Swain, M. B. Elowitz, and E. D. Siggia, “Intrinsic and extrinsic contributions to stochasticity in gene expression,” Proceedings of the National Academy of Sciences of the United States of America, vol. 99, no. 20, pp. 12795–12800, 2002.

K. M. Herbert, A. La Porta, B. J. Wong et al., “”Sequence resolved detection of pausing by single RNA polymerase molecules,” Cell, vol. 125, no. 6, pp. 1083–1094, 2006.

A. S. Ribeiro, O.-P. Smolander, T. Rajala, A. Häkkinen, and O. Yli-Harja, “Delayed stochastic model of transcription at the single nucleotide level,” Journal of Computational Biology, vol. 16, no. 4, pp. 539–553, 2009.

M. Voliotis, N. Cohen, C. Molina-París, and T. B. Liverpool, “Fluctuations, pauses and backtracking in DNA transcription,” Biophysical Journal, vol. 94, no. 2, pp. 334–348, 2008.

A. S. Ribeiro, “Stochastic and delayed stochastic models of gene expression and regulation,” Mathematical Biosciences, vol. 223, no. 1, pp. 1–11, 2010.

M. J. Dunlop, R. S. Cox III, J. H. Levine, R. M. Murray, and M. B. Elowitz, “Regulatory activity revealed by dynamic correlations in gene expression noise,” Nature Genetics, vol. 40, no. 12, pp. 1493–1498, 2008.

M. B. Elowitz and S. Leibier, “A synthetic oscillatory network of transcriptional regulators,” Nature, vol. 403, no. 6767, pp. 335–338, 2000.

M. B. Elowitz, A. J. Levine, E. D. Siggia, and P. S. Swain, “Stochastic gene expression in a single cell,” Science, vol. 297, no. 5584, pp. 1183–1186, 2002.

S. Rosenfeld, “Stochastic oscillations in genetic regulatory networks: application to microarray experiments,” Eurasip Journal on Bioinformatics and Systems Biology, vol. 2006, Article ID 59526, 2006.

M. Kærn, T. C. Elston, W. J. Blake, and J. J. Collins, “Stochasticity in gene expression: from theories to phenotypes,” Nature Reviews Genetics, vol. 6, no. 6, pp. 451–464, 2005.

D. Bratsun, D. Volfson, L. S. Tsimring, and J. Hasty, “Delay-induced stochastic oscillations in gene regulation,” Proceedings of the National Academy of Sciences of the United States of America, vol. 102, no. 41, pp. 14593–14598, 2005.

E. Zavala and T. T. Marquez-Lago, “Delays induce novel stochastic effects in negative feedback gene circuits,” Biophysical Journal, vol. 106, no. 2, pp. 467–478, 2014.

T. S. Gardner, C. R. Cantor, and J. J. Collins, “Construction of a genetic toggle switch in Escherichia coli,” Nature, vol. 403, no. 6767, pp. 339–342, 2000.

H. El Samad, M. Khammash, L. Petzold, and D. Gillespie, “Stochastic modelling of gene regulatory networks,” International Journal of Robust and Nonlinear Control, vol. 15, no. 15, pp. 691–711, 2005.

M. Barrio, K. Burrage, A. Leier, and T. Tian, “Oscillatory regulation of Hes1: discrete stochastic delay modelling and simulation,” PLoS Computational Biology, vol. 2, no. 9, article e117, 2006.

X. Cai, “Exact stochastic simulation of coupled chemical reactions with delays,” Journal of Chemical Physics, vol. 126, Article ID 124108, 2007.

M. Jansen and P. Pfaffelhuber, “Stochastic gene expression with delay,” Journal of Theoretical Biology, vol. 364, pp. 355–363, 2015.

N. Friedman, M. Linial, I. Nachman, and D. Pe'er, “Using Bayesian networks to analyze expression data,” Journal of Computational Biology, vol. 7, no. 3-4, pp. 601–620, 2000.

D. Heckerman, “A tutorial on learning with Bayesian networks,” Microsoft Research Technical Report MSR-TR-95-06, 1996.

A. V. Werhli and D. Husmeier, “Reconstructing gene regulatory networks with Bayesian networks by combining expression data with multiple sources of prior knowledge,” Statistical Applications in Genetics and Molecular Biology, vol. 6, article 15, 2007.

P. T. Spellman, G. Sherlock, M. Q. Zhang et al., “Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization,” Molecular Biology of the Cell, vol. 9, no. 12, pp. 3273–3297, 1998.

R. J. Prill, P. A. Iglesias, and A. Levchenko, “Dynamic properties of network motifs contribute to biological network organization,” PLoS Biology, vol. 3, no. 11, article e343, 2005.

Y. N. Kyrychko and K. B. Blyuss, “Global properties of a delayed SIR model with temporary immunity and nonlinear incidence rate,” Nonlinear Analysis: Real World Applications, vol. 6, no. 3, pp. 495–507, 2005.

H. B. Tiedemann, E. Schneltzer, S. Zeiser et al., “From dynamic expression patterns to boundary formation in the presomitic mesoderm,” PLoS Computational Biology, vol. 8, no. 6, Article ID e1002586, 2012.

C. Li and J. Wang, “Quantifying cell fate decisions for differentiation and reprogramming of a human stem cell network: landscape and biological paths,” PLoS Computational Biology, vol. 9, no. 8, Article ID e1003165, 2013.

J. D. Murray, Mathematical biology. I. An Introduction, Springer, New York, NY, USA, 2002.

K. B. Blyuss and Y. N. Kyrychko, “Stability and bifurcations in an epidemic model with varying immunity period,” Bulletin of Mathematical Biology, vol. 72, no. 2, pp. 490–505, 2010.

Y. N. Kyrychko, K. B. Blyuss, and E. Schöll, “Amplitude death in systems of coupled oscillators with distributed-delay coupling,” European Physical Journal B, vol. 84, no. 2, pp. 307–315, 2011.

Y. N. Kyrychko, K. B. Blyuss, and E. Schöll, “Amplitude and phase dynamics in oscillators with distributed-delay coupling,” Philosophical Transactions of the Royal Society of London A. Mathematical, Physical and Engineering Sciences, vol. 371, no. 1999, Article ID 20120466, 2013.

Y. N. Kyrychko, K. B. Blyuss, and E. Schöll, “Synchronization of networks of oscillators with distributed delay coupling,” Chaos, vol. 24, no. 4, Article ID 043117, 2014.

M. Lakshmanan and D. V. Senthilkumar, Dynamics of Nonlinear Time-Delay Systems, Springer Series in Synergetics, Springer, Berlin, Germany, 2010.

R. Suresh, D. V. Senthilkumar, M. Lakshmanan, and J. Kurths, “Transition to complete synchronization and global intermittent synchronization in an array of time-delay systems,” Physical Review E, vol. 86, no. 1, Article ID 016212, 2012.

J. Sieber, P. Kowalczyk, S. J. Hogan, and M. Di Bernardo, “Dynamics of symmetric dynamical systems with delayed switching,” Journal of Vibration and Control, vol. 16, no. 7-8, pp. 1111–1140, 2010.

Z. Li, P. Li, A. Krishnan, and J. Liu, “Large-scale dynamic gene regulatory network inference combining differential equation models with local dynamic Bayesian network analysis,” Bioinformatics, vol. 27, no. 19, Article ID btr454, pp. 2686–2691, 2011.

O. Folger, L. Jerby, C. Frezza, E. Gottlieb, E. Ruppin, and T. Shlomi, “Predicting selective drug targets in cancer through metabolic networks,” Molecular Systems Biology, vol. 7, article 501, 2011.

K. Raza and R. Jaiswal, “Reconstruction and analysis of cancer-specific gene regulatory networks from gene expression profiles,” International Journal on Bioinformatics and Biosciences, vol. 3, no. 2, pp. 25–34, 2013.

T.-T. Zhou, “Network systems biology for targeted cancer therapies,” Chinese Journal of Cancer, vol. 31, no. 3, pp. 134–141, 2012.

G. A. Banks, R. J. Roselli, R. Chen, and T. D. Giorgio, “A model for the analysis of nonviral gene therapy,” Gene Therapy, vol. 10, no. 20, pp. 1766–1775, 2003.

A. Tsygvintsev, S. Marino, and D. E. Kirschner, “A mathematical model of gene therapy for the treatment of cancer,” in Mathematical Methods and Models in Biomedicine, U. Ledzewicz, H. Schattler, A. Friedman, and E. Kashdan, Eds., Lecture Notes on Mathematical Modelling in the Life Sciences, pp. 367–385, Springer, New York, NY, USA, 2013.