Effects of noise and time delay on E2F's expression level in a bistable Rb-E2F gene's regulatory network.

The bistable Rb-E2F gene regulatory network plays a central role in regulating cellular proliferation-quiescence transition. Based on Gillespie's chemical Langevin method, the stochastic bistable Rb-E2F gene's regulatory network with time delays is proposed. It is found that under the moderate intensity of internal noise, delay in the Cyclin E synthesis rate can greatly increase the average concentration value of E2F. When the delay is considered in both E2F-related positive feedback loops, within a specific range of delay (3-13) hr , the average expression of E2F is significantly increased. Also, this range is in the scope with that experimentally given by Dong et al. [65]. By analysing the quasi-potential curves at different delay times, simulation results show that delay regulates the dynamic behaviour of the system in the following way: small delay stabilises the bistable system; the medium delay is conducive to a high steady-state, making the system fluctuate near the high steady-state; large delay induces approximately periodic transitions between high and low steady-state. Therefore, by regulating noise and time delay, the cell itself can control the expression level of E2F to respond to different situations. These findings may provide an explanation of some experimental result intricacies related to the cell cycle.

INTRODUCTION

Noise and time delay are ubiquitous in various biological processes, existing over multiple scales of biological systems from the population level to the molecular level. Noise might be perceived as a nuisance due to its effects on molecular events. However, it modulates the system's dynamics [1, 2]. Noise permeates biological networks at all levels, from the basic cellular to population levels, where it causes inevitable effects to the system [3, 4, 5, 6, 7, 8].

It is shown that cellular processes such as gene expression, signalling pathways, circadian rhythms, etc., when exposed to noise, deal with it differently [9, 10]. It makes it rather hard to give a general effect of noise precisely. Experimentally, bimodal population distributions of cells demonstrate that noise can cause switching on and off states in genetic toggle switches to determine the cell’s fate in a bistable system [11, 12].

In a stable genetic regulatory network (GRN) maintained by positive feedback, intrinsic noise impacts the unstable intermediate states [13]. Stochastic resonance (SR), in which some noise level is added to a weak signal in a non‐linear system, enhances weak signals' detection and transmission [14, 15]. The noise‐enhanced stability (NES) phenomenon has also been studied to show how noise can stabilize the system's fluctuating metastable states [16, 17, 18, 19, 20, 21]. Zeng and Wang [22] applied the NES with coloured noise and showed how tumour cell growth could be stabilised. Thus, it is a clear indication that noise [23] has a significant role in modulating a system's dynamics.

On the other hand, the time delay [24, 25, 26] is one of the essential aspects featured in GRNs [7, 8, 13, 27]. Time delays describe the propagation phenomena, material transfer in intercommoned systems, and data transmission in communicating systems in many natural and physical setups [28]. Typically, time delays deem to complicate the system’s dynamics by causing oscillations, instability, and poor control performances [29, 30]. In the cell cycle, delays are shown to affect the cell cycle [31] by inducing cell arrest and apoptosis through the p53 network upon sensing damaged DNA [32]. These and many other existing works, such as [33, 34, 35] to mention a few, have analysed delayed systems deterministically. However, there are increasing attempts to study the behaviour of stochastic systems with the time delay [36, 37, 38]. Predominantly, most of the indagations are based on modifications on the Chemical Langevin Equations (CLEs) to account for time delay [36, 38, 39, 40]. In fact, Bratsun et al., [38] analysed the combined effects of time delay and intrinsic noise on gene regulation and noticed that the two entities' variability was almost the same. Also, Brett and Galla [41] using the Brusselator model showed that in the presence of noise and delay, the system not only fluctuates with small amplitudes about the deterministic fixed point, but intermittent spikes are also invoked. By using a variety of delay models, Jiang et al. show that Markov’s model learning by the neural network can accurately reflect the stochastic dynamics across the parameter space [42]. Due to delays in gene regulatory networks in an oscillating system, noise can enhance robustness to maintain the oscillating expression pattern [13].

Human cell cycle progression, apoptosis, quiescence, and proliferation have been well‐studied to examine the cell's fate in the case of anti‐mitogens such as DNA damage that would otherwise result in detrimental conditions such as cancer [43, 44]. The cell can ensure the fidelity of DNA replication processes before entry and progression from the G1 to the S phase through several signalling pathways [45]. Signalling events controlled by the Rb/E2F pathway are conserved in many eukaryotes, including mammals, as key regulators of the G1/S checkpoint [46, 47]. E2F is a family of transcription factors that activates many growth target genes that drive entry and advancement through the cell cycle [48, 49]. In contrast, retinoblastoma (Rb) is a tumour suppressor and transcription factor regulating proliferation in almost all cell types [43, 50].

