Funneled potential and flux landscapes dictate the stabilities of both the states and the flow: Fission yeast cell cycle

Using fission yeast cell cycle as an example, we uncovered that the non-equilibrium network dynamics and global properties are determined by two essential features: the potential landscape and the flux landscape. These two landscapes can be quantified through the decomposition of the dynamics into the detailed balance preserving part and detailed balance breaking non-equilibrium part. While the funneled potential landscape is often crucial for the stability of the single attractor networks, we have uncovered that the funneled flux landscape is crucial for the emergence and maintenance of the stable limit cycle oscillation flow. This provides a new interpretation of the origin for the limit cycle oscillations: There are many cycles and loops existed flowing through the state space and forming the flux landscapes, each cycle with a probability flux going through the loop. The limit cycle emerges when a loop stands out and carries significantly more probability flux than other loops. We explore how robustness ratio (RR) as the gap or steepness versus averaged variations or roughness of the landscape, quantifying the degrees of the funneling of the underlying potential and flux landscapes. We state that these two landscapes complement each other with one crucial for stabilities of states on the cycle and the other crucial for the stability of the flow along the cycle. The flux is directly related to the speed of the cell cycle. This allows us to identify the key factors and structure elements of the networks in determining the stability, speed and robustness of the fission yeast cell cycle oscillations. We see that the non-equilibriumness characterized by the degree of detailed balance breaking from the energy pump quantified by the flux is the cause of the energy dissipation for initiating and sustaining the replications essential for the origin and evolution of life. Regulating the cell cycle speed is crucial for designing the prevention and curing strategy of cancer.


Introduction
The global stability and robustness are crucial for maintaining the function. They are also important for uncovering underlying mechanisms of the networks. [1][2][3][4][5][6][7] However, it is difficult to quantify them for dynamic systems and networks. This presents a challenge for the dynamical systems and the field of systems biology. [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23] In equilibrium systems, the global nature of the system is characterized by the underlying equilibrium potential landscape U which is directly linked to the equilibrium probability through the Boltzmann distribution law P * exp(−βU). The local dynamics is determined by the gradient of the equilibrium potential landscape. However, most dynamical systems do not typically have a gradient potential as in the equilibrium case. They are open systems usually not in isolations. Global natures of such systems are hard to address. In addition, for mesoscopic systems, the intrinsic fluctuations can also be significant. Under stochastic fluctuations, instead of following the dynamical trajectories which are stochastic and unpredictable, the evolution of the probabilistic distributions should be followed, which is inherently global as well as predictable due to its intrinsic linearity. The probabilistic evolution is governed by the master equations for discrete state space (more general) and Fokker-Planck equations for continuous state space.
In this work, we study the more general stochastic dynamics in discrete space of the nonequilibrium networks (Markov chains) governed by the probabilistic master equation. We found the network dynamics and global properties are determined by two features: the potential landscape and the probability flux landscape. While potential landscape quantifies the probabilities of different states forming hills and valleys, the probability flux landscape is composed of many flux loops flowing in state space. Therefore, statistics of the flux loops becomes important. These two landscapes can be quantitatively constructed through the decomposition of the dynamics into the detailed balance part and non-detailed balance part.
We found that while funneled landscape is crucial for the stability of the single attractor networks and stability of oscillation states. The funneled flux landscape is crucial for maintaining the stable limit cycle oscillation flow. The stability and the robustness of the networks can be quantified through a dimensionless ratio of the gap or steepness versus the averaged variations or roughness of the landscape (which measures the degree of funnel, we termed as robustness ratio RR), and explored under the changes of the network topologies and stochastic fluctuations.
This flux landscape picture provides a new interpretation of the origin of the limit cycle oscillations. The global oscillation only emerges when one specific loop stands out and carries much more probability flux, and therefore becomes more probable than the rest of the others.
We specifically studied the fission yeast cell cycle as an example to illustrate the idea. We found the flux landscape of the fission yeast cell cycle oscillations is funneled, which guarantees its stability and the robustness of the oscillation flows. The global stability is quantified by the robustness ratio RR of the funneled flux landscape and the robustness is quantified by RR against the changes in topology of the network (wirings) and stochastic fluctuations. The flux is directly related to the speed of the cell cycle. The landscape analysis here allows us to identify the key factors and structure elements of the networks in determining the global stability, speed, and robustness of the fission yeast cell cycle oscillations. We see that the non-equilibriumness characterized by the degree of detailed balance breaking from the energy pump quantified by the flux is the cause of the energy dissipation for initiating and sustaining the replications essential for the origin and evolution of life. The cell cycle speed is a hallmark of cancer. Regulating the cell cycle speed thus provides a possible strategy for preventing and curing strategy against cancer.

