Fractional Dynamics of Globally Slow Transcription and Its Impact on Deterministic Genetic Oscillation

In dynamical systems theory, a system which can be described by differential equations is called a continuous dynamical system. In studies on genetic oscillation, most deterministic models at early stage are usually built on ordinary differential equations (ODE). Therefore, gene transcription which is a vital part in genetic oscillation is presupposed to be a continuous dynamical system by default. However, recent studies argued that discontinuous transcription might be more common than continuous transcription. In this paper, by appending the inserted silent interval lying between two neighboring transcriptional events to the end of the preceding event, we established that the running time for an intact transcriptional event increases and gene transcription thus shows slow dynamics. By globally replacing the original time increment for each state increment by a larger one, we introduced fractional differential equations (FDE) to describe such globally slow transcription. The impact of fractionization on genetic oscillation was then studied in two early stage models – the Goodwin oscillator and the Rössler oscillator. By constructing a “dual memory” oscillator – the fractional delay Goodwin oscillator, we suggested that four general requirements for generating genetic oscillation should be revised to be negative feedback, sufficient nonlinearity, sufficient memory and proper balancing of timescale. The numerical study of the fractional Rössler oscillator implied that the globally slow transcription tends to lower the chance of a coupled or more complex nonlinear genetic oscillatory system behaving chaotically.


Introduction
The design and construction of genetic circuits are of great importance to the nascent field of synthetic biology [1]. A vital design principle for synthetic genetic circuits is to insure the capability of self-sustained oscillations for adapting the biological rhythms or environmental cycles. For example, the Goodwin oscillator which is considered to be the simplest genetic oscillator and the basis of the repressilators [1,2] has been used as a minimal model to interpret the circadian rhythms occurring in gene's negative autoregulation [3,4]. Usually, one may expect to describe a basic genetic circuit by using a minimal dynamical model (with as few equations as possible), for the purpose of simplicity. The variables in such model represent the quantities of several key products in the circuit. However, as we have known, even the simplest genetic regulation includes complex intermediate processes like transcription, transportation of RNA, RNA splicing, RNA capping, translation, transportation of mature protein and other steps of post-translational modification. Therefore, a minimal model (sometimes with only a single equation) often lacks power to cover such complex intermediate processes in a regular timescale. This can be seen from the case that a Goodwin oscillator which requires an unrealistic high Hill coefficient (larger than 8) for the destabilization of a fixed point and generating limit cycle oscillations [5,6].
However, by altering the timescale of gene transcription and introducing slow dynamics (e.g. considering the time lag in the protein transportation), one can readily obtain the desirable dynamical behaviors by using minimal models. A common method to achieve slow dynamics is introducing explicit time delay. In such way, sufficient time delay is considered to be one of the general requirements for sustained oscillations [7]. Another method is inserting an additional equation for lagging fast change in protein level, which plays a dynamical role similar to explicit time delays or to transport equations [8]. The two methods above display certain sorts of memory effects in gene transcription [7,8]. These optimizations in timescale are designed merely for a specific intermediate product (e.g. protein) or for a specific intermediate biochemical step. We thus regard the above timescale changes as non-global.
In this paper, we would like to introduce the globally slow transcription which can be described by fractional differential equations (FDE). This idea is based on the recent evidence that discontinuous transcription may be more common than continuous transcription. In the discontinuous transcription case, by inserting silent intervals into neighboring transcriptional events which are presupposed to be continuous, and appending an inserted silent interval to the end of the preceding event, we aim to establish that the running time for an intact transcriptional event increases and the transcription system thus shows the property of globally slow transcription. Effects of such globally slow transcription are investigated in a minimal Goodwin oscillator and a Rössler oscillator, both of which are well-known in genetic regulation.