Upon detection of DNA damages or insufficiencies, Rb suppresses cell cycle progression by inhibiting the proliferation activity of E2F [43, 51, 52, 53]. Moreover, CDK Inhibitors such as p27, p21 can also promote quiescence by blocking Rb's phosphorylation, allowing it to bind to E2F hence suppressing cell cycle procession [54, 55]. Myc, Cyclin D/CDK4, 6, and Cyclin E/CDK2 can mediate Rb's phosphorylation, preventing it from binding to E2F [56, 57]. Therefore, this will enhance the activation of E2F, leading to the progression of the cell cycle. Based on rather intricate regulatory networks comprising feedback loops, feedforward loops, or autoregulatory loops [58] involved, the E2F‐Rb homoeostatic mechanism in the cell cycle is dependent on external serum or external and internal random disturbance [59]. Yao et al. [60] described a bistable Rb‐E2F switch mechanism in which inhibitory and growth signals are converted into a binary OFF/ON signal by E2F in an ‘alright or none’ manner leading to quiescence or proliferation [43, 53].

Generally, the role of time delay is rather complex and system dependent. It can induce either on/off switching of oscillation states and create coherent bistability that can aid in the decision‐making of the system [61]. Since E2F plays a vital role in the cell cycle process, it is essential to know how noise and time delay regulate this gene expression.

In this article, the effects of internal noises and time delay on the expression of transcription factor E2F are studied in the bistable Rb‐E2F GRN. In section 2, based on Gillespie's CLEs method, the stochastic differential equations of bistable Rb‐E2F GRN with time delay are proposed. The effects of noise and time delay on the expression of transcription factor E2F are discussed in section 3. We end with the conclusions of the results and the possible biological significance in section 4.

MODEL AND METHOD

When the internal noise is ignored, the E2F regulated cell‐cycle network can be described by Yao's model [60]. However, noise inevitably exists in living cell systems. Therefore, stochastic equations are more valid than deterministic equations. Here we only consider the internal noise, which originates from the random birth‐death chemical species. In this article, Gillespie's stochastic simulation method is adopted, which is mainly based on the reaction rate to determine which reaction will happen next. According to this method, we introduce the number of intracellular E2F molecules as Mole . Similarly, Mole , Mole , Mole , Mole , Mole , Mole are the number of intracellular molecules of Myc, Rb, phosphorylated Rb, CycD‐Cdk4/6 complex, Rb/E2F complex, and CycE‐Cdk2 complex, respectively. Then, the concentrations of the genes involved in the E2F regulated cell‐cycle network are expressed in the following ways: where Ω is the volume of the whole compartment of the cell. There are 20 elementary processes for the current E2F network model, and the corresponding transition rates are proportional to the system size Ω. Figure 1 illustrates all the interactions between different species, and Table 1 describes the respective processes and transition rates. Here [S] represent the concentration of the external serum growth signal.

FIGURE 1

The schematic model of the Rb‐E2F network with the corresponding 20 reactions and time delays incorporated in the auto‐synthesis of E2F and in E2F‐dependent production of CycE‐Cdk2 complex

TABLE 1

Reaction kinetics, terms, and descriptions

Transition Number	Transition Processes	Transition Rates	Description
V 1	E2F→E2F+1	k_E([Myc]I_M+[Myc])([E2F]I_E+[E2F])Ω	E2F production modulated by cooperation between Myc and E2F.
V 2	E2F→E2F+1	k_b[Myc]I_M+[Myc]Ω	E2F production modulated by only Myc
V 3	E2F→E2F+1RP→RP+1RE→RE−1	k_p1[CD][RE]I_CD+[RE]Ω	E2F release from Rb‐E2F complex due to its phosphorylation mediated by the CycD‐Cdk4/6 complex
V 4	E2F→E2F+1RP→RP+1RE→RE−1	k_p2[CE][RE]I_CE+[RE]Ω	E2F release from Rb‐E2F complex due to its phosphorylation mediated by the CycE‐Cdk2 complex
V 5	E2F→E2F−1	d_E[E2F]Ω	E2F degradation
V 6	E2F→E2F−1R→R−1RE→RE+1	k_RE[R][E2F]Ω	Rb‐E2F complex formation
V 7	R→R+1	k_RΩ	Rb integrated production
V 8	R→R+1RP→RP−1	k_DP[RP]I_RP+[RP]Ω	Rb production due to the dephosphorylation of RP
V 9	R→R−1RP→RP+1	k_p1[CD][R]I_CD+[R]Ω	CycD‐Cdk4/6 complex phosphorylation of Rb
V 10	R→R−1RP→RP+1	k_p2[CE][R]I_CE+[R]Ω	CycE‐Cdk2 complex phosphorylation of Rb
V 11	R→R−1	d_RΩ	Rb degradation
V 12	RP→RP−1	d_RP[RP]Ω	RP degradation
V 13	RE→RE−1	d_RE[RE]Ω	RE degradation
V 14	Myc→Myc+1	k_M[S]K_S+[S]Ω	Myc production owing to an external serum stimulation
V 15	Myc→Myc−1	d_M[Myc]Ω	Myc degradation
V 16	CD→CD+1	k_CD[Myc]I_M+[Myc]Ω	Myc drives CycD‐Cdk4/6 complex activation and formation
V 17	CD→CD+1	k_CDS[S]K_S+[S]Ω	CycD‐Cdk4/6 complex activation and formation facilitated by the external serum signal (S)
V 18	CD→CD−1	d_CD[CD]Ω	CycD‐Cdk4/6 complex degradation
V 19	CE→CE+1	k_CE[E2F]I_E+[E2F]Ω	CycE‐Cdk2 complex activation and formation by E2F
V 20	CE→CE−1	d_CE[CE]Ω	CycE‐Cdk2 complex degradation