Model and methodology
The stochastic boolean model of fission yeast cell cycle We follow a boolean network model mainly built for fission yeast in [32]. As illustrated in the Table 1, the cell cycle of fission yeast is divided into several phases: START phase ! G1 phase ! S phase ! G2 phase ! M phase ! G1 phase. The network wiring diagram shown in Fig 1 consists of one check point and 9 gene nodes, that is CS (Cell Size), SK, Cdc2/Cdc13, Ste9, Rum1, Slp1, Cdc2/Cdc13 Ã , Wee1/Mik1, Cdc25, and PP. The check point of this cell cycle network is named as Cell Size (CS) as it mainly works as an indicator of mass of the cell [32,33]. In the global state space, there are 2 10 states. Each state is the combination of the "on" (s i = 1) and "off" (s i = 0) states, which can be represented by a state vector S = {s 1 , s 2 , s 3 , . . ., s 10 }. With this representation, it can form a state space of a complex gene regulation interaction networks [34]. In general, one can use the boolean dynamics model to explore the coarse grained dynamics with the information of the wiring topology of the networks [33], and one can find that there are robust biological pathway and a global robust G1 state inside the state space.
However, from the view of cell cycle, the stationary state G1 which converges the biological pathway is a temporal fixed point. To perform the function of cell cycle, one needs a positive feedback to push the G1 state to pass through the checkpoint of cell size to activate the node SK, and then go through the dynamic oscillatory trajectory, and finally get back to the G1 state. Therefore we add a kinetic excitation once the global G1 state is reached. For biological meaning, if there is no nutrition supply constantly pumping into the system, the cell cycle will stay in the stationary G1 state (so-called G0 state). While given enough nutrition supply, the cell staying at G1 gains a large transition probability to reach the cell size checkpoint and then activate the node SK to drive the cycle. Such a pumping or driving force has been explicitly discussed in a model system for limit cycle oscillation as the chemical energy supply [35]. There, the chemical energy supply, which may be generated from ATP or GTP, acted as a chemical pump or battery for initiating and maintaining the cycle. In dynamics, the chemical energy supply is in terms of flux driving the cycle flow.
This cycling conception can be emerged well in a stochastic Boolean network model. We have developed an approach to probabilistically describe the network dynamics as follow: Firstly, we calculate transition probabilities T(s i (t 0 )|S(t)) of gene node i jumping from state s i (t) to state s i (t 0 ), where t and t 0 are two close neighborhood time moments. Based on the Markovian process theory, all the cases of transition probability T(s i (t 0 )|S(t)) are shown in the black box in Table 2. [24,28,[36][37][38][39][40].
In Table 2, I ¼ P 10 j¼1 a ij s j ðtÞ is defined as the total input of the ith gene node at time t, which is the summation of each of the interaction strength a ij from jth gene node to ith gene node.
Just considering the simplifications of the interactions between two nodes (activated (+1) or repressed (-1)), one can obtain that: when the total input I > 0, gene node i has high probability to stay at state s i (t 0 ) = 1; when the total input I < 0, gene node i tends to be repressed at state s i (t 0 ) = 0. It is similar with the behavior of switching function, that is why we define the transition probability as T(s i (t 0 )|S(t)) = 1/2 ± 1/2tanh(μI). In this way, we can clearly see that the transition probability is mainly determined by the input. While in the case of I = 0, the ith gene node mainly stays at present state, just with a little probability of c to change the state due to the background production and self-degradation. There is only one exception for the case I = 0, that is when the cell cycle stays at the G1/G0 state, where we add a cycling activation strength of γ to activate the check point of CS, which corresponds to the condition with enough chemical supply as mentioned above.
The biological meaning of the three parameters μ, c and γ in Table 2 can be illustrated as follows: μ can be considered as a mean transition strength from the input to output of a gene or protein node, which is also related to the inverse of the fluctuation or noise strength.
c is a parameter to quantify the effect of perturbation from the background production or degradation with not input (I = 0). For example, gene node i without any input interactions has a probability of c to change the present state, with 0 ! 1 due to the background production, while 1 ! 0 due to the background self-degradation (direct degradation and growth dilution).    γ represents the probability of the positive feedback or stimulation from the stationary G1 state to the check point of cell size, that is from G1 phase to START phase. The cell cycle stays at G1/G0 at present implies the fact that there is no input interaction (I = 0) for all the node (so-called global stable steady state). When the nutrition supply is enough, the chemical energy pump will push the stationary G1/G0 state to activate the check point node CS (i = 1) and reach the state of START phase, with a probability of γ (with a probability of (1 − γ) to stay at G0).
Secondly, we figure out that the transition probability between two neighbor states in the Markov time series chain can be written as the product of each node transition [23,27] TðSðt 0 Þ ¼ fs 1 ðt 0 Þ; s 2 ðt 0 Þ; :::; s 10 ðt 0 ÞgjSðtÞÞ ¼ Finally, we obtain the evolution equation to guide the probabilistic dynamics, which is so called master equation [23,41,42]: where P i represents the probability of state i, and the transition probability T ij represents the transition probability from state i to state j. Here, we use T ij as a discrete transition probability to fit the master equation, whose meaning is equal to the continuous transition probability expression of T(S(t 0 )|S(t)). The physical meaning of the master equation is the conservation law of probability: the local change of the probability of a particular state i in time is equal to the probability flow (flux) from the other states to this state i given by ∑ j T ji P j subtracting the probability flow (flux) from the state i to other states ∑ j T ij P i . By solving the 2 10 = 1024 master equations numerically, we obtain both the time-dependent evolution and the steady-state probability of each state in the global state space.

Decompositions of boolean dynamics and probability flux loops
For steady state, we set the left term of the master eq (2) to zero, that is dP i dt ¼ 0, then we obtain the numerical steady state solution P ss i , which is the long time limit. Given the steady solution, we can define the steady state flux between state i and j as: F ij ss ¼ À T ij P ss i þ T ji P ss j , If for any i, j pair, F ij ss = 0, this Markov chain is detail-balance preserved, and the steady state of the system becomes the equilibrium state (without net local flux), since dP i /dt = ∑ j F ij ss = 0. However, in general the steady state probability can be obtained, but it does not have to satisfy the detailed balance condition(F ij ss 6 ¼ 0). In other words, the net flux does not have to be zero. The system is then in non-equilibrium steady state. Although the steady-state distribution is fixed and does not change in time, there can be an internal probability flow among states. Table 2. Transition probability T(s i (t 0 )|S(t)) of gene node i.

Input Output
s i (t 0 ) = 1 s i (t 0 ) = 0 In order to study the non-equilibrium steady states and characterize the global properties, one can separate the dynamical process into two parts, a detailed balance part and a pure irreversible non-detailed balance flux part by decomposing the transition probability matrix M [29]. The master equation can be rewritten as dP=dt ¼ M T P, where P is the vector of probability of all the discrete states, M is the transition probability matrix (or rate matrix) with M ij = T ij , i 6 ¼ j and M ii = (−1)∑ j T ij . We define a matrix C such that the ith row and jth column of it is given as C ij ¼ maxfT ij P ss i À T ji P ss j ; 0g=P ss i ; i 6 ¼ j and C ii = (−1)∑ j C ij , and matrix D whose ith row and jth column is given as D ij ¼ minfT ij P ss i ; T ji P ss j g=P ss i ; i 6 ¼ j and D ii = (−1)∑ j D ij . It follows that M ¼ C þ D and D T P ss ¼ 0. Since M T P ss ¼ ðC þ DÞ T P ss ¼ 0, C T P ss ¼ 0. By separating the transition probability matrix this way, two Markov processes are obtained [29]. The probability transition matrix (or rate matrix) M for characterizing the dynamics can be decomposed into two terms: C and D. Both C and D have the same steady-state(stationary) probability distribution, and one of the processes D satisfies detailed balance(D ij P i ss = D ji P j ss ), while the other C is non-detailed balanced and irreversible (if C ij P i ss > 0, C ji P j ss = 0). In this way, the dynamics is decomposed to detailed balance preserving part and detailed balanced breaking part. The non-equilibrium irreversible part can be termed as the circulation or flux part, since it can be further decomposed into flux circles or loops with a flux value on each cycle [29]. The prove of the circulation also provides a way to obtain all the circles and their corresponding flux values for the dynamic part. By definition of flux, we have We can keep on doing this, until a repeat is found: If J (1) 6 ¼ 0, repeat what we did above to find another cycle as well as its flux value, then subtract those fluxes from J (1) to get J (2) . Since the number of non-zero elements in J (i) is at least one less than that in J (i−1) , there exits an integer N such that J (N+1) = 0 (all the elements of the matrix J are zero). Therefore for the non-detailed balanced part of the dynamics, the flux can be decomposed to finite number of circles or closed loops, each with a flux value, J ¼ P i¼N i¼1 J ðiÞ [29]. Therefore, the decomposition and associated flux statistics can be directly carried out at the master equation level.
On the other hand, one can also find another way for decomposition and associated statistics of the fluxes from the stochastic boolean trajectories. One can follow the trajectory from one state S(t) at moment t to another state S(t 0 ) at the next moment t 0 . If one finds the same state at different moment, that is S(t@) = S(t), then all the states between these two same states can be defined as one loop. Then one should remove this loop from the trajectory and repeat the steps above to obtain all the loops. By making statistical analysis on all the flux loops, one can calculate the flux landscape using the formula below, U flux = −lnP flux , where [26][27][28][29]43] Therefore we have two quantitative features to characterize the system, one is the steady state probability and the other is the non-zero steady state probability flux which can be further decomposed into loops. The steady state probability obeys the evolution equation of the transition probability matrix (or rate matrix) at long time limit characterizing only the detailed balance part. The detailed balance condition allows one to identify the path independent probability measures [30]. This naturally leads to the potential. We can see how both potential and flux landscape influence the dynamics and stability of the system through an example on fission yeast cell cycle.
Entropy production rate and dissipation in the fission yeast boolean network As known, when an open system is under long time evolution, it can reach non-equilibrium steady state (NESS) [6, 7, 19-23, 31, 44]. The local steady state flux F ij ss ¼ À T ij P ss i þ T ji P ss j is not necessarily equal to zero (no detailed balance). In this condition we can define a generalized force referring to the generalized chemical potential (from j to i) A ji ¼ lnð [8, 10, 19-23, 31, 44]. There is a mapping between the cellular networks and electric circuits. The flux F ij corresponds to current I and chemical potential A ij corresponds to voltage V. The non-equilibrium cell network dissipates energy just as the electric circuits.
In the steady state, the heat loss rate is related to the entropy production rate. The entropy production or dissipation characterizes "time irreversibility" and provides a lower bound for the actual heat loss in Boolean network [8,10,[19][20][21][22][23]44]. The total entropy change is equal to the part from the system or source plus the part from the bath or sink (dissipation). Since in steady state the entropy change of the system is equal to zero, thus the total entropy change (source) is equal to the entropy change of the sink (dissipation). The total entropy change (source) = ∑F ij A ij is the entropy production and the sink term is dissipation. Therefore in steady state, knowing the entropy production, we know the dissipation quantitatively. The entropy S from the system part is defined as S = −∑ i P i lnP i and entropy production rate dS tot dt is given by: Entropy production is correlated with flux. When the steady state flux is zero, the entropy production or dissipation at steady state is zero. When the flux increases, the entropy production typically increases. Therefore, the entropy production or dissipation can also serve as a quantitative measure of how far away the system is from the equilibrium, or in other words, the degree of the detailed balance breaking.

