## Abstract

Recent studies on stochastic oscillations mostly focus on the power spectral analysis. However, the power spectrum yields information only on the frequency of oscillation and cannot differentiate between a stable limit cycle and a stable focus. The cycle flux, introduced by Hill (Hill 1989 *Free energy transduction and biochemical cycle kinetics*), is a quantitative measure of the net movement over a closed path, but it is impractical to compute for all possible cycles in systems with a large state space. Through simple examples, we introduce concepts used to quantify stochastic oscillation, such as the cycle flux, the Hill–Qian stochastic circulation and rotation number. We introduce a novel device, the Poincaré–Hill cycle map (PHCM), which combines the concept of Hill’s cycle flux with the Poincaré map from nonlinear dynamics. Applying the PHCM to a reversible extension of an oscillatory chemical system, the Schnakenberg model, reveals stable oscillations outside the Hopf bifurcation region in which the deterministic system contains a limit cycle. Bistable behaviour is found on the small volume scale with high probabilities around both the fixed point and the limit cycle. Convergence to the deterministic system is found in the thermodynamic limit.

## 1. Introduction

The chemical master equation (CME), first studied by Delbrück (1940), has become the leading mathematical theory for studying mesoscopic nonlinear chemical reaction systems with small volume of the order of a living cell (Beard & Qian 2008). The chemical reactions are characterized in terms of the stochastic copy numbers of the various dynamic chemical species. The CME provides the equation for the time-dependent joint probabilities of the number of molecules while the Gillespie (2007) algorithm gives the stochastic trajectories. The latter to a CME is the same as a stochastic differential equation to a Fokker–Planck equation.

We have recently investigated the bistability in connection with a simple nonlinear chemical reaction system: the Schlögl model (Vellela & Qian 2009). Mathematically, the Schlögl model is a one-dimensional stochastic dynamical system on , the space of integers. One is naturally interested in nonlinear stochastic dynamics on with limit cycles (Qian *et al.* 2002). In this paper, we study a simple two-dimensional chemical oscillator, the Schnakenberg (1979) model, which is one of the widely studied nonlinear chemical limit cycles with Hopf bifurcation (Murray 2002).

We introduce the Poincaré–Hill cycle map (PHCM) to study biased rotational random walks in , particularly those which exhibit some kind of sustained oscillation. The PHCM is a natural generalization and combination of both the deterministic Poincaré return map in nonlinear dynamics (Strogatz 2001) and Hill’s cycle fluxes in stochastic kinetics (Hill 1989).

Biochemical dynamics inside a single living cell has become increasingly quantitative in recent years owing to the technological advances in optical imaging, single cell analysis and green fluorescence proteins. Such dynamics are often nonlinear and stochastic, exhibiting complex behaviour. How to understand sustained chemical oscillations in the presence of stationary noise is a challenging problem for both experimentalists and theorists. Interest in oscillations of biochemical kinetic systems has given rise to the active research area of stochastic (or coherent) resonance (Gammaitoni *et al.* 1998); see Qian & Qian (2000), Zhu *et al.* (2006), Zia & Schmittmann (2007) and Wang *et al.* (2008) for recent discussions on this subject. The emerging consensus is that both stationary probability distribution and stationary flux are important in understanding stationary stochastic dynamics that exhibit nonlinear oscillation.

This paper is organized as follows. Section 2 contains important concepts and illustrates their use through examples. Section 3 introduces the PHCM and a reversible version of the Schnakenberg model. In §4, the PHCM is applied to this model and results are presented. Section 5 provides discussions.

## 2. Background concepts in stochastic oscillations