The direct stochastic simulation method is accurate, but it takes too much time when the system size is large. To overcome this challenge, based on the 'macro‐infinitesimal' timescale of the system, Gillespie found that CLEs [62] is instead a good approximation for studying the system's dynamics. It is, therefore, more convenient to qualitatively examine the influence of internal noise using CLEs. Then, according to the relationship between the concentration and molecular number, the CLEs corresponding to Yao's [60] macroscopic differential equations can be represented as:

In the above mesoscopic stochastic differential equations, internal noise intensity is proportional to Therefore, if the cell volume is big enough, that is, as Ω → ∞, internal noise will vanish in this thermodynamic limit.

A natural boundary at , , , , , and are used since the noise may drive the system variables to negative values. Artificial boundary conditions may affect the calculation accuracy [63]. Schnoerr et al. have shown that the chemical Fokker‐Planck equation is a good approximation of the CLEs even for a small number of molecules [64]. Fortunately, even in small volume Ω = 50, our calculation does not encounter the boundary.

There are time delays between any change in E2F's concentration and the response to that change in E2F synthesis and degradation processes. For example, using a dual E2F reporter system, Dong et al. suggest two delays in the G1 phase: an initial delay in the transcription of E2F and the delay between E2F's transcriptional dynamics and E2F's activity [65]. The later delay is due to the immediate sequestration of the newly synthesised E2F proteins by Rb. It not only affects the first delay but also regulates the whole cell cycle duration. Therefore, the delay between E2F's transcriptional dynamics and E2F's activity plays an essential role in the cell cycle process.

In this Rb‐E2F GRN, the reaction rates related to the concentration of E2F itself were V 1 and V 19. These two reaction rates are contained in two positive feedback loops. V 1 corresponds to the self‐regulated positive loop, and V 19 is related to the E2F‐CE‐E2F loop, as illustrated in Figure 1. Therefore, considering delays, equations V 1 and V 19 can be rewritten as follows:

To study the effects of noise and time delay on the expression level of [E2F], stochastic numerical simulations are needed, and the seven above mesoscopic equations are integrated by using a simple forward Euler algorithm. The time step is set at 0.001 h, and the initial transient dynamics were discarded in every simulation. Values of Ω, τ1 and τ2 were estimated and varied accordingly, while other parameter values in the model are adapted from reference [60]. The average expression level of [E2F], quasi‐potential function, and time evolutions are examined in four different cases: (i) when only noise is present; (ii) besides noise, when the time delay is only present in the production rate ; (iii) besides noise, when the time delay is only present in the production rate ; (iv) besides noise, when the time delay is considered in both production rates.

RESULTS AND DISCUSSIONS

Effects of internal noise on the system dynamics

A high expression level of E2F can push the cell cycle into progression [49]. Therefore, we can consider the high concentration state as an ‘on’ state, whereas the low concentration state is an ‘off’ state. Suppose the external stimulation [S] is chosen in the bistable region. The system is either in a high steady‐state or low steady‐state in the deterministic case, depending on the initial conditions. The time courses of [E2F] are plotted for different system volumes at [S] = 0.3, as shown in Figure 2, respectively. Here [S] = 0.3 is in the bistable region, but the E2F system is preferred to stay at a low stable state.

FIGURE 2

Time evolutions of [E2F] at different system volumes under a low serum signal [S] = 0.3. Noise is inversely proportional to cell volume, and at moderate volume size (Ω = 200), the transitions between high and low steady states are very obvious. When the volume size is big (Ω = 500), the noise level is too small to induce switching to the ‘on’ state. The [E2F] fluctuates around the ‘off’ state. The red line in the inset is the time evolution from the deterministic model

Under the disturbance of the internal noise, the state system is different. The smaller the system volume, the bigger the internal noise will be. At a small system volume Ω = 50, the noise is too large to induce approximately sustained random fluctuations with high amplitude of E2F. There is no obvious high‐low two‐state form of E2F. Even when the system volume increases to Ω = 100, oscillations in [E2F] begin to decrease in amplitude. When Ω = 200, the concentration of E2F is divided into evidently high and low bistable states, and each state has a long duration. There are very few transitions between the two states. However, when Ω = 500, the noise level is too small to induce switching to the ‘on’ state, so the system is completely maintained on the ‘off’ state. The red line in the inset is the time evolution from the deterministic model. The average value of [E2F] as a function of Ω is plotted in Figure 3, and it decreases monotonously with increasing of Ω.

FIGURE 3

E2F expression levels. The average value of [E2F] decreases with the increase of system volume under a low serum signal [S] = 0.3

