New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Potential landscape and flux framework of nonequilibrium networks: Robustness, dissipation, and coherence of biochemical oscillations

Edited by José N. Onuchic, University of California at San Diego, La Jolla, CA, and approved June 3, 2008 (received for review January 22, 2008)
Abstract
We established a theoretical framework for studying nonequilibrium networks with two distinct natures essential for characterizing the global probabilistic dynamics: the underlying potential landscape and the corresponding curl flux. We applied the idea to a biochemical oscillation network and found that the underlying potential landscape for the oscillation limit cycle has a distinct closed ring valley (Mexican hatlike) shape when the fluctuations are small. This global landscape structure leads to attractions of the system to the ring valley. On the ring, we found that the nonequilibrium flux is the driving force for oscillations. Therefore, both structured landscape and flux are needed to guarantee a robust oscillating network. The barrier height separating the oscillation ring and other areas derived from the landscape topography is shown to be correlated with the escaping time from the limit cycle attractor and provides a quantitative measure of the robustness for the network. The landscape becomes shallower and the closed ring valley shape structure becomes weaker (lower barrier height) with larger fluctuations. We observe that the period and the amplitude of the oscillations are more dispersed and oscillations become less coherent when the fluctuations increase. We also found that the entropy production of the whole network, characterizing the dissipation costs from the combined effects of both landscapes and fluxes, decreases when the fluctuations decrease. Therefore, less dissipation leads to more robust networks. Our approach is quite general and applicable to other networks, dynamical systems, and biological evolution. It can help in designing robust networks.
Biological rhythms exist on many levels in living organisms. The study of oscillation behavior in an integrated and coherent way is crucial to understanding how rhythms function biologically (1, 2).
Biological clock dynamics is often described by a network of deterministic nonlinear chemical reactions of the corresponding averaged protein concentrations in the bulk. In a cell, there is a finite number of molecules. Thus, intrinsic statistical fluctuations can be significant for dynamics. However, external fluctuations from highly dynamical and inhomogeneous cellular environments can also be important (3). It is therefore important to investigate the roles of statistical fluctuations in the robustness and stability of oscillation. Instead of the averaged deterministic dynamics, we need to develop a probabilistic description to model the corresponding cellular process. This can be realized by constructing a master equation for the intrinsic fluctuations or a diffusion equation for external fluctuations for probability evolution (4, 5). Even for intrinsic fluctuations, we can simplify the master equation into a diffusion equation in the weak noise or large number limit (5, 6).
By solving the diffusion equation, we can obtain the time evolution and longtime steady state of the probability distribution in protein concentrations of the network. In analogy to equilibrium systems, the generalized potential can be shown to be closely associated with the steadystate probability of the nonequilibrium network, with a few applications (6–19). Once the network problem is formulated in terms of potential landscape, the issue of the global stability or robustness is much easier to address (16–19). Although deterministic dynamics might be nonlinear and chaotic, the corresponding probabilistic distribution obeying linear evolution equations is usually ordered and can often be characterized globally. An interesting study is described in ref. 14 where the stability and emergence of competence cycles in Bacillus subtilis are controlled by the intrinsic statistical fluctuations and a nonconventional stochasticity from nonadiabatic conditions with finite binding/unbinding rates of the proteins to DNA. The resulting cycle dynamics between stable vegetation and competence state is incoherent (defined competent state duration but stochastic intervals between) in low and coherent in high Com K expression. In this work, we will focus on coherent dynamics of limit cycles with certain periodicity.
Landscape and Flux Framework for Nonequilibrium Networks
Landscape ideas were introduced for uncovering global principles in biology for protein dynamics (20), protein folding, and interactions (21, 22). All of these ideas were based on a quasiequilibrium assumption with known potentials. For a nonequilibrium open system constantly exchanging energies and information with outside environments, the potential landscape is not known a priori and needs to be uncovered. Even when the probability landscape could be discussed such as in population dynamics and developmental biology (23–25), it was not clear the relationship of landscape with dynamics. Furthermore, probability flux that is zero in equilibrium case now becomes significant. It is the purpose of this article to study the global robustness and physical mechanism of nonequilibrium network through the introduction of the concepts and quantifications of the potential landscape and the nonequilibrium probability flux, with an example of oscillations against fluctuations in the cell.
The conventional way of describing the dynamics of a network is to write down the underlying chemical rate equations:
Here, we provide a general way to find potential and flux. As mentioned, the cellular network is under intrinsic and external fluctuations (3). The dynamics of the network system is therefore more accurately described by the probabilistic approach:
Instead of deterministic trajectories, we will focus on the probabilistic evolution of diffusion equation covering the whole concentration space:
If the steady state exists that is true for many network systems (see SI Text for conditions),
For nonequilibrium systems in general, however, in steady state, ∇·J(x, t) = 0 does not necessarily mean that J has to vanish; there is no guarantee that the detailed balance condition J = 0 is satisfied. In general, the divergencefree nature implies flux J is a rotational curl or more precisely recurrent field: for example, in three dimensions, J = ∇ × A, with nonzero curl of vector A [for higher dimensions, see Hodge decomposition (28)]. This implies, from the definition of the flux, J_{ss} = FP_{ss} − D·
For nonequilibrium networks, the dynamics and global properties are therefore determined not only by gradient of potential landscape but also by the divergent free curl flux field. This may provide the missing link between the probability landscape and the underlying dynamics of networks (for example, in population evolution dynamics). The dynamics of a nonequilibrium network spirals (from flux) down the gradient (from potential) instead of only following the gradient as in the equilibrium case, just like electrons moving in both electric and magnetic field. As we shall see, the best example of illustrating the interplay of both potential landscape and flux in action is the oscillatory network.
Landscape and Flux of Biochemical Oscillation Network
To explore the nature of the oscillation mechanism, we will study a simplified yet important example of biochemical network of cell cycle: a periodic accumulation and degradation of two types of cyclins during the division cycle in budding yeast. The oscillations connected to dynamical interactions between CLNtype cyclins and CLBtype cyclins were found (29). CLN/CDC28s, which are CLNtype cyclins associated with Cdc28kinase, activate their own synthesis (“selfactivation”) and inhibit the degradation of CLB/CDC28s, which are CLBtype cyclins associated with CDC28 kinase. As the concentration of CLB/CDC28 becomes larger, it inhibits the synthesis of CLN/CDC28. The mutual interplay of CLN/CDC28 and CLB/CDC28 generates the periodic appearance of their associated kinase activities, which drive bud emergence, DNA synthesis, mitosis, and cell division of the budding yeast cell cycle. Fig. 1 shows the mechanism of CLN/CDC28 and CLB/CDC28 oscillation. CLN/CDC28 and CLB/CDC28 subunits are limited by cyclin availability, because kinase CDC28 is excessive (1).
For the protein network, based on the Michaelis–Menten enzyme kinetic equations, one can derive a set of differential equations that describe the variation rate of each component's concentration in the network. We have two independent simplified equations: In these equations, X_{1} and X_{2} are the average concentration of CLN/CDC28 and CLB/CDC28, respectively. The k′'s are the rate constants, K's are the equilibrium binding constants, and the J's are the Michaelis constants. The first term describes the synthesis, and the last term describes the decay. In terms of dimensionless variables (x_{1} = X_{1}/K_{m}, x_{2} = X_{2}/K_{n}, t′ = v_{1}t/K_{m}), these equations become where a = k_{2}K_{m}/v_{1}, b = k_{3}/(k_{4}K_{n}), c = (K_{m}/K_{j}) 2, and τ_{0} = v_{1}/(k_{4}K_{m}).
The corresponding diffusion equation for probability distributions of protein concentrations for x_{1} and x_{2} with noise due to fluctuations is
Here, D is the diffusion coefficient tensor (or matrix) assumed to be homogeneous and isotropic constant for simplicity (D_{11} = D_{22} = D and D_{12} = D_{21} = 0). The associated flux vector components in twodimensional protein concentration space are J_{1}(x_{1}, x_{2}, t) = F_{1}(x_{1}, x_{2})P − D
We fix all parameters except b and c. b represents the relative effectiveness of production of x_{2}, and c represents the relative effectiveness of inhibiting x_{2}'s degradation. The other parameter values are a = 0.1, ε = 0.1, and τ = 5.0. Fig. 2 shows the phase plane of b and c from the analysis of the deterministic equations (Eqs. 1 and 2). We can see that the system has three phase regions: an unstable limit cycle oscillation phase, a bistable phase, and a monostable phase. Large b and large c lead to effective inhibition of x_{1} production and leave only with degradation of x_{1}, and therefore yield a monostable decay. Smaller b and c can provide the balance between the activation and degradation of x_{1}. Therefore, when b is fixed to be small and c is large, oscillation emerges. In contrast, when c is fixed to be small and b is large, bistability emerges. When both c and b are small, there is neither effective production of x_{2} nor effective inhibition of x_{2}'s degradation. This leads to effective production of x_{1} and again monostability. We choose for convenience a set of specific parameters b = 0.1, c = 100, at which the fixed point is unstable and a limit cycle emerges.
To solve the diffusion equation, we impose the reflecting boundary condition, J = 0, in this work. We have also explored absorbing boundary conditions and obtained similar solutions. By giving certain initial conditions (either homogeneous or inhomogeneous) and taking the long time limit, we obtained the same steadystate solution using the finite difference method.
As we discussed before, from the steadystate probability distribution P, we can identify U(x) = −ln P(x, t → ∞) = −ln P_{ss} (when
The potential landscape becomes flatter as the diffusion constant D grows larger, as indicated by the shallower energies along the closed ring compared with both the inside and the outside of the closed ring. The landscape transforms from a distinct irregular inhomogeneous closed ring valley into a flatter structure as shown in Fig. 3B. It implies that, when the system is under larger fluctuations, there is more freedom to go to other states and therefore less attraction to the deterministic oscillation path. Less time is then spent on that path. The resulting underlying landscape departs from the clear closed ring valley structure of oscillations, reflecting the large fluctuations. When the diffusion coefficient increases, the attraction to the limit cycle will be weaker; conversely, the weaker the fluctuation, the more robust the oscillation.
We can also obtain the steadystate probability flux, another essential quantity for the network system once we get the steadystate solution of the diffusion equation. At the steady state, there is a circulating flow with nonzero curl as shown in Figs. 3 and 4. In Fig. 3, the blue arrows represent the steadystate probability flux and the white arrows represent the force from negative gradient of the potential landscape. Fig. 4A shows both the magnitude and direction of the flux; Fig. 4B shows only the direction of the flux flow. The magnitude of the flux is small inside and outside of the closed ring but significant along the ring (Figs. 3 and 4). The direction of the flux near the ring is parallel to the oscillation path. The forces from negative gradient of the potential landscape are insignificant along the closed ring and significant inside and outside of the ring. So, inside and outside of the closed ring valley, the network is attracted by the landscape toward the closed ring. Along the closed ring valley, the network is driven by the curl flux flow for oscillation.
The magnitude and direction of the residual curl flux force F′(x) = F + D·∇U is also shown in Fig. 4C. Fig. 4D shows the direction of the residual force F′(x). We see that the direction of curl flux J is parallel to that of the residual force F′(x). This is expected from the force decomposition discussed earlier: F + D·∇U = J_{ss}/P_{ss}. The residual force is thus parallel to the flux J and is the driving force for the curl field of probability flux.
Without the landscape's gradientpotential force (Figs. 3 and 4 E and F), the system will not be attracted to the oscillation ring (major stages such as G_{1}, S, G_{2}, and M of cell cycle). Without the curl flux driving the system (nutrition supply as the pump), the system will get stuck in low potential valleys on the ring without moving further(check points in cell cycle), and oscillation will not occur (Figs. 3 and 4 A–D). There is an important interplay between the dominant attractive force from the landscape inside as well as outside the closed ring and the dominant driving force from the flux along the closed ring. So both landscape and flux are necessary to characterize this kind of nonequilibrium system, and this oscillation (of cell cycle) provides an excellent illustration of that necessity.
We also notice that, when the diffusion coefficient D is small, the curl flux J_{ss} is almost parallel to real force F(x). This is because the gradient component of the force is proportional to D and the residual force gives dominant contribution to total force when D is small. In this case, dx/dt ∼ J_{ss}/P_{ss}, so the period of oscillation can be approximated through the loop integral of inverse flux along the oscillation path: T ∼ ∮dl/(J_{ss}^{l}/P_{ss}). This provides a possibility through the observation of oscillation period and local speed to explore the natures of the network flux.
Transition Time, Barrier Height, and Robustness
We now study the stability and robustness of the network. The stability is related to the escape time from the basins of attraction. Because the system is characterized by the basins of attractor with large weights, the easier it is to escape, the less stable is the system. For the probabilistic description of the network above with diffusion equation, the mean firstpassage time for escape τ(x_{1}, x_{2}) starting from the point (x_{1}, x_{2}) obeys (30) F(x)·∇_{x}τ + D·∇_{x}∇_{x}τ = −1, and in our case of two dimensions: It is essentially the average time it takes from a initial position to reach a given final position. The equation can be solved by an absorbing boundary condition at the given site and reflecting boundary conditions for the rest.
For an equilibrium system, the barrier height on the potential landscape is intimately related to the escape time by Arrhenius law. The question is, Will there be still a direct relationship between the escape time and barrier height for a nonequilibrium network? If so, the landscape topography will then provide a quantitative measure of the hardness of the system to escape from the limit cycle attractor to outside and therefore of the stability and robustness.
We define the barrier heights as barrier1 = U_{fix} − U_{min} and barrier2 = U_{fix} − U_{max}. U_{max} is the potential maximum along the limit cycle attractor. U_{min} is the potential minimum along the limit cycle attractor. U_{fix} is the potential at the local maximum point inside the limit cycle circle. In Fig. 5A, as the diffusion coefficient characterizing the fluctuations decreases, the barrier heights associated with escaping from the limit cycle attractor barrier1 and barrier2 are higher. In Fig. 5B, we see a direct relationship between the escape time and landscape barrier heights for nonequilibrium network: as the barrier for escape becomes higher, the escape time becomes longer. The resulting limit cycle attractor becomes more stable because it is harder to go from the ring to the outside. Therefore, small fluctuations and large barrier heights lead to robustness and stability in the oscillatory protein network.
Entropy Production, Barrier Height, and Robustness
The nonequilibrium network is often an open system exchanging information and energies with its surroundings. Therefore, the nonequilibrium steady state dissipates energy and causes entropy, just as electric circuits dissipate heat due to the action of both voltage (potential) and current (flux). Therefore, dissipation can be determined globally with both the landscape and flux. In the steady state, the heat loss rate is equivalent to the entropy production rate. We will explore the dissipation cost via entropy production rate at the steady state (details in SI Text) (18, 19, 31, 32). Fig. 6A shows the entropy production rate for different diffusion coefficients. We can see that the dissipation or the entropy production rate decreases when the diffusion coefficient decreases. This implies that when the fluctuations of the systems become smaller, the associated dissipation is smaller. Because fewer fluctuations lead to more robust oscillations as shown above in Fig. 5, less dissipation should be closely linked with fewer fluctuations and a more stable network. Indeed, we see that less dissipation leads to higher barrier heights barrier1 and barrier2, and therefore a more stable network (Fig. 6B). Because the entropy production is a global characterization of the network, minimization of the dissipation cost might serve a design principle for evolution of the network. It is intimately related to the robustness of the network.
Period, Amplitude, and Coherence of Oscillations Against Fluctuations
To address more of the robustness of the oscillations, we study the chemical reaction network equations under the fluctuating environments by simulating the stochastic dynamics for different values of D. That is, we follow the stochastic Brownian dynamics instead of the deterministic average dynamics. Fig. 7 A and B shows the distributions of the period of oscillations calculated for 1,400 successive cycles. We can see that the distribution becomes more spread out with a mean period μ that is still close to the deterministic period of oscillations (μ = 343.3) when the fluctuations increase. The standard deviation σ from the mean increases, and more other possible values of the period of oscillations can appear when the fluctuations increase (Fig. 7C) (2). This implies that fewer fluctuations lead to more coherent oscillations with single period instead of multiple periods. We also see that the period distribution becomes less dispersed when the entropy production rate is less. This shows that a less dissipated network can lead to a more coherent oscillation with a unique period instead of a distribution of periods. We see also that higher barrier heights lead to less dispersed period distribution (Fig. 7D). All of these show that more a robust network leads to more coherent oscillations focusing on a single rather than multiple periods.
We also show the distributions of the amplitude for x_{2} as D increases. The distribution becomes more dispersed but keeps the same mean value close to the deterministic amplitude of x_{2} [Amplitude(x_{2}) = 4.79], as the fluctuations increase in Fig. 8A. The standard deviation σ increases when D goes up in Fig. 8B. This also shows that less fluctuation leads to a more robust and coherent oscillation.
The robustness of the oscillation can be quantified further by the phase coherence ξ that measures the degree of periodicity of the time evolution of a given variable (33) (details in SI Text). In the presence of fluctuations, the more periodic the evolution, the larger the value of ξ. In Fig. 9A, ξ decreases when the diffusion coefficient increases. This means that larger fluctuations tend to destroy the coherence of the oscillations and therefore the robustness. In contrast, in Fig. 9B, ξ increases when barrier heights increase, showing that a more robust network leads to more coherent oscillations. ξ decreases when entropy production increases. Dissipation tends to destroy coherence.
Conclusions
We have uncovered the underlying potential landscape and flux of nonequilibrium networks, crucial for determining its dynamics and global robustness. The landscape of the oscillation network has a closed ring valley shape attracting the system down, and the curl flux along the ring is the driving force for oscillation. The potential barrier height, shown here to be correlated with escape time, provides a measure of the likelihood of escaping from the limit cycle attractor, which determines the robustness of the network.
We observe that when the fluctuations increase, global dissipation, period, and amplitude variations increase and barrier height and coherence decrease. The period and amplitude variations can be experimentally measured and compared with theoretical predications (refs. 34–37; similar topology as here in refs. 35 and 37). Furthermore, minimizing the dissipation costs may lead to a general design principle for robust networks. The framework and methods in this article can be applied to more complicated and realistic networks and dynamical systems to explore the underlying global potential landscape and flux for probabilistic population dynamics and biological evolution.
Acknowledgments
We acknowledge the use of VCell software for part of the calculations in this work. J.W. thanks Profs. H. Qian, P. Ao, and P. Mitra for discussions. J.W. thanks the National Science Foundation for a Career Award. L.X. and E.W. are supported by National Natural Science Foundation of China Grants 90713022 and 20735003.
Footnotes
 ^{‡}To whom correspondence may be addressed. Email: jin.wang.1{at}stonybrook.edu or ekwang{at}ciac.jl.cn

Author contributions: J.W. and E.W. designed research; J.W. and L.X. performed research; J.W. and E.W. contributed new reagents/analytic tools; J.W., L.X., and E.W. analyzed data; and J.W., L.X., and E.W. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/cgi/content/full/0800579105/DCSupplemental.
 © 2008 by The National Academy of Sciences of the USA
References
 ↵
 Fall CP,
 Marland ES,
 Wagner JM,
 Tyson JJ
 ↵
 Gonze D,
 Halloy J,
 Goldbeter A
 ↵
 Swain PS,
 Elowitz MB,
 Siggia ED
 ↵
 ↵
 ↵
 Qian H,
 Saffarian S,
 Elson EL
 ↵
 Haken H
 ↵
 Graham R
 Moss F,
 et al.
 ↵
 ↵
 Sasai M,
 Wolynes PG
 ↵
 Walczak AM,
 Onuchic JN,
 Wolynes PG
 ↵
 ↵
 ↵
 Schultz D,
 Jacob EB,
 Onuchic JN,
 Wolynes PG
 ↵
 ↵
 Wang J,
 Huang B,
 Xia XF,
 Sun ZR
 ↵
 ↵
 ↵
 Lapidus S,
 Han B,
 Wang J
 ↵
 Frauenfelder H,
 Sligar SG,
 Wolynes PG
 ↵
 Wolynes PG,
 Onuchic JN,
 Thirumalai D
 ↵
 ↵
 Fisher RA
 ↵
 Wright S
 ↵
 Waddington CH
 ↵
 Arfken GB,
 Weber HJ
 ↵
 ↵
 Bott R,
 Tu LW
 ↵
 ↵
 Risken H
 ↵
 Qian H
 ↵
 Han B,
 Wang J
 ↵
 ↵
 ↵
 ↵
 Rust MJ,
 Markson JS,
 Lane WS,
 Fisher DS,
 O'Shea EK
 ↵
Citation Manager Formats
Sign up for Article Alerts
Jump to section
 Article
 Abstract
 Landscape and Flux Framework for Nonequilibrium Networks
 Landscape and Flux of Biochemical Oscillation Network
 Transition Time, Barrier Height, and Robustness
 Entropy Production, Barrier Height, and Robustness
 Period, Amplitude, and Coherence of Oscillations Against Fluctuations
 Conclusions
 Acknowledgments
 Footnotes
 References
 Figures & SI
 Info & Metrics