Definitions of Fractional Calculus and Memory Weighted via the Convolution Kernel
In spite of the existence of different definitions for the fractional derivatives, the fractional integral is the common foundation of fractional calculus. The notion of the left-side fractional integral operator with order l[R z is in fact an extension of the Cauchy's formula for repeated integrals which replaces the l-fold integrals of a function f by a simple convolution: where tw0, g(t)~t l{1 C(l) , and C( : ) is the Gamma function. Moreover, under certain reasonable assumptions there exists lim we then obtain an identity operator denoted by [9] Considering that in real world applications, the evolution of a general dynamical system governed by the principle of causality is apriori timeirreversible, we use the left-side integral/derivative operators and the initial time t 0~0 throughout this paper. The derivative operator under the Caputo definition is expressed as follows [10]: (a) If l~n[N 0 , we have the integer order derivatives: where D 1 f (t)~lim Dt?0 , demonstrating the local property at a given time point t. In particular, there exists where n~qlr and n{1vlvn. In particular, when 0vlv1, we have From Eq. (3) and Eq. (4) we can see that the difference between the integer order derivative operator and the fractional (non-integer order) derivative operator is caused by the integral 0 J n{l t (0vn{lv1) which endows the operator C 0 D l t f (t) with the non-local property because the information of the entire integral interval 0,t ½ is involved. This is why we regard the fractional integral as the common foundation of fractional calculus.
In dynamical systems theory, memory effects are usually described by explicit non-local terms about the state variables. In Eq. (1), we know g : 0,? ½ Þ?R. If f : 0,? ½ Þ?R, according to the commutative property of convolution, we have Then, Eq. (1) can be written as It is obvious that the fractional integral Eq. (8) reflects a weighted average of delays f (t{s) via a specific convolution kernel function g(s)~s l{1 C(l) . The values of the weighting function g(s) change with s under different order l are illustrated in Fig. 1. In the circumstance where lw1, the memory is enhanced with s increasing, while states closed to present are given little weights. For 0vlv1, g(s) decreases with increasings, leading to the ''fading memory'' property with which the importance of the past state f (t{s) fades out. In this scenario, the larger values of l provide slower decay of g(s). In the limit as l approaches zero, g(s) approaches the Dirac delta function, leading to the identity operator (zero order calculus). In contrast, in the limit as l approaches 1, the weighting function g(s)~s

Transcriptional Discontinuity Builds the Link between Gene Transcription and Fractional Calculus
In early models of genetic regulation, the system of gene transcription is presupposed to be a continuous dynamical system and thus its long-time behavior can be described by using an ordinary differential equation (ODE) which contains a Michaelis-Menten (MM) mRNA synthesis term. For example, the onevariable Goodwin model is where a is the degradation coefficient and m is Hill coefficient [6,11]. By examining the two-variable extension of Eq. (9): and the three-variable extension: Griffith [6] found that the addition of intermediate steps (with the non-MM synthesis term for an intermediate product such like protein) tends to make the oscillation be achieved more easily. This finding makes us realize that the transcription equation which contains the MM mRNA synthesis term is the fundamental part of a basic autoregulation model, while other intermediate steps (with non-MM synthesis) play the role of introducing slow dynamics for certain intermediate variables, as described in [8]. The above suggests that one may expect to compact any Goodwin model to a single equation for gene transcription like Eq. (9).
Eq. (9) is a single ODE which reflects a continuous dynamical system. The increment of continuous state variable at the time point t can be expressed in difference form: or in differential form: where F (x,t) is the sum rate of synthesis rate minus degradation rate, and the time increment Dt or dt is an extremely tiny value. However, a recent molecular experiment shows that both continuous transcription and discontinuous transcription exist in yeast gene expression [12]. Moreover, more and more direct evidence show that the process of transcription for most genes is interspersed by ''gene-on'' and ''gene-off'' states by turns, in both prokaryotes and eukaryotes [13][14][15]. A silent interval (in which genes switch to ''off'' state and enter a refractory period) lying between two neighboring transcriptional events makes gene transcription not so continuous, and therefore, dynamics in such situation is described to be temporally discontinuous [15]. In this sense, if we remove those silent intervals, the presupposed continuous dynamical system of gene transcription is recovered. In order to describe the temporally discontinuous gene transcription by using differential equations, we impose a silent interval upon its preceding transcriptional event. The schedule of an intact transcriptional event is then extended with a refractory period being appended as the end part, implying a relative slow transcriptional dynamics for that the total running time of an intact transcriptional event increases. With this treatment, two neighboring transcriptional events can be connected smoothly and continuously without an interval (Fig. 2). The idea of inserting silent intervals into continuous events can be traced back to the case where a trapping event (the Brown particle is temporarily immobile) is inserted into two neighboring jumps of continuous time random walk (Fig. 2) and eventually leads to fractional dynamics [16][17][18][19]. Different from this stochastic dynamics whose solution is usually described by transition probability, in our deterministic model, imposing a time interval upon a transcriptional event and increasing the running time of an intact transcriptional event will make each variable increment dx occur within a larger time increment rather than within the original dt. Since the original dt is an extremely tiny value (e.g. dt~0:0001 in dimensionless form), we take Jumarie's notion of time increment (dt) l , where 0vlv1 leads to (dt) l wdt while lw1 leads to (dt) l vdt [20]. By setting (dt) l with 0vlv1 to be the larger time increment, we obtain the analog of Eq. (10b) which reads.
The relation between the fractional difference and the finite difference has been given as or in differential form: where l[(0, 1) [20][21][22]. Multiplying both sides of Eq. (11) by C(1zl) and taking account of Eq. (12b), we obtain Hence, with time evolving, we have a l-order differential equation: The symbol of l-order derivative d l x (dt) l makes one recall the ageold issue presented in the communications between L'Hôspital and Leibniz: what if the order is 1=2 (see [23] and the Preface of [10]). Mathematicians have been inspirited by this story for over 300 years and their endeavors have led to a variety of nonequivalent definitions for the fractional order derivatives. By taking Caputo's definition of the left-side fractional derivative and normalizing the constant coefficient C(1zl) in Eq. (14) to unity, we obtain the generalized form: is differentiable, by calculating (1{l)-order derivatives of the both sides of Eq. (15a), we can obtain an equivalent equation: where _ x x represents the first order derivative of x. Different from the traditional ordinary differential equation _ x x~F(x,t) in which the product rate _ x x is locally determined by F(x,t), Eq. (15) shows long-term memory because the fractional derivative operator is non-local. Take C 0 D 1{l t in Eq. (15b) for example, according to Eq. (5) and the commutative property of convolution, we know that Since F (x,t) is the sum of synthesis rate minus degradation rate, we define the first order derivative F (1) (x,t) as the ''acceleration of product gain''. The convolution kernel g(s)~s l{1 C(l) is used to weight F (1) (x,t{s) from nonce to remote history with s increasing from zero to current time point t. As depicted in Fig. 1, with normal timescale (l~1), weights are always equal, implying that each transcriptional event launches de novo. Under this condition, the traditional ordinary differential equation _ x x~F (x,t) is recovered. In the case where lw1, the farthest transcriptional event gives the highest impact while the weight given by the nearest transcriptional event is nearly close to zero. This case is in contradiction with regular physiological phenomena. Besides, l can be infinitely large under the condition of lw1. Therefore, we exclude lw1 in our study. In contrast, the case in which 0vlv1 shows the property of ''fading memory''. With such memory, a current event carries information form the preceding events, especially the nearest one. This property seems to be consistent with the observed phenomenon that transcriptional reinitiation is more common than de novo initiation in a discontinuouslytranscribed gene [13]. However, the mechanisms of such transcriptional memory are not very clear by now. A feasible explanation may be that the loop scaffold forming in gene transcription would retain certain enzymes (or regulatory factors) of the preceding transcriptional events [24,25] and sterically hinder the new recruited enzymes to take their place. Therefore, a subsequent transcriptional event remembers the preceding transcriptional events and takes a reinitiation rather than initiating de novo. The likelihood of reinitiation diminishes as the interval time elapses, and such relation can fit a simple exponential decay function [13]. Since l[(0, 1), we simply use a regular decay expression l~e {rDhD to represent the probability for a gene still ''surviving'' with the capability of transcriptional reinitiation after a time length DhD (time span of the silent interval), while 1{l is the probability for the gene decaying to de novo initiation. The unknown system constant r represents the probability per unit time for a still ''surviving'' gene to decay [26]. In this way, a relative small DhD indicates a relative large l, resulting in a relative small (dt) l . If the silent interval does not exist, namely, in the limit as DhD approaches zero, l will approaches one and (dt) l will approaches dt, leading to the recovery of the early stage models for the continuous transcription.

Numerical Studies of Fractional Dynamics in two Illustrative Examples -a minimal goodwin oscillator and a Rö ssler Oscillator
The first example is the Goodwin oscillator which was introduced originally in 1950s to simulate physiological oscillations in a closed loop with negative feedback.
The one-variable Goodwin model is expressed as Eq. (9). This model reflects by default a fast dynamics in which the products will be put into the feedback loop very quickly and then participate in the reactions immediately (Fig. 3A). However, this model does not generate sustained oscillation [6]. When the explicit time delay is introduced for generating sustained oscillation, a slow dynamics is achieved and the model becomes a delay Goodwin oscillator (Fig. 3B). If transcriptional discontinuity is considered, by involving explicit time delay and the deduced Eq. (15a), we obtain a ''fractional delay Goodwin oscillator'' which reads where hƒ0 is the time delay and the symbol C 0 D l t with l[(0, 1) denotes the Caputo fractional derivative operator. Since both the time delay and the fractional operator are non-local, Eq. (17) describes a ''dual memory'' system.
In this model, we assume that the silent interval retards the inhibitor produced in a transcriptional event to act on the gene promoter until the transcriptional reinitiation starts. Under this assumption, the time delay can be set to be equal to h (hƒ0) when DhD denotes the time span of the silent interval ( Fig. 2C and Fig. 3C).
In order to avoid negative solution when time delay is involved in the Goodwin oscillator [27], the conditions required for nonnegative solution should be established. Since the Goodwin oscillator can often be developed into several variants (e.g. the degradation term {ax(t) can be replaced by an MM expression  [7,28]), we investigate the existence and the uniqueness of the nonnegative solution for the generalized form of the fractional delay Goodwin oscillator (see Appendix S1).
The Simulink block diagram of Eq. (17) is depicted in Fig. 4. The fractional integrator used in this section is the well-established Oustaloup recursive filter which has been proved to fit well to the fractional operators and has been widely used in control systems [29][30][31]. It is suggested that because the same orders of the numerator and the denominator in the ordinary Oustaloup filter may cause algebraic loops in simulation, a low-pass filter must be appended to the Oustaloup filter to avoid such disadvantage; the stiff equation solver ode23tb is selected to ensure high efficiency and accuracy [32]. In simulation, we fix the parameter a~0:2 in the degradation term.
The second example is the Rössler oscillator which is originally used to study chaotic kinetics in biochemical system [33]. Novak and Tyson [7] used this model to describe the activator amplification with two negative feedback loops in parallel. From the viewpoint of mathematical modeling, the common phenomenon of coupling different genetic oscillator motifs in gene regulation would tend to cause complex non-linear oscillations like chaos. However, the deterministic chaos has not been detected more often in real data from experiments, and Novak and Tyson speculated that such chaos might be swamped by white noise and averaged out in large populations of cells [7]. This speculation is passable if stochastic factors are introduced. However, since all we have discussed so far are deterministic models, can we give an alternative explanation about the lack of chaos in gene regulation merely from the angle of determinism?
For convenience, we use the classic Rössler system (the ''model of a model''; [34]) to illustrate how the globally slow transcription affects the system's dynamical behavior. The parameters are set by a~0:15, b~0:2 and c~10 with sampling period Dt~0:1s for time series, according to the reference [35]. Because of the global property of (dt) l , all three individual equations of Eq. (18) share a commensurate order just like the way given by the traditional commensurate first order Rössler system (but with l[(0, 1)). The numerical method for fractional calculus is the Oustaloup recursive filter mentioned above. The largest Lyapunov exponent (LLE) is calculated by using the method of Rosenstein et al. [35], and the embedded delay and dimension are evaluated by using the methods of Kim et al. [36] and Kugiumtzis [37].

Results and Discussion
By inserting silent intervals into consecutive transcriptional events and globally replacing dt by a larger increment (dt) l , we make the early stage models change from ODE to FDE. The impact of the order fractionization on genetic oscillation can be shown by the Goodwin oscillator and the Rössler oscillator.
For the fractional delay Goodwin oscillator, the numerical results show that no sustained oscillations occur when m~1; however, with fixed l and h, the increase of m leads to the destabilization of the steady state (data not shown). In the case of m~2 (Fig. 5A and Fig. 5B), when h is specified, the increase of l leads to the stability loss. The impact of time delay is clearly shown in Fig. 5B-5D. When m~2 and l~0:9, setting the delay time h~{15 generates damped oscillation which is on the way to Step block is used for assigning the initial value in the beginning step, and the Outport block is used for collecting the simulation results. doi:10.1371/journal.pone.0038383.g004 steady state, while h~{20 produces periodical sustained oscillation. In this case, by plotting the rate of degradation or synthesis versus the present values x(t), the trajectory of the time-delayed synthesis rate with h~{15 is attracted to the intersection point (the steady state of x) of the degradation rate and the non-delayed synthesis rate (Fig. 5C). With an extending delay of h~{20, the time-delayed loop overshoots and undershoots the steady state, indicating the periodical sustained oscillation of x(t) (Fig. 5D).
Novak and Tyson [7] proposed four general requirements for oscillation in gene regulation: (1) negative feedback loop; (2) sufficient nonlinearity; (3) sufficient time delay; (4) proper balancing of timescale (namely, a in the degradation term of Eq. (17) must not be too large). However, by introducing fractional dynamics, the decrease of the order l would tend to make the oscillation more stable, even though a sufficient time delay is reached. This can be clearly shown in Fig. 5A and Fig. 5B: when sufficient nonlinearity (m~2), sufficient time delay (h~{20) and proper balancing of timescale (a~0:2) are satisfied, the decrease of l (from 0.9 to 0.7) makes the oscillatory limit cycle become a fixed point attractor. Since the order l (as shown in Fig. 1) can be used to indicate the strength of memory in fractional dynamics, we suggest the use of the term ''sufficient memory'' for the 3rd requirement rather than using only ''sufficient time delay'' when fractional dynamics is involved.
For the Rössler oscillator, the traditional commensurate first order model shows the classic unimodal folded chaos with LLE of +0.0995 (Fig. 6F). When fractional dynamics is involved, it is clearly shown that the dynamical behavior changes from nonperiodic (chaotic) motion (LLE.0) to periodic motion (LLE = 0) with l decreasing (Fig. 6F). Some points that seem to be outliers at around l&0:993 (Fig. 5F) are attributed to the meeting of the period-three window which is embedded in the chaotic region ( Fig. 6C and Fig. 6E). The existence of the period three window in the diagram of period doubling bifurcation (Fig. 6E) implies that the chaotic dynamics of the fractional Rössler model can be interpreted through the Sharkovskii order [38] or Li-Yorke theorem [39]. The critical value for the route to chaos is l C &0:9845, and the region corresponding to the order interval 0:9845ƒlƒ1 is defined as the chaotic region. Therefore, if we simply consider a uniform distribution of l within the order  17) with h~{15 indicates a damped oscillation on the way to steady state, while the case with h~{20 corresponds to sustained oscillation. C. Plots of the synthesis rate (without delay), the degradation rate and the rate of synthesis containing delay, against the present value x(t). The blue curve indicates the damped oscillation. D. Plots of the synthesis rate (without delay), the degradation rate and the rate of synthesis containing delay, against the present value x(t). The blue circle indicates sustained oscillation. doi:10.1371/journal.pone.0038383.g005 interval (0,1), the fractional Rössler oscillator with originally specific parameters will show a probability of more than 0.98 to behave non-chaotically. This deterministic fractional dynamics may provide an alternative explanation for understanding why the deterministic chaos in gene regulation has not been detected more often in real data from experiments [7].
In conclusion, a transcription equation with a term representing MM mRNA synthesis is the core of a basic autoregulation model. Therefore, any basic autoregulation circuit can be minimized to such a single equation. Traditionally, slow dynamics can be achieved via setting explicit time delay to the state variable of a minimal model. In this study, we propose a globally slow transcription on the basis of the recent observation that discontinuous transcription may be more common than continuous transcription. By inserting silent intervals into neighboring transcriptional events which are presupposed to be continuous, and appending an inserted silent interval to the end of the preceding event, we establish that the running time for an intact transcriptional event increases. By globally replacing the original time increment for each state increment by a larger one, we obtain a fractional model for gene transcription. With the assumption that the silent interval hinders the product to trigger the reinitiation of gene transcription and hence causes time delay, we construct a fractional delay Goodwin oscillator. Since both time delay and fractional operator are non-local, such new type of Goodwin oscillator is a ''dual memory'' system. To avoid negative solution when time delay is involved, the existence and the uniqueness of the non-negative solution for the generalized form of the fractional delay Goodwin oscillator are also studied. The numerical studies show that the explicit time delay tends to destabilize the steady state, while the fractionization of the order tends to make the system stable. This result makes us realize that the requirement ''sufficient time delay'' for genetic oscillation is not sufficient and should be changed to ''sufficient memory'' when fractional dynamics is involved. When we examine another wellknown genetic oscillator -the Rössler oscillator which describes the activator amplification coupled with two negative feedback loops in parallel, the diagram of period doubling bifurcation against the orders reveals that the globally slow dynamics induced via discontinuous transcription tends to lower the chance of a coupled or more complex nonlinear genetic oscillatory system behaving chaotically.

Supporting Information
Appendix S1 Existence and uniqueness of the nonnegative solution for the generalized form of the fractional delay Goodwin oscillator. (DOC)