Results and discussions Potential landscape of fission yeast cell cycle
Quantifying the potential landscape. To explore the global quantification of dynamical systems of this fission yeast networks, we define the underlying potential landscape U from the steady state probability U = −ln(P ss ). To visualize the potential landscape U over the 2 10 states space, we firstly draw the landscape spectrum, which is shown in Fig 2. Then we define the Robustness Ratio (RR) as the ratio of potential energy gap between the lowest potential energy value and the average potential energy, with the average variations measured by the standard deviation of potentials.
. The gap represents the separation between the lowest potential state with the rest of the decoys. This definition of RR holds for a single attractor. When we deal with cycle oscillations as fission yeast, we are interested in the stability of not only one state but the whole oscillating cycle states. Therefore we extend our definition to include all the states (10 states in our example of fission yeast cell cycle) on the oscillation path as the "native" states. So, for the cell cycle, the gap quantifies the bias or slope of the underlying landscape towards the native (potentials of all the states on the oscillation paths basin of attraction) while the average variations measure the fluctuations or roughness of the underlying landscape. In this way, we can quantify the topography of underlying potential landscape. Fig 2 was shown with the parameters μ = 5, c = 0.001 and γ = 60% for the fission yeast cell cycle model (Table 2). In the potential landscape spectrum there are 10 states representing the biological cell cycle phases (G1, S, G2, M) which are marked as green lines. We notice that they all settle at the bottom. This means the 10 states are sitting on the most stable cell cycle pathway. We then calculate RR = 2.29813. RR is significantly larger than 1, which indicates the 10 states along the biological cell cycle is separated from the others. So we state that this can be the reason that the cell cycle states are stable since the states of the biological path are all at low potentials and high probabilities sufficiently separated from others.
Furthermore, we draw all those states (2 10 = 1024) on a 2-dimensional surface by minimizing the distances between two states which have the strongest connections, and using the underlying potential energy to be the vertical axis. The color represents the potential level of each state on both surface and the bottom contour map. We obtain Fig 3, in which we notice that the potential landscape shows a distinct topology with a Mexican hat like shape. This gives us a visualized 3D picture for the potential landscape. The valleys of this landscape correspond to exactly the biological cell cycle with low potentials more easily seen at the bottom contour map.
In our earlier study of budding yeast system without explicitly putting in the excitations from G1 ground state to the start signal [22,23], we see quite different dynamics and landscape. There the potential landscape has a funneled shape. The system has one dominant basin of attraction pointing towards G1. The model can explain the dynamical process of yeast cell cycle once the start signal kicks off. Since the end state is always the G1 state (bottom of the funnel or basin of attraction), it does not contain the cycle part. The excitation from G1 state here giving the cycle is triggered by the nutrition supply or energy pump [17][18][19]35].
With the excitation explicitly implemented in the fission yeast cell cycle, we can see in Fig 3 the potential landscape changes the shape from single attractor funnel to Mexican hat shape which can guarantee the stabilities of the oscillation states. If we see Mexican hat landscape a global quantification of a closed loop line attractor, an effective "funneled" potential landscape towards the oscillation path emerges. The degree of the funneless is quantified by RR. The "funneled" potential landscape towards the oscillation guarantees the global stabilities of the states on oscillation paths. However, this potential landscape can not guarantee the stable directional flows for the oscillations. The state transitions or switching along the biological path from the energy pump provided by the nutrition supply directs the oscillation cycle of the fission yeast cycle (G1 ! S ! G2 ! M ! G1). This is globally manifested by the underlying flux landscape, as we will show later in this study.
Robustness of potential landscape against changes in sharpness of response, self degradation, and stimulations. To further explore the robustness of the networks with internal and external perturbations, we calculate the Robustness Ratio (RR), probability of the cell cycle path, entropy production against the variations of different parameters.
In Fig 4(a), we have shown the stationary probabilities for both G1 state and cell cycle path by fixing c = 0.001, γ = 60% and changing μ. As we mentioned above, the parameter μ refers to the mean transition rate of gene or protein switching and can be considered as the inverse of the noise strength or environment temperature [23,27]. So as the transition rate μ increases or the fluctuation decreases, the whole biological cell cycle, as well as the G1 state, becomes more stable monotonously. When comparing the biological cell cycle with the G1 phase by changing μ, the biological cell cycle is more statable than the G1 phase with higher probabilities.
The changing of RR while switching the parameter μ is shown in Fig 4(b). As we see at first the Robustness Ratio increases when μ increases or the fluctuation decreases. This means a large transition from input to output or smaller fluctuations makes a more robust network. Notice that the RR goes down to reach a certain value when μ is large enough. This result is similar to the behavior shown in the budding yeast cell cycle without the excitation from G1 ground state [23]. We have stated that there are many traps beyond the G1 state and circle path in the state space, which have low potentials and high weights at "low temperatures" (small fluctuations or high transition rate μ). The presence of these traps leads the potential landscape spectrum less separated from the lowest potential state, and therefore less RR. One needs to increase the "temperature" with more fluctuations or decrease transition (or switching) μ (at high μ) in order to "kick" the system out of the traps and therefore increase RR. That is why there is a peak of RR against the variation of the parameter μ.
In Fig 4(c) and 4(d), we plotted the entropy production rate or the dissipation cost of the network, dS dt , against μ and RR. We can see that the sharper the switching is, and therefore the more stable the oscillation is, the more dissipation cost is. The stable oscillation requires more energy consumptions to maintain it. The entropy production rate is the accumulated effects from the combination of both landscape and flux. Therefore the entropy production rate is in general a nonlinear function of these accumulated effects of landscape and flux. In addition, we see that traps consumes more energy and this leads to less stability of oscillation states through RR. Fig 5(a) shows the steady state probability of the biological cell cycle versus the variation of the perturbation parameter c (fixing μ = 5, γ = 60%). We notice that large (small) c indicates large (small) perturbation, and then the probability of biological cell cycle is decreasing with respect to c. It indicates that less perturbation gives more stable biological cell cycle, and therefore a more robust network which is shown in Fig 5(b). Fig 5(c) and 5(d) show the energy dissipation decreases as the perturbation effect increases and as the RR decreases, which shows that the robust system needs more energy consumption against the larger perturbation effect.
γ can be considered as the jumping probability which represents that the state of the fission yeast receives a start signal to begin from the G1 phase. We see in Fig 6(a) that the weights or occupational probabilities of the states on the oscillation path do not change significantly with respect to cycling activation strength γ. Although γ does not change the stabilities of these oscillating states much, it does lead to the directed flow and therefore the stable oscillations as seen in the later flux landscape discussions. Consequently, stimulations through the cycling activation strength γ does not change significantly the shape of the potential landscape as shown in Fig 6(b). Fig 6(c) and 6(d) show that the recycling probability of the cell cycle costs more energy and maintaining a more stable oscillating system requires more energy consumption. Fig 7(a) and 7(b) show the steady state probability of G1 state and the steady state probability of the biological cell cycle path under various stimulation levels (γ) at different switching or fluctuations μ. It has been studied above that μ should characterize the inverse noise strength and increase of which enhances the stability of both the G1 phase and the biological cell cycle path. When considering the changes of jumping probability γ, one can see in Fig 7(a) and 7(b) that the larger the stimulation probability from G1 phase is,the less stable G1 state becomes. The weights or the occupation probabilities of the states on the oscillation paths are also decreased slightly with respect to the cycling activation strength γ, but not significantly compared to that of G1. This means that although the individual states or phases of the cell cycle such as G1 becomes less stable, the stimulation does not change significantly the overall occupations of the states on the oscillation paths and therefore the associated stability.
Power spectrum under different stimulation and sharpness of responses. To illustrate the oscillation cycle, one can introduce the correlation functions of the dynamical observables. For the fission yeast cell cycle, the observables are the expressions of individual genes. This gives the measure on the dynamical response and fluctuations of the systems. The Fourier transform of the autocorrelation function in time of gene expression variables gives the power spectrum. The power spectrum analysis is widely used in the periodical signal transduction mixed with noise interference [45]. The gene expression time series of fission yeast cell cycle can be viewed as the noisy signals. We try to calculate the power spectrum of the fission yeast cell cycle to find how oscillation is influenced by different parameters, and then study their relationship with the flux landscape in the global state space.  The results from Fig 8(a) to 8(f) show that a main peak of the power spectrum emerges and becomes more prominent with the cycling activation strength γ being larger (0 < γ < 0.4) when fixing c = 0.001, μ = 5. This peak corresponds to exactly the speed or frequency of the cell cycle oscillation path. This shows that the intrinsic frequency of the oscillations from the loop flux (the loop flux is proportional to the speed and inversely related to the period of oscillations) coincides with the external response frequency at power spectrum peak. In other words, this results the resonance from the system's external response by power spectrum to reflect the intrinsic frequency of the cell cycle oscillations. When γ becomes larger than 0.4, we also find that the power spectrum will have even more significate peaks. As one peak can map into one periodical oscillation, we can see that with γ becoming higher, there are more and more smaller peaks representing smaller loops emerging, and the biological cell cycle path tends towards the most prominent loop, which can be considered as the dominant loop compared with others.
While by setting c = 0.001, γ = 60% with changing μ , Fig 9(a) to 9(h) give us a clear view that the increasing μ leads to one prominent peak and several other smaller peaks, which again shows that the fission yeast cell cycle contains one main dominant periodical oscillation, and many other much smaller cycles. At even larger stimulations chaos might emerge.