In this section, we first give some background that is pertinent to the concept of the PHCM. Using the harmonic oscillator as an example of a nonlinear dynamics on the circle, we will illustrate the existence of stochastic oscillations in both underdamped and overdamped conditions. We also define one-way cycle fluxes (Hill 1989) in terms of an exit problem in a single-cycle Markov chain, and show that the cycle fluxes are a measure of stochastic oscillations. We then introduce a concept from the mathematical theory of irreversible Markov processes called *stochastic circulation* (Jiang *et al.* 2004). This term first appeared in Tomita & Tomita (1974) and Hill & Chen (1975) and its theory has been extensively developed by Hill and Qian (Hill 1989; Jiang *et al.* 2004). The Hill–Qian circulation is shown to be equivalent to the cycle flux in a two-dimensional random walk problem with rotational symmetry.

### (a) Nonlinear dynamics on the circle

The most widely studied oscillator is the damped harmonic system. In polar coordinates, it is represented by the system of differential equations
2.1
in which and *a*=*η*/2*m*. The parameters *m*, *k* and *η* are the mass, spring constant and frictional coefficient. Standard textbooks show that *η*^{2}/4*k**m*=(*a*/*ω*)^{2}>1 corresponds to the overdamped (non-oscillatory) case and (*a*/*ω*)^{2}<1 corresponds to the underdamped (oscillatory) case.

There is a different perspective. As the equations in (2.1) are not coupled, we can study the angular movement, *θ*(*t*), as nonlinear dynamics on a circle, independent of *r*. One can find a saddle–node bifurcation at |*a*|=|*ω*|. This leads to two scenarios.

(i) When |

*a*|>|*ω*|, there are two ‘fixed points’ in [0,*π*): one stable and one unstable. There is no continuous rotation. One can compute the*rotation number*2.2 and find*α*=0. In terms of Cartesian coordinates, the stable and unstable fixed points represent the directions of the eigenvectors corresponding to the larger and smaller eigenvalues, respectively.(ii) When |

*a*|<|*ω*|, there are no fixed points and*θ*(*t*) is a continuous rotation on the circle. In fact, in this case the rotation number is*α*=(*ω*^{2}−*a*^{2})^{1/2}/2*π*, which is precisely the frequency of the underdamped system.

We see a mutual exclusion between the fixed points and the rotation number: when there is rotation, there are no fixed points, and when there are fixed points, there is no rotation. It turns out that this mutual exclusion is not valid when noise is present. If we associate the fixed points with the extrema of the probability distribution, and rotation number with cycle flux, then the stochastic dynamics on the circle can have a non-uniform probability distribution while at the same time a non-zero cycle flux (Qian *et al.* 2000).

### (b) One-way cycle fluxes of stochastic dynamics on the circle

Let us now consider an *n*-state kinetic cycle
2.3
Starting at state 1, one is interested in the respective rates of the system completing a forward cycle 1,2,…*n*,1 and a backward cycle 1,*n*,…2,1. By completing a cycle, we mean the system can go back and forth, as long as it first returns to 1 by finishing the cycle in the appropriate direction. In single-molecule enzymology, the forward and backward cycles precisely correspond to the conversion of one reactant to a product, and *vice versa* (Qian & Xie 2006).

The one-way forward and backward cycle fluxes are the number of cycles completed per unit time on average. Hill was among the first to formulate the cycle flux problem into the canonical exit problem in probability in which completing a forward or backward cycle is represented by reaching the state 1* or 1′, respectively. The kinetic scheme is 2.4 In terms of the scheme in equation (2.4), the following has been shown.

(i) Starting in the central state , there are exit (absorbing) probabilities

*p*′ and*p**, where*p*′+*p**=1 and 2.5(ii) The normalized probability density functions for the time of absorption at 1′ and 1* are identical, with a common mean time 〈

*T*〉.

Then, for the cycle kinetics in equation (2.3), one has the forward and backward cycle fluxes (Qian 2008)
2.6
and the net forward cycle flux *J*=*J*^{+}−*J*^{−}. When the cycle kinetics is in detailed balance, the right-hand side of equation (2.5) is 1. Hence, *J*^{+}=*J*^{−} and the net forward cycle flux is zero. This is the only case in which there is no net rotational motion. If all the *u*_{i}=*u*, *w*_{i}=*w* and *u*>*w*, then *p*′/*p**=(*w*/*u*)^{n} is exponentially small for large *n*. Thus, in the large system limit (), one can essentially neglect the backward cycle.