These results here all together demonstrate that there exists a moderate noise intensity above which E2F is in the high‐frequency oscillation to elevate [E2F] (hence proliferation) and below which E2F switches off to quiescence. Our results rhyme with previous findings that if the volume is excessively increased, it results in reduced E2F expression levels (Figure 3), which consequently leads to quiescence, and if this is prolonged, it might result in senescence [66].

Effects of time delay in E2F itself

The cumulative time of E2F under the regulation of time delay was firstly explored in Figure 4. The serum signal [S] = 0.8 μM was chosen, since this concentration is significant enough to elevate the E2F to the ‘on’ state at a large cell volume that would otherwise cause quiescence when at [S] = 0.3 μM in Figure 2. To emphasize this delay's effect, we consider a very small internal noise intensity at Ω = 500. It is clearly shown in Figure 4 that increasing time delays postpone the rising time of E2F with all the trajectories raised to approximately the same maximum [E2F] level. It indicates how time delay regulates the E2F activation time to enable the cell exit quiescence, intimating on cellular quiescence depths. For deeper quiescence, longer activation times of E2F will be required to proliferate the cells from this state [43, 67]. As previously demonstrated, switching E2F from ‘off’ to ‘on’ positively correlates with transition from deep or shallow quiescence. Also, cells have a higher E2F switching threshold during deeper quiescence under the same growth conditions [53]. We can ratify that deeper quiescent cells also exhibit delayed E2F activation, hence delayed DNA replication is observed compared to cells in shallow quiescence. Delay times could facilitate the lock‐in strategy proposed in cancer therapy where cancer stem cells can be maintained in the G0 state to prevent tumour growth and metastasis [68]. But the final concentration of E2F changes not with the increase in delays. Therefore, time delay does not change the final E2F concentration but in a high steady‐state.

FIGURE 4

The time evolution of E2F for different time delays for Ω = 500 (low noise intensity) at constant serum signal concentration [S] = 0.8. Longer time delays stretch the rising time to shoot up the E2F expression to the ‘on’ state without altering the final E2F concentration level

Then, we sought to examine the influence of the delay time on the average value of [E2F] at a low steady‐state for different intrinsic noise levels. The serum growth signal chosen is low beyond the bistability region that is, [S] = 0.3 μM. Hence, as shown in Figure 2, without involvement of a stimulant such as some noise level, the expression level of E2F stays at the ‘off’ state. From Figure 5, it is clearly seen that the average value of [E2F] changes slightly with the increase of time delay at small (Ω = 50, Ω = 100) and big (Ω = 300) system volumes. These volumes correspond to high and low noise levels respectively. However, the trend changes remarkably with increasing delay time at medium noise level that is, system volume at Ω = 200. Firstly, E2F′s average value was slightly increased with time delay to a maximum by approximately 7 h. Then, it dramatically decreases to a minimum value by around 15 h. Finally, it increases again to almost a constant value (τ 1 > 20 h). These results concretise the fact that only at moderate noise intensity, time delay causes significant changes in the dynamics of average E2F concentrations.

FIGURE 5

The average value of [E2F] as a function of time delay under different system volumes/noise levels at a low serum signal concentration [S] = 0.3. At moderate noise intensity when Ω = 200, the average E2F changed remarkably with the increasing delay time

Next, we focus on the dynamic behaviour of the system when the volume is Ω = 200. Time courses of [E2F] are plotted for different time delay values, as shown in Figure 6. The residence times for which the system stays at high and low stable states are longer at τ 1 = 7 than at τ 1 = 0. Increasing the time delay to τ 1 = 10 does not increase the residence time further but initiates the appearance of pulse‐like oscillations. Low steady‐state and pulse oscillation state coexist at τ 1 = 15, and the amplitude of the oscillations begins to increase. From τ 1 = 20 to τ 1 = 50, the system is in an oscillation state. With the increase of time delay, the oscillation becomes more regular and the oscillation period becomes larger. Comparison of the state change from τ1 = 0 to τ 1 = 15 and from τ 1 = 20 to τ 1 = 50, are obvious differently. Therefore, we choose τ 1 = 15 as the critical delay value and discuss the behaviour of the system in the case of small‐time delay (≤15) and large delay (>15). From our results, we can demystify the existence of a certain delay time value below which E2F expression levels in off and on states are sustainably distinct and above which, the system goes into absolute oscillation. However, the two steady states coexist at this critical value.

FIGURE 6

Time evolutions of [E2F] at different time delay at system volume Ω = 200, at a constant low serum signal concentration [S] = 0.3. From τ 1 = 0 to τ 1 = 10, there are transitions between ‘on’ and ‘off’ states. At τ 1 = 15, oscillation state begins to emerge. From τ 1 = 20 to τ 1 = 50, typical oscillations only with changes in amplitude