Funneled flux landscape leads to robust limit cycle oscillations
Quantifying the flux landscape. To give a whole picture of the probabilistic flux, we further calculate the flux landscape of this fission yeast cell cycle, and analyze how the flux landscape is influenced by different parameter variations, such as μ, c and γ. By comparing the characteristics of landscape of flux with that of the landscape of potential, we expect to gain unique insights on the non-equilibrium biological cell cycle.
After the decomposition of the driving force into the probability landscape and probability flux, we further decompose the probability flux into the cycle loops. This forms the flux landscape with different loops. The results are shown in Fig 10 (we did not show all the loops for the purpose of clear view on the figure). The thickness of the arrows represents the magnitude of the probability fluxes and the node size represents the steady state probability of that state. We can see that the blue loop can be considered as a dominant loop of the biological cell cycle path with dominant flux flowing along the 10 states of the cell cycle. There are also other secondary loops which can map into those secondary peaks in Figs 8(f) or 9(c). We notice that there is a small number of states clustered together but separated from the major state cluster. These are possible traps of states.
To further quantifying the flux landscape, we show that the flux landscape spectrum in Fig  11. We can see that there is a clear separation seen from the large gap between the flux from "native" cycle and the rest of the others under the chosen parameter set. This means the cell cycle loop stands out from the sea of many loops and becomes dominant. As a result, the flux landscape is funneled towards the dominant loop state. This provides a physical picture for the origin of the limit cycle.
Robustness of flux landscapes against changes in sharpness of the response, self degradations and stimulations. We plot the dominant flux flowing along the biological cell cycle,   {0.1, 0.5, 0.9, 1.0, 2.0, 3.0, 5.0, 9.0} for (a)-(h)  by varying response μ and fixing parameters c = 0.001 and γ = 60% in Fig 12(a). We can see that the increasing μ can directly lead to the increase of the flux. Reminding the discussions above on the robustness of potential landscape (Fig 4(a)), we can state that high probability of switching between states or less environment fluctuations will lead to higher flux of the cell cycle.
We explore the robustness RR of the flux landscape in Fig 12(b) and 12(c) to show how it relates to μ and entropy production rate dS/dt. As we have mentioned, increasing μ can lead to higher probability of the G1 ground state, as well as the cycle flux. Here we see the robustness RR also increases to certain level. The excessively large μ can lead to higher probabilities of other states as traps near G1. This results to the decrease of RR of both potential landscape and flux landscape. Larger fluctuations (smaller μ) from the large μ side will lead to the "escape" from the traps and therefore larger RR as shown in Fig 12(b). Fig 12(c) indicates that more robust oscillating cell cycle often costs more energy to maintain, while traps can consume even more energies. Comparing with the potential landscape shown in Fig 4(b) and 4(d), we obtain Fig 12(d), which shows the consistency of the robustness measure of the potential landscape using states on cell cycle as "native" states and the robustness measure of the flux landscape.
We also calculate how flux landscape is influenced by the perturbation parameter c (by fixing μ = 5, γ = 60%). From Fig 13(a)-13(c), when the perturbation parameter c decreases, the cycle flux increases and the funneled flux landscape becomes more robust, which costs more energy to maintain. Fig 13(d) also shows that the robustness of the flux landscape is monotonic related to the robustness of the potential landscape upon varying of c.
We study the relationship between the stimulation or activation strength from the G1 state to START state and the flux landscape. The flux originates from the the nutrition supply which provides the energy pump (for example, thorough the release of the ATP production). We show in Fig 14(a)-14(d) how the cycling activation strength γ influences the shape or topography of the underlying flux landscape. Before the pumping, the biological cell cycle is an oneway stable pathway to G1 [23,33]. At this stage, even the occupations of the states of the oscillation path including G1 are higher with high RR for potential landscape, the robust directional flow along the cycle has not been formed yet. As a result, no oscillations are emerged. The increase of activation strength can be considered as a switch to form the flux landscape with cycle loops from a deep stable one basin potential landscape. When the pumping increases, many cycle loops start to form and flux landscape starts to emerge. When activation strength γ is large enough the flux and RR reaches saturation. Further pumping will not be effective since robust cycle and flux landscape has already been formed. The average cycle flux increases and the robustness of the flux landscape increases. This indicates that a single flux loop dominates and stands out from the rest. It leads to robust cell cycle. The robust cycle costs more energy to maintain. Fig 14(d) shows that although the increase of the stimulation or pumping from the G1 to the start of the cell cycle leads to a slight decrease of the occupations of the states on the oscillating path (due to the excitations), the separation between the dominant flux loop and the rest of the decoys increases. As a result, the dominant flux loop stands out and forms the directional flow along the yeast cell cycle.
Energy pump, curl flux, dissipation, speed of the cell cycle, and origin of life. To further explore the intrinsic mechanism of flux landscape formation, we study the relationship between the oscillation speed, entropy production and flux. From Fig 8, we can find that when the activation strength γ increase from 1% to 30%, the frequency of most prominent oscillation also increases. We can associate the most prominent frequency and flux. We focus on these relationships in Fig 14(a), 14(b) and 14(d), in which we can see clearly that larger flux leads to the faster cell cycle oscillation. Therefore, the energy pump through γ is the origin of the flux and the flux drives the limit cycle and determines the associated speed or period.
Furthermore, we also explore the relationship between entropy production and flux under activation strength γ changes in Fig 15(c), 15(e) and 15(f). The entropy production rate here represents the energy dissipation in steady state. These figures show that larger flux gives larger energy dissipation and the total entropy production becomes larger. This demonstrates that the degree of detailed balance breaking from the energy pump measured by the flux is the cause of the energy dissipation for sustaining the oscillation.
A fundamental signature of living is the biological replications which can be described by the limit cycle oscillations. From this quantitative study, we can see that energy pump is required for the replications to emerge and survive. Therefore energy supply is necessary for life. This also has evolution implications for the origin of life being initiated and sustained by energy supply.
As a summary, we can state that the limit cycle oscillation is maintained to be stable due to two driving forces: the funneled potential landscape which tends to attract the system down to the close ring valley, leading to high occupations of the states along the oscillation path. The directional flow along the oscillation path is driven by the probability flux originated from the nutrition supply manifested as the stimulation or excitation from G1 to the start of cell cycle. The funneled flux landscape guarantees the clear separation between dominant flux loop and the rest of the other flux loops. Consequently, the dominant flux loop stands out and forms the yeast cell (limit) cycle. Both forces from potential landscape and flux landscape are essential for the stability of the fission yeast cell cycle. Backbone subnetwork of fission yeast from the evaluation of global stability and robustness of both potential landscape and flux landscape We have uncovered that both the potential landscape and the flux landscape are crucial for the stability and the robustness of the limit cycle oscillation of the fission yeast cell cycle. As a practical application, we will perform global sensitivity analysis based on the two landscapes, to explore a backbone subnetwork to carry out the biological functions. We firstly perform perturbations through adding, deleting or repressing arrows between nodes in the wiring diagram in Fig 1, or replacing an active arrow with an inactive arrow, or deleting an individual node. And then we try to analyze the variation of the important characteristic of the two landscapes, such as RR, P G1 , P Circle , RR of flux and so on. Finally, we try to work out which key links or nodes are responsible for the stability, speed (function), and robustness of the cell cycle.
Global stability and robustness of potential landscape under perturbations of mutations and regulation strengths. Fig 16(a) shows the RR versus the probability of the biological cell cycle (fixing μ = 5, γ = 60%) against various perturbations. The perturbations are through adding, deleting or repressing the arrows between the nodes in the wiring diagram in Fig 1, or replacing an activating arrow with an inactivating arrow, or deleting an individual node. We see the larger the RR is, the higher the occupation is of the states on oscillation cell cycle upon perturbations of links and nodes. This indicates the more stabilities of the states on the cell cycle. This provides the rational of using RR as a robustness measure for the cell cycle network. Fig 16(b) shows again that more stable oscillations requires more energy consumption. The points of low dissipations with large RR values are ignored, as the corresponding probabilities of the biological cycle path are low. Fig 16(c) shows the relationship between the steady state probability of the biological cell cycle path and the steady state probability of the G1 state upon perturbations of links and nodes. The stable G1 often correlates with higher occupations of the states on oscillating cycle paths, as G1 is the starting point in the biological cycle path. However the higher occupations of the biological cycle do not always imply stable G1 state, as quite a few perturbations will disable the stability of state G1 and enhance the stabilities of some other states along the cell cycle path. We see that upon perturbations on links and nodes of the fission yeast cell cycle network, higher probabilities of states on oscillation path often accompany with higher robustness ratio of the flux landscape separating the cell cycle loop from the rest as shown in Fig 17. There are some exceptions where very high and very low probabilities of the states on oscillation cell cycle path have varying robustness ratios of the flux landscapes. We can ignore those due to the insignificant cell cycle or insignificant decoys.