In the context of the stochastic dynamics on the circle, when *u*_{i}≥*w*_{i} (or *u*_{i}≤*w*_{i}) ∀*i*, it corresponds to the case of an underdamped system and there are no fixed points. Otherwise, at least one ‘fixed point’ exists, and it is overdamped. An underdamped system has more pronounced oscillations in one direction and exponentially small backward (or forward, depending on the ratios *u*_{n}/*w*_{n}) cycle flux as *n* increases. Therefore, in the large system limit, it has a sustained rotation. An overdamped system with increasing *n*, however, has more and more balanced forward and backward cycle fluxes. Thus, in the limit the rotation becomes zero.

### (c) Hill–Qian stochastic circulation and rotation number

Hill’s one-way cycle flux has been generalized to the mathematical concept of stochastic circulation in the theory of Markov processes (Jiang *et al.* 2004). When applied to the dynamics on the circle in §2*a*, it turns out to be the natural generalization of the deterministic rotation number given in equation (2.2), divided by 2*π* (Wang *et al.* 1997). Consider the stochastic counterpart of the angular dynamics in equation (2.1)
2.7
where *B*(*t*) is the standard Brownian motion. One considers the stationary process *Θ*(*t*) and counts the number of times it completes the forward cycle, *N*^{+}(*t*), and backward cycle, *N*^{−}(*t*), up to time *t*. Then, the stochastic circulations are defined as
2.8
Although *N*^{+}(*t*) and *N*^{−}(*t*) are random variables, the limits in equation (2.8) exist.

We can consider the dynamics on the circle as movement on a line with the endpoints being joined. Then, the and represent the probability of exiting the line through the left end (−2*π*), or the right end (2*π*), respectively, divided by the mean wait time, 〈*T*〉. The exit probabilities can be expressed in terms of the function (Gardiner 1985)
2.9
The probability of exit through either endpoint starting at an initial point *Θ*_{0} is then
2.10
Starting initially at *Θ*_{0}=0, the ratio and difference between the two probabilities are
2.11
Note that, when *A*=0, we have *π*^{+}(0)=0 and *π*^{−}(0)=1. Thus, deterministically there is a clear direction of rotation.

The mean wait time to either endpoint is
2.12
The Hill–Qian circulation numbers are *ϕ*^{+}_{HQ}=*π*^{+}(0)/〈*T*〉 and *ϕ*^{−}_{HQ}=*π*^{−}(0)/〈*T*〉. The net Hill–Qian circulation, *ϕ*^{+}_{HQ}−*ϕ*^{−}_{HQ}, also called the stochastic rotation number, can also be obtained as the stationary flux on the circle (Wang *et al.* 1997).

Figure 1 shows a plot of the net Hill–Qian circulation versus the damping parameter *a*. When *a*=*ω*=1, we see the disappearance of oscillations in the deterministic case, *A*=0. For positive *A*, there is always net rotation.

### (d) Rotational random walk in two dimensions

We now consider a two-dimensional master equation system (Q-process) in polar coordinates. The symbols *n* and *m* in this section bear no relation to the *n* and *m* in the previous sections. We assume that the angular coordinate, *j*, is uniformly divided into *n* states, as shown in figure 2. Furthermore, we assume that all the rate constants are independent of *j*
2.13
The system has a rotational symmetry. With this set of rate constants, the motion in the radial direction is uncoupled to the motion in the angular direction.

Taking advantage of the rotational symmetry, we can solve the stationary probability problem of the master equation
2.14
The net cycle flux on each ring with radius *i* is
2.15
in which the term in is the total probability of the cycle with radius *i*. It is clear that, while the can be large, whether there is a significant net cycle flux also crucially depends on the term (*u*_{i}−*w*_{i})/*n*. The forward and backward one-way cycle fluxes obey the relations
2.16
Therefore,
2.17
We see that if *u*_{i}>*w*_{i} and *n* is large, then .

The above master equation system is the stochastic counterpart of the nonlinear dynamics in polar coordinates known as the *λ*–*ω* system in Murray (2002)
2.18
When *r*_{i}=*s*_{i}, we have *Ω*(*i*)=0. In this case, the master equation system has detailed balance; there is no net cycle flux. For a dynamical system with *Ω*≠0∀*i*, irrespective of the nature of *f*(*r*), there are complex eigenvalues, which correspond to either a stable fixed point of spiral type, or a stable limit cycle. In both cases, the system exhibits oscillatory dynamics with off-zero power spectral peak.

However, oscillatory dynamics does not necessarily imply that there is a stable limit cycle. A limit cycle has to have a certain ‘attractivity’, which cannot be represented through the power spectrum alone. The present paper focuses on this aspect of stochastic limit cycles. Essentially, all the previous studies on stochastic oscillations have been focused on the power spectrum (Xiao *et al.* 2007), with the only exception being Treutlein & Schulten (1985).

Recently, Ge and his co-workers (Ge *et al.* 2008; Ge & Qian 2009) have introduced the stochastic limit cycle based on the theory of circulation of Markov chain (Jiang *et al.* 2004). The Hill–Qian circulation decomposition theorem provides each and every directed circuit with a *circulation*. It is defined as the ratio of the net number of times the stationary process completes a circuit to the total run time. The Hill–Qian circulation along the ring with radius *i* is
2.19
Thus, with the independence of the motion in radial and angular directions, the Hill–Qian circulation along a ring is precisely the product of the total probability of the ring times the angular velocity. It is exactly the net cycle flux in equation (2.15).

These examples emphasize the importance of the cycle flux, as opposed to the probability steady state or the frequency (power spectrum) alone in understanding stochastic circulation. Without radial or rotational symmetry, the number of cycles in a two-dimensional system becomes very difficult to enumerate, and thus Hill’s cycle flux is almost impossible to obtain. In the next section, we introduce the PHCM, a natural extension of the Poincaré map into stochastic systems. The PHCM is an attempt to look into stochastic oscillation without necessarily calculating the cycle flux.

## 3. The PHCM and reversible Schnakenberg model

In this section, we define the PHCM and the reversible Schnakenberg model of oscillatory chemical kinetics, to which we will apply the map.

### (a) Introducing PHCM

We now combine the idea of cycle fluxes on a circle with the Poincaré return map from nonlinear dynamical systems theory. In a two-dimensional system, the Poincaré return map is a device for studying limit cycle oscillations. In the following construction, we assume that the system contains only one fixed point, and, when this fixed point is unstable, it is surrounded by a stable limit cycle.

The PHCM is constructed as follows.

(i) Define a horizontal line through the point (*n**,*m**)=*V* (*x**,*y**), where (*x**,*y**) is the fixed point of the corresponding ordinary differential equation (ODE) system and *V* is the volume size. The half line to the right of the fixed point will be referred to as the Poincaré cut.1

(ii) Starting with an initial value *n*_{0} on the Poincaré cut, simulate the CME using the Gillespie algorithm until the system returns to the line *after completing one cycle* in either the forward or the backward direction (see figure 3*a*).

(iii) If the cycle was completed in the counterclockwise direction, the point at which the trajectory intersects the Poincaré cut is called *P*^{+}(*n*_{0}). If the cycle was completed in the clockwise direction, the point is called *P*^{−}(*n*_{0}) (see figure 3*b*).

(iv) Because the system is stochastic, *P*^{+}(*n*_{0}) and *P*^{−}(*n*_{0}) are random variables. By repeating the simulation for the duration of one cycle with the same initial condition (*n*_{0}) many times, we can form the conditional histograms, *p*^{+}(*n*|*n*_{0}) and *p*^{−}(*n*|*n*_{0}), which represent the unnormalized probability density functions of *P*^{+}(*n*_{0}) and *P*^{−}(*n*_{0}), respectively.

(v) We are interested in the net circulation (as in §2*c*). The distribution for the PHCM with initial condition *n*_{0}, *p*(*n*|*n*_{0}), is the difference *p*(*n*|*n*_{0})=*p*^{+}(*n*|*n*_{0})−*p*^{−}(*n*|*n*_{0}) (or *p*^{−}(*n*|*n*_{0})−*p*^{+}(*n*|*n*_{0}), depending on the direction of deterministic rotation). Note that *g*(*n*|*n*_{0})=*p*^{+}(*n*|*n*_{0})+*p*^{−}(*n*|*n*_{0}) defines a Markov map, but *p*(*n*|*n*_{0}) does not.

In §4, we will apply this map to the reversible Schnakenberg model and discuss methods of analysis. First, we introduce the model and its behaviour as a system of ODEs.

### (b) The reversible Schnakenberg model

Schnakenberg’s original model is a simple trimolecular system able to support limit cycle behaviour (Schnakenberg 1979); see Murray (2002) and Qian *et al.* (2002) for further discussion. Macroscopic oscillation is intimately related to irreversibility of the mesoscopic dynamics; a consequence of external driving force in chemical dynamics (Qian 2007). We shall study an extension of the Schnakenberg model in which all the individual reactions are microscopically reversible, hence the driving force is well defined. This modified model is
3.1
in which specices *A* and *B* are fixed while *X* and *Y* are dynamic. In order to be consistent with the original model, the backward rates used in simulation will be relatively small compared with the forward rates. The consequence of the extension is that one can quantify the chemical driving force of the overall reaction, , with chemical potential difference .

We form a two-dimensional system of ODEs using the law of mass action. Here, *x*(*t*) and *y*(*t*) represent the concentrations of molecules *X* and *Y* , and *a* and *b* represent the concentrations of molecules *A* and *B*, which are held constant. The ODEs are
3.2
It can be shown that, for reasonable parameter values, there will be only one positive fixed point.

Figure 4*a* shows a phase plane solution to the ODE with reaction rates and fixed concentrations
3.3
Under these conditions, a stable limit cycle exists and a Poincaré map analysis is appropriate. The dotted line in figure 4*a* represents the line used for the Poincaré map (chosen to match the Poincaré cut described in §3) and the star represents the unstable fixed point, (*x**,*y**). Figure 4*b* plots this map, *P*_{ODE}(*x*_{0}), for initial conditions *x*_{0}≥*x**.

The Poincaré map has fixed points when *P*_{ODE}(*x*)=*x*, i.e. where the map intersects the diagonal. In figure 4*b*, the diagonal is drawn as a dotted line and there are two fixed points. Stability of each fixed point is determined by the derivative: if |*P*′_{ODE}(*x*)|>1 at the fixed point, it is unstable, and if |*P*′_{ODE}(*x*)|<1, it is stable. Thus, we see that the left fixed point, which represents the unstable fixed point in the ODE, is an unstable fixed point of the Poincaré map. The right fixed point is stable, and represents the intersection of the Poincaré line (the dotted line in figure 4*a*) and the limit cycle.

The Poincaré map in figure 4*b* is the infinite volume limit of the PHCM. There is no counterclockwise map *P*^{+}(*x*_{0}) because it is impossible for a deterministic system to flow against the direction of a vector field. In the infinite volume limit, the distributions for each *P*(*n*_{0}) have evolved to delta functions over the deterministic values *P*_{ODE}(*x*_{0}).

## 4. Application of the PHCM

In this section, we discuss results and insight from applying the PHCM to the Schnakenberg model. The CME formulation of the model, including the birth and death rates used in the Gillespie algorithm, can be found in the electronic supplementary material. Because the deterministic flow of the Schnakenberg model is in the clockwise direction, the net conditional distributions of the PHCM will be defined as *p*(*n*|*n*_{0})=*p*^{−}(*n*|*n*_{0})−*p*^{+}(*n*|*n*_{0}). Unless otherwise stated, steps 2 and 3 of the algorithm described in §3*a* were repeated 10 000 times to produce each distribution.

### (a) Steady-state distribution of the PHCM

For a Poincaré map, one is interestd in its fixed points. For the stochastic cycling, we now have *p*(*n*|*n*_{0})=*p*^{+}(*n*|*n*_{0})−*p*^{−}(*n*|*n*_{0}) and *g*(*n*|*n*_{0})=*p*^{+}(*n*|*n*_{0})+*p*^{−}(*n*|*n*_{0}). As *g*(*n*|*n*_{0}) is Markovian, it has a stationary distribution . Then, the steady-state distribution of the PHCM is defined as . Another way to compare the PHCM with the deterministic Poincaré map is to consider the positions *n**(*n*_{0}) at which *p*(*n*|*n*_{0}) is at its peak value when *n*=*n**. In general, there are multiple *n**(*n*_{0}) for a given *n*_{0}.

We applied this approach to the reversible Schnakenberg model, using the parameter values in equation (3.3), which lead to a stable limit cycle in the ODE. We found that, for each distribution *p*(*n*|*n*_{0}) with a given *n*_{0}, there was a ‘peak’ near the location of the limit cycle. Thus, we construct a curve: the location of the peak of the distribution for each initial condition. We then plot this curve against a diagonal line, and define the fixed point(s) as the intersection(s). In this way, results from the PHCM can be directly compared with those of the deterministic Poincaré map in figure 4*b*. Figure 5*a* plots the location of the peak of each distribution for several initial conditions. Results for small volume (*V* =100), large volume (*V* =900) and infinite volume (from figure 4*b*) are all plotted together in units of concentration.

Figure 5*b* is a plot of the distribution *p*(*n*|*n*_{0}) for the *V* =100, with *n*_{0}=30. It has a peak at *n*=30. To systematically locate the peak, a smoothing spline was used on each distribution and its critical points were numerically calculated. The spline parameters were kept constant for each volume size. The spline fit is shown in figure 5*b* as a dotted line.

In the small volume cases, we observed a ‘spike’ near 0 in addition to the smoother hump near the limit cycle. Throughout the following discussion, we refer to the former maxima as the ‘spike’ and latter maxima as the ‘peak’. The spline fit was used only to locate the peak; in our code we omitted the spike during fitting. We will discuss the meaning of the spike in §4*b*.

#### (i) p(n|n_{0}) with small and large V

For small *V* , the *p*(*n*|*n*_{0}) with different initial value *n*_{0} look very similar. As seen in figure 5*a*, the location of the peak increases slowly with increasing *n*_{0}. However, as *V* increases, so do the differences in the peak locations for different *n*_{0}. The peak values for *V* =900 are very close to the values of the deterministic Poincaré map.

Figure 6 illustrates the difference between small and large *V* with schematics. For small *V* , a distribution very similar to the steady-state distribution is obtained for a wide range of *n*_{0}. For large *V* , the peaks of the distributions change with *n*_{0}; distributions similar to the steady-state distribution can be obtained only for a restricted set of initial conditions.

The reason for this difference is as follows: in systems with small *V* , movement perpendicular to the limit cycle is much faster than the rotational movement. This leads to a ‘memoryless’ dynamics, where the system behaves the same as it cycles across the Poincaré line regardless of where it started. This indicates a rapid approach to steady-state distribution. Looking at figure 5*a*, we see that, at the point where the line intersects the diagonal, the slope is smaller for small *V* than for large *V* and ODE cases. The slope at the intersection point with the diagonal is a measure of the rate of convergence to the fixed point.

Thus, in small volume systems, there is a ‘fuzzy’ attractor, which is reached within one cycle regardless of the initial condition, but has a large variance around its peak. In large volume systems, multiple cycles are needed to approach the attractor, but, once it is reached, there is a very small variance around the peak.

### (b) Bistability in small volume systems

As mentioned in §4*a*, for small volume systems and parameters leading to a stable limit cycle, there are two places of high probability in the PHCM distributions: a spike near the unstable fixed point (represented by 0), and a second peak closer to the limit cycle. We expect the limit cycle to be attracting, but find that the unstable fixed point is also attracting, even when the initial condition is to the right of the limit cycle (i.e. the simulation is started outside the limit cycle). When the volume is increased, the spike near the unstable fixed point appears only for cases of initial conditions inside the limit cycle. As the volume is increased further, the spike disappears completely and the peak becomes more pronounced.

Figure 7 plots the log of the ratio of the spike height to the peak height against volume size for steady-state distributions of the PHCM, using the same parameters as in §4*a*. For small volume sizes, there is a kind of bistability in the PHCM; both the limit cycle and the fixed point are attracting. The appearance of bistability in small volume systems where there is none in the large volume equivalent has been observed in other models as well (Bishop & Qian in press). There is a clear exponential relationship between the spike to peak ratio and volume size (a linear fit is shown as a dotted line in figure 7). The spike is strictly a small volume phenomenon.

### (c) Case of a stable fixed point

We also applied the PHCM to the Schnakenberg model using the same reaction rates as in equation (3.3), but with new pump parameter values: *a*=0.25 and *b*=0.4. These parameters lead to a stable fixed point in the ODE, but are near the Hopf bifurcation.

We saw that the PHCM (with *V* =100) yields distributions which again have a spike and a peak. For these parameters, the spike near 0 is expected, but there is also a clear peak to the right of the stable fixed point. This peak represents a limit cycle which is not present in the ODE. It cannot be seen in the deterministic equations, but the vector field is such that cycles are occurring away from the stable fixed point. Through collaboration with Cao & Liang (2008), we obtained the two-dimensional steady-state distribution for this system with *V* =100. The steady-state distribution shows a single maximum around the stable fixed point; there is no evidence of a ‘limit cycle’. This illustrates the importance of the PHCM, which takes into account steady-state behaviour and cyclic movement.

There is a clear Hopf bifurcation in the ODE in which a limit cycle grows out of a single fixed point. When applying the PHCM with a certain initial condition, the analogous bifurcation in the CME is the point at which a peak appears in the distribution. As one or more parameters are varied, the distribution changes from having only one maximum (near 0) to having an additional maximum to the right of 0. The latter represents a stochastic limit cycle and can appear even for parameters which lead to a stable fixed point in the ODE.

Figure 8*a* shows the *a, b* parameter plane. The bold line indicates the location of the Hopf bifurcation, which can be found analytically. The dotted line shows the approximate location of the analogous bifurcation in the PHCM with *V* =100. This bifurcation is defined as the parameter values at which there is no longer a peak in the spline fit of the distribution with initial condition *x*_{0}=50. As discussed before, when the volume is small, the distribution will be the same for most initial conditions, so the choice of *x*_{0}=50 is arbitrary.

Figure 8*b* shows two distributions on either side of the dotted line bifurcation in figure 8*a* (the distributions were normalized for better comparison). As the parameters near the bifurcation, the peak begins to flatten without necessarily moving closer to 0. This is another small volume phenomenon: when the limit cycle is born, it is relatively far away from the fixed point. This is unlike large volume (or ODE) behaviour in which the peak (or limit cycle) begins very close to the fixed point and moves outward as the parameters move away from the bifurcation point.

### (d) Convergence to the deterministic map

We applied the PHCM using the parameters in equation (3.3), with *b*=0.4 and allowing *a* to vary across the Hopf bifurcation, from *a*=0.1 to 0.3. We calculated the steady-state distribution for each parameter pair and located the peak. Figure 9 shows a plot of the peak in the PHCM steady-state distribution as the parameter *a* is varied. The peak is plotted for three different volume sizes (in terms of concentration) along with the location of the ODE limit cycle.

Figure 9 is reminiscent of figure 1 because it illustrates the presence of rotation in stochastic systems where there is none in the corresponding deterministic system. There are peaks which appear in the steady-state distributions for *a*>0.2 (when the ODE fixed point is stable) for all volume sizes, but as the volume increases, the peak moves closer to 0. This convergence is similar to the manner in which the net stochastic circulation in figure 1 approaches 0 as the noise decreases.

## 5. Discussion

A two-dimensional stochastic process with a stationary distribution will exhibit sustained rotations when the deterministic counterparts do not. However, the stochastic limit cycle is difficult to define, as noise allows the system to complete closed paths in a large number of ways. The mathematical theory of stochastic circulation for a Markov chain requires enumerating all the possible cycles and quantifying the net flux on each and every cycle (Hill 1989; Jiang *et al.* 2004). While this approach has generated interesting insights in systems with a small number of nodes (Ge *et al.* 2008; Ge & Qian 2009), the application of this theory to the CME systems in high dimension (*n*≥2), unfortunately, is not realistic.

The state space of a CME, however, has a well-defined geometry. Thus, one can introduce the PHCM, a generalization of the Poincaré return map from dynamical systems theory and Hill’s cycle flux theory. It is a concept which can be algorithmically applied to any stochastic trajectory. Using this tool, we can observe stochastic oscillations which are unseen by the steady-state probability distribution function. The PHCM differs from a power spectrum analysis; it focuses on the location of the oscillatory behaviour, not the frequency. The ability to locate the limit cycle gives additional information such as the amplitude of oscillations and can distinguish between sustained (limit cycle) and dissipating (focus) oscillations.

We applied the PHCM to a reversible extension of a classic model of chemical limit cycle, the Schnakenberg (1979) model. We used the PHCM to locate the limit cycle. Because Boolean networks possess no geometry, only topology, there is not a metric and therefore no concept of the geometry for a limit cycle (Ge *et al.* 2008; Ge & Qian 2009), as defined here.

When the volume is very small, the PHCM shows evidence of two attractors in the Schnakenberg system even though its deterministic counterpart has a single stable limit cycle. The two attractors, one very close to the unstable fixed point and one near the limit cycle of the ODE, represent the maxima in the steady-state distribution of the PHCM. As the volume increases, the relative heights of the two maxima shrink exponentially and only the limit cycle remains. This observation is akin to the stochastic bistability recently observed (Bishop & Qian in press).

In the case when ODE has a stable fixed point, the PHCM again reveals bistability when the volume is small. There is an attractive cycle in addition to the attractive fixed point. Unlike the birth of a limit cycle in a deterministic Hopf bifurcation, when the small volume cycle appears, it is already away from the fixed point. The PHCM shows evidence of circulation in stochastic systems over a region of the parameter space which is larger than the limit cycle region for the ODE system.

As the volume increases, we see an agreement with the ODE model: bistability disappears and the peak of the distribution closely matches the deterministic limit cycle. The distributions obtained from the PHCM become more dependent on the initial condition and more cycles are required to reach the steady state.

## Acknowledgements

The authors would like to thank Jie Liang and Cao Youfang for helpful discussions and computation of probability steady-state solutions.

## Footnotes

↵1 Although we choose a horizontal line here, as with the traditional Poincaré map, any line through the fixed point may be used. Furthermore, for systems with dimension

*n*, the cut is a*n*−1 dimensional hyperplane.- Received July 3, 2009.
- Accepted October 13, 2009.

- © 2009 The Royal Society