Since the toggle between different states corresponds to different phenotypes, it is expedient to obtain relative depth or order of the metastable states through a quasi‐potential function U(x), where the variable x here refers to [E2F] [69]. To calculate the potential U, CLEs were run for 50,000 h to get a steady‐state histogram of [E2F]. The simulation time is long enough, so we get a smooth steady‐state probability distribution (P ). Then the quasi‐potential as per Kwon et.al approach [53] is defined as: Where P is the steady‐state distribution at the saddle point. It is clearly shown in Figure 7 that there is a typical peak between two troughs in the landscape, which corresponds to an unstable state of E2F, the saddle point. With such a function, it is possible to explicitly illuminate the probability and direction of transition between the stable states in a noisy or stochastic system [69, 70, 71]. There are two troughs in the E2F landscape, one corresponding to E2F 'off' state (left U ) and the other one corresponding to the 'on' state (alright U ). Quiescence is associated with E2F 'off' state, while the activation and cell cycle entry are associated with a high concentration level of E2F, that is, E2F 'on' state. Besides, there is a common peak in the landscape, which corresponds to an unstable state of E2F. This unstable state can be considered as the saddle point of the potential U . Then, the potential barrier for E2F activation from the ‘off’ state is given by ΔU = U −U . The potential barrier for E2F inactivation from the ‘on’ state is given by ΔU = U −U . The deeper the potential trough, or the higher the potential barrier, the harder the steady state transition [53]. Kwon et al., 2017 illustrated that the depth of quiescence regulated by Cdk inhibition rate affects cell cycle re‐entry and correlates with the relative change in quasi‐potential [53].

FIGURE 7

The system's quasi‐potential landscape at different delays for system volume Ω = 200, [S] = 0.3. (a) Small time delay cases: τ 1 = 0 (black line) has the deepest quiescent well (corresponding to ‘off’ state). Depth of other delay time does not differ so much in the ‘off’ state. For τ 1 = 10 (orange line) and τ 1 = 15 (blue line), the depth in the ‘on’ state distinctly differs. Deeper potential barriers in the E2F‐on states than the ‘off’ state for all time delays except for τ 1 = 15. (b) Big time delay cases: No palpable difference in depths between the steady states corresponding to different time delay values. Inset: enlargement of the ‘on’ troughs, the location of the ‘on’ trough shift to the high value of E2F from delay time τ 1 = 15 (blue line) to τ 1 = 30 (magenta line)

Here, we categorised delay times as small (i.e. τ 1 ≤ 15) or big (i.e. τ 1 > 15). In the case of small delays, from τ 1 = 0 (black line) to τ 1 = 4 (red line) and then to τ 1 = 7 (green line) in Figure 7(a), it can be seen that with the increase of time delay, the depth of low steady‐state potential trough becomes shallower, while the depth of high steady‐state potential trough remains unchanged. ΔU is smaller when τ 1 = 7, so the system is easier to be activated and thus increased E2F activation readiness towards cell proliferation. Therefore, the average value of E2F in Figure 5(c) reaches the maximum value at τ 1 = 7 h.

When the time delay approaches the critical delay, see τ 1 = 10 (orange line) and τ 1 = 15 (blue line) in Figure 7(a), the depths of the ‘off’ state are almost similar to that of τ 1 = 7. While the trough at the ‘on’ state is becoming shallow, with τ 1 = 15 reaching the most superficial state, it results in the decrease of the inactivation barrier potential, ΔU . It signifies an increase in the abilities of E2F switching back to the ‘off’ state. Therefore, the probability that E2F stays at the ‘on’ state decreases, and the average value of E2F reaches the lowest value in Figure 5(c). Besides, the depth of the two troughs is shallow and is very similar at τ 1 = 15. It indicates that switching between the ‘on’ and ‘off’ states is easy under this delay, and oscillations began to appear. Thus, the system dynamics change from longer residence times in stable states to oscillation with increased time delay.

In the case of big time delays, there is no obvious difference in depth of the troughs on either the 'off' or 'on' states with different delay times in Figure 7(b). The depth of the two troughs of all quasi‐potential curves is shallow and approximately the same. Therefore, switching E2F between 'on' and 'off' states is easy, and the whole system is in the oscillation state. But analysing the ‘on’ state at an enlarged scale in Figure 7(b) inset, it is evident that the location of the 'on' trough shifts to a higher value of E2F at τ 1 = 30 (magenta line) than at τ 1 = 15 (blue line). It explains why the average of E2F goes up from τ 1 = 15 to τ 1 = 30 in Figure 5 (green line with stars). All put together, the delay in E2F itself does not destroy the bistability of E2F but fine‐tune its expression level.

To qualitatively determine the value of delay at which the system is set into periodic oscillation, we calculated the power spectra of the time series for the different τ 1 values. The power spectra were computed from 500 averages of the power spectra for the time series of [E2F] with length 16,384 using a fast Fourier transform. It is shown in Figure 8(a) that at small τ 1 values (0‐10), no peak was obtained. Hence no oscillations are invoked. As the delay increased to τ 1 = 15 (blue line in Figure 8(b)), a small peak emerges, which characterises the onset of periodic oscillations. The height of the peak grows with an increasing value of τ 1 in Figure 8(b). Actually, from the power spectrum, τ 1 = 50 (brown line) forms the highest peak. In addition, by comparing the power spectra at τ 1 = 20 (magenta line) with τ 1 = 50 (brown line) in Figure 8(b), it is found that the frequency corresponding to the highest peak is decreased.