Backbone subnetwork contained in the fission yeast networks
Through the global stability analysis for the key wirings of the networks upon perturbations of links and nodes, one can identify the key network structure elements or motifs responsible for the stability and biological function. To further identify the stable and functional backbone subnetwork, we choose to delete the link one by one and find out which will make the network more unstable. Based on the discussion in the text above, we select three essential elements to measure the importance of each network edge, that is RR, P G1 , and RRflux (RR of flux landscape), which can represent not only the global stability of G1 state but also the robustness of the biological path.
To measure the robustness of each edge, we performed the orders as follow: First of all, for each edge, we respectively calculate the difference of the three essential elements between mutated network and original network. Then we rank the difference of the 27 edges and give them score from 1 to 27. The larger is the difference, the larger the score is. Then we rank the summation of all the three scores as a total evaluation score(TES), and obtain the rank of robustness of all the 27 edges of the fission yeast network. (One can see the result in Table 3) Based on the edge robustness, we attempt to reconstruct a minimal but most robust or stable backbone network of the fission yeast. However, nature does not necessarily use the robust edges to build the network, as the aim of a network is to perform certain biological function. Therefore, we need to reconstruct a biological meaningful network based on the major functional biological path for cell cycle. There are several strategies to do so.
In an approach [46], the backbone network is obtained directly from the criterion of emergence of major biological path for cell cycle. However, in our approach, we chose the backbone subnetwork based on the global sensitivity and robustness from potential and flux landscape for quantitatively guarantee a stable biological path. We have labeled each edge of this subnetwork both in the Table 3 with italic type and the original network wiring in Fig 18 with red color.
Notice that our backbone subnetwork has large overlaps with the one obtained by Zeng et al. [46] through biological path requirement. This shows that the backbone network chosen from biological path is also likely the stable one as we have quantitatively shown here for this case of fission yeast cell cycle. This result seems not so intuitive. As the backbone subnetwork directs to the minimal network to perform the biological function, while the robust network means that it has strong capability to persist in the mutation or fluctuation environments. Therefore, the minimal network should not always lead to the most robustness one.
In our study, we explore parameters listed in the Table 3 of both potential and flux landscapes. The three elements which contribute the TES are the key features representing the potential and flux landscape topography. Therefore, the rank of the TES is calculated in a quasi-quantitative way to gain insights on both potential stability and period persistence for each edge. It suggests that due to the consideration of flux landscape in the periodical dynamics, the minimal backbone subnetwork with highest rank of TES tends to be a global robust one to perform the cell cycle function. Therefore, we state that this study gives a physical principle and basis in terms of the potential and flux landscape for the backbone finding.
Furthermore, based on the global sensitivity analysis, we can identify key links that change significantly the occupation probabilities of the states on the cell cycle path and robustness ratio separating the dominant cycle loop from the rest compared to the wild type. These links and nodes are responsible for biological function of the cell cycle.
As we have stated, for this case of oscillation network, the flux landscape have a large impact on the dynamical behavior of the network. In Table 3, if we delete the index of flux landscape, Table 3. Robustness ranking of 27 links of the fission yeast networks. RRflux  Fig 18. The minimal robustness backbone subnetworks of the fission yeast networks. All the links marked by red color contribute a backbone subnetwork, which is generated from the evaluation of robustness in Table 3 (marked by italic type). The remaining links form a residual auxiliary subnetwork. All the signs of the links are the same as in Fig 1. https://doi.org/10.1371/journal.pcbi.1005710.g018

RR P G1
i.e. RRflux, there are several edges in the original network acutely changing their orders of importance. The edges such as SK to Rum1, and Cdc25 to Cdc2/Cdc13 Ã will lose their high ranks due to the missing of flux index, while the remaining edges such as Cdc2/Cdc13 to Rum1 and to Ste9 and others tend to squeeze into the subnetwork. This leads to the conclusion that the former edges tend to perform the biological cycle function while the remaining edges tend to keep the robustness of G1 state. It is in this perspective these results provide a strong support for the potential and flux landscape theory in the study of cell cycle. Cell cycle is a hallmark of cancer. Cancer cells has a much faster speed of cell cycle than normal cells. Therefore, regulating the cell cycle speed is crucial for preventing and curing the cancer. From above global sensitivity analysis, we can identify the key nodes and links in the fission yeast cell cycle network for regulating the cell cycle speed. These key links and nodes form the backbone network of the cell cycle. Therefore, we can based on this to do network design and network medicine discovery targeting the cancer.

Conclusions
We explore the global natures of the networks. We found the network dynamics and global properties are determined by two essential features: the potential landscape and the flux landscape. While potential landscape quantifies the probabilities of different states forming hills and valleys, the flux landscape quantifies the probability fluxes of different loops flowing through states. These two landscapes can be quantified through the decomposition of the dynamics into the detailed balance preserving part and detailed balance breaking part. While funneled landscape is crucial for the stability of the single attractor networks, the argument can be extended to the stabilities of the states on the oscillation paths by including them in the same (line) basin of attraction. Importantly, we have uncovered that the funneled flux landscape is crucial for the stable and robust oscillation flow.
This provides a new interpretation of the origin of the limit cycle oscillations: There are always many cycles and loops forming the flux landscapes, each with a probability flux going through the loop. The oscillation only emerges when one specific loop stands out and carries much more probability flux than the rest of the others.
We studied the fission yeast cell cycle as an example to illustrate the idea. We found both the potential landscape and the flux landscape of the fission yeast cell cycle oscillations are funneled, which guarantees the global stability. While the funneled potential landscape guarantees the stabilities of the states on the oscillating path, the funneled flux landscape guarantees the directional flow of the oscillations which breaks the detailed balance and time reversal symmetry, leading to the stand out of the dominant flux loop against others. The stability and robustness of the oscillations are quantified through a dimensionless ratio of the steepness or gap versus the averaged variations or roughness of the landscape (measuring funnelness as we termed as robustness ratio RR).
We explore how RR changes with respect to the stimulations, self degradations, state switching rate or fluctuations, and changes in topology of the network (wirings). This allows us to identify the key factors and structure elements of the networks in determining the stability, speed and robustness of the fission yeast cell cycle oscillations.
Based on the global sensitivity analysis, we obtain that our most robust subnetwork is nearly the same as the minimal biological functional network, and by setting the cell cycle period as the evolution goal, we suggest the fission yeast should follow this evolution goal to form a 27-link network with faster period but not using minimal backbone network. We see that the non-equilibriumness characterized by the degree of detailed balance breaking from the energy pump quantified by the flux is the cause of the energy dissipation for initiating and sustaining the replications essential for the origin and evolution of life. Finally we are looking forward to the good future by controlling the speed of the cell cycle as an important in designing targeting drugs for preventing and curing the cancer.