FIGURE 8

Power spectra of [E2F] at different delays for system volume Ω = 200, [S] = 0.3. Small τ 1 values (0‐10) no peak observed. Big τ 1 (15‐50), peak starts to appear and its value increases steadily with delay value, a manifestation of increase in periodic oscillation

Effects of time delay in CE

Next, we sought to examine the influences of the delay time on the average value of [E2F]. Similarly, the serum growth signal is chosen low beyond the bistability region at [S] = 0.3 μM. From Figure 9, it can be clearly seen that the average value of [E2F] decreases with the increase in time delay at small Ω = 50 (black line with rectangles) and big Ω = 300 (blue line with triangles) system volumes. Moreover, when Ω = 300, the average value of [E2F] is very close to 0, indicating that it has been in the low steady‐state and cannot transit to the high steady‐state. It increases slightly with the increase of delay time at Ω = 100 (red line with open circles). Interestingly, it increases remarkably at medium system volume, see at Ω = 200 (green line with stars). From to , the average value of [E2F] increased rapidly with the delay, and then, when , it remained in a high expression level almost unchanged. Medium system volume Ω = 200 corresponds to moderate intensity of internal noise. Therefore, there is a synergistic effect between noise and delay, making the system stay in the ‘on’ state.

FIGURE 9

The average value of [E2F] as a function of time delay under different system volumes, [S] = 0.3. At moderate volume, Ω = 200, the average E2F increased markedly with delay time

We focus on the dynamic behaviour of the system when the volume is Ω = 200. Time courses of [E2F] are plotted for different time delays, as shown in Figure 10. In the comparison of time evolutions, at τ 2 = 2 and τ 2 = 0, it can be seen that the addition of a small delay can significantly increase the residence times. It is consistent with the conclusion of reference [72] that the increasing delay can stabilise the bistable gene network. Increased time delay to τ 2 = 5, to τ 2 = 10, and then to τ 2 = 20 does not increase the residence time further but destabilised the bistability. The system is preferred to stay at a high steady‐state. A big time delay τ 2 = 50, approximate regular transitions between ‘on’ and ‘off’ states but not oscillations as in the τ 1 case appear. However, it is obvious that the high steady‐state's residence time is longer than for the low steady‐state.

FIGURE 10

Time evolutions of [E2F] at different τ 2 at system volume Ω = 200, [S] = 0.3. From τ2 = 0 to τ2 = 2, the transitions between ‘on’ and ‘off’ states are decreased, that is, increasing delay increases the residence times in ‘on’ and ‘off’ stable states. From τ2 = 5 to τ2 = 20, the value of [E2F] fluctuates randomly near the high steady state. At τ2  = 50, it begins to transit a little regularly between high and low steady‐states

The system's quasi‐potential landscapes under the regulation of different time delays for four system volume are shown in Figure 10. Though at small (Ω = 50) and big (Ω = 300) system volumes, increasing time delay will reduce the average [E2F] similarly (as shown in Figure 9). The shapes of the quasi‐potential function U(x) change differently with the increasing delays in Figure 11(a), (b), and (d). When τ 2 increases from 0 to 5 in Figure 11(a), the trough for the low steady‐state disappears. It means that an increase in τ 2 causes the system to lose bistability, and E2F is always at a high expression level. While for the big system cases in Figure 11(d), an increase in τ 2 causes the system to lose bistability, and E2F is always at a low expression level.

FIGURE 11

The system's quasi‐potential landscapes at different time delays for four system volumes, [S] = 0.3. (a) Small system volume cases: Ω = 50, increasing the delay destroys the bistability of the system, and only the high steady‐state exists at (green line). (b) system volume: Ω = 100, similar to small volume cases, and only the high steady‐state exists at (green line). (c) medium volume cases: Ω = 200, increasing the delay to (red line) stabilised the bistability of the system (Inset: enlargement of (black line) and ). However, big delays stabilised the high steady‐state, see (green line), (blue line), (claret line) and (orange line). (d) Big system volume cases: Ω = 300, increasing the delay destroys the bistability of the system, and only the low steady‐state exists at (red line)

The shape of U(x) changes complicatedly as time delay changes at medium system volume Ω = 200 in Figure 11(c). Slightly increasing time delay to τ 2 = 2 (red line in the inset of Figure 11(c)) causes the trough depth to deepen and widen. It implies that an increase in time delay to τ 2 = 2 makes the system's transition between the potential troughs more difficult. Time delay stabilised the bistability of the E2F gene regulate network [72]. Increasing the delay further to τ 2 = 5, to τ 2 = 10, to τ 2 = 20, and then to τ 2 = 50 facilitates the high steady‐state. Consequently, the average value of [E2F] is maintained at a high expression level in Figure 9 (green line with stars).

Similarly, we calculated the power spectra of the time series for the different τ 2 values. It is clearly shown in Figure 12 that for all the different τ 2 values, the power spectra peaks are either very minimal or nothing at all. Therefore, increasing this type of delay will not induce oscillation. Anyway, as observed in Figure 10, no oscillations were created by this kind of delay.

FIGURE 12

Power spectra of [E2F] at different delays for system volume Ω = 200, [S] = 0.3. With the same scale along the y‐axis as in Figure 8. Generally, almost no peak or very minimal peaks formed at all τ 2 values, which implies no or marginal periodic oscillations occurs

Effects of time delay in E2F and CE

Next, we consider the two delays simultaneously, that is, consider the delay between E2F transcriptional dynamics and activity dynamics. The average value of [E2F] versus delay time τ 1&τ 2 under different system volumes are shown in Figure 13. In the case of small volume Ω = 50, Ω = 100 and large volume Ω = 300, the variation of average [E2F] with two delays is similar to that of τ 2 considered alone. However, at moderate noise level, that is, Ω = 200, the curve is increased first and then decreased. The change of average [E2F] with delay is similar to resonance. Moreover, the range of delays (3‐13 h), which induced the high expression level of E2F, coincides with the experimental range by Dong et al., (for the transcriptional delay is in the range 6‐20 h and for the delay between E2F's transcriptional dynamics and E2F's activity is in the range of 0‐10 h) [65]. Therefore, this range of delays is physiologically meaningful.

FIGURE 13

The average value of [E2F] as a function of two delays under different system volumes, [S] = 0.3. At moderate system volume, Ω = 200, the average E2F increased remarkably with delay time

Time courses of [E2F] are plotted for different time delays at system volume Ω = 200, as shown in Figure 14. The residence times in ‘on’ and ‘off’ states are extended by adding small delays, see τ 1&τ 2 = 0 and τ 1&τ 2 = 2. Increased time delays to τ 1&τ 2 = 10 do not increase the residence time further but destroyed the bistability. The system is preferred to stay at a high steady‐state. These [E2F] changes with the increasing delay times τ 1&τ 2 were similar to those of τ 2 alone. However, when τ 1&τ 2 = 15, the situation is different. There are fast and nearly periodic transitions between high and low steady‐states that begin to appear. These are transitions but not oscillation as in the case where τ 1 regulates alone in Figure 6 for τ 1 = 15. Increased time delay to τ 1&τ 2 = 20 and to τ 1&τ 2 = 50 increases the period of the transiting, and the residence time of high steady‐state is increased.

FIGURE 14

Time evolutions of [E2F] at different τ 1 & τ 2 at system volume Ω = 200, [S] = 0.3. From τ 1 & τ 2 = 0 to τ 1 & τ 2 = 2, the transitions between ‘on’ and ‘off’ states are decreased, that is, increasing delay increases the residence times in ‘on’ and ‘off’ stable states. At τ 1 & τ 2 = 10, [E2F] are occasionally transmitted to low steady‐state. At τ 1 & τ 2 = 15, [E2F] are transited frequently between ‘on’ and ‘off’ states. At τ 1 & τ 2 = 20, it is distinct that the frequency of transiting is decreased. At τ 1 & τ 2 = 50, the residence time of high steady‐state is increased, and the frequency of transiting is decreased further

Figure 15 depicts the quasi‐potential U(x) for four system volumes and different delay times. Small (Ω = 50, 100) and big (Ω = 300) system volumes, and the variation of quasi‐potential U with the increase of delay is similar to that of τ 2 alone. At medium system volume Ω = 200, increasing the delays to τ 1 & τ 2 = 2(red line) deepens both troughs. Delays stabilised the bistability of the system. Increasing the delay to τ 1 & τ 2 = 5 h (Inset: green line) and τ 1 & τ 2 = 10 h (blue line) makes the trough of high steady‐state more profound. It means that increasing the delay is beneficial to a high steady‐state. However, further increasing delay times to τ 1 & τ 2 = 15 h (magenta line) and τ 1 & τ 2 = 50 h (orange line) makes the two potential troughs change significantly, the low steady‐state potential trough becomes deeper, and the high steady‐state potential trough becomes shallower. At τ 1 & τ 2 = 50 h, the depth of the two troughs is very close. It implies that the transitions between high and low steady‐states become easier. This corresponds to Figure 13, when τ 1 & τ 2 is greater than 15 h, the average value of E2F begins to decrease.

FIGURE 15

The system's quasi‐potential landscapes at different time delays τ 1 & τ 2 for four system volumes, [S] = 0.3. (a) and (b) Small system volume cases: increasing the delay destroys the bistability of the system, and only the high steady‐state exists at τ 1 & τ 2 = 5hr (green line). (c) Medium volume cases: Ω = 200, increasing the delay to τ 1 & τ 2 = 2hr (red line) stabilised the bistability of the system (Inset: enlargement of τ 1 & τ 2 = 0 (black line) and ). Big delays stabilised the high steady‐state, see (Inset: green line), τ 1 & τ 2 = 10 hr (blue line). However, further increasing delay times to (magenta line) and to τ 1 & τ 2 = 50 hr (orange line) induce transitions between ‘on’ and ‘off’ steady‐states. (d) Big system volume cases: Ω = 300, increasing the delay destroys the bistability of the system, and only the low steady‐state exists at τ 1 & τ 2 = 2hr (red line)

The power spectra of the time series for the different τ 1 & τ 2 values are shown in Figure 16. It was evident that no peaks were created for small τ 1 & τ 2 values (0‐10), which explains the absence of oscillations for such small delays seen in Figure 14. However, an increase in delay values say at τ 1 & τ 2 = 15(blue) is enough to cause periodic oscillations. A delay value of τ 1 & τ 2 = 20 (magenta) shown in Figure 16 produced the highest peak, which is very similar to periodic oscillation. For the biggest delay τ 1 & τ 2 = 50 (brown), the peak value of the power spectra decreases, and the periodicity decreases. At the same time, the frequency corresponding to the peak is decreasing, as observed in Figure 14.

FIGURE 16

Power spectra of [E2F] at different delays for system volume Ω = 200, [S] = 0.3. For small delay values (τ 1 & τ 2 = 0‐10), no peaks of the power spectra were observed. For the case of big delays, peaks form in the spectra such that the highest is at τ 1 & τ 2 = 20 (magenta)

CONCLUSIONS AND DISCUSSION

Noise and delays can play pivotal roles in regulating the functions of genes. Based on the stochastic Rb‐E2F model, the effects of internal noise and time delay on the E2F's dynamics related cell cycle system have been investigated. It was found that the internal noise will induce the switching between 'on' and 'off' states when the external stimulation is set in the bistable region. The bigger the noise, the more frequent are the switches, and the higher will be the average [E2F].

At the moderate intensity of internal noise, when Ω = 200, adding different time delays in the system can significantly change the average value of [E2F]. Especially when the delay is added to the Cyclin E synthesis rate, increasing the delay to can greatly increase the average concentration value of E2F. When the delay is considered in both E2F‐related positive feedback loops, within a specific range of delay (3‐13), the average expression of E2F was significantly increased, also this range is in the scope given experimentally by Dong et al. [65]. By analysing the quasi‐potential curves at different delay times, it was found that whether it is under the regulation of or under the regulation of τ 1 & τ 2 simultaneously, increasing the time delays will experience: small delay which stabilises the bistable system; the medium delay that favours high steady‐state conditions, thereby making the system fluctuate near the high steady‐state; large delay induces approximately periodic transitions between high and low steady‐state.

Our exploration of the effects of noise and delay is based on Yao's model [60], which uses the Michaelis–Menten (MM) rate law to describe the key interactions of E2F during the G1/S cell cycle transition. Results from Yao's model are qualitatively in good agreement with those in their experimental work. Besides, Yao et al. also give a minimum network structure with which the system bistability can be maintained [73]. Also, some other study points out the critical parameters that control the activation threshold of this kind of network switch [53]. However, Segel and Slemrod [74] spotlighted that only when the enzyme concentration is low enough to the substrate, quasi‐steady state approximation (QSSA) [75] can be used to simplify the full mass‐action rate laws to the MM rate laws. Recent studies [76, 77, 78, 79] have shown that the concentrations of transcription factors (TFs), kinases, phosphatases, inhibitors, etc., in protein interaction networks are usually comparable. Therefore, it is risky to use an MM rate law to describe the interactions between TFs and gene promoter regions and describe the interaction between ligands and receptors. Instead, Kim et al., 2020 [80] suggested that the total quasi‐steady state approximation (tQSSA) model is accurate for any combination of substrate and enzyme concentrations. Future work is best to confirm the experimental basis of MM rate law's validity condition in the E2F‐Rb cell cycle network.

Then, by converting the protein concentration to its number of molecules, it gives a good approximation of Gillespie's stochastic simulation method, so CLEs are adopted in our paper. However, Kim et al. [81] show that the accuracy of such stochastic QSSA is determined by timescale separation and by the sensitivity of the QSSA solution at the same time. Furthermore, by using both theory and simulations, James and Roman show that Hill functions in describing transcriptional regulation are only valid in a fast promoter switching condition [82]. Stochastic (tQSSA) is more accurate than the stochastic QSSA. It is difficult to find an analytical solution of the tQSSA of the complex E2F‐Rb network, Kim and partners' [80] excellent work on a single substrate enzyme‐catalysed reaction system hints that in the future, we can make use of the tQSSA to get the simplified model from the full detailed model based on mass‐action kinetics.

In general, the previous work showed that cell quiescence depth is controlled by the activation threshold of an Rb‐E2F gene switch [53], that is, the barrier of activation in the quasi‐potential landscape. Different components (parameters) in the network have different effects on the threshold of the switch. However, in our work, we find that the time delay also adjusts the troughs' depth remarkably. Although the delay in our simulation is artificially added, an experiment by Dong et al. pointed out that there is a delay between E2F transcriptional and activity dynamics [65]. Our work may be helpful to explain some future experimental results intricacies related to the cell cycle.